cloudy trunk
Loading...
Searching...
No Matches
atom_hyperfine.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3/* HyperfineCreat establish space for hf arrays, reads atomic data from hyperfine.dat */
4/* HyperfineCS - returns collision strengths for hyperfine struc transitions */
5/*H21cm computes rate for H 21 cm from upper to lower excitation by atomic hydrogen */
6/*h21_t_ge_20 compute rate for H 21 cm from upper to lower excitation by atomic hydrogen */
7/*h21_t_lt_20 compute rate for H 21 cm from upper to lower excitation by atomic hydrogen */
8/*H21cm_electron compute H 21 cm rate from upper to lower excitation by electrons - call by CoolEvaluate */
9/*H21cm_H_atom - evaluate H atom spin changing collision rate, called by CoolEvaluate */
10/*H21cm_proton - evaluate proton spin changing H atom collision rate, */
11#include "cddefines.h"
12#include "conv.h"
13#include "lines_service.h"
14#include "phycon.h"
15#include "dense.h"
16#include "rfield.h"
17#include "taulines.h"
18#include "iso.h"
19#include "trace.h"
20#include "hyperfine.h"
21#include "physconst.h"
22
23/* H21_cm_pops - fine level populations for 21 cm with Lya pumping included
24 * called in CoolEvaluate */
25void H21_cm_pops( void )
26{
27 /*atom_level2( HFLines[0] );*/
28 /*return;*/
29 /*
30 things we know on entry to this routine:
31 total population of 2p: iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH2p].Pop
32 total population of 1s: iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop
33 continuum pumping rate (lo-up) inside 21 cm line: HFLines[0].pump()
34 upper to lower collision rate inside 21 cm line: HFLines[0].cs*dense.cdsqte
35 occupation number inside Lya: OccupationNumberLine( &iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s) )
36
37 level populations (cm-3) must be computed:
38 population of upper level of 21cm: HFLines[0].Hi->Pop
39 population of lower level of 21cm: (*HFLines[0].Lo()).Pop
40 stimulated emission corrected population of lower level: HFLines[0].Emis->PopOpc()
41 */
42
43 double x,
44 PopTot;
45 double a32,a31,a41,a42,a21, occnu_lya ,
46 rate12 , rate21 , pump12 , pump21 , coll12 , coll21,
47 texc , occnu_lya_23 , occnu_lya_13,occnu_lya_24,occnu_lya_14, texc1, texc2;
48
49 PopTot = iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop();
50 /* population can be zero in certain tests where H is turned off,
51 * also if initial solver does not see any obvious source of ionization
52 * also possible to set H0 density to zero with element ionization command,
53 * as is done in func_set_ion test case */
54 if( PopTot <0 )
56 else if( PopTot == 0 )
57 {
58 /*return after zeroing local variables */
59 (*HFLines[0].Hi()).Pop() = 0.;
60 (*HFLines[0].Lo()).Pop() = 0.;
61 HFLines[0].Emis().PopOpc() = 0.;
62 HFLines[0].Emis().phots() = 0.;
63 HFLines[0].Emis().xIntensity() = 0.;
64 HFLines[0].Emis().ColOvTot() = 0.;
65 hyperfine.Tspin21cm = 0.;
66 return;
67 }
68
69 a31 = 2.08e8; /* Einstein co-efficient for transition 1p1/2 to 0s1/2 */
70 a32 = 4.16e8; /* Einstein co-efficient for transition 1p1/2 to 1s1/2 */
71 a41 = 4.16e8; /* Einstein co-efficient for transition 1p3/2 to 0s1/2 */
72 a42 = 2.08e8; /* Einstein co-efficient for transition 1p3/2 to 1s1/2 */
73 /* These A values are determined from eqn. 17.64 of "The theory of Atomic structure
74 * and Spectra" by R. D. Cowan
75 * A hyperfine level has degeneracy Gf=(2F + 1)
76 * a2p1s = 6.24e8; Einstein co-efficient for transition 2p to 1s */
77 a21 = 2.85e-15; /* Einstein co-efficient for transition 1s1/2 to 0s1/2 */
78
79 /* above is spontaneous rate - the net rate is this times escape and destruction
80 * probabilities */
81 a21 *= (HFLines[0].Emis().Pdest() + HFLines[0].Emis().Pesc() + HFLines[0].Emis().Pelec_esc());
82 ASSERT( a21>0. );
83
84 /* hyperfine.lgLya_pump_21cm is option to turn off Lya pump
85 * of 21 cm, with no 21cm lya pump command - note that this
86 * can be negative if Lya mases - can occur during search phase */
88 hyperfine.lgLya_pump_21cm;
89 if( occnu_lya<0 )
90 {
91 static bool lgCommentDone = false;
92 /* Lya is masing - could be due to very bad solution in search phase */
93 if( !conv.lgSearch && !lgCommentDone )
94 {
95 fprintf(ioQQQ,
96 "NOTE Lya masing will invert 21 cm, occupation number set zero\n");
97 lgCommentDone = true;
98 }
99 occnu_lya = 0.;
100 }
101
102 /* Lya occupation number for the hyperfine levels 0S1/2 and 1S1/2 are different*/
103 texc = TexcLine( iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s) );
104 /* Energy difference between 2p1/2 and 2p3/2 taken from NSRDS */
105 if( texc > 0. )
106 {
107 /* convert to Boltzmann factor, which will applied to occupation
108 * number of higher energy transition */
109 texc1 = sexp(0.068/texc);
110 texc2 = sexp(((82259.272-82258.907)*T1CM)/texc);
111 }
112 else
113 {
114 texc1 = 0.;
115 texc2 = 0.;
116 }
117
118 /* the continuum within Lya seen by the two levels is not exactly the same brightness. They
119 * differ by the exp when Lya is on Wein tail of black body, which must be true if 21 cm is important */
120
121 occnu_lya_23 = occnu_lya;
122 occnu_lya_13 = occnu_lya*texc1;
123 occnu_lya_14 = occnu_lya_13*texc2;
124 occnu_lya_24 = occnu_lya*texc2;
125
126 /* this is the 21 cm upward continuum pumping rate [s-1] for the attenuated incident and
127 * local continuum and including line optical depths */
128 pump12 = HFLines[0].Emis().pump();
129 pump21 = pump12 * (*HFLines[0].Lo()).g() / (*HFLines[0].Hi()).g();
130
131 /* collision rates s-1 within 1s,
132 * were multiplied by collider density when evaluated in CoolEvaluate */
133 /* ContBoltz is Boltzmann factor for wavelength of line */
134 ASSERT( HFLines[0].Coll().col_str()>0. );
135 coll12 = HFLines[0].Coll().col_str()*dense.cdsqte/(*HFLines[0].Lo()).g()*rfield.ContBoltz[HFLines[0].ipCont()-1];
136 coll21 = HFLines[0].Coll().col_str()*dense.cdsqte/(*HFLines[0].Hi()).g();
137
138 /* set up rate (s-1) equations
139 * all process out of 1 that eventually go to 2 */
140 rate12 =
141 /* collision rate (s-1) from 1 to 2 */
142 coll12 +
143 /* direct external continuum pumping (s-1) in 21 cm line - usually dominated by CMB */
144 pump12 +
145 /* pump rate (s-1) up to 3, times fraction that decay to 2, hence net 1-2 */
146 3.*a31*occnu_lya_13 *a32/(a31+a32)+
147 /* pump rate (s-1) up to 4, times fraction that decay to 2, hence net 1-2 */
148 /* >>chng 05 apr 04, GS, degeneracy corrected from 6 to 3 */
149 3.*a41*occnu_lya_14 *a42/(a41+a42);
150
151 /* set up rate (s-1) equations
152 * all process out of 2 that eventually go to 1 */
153 /* spontaneous + induced 2 -> 1 by external continuum inside 21 cm line */
154 /* >>chng 04 dec 03, do not include spontaneous decay, for numerical stability */
155 rate21 =
156 /* collisional deexcitation */
157 coll21 +
158 /* net spontaneous decay plus external continuum pumping in 21 cm line */
159 pump21 +
160 /* rate from 2 to 3 time fraction that go back to 1, hence net 2 - 1 */
161 /* >>chng 05 apr 04,GS, degeneracy corrected from 2 to unity */
162 occnu_lya_23*a32 * a31/(a31+a32)+
163 occnu_lya_24*a42*a41/(a41+a42);
164
165 /* x = (*HFLines[0].Hi()).Pop/(*HFLines[0].Lo()).Pop */
166 x = rate12 / SDIV(a21 + rate21);
167
168 /* the Transitions term is the total population of 1s */
169 (*HFLines[0].Hi()).Pop() = (x/(1.+x))* PopTot;
170 (*HFLines[0].Lo()).Pop() = (1./(1.+x))* PopTot;
171 ASSERT( (*HFLines[0].Hi()).Pop() >0. );
172 ASSERT( (*HFLines[0].Lo()).Pop() >0. );
173
174 /* the population with correction for stimulated emission */
175 HFLines[0].Emis().PopOpc() = (*HFLines[0].Lo()).Pop()*((3*rate21- rate12) + 3*a21)/SDIV(3*(a21+ rate21));
176
177 /* number of escaping line photons, used elsewhere for outward beam */
178 HFLines[0].Emis().phots() = (*HFLines[0].Hi()).Pop() * HFLines[0].Emis().Aul() *
179 (HFLines[0].Emis().Pesc() + HFLines[0].Emis().Pelec_esc());
180 ASSERT( HFLines[0].Emis().phots() >= 0. );
181 /* intensity of line */
182 HFLines[0].Emis().xIntensity() = HFLines[0].Emis().phots()*HFLines[0].EnergyErg();
183
184 /* ratio of collisional to total (collisional + pumped) excitation */
185 HFLines[0].Emis().ColOvTot() = coll12 / rate12;
186
187 /* finally save the spin temperature */
188 if( (*HFLines[0].Hi()).Pop() > SMALLFLOAT )
189 {
190 hyperfine.Tspin21cm = TexcLine( HFLines[0] );
191 /* this line must be non-zero - it does strongly mase in limit_compton_hi_t sim -
192 * in that sim pop ratio goes to unity for a float and TexcLine ret zero */
193 if( hyperfine.Tspin21cm == 0. )
194 hyperfine.Tspin21cm = phycon.te;
195 }
196 else
197 {
198 hyperfine.Tspin21cm = phycon.te;
199 }
200
201 return;
202}
203
204/*H21cm_electron computes rate for H 21 cm from upper to lower excitation by electrons - call by CoolEvaluate
205 * >>refer H1 CS Smith, F. J. 1966, Planet. Space Sci., 14, 929 */
206double H21cm_electron( double temp )
207{
208 double hold;
209 temp = MIN2(1e4 , temp );
210 /* following fit is from */
211 /* >>refer H1 21cm Liszt, H. 2001, A&A, 371, 698 */
212
213 hold = -9.607 + log10( sqrt(temp)) * sexp( pow(log10(temp) , 4.5 ) / 1800. );
214 hold = pow(10.,hold );
215 return( hold );
216}
217
218/* computes rate for H 21 cm from upper to lower excitation by atomic hydrogen
219 * from
220 * >>refer H1 CS Allison, A. C., & Dalgarno A. 1969, ApJ 158, 423 */
221/* the following is the best current survey of 21 cm excitation */
222/* >>refer H1 21cm Liszt, H. 2001, A&A, 371, 698 */
223#if 0
224STATIC double h21_t_ge_20( double temp )
225{
226 double y;
227 double x1,
228 teorginal = temp;
229 /* data go up to 1,000K must not go above this */
230 temp = MIN2( 1000.,temp );
231 x1 =1.0/sqrt(temp);
232 y =-21.70880995483007-13.76259674006133*x1;
233 y = exp(y);
234
235 /* >>chng 02 feb 14, extrapolate above 1e3 K as per Liszt 2001 recommendation
236 * page 699 of */
237 /* >>refer H1 21cm Liszt, H. 2001, A&A, 371, 698 */
238 if( teorginal > 1e3 )
239 {
240 y *= pow(teorginal/1e3 , 0.33 );
241 }
242
243 return( y );
244}
245
246/* this branch for T < 20K, data go down to 1 K */
247STATIC double h21_t_lt_20( double temp )
248{
249 double y;
250 double x1;
251
252 /* must not go below 1K */
253 temp = MAX2( 1., temp );
254 x1 =temp*log(temp);
255 y =9.720710314268267E-08+6.325515312006680E-08*x1;
256 return(y*y);
257}
258#endif
259
260/* >> chng 04 dec 15, GS. The fitted rate co-efficients (cm3s-1) in the temperature range 1K to 300K is from
261 * >>refer H1 CS Zygelman, B. 2005, ApJ, 622, 1356
262 * The rate is 4/3 times the Dalgarno (1969) rate for the
263 temperature range 300K to 1000K. Above 1000K, the rate is extrapolated according to Liszt 2001.*/
264STATIC double h21_t_ge_10( double temp )
265{
266 double y;
267 double x1,x2,x3,
268 teorginal = temp;
269 /* data go up to 300K */
270 temp = MIN2( 300., temp );
271 x1 =temp;
272 y =1.4341127e-9+9.4161077e-15*x1-9.2998995e-9/(log(x1))+6.9539411e-9/sqrt(x1)+1.7742293e-8*(log(x1))/pow2(x1);
273 if( teorginal > 300. )
274 {
275 /* data go up to 1000*/
276 x3 = MIN2( 1000., teorginal );
277 x2 =1.0/sqrt(x3);
278 y =-21.70880995483007-13.76259674006133*x2;
279 y = 1.236686*exp(y);
280
281 }
282 if( teorginal > 1e3 )
283 {
284 /*data go above 1000*/
285 y *= pow(teorginal/1e3 , 0.33 );
286 }
287 return( y );
288}
289/* this branch for T < 10K, data go down to 1 K */
290STATIC double h21_t_lt_10( double temp )
291{
292 double y;
293 double x1;
294
295 /* must not go below 1K */
296 temp = MAX2(1., temp );
297 x1 =temp;
298 y =8.5622857e-10+2.331358e-11*x1+9.5640586e-11*pow2(log(x1))-4.6220869e-10*sqrt(x1)-4.1719545e-10/sqrt(x1);
299 return(y);
300}
301
302/*H21cm_H_atom - evaluate H atom spin changing H atom collision rate,
303 * called by CoolEvaluate
304 * >>refer H1 CS Allison, A. C. & Dalgarno, A. 1969, ApJ 158, 423
305 */
306double H21cm_H_atom( double temp )
307{
308 double hold;
309 if( temp >= 10. )
310 {
311 hold = h21_t_ge_10( temp );
312 }
313 else
314 {
315 hold = h21_t_lt_10( temp );
316 }
317
318 return hold;
319}
320
321/*H21cm_proton - evaluate proton spin changing H atom collision rate,
322* called by CoolEvaluate */
323double H21cm_proton( double temp )
324{
325 /*>>refer 21cm p coll Furlanetto, S. R. & Furlanetto, M. R. 2007, MNRAS, 379, 130
326 * previously had used proton rate, which is 3.2 times H0 rate according to
327 *>>refer 21cm CS Liszt, H. 2001, A&A, 371, 698 */
328 /* fit to table 1 of first paper */
329 /*--------------------------------------------------------------*
330 TableCurve Function: c:\storage\litera~1\21cm\atomic~1\p21cm.c Jun 20, 2007 3:37:50 PM
331 proton coll deex
332 X= temperature (K)
333 Y= rate coefficient (1e-9 cm3 s-1)
334 Eqn# 4419 y=a+bx+cx^2+dx^(0.5)+elnx/x
335 r2=0.9999445384690351
336 r2adj=0.9999168077035526
337 StdErr=5.559328579039901E-12
338 Fstat=49581.16793656295
339 a= 9.588389834316704E-11
340 b= -5.158891920816405E-14
341 c= 5.895348443553458E-19
342 d= 2.05304960232429E-11
343 e= 9.122617940315725E-10
344 *--------------------------------------------------------------*/
345
346 double hold;
347 /* only fit this range, did not include T = 1K point which
348 * causes an inflection */
349 temp = MAX2( 2. , temp );
350 temp = MIN2( 2e4 , temp );
351
352 /* within range of fitted rate coefficients */
353 double x1,x2,x3,x4;
354 x1 = temp;
355 x2 = temp*temp;
356 x3 = sqrt(temp);
357 x4 = log(temp)/temp;
358 hold =9.588389834316704E-11 - 5.158891920816405E-14*x1
359 +5.895348443553458E-19*x2 + 2.053049602324290E-11*x3
360 +9.122617940315725E-10*x4;
361
362 return hold;
363}
364
365/*
366 * HyperfineCreate, HyperfineCS written July 2001
367 * William Goddard for Gary Ferland
368 * This code calculates line intensities for known
369 * hyperfine transitions.
370 */
371
372/* two products, the transition structure HFLines, which contains all information for the lines,
373 * and nHFLines, the number of these lines.
374 *
375 * these are in taulines.h
376 *
377 * info to create them contained in hyperfine.dat
378 *
379 * abundances of nuclei are also in hyperfine.dat, stored in
380 */
381
382/* Ion contains twelve varying temperatures, specified above, used for */
383/* calculating collision strengths. */
384typedef struct
385{
386 double strengths[12];
387} Ion;
388
390
391/* HyperfineCreate establish space for hf arrays, reads atomic data from hyperfine.dat */
393{
394 FILE *ioDATA;
395 char chLine[INPUT_LINE_LENGTH];
396 bool lgEOL;
397 /*double c, h, k, N, Ne, q12, q21, upsilon, x;*/
398 realnum spin, wavelength;
399 long int i, j, mass, nelec, ion, nelem;
400
401 DEBUG_ENTRY( "HyperfineCreate()" );
402
403 /* list of ion collision strengths for the temperatures listed in table */
404 /* HFLines containing all the data in Hyperfine.dat, and transition is */
405 /* defined in cddefines.h */
406
407 /*transition *HFLines;*/
408
409 /* get the line data for the hyperfine lines */
410 if( trace.lgTrace )
411 fprintf( ioQQQ," Hyperfine opening hyperfine.dat:");
412
413 ioDATA = open_data( "hyperfine.dat", "r" );
414
415 /* first line is a version number and does not count */
416 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
417 {
418 fprintf( ioQQQ, " Hyperfine could not read first line of hyperfine.dat.\n");
420 }
421 /* count how many lines are in the file, ignoring all lines
422 * starting with '#' */
423 nHFLines = 0;
424 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
425 {
426 /* we want to count the lines that do not start with #
427 * since these contain data */
428 if( chLine[0] != '#')
429 ++nHFLines;
430 }
431
432 /* allocate the transition HFLines array */
433 HFLines.resize(nHFLines);
434 AllTransitions.push_back(HFLines);
435 /* initialize array to impossible values to make sure eventually done right */
436 for( i=0; i< nHFLines; ++i )
437 {
438 HFLines[i].Junk();
439 HFLines[i].AddHiState();
440 HFLines[i].AddLoState();
441 HFLines[i].AddLine2Stack();
442 }
443
444 Strengths = (Ion *)MALLOC( (size_t)(nHFLines)*sizeof(Ion) );
445 hyperfine.HFLabundance = (realnum *)MALLOC( (size_t)(nHFLines)*sizeof(realnum) );
446
447 /* now rewind the file so we can read it a second time*/
448 if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
449 {
450 fprintf( ioQQQ, " Hyperfine could not rewind hyperfine.dat.\n");
452 }
453
454 /* check that magic number is ok, read the line */
455 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
456 {
457 fprintf( ioQQQ, " Hyperfine could not read first line of hyperfine.dat.\n");
459 }
460
461 /* check that magic number is ok, scan it in */
462 i = 1;
463 nelem = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
464 nelec = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
465 ion = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
466
467 {
468 /* the following is the set of numbers that appear at the start of hyperfine.dat 01 07 12 */
469 static int iYR=2, iMN=10, iDY=28;
470 if( ( nelem != iYR ) || ( nelec != iMN ) || ( ion != iDY ) )
471 {
472 fprintf( ioQQQ,
473 " Hyperfine: the version of hyperfine.dat in the data directory is not the current version.\n" );
474 fprintf( ioQQQ,
475 " I expected to find the number %i %i %i and got %li %li %li instead.\n" ,
476 iYR, iMN , iDY ,
477 nelem , nelec , ion );
478 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
480 }
481 }
482
483 /*
484 * scan the string taken from Hyperfine.dat, parsing into
485 * needed variables.
486 * nelem is the atomic number.
487 * IonStg is the ionization stage. Atom = 1, Z+ = 2, Z++ = 3, etc.
488 * Aul is used to find the einstein A in the function GetGF.
489 * most of the variables are floats.
490 */
491
492 /* this will count the number of lines */
493 j = 0;
494
495 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
496 {
497 /* only look at lines without '#' in first col */
498 /* make sure line in data file is < 140 char */
499 if( chLine[0] != '#')
500 {
501 double Aul;
502 ASSERT( j < nHFLines );
503
504 double help[3];
505 sscanf(chLine, "%li%i%le%i%le%le%le%le%le%le%le%le%le%le%le%le%le%le%le",
506 &mass,
507 &(*HFLines[j].Hi()).nelem(),
508 &help[0],
509 &(*HFLines[j].Hi()).IonStg(),
510 &help[1],
511 &help[2],
512 &Aul,
513 &Strengths[j].strengths[0],
514 &Strengths[j].strengths[1],
515 &Strengths[j].strengths[2],
516 &Strengths[j].strengths[3],
517 &Strengths[j].strengths[4],
518 &Strengths[j].strengths[5],
519 &Strengths[j].strengths[6],
520 &Strengths[j].strengths[7],
521 &Strengths[j].strengths[8],
522 &Strengths[j].strengths[9],
523 &Strengths[j].strengths[10],
524 &Strengths[j].strengths[11]);
525 spin = (realnum)help[0];
526 hyperfine.HFLabundance[j] = (realnum)help[1];
527 wavelength = (realnum)help[2];
528 HFLines[j].Emis().Aul() = (realnum)Aul;
529 HFLines[j].Emis().damp() = 1e-20f;
530 (*HFLines[j].Hi()).g() = (realnum) (2*(spin + .5) + 1);
531 (*HFLines[j].Lo()).g() = (realnum) (2*(spin - .5) + 1);
532 HFLines[j].WLAng() = wavelength * 1e8f;
533 HFLines[j].EnergyWN() = (realnum) (1. / wavelength);
534 HFLines[j].Emis().gf() = (realnum)(GetGF(HFLines[j].Emis().Aul(), HFLines[j].EnergyWN(), (*HFLines[j].Hi()).g()));
535 /*fprintf(ioQQQ,"HFLinesss1\t%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
536 j,HFLines[j].opacity() , HFLines[j].Emis().gf() , HFLines[j].Emis().Aul() , HFLines[j].EnergyWN,HFLines[j].Lo()->g());*/
537 (*HFLines[j].Lo()).nelem() = (*HFLines[j].Hi()).nelem();
538 (*HFLines[j].Lo()).IonStg() = (*HFLines[j].Hi()).IonStg();
539
540
541 ASSERT(HFLines[j].Emis().gf() > 0.);
542 /* increment the counter */
543 j++;
544 }
545
546 }
547 ASSERT( j==nHFLines );
548
549 /* closing the file */
550 fclose(ioDATA);
551
552# if 0
553 /* for debugging and developing only */
554 /* calculating the luminosity for each isotope */
555 for(i = 0; i < nHFLines; i++)
556 {
557 N = dense.xIonDense[(*HFLines[i].Hi()).nelem()-1][(*HFLines[i].Hi()).IonStg()-1];
558 Ne = dense.eden;
559
560 h = 6.626076e-27; /* erg * sec */
561 c = 3e10; /* cm / sec */
562 k = 1.380658e-16; /* erg / K */
563
564 upsilon = HyperfineCS(i);
565 /*statistical weights must still be identified */
566 q21 = COLL_CONST * upsilon / (phycon.sqrte * (*HFLines[i].Hi()).g());
567
568 q12 = (*HFLines[i].Hi()).g()/ (*HFLines[i].Lo()).g() * q21 * exp(-1 * h * c * HFLines[i].EnergyWN / (k * phycon.te));
569
570 x = Ne * q12 / (HFLines[i].Emis().Aul() * (1 + Ne * q21 / HFLines[i].Aul()));
571 HFLines[i].xIntensity() = N * HFLines[i].Emis().Aul() * x / (1.0 + x) * h * c / (HFLines[i].WLAng() / 1e8);
572
573 }
574# endif
575
576 return;
577}
578
579
580/*HyperfineCS returns interpolated collision strength for element nelem and ion ion */
581double HyperfineCS( long i )
582{
583
584 /*Search HFLines to find i of ion */
585 int /*i = 0,*/ j = 0;
586# define N_TE_TABLE 12
587 double slope, upsilon, TemperatureTable[N_TE_TABLE] = {.1e6, .15e6, .25e6, .4e6, .6e6,
588 1.0e6, 1.5e6, 2.5e6, 4e6, 6e6, 10e6, 15e6};
589
590 DEBUG_ENTRY( "HyperfineCS()" );
591
592 /*while(i < nHFLines+1 && (*HFLines[i].Hi()).nelem() != nelem && (*HFLines[i].Hi()).IonStg() != ion)
593 ++i;*/
594 ASSERT( i >= 0. && i <= nHFLines );
595
596 /*j = temperature */
597 /* calculate actual cooling rate for isotope. */
598 /* The temperature-dependent collision strength must first be calculated. */
599 /* phycon.te is compared to the first, last and intermediate values of strengths[]*/
600 /* and interpolated by taking the log of the collision strength */
601 /* and temperature. */
602
603 if( phycon.te <= TemperatureTable[0])
604 {
605 /* temperature below bounds of table */
606 j = 0;
607 slope = (log10(Strengths[i].strengths[j+1]) - log10(Strengths[i].strengths[j])) /
608 (log10(TemperatureTable[j+1]) - log10(TemperatureTable[j]));
609 upsilon = (log10((phycon.te)) - (log10(TemperatureTable[j])))*slope + log10(Strengths[i].strengths[j]);
610 upsilon = pow(10., upsilon);
611 }
612 else if( phycon.te >= TemperatureTable[N_TE_TABLE-1])
613 {
614 /* temperature above bounds of table */
615 j = N_TE_TABLE - 1;
616 slope = (log10(Strengths[i].strengths[j-1]) - log10(Strengths[i].strengths[j])) /
617 (log10(TemperatureTable[j-1]) - log10(TemperatureTable[j]));
618 upsilon = (log10((phycon.te)) - (log10(TemperatureTable[j])))*slope + log10(Strengths[i].strengths[j]);
619 upsilon = pow(10., upsilon);
620 }
621 else
622 {
623 j = 1;
624 /* want Table[j-1] < te < Table[j] */
625 /*while ( phycon.te >= TemperatureTable[j] && j < N_TE_TABLE )*/
626 /*while ( TemperatureTable[j] < phycon.te && j < N_TE_TABLE )*/
627 while ( j < N_TE_TABLE && TemperatureTable[j] < phycon.te )
628 j++;
629
630 ASSERT( j >= 0 && j < N_TE_TABLE);
631 ASSERT( (TemperatureTable[j-1] <= phycon.te ) && (TemperatureTable[j] >= phycon.te) );
632
633 if( fp_equal( phycon.te, TemperatureTable[j] ) )
634 {
635 upsilon = Strengths[i].strengths[j];
636 }
637 /*phycon.te must be less than TemperatureTable[j], greater than TemperatureTable[j-1][0]*/
638 else if( phycon.te < TemperatureTable[j])
639 {
640 slope = (log10(Strengths[i].strengths[j-1]) - log10(Strengths[i].strengths[j])) /
641 (log10(TemperatureTable[j-1]) - log10(TemperatureTable[j]));
642
643 upsilon = (log10((phycon.te)) - (log10(TemperatureTable[j-1])))*slope + log10(Strengths[i].strengths[j-1]);
644 upsilon = pow(10., upsilon);
645 }
646 else
647 {
648 upsilon = Strengths[i].strengths[j-1];
649 }
650 }
651
652 return upsilon;
653}
static double x2[63]
static double x1[83]
#define N_TE_TABLE
STATIC double h21_t_lt_10(double temp)
double H21cm_electron(double temp)
double H21cm_H_atom(double temp)
double H21cm_proton(double temp)
double HyperfineCS(long i)
void H21_cm_pops(void)
void HyperfineCreate(void)
STATIC double h21_t_ge_10(double temp)
static Ion * Strengths
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
sys_float sexp(sys_float x)
Definition service.cpp:914
#define STATIC
Definition cddefines.h:97
T pow2(T a)
Definition cddefines.h:931
#define MIN2
Definition cddefines.h:761
#define MALLOC(exp)
Definition cddefines.h:501
const int INPUT_LINE_LENGTH
Definition cddefines.h:254
#define EXIT_FAILURE
Definition cddefines.h:140
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition service.cpp:381
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition service.cpp:70
NORETURN void TotalInsanity(void)
Definition service.cpp:886
sys_float SDIV(sys_float x)
Definition cddefines.h:952
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
t_conv conv
Definition conv.cpp:5
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition cpu.cpp:625
const realnum SMALLFLOAT
Definition cpu.h:191
t_dense dense
Definition dense.cpp:24
t_hyperfine hyperfine
Definition hyperfine.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
const int ipH1s
Definition iso.h:27
const int ipH2p
Definition iso.h:29
const int ipH_LIKE
Definition iso.h:62
double GetGF(double trans_prob, double enercm, double gup)
static realnum * wavelength
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double T1CM
Definition physconst.h:167
UNUSED const double COLL_CONST
Definition physconst.h:229
t_rfield rfield
Definition rfield.cpp:8
static double ** col_str
Definition species2.cpp:29
static double * g
Definition species2.cpp:28
double strengths[12]
TransitionList HFLines("HFLines", &AnonStates)
long int nHFLines
Definition taulines.cpp:31
vector< TransitionList > AllTransitions
Definition taulines.cpp:8
static const int N
t_trace trace
Definition trace.cpp:5
double OccupationNumberLine(const TransitionProxy &t)
double TexcLine(const TransitionProxy &t)