00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "cddefines.h"
00011 #include "lines_service.h"
00012 #include "path.h"
00013 #include "phycon.h"
00014 #include "dense.h"
00015 #include "rfield.h"
00016 #include "taulines.h"
00017 #include "iso.h"
00018 #include "trace.h"
00019 #include "hyperfine.h"
00020 #include "physconst.h"
00021
00022
00023
00024 void H21_cm_pops( void )
00025 {
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 double x,
00043 PopTot;
00044 double a32,a31,a41,a42,a21, occnu_lya ,
00045 rate12 , rate21 , pump12 , pump21 , coll12 , coll21,
00046 texc , occnu_lya_23 , occnu_lya_13,occnu_lya_24,occnu_lya_14, texc1, texc2;
00047 a31 = 2.08e8;
00048 a32 = 4.16e8;
00049 a41 = 4.16e8;
00050 a42 = 2.08e8;
00051
00052
00053
00054
00055 a21 = 2.85e-15;
00056
00057
00058
00059 a21 *= (HFLines[0].Pdest + HFLines[0].Pesc + HFLines[0].Pelec_esc);
00060
00061
00062
00063 occnu_lya = OccupationNumberLine( &EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s] ) *
00064 hyperfine.lgLya_pump_21cm;
00065
00066
00067
00068 texc = TexcLine( &EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s] );
00069
00070 if( texc > 0. )
00071 {
00072
00073 texc1 = sexp(0.068/texc);
00074 texc2 = sexp(((82259.272-82258.907)*T1CM)/texc);
00075 }
00076 else
00077 {
00078 texc1 = 0.;
00079 texc2 = 0.;
00080 }
00081
00082
00083
00084
00085 occnu_lya_23 = occnu_lya;
00086 occnu_lya_13 = occnu_lya*texc1;
00087 occnu_lya_14 = occnu_lya_13*texc2;
00088 occnu_lya_24 = occnu_lya*texc2;
00089
00090
00091
00092 pump12 = HFLines[0].pump;
00093 pump21 = pump12 * HFLines[0].gLo / HFLines[0].gHi;
00094
00095
00096
00097
00098 coll12 = HFLines[0].cs*dense.cdsqte/HFLines[0].gLo*rfield.ContBoltz[HFLines[0].ipCont-1];
00099 coll21 = HFLines[0].cs*dense.cdsqte/HFLines[0].gHi;
00100
00101
00102
00103 rate12 =
00104
00105 coll12 +
00106
00107 pump12 +
00108
00109 3.*a31*occnu_lya_13 *a32/(a31+a32)+
00110
00111
00112 3.*a41*occnu_lya_14 *a42/(a41+a42);
00113
00114
00115
00116
00117
00118 rate21 =
00119
00120 coll21 +
00121
00122 pump21 +
00123
00124
00125 occnu_lya_23*a32 * a31/(a31+a32)+
00126 occnu_lya_24*a42*a41/(a41+a42);
00127
00128
00129 x = rate12 / SDIV(a21 + rate21);
00130
00131 PopTot = iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s]*dense.xIonDense[ipHYDROGEN][1];
00132
00133 HFLines[0].PopHi = (x/(1.+x))* PopTot;
00134 HFLines[0].PopLo = (1./(1.+x))* PopTot;
00135
00136
00137
00138
00139
00140
00141 HFLines[0].PopOpc = HFLines[0].PopLo*((3*rate21- rate12) + 3*a21)/SDIV(3*(a21+ rate21));
00142
00143
00144
00145
00146
00147
00148
00149 if( HFLines[0].PopHi > SMALLFLOAT )
00150 {
00151 hyperfine.Tspin21cm = TexcLine( &HFLines[0] );
00152
00153
00154 if( hyperfine.Tspin21cm == 0. )
00155 hyperfine.Tspin21cm = phycon.te;
00156 }
00157 else
00158 {
00159 hyperfine.Tspin21cm = phycon.te;
00160 }
00161
00162 return;
00163 }
00164
00165
00166
00167 double H21cm_electron( double temp )
00168 {
00169 double hold;
00170 temp = MIN2(1e4 , temp );
00171
00172
00173
00174 hold = -9.607 + log10( sqrt(temp)) * sexp( pow(log10(temp) , 4.5 ) / 1800. );
00175 hold = pow(10.,hold );
00176 return( hold );
00177 }
00178
00179
00180
00181
00182
00183
00184 #if 0
00185 static double h21_t_ge_20( double temp )
00186 {
00187 double y;
00188 double x1,
00189 teorginal = temp;
00190
00191 temp = MIN2( 1000.,temp );
00192 x1 =1.0/sqrt(temp);
00193 y =-21.70880995483007-13.76259674006133*x1;
00194 y = exp(y);
00195
00196
00197
00198
00199 if( teorginal > 1e3 )
00200 {
00201 y *= pow(teorginal/1e3 , 0.33 );
00202 }
00203
00204 return( y );
00205 }
00206
00207
00208 static double h21_t_lt_20( double temp )
00209 {
00210 double y;
00211 double x1;
00212
00213
00214 temp = MAX2( 1., temp );
00215 x1 =temp*log(temp);
00216 y =9.720710314268267E-08+6.325515312006680E-08*x1;
00217 return(y*y);
00218 }
00219 #endif
00220
00221
00222
00223
00224
00225 static double h21_t_ge_10( double temp )
00226 {
00227 double y;
00228 double x1,x2,x3,
00229 teorginal = temp;
00230
00231 temp = MIN2( 300., temp );
00232 x1 =temp;
00233 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);
00234 if( teorginal > 300. )
00235 {
00236
00237 x3 = MIN2( 1000., teorginal );
00238 x2 =1.0/sqrt(x3);
00239 y =-21.70880995483007-13.76259674006133*x2;
00240 y = 1.236686*exp(y);
00241
00242 }
00243 if( teorginal > 1e3 )
00244 {
00245
00246 y *= pow(teorginal/1e3 , 0.33 );
00247 }
00248 return( y );
00249 }
00250
00251 static double h21_t_lt_10( double temp )
00252 {
00253 double y;
00254 double x1;
00255
00256
00257 temp = MAX2(1., temp );
00258 x1 =temp;
00259 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);
00260 return(y);
00261 }
00262
00263 #if 0
00264 double H21cm_H_atom( double temp )
00265 {
00266 double hold;
00267 if( temp >= 20. )
00268 {
00269 hold = h21_t_ge_20( temp );
00270 }
00271 else
00272 {
00273 hold = h21_t_lt_20( temp );
00274 }
00275
00276 return hold;
00277 }
00278 #endif
00279
00280 double H21cm_H_atom( double temp )
00281 {
00282 double hold;
00283 if( temp >= 10. )
00284 {
00285 hold = h21_t_ge_10( temp );
00286 }
00287 else
00288 {
00289 hold = h21_t_lt_10( temp );
00290 }
00291
00292 return hold;
00293 }
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316 typedef struct
00317 {
00318 double strengths[12];
00319 } Ion;
00320
00321 static Ion *Strengths;
00322
00323
00324 void HyperfineCreate(void)
00325 {
00326 FILE *ioDATA;
00327 char chFilename[FILENAME_PATH_LENGTH_2] ,chLine[INPUT_LINE_LENGTH];
00328 bool lgEOL;
00329
00330 float spin, wavelength;
00331 long int i, j, mass, nelec, ion, nelem;
00332
00333 DEBUG_ENTRY( "HyperfineCreate()" );
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345 if( lgDataPathSet == true )
00346 {
00347
00348 strcpy( chFilename , chDataPath );
00349 strcat( chFilename , "hyperfine.dat" );
00350 }
00351 else
00352 {
00353
00354 strcpy( chFilename , "hyperfine.dat" );
00355 }
00356
00357 if( trace.lgTrace )
00358 fprintf( ioQQQ," Hyperfine opening hyperfine.dat:");
00359
00360 if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
00361 {
00362 fprintf( ioQQQ, " Hyperfine could not open hyperfine.dat.\n" );
00363 if( lgDataPathSet == true )
00364 fprintf( ioQQQ, " even tried path\n" );
00365
00366 if( lgDataPathSet == true )
00367 {
00368 fprintf( ioQQQ, " Hyperfine could not open hyperfine.dat.\n");
00369 fprintf( ioQQQ, " path is ==%s==.\n",chDataPath );
00370 fprintf( ioQQQ, " final path is ==%s==.\n",chFilename );
00371 }
00372
00373
00374 path_not_set();
00375
00376 puts( "[Stop in Hyperfine]" );
00377 cdEXIT(EXIT_FAILURE);
00378 }
00379
00380
00381 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00382 {
00383 fprintf( ioQQQ, " Hyperfine could not read first line of hyperfine.dat.\n");
00384 puts( "[Stop in Hyperfine]" );
00385 cdEXIT(EXIT_FAILURE);
00386 }
00387
00388
00389 nHFLines = 0;
00390 while( fgets( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
00391 {
00392
00393
00394 if( chLine[0] != '#')
00395 ++nHFLines;
00396 }
00397
00398
00399 HFLines = (EmLine *)MALLOC( (size_t)(nHFLines)*sizeof(EmLine) );
00400
00401 for( i=0; i< nHFLines; ++i )
00402 {
00403 EmLineJunk( &HFLines[i] );
00404 }
00405
00406 Strengths = (Ion *)MALLOC( (size_t)(nHFLines)*sizeof(Ion) );
00407 hyperfine.HFLabundance = (float *)MALLOC( (size_t)(nHFLines)*sizeof(float) );
00408
00409
00410 if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
00411 {
00412 fprintf( ioQQQ, " Hyperfine could not rewind hyperfine.dat.\n");
00413 puts( "[Stop in Hyperfine]" );
00414 cdEXIT(EXIT_FAILURE);
00415 }
00416
00417
00418 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00419 {
00420 fprintf( ioQQQ, " Hyperfine could not read first line of hyperfine.dat.\n");
00421 puts( "[Stop in Hyperfine]" );
00422 cdEXIT(EXIT_FAILURE);
00423 }
00424
00425
00426 i = 1;
00427 nelem = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00428 nelec = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00429 ion = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00430
00431 {
00432
00433 static int iYR=2, iMN=10, iDY=28;
00434 if( ( nelem != iYR ) || ( nelec != iMN ) || ( ion != iDY ) )
00435 {
00436 fprintf( ioQQQ,
00437 " Hyperfine: the version of hyperfine.dat in the data directory is not the current version.\n" );
00438 fprintf( ioQQQ,
00439 " I expected to find the number %i %i %i and got %li %li %li instead.\n" ,
00440 iYR, iMN , iDY ,
00441 nelem , nelec , ion );
00442 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
00443 puts( "[Stop in Hyperfine]" );
00444 cdEXIT(EXIT_FAILURE);
00445 }
00446 }
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458 j = 0;
00459
00460 while( fgets( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
00461 {
00462
00463
00464 if( chLine[0] != '#')
00465 {
00466 double Aul;
00467 ASSERT( j < nHFLines );
00468
00469
00470 sscanf(chLine, "%li%i%e%i%e%e%le%le%le%le%le%le%le%le%le%le%le%le%le",
00471 &mass,
00472 &HFLines[j].nelem,
00473 &spin,
00474 &HFLines[j].IonStg,
00475 &hyperfine.HFLabundance[j],
00476 &wavelength,
00477 &Aul,
00478 &Strengths[j].strengths[0],
00479 &Strengths[j].strengths[1],
00480 &Strengths[j].strengths[2],
00481 &Strengths[j].strengths[3],
00482 &Strengths[j].strengths[4],
00483 &Strengths[j].strengths[5],
00484 &Strengths[j].strengths[6],
00485 &Strengths[j].strengths[7],
00486 &Strengths[j].strengths[8],
00487 &Strengths[j].strengths[9],
00488 &Strengths[j].strengths[10],
00489 &Strengths[j].strengths[11]);
00490 HFLines[j].Aul = (float)Aul;
00491 HFLines[j].gHi = (float) (2*(spin + .5) + 1);
00492 HFLines[j].gLo = (float) (2*(spin - .5) + 1);
00493 HFLines[j].WLAng = wavelength * 1e8f;
00494 HFLines[j].EnergyWN = (float) (1. / wavelength);
00495 HFLines[j].gf = (float)(GetGF(HFLines[j].Aul, HFLines[j].EnergyWN, HFLines[j].gHi));
00496
00497
00498
00499 ASSERT(HFLines[j].gf > 0.);
00500
00501 j++;
00502 }
00503
00504 }
00505
00506
00507 fclose(ioDATA);
00508
00509 # if 0
00510
00511
00512 for(i = 0; i < nHFLines; i++)
00513 {
00514 N = dense.xIonDense[HFLines[i].nelem-1][HFLines[i].IonStg-1];
00515 Ne = dense.eden;
00516
00517 h = 6.626076e-27;
00518 c = 3e10;
00519 k = 1.380658e-16;
00520
00521 upsilon = HyperfineCS(i);
00522
00523 q21 = COLL_CONST * upsilon / (phycon.sqrte * HFLines[i].gHi);
00524
00525 q12 = HFLines[i].gHi/ HFLines[i].gLo * q21 * exp(-1 * h * c * HFLines[i].EnergyWN / (k * phycon.te));
00526
00527 x = Ne * q12 / (HFLines[i].Aul * (1 + Ne * q21 / HFLines[i].Aul));
00528 HFLines[i].xIntensity = N * HFLines[i].Aul * x / (1.0 + x) * h * c / (HFLines[i].WLAng / 1e8);
00529
00530 }
00531 # endif
00532
00533 DEBUG_EXIT( "HyperfineCreate()" );
00534
00535 return;
00536 }
00537
00538
00539
00540 double HyperfineCS( long i )
00541 {
00542
00543
00544 int j = 0;
00545 # define N_TE_TABLE 12
00546 double slope, upsilon, TemperatureTable[N_TE_TABLE] = {.1e6, .15e6, .25e6, .4e6, .6e6,
00547 1.0e6, 1.5e6, 2.5e6, 4e6, 6e6, 10e6, 15e6};
00548
00549 DEBUG_ENTRY( "HyperfineCS()" );
00550
00551
00552
00553 ASSERT( i >= 0. && i <= nHFLines );
00554
00555
00556
00557
00558
00559
00560
00561
00562 if( phycon.te <= TemperatureTable[0])
00563 {
00564
00565 upsilon = Strengths[i].strengths[0];
00566 }
00567 else if( phycon.te >= TemperatureTable[N_TE_TABLE-1])
00568 {
00569
00570 upsilon = Strengths[i].strengths[N_TE_TABLE-1];
00571 }
00572 else
00573 {
00574 j = 1;
00575
00576
00577
00578 while ( j < N_TE_TABLE && TemperatureTable[j] < phycon.te )
00579 j++;
00580
00581 ASSERT( j >= 0 && j < N_TE_TABLE);
00582 ASSERT( (TemperatureTable[j-1] <= phycon.te ) && (TemperatureTable[j] >= phycon.te) );
00583
00584
00585 if( phycon.te == TemperatureTable[j])
00586
00587 {
00588 upsilon = Strengths[i].strengths[j];
00589 }
00590
00591 else if( phycon.te < TemperatureTable[j])
00592 {
00593 slope = (log10(Strengths[i].strengths[j-1]) - log10(Strengths[i].strengths[j])) /
00594 (log10(TemperatureTable[j-1]) - log10(TemperatureTable[j]));
00595
00596 upsilon = (log10((phycon.te)) - (log10(TemperatureTable[j-1])-6.))*slope + log10(Strengths[i].strengths[j-1]);
00597 upsilon = pow(10., upsilon);
00598 }
00599 else
00600 {
00601 upsilon = Strengths[i].strengths[j-1];
00602 }
00603 }
00604
00605 DEBUG_EXIT( "HyperfineCS()" );
00606
00607 return upsilon;
00608 }