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;
56 else if( PopTot == 0 )
61 HFLines[0].Emis().PopOpc() = 0.;
63 HFLines[0].Emis().xIntensity() = 0.;
64 HFLines[0].Emis().ColOvTot() = 0.;
91 static bool lgCommentDone =
false;
93 if( !
conv.lgSearch && !lgCommentDone )
96 "NOTE Lya masing will invert 21 cm, occupation number set zero\n");
109 texc1 =
sexp(0.068/texc);
110 texc2 =
sexp(((82259.272-82258.907)*
T1CM)/texc);
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;
128 pump12 =
HFLines[0].Emis().pump();
129 pump21 = pump12 * (*
HFLines[0].Lo()).g() / (*
HFLines[0].Hi()).g();
146 3.*a31*occnu_lya_13 *a32/(a31+a32)+
149 3.*a41*occnu_lya_14 *a42/(a41+a42);
162 occnu_lya_23*a32 * a31/(a31+a32)+
163 occnu_lya_24*a42*a41/(a41+a42);
166 x = rate12 /
SDIV(a21 + rate21);
169 (*
HFLines[0].Hi()).Pop() = (x/(1.+x))* PopTot;
170 (*
HFLines[0].Lo()).Pop() = (1./(1.+x))* PopTot;
175 HFLines[0].Emis().PopOpc() = (*
HFLines[0].Lo()).Pop()*((3*rate21- rate12) + 3*a21)/
SDIV(3*(a21+ rate21));
185 HFLines[0].Emis().ColOvTot() = coll12 / rate12;
209 temp =
MIN2(1e4 , temp );
213 hold = -9.607 + log10( sqrt(temp)) *
sexp( pow(log10(temp) , 4.5 ) / 1800. );
214 hold = pow(10.,hold );
224STATIC double h21_t_ge_20(
double temp )
230 temp =
MIN2( 1000.,temp );
232 y =-21.70880995483007-13.76259674006133*
x1;
238 if( teorginal > 1e3 )
240 y *= pow(teorginal/1e3 , 0.33 );
247STATIC double h21_t_lt_20(
double temp )
253 temp =
MAX2( 1., temp );
255 y =9.720710314268267E-08+6.325515312006680E-08*
x1;
270 temp =
MIN2( 300., 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. )
276 x3 =
MIN2( 1000., teorginal );
278 y =-21.70880995483007-13.76259674006133*
x2;
282 if( teorginal > 1e3 )
285 y *= pow(teorginal/1e3 , 0.33 );
296 temp =
MAX2(1., 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);
349 temp =
MAX2( 2. , temp );
350 temp =
MIN2( 2e4 , temp );
358 hold =9.588389834316704E-11 - 5.158891920816405E-14*
x1
359 +5.895348443553458E-19*
x2 + 2.053049602324290E-11*x3
360 +9.122617940315725E-10*x4;
399 long int i, j, mass, nelec, ion, nelem;
411 fprintf(
ioQQQ,
" Hyperfine opening hyperfine.dat:");
413 ioDATA =
open_data(
"hyperfine.dat",
"r" );
418 fprintf(
ioQQQ,
" Hyperfine could not read first line of hyperfine.dat.\n");
424 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL )
428 if( chLine[0] !=
'#')
448 if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
450 fprintf(
ioQQQ,
" Hyperfine could not rewind hyperfine.dat.\n");
457 fprintf(
ioQQQ,
" Hyperfine could not read first line of hyperfine.dat.\n");
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);
469 static int iYR=2, iMN=10, iDY=28;
470 if( ( nelem != iYR ) || ( nelec != iMN ) || ( ion != iDY ) )
473 " Hyperfine: the version of hyperfine.dat in the data directory is not the current version.\n" );
475 " I expected to find the number %i %i %i and got %li %li %li instead.\n" ,
477 nelem , nelec , ion );
478 fprintf(
ioQQQ,
"Here is the line image:\n==%s==\n", chLine );
495 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL )
499 if( chLine[0] !=
'#')
505 sscanf(chLine,
"%li%i%le%i%le%le%le%le%le%le%le%le%le%le%le%le%le%le%le",
529 HFLines[j].Emis().damp() = 1e-20f;
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);
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};
603 if(
phycon.te <= TemperatureTable[0])
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);
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);
631 ASSERT( (TemperatureTable[j-1] <=
phycon.te ) && (TemperatureTable[j] >=
phycon.te) );
638 else if(
phycon.te < TemperatureTable[j])
641 (log10(TemperatureTable[j-1]) - log10(TemperatureTable[j]));
643 upsilon = (log10((
phycon.te)) - (log10(TemperatureTable[j-1])))*slope + log10(
Strengths[i].strengths[j-1]);
644 upsilon = pow(10., upsilon);
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 HyperfineCreate(void)
STATIC double h21_t_ge_10(double temp)
sys_float sexp(sys_float x)
const int INPUT_LINE_LENGTH
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
bool fp_equal(sys_float x, sys_float y, int n=3)
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
NORETURN void TotalInsanity(void)
sys_float SDIV(sys_float x)
#define DEBUG_ENTRY(funcname)
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
t_iso_sp iso_sp[NISO][LIMELM]
double GetGF(double trans_prob, double enercm, double gup)
static realnum * wavelength
UNUSED const double COLL_CONST
TransitionList HFLines("HFLines", &AnonStates)
vector< TransitionList > AllTransitions
double OccupationNumberLine(const TransitionProxy &t)
double TexcLine(const TransitionProxy &t)