/*Hydrogenic main routine to call HydroLevel and determine model hydrogen atom level balance */ #include "cddefines.h" #include "taulines.h" #include "hydrogenic.h" #include "hphoto.h" #include "trace.h" #include "ionfracs.h" #include "secondaries.h" #include "h.h" #include "phycon.h" #include "photn2.h" #include "hmi.h" #include "hionfraccom.h" #include "hcolrec.h" #include "colion.h" #include "elmton.h" #include "ionrange.h" #include "hcolrt.h" #include "hcaseb.h" #include "printe82.h" #include "hmole.h" void Hydrogenic(void) { long int i , ipLo , ipHi , ipZ; double coltot, gamtot, grddst , photex; float hmovh1 = 0.; # ifdef DEBUG_FUN fputs( "<+>Hydrogenic()\n", debug_fp ); # endif for( ipZ=0; ipZ < LIMELM; ipZ++ ) { /* do not consider elements that have been turned off */ if( elmton.lgElmtOn[ipZ] ) { /* note that ipZ scale is totally on c not physical scale, so 0 is h */ /* evaluate hydrogenic balance if ionization reaches this high */ if( (IonRange.IonHigh[ipZ] == ipZ + 2) ) { /* evaluate collisional rates */ HydroCollid(ipZ); /* evaluate photoionization rates */ HydroPhoto(ipZ); /* evaluate recombination rates */ HydroRecom(ipZ); /* form radiative sums needed for HydroLevel. escape and destruction * probs were done in RTMake */ HydroSumTrans(ipZ); /* solve for the level populations */ HydroLevel(ipZ); if( trace.lgTrace ) { fprintf( ioQQQ, " Hydrogenic Z=%2ld H2OVH1=",ipZ); fprintf( ioQQQ,PrintEfmt("%10.3e", HIonFraccom.HIonFrac[ipZ])); fprintf( ioQQQ, " simple="); fprintf( ioQQQ,PrintEfmt("%10.3e", HIonFraccom.HIonSimple[ipZ])); fprintf( ioQQQ, " H-ovH1:"); fprintf( ioQQQ,PrintEfmt("%10.2e", hmovh1 ) ); fprintf( ioQQQ, "\n"); } } else { HIonFraccom.HIonFrac[ipZ] = 0.; HIonFraccom.HIonSimple[ipZ] = 0.; /* zero it out since no population*/ for( ipHi=IP2S; ipHi <= nhlevel; ipHi++ ) { for( ipLo=IP1S; ipLo < ipHi; ipLo++ ) { /* population of lower level rel to ion */ HydroLines[ipZ][ipHi][ipLo].PopLo = 0.; /* population of upper level rel to ion */ HydroLines[ipZ][ipHi][ipLo].PopHi = 0.; /* population of lower level rel to ion, corrected for stim em */ HydroLines[ipZ][ipHi][ipLo].PopOpc = 0.; } } } } } /* rest is for hydrogen only */ /* do molecular balance * hmovh1 will be ratio of molecular to atomic hydrogen * HIonFrac is fration of H that is ionized, ratio of ion to atom */ hmole(&hmovh1,HIonFraccom.HIonFrac[0]); /* get ionization balance */ h.hii = (float)(phycon.hden*HIonFraccom.HIonFrac[0]/ (1. + HIonFraccom.HIonFrac[0] + hmovh1)); /* save it into the master ionization balance array */ xIonFracs[2][0] = h.hii; h.hi = (float)(phycon.hden/(1. + HIonFraccom.HIonFrac[0] + hmovh1)); xIonFracs[1][0] = h.hi; /* is hydrogen atom in condition of collisional rec? if so then * only very small dT's can be tolerated * HRecEffec(ipZ) is sum of totcap(i) * htotrc(ipZ) is sum of radiative recombination */ if( hcaseb.HRecCol[0]/hcaseb.htotrc[0] > 100. ) { HColRec.lgHColRec = TRUE; } else { HColRec.lgHColRec = FALSE; } /* this is test whether collisions are important first define ratio of exits * from ground due to coll, rel to everthing else, then flag is large */ if( hcolrt.hlte[0][IP1S] > 1e-30 ) { coltot = hydro.hcol[0][IP1S][IP2P]* hcolrt.hlte[0][IP2P]/hcolrt.hlte[0][IP1S]* phycon.eden*(HPhoto.hgamnc[0][IP2P] + hydro.hcolnc[0][IP2P] + hydro.hcol[0][IP2P][3]*hcolrt.hlte[0][3]/hcolrt.hlte[0][IP2P])/ (hydro.hcol[0][IP1S][IP2P]*phycon.eden + hydro.HRadNet[0][IP1S][IP2P]); /* caution that collisions are important (and so only small changes in * temperature should happen) if more than 20% of the total */ if( coltot > 0.2 ) { colion.lgHColionImp = TRUE; } } else { colion.lgHColionImp = FALSE; } /* remember the ratio of pops of 2p to 1s for possible printout in prtComment * and to obtain Lya excitation temps. the pop of ground is not defined if * NO ionization at all since these pops are relative to ion */ /* >>chng 99 Jun 03, added MAX2 to protect against totally neutral gas */ if( hydro.hn[0][IP2P]/MAX2(SMALLFLOAT,hydro.hn[0][IP1S]) > 0.1 ) { hydro.lgHiPop2 = TRUE; hydro.pop2mx = (float)MAX2(hydro.hn[0][IP2P]/hydro.hn[0][IP1S], hydro.pop2mx); } photex = 0.; for( i=IP2S; i <= nhlevel; i++ ) { photex += hydro.hn[0][i]*HPhoto.hgamnc[0][i]; } gamtot = HPhoto.hgamnc[0][IP1S] + Secondaries.csupra; coltot = hydro.hcolnc[0][IP1S] + hydro.hcol[0][IP1S][IP2P]*4.* hydro.hlbolt[0][IP1S][IP2P]; /* ground state destruction rate*/ grddst = gamtot + coltot*phycon.eden + Secondaries.x12tot + photex/MAX2(1e-37,hydro.hn[0][IP1S] ) ; if( grddst > 0. ) { photn2Com.photn2 = (float)(HPhoto.hgamnc[0][IP1S]/grddst ); } else { photn2Com.photn2 = 0.; } /* this flag is used in ionize to decide whether we * really need to converge the secondary ionization rates */ if( grddst > 0. ) { Secondaries.sec2total = (float)(Secondaries.csupra / grddst) ; } else { Secondaries.sec2total = 0. ; } /* now do 21 cm, will define populations here */ TauLines[ipH21cm].PopHi = hydro.hn[0][IP1S]*0.75f; TauLines[ipH21cm].PopLo = hydro.hn[0][IP1S]*0.25f; /* assume spin temperature is gas kinetic temperature */ TauLines[ipH21cm].PopOpc = hydro.hn[0][IP1S] * 0.25f * TauLines[ipH21cm].EnergyK / phycon.te; if( trace.lgTrace ) { fprintf( ioQQQ, " Hydrogenic retrn H,2,NE="); fprintf(ioQQQ,PrintEfmt("%9.2e", h.hi)); fprintf(ioQQQ,PrintEfmt("%9.2e", h.hii)); fprintf(ioQQQ,PrintEfmt("%9.2e", phycon.eden)); fprintf( ioQQQ, " REC, COL, GAMT= "); fprintf(ioQQQ,PrintEfmt("%9.2e", hcaseb.htotrc[0] )); fprintf(ioQQQ,PrintEfmt("%9.2e", coltot)); fprintf(ioQQQ,PrintEfmt("%9.2e", gamtot)); fprintf( ioQQQ, " CSUP="); PrintE82( ioQQQ, Secondaries.csupra); fprintf( ioQQQ, " Htwo="); fprintf(ioQQQ,PrintEfmt("%9.2e", hmi.htwo)); fprintf( ioQQQ, "\n"); } # ifdef DEBUG_FUN fputs( " <->Hydrogenic()\n", debug_fp ); # endif return; }