44 double a32,a31,a41,a42,a21, occnu_lya ,
45 rate12 , rate21 , pump12 , pump21 , coll12 , coll21,
46 texc , occnu_lya_23 , occnu_lya_13,occnu_lya_24,occnu_lya_14, texc1, texc2;
72 texc1 =
sexp(0.068/texc);
73 texc2 =
sexp(((82259.272-82258.907)*
T1CM)/texc);
84 occnu_lya_23 = occnu_lya;
85 occnu_lya_13 = occnu_lya*texc1;
86 occnu_lya_14 = occnu_lya_13*texc2;
87 occnu_lya_24 = occnu_lya*texc2;
108 3.*a31*occnu_lya_13 *a32/(a31+a32)+
111 3.*a41*occnu_lya_14 *a42/(a41+a42);
124 occnu_lya_23*a32 * a31/(a31+a32)+
125 occnu_lya_24*a42*a41/(a41+a42);
128 x = rate12 /
SDIV(a21 + rate21);
142 HFLines[0].Emis->xIntensity = HFLines[0].Emis->phots*HFLines[0].EnergyErg;
145 HFLines[0].Emis->ColOvTot = coll12 / rate12;
169 temp =
MIN2(1e4 , temp );
173 hold = -9.607 + log10( sqrt(temp)) *
sexp( pow(log10(temp) , 4.5 ) / 1800. );
174 hold = pow(10.,hold );
184 STATIC double h21_t_ge_20(
double temp )
190 temp =
MIN2( 1000.,temp );
192 y =-21.70880995483007-13.76259674006133*
x1;
198 if( teorginal > 1e3 )
200 y *= pow(teorginal/1e3 , 0.33 );
207 STATIC double h21_t_lt_20(
double temp )
213 temp =
MAX2( 1., temp );
215 y =9.720710314268267E-08+6.325515312006680E-08*
x1;
230 temp =
MIN2( 300., temp );
232 y =1.4341127e-9+9.4161077e-15*x1-9.2998995e-9/(log(x1))+6.9539411e-9/sqrt(x1)+1.7742293e-8*(log(x1))/pow(x1,2);
233 if( teorginal > 300. )
236 x3 =
MIN2( 1000., teorginal );
238 y =-21.70880995483007-13.76259674006133*
x2;
242 if( teorginal > 1e3 )
245 y *= pow(teorginal/1e3 , 0.33 );
256 temp =
MAX2(1., temp );
258 y =8.5622857e-10+2.331358e-11*x1+9.5640586e-11*pow((log(x1)),2)-4.6220869e-10*sqrt(x1)-4.1719545e-10/sqrt(x1);
310 temp =
MAX2( 2. , temp );
311 temp =
MIN2( 2e4 , temp );
319 hold =9.588389834316704E-11 - 5.158891920816405E-14*x1
320 +5.895348443553458E-19*x2 + 2.053049602324290E-11*x3
321 +9.122617940315725E-10*x4;
347 double strengths[12];
360 long int i, j, mass, nelec, ion, nelem;
372 fprintf(
ioQQQ,
" Hyperfine opening hyperfine.dat:");
374 ioDATA =
open_data(
"hyperfine.dat",
"r" );
379 fprintf(
ioQQQ,
" Hyperfine could not read first line of hyperfine.dat.\n");
385 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL )
389 if( chLine[0] !=
'#')
408 if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
410 fprintf(
ioQQQ,
" Hyperfine could not rewind hyperfine.dat.\n");
417 fprintf(
ioQQQ,
" Hyperfine could not read first line of hyperfine.dat.\n");
429 static int iYR=2, iMN=10, iDY=28;
430 if( ( nelem != iYR ) || ( nelec != iMN ) || ( ion != iDY ) )
433 " Hyperfine: the version of hyperfine.dat in the data directory is not the current version.\n" );
435 " I expected to find the number %i %i %i and got %li %li %li instead.\n" ,
437 nelem , nelec , ion );
438 fprintf(
ioQQQ,
"Here is the line image:\n==%s==\n", chLine );
455 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL )
459 if( chLine[0] !=
'#')
465 sscanf(chLine,
"%li%i%le%i%le%le%le%le%le%le%le%le%le%le%le%le%le%le%le",
523 q12 = HFLines[i].Hi->g/ HFLines[i].Lo->g * q21 * exp(-1 * h * c * HFLines[i].EnergyWN / (k *
phycon.
te));
525 x = Ne * q12 / (HFLines[i].Emis->Aul * (1 + Ne * q21 / HFLines[i].Aul));
526 HFLines[i].xIntensity =
N * HFLines[i].Emis->Aul * x / (1.0 + x) * h * c / (HFLines[i].WLAng / 1e8);
541 # define N_TE_TABLE 12
542 double slope, upsilon, TemperatureTable[
N_TE_TABLE] = {.1e6, .15e6, .25e6, .4e6, .6e6,
543 1.0e6, 1.5e6, 2.5e6, 4e6, 6e6, 10e6, 15e6};
558 if(
phycon.
te <= TemperatureTable[0])
585 else if(
phycon.
te < TemperatureTable[j])
587 slope = (log10(Strengths[i].strengths[j-1]) - log10(Strengths[i].strengths[j])) /
588 (log10(TemperatureTable[j-1]) - log10(TemperatureTable[j]));
590 upsilon = (log10((
phycon.
te)) - (log10(TemperatureTable[j-1])-6.))*slope + log10(Strengths[i].strengths[j-1]);
591 upsilon = pow(10., upsilon);