00001
00002
00003
00004 #include "cddefines.h"
00005 #include "atmdat.h"
00006 #include "dense.h"
00007 #include "hydrogenic.h"
00008 #include "iso.h"
00009 #include "rfield.h"
00010 #include "geometry.h"
00011 #include "lines.h"
00012 #include "lines_service.h"
00013 #include "phycon.h"
00014 #include "radius.h"
00015 #include "secondaries.h"
00016 #include "taulines.h"
00017 #include "trace.h"
00018
00019 void lines_hydro(void)
00020 {
00021 long ipISO = ipH_LIKE;
00022 long int i, nelem, ipHi, ipLo;
00023 char chLabel[5]=" ";
00024
00025 double hbetab,
00026 em ,
00027 pump ,
00028 caseb;
00029
00030 DEBUG_ENTRY( "lines_hydro()" );
00031
00032 if( trace.lgTrace )
00033 {
00034 fprintf( ioQQQ, " lines_hydro called\n" );
00035 }
00036
00037 i = StuffComment( "H-like iso-sequence" );
00038 linadd( 0., (realnum)i , "####", 'i',
00039 " start H -like iso sequence ");
00040
00041 linadd(MAX2(0.,iso.xLineTotCool[ipH_LIKE][ipHYDROGEN]),912,"Clin",'c',
00042 " total collisional cooling due to all hydrogen lines ");
00043
00044 linadd(MAX2(0.,-iso.xLineTotCool[ipH_LIKE][ipHYDROGEN]),912,"Hlin",'h' ,
00045 " total collisional heating due to all hydrogen lines ");
00046
00047
00048
00049
00050
00051
00052 linadd(MAX2(0.,iso.cLya_cool[ipH_LIKE][ipHYDROGEN]),1216,"Cool",'i',
00053 "collisionally excited La cooling ");
00054
00055 linadd(MAX2(0.,-iso.cLya_cool[ipH_LIKE][ipHYDROGEN]),1216,"Heat",'i',
00056 " collisionally de-excited La heating ");
00057
00058 linadd(MAX2(0.,iso.cLyrest_cool[ipH_LIKE][ipHYDROGEN]),960,"Crst",'i',
00059 " cooling due to n>2 Lyman lines ");
00060
00061 linadd(MAX2(0.,-iso.cLyrest_cool[ipH_LIKE][ipHYDROGEN]),960,"Hrst",'i',
00062 " heating due to n>2 Lyman lines ");
00063
00064 linadd(MAX2(0.,iso.cBal_cool[ipH_LIKE][ipHYDROGEN]),4861,"Crst",'i',
00065 " cooling due to n>3 Balmer lines ");
00066
00067 linadd(MAX2(0.,-iso.cBal_cool[ipH_LIKE][ipHYDROGEN]),4861,"Hrst",'i',
00068 " heating due to n>3 Balmer lines ");
00069
00070 linadd(MAX2(0.,iso.cRest_cool[ipH_LIKE][ipHYDROGEN]),0,"Crst",'i',
00071 " cooling due to higher Paschen lines ");
00072
00073 linadd(MAX2(0.,-iso.cRest_cool[ipH_LIKE][ipHYDROGEN]),0,"Hrst",'i',
00074 " heating due to higher Paschen lines ");
00075
00076
00077 secondaries.SecHIonMax = MAX2( secondaries.SecHIonMax , secondaries.sec2total );
00078
00079
00080 atmdat.HIonFracMax = MAX2( atmdat.HIonFracMax, atmdat.HIonFrac);
00081
00082
00083 hydro.HCollIonMax =
00084 (realnum)MAX2( hydro.HCollIonMax , hydro.H_ion_frac_collis );
00085
00086 linadd(secondaries.x12tot*dense.xIonDense[ipHYDROGEN][1]*StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop*1.634e-11,1216,"LA X" ,'i',
00087 "Lyaa contribution from suprathermal secondaries from ground ");
00088
00089
00090
00091
00092 pump = (double)(Transitions[ipH_LIKE][ipHYDROGEN][ipH4p][ipH1s].Emis->pump*StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop*
00093 dense.xIonDense[ipHYDROGEN][1]*4.09e-12*0.4836);
00094 linadd(pump,4861,"Pump",'r',
00095 "part of Hbeta formed by continuum pumping");
00096
00097 linadd(MAX2(0.,iso.coll_ion[ipH_LIKE][ipHYDROGEN]),0,"CION",'c',
00098 "collision ionization cooling of hydrogen ");
00099
00100 linadd(MAX2(-iso.coll_ion[ipH_LIKE][ipHYDROGEN],0.),0,"3bHt",'h',
00101 " this is the heating due to 3-body recombination ");
00102
00103 fixit();
00104 linadd(dense.xIonDense[ipHYDROGEN][1]*StatesElem[ipH_LIKE][ipHYDROGEN][ipH2p].Pop*0.*iso.pestrk[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s]*1.634e-11,1216,"Strk",'i',
00105 " Stark broadening contribution to line ");
00106
00107 linadd(dense.xIonDense[ipHYDROGEN][1]*StatesElem[ipH_LIKE][ipHYDROGEN][ipH3s].Pop*iso.pestrk[ipH_LIKE][ipHYDROGEN][ipH3s][ipH2p]*3.025e-12,
00108 6563,"Strk",'i',
00109 " Stark broadening contribution to line ");
00110
00111 linadd(dense.xIonDense[ipHYDROGEN][1]*StatesElem[ipH_LIKE][ipHYDROGEN][ipH4s].Pop*iso.pestrk[ipH_LIKE][ipHYDROGEN][ipH4s][ipH2p]*4.084e-12,
00112 4861,"Strk",'i',
00113 "Stark broadening contribution to line ");
00114
00115 linadd(dense.xIonDense[ipHYDROGEN][1]*StatesElem[ipH_LIKE][ipHYDROGEN][ipH4p].Pop*iso.pestrk[ipH_LIKE][ipHYDROGEN][ipH4p][ipH3s]*1.059e-12,
00116 18751,"Strk",'i',
00117 " Stark broadening contribution to line ");
00118
00119
00120
00121
00122
00123 if( iso.numLevels_max[ipH_LIKE][ipHYDROGEN] > 5 )
00124 linadd(dense.xIonDense[ipHYDROGEN][1]*StatesElem[ipH_LIKE][ipHYDROGEN][ipH5p].Pop*iso.pestrk[ipH_LIKE][ipHYDROGEN][ipH5p][ipH4s]*4.900e-13,40512,"Strk",'i',
00125 "Stark broadening part of line");
00126
00127
00128
00129 ASSERT( LineSave.ipass <1 ||
00130 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->ots>= 0.);
00131
00132 linadd(Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->ots*Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].EnergyErg, 1216,"Dest",'i',
00133 " portion of line lost due to absorp by background opacity ");
00134
00135
00136
00137
00138 if( iso.numLevels_max[ipH_LIKE][ipHYDROGEN] > ipH3p )
00139 linadd(Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH2s].Emis->ots*Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH2s].EnergyErg, 6563,"Dest",'i',
00140 "Ha destroyed by background opacity");
00141
00142
00143
00144
00145 if( iso.numLevels_max[ipH_LIKE][ipHYDROGEN] > ipH5p )
00146 linadd(Transitions[ipH_LIKE][ipHYDROGEN][ipH5p][ipH4s].Emis->ots*Transitions[ipH_LIKE][ipHYDROGEN][ipH5p][ipH4s].EnergyErg,40516, "Dest",'i',
00147 "portion of line lost due to absorb by background opacity");
00148
00149
00150
00151
00152 if( iso.numLevels_max[ipH_LIKE][ipHYDROGEN] > ipH4p )
00153 linadd(Transitions[ipH_LIKE][ipHYDROGEN][ipH4p][ipH2s].Emis->ots*Transitions[ipH_LIKE][ipHYDROGEN][ipH4p][ipH2s].EnergyErg, 4861,"Dest",'i',
00154 "portion of line lost due to absorb by background opacity");
00155
00156
00157
00158
00159 if( iso.numLevels_max[ipH_LIKE][ipHYDROGEN] > ipH4p )
00160 linadd(Transitions[ipH_LIKE][ipHYDROGEN][ipH4p][ipH3s].Emis->ots*Transitions[ipH_LIKE][ipHYDROGEN][ipH4p][ipH3s].EnergyErg ,18751, "Dest",'i',
00161 "portion of line lost due to absorb by background opacity");
00162
00163 linadd(StatesElem[ipH_LIKE][ipHYDROGEN][ipH2p].Pop*dense.xIonDense[ipHYDROGEN][1]*Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Aul*
00164 hydro.dstfe2lya*Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].EnergyErg , 1216 , "Fe 2" , 'i',
00165 "Ly-alpha destroyed by overlap with FeII " );
00166
00167 linadd(iso.RadRec_caseB[ipH_LIKE][ipHYDROGEN]*dense.xIonDense[ipHYDROGEN][1]*dense.eden * 1.64e-11,1216,"Ca B",'i',
00168 " simple high-density case b intensity of Ly-alpha, no two photon ");
00169
00170
00171 if( nzone == 1 )
00172 {
00173
00174 caseb = rfield.qhtot*
00175 atmdat_HS_caseB( 4 , 2 , 1 , phycon.te , dense.eden, 'b' ) / iso.RadRec_caseB[ipH_LIKE][ipHYDROGEN];
00176
00177
00178 if( caseb < 0 )
00179 {
00180 caseb = rfield.qhtot*4.75e-13;
00181 }
00182 LineSv[LineSave.nsum].sumlin[LineSave.lgLineEmergent] = 0.;
00183 }
00184 else
00185 {
00186 caseb = 0.;
00187 }
00188
00189
00190 linadd( caseb/radius.dVeff*geometry.covgeo , 4861 , "Q(H)" , 'i' ,
00191 "Case B H-beta computed from Q(H) and specified covering factor");
00192
00193 if( nzone == 1 )
00194 {
00195 caseb = rfield.qhtot*Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].EnergyErg;
00196 LineSv[LineSave.nsum].sumlin[LineSave.lgLineEmergent] = 0.;
00197 }
00198 else
00199 {
00200 caseb = 0.;
00201 }
00202
00203 linadd( caseb/radius.dVeff*geometry.covgeo , 1216 , "Q(H)" , 'i',
00204 "Ly-alpha from Q(H), high-dens lim, specified covering factor" );
00205
00206
00207 for( nelem=ipISO; nelem < LIMELM; nelem++ )
00208 {
00209 if( dense.lgElmtOn[nelem] )
00210 {
00211 for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
00212 {
00213 for( ipLo=0; ipLo < ipHi; ipLo++ )
00214 {
00215 if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA )
00216 continue;
00217
00218
00219 Transitions[ipISO][nelem][ipHi][ipLo].Emis->phots =
00220 Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul*
00221 StatesElem[ipISO][nelem][ipHi].Pop*
00222 Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pesc*
00223 dense.xIonDense[nelem][nelem+1-ipISO];
00224
00225
00226 Transitions[ipISO][nelem][ipHi][ipLo].Emis->xIntensity =
00227 Transitions[ipISO][nelem][ipHi][ipLo].Emis->phots*
00228 Transitions[ipISO][nelem][ipHi][ipLo].EnergyErg;
00229 }
00230 }
00231 }
00232 }
00233
00234
00235
00236 for( nelem=0; nelem < LIMELM; nelem++ )
00237 {
00238 if( dense.IonHigh[nelem] == nelem + 1 )
00239 {
00240
00241 for( ipHi=3; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
00242 {
00243 long index_of_nHi_P;
00244
00245
00246 if( N_(ipHi) > iso.n_HighestResolved_max[ipH_LIKE][nelem] )
00247 index_of_nHi_P = ipHi;
00248 else
00249 index_of_nHi_P = iso.QuantumNumbers2Index[ipH_LIKE][nelem][ N_(ipHi) ][1][2];
00250
00251
00252 for( ipLo=0; ipLo < ipHi; ipLo++ )
00253 {
00254 long index_of_nLo_S = iso.QuantumNumbers2Index[ipH_LIKE][nelem][ N_(ipLo) ][0][2];
00255
00256
00257
00258 if( N_(ipLo) > iso.n_HighestResolved_local[ipH_LIKE][nelem] || N_(ipLo) == N_(ipHi) )
00259 break;
00260
00261 if( Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA )
00262 continue;
00263
00264
00265 if( ipHi == index_of_nHi_P && ipLo == index_of_nLo_S )
00266 continue;
00267 else
00268 {
00269
00270 Transitions[ipH_LIKE][nelem][index_of_nHi_P][index_of_nLo_S].Emis->xIntensity +=
00271 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->xIntensity;
00272
00273 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->xIntensity = 0.;
00274
00275 }
00276 }
00277 }
00278 }
00279 }
00280
00281
00282 hbetab = (double)((pow(10.,-20.89 - 0.10612*POW2(phycon.alogte - 4.4)))/
00283 phycon.te);
00284
00285
00286
00287 ASSERT( iso.numLevels_max[ipH_LIKE][ipHYDROGEN] > 4 );
00288 hbetab *= dense.xIonDense[ipHYDROGEN][1]*dense.eden;
00289
00290 lindst(hbetab, -4861 ,"CaBo",
00291 1 ,'i',false,
00292 " this is old case b based on Ferland (1980) PASP ");
00293
00294 if( dense.lgElmtOn[ipHELIUM] )
00295 {
00296
00297
00298
00299 ASSERT( iso.numLevels_max[ipH_LIKE][ipHELIUM] > 4 );
00300
00301 em = 2.03e-20/(phycon.te70*phycon.te10*phycon.te03);
00302 em *= dense.xIonDense[ipHELIUM][2]*dense.eden;
00303
00304 lindst(em,-1640,"CaBo",
00305 1,'i',false,
00306 " old prediction of He II 1640, Case B at low densities");
00307
00308
00309
00310 em = 2.52e-20/(pow(phycon.te,1.05881));
00311 em *= dense.xIonDense[ipHELIUM][2]*dense.eden;
00312
00313 lindst(em,-4686,"CaBo", 1,'i',false,
00314 " old prediction of He II 4686, Case B at low densities");
00315 }
00316
00317
00318
00319 if( LineSave.ipass <= 0 )
00320 {
00321 for(nelem=0; nelem<HS_NZ; ++nelem )
00322 {
00323 atmdat.lgHCaseBOK[0][nelem] = true;
00324 atmdat.lgHCaseBOK[1][nelem] = true;
00325 }
00326 }
00327
00328 for( nelem=0; nelem < LIMELM; nelem++ )
00329 {
00330 if( dense.lgElmtOn[nelem] )
00331 {
00332
00333
00334
00335 if( nelem < HS_NZ && (nelem<2 || nelem>4) )
00336 {
00337 int iCase;
00338 for( iCase=0; iCase<2; ++iCase )
00339 {
00340 char chAB[2]={'A','B'};
00341 char chLab[5]="Ca ";
00342
00343
00344
00345
00346
00347 for( ipLo=1+iCase; ipLo<MIN2(10,iso.n_HighestResolved_max[ipH_LIKE][nelem] + iso.nCollapsed_max[ipH_LIKE][nelem]); ++ipLo )
00348 {
00349 for( ipHi=ipLo+1; ipHi< MIN2(25,iso.n_HighestResolved_max[ipH_LIKE][nelem] + iso.nCollapsed_max[ipH_LIKE][nelem]+1); ++ipHi )
00350 {
00351 realnum wl;
00352 double case_b_Intensity;
00353 long int ipCHi , ipCLo;
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364 case_b_Intensity = atmdat_HS_caseB( ipHi,ipLo , nelem+1, phycon.te , dense.eden, chAB[iCase] );
00365 if( case_b_Intensity<=0. )
00366 {
00367 atmdat.lgHCaseBOK[iCase][nelem] = false;
00368 case_b_Intensity = 0.;
00369 }
00370
00371 case_b_Intensity *= dense.xIonDense[nelem][nelem+1-ipISO]*dense.eden;
00372
00373 if( iCase==0 && ipLo==1 )
00374 {
00375
00376 ipCHi = ipHi;
00377 ipCLo = 0;
00378 }
00379 else
00380 {
00381
00382 ipCHi = ipHi;
00383 ipCLo = ipLo;
00384 }
00385
00386
00387 chLab[3] = chAB[iCase];
00388
00389
00390 if( ipCHi > 2 )
00391 {
00392 if( ipCLo > 2 )
00393 {
00394
00395 ipCHi = iso.QuantumNumbers2Index[ipH_LIKE][nelem][ipCHi][1][2];
00396 ipCLo = iso.QuantumNumbers2Index[ipH_LIKE][nelem][ipCLo][0][2];
00397 }
00398 else if( ipCLo == 2 )
00399 {
00400
00401 ipCHi = iso.QuantumNumbers2Index[ipH_LIKE][nelem][ipCHi][0][2];
00402 }
00403 else if( ipCLo == 1 || ipCLo == 0 )
00404 {
00405
00406 ipCHi = iso.QuantumNumbers2Index[ipH_LIKE][nelem][ipCHi][1][2];
00407 }
00408 }
00409
00410
00411 wl = Transitions[ipH_LIKE][nelem][ipCHi][ipCLo].WLAng;
00412 atmdat.WaveLengthCaseB[nelem][ipHi][ipLo] = wl;
00413
00414 lindst(case_b_Intensity,wl,chLab,Transitions[ipH_LIKE][nelem][ipCHi][ipCLo].ipCont,'i',false,
00415 " case a or case b from Hummer & Storey tables" );
00416 }
00417 }
00418 }
00419 }
00420
00421
00422
00423
00424
00425 if( LineSave.ipass == 0 )
00426 {
00427
00428
00429
00430 chIonLbl(chLabel, &Transitions[ipH_LIKE][nelem][ipH2s][ipH1s]);
00431 }
00432 linadd( Transitions[ipH_LIKE][nelem][ipH2s][ipH1s].Emis->xIntensity , 0,chLabel,'r',
00433 "two-photon emission");
00434
00435 linadd(
00436 StatesElem[ipH_LIKE][nelem][ipH2s].Pop*
00437 dense.xIonDense[nelem][nelem+1-ipISO]*
00438 iso.TwoNu_induc_dn[ipH_LIKE][nelem]*
00439 Transitions[ipH_LIKE][nelem][ipH2s][ipH1s].EnergyErg,
00440 22, chLabel ,'i',
00441 "induced two-photon emission ");
00442
00443
00444
00445 for( ipLo=ipH1s; ipLo < iso.numLevels_max[ipH_LIKE][nelem]-1; ipLo++ )
00446 {
00447
00448 if( ipLo==ipH2p )
00449 continue;
00450
00451 for( ipHi=ipLo+1; ipHi < iso.numPrintLevels[ipH_LIKE][nelem]; ipHi++ )
00452 {
00453
00454
00455 if( Transitions[ipH_LIKE][nelem][ipHi][ipLo].ipCont < 1 )
00456 continue;
00457
00458 #if 0
00459 char chCommentTemp[23];
00460 strcpy( chCommentTemp, StatesElem[ipH_LIKE][nelem][ipLo].chConfig );
00461 strcat( chCommentTemp, " - " );
00462 strcat( chCommentTemp, StatesElem[ipH_LIKE][nelem][ipHi].chConfig );
00463 const char* chComment = chCommentTemp;
00464 PutLine(&Transitions[ipH_LIKE][nelem][ipHi][ipLo], chComment );
00465 #else
00466 PutLine(&Transitions[ipH_LIKE][nelem][ipHi][ipLo],
00467 "predicted line, all processes included");
00468 #endif
00469
00470 }
00471 }
00472 }
00473 }
00474
00475 if( trace.lgTrace )
00476 {
00477 fprintf( ioQQQ, " lines_helium returns\n" );
00478 }
00479 return;
00480 }