/*DoPunch produce punch output during calculation */ /* chTime is 'MIDL' during calculation, 'LAST' at the end */ #include #include "cddefines.h" #include "showme.h" #include "physconst.h" #include "taulines.h" #include "hydrogenic.h" #include "nhe1lvl.h" #include "secondaries.h" #include "nhe3lvl.h" #include "grainvar.h" #include "linesave.h" #include "pnunit.h" #include "opacpoint.h" #include "tepowers.h" #include "mappar.h" #include "colden.h" #include "itercnt.h" #include "hevmolec.h" #include "flxsav.h" #include "converge.h" #include "heating.h" #include "sphere.h" #include "called.h" #include "compton.h" #include "opac.h" #include "rfield.h" #include "phycon.h" #include "timesc.h" #include "radius.h" #include "pop371.h" #include "assertresults.h" #include "h.h" #include "he.h" #include "zonecnt.h" #include "heat.h" #include "nh.h" #include "phe1lv.h" #include "heopfr.h" #include "dclgas.h" #include "maxcor.h" #include "wind.h" #include "deptau.h" #include "hmi.h" #include "levery.h" #include "ipvfil.h" #include "tlast.h" #include "punchskip.h" #include "linelabl.h" #include "heatcool.h" #include "fourpi.h" #include "hcaseb.h" #include "ionfracs.h" #include "pressure.h" #include "tpred.h" #include "he3rate.h" #include "he1rate.h" #include "herec.h" #include "he3gams.h" #include "nhe1.h" #include "nsshells.h" #include "ph1com.h" #include "elementnames.h" #include "shell.h" #include "sumoutinc.h" #include "plwidth.h" #include "hphoto.h" #include "densty.h" #include "ipoint.h" #include "printe82.h" #include "punlin.h" #include "pfeii.h" #include "puncool.h" #include "pungamma.h" #include "pgaunt.h" #include "punheat.h" #include "gammas.h" #include "pllabels.h" #include "pldata.h" #include "plinin.h" #include "result.h" #include "punspec.h" #include "punopac.h" #include "dopunch.h" #include "hchargtran.h" #include "punt.h" #include "anuunit.h" #include "pmeani.h" #include "plankf.h" void DoPunch( char *chTime) /* chTime is null terminated 4 char string, either "MIDL" or "LAST" */ { char chEner[LIMELM][7]; int lgOK, lgTlkSav; FILE * isav; long int i, i1, ion, ipConMax, ipLinMax, ips, j, jj, limit, n, nd, nelem; double ConMax, RateInter, av, conem, damp, eps, flxatt, flxcor, flxin, flxref, flxtrn, fout, fref, fsum, opConSum, opLinSum, pop, refac, stage, sum, texc, xLinMax; # ifdef DEBUG_FUN fputs( "<+>DoPunch()\n", debug_fp ); # endif /* * the "last" option on punch command, to punch on last iteration, * is parsed at the top of the loop in only one place. * no further action is needed at all for punch last to work * ok throughout this routine */ /* * each branch can have a test whether chTime is or is not "LAST" * if( strcmp(chTime,"LAST") == 0 ) * if "LAST" then this is last call to routine after calc complete * punch only if "LAST" when results at end of cal are needed * test for .not."LAST" is for every zone result, where you do not * want to punch last zone twice */ /* punch results on unit ipPnunit */ if( pnunit.npunch < 1 ) { # ifdef DEBUG_FUN fputs( " <->DoPunch()\n", debug_fp ); # endif return; } for( i=0; i < pnunit.npunch; i++ ) { /* * lgPunLstIter set true if 'last' key occured on punch command * normally is false. This will skip punching if last set and * this is not last iteration */ if( IterCnt.lgLastIt || (!pnunit.lgPunLstIter[i]) ) { if( strcmp(pnunit.chPunch[i],"ABUN") == 0 ) { /* punch abundances vs depth */ if( strcmp(chTime,"LAST") != 0 ) { for( j=1; j <= LIMELM; j++ ) { fprintf( pnunit.ipPnunit[i], "%6.2f", log10(xIonFracs[0][j-1]) ); } fprintf( pnunit.ipPnunit[i], "\n" ); } } else if( strcmp(pnunit.chPunch[i],"AGES") == 0 ) { /* punch timescales vs depth */ if( strcmp(chTime,"LAST") != 0 ) { fprintf( pnunit.ipPnunit[i], "%10.2e%10.2e%10.2e%10.2e%10.2e\n", /* depth, cm */ radius.depth, /* cooling timescale */ densty.pden*BOLTZMANN*1.5*phycon.te/ heat.htot, /* H2 formation timescale */ timesc.AgeH2MoleDest, /* CO formation timescale */ timesc.AgeCOMoleDest, /* H recombination timescale */ 1./(phycon.eden*2.90e-10/(TePowers.te70*TePowers.te10/TePowers.te03)) ); } } else if( strcmp(pnunit.chPunch[i],"ASSE") == 0 ) { if( strcmp(chTime,"LAST") == 0 ) { /* punch the assert output */ /*lint -e534 ignore return value of following - nothing to do with it */ lgCheckAsserts(pnunit.ipPnunit[i]); /*lint +e534 ignore return value of following - nothing to do with it */ } } else if( strcmp(pnunit.chPunch[i],"CHAR") == 0 ) { /* charge exchange rate coefficients, entered with * punch charge transfer command. Queries Jim Kingdon's * CT fits and routines to get H - He and higher CT rates * rates are evaluated at the current temperature, which can * be specified with the constant temperature commmand */ if( strcmp(chTime,"LAST") != 0 ) { /* first group is charge transfer recombination, * the process A+x + H => A+x-1 + H+ */ fprintf( pnunit.ipPnunit[i], "recombination rates, atomic number\n"); for( j=1; j < LIMELM; j++ ) { fprintf( pnunit.ipPnunit[i], "%3ld", j ); /* only do first 6 unless element very light, then do less */ limit = MIN2(6,j); for( jj=0; jj < limit; jj++ ) { fprintf( pnunit.ipPnunit[i], "%10.2e", HCTRecom(jj,j) ); } fprintf( pnunit.ipPnunit[i], "\n" ); } /* second group is charge transfer ionization, * the process A+x + H+ => A+x+1 + H */ fprintf( pnunit.ipPnunit[i], "\nionization rates, atomic number\n"); for( j=1; j < LIMELM; j++ ) { limit = MIN2(6,j); fprintf( pnunit.ipPnunit[i], "%3ld", j ); for( jj=0; jj < limit; jj++ ) { fprintf( pnunit.ipPnunit[i], "%10.2e", HCTIon(jj,j) ); } fprintf( pnunit.ipPnunit[i], "\n" ); } } } else if( strcmp(pnunit.chPunch[i],"COOL") == 0 ) { if( strcmp(chTime,"LAST") != 0 ) PunCool(pnunit.ipPnunit[i]); } else if( strcmp(pnunit.chPunch[i],"COMP") == 0 ) { /* compton energy exchange coefficients */ if( strcmp(chTime,"LAST") != 0 ) { for( jj=0; jj>chng 96 oct 22, format of anu to 11.5 to resolve grid near 1 ryd */ fprintf( pnunit.ipPnunit[i], "%12.4e %10.3e\n", AnuUnit(j+1), flxin ); } } } else if( strcmp(pnunit.chPunch[i],"CONR") == 0 ) { /* punch reflected continuum */ /* set pointer for possible change in units of energy in continuum * AnuUnit will give anu in whatever units were set with punch units */ pnunit.ipConPunEnr = i; if( strcmp(chTime,"LAST") == 0 ) { if( sphere.lgSphere ) { fprintf( pnunit.ipPnunit[i], " Reflected continuum not predicted when SPHERE is set.\n" ); puts( "[Stop in DoPunch]" ); exit(1); } for( j=0; j 1e-25 ) { av = rfield.reflux[j]/flxsav.FluxSave[j]; } else { av = 0.; } /* >>chng 96 oct 22, format of anu to 12.5 to resolve grid near 1 ryd */ fprintf( pnunit.ipPnunit[i], "%12.5e%12.5e%12.5e%12.5e%12.5e %4.4s\n", AnuUnit(j+1), flxref, fref, flxref + fref, av, LineLabl.chContLabel[j] ); } } } else if( strcmp(pnunit.chPunch[i],"CONB") == 0 ) { /* punch continuum binning */ /* set pointer for possible change in units of energy in continuum * AnuUnit will give anu in whatever units were set with punch units */ pnunit.ipConPunEnr = i; if( strcmp(chTime,"LAST") != 0 ) { for( j=1; j<=rfield.nupper; j = j + PunchSkip.ncPunchSkip) { fprintf( pnunit.ipPnunit[i], "%14.5e%14.5e\n", AnuUnit(j), rfield.widflx[j-1] ); } } } else if( strcmp(pnunit.chPunch[i],"COND") == 0 ) { /* punch diffuse continuum */ /* set pointer for possible change in units of energy in continuum * AnuUnit will give anu in whatever units were set with punch units */ pnunit.ipConPunEnr = i; if( strcmp(chTime,"LAST") == 0 ) { /* this option to punch diffuse continuum */ for( j=0; j>chng 96 oct 22, format of anu to 12.5 to resolve grid near 1 ryd */ fprintf( pnunit.ipPnunit[i], "%12.5e %12.5e\n", AnuUnit(j+1), rfield.DiffLocal[j]*rfield.anu2[j]* EN1RYD/rfield.widflx[j] ); } } } else if( strcmp(pnunit.chPunch[i],"CONE") == 0 ) { /* punch emitted continuum */ /* set pointer for possible change in units of energy in continuum * AnuUnit will give anu in whatever units were set with punch units */ pnunit.ipConPunEnr = i; if( strcmp(chTime,"LAST") == 0 ) { /* punch emitted continuum */ if( sphere.lgSphere ) { refac = 0.; } else { refac = 1.; } for( j=0; j>chng 96 oct 22, format of anu to 11.5 to resolve grid near 1 ryd */ fprintf( pnunit.ipPnunit[i], "%11.5e%11.3e%11.3e%11.3e %4.4s %4.4s \n", AnuUnit(j+1), flxref, conem, flxref + conem, LineLabl.chLineLabel[j], LineLabl.chContLabel[j] ); } } } else if( strcmp(pnunit.chPunch[i],"CONi") == 0 ) { /* punch continuum interactions */ /* set pointer for possible change in units of energy in continuum * AnuUnit will give anu in whatever units were set with punch units */ pnunit.ipConPunEnr = i; /* continuum interactions */ if( strcmp(chTime,"LAST") != 0 ) { /* this is option to set lowest energy */ if( pnunit.punarg[0][i] <= 0. ) { i1 = 1; } else if( pnunit.punarg[0][i] < 100. ) { i1 = ipoint(pnunit.punarg[0][i]); } else { i1 = (long int)pnunit.punarg[0][i]; } fref = 0.; fout = 0.; fsum = 0.; sum = 0.; flxin = 0.; for( j=i1-1; j < rfield.nflux; j++ ) { fref += rfield.flux[j]*opac.opac[j]; fout += rfield.otslin[j]*opac.opac[j]; fsum += rfield.otscon[j]*opac.opac[j]; sum += rfield.outcon[j]*opac.opac[j]; flxin += rfield.outlin[j]*opac.opac[j]; } fprintf( pnunit.ipPnunit[i], "%10.2e%10.2e%10.2e%10.2e%10.2e\n", fref, fout, fsum, sum, flxin ); } } else if( strcmp(pnunit.chPunch[i],"CONI") == 0 ) { /* punch ionizing continuum */ /* set pointer for possible change in units of energy in continuum * AnuUnit will give anu in whatever units were set with punch units */ pnunit.ipConPunEnr = i; if( strcmp(chTime,"LAST") == 0 ) { /* punch ionizing continuum command * this is option to set lowest energy */ if( pnunit.punarg[0][i] <= 0. ) { i1 = 1; } else if( pnunit.punarg[0][i] < 100. ) { i1 = ipoint(pnunit.punarg[0][i]); } else { i1 = (long int)pnunit.punarg[0][i]; } sum = 0.; for( j=i1-1; j < rfield.nflux; j++ ) { flxcor = rfield.flux[j] + rfield.otslin[j] + rfield.otscon[j] + rfield.outcon[j] + rfield.outlin[j]; sum += flxcor*opac.opac[j]; } if( sum > 0. ) { sum = 1./sum; } else { sum = 1.; } fsum = 0.; for( j=i1-1; j pnunit.punarg[1][i] ) { /* >>chng 96 oct 22, format of anu to 11.5 to resolve grid near 1 ryd */ fprintf( pnunit.ipPnunit[i], "%11.5e%10.2e%10.2e%5.2f%5.2f%5.2f%5.2f%10.2e%10.2e %4.4s %4.4s \n", AnuUnit(j+1), flxcor, flxcor*opac.opac[j], rfield.flux[j]/flxcor, rfield.otslin[j]/ flxcor, rfield.otscon[j]/flxcor, (rfield.outcon[j] + rfield.outlin[j])/flxcor, RateInter, fsum*sum, LineLabl.chLineLabel[j], LineLabl.chContLabel[j] ); } } } } else if( strcmp(pnunit.chPunch[i],"CORA") == 0 ) { /* punch raw continuum */ /* set pointer for possible change in units of energy in continuum * AnuUnit will give anu in whatever units were set with punch units */ pnunit.ipConPunEnr = i; if( strcmp(chTime,"LAST") == 0 ) { /* this option to punch all raw ionizing continuum */ for( j=0;j>chng 97 jul 10, remove PunchLWidth from this one only since * we must conserve energy even in lines */ flxatt = rfield.flux[j]*rfield.anu2[j]*EN1RYD/ rfield.widflx[j]*radius.r1r0sq; conem = (rfield.ConOutNoInter[j] + rfield.outcon[j]+rfield.outlin[j])* rfield.anu2[j] ; conem *= radius.r1r0sq*EN1RYD*sphere.covgeo; flxtrn = conem + flxatt; /* use AnuOrg here instead of anu since probably * going to be used by table read * and we want good anu array for sanity check*/ fprintf( pnunit.ipPnunit[i],PrintEfmt("%10.3e", AnuUnit(j+1) )); fprintf( pnunit.ipPnunit[i],PrintEfmt("%10.3e", flxtrn )); fprintf( pnunit.ipPnunit[i], "\n"); } } } else if( strcmp(pnunit.chPunch[i],"DUSO") == 0 ) { /* punch grain opacity */ if( strcmp(chTime,"LAST") == 0 ) { for( j=1; j <= rfield.nflux; j++ ) { fprintf( pnunit.ipPnunit[i], "%12.4e%10.2e%10.2e%10.2e\n", rfield.anu[j-1], GrainVar.dstab[j-1] + GrainVar.dstsc[j-1], GrainVar.dstab[j-1], GrainVar.dstsc[j-1] ); } } } else if( strcmp(pnunit.chPunch[i],"DUSP") == 0 ) { /* punch grain physical conditions */ fprintf( pnunit.ipPnunit[i], " %4ld", ZoneCnt.nzone ); for( j=1; j <= GrainVar.nDustOn; j++ ) { fprintf( pnunit.ipPnunit[i], "%10.2e", GrainVar.tedust[GrainVar.ipdon[j-1]-1] ); } for( j=1; j <= GrainVar.nDustOn; j++ ) { fprintf( pnunit.ipPnunit[i], "%10.2e", GrainVar.dstpot[GrainVar.ipdon[j-1]-1]*EVRYD ); } for( j=1; j <= GrainVar.nDustOn; j++ ) { fprintf( pnunit.ipPnunit[i], "%10.2e", GrainVar.DustDftVel[GrainVar.ipdon[j-1]-1] ); } fprintf( pnunit.ipPnunit[i], "%10.2e%10.2e\n", HeatingCom.heating[13][0]/heat.htot, dclgas.DustCoolGas/heat.htot ); } else if( strcmp(pnunit.chPunch[i],"DUSQ") == 0 ) { /* punch grain Qs */ if( strcmp(chTime,"LAST") == 0 ) { for( j=0; j < rfield.nflux; j++ ) { fprintf( pnunit.ipPnunit[i], "%11.4e", rfield.anu[j] ); for( nd=0; nd < GrainVar.nDustOn; nd++ ) { fprintf( pnunit.ipPnunit[i], "%9.1e%9.1e", GrainVar.dqabs[GrainVar.ipdon[nd]-1][j], GrainVar.dqscat[GrainVar.ipdon[nd]-1][j] ); } fprintf( pnunit.ipPnunit[i], "\n" ); } } } else if( strcmp(pnunit.chPunch[i],"ELEM") == 0 ) { if( strcmp(chTime,"LAST") != 0 ) { /* this is the pointer to atomic number */ jj = (long int)pnunit.punarg[0][i]; fprintf( pnunit.ipPnunit[i], " %11.4e", radius.depth ); for( j=1; j <= (jj + 1); j++ ) { fprintf( pnunit.ipPnunit[i], "%5ld", (long)(-1000.*log10(MAX2(5e-10,xIonFracs[j][jj-1]/ xIonFracs[0][jj-1]))) ); } fprintf( pnunit.ipPnunit[i], "\n" ); } } else if( strcmp(pnunit.chPunch[i],"RECE") == 0 ) { /* punch recombination efficiencies, * option turned on with the "punch recombination efficiencies" command * output for the punch recombination coefficients command is actually * produced by a series of routines, as they generate the recombination * coefficients. these include * dielsupres, helium, hydrorecom, iibod, and makerecomb*/ fprintf( pnunit.ipPnunit[i], "%12.4e %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e\n", hydro.hrec[0][0][ipRecRad], hydro.hrec[0][0][ipRecNetEsc] , hydro.hrec[0][2][ipRecRad], hydro.hrec[0][2][ipRecNetEsc], phe1lv.he1rec[1][0], phe1lv.he1rec[2][0], heopfr.ophe1f[0] ); } else if( strcmp(pnunit.chPunch[i],"FE2f") == 0 ) { /* set with punch feii fred command, this punches some stuff from * fred hamann's feii atom */ if( strcmp(chTime,"LAST") != 0 ) pfeii(pnunit.ipPnunit[i]); } else if( strcmp(pnunit.chPunch[i],"FE2d") == 0 ) { /* punch some departure cefficients for large feii atom */ if( strcmp(chTime,"LAST") != 0 ) FeIIPunDepart(pnunit.ipPnunit[i],FALSE); } else if( strcmp(pnunit.chPunch[i],"FE2D") == 0 ) { /* punch all departure cefficients for large feii atom */ if( strcmp(chTime,"LAST") != 0 ) FeIIPunDepart(pnunit.ipPnunit[i],TRUE); } else if( strcmp(pnunit.chPunch[i],"FE2p") == 0 ) { /* punch some level populations for large feii atom */ if( strcmp(chTime,"LAST") != 0 ) FeIIPunPop(pnunit.ipPnunit[i],FALSE); } else if( strcmp(pnunit.chPunch[i],"FE2P") == 0 ) { /* punch all level populations for large feii atom */ if( strcmp(chTime,"LAST") != 0 ) FeIIPunPop(pnunit.ipPnunit[i],TRUE); } else if( strcmp(pnunit.chPunch[i],"GAMM") == 0 ) { if( strcmp(chTime,"LAST") != 0 ) PunGamma(pnunit.ipPnunit[i]); } else if( strcmp(pnunit.chPunch[i],"GAUN") == 0 ) { if( strcmp(chTime,"LAST") != 0 ) pgaunt(pnunit.ipPnunit[i]); } else if( strcmp(pnunit.chPunch[i],"HEAT") == 0 ) { if( strcmp(chTime,"LAST") != 0 ) PunHeat(pnunit.ipPnunit[i]); /* punch helium */ } /* informatin ou the helium singlets */ else if( strcmp(pnunit.chPunch[i],"HEL1") == 0 ) { if( strcmp(chTime,"LAST") != 0 ) { if( he1rate.smhe2ov1 > 1e-30 ) { /* this is simple estimate of singlet to total he */ fref = (xIonFracs[2][ipHE]/he1rate.smhe2ov1)/xIonFracs[0][ipHE]; } else { fref = 0.; } /* second number is ratio of predicted HeII/AheII to simple * clhe2ov1 is computed He1/He1 */ fprintf( pnunit.ipPnunit[i], "he1%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e\n", xIonFracs[1][ipHE]/xIonFracs[0][ipHE], fref, he1rate.clhe2ov1, he1rate.smhe2ov1, herecCom.reci*phycon.eden, phe1lv.he1gam[0], he.hei3/xIonFracs[0][ipHE] ); /* print out photoionization rates */ GammaPrt(nhe1Com.nhe1[0],rfield.nflux,OpacPoint.iophe1[0], pnunit.ipPnunit[i],phe1lv.he1gam[0],phe1lv.he1gam[0]* 0.05); } } /* information about the helium triplets */ else if( strcmp(pnunit.chPunch[i],"HEL3") == 0 ) { if( strcmp(chTime,"LAST") != 0 ) { /* this is simple estimate of triplet to total he */ if( he3rate.she2ov3 > 0. ) { fref = (xIonFracs[2][ipHE]/he3rate.she2ov3)/xIonFracs[0][ipHE]; } else { fref = 0.; } fprintf( pnunit.ipPnunit[i], "HeTrip%10.2e%10.2e%10.2e%10.2e%10.2e\n", he.hei3/xIonFracs[0][ipHE], fref, herecCom.rectri* phycon.eden, he3gams.he3gam[0], he3rate.collhe3 ); /* print out photoionization rates */ GammaPrt(he.nhei3,nhe1Com.nhe1[0],OpacPoint.ioptri, pnunit.ipPnunit[i] , he3gams.he3gam[0] , he3gams.he3gam[0]*0.05); } } /* punch hummer, results needed for Lya transport, to feed into David's routine */ else if( strcmp(pnunit.chPunch[i],"HUMM") == 0 ) { damp = HydroLines[0][IP2P][IP1S].damp; eps = HydroLines[0][IP2P][IP1S].Aul/ (hydro.hcol[0][IP1S][IP2P]*phycon.eden); fprintf( pnunit.ipPnunit[i], " %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n", radius.depth, HydroLines[0][IP2P][IP1S].TauIn, hydro.hn[0][IP1S]*h.hii, hydro.hn[0][IP2P]*h.hii, phycon.te, damp, eps ); } else if( strcmp(pnunit.chPunch[i],"HYDc") == 0 ) { /* punch hydrogen physical conditions */ fprintf( pnunit.ipPnunit[i], " %10.2e%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e\n", radius.depth, phycon.te, phycon.hden, phycon.eden, h.hi/phycon.hden, h.hii/phycon.hden, hmi.htwo/phycon.hden, hmi.h2plus/phycon.hden, hmi.h3plus/phycon.hden, hmi.hminus/phycon.hden ); } else if( strcmp(pnunit.chPunch[i],"HYDi") == 0 ) { if( strcmp(chTime,"LAST") != 0 ) { /* punch hydrogen ionization * this will be total decays to ground */ RateInter = 0.; stage = hydro.hcolnc[0][0]*phycon.eden*hydro.hn[0][IP1S]; fref = HPhoto.hgamnc[0][IP1S]*hydro.hn[0][IP1S]; fout = hydro.hn[0][IP1S]; for( ion=IP2S; ion <= nhlevel; ion++ ) { /* this is total decays to ground */ RateInter += hydro.HRadNet[0][IP1S][ion]; /* total photo from all levels */ fref += HPhoto.hgamnc[0][ion]*hydro.hn[0][ion]; /* total col ion from all levels */ stage += hydro.hcolnc[0][ion]*phycon.eden* hydro.hn[0][ion]; fout += hydro.hn[0][ion]; } fprintf( pnunit.ipPnunit[i], "hion\t%4ld\t%10.2e\t%10.2e", ZoneCnt.nzone, HPhoto.hgamnc[0][IP1S], hcaseb.HRecEffec[0] ); for(j=0; j < LIMELM; j++) { fprintf( pnunit.ipPnunit[i], "\t%10.2e", hcaseb.HRecRad[j] ); } fprintf( pnunit.ipPnunit[i], "\t%10.2e\t%10.2e\t%10.2e\t%10.2e\t%10.2e\t%10.2e\t%10.2e\t%10.2e\n", h.hii/h.hi, HPhoto.hgamnc[0][IP1S]/(phycon.eden*hcaseb.HRecEffec[0]), hydro.hrec[0][1][ipRecEsc], RateInter, fref/MAX2(1e-37,fout), stage/MAX2(1e-37,fout), HPhoto.hgamnc[0][IP1S]*h.hi/(phycon.eden* h.hii) , Secondaries.csupra); GammaPrt(nh.ipHn[0][0],rfield.nflux,OpacPoint.ipHOpac[0][IP1S], pnunit.ipPnunit[i],HPhoto.hgamnc[0][IP1S],HPhoto.hgamnc[0][IP1S]* 0.05); } } else if( strcmp(pnunit.chPunch[i],"HYDp") == 0 ) { /* punch hydrogen populations */ fprintf( pnunit.ipPnunit[i], " %10.2e%10.2e%10.2e", radius.depth, h.hi, h.hii ); for( j=IP1S; j <= 6; j++ ) { fprintf( pnunit.ipPnunit[i], "%10.2e", hydro.hn[0][j] ); } fprintf( pnunit.ipPnunit[i], "%10.2e%10.2e\n", hydro.hn[0][IP2S], hydro.hn[0][IP2P] ); } else if( strcmp(pnunit.chPunch[i],"IONI") == 0 ) { if( strcmp(chTime,"LAST") == 0 ) { /* punch mean ionization distribution */ pmeani( 'i', pnunit.ipPnunit[i] ); } } else if( strcmp(pnunit.chPunch[i]," IP ") == 0 ) { if( strcmp(chTime,"LAST") == 0 ) { /* punch valence shell ip's */ for( nelem=0; nelem < LIMELM; nelem++ ) { for( n=0; n <= 15; n += 15 ) { if( n + 1 <= nelem+1 ) { fprintf( pnunit.ipPnunit[i], " \n" ); limit = MIN2(n+15,nelem+1); fprintf( pnunit.ipPnunit[i], "%2.2s", elementnames.chElementSym[nelem]); for( ion=n + 1; ion <= limit; ion++ ) { fprintf( pnunit.ipPnunit[i], "%4ld", ion ); } fprintf( pnunit.ipPnunit[i], "\n" ); } if( n + 1 <= nelem+1 ) { for( ips=0; ips < nsShellsCom.nsShells[n][MIN2(n+15,nelem+1)-1]; ips++ ) { limit = MIN2(n+15,nelem+1); for( j=n; j < limit; j++ ) { sum = PH1COM.PH1[ips][nelem+1-j][nelem][0]; if( sum < 2. ) { strcpy( chEner[j], " " ); } else if( sum < 10. ) { sprintf( chEner[j], "%6.3f", sum ); } else if( sum < 100. ) { sprintf( chEner[j], "%6.2f", sum ); } else if( sum < 1000. ) { sprintf( chEner[j], "%6.1f", sum ); } else { sprintf( chEner[j], "%6ld", (long)(sum) ); } } limit = MIN2(n+15,nelem+1); fprintf( pnunit.ipPnunit[i], "%2.2s", Shell.chShell[ips] ); for( j=n; j < limit; j++ ) { fprintf( pnunit.ipPnunit[i], "%6.6s", chEner[j] ); } fprintf( pnunit.ipPnunit[i], "\n" ); } } } } } } else if( strcmp(pnunit.chPunch[i],"LINC") == 0 ) { /* punch line cumulative */ /* lgOK not used, placdholder */ if( strcmp(chTime,"LAST") != 0 ) { /* not used for anything here, but should not pass baloney*/ lgOK = TRUE; punlin(pnunit.ipPnunit[i],"PUNC",lgOK); } } else if( strcmp(pnunit.chPunch[i],"LINL") == 0 ) { /* punch line labels, then stop */ pllabels(pnunit.ipPnunit[i]); puts( "[Stop in DoPunch]" ); fprintf( ioQQQ, " Cloudy ends:\n" ); exit(1); } else if( strcmp(pnunit.chPunch[i],"LIND") == 0 ) { /* punch line data, then stop */ pldata(pnunit.ipPnunit[i]); puts( "[Normal stop in DoPunch]" ); fprintf( ioQQQ, " Cloudy ends:\n" ); exit(1); } else if( strcmp(pnunit.chPunch[i],"LINS") == 0 ) { /* punch line structure */ if( strcmp(chTime,"LAST") != 0 ) { /* not acutally used, but should not pass baloney*/ lgOK = TRUE; punlin(pnunit.ipPnunit[i],"PUNS",lgOK); } } else if( strcmp(pnunit.chPunch[i],"LINA") == 0 ) { if( strcmp(chTime,"LAST") == 0 ) { /* punch out all lines with energies */ for( j=0; j < LineSave.nsum; j++ ) { if( LineSv[j].xLineEnergy > 0. && LineSv[j].sumlin > 0. ) { fprintf( pnunit.ipPnunit[i], "%12.5e%8.3f %4.4s\n", LineSv[j].xLineEnergy, log10(LineSv[j].sumlin) + fourpi.pirsq, LineSv[j].chALab ); } } } } else if( strcmp(pnunit.chPunch[i],"LINI") == 0 ) { if( strcmp(chTime,"LAST") == 0 && (ZoneCnt.nzone/levery.LinEvery)*levery.LinEvery != ZoneCnt.nzone ) { /* this is the last zone * punch line intensities - but do not do last zone twice */ plinin(pnunit.ipPnunit[i]); } else if( strcmp(chTime,"LAST") != 0 ) { /* following so we only punch first zone if LinEvery reset */ if( (levery.lgLinEvery && ZoneCnt.nzone == 1) || (ZoneCnt.nzone/levery.LinEvery)*levery.LinEvery == ZoneCnt.nzone ) { /* this is middle of calculation * punch line intensities */ plinin(pnunit.ipPnunit[i]); } } } /* punch Lya - some details about Lya */ else if( strcmp(pnunit.chPunch[i],"LYMA") == 0 ) { pop = (hydro.hn[0][IP2P]/6.)/(hydro.hn[0][IP1S]/2.); if( pop > 0. ) { texc = -1.183e5/log(pop); } else { texc = 0.; } fprintf( pnunit.ipPnunit[i], " %.2e %.2e %.2e %.2e %.2e %.2e %.2e %.2e %.2e\n", radius.depth, HydroLines[0][IP2P][IP1S].TauIn, pop, texc, phycon.te, texc/phycon.te , HydroLines[0][IP2P][IP1S].Pdest, opac.opac[HydroLines[0][IP2P][IP1S].ipCont], opac.albedo[HydroLines[0][IP2P][IP1S].ipCont] ); } else if( strcmp(pnunit.chPunch[i],"MAP ") == 0 ) { /* do the map now if we are at the zone, or if this * is the LAST call to this routine and map not done yet */ if( (MapPar.lgMapDone==FALSE ) && (ZoneCnt.nzone == MapPar.MapZone || strcmp(chTime,"LAST") == 0 ) ) { lgTlkSav = called.lgTalk; called.lgTalk = TRUE; isav = ioQQQ; ioQQQ = pnunit.ipPnunit[i]; punt(" map"); ioQQQ = isav; called.lgTalk = lgTlkSav; } } else if( strcmp(pnunit.chPunch[i],"MOLE") == 0 ) { fprintf( pnunit.ipPnunit[i], " %10.2e%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e\n", radius.depth, hmi.htwo/phycon.hden, hevmolec.hevmol[IPCO-1]/ phycon.hden, hevmolec.hevmol[IPH2O-1]/phycon.hden, hevmolec.hevmol[IPOH-1]/phycon.hden, hevmolec.hevmol[IPCH-1]/phycon.hden, hevmolec.hevmol[IPOTWO-1]/phycon.hden ); } else if( strcmp(pnunit.chPunch[i],"OPAC") == 0 ) { /* punch opacity - routine will parse which type of opacity punch to do */ if( strcmp(chTime,"LAST") == 0 ) PunOpac(pnunit.ipPnunit[i],i); } else if( strcmp(pnunit.chPunch[i],"OPTI") == 0 ) { if( strcmp(chTime,"LAST") != 0 ) { for( j=0; j < rfield.nflux; j++ ) { fprintf( pnunit.ipPnunit[i], "%12.4e%12.4e%12.4e%12.4e\n", rfield.anu[j], opac.TauTotal[0][j], opac.tauabs[0][j], opac.tausct[0][j] ); } } } else if( strcmp(pnunit.chPunch[i]," OTS") == 0 ) { ConMax = 0.; xLinMax = 0.; opConSum = 0.; opLinSum = 0.; ipLinMax = 1; ipConMax = 1; for( j=0; j < rfield.nflux; j++ ) { opConSum += rfield.otscon[j]*opac.opac[j]; opLinSum += rfield.otslin[j]*opac.opac[j]; if( rfield.otslin[j]*opac.opac[j] > xLinMax ) { xLinMax = rfield.otslin[j]*opac.opac[j]; ipLinMax = j+1; } if( rfield.otscon[j]*opac.opac[j] > ConMax ) { ConMax = rfield.otscon[j]*opac.opac[j]; ipConMax = j+1; } } fprintf( pnunit.ipPnunit[i], "tot con lin=%10.2e%10.2e lin=%4.4s%11.4e%10.2e con=%4.4s%11.4e%10.2e\n", opConSum, opLinSum, LineLabl.chLineLabel[ipLinMax-1] , rfield.anu[ipLinMax-1], xLinMax, LineLabl.chContLabel[ipConMax-1] , rfield.anu[ipConMax-1], ConMax ); } else if( strcmp(pnunit.chPunch[i],"OVER") == 0 ) { /* overview of model results, * depth, te, hden, eden, ion fracs H, He, c, O */ if( strcmp(chTime,"LAST") != 0 ) { fprintf( pnunit.ipPnunit[i],PrintEfmt("%10.3e", radius.depth )); fprintf( pnunit.ipPnunit[i], "%6.3f%8.3f%7.3f%7.3f%6.3f%6.3f", log10(phycon.te), log10(heat.htot), log10(phycon.hden), log10(phycon.eden), -log10(MAX2(5e-10, h.hi/phycon.hden)), -log10(MAX2(5e-10,h.hii/ phycon.hden)) ); for( j=1; j <= 3; j++ ) { fprintf( pnunit.ipPnunit[i], "%6.3f", -log10(MAX2(5e-10,xIonFracs[j][1]/ xIonFracs[0][1])) ); } for( j=1; j <= 4; j++ ) { fprintf( pnunit.ipPnunit[i], "%6.3f", -log10(MAX2(5e-10,xIonFracs[j][5]/ xIonFracs[0][5])) ); } for( j=1; j <= 6; j++ ) { fprintf( pnunit.ipPnunit[i], "%6.3f", -log10(MAX2(5e-10,xIonFracs[j][7]/ xIonFracs[0][7])) ); } fprintf( pnunit.ipPnunit[i], "\n" ); } } else if( strcmp(pnunit.chPunch[i]," PDR") == 0 ) { if( strcmp(chTime,"LAST") != 0 ) { av = opac.TauTotal[0][ipvfil.ipVFilter-1]*1.08574; fprintf( pnunit.ipPnunit[i], "%8.2e%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e\n", radius.depth, coldenCom.colden[IPCOLUMNDENSITY-1], av, phycon.te, h.hi/phycon.hden, hmi.htwo/phycon.hden, xIonFracs[1][5]/xIonFracs[0][5], hevmolec.hevmol[IPCO-1]/xIonFracs[0][5], hevmolec.hevmol[IPH2O-1]/xIonFracs[0][7] ); } } else if( strcmp(pnunit.chPunch[i],"PHYS") == 0 ) { if( strcmp(chTime,"LAST") != 0 ) { /* punch physical conditions */ fprintf( pnunit.ipPnunit[i], "%14.6e%12.4e%11.3e%11.3e%11.3e%11.3e\n", radius.depth, phycon.te, phycon.hden, phycon.eden, heat.htot, wind.AccelTot ); } } else if( strcmp(pnunit.chPunch[i],"PRES") == 0 ) { if( strcmp(chTime,"LAST") != 0 ) { fprintf( pnunit.ipPnunit[i], " %10.2e%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e %10.2e %10.2e %c\n", radius.depth, pressure.pnow, pressure.PressureInit + pressure.PresInteg, pressure.PressureInit, pressure.GasPres, pressure.RamPres, pressure.RadPres, /* subtract continuum rad pres which has already been added on */ pressure.PresInteg - pressure.pinzon, wind.windv/1e5, TorF(conv.lgConvPres) ); } } else if( strcmp(pnunit.chPunch[i],"RADI") == 0 ) { if( strcmp(chTime,"LAST") != 0 ) { fprintf( pnunit.ipPnunit[i], " %4ld%12.4e%12.4e%12.4e\n", ZoneCnt.nzone, radius.Radius, radius.depth, radius.drad ); } } else if( strcmp(pnunit.chPunch[i],"RESU") == 0 ) { /* punch results of the calculation */ if( strcmp(chTime,"LAST") == 0 ) result(pnunit.ipPnunit[i]); } else if( strcmp(pnunit.chPunch[i],"SOUS") == 0 ) { /* full spectrum of continuum source function at 1 depth * command was "punch source spectrum" */ if( strcmp(chTime,"LAST") != 0 ) { limit = MIN2(maxcor.ipMaxBolt,rfield.nflux); for( j=0; j < limit; j++ ) { fprintf( pnunit.ipPnunit[i], "%12.4e%12.4e%12.4e%12.4e%12.4e\n", rfield.anu[j], rfield.DiffLocal[j]/rfield.widflx[j], opac.opac[j], rfield.DiffLocal[j]/rfield.widflx[j]/MAX2(1e-35,opac.opac[j]), rfield.DiffLocal[j]/rfield.widflx[j]/MAX2(1e-35,opac.opac[j])/plankf(j) ); } } } else if( strcmp(pnunit.chPunch[i],"SOUD") == 0 ) { /* parts of continuum source function vs depth * command was source depth */ j = nh.ipHn[0][IP1S] + 2; fprintf( pnunit.ipPnunit[i], "%12.4e%12.4e%12.4e%12.4e\n", deptau.depabs[j-1], rfield.DiffLocal[j-1]/rfield.widflx[j-1]/MAX2(1e-35,opac.opac[j-1]), rfield.otscon[nh.ipHn[0][IP1S]-1], rfield.otscon[nh.ipHn[0][0]-1]/opac.opac[nh.ipHn[0][IP1S]-1] ); } /* this is punch special option */ else if( strcmp(pnunit.chPunch[i],"SPEC") == 0 ) { PunSpec(pnunit.ipPnunit[i],chTime); } else if( strcmp(pnunit.chPunch[i],"TEGR") == 0 ) { /* punch history of heating and cooling */ if( strcmp(chTime,"LAST") != 0 ) { fprintf( pnunit.ipPnunit[i], " " ); for( j=0; j < NGRID; j++ ) { fprintf( pnunit.ipPnunit[i], "%4ld%11.3e%11.3e%11.3e\n", HeatCool.nZonGrid[j], HeatCool.TeGrid[j], HeatCool.HtGrid[j], HeatCool.ClGrid[j] ); } } } else if( strcmp(pnunit.chPunch[i],"TEMP") == 0 ) { /* temperature and its derivative */ fprintf( pnunit.ipPnunit[i], " %4ld%12.4e%12.4e\n", ZoneCnt.nzone, phycon.te, (phycon.te - tlastCom.tlast)/ radius.drad ); } else if( strcmp(pnunit.chPunch[i],"TPRE") == 0 ) { /* temperature and its predictors, turned on with punch tprid */ fprintf( pnunit.ipPnunit[i], "%5ld %11.4e %11.4e %11.4e %g\n", ZoneCnt.nzone, tpred.TeInit, tpred.TeProp, phycon.te, (tpred.TeProp- phycon.te)/phycon.te ); } else if( strcmp(pnunit.chPunch[i],"VERN") == 0 ) { if( (strcmp(chTime,"LAST") == 0) && FeII.lgFeIION) { /*FeIILines("pun",norm.ScaleNormLine/LineSv[norm.ipNormWavL].sumlin);*/ /*FeIIPunchLines(FeII.ioVerner);*/ FeIIPunchLines( pnunit.ipPnunit[i] ); } } else if( strcmp(pnunit.chPunch[i],"WIND") == 0 ) { /* wind velocity, radiative accel, and ratio total to e- accel */ fprintf( pnunit.ipPnunit[i], "%12.4e%12.4e%12.4e%12.4e%12.4e\n", radius.Radius, radius.depth, wind.windv, wind.AccelTot, wind.fmul ); } else { /* this is insanity, internal flag set in getpunch not seen here */ fprintf( ioQQQ, " DoPunch does not recognize flag %4.4s set by GETPUNCH. This is impossible.\n", pnunit.chPunch[i] ); ShowMe(); puts( "[Stop in DoPunch]" ); exit(1); } /* at start of loop, broke to here if not last iteration, * and punch last was set * throw special line if punching every iteration */ if( strcmp(chTime,"LAST") == 0 && !IterCnt.lgLastIt ) { fprintf( pnunit.ipPnunit[i], " ###########################\n" ); } } } # ifdef DEBUG_FUN fputs( " <->DoPunch()\n", debug_fp ); # endif return; }