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