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 }