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