/*tauout set initial outward optical depths at start of first iteration, * it is only called by cloudy one time per complete calculation, just after * continuum set up and before start of convergence attempts. * * after first iteration, update updates the optical depths, mirroring this * routine but with the previous iteration's variables */ #include "cddefines.h" #include "physconst.h" #define TAULIM 1e8 #include "taulines.h" #include "nhe1lvl.h" #include "hydrogenic.h" #include "trace.h" #include "rfield.h" #include "turb.h" #include "showme.h" #include "xry.h" #include "taumin.h" #include "opac.h" #include "nh.h" #include "he1tau.h" #include "he3tau.h" #include "he1as.h" #include "he1nxt.h" #include "phycon.h" #include "qs.h" #include "elmton.h" #include "rt.h" #include "caseb.h" #include "sphere.h" #include "tseton.h" #include "double.h" #include "stopcalc.h" #include "taddlya.h" #include "solar.h" #include "converge.h" #include "tfidle.h" #include "printe82.h" #include "pop371.h" #include "emlinejunk.h" #include "hlife.h" #include "tauout.h" void tauout(void) { long int i, ipZ, ipHi, ipLo, limit ; long lgHit; /* used for logic check */ double AbunRatio, Vel, balc, coleff, f, factor , z4 ,/* used for saving ipZ+1**4 */ BalmerTauOn; # ifdef DEBUG_FUN fputs( "<+>tauout()\n", debug_fp ); # endif /* set initial and total optical depths in arrays * TAUNU is set when lines read in, sets stopping radius */ if( StopCalc.taunu > 0. ) { /* an optical depth has been specified */ if( StopCalc.iptnu >= nh.ipHn[0][IP1S] ) { /* at ionizing energies */ for( i=0; i < (nh.ipHn[0][IP1S] - 1); i++ ) { /* taumin set to 1e-20, can be reset with taumin command */ opac.tauabs[1][i] = tauminCom.taumin; opac.tausct[1][i] = tauminCom.taumin; opac.TauTotal[1][i] = tauminCom.taumin; } for( i=nh.ipHn[0][IP1S]-1; i < rfield.nupper; i++ ) { /* tauabs(i,2) = tauend * (anu(i)/anu(iptnu))**(-2.43) */ opac.tauabs[1][i] = StopCalc.tauend; opac.tausct[1][i] = tauminCom.taumin; opac.TauTotal[1][i] = opac.tauabs[1][i] + opac.tausct[1][i]; } } else { /* not specified at ionizing energies */ for( i=0; i < (nh.ipHn[0][IP1S] - 1); i++ ) { opac.tauabs[1][i] = StopCalc.tauend; opac.tausct[1][i] = StopCalc.tauend; opac.TauTotal[1][i] = StopCalc.tauend; } for( i=nh.ipHn[0][IP1S]-1; i < rfield.nupper; i++ ) { opac.tauabs[1][i] = (float)(TAULIM*pow(rfield.anu[i],-2.43)); opac.tausct[1][i] = tauminCom.taumin; opac.TauTotal[1][i] = opac.tauabs[1][i] + opac.tausct[1][i]; } } } else { /* ending optical depth not specified, assume 1E8 at 1 Ryd */ for( i=0; i < (nh.ipHn[0][IP1S] - 1); i++ ) { opac.tauabs[1][i] = tauminCom.taumin; opac.tausct[1][i] = tauminCom.taumin; opac.TauTotal[1][i] = tauminCom.taumin; } for( i=nh.ipHn[0][IP1S]-1; i < rfield.nupper; i++ ) { opac.tauabs[1][i] = (float)(TAULIM*pow(rfield.anu[i],-2.43)); opac.tausct[1][i] = tauminCom.taumin; opac.TauTotal[1][i] = opac.tauabs[1][i] + opac.tausct[1][i]; } } /* if lgSphere then double outer, set inner to half * assume will be optically thin at sub-ionizing energies */ if( sphere.lgSphere || caseb.lgCaseb ) { for( i=0; i < (nh.ipHn[0][IP1S] - 1); i++ ) { opac.tauabs[0][i] = tauminCom.taumin; opac.tauabs[1][i] = tauminCom.taumin*2.f; opac.tausct[0][i] = tauminCom.taumin; opac.tausct[1][i] = tauminCom.taumin*2.f; opac.TauTotal[0][i] = 2.f*tauminCom.taumin; opac.TauTotal[1][i] = 4.f*tauminCom.taumin; } for( i=nh.ipHn[0][IP1S]-1; i < rfield.nupper; i++ ) { opac.tauabs[0][i] = opac.tauabs[1][i]; opac.tauabs[1][i] *= 2.; opac.tausct[0][i] = opac.tausct[1][i]; opac.tausct[1][i] *= 2.; opac.TauTotal[0][i] = opac.TauTotal[1][i]; opac.TauTotal[1][i] *= 2.; } if( StopCalc.taunu > 0. ) { /* ending optical depth specified, and lgSphere also */ StopCalc.tauend *= 2.; } } /* fix escape prob for He, metals, first set log of Te, needed by RECEFF * do not do this if temperature has been set by command */ if( !tseton.lgTSetOn ) { phycon.te = 2e4; } /* propagate this temperature through the code */ tfidle(); /* set inward optical depths for hydrogenic ions to small number proportional to abundance */ for( ipZ=0; ipZ < LIMELM; ipZ++ ) { /* will need this for some scaling laws - physical Z to the 4th power*/ z4 = POW2(ipZ+1.); z4 *= z4; if( elmton.lgElmtOn[ipZ] ) { /* the actual branching ratio from high levels down to * 2s and 2p depend on the density. the code goes to * the correct low and high density limits - I know of * nothing better than can be done. */ factor = 0.32 - 0.07*phycon.eden/(phycon.eden + 1e7); /* upper limit for HydroBranch */ limit = MIN2(15,nhlevel); /* first reset line opacities and branching ratios */ for( ipHi=4; ipHi <= limit; ipHi++ ) { HydroLines[ipZ][ipHi][IP2S].Aul = (float)(hlife.HyLife[ipHi]* factor*HydroBranch(ipHi,2,ipZ+1)*z4); assert(HydroLines[ipZ][ipHi][IP2S].Aul > 0.); /* do 2p in terms of 2s */ HydroLines[ipZ][ipHi][IP2P].Aul = (float)( HydroLines[ipZ][ipHi][IP2S].Aul / factor * (1. - factor ) ); assert(HydroLines[ipZ][ipHi][IP2P].Aul > 0.); /* make self-consistent opaciity */ HydroLines[ipZ][ipHi][IP2S].opacity = (float)(HydroLines[ipZ][ipHi][IP2S].Aul* 2.2448e-26*HydroLines[ipZ][ipHi][IP2S].gHi/ HydroLines[ipZ][ipHi][IP2S].gLo* POW3(RYDLAM/HydroLines[ipZ][ipHi][IP2S].EnergyRyd)); assert(HydroLines[ipZ][ipHi][IP2S].opacity > 0.); HydroLines[ipZ][ipHi][IP2P].opacity = (float)(HydroLines[ipZ][ipHi][IP2P].Aul* 2.2448e-26*HydroLines[ipZ][ipHi][IP2P].gHi/ HydroLines[ipZ][ipHi][IP2P].gLo* POW3(RYDLAM/HydroLines[ipZ][ipHi][IP2P].EnergyRyd)); assert(HydroLines[ipZ][ipHi][IP2P].opacity > 0. ); for( ipLo=3; ipLo < ipHi; ipLo++ ) { HydroLines[ipZ][ipHi][ipLo].Aul = (float)(hlife.HyLife[ipHi]* HydroBranch(ipHi,ipLo,ipZ+1)*z4); assert(HydroLines[ipZ][ipHi][ipLo].Aul > 0. ); /* make self-consistent opacity, convert new As back into opacities */ HydroLines[ipZ][ipHi][ipLo].opacity = (float)(HydroLines[ipZ][ipHi][ipLo].Aul* 2.2448e-26*HydroLines[ipZ][ipHi][ipLo].gHi/ HydroLines[ipZ][ipHi][ipLo].gLo* POW3(RYDLAM/HydroLines[ipZ][ipHi][ipLo].EnergyRyd)); assert(HydroLines[ipZ][ipHi][ipLo].opacity > 0.); } } /* now get actual optical depths */ AbunRatio = solarCom.solar[ipZ]/solarCom.solar[0]; for( ipLo=IP1S; ipLo <= (nhlevel - 1); ipLo++ ) { for( ipHi=ipLo + 1; ipHi <= nhlevel; ipHi++ ) { /* set all inward optical depths to taumin, regardless of abundance * this is a very small number, 1e-20 */ HydroLines[ipZ][ipHi][ipLo].TauIn = tauminCom.taumin; } } /* La may be case B, tlamin set to 1e9 by default with case b command */ HydroLines[ipZ][IP2P][IP1S].TauIn = tauminCom.tlamin; /* scale factor so that alll other lyman lines are appropriate for this Lya optical depth*/ f = tauminCom.tlamin/HydroLines[ipZ][IP2P][IP1S].opacity; for( ipHi=3; ipHi <= nhlevel; ipHi++ ) { HydroLines[ipZ][ipHi][IP1S].TauIn = (float)(f*HydroLines[ipZ][ipHi][IP1S].opacity); assert( HydroLines[ipZ][ipHi][IP1S].TauIn > 0. ); } /* account for turbulence */ Vel = 1.28e6/sqrt(1.65e12+turb.TurbVel*turb.TurbVel); /* after this set of if's the total lya optical depth will be known, * common code later on will set rest of lyman lines * if case b then set total lyman to twice inner */ if( caseb.lgCaseb ) { /* force outer optical depth to twice inner if case B */ HydroLines[ipZ][IP2P][IP1S].TauTot = (float)(2.*HydroLines[ipZ][IP2P][IP1S].TauIn); /* force off Balmer et al optical depths */ BalmerTauOn = 0.; } else { /* usual case for H LyA; try to guess outer optical depth */ if( StopCalc.colnut < 6e29 ) { /* neutral column is defined */ HydroLines[ipZ][IP2P][IP1S].TauTot = (float)(StopCalc.colnut* double_.DoubleTau*7.6e-14*Vel*AbunRatio); assert( HydroLines[ipZ][IP2P][IP1S].TauTot > 0. ); } else if( StopCalc.taunu < 3. && StopCalc.taunu >= 0.99 ) { /* Lyman continuum optical depth defined */ HydroLines[ipZ][IP2P][IP1S].TauTot = (float)(StopCalc.tauend* 1.2e4*Vel*double_.DoubleTau*AbunRatio); assert( HydroLines[ipZ][IP2P][IP1S].TauTot > 0. ); } else if( StopCalc.HColStop < 6e29 ) { /* neutral column not defined, guess from total col and U */ coleff = StopCalc.HColStop - MIN2(MIN2(qs.qhtot/phycon.eden,1e24)/2.6e-13,0.6*StopCalc.HColStop); HydroLines[ipZ][IP2P][IP1S].TauTot = (float)(coleff* 7.6e-14*Vel*AbunRatio); assert( HydroLines[ipZ][IP2P][IP1S].TauTot > 0. ); } else { /* no way to estimate 912 optical depth, set to large number */ HydroLines[ipZ][IP2P][IP1S].TauTot = (float)(1e20* Vel*AbunRatio); assert( HydroLines[ipZ][IP2P][IP1S].TauTot > 0. ); } /* allow Balmer et al. optical depths */ BalmerTauOn = 1.; } /* Lya total optical depth now known, is it optically thick?*/ if( HydroLines[ipZ][IP2P][IP1S].TauTot > 1. ) { TAddLya.TAddHLya = (float)MIN2(1.,HydroLines[ipZ][IP2P][IP1S].TauTot/ 1e4); HydroLines[ipZ][IP2P][IP1S].TauIn += TAddLya.TAddHLya; } else { TAddLya.TAddHLya = tauminCom.tlamin; } /* this scale factor is to set other lyman lines, given the Lya optical depth */ f = HydroLines[ipZ][IP2P][IP1S].TauTot/ HydroLines[ipZ][IP2P][IP1S].opacity; for( ipHi=3; ipHi <= nhlevel; ipHi++ ) { /* set total optical depth for higher lyman lines */ HydroLines[ipZ][ipHi][IP1S].TauTot = (float)(HydroLines[ipZ][ipHi][IP1S].opacity* f); assert( HydroLines[ipZ][ipHi][IP1S].TauTot > 0.); /* increment inward optical depths by TAddLya all lyman lines, inward * optical depth was set above, this adds to it. this was originally * set some place above */ HydroLines[ipZ][ipHi][IP1S].TauIn += TAddLya.TAddHLya* HydroLines[ipZ][ipHi][IP1S].opacity/ HydroLines[ipZ][IP2P][IP1S].opacity; assert( HydroLines[ipZ][ipHi][IP1S].TauIn > 0.); } /* try to guess what Balmer cont optical guess, * first branch is case where we will stop at a balmer continuum optical * depth - taunu is energy where tauend was set */ if( StopCalc.taunu > 0.24 && StopCalc.taunu < 0.7 ) { HydroLines[ipZ][3][2].TauTot = (float)(StopCalc.tauend* 3.7e4*Vel*BalmerTauOn*AbunRatio + 1e-20); } else { /* this is a guess based on Ferland&Netzer 1979, but it gets very large */ balc = qs.qhtot*2.1e-19*BalmerTauOn*AbunRatio + 1e-20; /* Vel if the factor that account for turbulence, =1 is TurbVel=0 */ HydroLines[ipZ][3][2].TauTot = (float)(MIN2(2e5, balc*3.7e4*Vel*BalmerTauOn+1e-20)); assert( HydroLines[ipZ][3][2].TauTot > 0.); } /* 2s - 2p strongly forbidden */ HydroLines[ipZ][IP2P][IP2S].TauTot = 2.f*tauminCom.taumin; HydroLines[ipZ][IP2P][IP2S].TauIn = tauminCom.taumin; /* 2s - 1s strongly forbidden */ HydroLines[ipZ][IP2S][IP1S].TauTot = 2.f*tauminCom.taumin; HydroLines[ipZ][IP2S][IP1S].TauIn = tauminCom.taumin; /* fill in rest of the alpha transitions */ HydroLines[ipZ][4][3].TauTot = (float)(HydroLines[ipZ][3][2].TauTot*0.01*Vel*BalmerTauOn + tauminCom.taumin); HydroLines[ipZ][5][4].TauTot = (float)(HydroLines[ipZ][3][2].TauTot* 0.0005*Vel*BalmerTauOn + tauminCom.taumin); HydroLines[ipZ][6][5].TauTot = (float)(HydroLines[ipZ][3][2].TauTot*7.7e-5*Vel*BalmerTauOn + tauminCom.taumin); for( ipHi=7; ipHi <= nhlevel; ipHi++ ) { HydroLines[ipZ][ipHi][ipHi-1].TauTot = (float)(HydroLines[ipZ][6][5].TauTot/ipHi*BalmerTauOn + tauminCom.taumin); } /* transitions down to 2s are special since 'alpha' (2s-2p) has * small optical depth, so use 3 to 2p instead */ f = HydroLines[ipZ][3][IP2P].TauTot/ HydroLines[ipZ][3][IP2P].opacity* BalmerTauOn; ipLo = IP2S; for( ipHi=ipLo + 1; ipHi <= nhlevel; ipHi++ ) { HydroLines[ipZ][ipHi][ipLo].TauTot = (float)(tauminCom.taumin + f* HydroLines[ipZ][ipHi][ipLo].opacity); assert(HydroLines[ipZ][ipHi][ipLo].TauTot > 0.); } /* this is for rest of lines, scale from opacity */ for( ipLo=IP2P; ipLo <= (nhlevel - 1); ipLo++ ) { f = HydroLines[ipZ][ipLo+1][ipLo].TauTot/ HydroLines[ipZ][ipLo+1][ipLo].opacity* BalmerTauOn; for( ipHi=ipLo + 1; ipHi <= nhlevel; ipHi++ ) { HydroLines[ipZ][ipHi][ipLo].TauTot = (float)(tauminCom.taumin + f* HydroLines[ipZ][ipHi][ipLo].opacity); assert(HydroLines[ipZ][ipHi][ipLo].TauTot > 0.); } } /* this loop is over all possible H lines, do some final cleanup */ for( ipLo=IP1S; ipLo <= (nhlevel - 1); ipLo++ ) { for( ipHi=ipLo + 1; ipHi <= nhlevel; ipHi++ ) { /* TauCon is line optical depth to inner face used for continuum pumping rate, * not equal to TauIn in case of static sphere since TauCon will never * count the far side line optical depth */ HydroLines[ipZ][ipHi][ipLo].TauCon = HydroLines[ipZ][ipHi][ipLo].TauIn; /* make sure inward optical depth is not larger than half total */ HydroLines[ipZ][ipHi][ipLo].TauIn = MIN2( HydroLines[ipZ][ipHi][ipLo].TauIn , HydroLines[ipZ][ipHi][ipLo].TauTot/2.f ); assert(HydroLines[ipZ][ipHi][ipLo].TauIn > 0.); /* this is fraction of line that goes inward, not known until second iteration*/ HydroLines[ipZ][ipHi][ipLo].FracInwd = 0.5; } } } } /* begin sanity check, total greater than 1/0.9 time inward */ lgHit = FALSE; for( ipZ=0; ipZ < LIMELM; ipZ++ ) { if( elmton.lgElmtOn[ipZ] ) { for( ipLo=IP1S; ipLo <= (nhlevel - 1); ipLo++ ) { for( ipHi=ipLo + 1; ipHi <= nhlevel; ipHi++ ) { if( HydroLines[ipZ][ipHi][ipLo].TauTot < 0.9*HydroLines[ipZ][ipHi][ipLo].TauIn ) { fprintf(ioQQQ, "tauout insanity for h line, Z=%li lo=%li hi=%li tot=%g in=%g \n", ipZ , ipLo, ipHi , HydroLines[ipZ][ipHi][ipLo].TauTot , HydroLines[ipZ][ipHi][ipLo].TauIn ); lgHit = TRUE; } } } } } if( lgHit ) { fprintf( ioQQQ," stopping due to insanity in tauout\n"); ShowMe(); # ifdef DEBUG_FUN fputs( " <->tauout()\n", debug_fp ); # endif exit(1); } /*end sanity check */ /* fix offset for effective column density optical depth */ xry.tauxry = opac.tauabs[0][xry.ipxry-1]; for( i=0; i < (NHE1LVL - 1); i++ ) { for( ipLo=0; ipLo < NHE1LVL; ipLo++ ) { he1nxtCOM.he1nxt[i][ipLo] = tauminCom.taumin; he1tauCOM.he1tau[i][ipLo] = tauminCom.taumin; } } /* La may be case B */ he1tauCOM.he1tau[0][1] = tauminCom.tlamin; he1nxtCOM.he1nxt[0][1] = tauminCom.tlamin; f = tauminCom.tlamin/he1tauCOM.he1opc[0][1]; for( i=2; i < NHE1LVL; i++ ) { he1tauCOM.he1tau[0][i] = (float)(f*he1tauCOM.he1opc[0][i]); he1nxtCOM.he1nxt[0][i] = he1tauCOM.he1tau[0][i]; } /* fix LA in ways of increasing accuracy */ /* account for turbulence */ Vel = 1.28e6/sqrt(1.65e12+turb.TurbVel*turb.TurbVel); if( caseb.lgCaseb ) { /* force outer optical depth to twice inner if case B */ he1tauCOM.he1lim[0][1] = (float)(2.*he1nxtCOM.he1nxt[0][1]); /* force off Balmer et al optical depths */ BalmerTauOn = 0.; } else { he1tauCOM.he1lim[0][1] = HydroLines[0][IP2P][IP1S].TauTot; /* allow Balmer et al. optical depths */ BalmerTauOn = 1.; } if( he1tauCOM.he1lim[0][1] > 1. ) { TAddLya.TAddHeI = (float)MIN2(1.,he1tauCOM.he1lim[0][1]/1e4); he1tauCOM.he1tau[0][1] += TAddLya.TAddHeI; he1nxtCOM.he1nxt[0][1] = he1tauCOM.he1tau[0][1]; } else { TAddLya.TAddHeI = tauminCom.tlamin; } /* min triggers if case b command entered */ TAddLya.TAddHeI = (float)MIN2(he1tauCOM.he1lim[0][1]/10.,he1tauCOM.he1tau[0][1]); f = he1tauCOM.he1lim[0][1]/he1tauCOM.he1opc[0][1]; for( i=2; i < NHE1LVL; i++ ) { he1tauCOM.he1lim[0][i] = (float)(he1tauCOM.he1opc[0][i]*f); he1tauCOM.he1tau[0][i] += TAddLya.TAddHeI*he1tauCOM.he1opc[0][i]/ he1tauCOM.he1opc[0][1]; he1nxtCOM.he1nxt[0][i] = he1tauCOM.he1tau[0][i]; } for( i=1; i < (NHE1LVL - 1); i++ ) { for( ipLo=i + 1; ipLo <= NHE1LVL; ipLo++ ) { he1tauCOM.he1lim[i][ipLo] = tauminCom.taumin; } } /* set helium triplet optical depths */ for( i=0; i < NHE3TAU; i++ ) { /* these are line optical depth arrays * inward optical depth */ he3tau[i].TauIn = tauminCom.taumin; he3tau[i].TauCon = tauminCom.taumin; he3tau[i].ColOvTot = 0.; /* outward optical depth */ he3tau[i].TauTot = 1e30f; /* escape probability */ he3tau[i].Pesc = 1.; /* inward part of line */ he3tau[i].FracInwd = 1.; /* destruction probability */ he3tau[i].Pdest = 0.; /* line pumping rate */ he3tau[i].pump = 0.; /* population of lower level */ he3tau[i].PopLo = 0.; /* >>chng 97 july 21, added following zero * population of upper level */ he3tau[i].PopHi = 0.; /* population of lower level with correction for stim emission */ he3tau[i].PopOpc = 0.; /* following two heat exchange excitation, deexcitation */ he3tau[i].cool = 0.; he3tau[i].heat = 0.; /* intensity of line */ he3tau[i].xIntensity = 0.; /* number of photons emitted in line */ he3tau[i].phots = 0.; /* opacity in line */ he3tau[i].dTau = 0.; /* indicates that this is a high quality transition */ he3tau[i].cs1 = 0.; } /* HeI 10830 when static */ if( sphere.lgStatic ) { he3tau[IPT10830-1].TauIn = 100.; he3tau[IPT10830-1].TauTot = 200.; } /* initialize heavy element line array */ for( i=1; i <= nLevel1; i++ ) { /* these are line optical depth arrays * inward optical depth */ TauLines[i].TauIn = tauminCom.taumin; TauLines[i].TauCon = tauminCom.taumin; TauLines[i].ColOvTot = 0.; /* outward optical depth */ TauLines[i].TauTot = 1e30f; /* escape probability */ TauLines[i].Pesc = 1.; /* inward part of line */ TauLines[i].FracInwd = 1.; /* destruction probability */ TauLines[i].Pdest = 0.; /* line pumping rate */ TauLines[i].pump = 0.; /* population of lower level */ TauLines[i].PopLo = 0.; /* >>chng 97 july 21, added following zero * population of upper level */ TauLines[i].PopHi = 0.; /* population of lower level with correction for stim emission */ TauLines[i].PopOpc = 0.; /* following two heat exchange excitation, deexcitation */ TauLines[i].cool = 0.; TauLines[i].heat = 0.; /* intensity of line */ TauLines[i].xIntensity = 0.; /* number of photons emitted in line */ TauLines[i].phots = 0.; /* opacity in line */ TauLines[i].dTau = 0.; /* indicates that this is a high quality transition */ TauLines[i].cs1 = 0.; } /* this is a dummy optical depth array for non-existant lines * when this goes over to struc, make sure all are set to zero here since * init in grid may depend on it */ /*for( i=0; i < NTA; i++ )*/ /*{*/ /*TauDummy.t[i] = 0.;*/ /*}*/ EmLineZero( &TauDummy ); /* lines in cooling function with Mewe approximate collision strengths */ for( i=0; i < nWindLine; i++ ) { /* these are line optical depth arrays * inward optical depth */ TauLine2[i].TauIn = tauminCom.taumin; TauLine2[i].TauCon = tauminCom.taumin; TauLine2[i].ColOvTot = 0.; /* outward optical depth */ TauLine2[i].TauTot = 1e20f; /* escape probability */ TauLine2[i].Pesc = 1.f; /* inward part of line */ TauLine2[i].FracInwd = 1.f; /* destruction probability */ TauLine2[i].Pdest = 0.; /* line pumping rate */ TauLine2[i].pump = 0.; /* population of lower level */ TauLine2[i].PopLo = 0.; /* population of upper level */ TauLine2[i].PopHi = 0.; /* population of lower level with correction for stim emission */ TauLine2[i].PopOpc = 0.; /* following two heat exchange excitation, deexcitation */ TauLine2[i].cool = 0.; TauLine2[i].heat = 0.; /* intensity of line */ TauLine2[i].xIntensity = 0.; /* number of photons emitted in line */ TauLine2[i].phots = 0.; /* opacity in line */ TauLine2[i].dTau = 0.; } /* initialize Dima and Katya Verner's FeII arrays */ if( FeII.lgFeIION ) { FeIITauInit(); } /* this is flag detected by dest prob routines to see whether ots rates are * oscaillating - damp them out if so */ conv.lgOscilOTS = FALSE; /* now that optical depths have been incremented, do escape prob, this * is located here instead on in cloudy.c (where it belongs) because * information generated by RTMake is needed for the following printout */ RTMake(TRUE); /* rest of routine is printout in case of trace */ if( trace.lgTrace ) { if( trace.lgHBug && trace.lgHBugFull ) { fprintf( ioQQQ, "\n\n up HydroLines TauTot array as setin TAUOUT ipZTrace=%3ld\n", trace.ipZTrace ); for( ipHi=2; ipHi <= nhlevel; ipHi++ ) { fprintf( ioQQQ, " %3ld", ipHi ); for( ipLo=IP1S; ipLo < ipHi; ipLo++ ) { fprintf( ioQQQ,PrintEfmt("%9.2e", HydroLines[trace.ipZTrace][ipHi][ipLo].TauTot )); } fprintf( ioQQQ, "\n" ); } fprintf( ioQQQ, "\n\n TauIn array\n" ); for( ipHi=1; ipHi <= nhlevel; ipHi++ ) { fprintf( ioQQQ, " %3ld", ipHi ); for( ipLo=IP1S; ipLo < ipHi; ipLo++ ) { fprintf( ioQQQ,PrintEfmt("%9.2e", HydroLines[trace.ipZTrace][ipHi][ipLo].TauIn )); } fprintf( ioQQQ, "\n" ); } fprintf( ioQQQ, "\n\n Aul As array\n" ); for( ipHi=1; ipHi <= nhlevel; ipHi++ ) { fprintf( ioQQQ, " %3ld", ipHi ); for( ipLo=IP1S; ipLo < ipHi; ipLo++ ) { fprintf( ioQQQ,PrintEfmt("%9.2e", HydroLines[trace.ipZTrace][ipHi][ipLo].Aul) ); } fprintf( ioQQQ, "\n" ); } fprintf( ioQQQ, "\n\n Aul*esc array\n" ); for( ipHi=1; ipHi <= nhlevel; ipHi++ ) { fprintf( ioQQQ, " %3ld", ipHi ); for( ipLo=IP1S; ipLo < ipHi; ipLo++ ) { fprintf( ioQQQ,PrintEfmt("%9.2e", HydroLines[trace.ipZTrace][ipHi][ipLo].Aul* (HydroLines[trace.ipZTrace][ipHi][ipLo].Pdest + HydroLines[trace.ipZTrace][ipHi][ipLo].Pesc) )); } fprintf( ioQQQ, "\n" ); } fprintf( ioQQQ, "\n\n H opac array\n" ); for( ipHi=1; ipHi <= nhlevel; ipHi++ ) { fprintf( ioQQQ, " %3ld", ipHi ); for( ipLo=IP1S; ipLo < ipHi; ipLo++ ) { fprintf( ioQQQ,PrintEfmt("%9.2e", HydroLines[trace.ipZTrace][ipHi][ipLo].opacity )); } fprintf( ioQQQ, "\n" ); } } else if( trace.lgHe1Bug ) { fprintf( ioQQQ, "\n\n HE1LIM array as set in TAUOUT\n" ); for( ipHi=1; ipHi < NHE1LVL; ipHi++ ) { fprintf( ioQQQ, " %3ld", ipHi ); for( ipLo=0; ipLo < (NHE1LVL - 1); ipLo++ ) { fprintf( ioQQQ,PrintEfmt("%9.2e", he1tauCOM.he1lim[ipLo][ipHi] )); } fprintf( ioQQQ, "\n" ); } fprintf( ioQQQ, "\n\n HE1TAU array\n" ); for( ipHi=1; ipHi < NHE1LVL; ipHi++ ) { fprintf( ioQQQ, " %3ld", ipHi ); for( ipLo=0; ipLo < (NHE1LVL - 1); ipLo++ ) { fprintf( ioQQQ,PrintEfmt("%9.2e", he1tauCOM.he1tau[ipLo][ipHi] )); } fprintf( ioQQQ, "\n" ); } fprintf( ioQQQ, "\n\n HE1 As array\n" ); for( ipHi=1; ipHi < NHE1LVL; ipHi++ ) { fprintf( ioQQQ, " %3ld", ipHi ); for( ipLo=0; ipLo < (NHE1LVL - 1); ipLo++ ) { fprintf( ioQQQ,PrintEfmt("%9.2e", he1as.He1EinA[ipLo][ipHi] )); } fprintf( ioQQQ, "\n" ); } fprintf( ioQQQ, "\n\n HE1 opac array\n" ); for( ipHi=1; ipHi < NHE1LVL; ipHi++ ) { fprintf( ioQQQ, " %3ld", ipHi ); for( ipLo=0; ipLo < (NHE1LVL - 1); ipLo++ ) { fprintf( ioQQQ,PrintEfmt("%9.2e", he1tauCOM.he1opc[ipLo][ipHi] )); } fprintf( ioQQQ, "\n" ); } } else { fprintf( ioQQQ, " TAUOUT called.\n" ); } } # ifdef DEBUG_FUN fputs( " <->tauout()\n", debug_fp ); # endif return; }