52 else if( PopTot == 0 )
92 "NOTE Lya masing will invert 21 cm, occupation number set zero\n");
106 double texc, texc1 = 0., texc2 = 0.;
116 enum { DEBUG_SPEC =
false };
119 fprintf(
ioQQQ,
"DEBUG texc %12.3e excitation %12.3e kinetic %12.3e\n",
128 texc1 =
sexp(0.068/texc);
129 texc2 =
sexp(((82259.272-82258.907)*T1CM)/texc);
135 double occnu_lya_23 = occnu_lya;
136 double occnu_lya_13 = occnu_lya*texc1;
137 double occnu_lya_14 = occnu_lya_13*texc2;
138 double occnu_lya_24 = occnu_lya*texc2;
143 double pump21 = pump12 * (*
HFLines[0].Lo()).g() / (*
HFLines[0].Hi()).g();
160 3.*a31*occnu_lya_13 *a32/(a31+a32)+
163 3.*a41*occnu_lya_14 *a42/(a41+a42);
176 occnu_lya_23*a32 * a31/(a31+a32)+
177 occnu_lya_24*a42*a41/(a41+a42);
180 double x = rate12 /
SDIV(a21 + rate21);
184 (*
HFLines[0].Hi()).Pop() = (x/(1.+x))* PopTot;
185 (*
HFLines[0].Lo()).Pop() = (1./(1.+x))* PopTot;
217 temp =
MIN2(1e4 , temp );
221 return exp10( -9.607 + log10( sqrt(temp)) *
sexp(
powpq(log10(temp),9,2) / 1800. ));
230 STATIC double h21_t_ge_20(
double temp )
236 temp =
MIN2( 1000.,temp );
238 y =-21.70880995483007-13.76259674006133*
x1;
244 if( teorginal > 1e3 )
246 y *=
pow(teorginal/1e3 , 0.33 );
253 STATIC double h21_t_lt_20(
double temp )
259 temp =
MAX2( 1., temp );
261 y =9.720710314268267E-08+6.325515312006680E-08*
x1;
272 double teorginal = temp;
275 temp =
MIN2( 300., temp );
277 double y = 1.4341127e-9
278 + 9.4161077e-15 * temp
279 - 9.2998995e-9 / log(temp)
280 + 6.9539411e-9 / sqrt(temp)
281 + 1.7742293e-8 * log(temp)/
pow2(temp);
282 if( teorginal > 300. )
285 temp =
MIN2( 1000., teorginal );
286 y = -21.70880995483007 - 13.76259674006133 / sqrt(temp);
289 if( teorginal > 1e3 )
292 y *=
pow( teorginal/1e3 , 0.33 );
300 temp =
MAX2(1., temp );
302 + 2.331358e-11 * temp
303 + 9.5640586e-11 *
pow2(log(temp))
304 - 4.6220869e-10 * sqrt(temp)
305 - 4.1719545e-10 / sqrt(temp);
354 temp =
MAX2( 2. , temp );
355 temp =
MIN2( 2e4 , temp );
358 return 9.588389834316704E-11
359 - 5.158891920816405E-14 * temp
360 + 5.895348443553458E-19 * temp * temp
361 + 2.053049602324290E-11 * sqrt(temp)
362 + 9.122617940315725E-10 * log(temp) / temp;
416 ioDATA =
open_data(
"hyperfine.dat",
"r" );
421 fprintf(
ioQQQ,
" Hyperfine could not read first line of hyperfine.dat.\n");
428 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL )
430 if( chLine[0] ==
'#' )
434 else if( strstr(chLine,
"TDATA:") == NULL)
437 int Aiso = atoi(data[0].c_str());
438 int nelem = atoi(data[1].c_str());
449 if (data.size() <= 1)
451 fprintf(
ioQQQ,
"HyperfineCreate: Error: Invalid number of temperatures in 'TDATA:': %d\n",
456 Ntemp = data.size() - 1;
457 csTemp = (
double *)
MALLOC( (
size_t)Ntemp*
sizeof(double) );
460 for (std::vector<string>::const_iterator it = data.begin(); it != data.end() && i <=
Ntemp; it++, i++)
464 csTemp[i-1] = atof((*it).c_str());
471 ASSERT(nHFLines > 0 && Ntemp > 0);
472 for(
int i = 0; i <
Ntemp; i++)
477 ASSERT(csTemp[i] > csTemp[i-1]);
497 colstr[j].cs = (
double *)
MALLOC( (
size_t)(
Ntemp)*
sizeof(
double) );
498 colstr[j].cs2d= (
double *)
MALLOC( (
size_t)(
Ntemp)*
sizeof(
double) );
504 if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
506 fprintf(
ioQQQ,
" Hyperfine could not rewind hyperfine.dat.\n");
513 fprintf(
ioQQQ,
" Hyperfine could not read first line of hyperfine.dat.\n");
521 int year = (int)
FFmtRead(chLine,&j,
sizeof(chLine),&lgEOL);
522 int month = (int)
FFmtRead(chLine,&j,
sizeof(chLine),&lgEOL);
523 int day = (int)
FFmtRead(chLine,&j,
sizeof(chLine),&lgEOL);
526 static int iYR=13, iMN=10, iDY=18;
527 if( ( year != iYR ) || ( month != iMN ) || ( day != iDY ) )
530 " Hyperfine: the version of hyperfine.dat in the data directory is not the current version.\n" );
532 " I expected to find the number %i %i %i and got %i %i %i instead.\n" ,
534 year , month , day );
535 fprintf(
ioQQQ,
"Here is the line image:\n==%s==\n", chLine );
553 if( chLine[0] ==
'#' || strstr(chLine,
"TDATA:") != NULL)
557 int Aiso = atoi(data[0].c_str());
558 int nelem = atoi(data[1].c_str());
569 (*
HFLines[j].Hi()).nelem() = nelem;
572 (*
HFLines[j].Hi()).IonStg() = atoi(data[2].c_str());
587 double tmp = (*
HFLines[j].Hi()).g();
592 double fenergyWN =
MAX2(ENERGY_MIN_WN, 1.0 / wavelength);
612 if( data.size() > 6 )
615 for (
int ij = 6, ii = 0; ij < (int) data.size() && ii <
Ntemp; ij++, ii++)
617 colstr[j].cs[ii] = atof(data[ij].c_str());
618 ASSERT(colstr[j].cs[ii] >= 0.0);
622 spline(csTemp, colstr[j].cs, Ntemp, 2e31, 2e31, colstr[j].cs2d);
627 free( colstr[j].cs ); colstr[j].cs = NULL;
628 free( colstr[j].cs2d ); colstr[j].cs2d = NULL;
640 for(
long nelem = 0; nelem <
LIMELM; nelem++ )
658 q21 = COLL_CONST * upsilon / (
phycon.
sqrte * (*HFLines[i].Hi()).g());
660 q12 = (*HFLines[i].Hi()).g()/ (*HFLines[i].Lo()).g() * q21 * exp(-1 * h * c * HFLines[i].EnergyWN / (k *
phycon.
te));
662 x = Ne * q12 / (HFLines[i].Emis().Aul() * (1 + Ne * q21 / HFLines[i].Aul()));
663 HFLines[i].xIntensity() =
N * HFLines[i].Emis().Aul() * x / (1.0 + x) * h * c / (HFLines[i].EnergyAng() / 1e8);
681 if( colstr[i].cs == NULL )
682 return HFLines[i].Coll().col_str();
687 upsilon = colstr[i].cs[0];
689 else if(
phycon.
te >= csTemp[Ntemp-1])
693 double slope = log10(colstr[i].cs[j-1]/colstr[i].cs[j]) / log10(csTemp[j-1]/csTemp[j]);
694 upsilon = log10(
phycon.
te/csTemp[j])*slope + log10(colstr[i].cs[j]);
695 upsilon =
exp10( upsilon);
699 splint( csTemp, colstr[i].cs, colstr[i].cs2d, Ntemp,
phycon.
te, &upsilon );
double HyperfineCS(size_t i)
double TexcLine(const TransitionProxy &t)
STATIC double h21_t_ge_10(double temp)
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
STATIC double h21_t_lt_10(double temp)
NORETURN void TotalInsanity(void)
double H21cm_H_atom(double temp)
void set_xIntensity(const TransitionProxy &t)
TransitionList HFLines("HFLines",&AnonStates)
sys_float sexp(sys_float x)
void HyperfineCreate(void)
void splint(const double xa[], const double ya[], const double y2a[], long int n, double x, double *y)
t_iso_sp iso_sp[NISO][LIMELM]
bool lgLyaMaseCommentDone
const double TEMP_LIMIT_LOW
double xIonDense[LIMELM][LIMELM+1]
double H21cm_proton(double temp)
double H21cm_electron(double temp)
void resize(size_t newsize)
void spline(const double x[], const double y[], long int n, double yp1, double ypn, double y2a[])
realnum getMagMom(int Anum)
realnum getSpin(int Anum)
const double TEMP_LIMIT_HIGH
const int INPUT_LINE_LENGTH
double OccupationNumberLine(const TransitionProxy &t)
LyaSourceFunctionShape LyaSourceFunctionShape_assumed
double GetGF(double trans_prob, double enercm, double gup)
#define DEBUG_ENTRY(funcname)
double powpq(double x, int p, int q)
static vector< realnum > wavelength
const double ENERGY_MIN_WN
int fprintf(const Output &stream, const char *format,...)
sys_float SDIV(sys_float x)
double pow(double x, int i)
void Split(const string &str, const string &sep, vector< string > &lst, split_mode mode)
vector< TransitionList > AllTransitions
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
void MakeCS(const TransitionProxy &t)
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)