/*PrtComment analyze model, generating comments on its features */ #include #include "cddefines.h" #include "cddrive.h" #include "showme.h" #include "physconst.h" #include "nhe1lvl.h" #include "grainvar.h" #include "taulines.h" #include "hydrogenic.h" #include "linesave.h" #include "fe2neg.h" #include "lamase.h" #include "lyaext.h" #include "trace.h" #include "mappar.h" #include "elecden.h" #include "sdrmin.h" #include "dtmase.h" #include "insane.h" #include "converge.h" #include "elmton.h" #include "insanity.h" #include "deplon.h" #include "opac.h" #include "powi.h" #include "zonecnt.h" #include "caseb.h" #include "hphoto.h" #include "plasnu.h" #include "hionrad.h" #include "testit.h" #include "elementnames.h" #include "ca2mly.h" #include "broke.h" #include "tseton.h" #include "badstp.h" #include "totlum.h" #include "turhet.h" #include "cmfrac.h" #include "jmin.h" #include "difher.h" #include "nhe1.h" #include "pressure.h" #include "dphoht.h" #include "codish.h" #include "h2dish.h" #include "con0.h" #include "negoi.h" #include "abnslr.h" #include "comneg.h" #include "negdrg.h" #include "he.h" #include "rfield.h" #include "occmax.h" #include "h2max.h" #include "nustbl.h" #include "colden.h" #include "phycon.h" #include "hlmax.h" #include "timesc.h" #include "tenden.h" #include "stimax.h" #include "hextra.h" #include "hmihet.h" #include "nh.h" #include "tcool.h" #include "heat.h" #include "compton.h" #include "tehigh.h" #include "he1tau.h" #include "dalerr.h" #include "radius.h" #include "the1lo.h" #include "failz.h" #include "itercnt.h" #include "taucon.h" #include "sphere.h" #include "drminu.h" #include "fudgec.h" #include "lyamas.h" #include "hemase.h" #include "edgtcm.h" #include "physok.h" #include "occ1hi.h" #include "bndsok.h" #include "frcind.h" #include "called.h" #include "habing.h" #include "reason.h" #include "warngs.h" #include "itagan.h" #include "cyclot.h" #include "he3pht.h" #include "tff.h" #include "filfac.h" #include "wind.h" #include "geopp.h" #include "h2ozer.h" #include "secondaries.h" #include "struc.h" #include "max2pht.h" #include "heatlmax.h" #include "extramax.h" #include "taddlya.h" #include "double.h" #include "turb.h" #include "chmax.h" #include "mgexc.h" #include "o3exc.h" #include "betaver.h" #include "noconverge.h" #include "autoit.h" #include "input.h" #include "forcte.h" #include "charexc.h" #include "chgetlbl.h" #include "totlin.h" #include "badprt.h" #include "warnings.h" #include "outsum.h" #include "chkcaheps.h" #include "chlinelbl.h" #include "agecheck.h" #include "prtcomment.h" void PrtComment(void) { char chLbl[11], chLine[133]; int lgHe1thk, lgHe2thk, lgThick; long int _o, i, imas, ipLo , ipHi , ipZ, isav, j, nc, nd, nline, nn, nneg, ns, nw; double BigJump, absint , aj, alpha, big, bigm, collin, comfrc, differ, error, flur, freqn, rate, ratio, reclin, rel, small, tauneg, ts , HBeta, relfl , relhm, fedest, GetHeat, SumNeg, c4363, t4363, outin, outout, totwid, outtot; double VolComputed , VolExpected , ConComputed , ConExpected; # ifdef DEBUG_FUN fputs( "<+>PrtComment()\n", debug_fp ); # endif if( trace.lgTrace ) { fprintf( ioQQQ, " PrtComment called.\n" ); } /* * enter all comments cautions warnings and surprises into a large * stack of statements * at end of this routine is master call to printing routines */ itagan.lgIterAgain = FALSE; /* initialize the line saver */ wcnint(); if( BetaVer.nBetaVer > 0 ) { sprintf( chLine, " !This is beta test version%4ld and is intended for testing only.", BetaVer.nBetaVer ); bangin(chLine); } /* this flag set by call to fixit() routine, * and show that parts of the code needs repair. */ if( broke.lgFixit ) { sprintf( chLine, " !The code needs to be fixed." ); bangin(chLine); } /* say why calculation stopped */ if( badstp.lgBadStop ) { /* this case stop probably was not intended */ sprintf( warnings.chRgcln[0], " W-Calculation stopped because %21.21s Iteration%3ld of%3ld", reason.chReason, IterCnt.iter, IterCnt.itermx + 1 ); sprintf( chLine, " W-Calculation stopped because %21.21s", reason.chReason ); warnin(chLine); sprintf( chLine, " W-This was not intended." ); warnin(chLine); } else { if( (autoit.lgAutoIt && IterCnt.iter != IterCnt.itermx + 1) && IterCnt.iter > 2 ) { sprintf( warnings.chRgcln[0], " Calculation stopped because %21.21s Iteration%3ld of%3ld, not converged due to %10.10s", reason.chReason, IterCnt.iter, IterCnt.itermx + 1, NoConverge.chNotConverged ); } else { sprintf( warnings.chRgcln[0], " Calculation stopped because %21.21s Iteration%3ld of%3ld", reason.chReason, IterCnt.iter, IterCnt.itermx + 1 ); } } /* check whether stopped because default number of zones hit, * not intended?? */ if( (!ZoneCnt.lgZoneSet) && ZoneCnt.lgZoneTrp ) { badstp.lgBadStop = TRUE; sprintf( chLine, " W-Calculation stopped because default number of zones reached. Was this intended???" ); warnin(chLine); sprintf( chLine, " W-Default limit can be increased while retaining this check with the SET NEND command." ); warnin(chLine); } /* check whether stopped because zones too thin - only happens when glob set * and depth + dr = depth * not intended */ if( drminu.lgDrMinUsed || radius.lgdR2Small ) { badstp.lgBadStop = TRUE; sprintf( chLine, " W-Calculation stopped zone thickness became too thin. This was not intended." ); warnin(chLine); sprintf( chLine, " W-The most likely reason was an uncontrolled oscillation." ); warnin(chLine); ShowMe(); } if( radius.lgdR2Small ) { sprintf( chLine, " W-This happened because the globule scale became very small relative to the depth." ); warnin(chLine); sprintf( chLine, " W-This problem is described in Hazy." ); warnin(chLine); } /* did insanity occur? */ if( insanity.lgInsane ) { badstp.lgBadStop = TRUE; sprintf( chLine, " W-Insanity occurred." ); warnin(chLine); ShowMe(); warnin(chLine); } /* possible that last zone does not have stored temp - if so * then make it up - this happens for some bad stops */ if( struc.testr[ZoneCnt.nzone-1] == 0. ) struc.testr[ZoneCnt.nzone-1] = struc.testr[ZoneCnt.nzone-2]; if( struc.ednstr[ZoneCnt.nzone-1] == 0. ) struc.ednstr[ZoneCnt.nzone-1] = struc.ednstr[ZoneCnt.nzone-2]; /* give indication of geometry */ rel = radius.depth/radius.rinner; if( rel < 0.1 ) { sprintf( warnings.chRgcln[1], " The geometry is plane-parallel." ); GeoPP.lgGeoPP = TRUE; } else if( rel >= 0.1 && rel < 3. ) { GeoPP.lgGeoPP = FALSE; sprintf( warnings.chRgcln[1], " The geometry is a thick shell." ); } else { GeoPP.lgGeoPP = FALSE; sprintf( warnings.chRgcln[1], " The geometry is spherical." ); } /* levels of warnings: Warning (possibly major problems) * Caution (not likely to invalidate the results) * [nothing] interesting note * [ !] surprise, but not disaster */ if( fudgec.nfudge > 0 ) { sprintf( chLine, " W-Fudge factors in place. Why?" ); warnin(chLine); } /* this flag set by call to broken(), * and show that the code is broken. */ if( broke.lgBroke ) { sprintf( chLine, " W-The code is broken." ); warnin(chLine); } /* incorrect electron density detected in radinc during the calculation */ if( ElecDen.lgEdenBad ) { if( ElecDen.nzEdenBad == ZoneCnt.nzone ) { sprintf( chLine, " C-The assumed electron density was incorrect for the last zone." ); caunin(chLine); sprintf( chLine, " C-Did a temperature discontinuity occur??" ); caunin(chLine); } else { sprintf( chLine, " W-The assumed electron density was incorrect during the calculation. This is bad." ); warnin(chLine); ShowMe(); warnin(chLine); } } if( called.lgBusted ) { sprintf( chLine, " W-The code returned with BUSTED set. Something REALLY went wrong!" ); warnin(chLine); } /* thermal map was done but results were not ok */ if( MapPar.lgMapDone && !MapPar.lgMapOK ) { sprintf( chLine, " W-A thermal map was done and problems were detected. Check map output." ); warnin(chLine); } if( phycon.hden > 1.1e13 ) { if( phycon.hden > 1e15 ) { sprintf( chLine, " C-Density greater than 10**15, heavy elements are very uncertain." ); caunin(chLine); } else { sprintf( chLine, " C-Density greater than 10**13" ); caunin(chLine); } } /* HBeta is used later in the code to check on line intensites */ if( !cdLine("Pump",4861,&relfl,&absint) ) { fprintf( ioQQQ, " Did not find flur H-beta\n" ); puts( "[Stop in PrtComment]" ); ShowMe(); exit(1); } /* now find total Hbeta */ if( !cdLine( "totl",4861,&HBeta,&absint) ) { fprintf( ioQQQ, " Did not find totl H-beta\n" ); puts( "[Stop in PrtComment]" ); ShowMe(); exit(1); } if( HBeta > 0. ) { flur = relfl/HBeta; if( flur > 0.01 ) { sprintf( chLine, " Continuum fluorescent production of H-beta was significant." ); notein(chLine); } } /* check if there were problems with initial wind velocity */ if( wind.windv0 > 0. && ((!wind.lgWindOK) || wind.windv < 1e6) ) { sprintf( chLine, " C-Wind velocity below sonic point; solution is not valid." ); caunin(chLine); } /* now confirm that mass flux here is correct */ if( wind.windv0 != 0. ) { rel = wind.emdot/(wind.windv*phycon.hden)/radius.r1r0sq; if( fabs(1.-rel)> 0.02 ) { ShowMe(); sprintf( chLine, " C-Wind mass flux error is %g%%",fabs(1.-rel)*100. ); caunin(chLine); } } /* check that we didn't overrun zone scale */ if( ZoneCnt.nzone >= NZLIM ) { sprintf( chLine, " C-Number of zones exceeded NZLIM - several structure quantities are not valid." ); caunin(chLine); sprintf( chLine, " C-Number of zones exceeded NZLIM - The emission line spectrum is ok though." ); caunin(chLine); } /* check whether energy is conserved * following is outward continuum */ outsum(&outtot,&outin,&outout); /* reclin is the sum of all recombination line energy */ reclin = totlin('r'); if( reclin == 0. ) { sprintf( chLine, " !Total recombination line energy is 0." ); bangin(chLine); } /* collin is the sum of all coolants */ collin = totlin('c'); if( collin == 0. ) { sprintf( chLine, " !Total cooling is zero." ); bangin(chLine); } /* TotalLumin is total incident photon luminosity, evaluated in setup * POWER is incremented in LINES and is the total heating integrated * over the computed volume * COLLIN is evaluated here, is sum of all coolants */ if( (collin + reclin) > totlum.TotalLumin && (!tseton.lgTSetOn) ) { if( (hextra.cryden == 0. && turhet.TurbHeat == 0.) && !Secondaries.lgCSetOn ) { sprintf( chLine, " W-Radiated luminosity (cool+rec=%8.2e+%10.2e) is greater than that in incident radiation (TotalLumin=%8.2e). Power=%8.2e", collin, reclin, totlum.TotalLumin, heat.power ); warnin(chLine); /* write same thing directly to output (above will be sorted later) */ fprintf( ioQQQ, " This calculation DID NOT CONSERVE ENERGY!\n" ); /* print out the offenders */ badprt(totlum.TotalLumin); sprintf( chLine, " W-Something is really wrong: the ratio of radiated to incident luminosity is %10.2e.", (collin + reclin)/totlum.TotalLumin ); warnin(chLine); /* this can be caused by the force te command */ if( forcte.ForceTemp > 0. ) { fprintf( ioQQQ, " This may have been caused by the FORCE TE command.\n" ); fprintf( ioQQQ, " Remove it and run again.\n" ); } else { ShowMe(); } warnin(chLine); } } /* warn if cosmic rays and bmag field both present */ if( hextra.cryden*cyclot.CycloCoolCoef > 0. ) { sprintf( chLine, " W-Magnetic field and cosmic rays are both present. This is not correctly treated." ); warnin(chLine); } /* PrtComment if test code is in place */ if( Testit.lgTestIt ) { sprintf( chLine, " !Test code is in place." ); bangin(chLine); } /* lgComUndr set to .true. if compton cooling rate underflows to 0 */ if( Compton.lgComUndr ) { sprintf( chLine, " !Compton cooling rate underflows to zero. Is this important?" ); bangin(chLine); } /* make note if input stream extended beyond col 80, extra char were ignored */ if( lgCardSavTooLong ) { sprintf( chLine, " !Some input lines were longer than 80 characters. Extra characters were ignored." ); bangin(chLine); } /* lgHionRad set to .true. if no hydrogen ionizing radiation */ if( HionRad.lgHionRad ) { sprintf( chLine, " !There is no Hydrogen-ionizing radiation. Was this intended?" ); bangin(chLine); } /* check whether certain zones were thermally unstable */ if( nustbl.nUnstable > 0 ) { sprintf( chLine, " Derivative of net cooling negative and so possibly thermally unstable in%4ld zones.", nustbl.nUnstable ); notein(chLine); } /* generate a bang if a large fraction of the zones were unstable */ if( ZoneCnt.nzone > 1 && (float)(nustbl.nUnstable)/(float)(ZoneCnt.nzone) > 0.25 ) { sprintf( chLine, " !A large fraction of the zones were possibly thermally unstable,%4ld out of%4ld", nustbl.nUnstable, ZoneCnt.nzone ); bangin(chLine); } /* comment if negative coolants were ever significant */ if( CHMax.CoolHeatMax > 0.2 ) { sprintf( chLine, " !Negative cooling reached %6.1f%% of the local heating, due to %4.4s%6ld", CHMax.CoolHeatMax*100., CHMax.chCoolHeatMax, CHMax.ipCoolHeatMax ); bangin(chLine); } else if( CHMax.CoolHeatMax > 0.05 ) { sprintf( chLine, " Negative cooling reached %6.1f%% of the local heating, due to %4.4s%6ld", CHMax.CoolHeatMax*100., CHMax.chCoolHeatMax, CHMax.ipCoolHeatMax ); notein(chLine); } /* check if charge transfer heating cooling was important */ if( CharExc.HCharHeatMax > 0.05 ) { sprintf( chLine, " !Charge transfer heating reached%10.2e%% of the local heating.", CharExc.HCharHeatMax*100. ); bangin(chLine); } else if( CharExc.HCharHeatMax > 0.005 ) { sprintf( chLine, " Charge transfer heating reached%10.2e%% of the local heating.", CharExc.HCharHeatMax*100. ); notein(chLine); } if( CharExc.HCharCoolMax > 0.05 ) { sprintf( chLine, " !Charge transfer cooling reached%10.2e%% of the local heating.", CharExc.HCharCoolMax*100. ); bangin(chLine); } else if( CharExc.HCharCoolMax > 0.005 ) { sprintf( chLine, " Charge transfer cooling reached%10.2e%% of the local heating.", CharExc.HCharCoolMax*100. ); notein(chLine); } /* did Lya mase? */ if( LaMase.xLaMase < 0. ) { if( LaMase.xLaMase < -0.01 ) { sprintf( chLine, " !H Lya mased!! Smallest optical depth was%10.2e", LaMase.xLaMase ); bangin(chLine); } else { sprintf( chLine, " H Lya mased. Smallest optical depth was%10.2e", LaMase.xLaMase ); notein(chLine); } } /* check whether photo from up level of Mg2 2798 ever important */ if( mgexc.xMg2Max > 0.1 ) { sprintf( chLine, " !Photoionization of upper level of Mg II 2798 reached %5.1f%% of the total Mg+ photo rate.", mgexc.xMg2Max*100. ); bangin(chLine); } else if( mgexc.xMg2Max > 0.01 ) { sprintf( chLine, " Photoionization of upper level of Mg II 2798 reached %5.1f%% of the total Mg+ photo rate.", mgexc.xMg2Max*100. ); notein(chLine); } /* check whether photo from up level of [OI] 6300 ever important */ if( o3exc.poimax > 0.1 ) { sprintf( chLine, " !Photoionization of upper levels of [OI] reached %5.1f%% of the total O destruction rate.", o3exc.poimax*100. ); bangin(chLine); } else if( o3exc.poimax > 0.01 ) { sprintf( chLine, " Photoionization of upper levels of [OI] reached %5.1f%% of the total O destruction rate.", o3exc.poimax*100. ); notein(chLine); } /* check whether photo from up level of [OIII] 5007 ever important */ if( (o3exc.poiii2Max + o3exc.poiii3Max) > 0.1 ) { sprintf( chLine, " !Photoionization of upper levels of [OIII] reached %5.1f%% of the total O++ photo rate.", (o3exc.poiii2Max + o3exc.poiii3Max)*100. ); bangin(chLine); } else if( (o3exc.poiii2Max + o3exc.poiii3Max) > 0.01 ) { sprintf( chLine, " Photoionization of upper levels of [OIII] reached %5.1f%% of the total O++ photo rate.", (o3exc.poiii2Max + o3exc.poiii3Max)*100. ); notein(chLine); } /* check whether photoionization of He 2trip S was important */ if( he3pht.he3photo > 0.01 ) { sprintf( chLine, " Photoionization of He 2TriS reached %5.1f%% of the total rate" " out at zone %li, %5.1f%% of that was Lya.", he3pht.he3photo*100, he3pht.nzone, he3pht.he3lya*100. ); notein(chLine); } /* frequency array may not have been defined for all energies */ if( !bndsok.lgMMok ) { sprintf( chLine, " C-Continuum not defined in extreme infrared - Compton scat, grain heating, not treated.?" ); caunin(chLine); } if( !bndsok.lgHPhtOK ) { sprintf( chLine, " C-Continuum not defined at photon energies which ionize excited states of H, important for H- and ff heating." ); caunin(chLine); } if( !bndsok.lgXRayOK ) { sprintf( chLine, " C-Continuum not defined at X-Ray energies - Compton scattering and Auger ionization wrong?" ); caunin(chLine); } if( !bndsok.lgGamrOK ) { sprintf( chLine, " C-Continuum not defined at gamma-ray energies - pair production and Compton scattering OK?" ); caunin(chLine); } if( con0.lgCon0 ) { sprintf( chLine, " C-Continuum zero at some energies." ); caunin(chLine); } if( con0.lgCoStarInterpolationCaution ) { sprintf( chLine , " C-DoCoStar interpolated between non-adjoining tracks, this may not be valid." ); caunin(chLine); } if( occ1hi.lgOcc1Hi ) { sprintf( chLine, " !The continuum occupation number at 1 Ryd is greater than unity." ); bangin(chLine); } /* this flag set true it set dr forced first zone to be too big */ if( sdrminCom.lgDR2Big ) { sprintf( chLine, " C-The thickness of the first zone was set larger than optimal by a SET DR command." ); caunin(chLine); sprintf( chLine, " C-Consider using the STOP THICKNESS command instead." ); caunin(chLine); } /* check whether non-col excit of 4363 was important */ if( !cdLine("TOTL",4363,&t4363,&absint) ) { fprintf( ioQQQ, " PrtComment could not find total OIII 4363 with cdLine.\n" ); puts( "[Stop in PrtComment]" ); ShowMe(); exit(1); } if( !cdLine("Coll",4363,&c4363,&absint) ) { fprintf( ioQQQ, " PrtComment could not find collis OIII 4363 with cdLine.\n" ); puts( "[Stop in PrtComment]" ); ShowMe(); exit(1); } /* only print this comment if 4363 is significant and density low */ if( HBeta > 1e-35 ) { if( t4363/HBeta > 1e-4 && phycon.hden < 1e7 ) { ratio = (t4363 - c4363)/t4363; if( ratio > 0.01 ) { sprintf( chLine, " !Non-collisional excitation of [OIII] 4363 reached %6.2f%% of the total.", ratio*100. ); bangin(chLine); } else if( ratio > 0.001 ) { sprintf( chLine, " Non-collisional excitation of [OIII] 4363 reached %6.2f%% of the total.", ratio*100. ); notein(chLine); } } } /* check for plasma shielding */ if( plasnu.lgPlasNu ) { sprintf( chLine, " !The largest plasma frequency was%10.2e Ryd. The incident continuum is set to 0 below this.", plasnu.plsfrqmax ); bangin(chLine); } if( occmaxCom.occmax > 0.1 ) { if( occmaxCom.occmnu > 1e-4 ) { sprintf( chLine, " !The largest continuum occupation number was %10.3e at %10.3e Ryd.", occmaxCom.occmax, occmaxCom.occmnu ); bangin(chLine); } else { /* not surprising if occ num bigger than 1 around 1e-5 Ryd, * since this is the case for 3K background */ sprintf( chLine, " The largest continuum occupation number was %10.3e at%10.3e Ryd.", occmaxCom.occmax, occmaxCom.occmnu ); notein(chLine); } } if( occmaxCom.occmax > 1e4 && occmaxCom.occ1nu > 0. ) { /* occ1nu is energy (ryd) where continuum occ num falls below 1 */ if( occmaxCom.occ1nu < 0.0912 ) { sprintf( chLine, " The continuum occupation number fell below 1 at %10.3e microns.", 0.0912/occmaxCom.occ1nu ); notein(chLine); } else if( occmaxCom.occ1nu < 1. ) { sprintf( chLine, " The continuum occupation number fell below 1 at %10.3e Angstroms.", 912./occmaxCom.occ1nu ); notein(chLine); } else { sprintf( chLine, " The continuum occupation number fell below 1 at %10.3e Ryd.", occmaxCom.occ1nu ); notein(chLine); } } if( occmaxCom.tbrmax > 1e3 ) { sprintf( chLine, " !The largest continuum brightness temperature was %10.3eK at %10.3e Ryd.", occmaxCom.tbrmax, occmaxCom.tbrmnu ); bangin(chLine); } if( occmaxCom.tbrmax > 1e4 ) { /* tbr4nu is energy (ryd) where continuum bright temp falls < 1e4 */ if( occmaxCom.tbr4nu < 0.0912 ) { sprintf( chLine, " The continuum brightness temperature fell below 10,000K at %10.3e microns.", 0.0912/occmaxCom.tbr4nu ); notein(chLine); } else if( occmaxCom.tbr4nu < 1. ) { sprintf( chLine, " The continuum brightness temperature fell below 10,000K at %10.3e Angstroms.", 912./occmaxCom.tbr4nu ); notein(chLine); } else { sprintf( chLine, " The continuum brightness temperature fell below 10,000K at %10.3e Ryd.", occmaxCom.tbr4nu ); notein(chLine); } } /* turbulence AND constant pressure do not make sense */ if( turb.TurbVel > 0. && strcmp(pressure.chCPres,"CPRE") == 0 ) { sprintf( chLine, " !Both constant pressure and turbulence makes no physical sense???" ); bangin(chLine); } /* filling factor AND constant pressure do not make sense */ if( filfac.FillFac < 1. && strcmp(pressure.chCPres,"CPRE") == 0 ) { sprintf( chLine, " !Both constant pressure and a filling factor makes no physical sense???" ); bangin(chLine); } /* grains and solar abundances do not make sense */ if( GrainVar.lgDustOn && abnslr.lgAbnSolar ) { sprintf( chLine, " !Grains are present, but the gas phase abundances were left at the solar default. This is not physical." ); bangin(chLine); } /* check if depletion command set but no grains, another silly thing to do */ if( deplon.lgDepln && !GrainVar.lgDustOn ) { sprintf( chLine, " !Grains are not present, but the gas phase abundances were depleted. This is not physical." ); bangin(chLine); } /* constant temperature greater than continuum energy density temperature */ if( tseton.lgTSetOn && forcte.ForceTemp*1.0001 < tenden.TEnerDen ) { sprintf( chLine, " C-The continuum energy density temperature (%g)" " is greater than the electron temperature (%g).", tenden.TEnerDen , forcte.ForceTemp); caunin(chLine); sprintf( chLine, " C-This is unphysical." ); caunin(chLine); } /* remark that grains not present but energy density was low */ if( !GrainVar.lgDustOn && tenden.TEnerDen < 800. ) { sprintf( chLine, " Grains were not present but might survive in this environment (energy density temperature was %10.2eK)", tenden.TEnerDen ); notein(chLine); } /* call routine that will check age of cloud */ AgeCheck(); /* check on CaH and H-epsilon overlapping * need to do this since need to refer to lines arrays */ chkCaHeps(&totwid); if( totwid > 121. ) { sprintf( chLine, " H-eps and CaH overlap." ); notein(chLine); } /* warning that something was turned off */ if( !physok.lgPhysOK ) { sprintf( chLine, " !A physical process has been disabled." ); bangin(chLine); } /* check on lifetimes of [OIII] against photoionization, only for low den */ if( phycon.hden < 1e8 ) { if( o3exc.r5007Max > 0.0263 ) { sprintf( chLine, " !Photoionization of upper level of [OIII] 5007 reached %10.2e%% of the radiative lifetime.", o3exc.r5007Max*100. ); bangin(chLine); } else if( o3exc.r5007Max > 0.0263/10. ) { sprintf( chLine, " Photoionization of upper level of [OIII] 5007 reached %10.2e%% of the radiative lifetime.", o3exc.r5007Max*100. ); notein(chLine); } if( o3exc.r4363Max > 1.78 ) { sprintf( chLine, " !Photoionization of upper level of [OIII] 4363 reached %10.2e%% of the radiative lifetime.", o3exc.r4363Max*100. ); bangin(chLine); } else if( o3exc.r4363Max > 1.78/10. ) { sprintf( chLine, " Photoionization of upper level of [OIII] 4363 reached %10.2e%% of the radiative lifetime.", o3exc.r4363Max*100. ); notein(chLine); } } /* check whether total heating and cooling matched * >>chng 97 mar 28, added GrossHeat, heat in terms normally heat-cool */ error = fabs(heat.power-tcool.totcol)/((heat.power + tcool.totcol + tcool.GrossHeat)/2.); if( tseton.lgTSetOn ) { if( error > 0.05 ) { sprintf( chLine, " !Heating - cooling mismatch =%5.1f%%. Caused by constant temperature assumption. ", error*100. ); bangin(chLine); } } else { if( error > 0.05 && error < 0.2 ) { sprintf( chLine, " C-Heating - cooling mismatch =%5.1f%%. Whats wrong?", error*100. ); caunin(chLine); } else if( error >= 0.2 ) { sprintf( chLine, " W-Heating - cooling mismatch =%10.2e Whats wrong????", error ); warnin(chLine); } } /* say if Ly-alpha photo of Ca+ excited levels was important */ if( ca2mly.Ca2RmLya > 0.01 ) { sprintf( chLine, " Photoionization of Ca+ 2D level by Ly-alpha reached %6.1f%% of the total rate out.", ca2mly.Ca2RmLya*100. ); notein(chLine); } /* check if Lya alpha ever hotter than gas */ if( LyaExT.nLyaHot > 0 ) { if( LyaExT.TLyaMax/LyaExT.TeLyaMax > 1.05 ) { sprintf( chLine, " !The excitation temp of Lya exceeded the electron temp, largest value was%10.2eK (gas temp there was%10.2eK, zone%4ld)", LyaExT.TLyaMax, LyaExT.TeLyaMax, LyaExT.nZTLaMax ); bangin(chLine); } } /* was line photoionization of hydrogen important */ hlmaxCom.hlmax[0] = 0.; isav = 0; for( i=2; i <= nhlevel; i++ ) { if( hlmaxCom.hlmax[i] > 0.001 ) { if( hlmaxCom.hlmax[i] > hlmaxCom.hlmax[1] ) { hlmaxCom.hlmax[1] = hlmaxCom.hlmax[i]; isav = i; } } } if( isav > 0 ) { sprintf( chLine, " Hydrogen self ionization by alpha transitions reached %6.1f%% of the destruction rate for level%3ld.", hlmaxCom.hlmax[1]*100., isav ); notein(chLine); } /* check if line absorption heating was important */ /* get all negative lines, check if line absorption significant heat source * this is used in "final" for energy budget print out */ if( !cdLine("Line",0,&SumNeg,&absint) ) { fprintf( ioQQQ, " did not get sumneg cdLine\n" ); puts( "[Stop in PrtComment]" ); ShowMe(); exit(1); } /* this is total heating */ if( !cdLine("TotH",0,&GetHeat,&absint) ) { fprintf( ioQQQ, " did not get GetHeat cdLine\n" ); puts( "[Stop in PrtComment]" ); ShowMe(); exit(1); } if( GetHeat > 0. ) { SumNeg /= GetHeat; if( SumNeg > 0.1 ) { sprintf( chLine, " !Line absorption heating reached%10.2e%% of the global heating.", SumNeg*100. ); bangin(chLine); } else if( SumNeg > 0.01 ) { sprintf( chLine, " Line absorption heating reached%10.2e%% of the global heating.", SumNeg*100. ); notein(chLine); } } /* this is check of extra lines added with g-bar */ if( ExtraMax.GBarMax > 0.1 ) { strcpy( chLbl, chLineLbl(&TauLine2[ExtraMax.ipMaxExtra-1]) ); sprintf( chLine, " !G-bar cooling lines reached%10.2e%% of the local cooling. Line=%10.10s", ExtraMax.GBarMax*100., chLbl ); bangin(chLine); } else if( ExtraMax.GBarMax > 0.01 ) { strcpy( chLbl, chLineLbl(&TauLine2[ExtraMax.ipMaxExtra-1]) ); sprintf( chLine, " G-bar cooling lines reached%10.2e%% of the local cooling. Line=%10.10s", ExtraMax.GBarMax*100., chLbl ); notein(chLine); } /* line absorption heating reached more than 10% of local heating? * HeatLineMax is largest heating(1,23)/htot */ if( heatlmax.HeatLineMax > 0.1 ) { if( heatlmax.levlmax == 1 ) { /* main block of lines */ strcpy( chLbl, chGetLbl(heatlmax.ipHeatlmax) ); } else if( heatlmax.levlmax == 2 ) { /* level 2 lines */ strcpy( chLbl, chLineLbl(&TauLine2[heatlmax.ipHeatlmax-1]) ); } else { fprintf( ioQQQ, " PrtComment has insane levlmax,=%5ld\n", heatlmax.levlmax ); } sprintf( chLine, " !Line absorption heating reached%10.2e%% of the local heating - largest by level%2ld line %10.10s", heatlmax.HeatLineMax*100., heatlmax.levlmax, chLbl ); bangin(chLine); } else if( heatlmax.HeatLineMax > 0.01 ) { sprintf( chLine, " Line absorption heating reached%10.2e%% of the local heating.", heatlmax.HeatLineMax*100. ); notein(chLine); } /* did Ly-alpha mase */ if( lyamas.lgHLyaMased ) { sprintf( chLine, " Ly-alpha mased." ); notein(chLine); } /* did HeII mase? */ if( hemase.lgHeMase ) { sprintf( chLine, " Helium continuum mased." ); notein(chLine); } /* check whether excited states of helium were optically thick */ lgHe1thk = FALSE; lgHe2thk = FALSE; for( i=1; i < 5; i++ ) { for( j=i + 1; j < 6; j++ ) { if( he1tauCOM.he1tau[i][j] > 1. ) lgHe1thk = TRUE; } } if( lgHe1thk ) { sprintf( chLine, " Some excited state HeI lines are thick." ); notein(chLine); } if( lgHe2thk ) { sprintf( chLine, " Some excited state HeII lines are thick." ); notein(chLine); } if( phycon.hden < 1e7 ) { /* check on IR fine structure lines */ lgThick = FALSE; tauneg = 0.; alpha = 0.; for( i=1; i <= nLevel1; i++ ) { /* define IR as anything longward of 1 micron */ if( TauLines[i].EnergyWN < 10000. ) { if( TauLines[i].TauIn > 1. ) { lgThick = TRUE; alpha = MAX2(alpha,TauLines[i].TauIn); } else if( TauLines[i].TauIn < tauneg ) { tauneg = tauneg; strcpy( chLbl, chLineLbl(&TauLines[i]) ); } } } /* now print results, were any fine structure lines optically thick? */ if( lgThick ) { sprintf( chLine, " !Some infrared fine structure lines are optically thick: largest tau was%10.2e", alpha ); bangin(chLine); } /* did any fine structure lines mase? */ if( tauneg < -0.01 ) { sprintf( chLine, " !Some fine structure lines mased:ion with smallest tau was%10.10s%10.2e", chLbl, tauneg ); bangin(chLine); } } /* were any other lines masing? */ tauneg = 0.; alpha = 0.; for( i=1; i <= nLevel1; i++ ) { /* define UV as anything shortward of 1 micron */ if( TauLines[i].EnergyWN >= 10000. ) { if( TauLines[i].TauIn < tauneg ) { tauneg = TauLines[i].TauIn; strcpy( chLbl, chLineLbl(&TauLines[i]) ); } } } /* did any lines mase? */ if( tauneg < -0.01 ) { sprintf( chLine, " !Some lines mased: most negative ion and tau were:%10.10s%10.2e", chLbl, tauneg ); bangin(chLine); } /* this is check that at least a second iteration was done with sphere static, * the test is orverridden with the (OK) option on the sphere static command, * which sets sphere.lgStaticNoIt TRUE */ if( (sphere.lgStatic && IterCnt.lgLastIt) && (IterCnt.iter == 1) && !sphere.lgStaticNoIt) { sprintf( chLine, " C-I must iterate when SPHERE STATIC is set." ); caunin(chLine); itagan.lgIterAgain = TRUE; } /* how important was induced two photon?? */ if( Max2Pht.xMax2Pht > 1. ) { sprintf( chLine, " !Rate of induced H 2-photon emission reached %10.2e s^-1", Max2Pht.xMax2Pht ); bangin(chLine); } else if( Max2Pht.xMax2Pht > 0.05 ) { sprintf( chLine, " Rate of induced H 2-photon emission reached %10.2e s^-1", Max2Pht.xMax2Pht ); notein(chLine); } /* how important was induced recombination? */ if( frcind.FracInd > 0.01 ) { sprintf( chLine, " Induced recombination was %5.1f%% of the total for H level%3ld", frcind.FracInd*100., frcind.ndclev ); notein(chLine); } if( frcind.fbul > 0.01 ) { sprintf( chLine, " Stimulated emission was%6.1f%% of the total for H transition%3ld -%3ld", frcind.fbul*100., frcind.nbul + 1, frcind.nbul ); notein(chLine); } /* check whether FeII destruction of La was important */ if( !cdLine("Fe 2",1216,&fedest,&absint) ) { fprintf( ioQQQ, " Did not find Fe ii Lya\n" ); puts( "[Stop in PrtComment]" ); ShowMe(); exit(1); } /* find total Lya for comparison */ if( !cdLine("TOTL",1216,&relhm,&absint) ) { fprintf( ioQQQ, " Did not find Lya\n" ); puts( "[Stop in PrtComment]" ); ShowMe(); exit(1); } if( relhm > 0. ) { ratio = fedest/(fedest + relhm)*100.; if( ratio > 10. ) { sprintf( chLine, " !FeII destruction of Ly-a removed %5.1f%% of the line.", ratio ); bangin(chLine); } else if( ratio > 0.1 ) { sprintf( chLine, " FeII destruction of Ly-a removed %5.1f%% of the line.", ratio ); notein(chLine); } } if( !cdLine("H-CT",6563,&relhm,&absint) ) { fprintf( ioQQQ, " Comment did not find H-CT H-alpha\n" ); puts( "[Stop in PrtComment]" ); ShowMe(); exit(1); } if( HBeta > 0. ) { if( relhm/HBeta > 0.01 ) { sprintf( chLine, " !Mutual neutralization production of H-alpha was significant." ); bangin(chLine); } } /* note about very high population in H n=2 rel to ground, set in hydrogenic */ if( hydro.lgHiPop2 ) { sprintf( chLine, " The population of H n=2 reached%10.2e relative to the ground state.", hydro.pop2mx ); notein(chLine); } /* check where diffuse emission error */ for( ipZ=0; ipZ < LIMELM; ipZ++ ) { if( elmton.lgElmtOn[ipZ] ) { if( difher.CaseBCheck[ipZ] > 1.5 ) { sprintf( chLine, " Ratio of computed diffuse emission to case B reached %g for element %li", difher.CaseBCheck[ipZ] , ipZ+1 ); notein(chLine); } } } if( tehigh.thist > 1e9 ) { if( tehigh.thist > 5e9 ) { sprintf( chLine, " W-Electrons were relativistic; High TE=%10.2e", tehigh.thist ); warnin(chLine); } else { sprintf( chLine, " C-Electrons were mildly relativistic; High TE=%10.2e", tehigh.thist ); caunin(chLine); } } rate = timesc.TimeErode*2e-26; if( rate > 1e-35 ) { /* 2E-26 is roughly cross section for photoerosion * see Boyd and Ferland ApJ Let 318, L21. */ ts = (1./rate)/3e7; if( ts < 1e3 ) { sprintf( chLine, " !Timescale-photoerosion of Fe=%.2e yr", ts ); bangin(chLine); } else if( ts < 1e9 ) { sprintf( chLine, " Timescale-photoerosion of Fe=%.2e yr", ts ); notein(chLine); } } /* check whether Compton heating was significant */ comfrc = Compton.comtot/heat.power; if( comfrc > 0.01 ) { sprintf( chLine, " Compton heating was %5.1f%% of the total.", comfrc*100. ); notein(chLine); } /* check on relative importance of induced compton heating */ if( comfrc > 0.01 && Compton.cinrat > 0.05 ) { sprintf( chLine, " !Induced Compton heating was %10.2e of the total Compton heating.", Compton.cinrat ); bangin(chLine); } /* check whether equilibrium timescales are short rel to Hubble time */ if( timesc.tcmptn > 5e17 ) { if( comfrc > 0.05 ) { sprintf( chLine, " C-Compton cooling is significant and the equilibrium timescale (%10.2es) is longer than the Hubble time.", timesc.tcmptn ); caunin(chLine); } else { sprintf( chLine, " Compton cooling equilibrium timescale (%10.2es) is longer than Hubble time.", timesc.tcmptn ); notein(chLine); } } if( timesc.ttherm > 5e17 ) { sprintf( chLine, " C-Thermal equilibrium timescale longer than Hubble time, =%10.2es", timesc.ttherm ); caunin(chLine); } /* check whether model large relative to Jeans length * DMEAN is mean density (gm per cc) * mean temp is weighted by mass density */ if( log10(radius.depth) > jmin.rjnmin ) { /* AJMIN is minimum Jeans mass, log in grams */ aj = pow(10.,jmin.ajmmin - log10(SOLAR_MASS)); if( strcmp(pressure.chCPres,"CPRE") == 0 ) { sprintf( chLine, " C-Cloud thicker than smallest Jeans length=%8.2ecm; stability problems? (smallest Jeans mass=%8.2eMo)", pow(10.,jmin.rjnmin), aj ); caunin(chLine); } else { sprintf( chLine, " Cloud thicker than smallest Jeans length=%8.2ecm; stability problems? (smallest Jeans mass=%8.2eMo)", pow(10.,jmin.rjnmin), aj ); notein(chLine); } } /* check whether grains too hot to survive */ for( nd=0; nd < NDUST; nd++ ) { if( GrainVar.lgDustOn1[nd] ) { if( GrainVar.TeGrainMax[nd] > GrainVar.Tsublimat[nd] ) { sprintf( chLine, " W-Maximum temperature of grain %10.10s was%10.2eK, above its sublimation temperature,%10.2eK.", chGrainVar.chDstLab[nd], GrainVar.TeGrainMax[nd], GrainVar.Tsublimat[nd] ); warnin(chLine); } else if( GrainVar.TeGrainMax[nd] > GrainVar.Tsublimat[nd]* 0.9 ) { sprintf( chLine, " C-Maximum temperature of grain %10.10s was%10.2eK, near its sublimation temperature,%10.2eK.", chGrainVar.chDstLab[nd], GrainVar.TeGrainMax[nd], GrainVar.Tsublimat[nd] ); caunin(chLine); } } } if( negdrg.lgNegGrnDrg ) { sprintf( chLine, " !Grain drag force <0." ); bangin(chLine); } /* is photoelectric heating of gas by photioniz of grains important */ if( dphoht.dphmax > 0.05 ) { sprintf( chLine, " Local grain-gas photoelectric heating rate reached %5.1f%% of the total.", dphoht.dphmax*100. ); notein(chLine); } if( dphoht.TotalDustHeat/heat.power > 0.01 ) { sprintf( chLine, " Global grain photoelectric heating of gas was%5.1f%% of the total.", dphoht.TotalDustHeat/heat.power*100. ); notein(chLine); if( dphoht.TotalDustHeat/heat.power > 0.25 ) { sprintf( chLine, " !Grain photoelectric heating is VERY important." ); bangin(chLine); } } /* grain-gas collis cooling of gas */ if( dphoht.dclmax > 0.05 ) { sprintf( chLine, " Local grain-gas cooling of gas rate reached %5.1f%% of the total.", dphoht.dclmax*100. ); notein(chLine); } /* check whether photodissociation of H_2^+ molecular ion was important */ if( dphoht.h2pmax > 0.10 ) { sprintf( chLine, " !The local H2+ photodissociation heating rate reached %5.1f%% of the total heating.", dphoht.h2pmax*100. ); bangin(chLine); } else if( dphoht.h2pmax > 0.01 ) { sprintf( chLine, " The local H2+ photodissociation heating rate reached %5.1f%% of the total heating.", dphoht.h2pmax*100. ); notein(chLine); } /* check whether photodissociation of molecular hydrogen (H2)was important */ if( h2dish.h2dfrc > 0.05 ) { sprintf( chLine, " The local H2 photodissociation heating rate reached %5.1f%% of the total.", h2dish.h2dfrc*100. ); notein(chLine); } if( h2dish.h2dtot/heat.power > 0.01 ) { sprintf( chLine, " Global H2 photodissociation heating of gas was%5.1f%% of the total.", h2dish.h2dtot/heat.power*100. ); notein(chLine); if( h2dish.h2dtot/heat.power > 0.25 ) { sprintf( chLine, " H2 photodissociation heating is VERY important." ); notein(chLine); } } /* check whether photodissociation of carbon monoxide (co) was important */ if( codish.codfrc > 0.25 ) { sprintf( chLine, " !Local CO photodissociation heating rate reached %5.1f%% of the total.", codish.codfrc*100. ); bangin(chLine); } else if( codish.codfrc > 0.05 ) { sprintf( chLine, " Local CO photodissociation heating rate reached %5.1f%% of the total.", codish.codfrc*100. ); notein(chLine); } if( codish.codtot/heat.power > 0.01 ) { sprintf( chLine, " Global CO photodissociation heating of gas was%5.1f%% of the total.", codish.codtot/heat.power*100. ); notein(chLine); if( codish.codtot/heat.power > 0.25 ) { sprintf( chLine, " CO photodissociation heating is VERY important." ); notein(chLine); } } if( edgtcm.lgEdnGTcm ) { sprintf( chLine, " Energy density of radiation field was greater than the Compton temperature. Is this physical?" ); notein(chLine); } /* was cooling due to induced recombination important? */ if( HPhoto.cintot/heat.power > 0.01 ) { sprintf( chLine, " Induced recombination cooling was%5.1f%% of the total.", HPhoto.cintot/heat.power*100. ); notein(chLine); } /* check whether free-free heating was significant */ if( hydro.FreeFreeTotHeat/heat.power > 0.1 ) { sprintf( chLine, " !Free-free heating was%5.1f%% of the total.", hydro.FreeFreeTotHeat/heat.power*100. ); bangin(chLine); } else if( hydro.FreeFreeTotHeat/heat.power > 0.01 ) { sprintf( chLine, " Free-free heating was%5.1f%% of the total.", hydro.FreeFreeTotHeat/heat.power*100. ); notein(chLine); } /* was heating due to H- absorption important? */ if( hmihetCom.hmitot/heat.power > 0.01 ) { sprintf( chLine, " H- absorption heating was%5.1f%% of the total.", hmihetCom.hmitot/heat.power*100. ); notein(chLine); } /* water desruction rate was zero */ if( h2ozer.lgH2Ozer ) { sprintf( chLine, " Water destruction rate zero." ); notein(chLine); } /* numerical instability in matrix inverson routine */ if( negoi.nNegOI > 0 ) { sprintf( chLine, " C-OI negative level populations%4ld times.", negoi.nNegOI ); caunin(chLine); } /* numerical instability in matrix inverson routine */ if( Fe2Neg.nFe2Neg > 0 ) { sprintf( chLine, " C-FeII negative level populations%4ld times.", Fe2Neg.nFe2Neg ); caunin(chLine); } /* check for negative optical depths, * optical depth in excited state helium lines */ small = 0.; imas = 0; isav = 0; i = 0; for( ipZ=0; ipZ 0 ) { sprintf( chLine, " !Some continuous optical depths <0. The lowest freq was%10.3e Ryd, and a total of%4ld", freqn, nneg ); bangin(chLine); sprintf( chLine, " !The smallest optical depth was%10.2e", tauneg ); bangin(chLine); } /* say if Balmer continuum optically thick */ if( opac.tauabs[0][nh.ipHn[0][2]-1] > 0.05 ) { sprintf( chLine, " The Balmer continuum optical depth was%10.2e.", opac.tauabs[0][nh.ipHn[0][2]-1] ); notein(chLine); } /* was correction for stimulated emission significant? */ if( stimaxCom.stimax[0] > 0.02 && opac.tauabs[0][nh.ipHn[0][IP1S]-1] > 0.2 ) { sprintf( chLine, " The Lyman continuum stimulated emission correction to optical depths reached%10.2e.", stimaxCom.stimax[0] ); notein(chLine); } else if( stimaxCom.stimax[1] > 0.02 && opac.tauabs[0][nh.ipHn[0][2]-1] > 0.1 ) { sprintf( chLine, " The Balmer continuum stimulated emission correction to optical depths reached%10.2e.", stimaxCom.stimax[1] ); notein(chLine); } /* say if Paschen continuum optically thick */ if( opac.tauabs[0][nh.ipHn[0][3]-1] > 0.2 ) { sprintf( chLine, " The Paschen continuum optical depth was%10.2e.", opac.tauabs[0][nh.ipHn[0][3]-1] ); notein(chLine); } /* some comments about near IR total optical depth */ if( opac.tauabs[0][0] > 1. ) { sprintf( chLine, " The continuum optical depth at the lowest energy considered (%10.3e Ryd) was %10.3e.", rfield.anu[0], opac.tauabs[0][0] ); notein(chLine); } /* comment if optical depth to Rayleigh scattering is big * cs from VAL 76 */ if( coldenCom.colden[IPCHI-1]*7e-24 > 0.01 ) { sprintf( chLine, " The optical depth to Rayleigh scattering at 1300A is%10.2e", coldenCom.colden[IPCHI-1]*6.71e-24 ); notein(chLine); } if( coldenCom.colden[IPCH2PLS-1]*7e-18 > 0.01 ) { sprintf( chLine, " The optical depth to the H2+ molecular ion is%10.2e", coldenCom.colden[IPCH2PLS-1]*7e-18 ); notein(chLine); } /* warn if optically thick to H- absorption */ if( taucon.thmin > 0.01 ) { sprintf( chLine, " Optical depth to negative hydrogen ion is%10.2e", taucon.thmin ); notein(chLine); } /* say if 3-dody rec coef function outside range of validity * tripped if te/z**2 < 100 or approx 10**13: => effect >50% of radiative * other intergers defined in source for da */ if( dalerr.ifail > 0 && dalerr.ifail <= 10 ) { sprintf( chLine, " 3 body rec coef outside range%4ld", dalerr.ifail ); notein(chLine); } else if( dalerr.ifail > 10 ) { sprintf( chLine, " C-3 body rec coef outside range%4ld", dalerr.ifail ); caunin(chLine); } /* check whether energy density less than background */ if( tenden.TEnerDen < 2.6 ) { sprintf( chLine, " !Incident radiation field energy density is less than 2.7K. Add background with FIREBALL command." ); bangin(chLine); } /* incident radiation field is less than Habing ISM field */ if( habing.lgHabing ) { sprintf( chLine, " !The intensity of the incident radiation field is less than 10 times the Habing diffuse ISM field. Is this OK?" ); bangin(chLine); } /* check whether cosmic rays on, but model thick to them */ if( hextra.cryden > 0. && (coldenCom.colden[IPCHI-1]/10. + coldenCom.colden[IPCHII-1]) > 1e23 ) { sprintf( chLine, " C-Model is thick to cosmic rays, which are on." ); caunin(chLine); } /* was ionization rate less than cosmic ray ionization rate in ISM? */ if( hextra.cryden == 0. && HPhoto.hgamnc[0][IP1S] < 1e-17 ) { sprintf( chLine, " Ionization rate fell below background cosmic ray ionization rate. Should this be added too?" ); notein(chLine); } /* more than 10% of H is in the H2 molecule */ if( h2max.BiggestH2 > 0.1 ) { sprintf( chLine, " !A large amount of the hydrogen is molecular, the fraction reached %5.1f%%.", MIN2(1.,h2max.BiggestH2)*100. ); bangin(chLine); } else if( h2max.BiggestH2 > 1e-3 ) { sprintf( chLine, " The H2/H molecular fraction reached%10.2e%%.", h2max.BiggestH2*100. ); notein(chLine); } /* how much carbon is molecular? */ if( cmfrac.CarMolFrac > 0.1 ) { sprintf( chLine, " !A large amount of the carbon is molecular, fraction reached %5.1f%%", MIN2(1.,cmfrac.CarMolFrac)*100. ); bangin(chLine); } else if( cmfrac.CarMolFrac > 1e-3 ) { sprintf( chLine, " The carbon molecular fraction reached%10.2e%%", cmfrac.CarMolFrac*100. ); notein(chLine); } /* brems optical depth */ if( tffCom.tff > 0.09 ) { sprintf( chLine, " !The cloud is optically thick at optical wavelengths, extending to%10.3e Ryd =%10.3eA", tffCom.tff, RYDLAM/tffCom.tff ); bangin(chLine); } else if( tffCom.tff > 0.009 ) { sprintf( chLine, " The continuum of the computed structure may be optically thick in the near infrared." ); notein(chLine); } /* did model run away to very large radius? */ if( radius.Radius > 1e23 && radius.Radius/radius.rinner > 10. ) { sprintf( chLine, " Is an outer radius of%10.2e reasonable?", radius.Radius ); notein(chLine); } /* following set true in LineTauInc if maser capped at tau = -1 */ if( dtMase.lgMaserCapHit ) { sprintf( chLine, " Laser maser optical depths capped in LineTauInc." ); notein(chLine); } /* lgPradCap is true if radiation pressure was capped on first iteration * also check that this is a constant total pressure model */ if( (pressure.lgPradCap && (strcmp(pressure.chCPres,"CPRE") == 0)) && pressure.lgRadPresON ) { sprintf( chLine, " Radiation pressure kept below gas pressure on this iteration." ); notein(chLine); } if( pressure.RadBetaMax > 0.5 ) { if( pressure.ipPradMax == 0 ) { sprintf( chLine, " The ratio of radiation to gas pressure reached %10.2e. Caused by Lyman alpha.", pressure.RadBetaMax ); notein(chLine); } else { sprintf( chLine, " The ratio of radiation to gas pressure reached %10.2e. " "Caused by line number%3ld, label %s", pressure.RadBetaMax, pressure.ipPradMax,pressure.chLineRadPres ); notein(chLine); } } if( taucon.telec > 0.05 && taucon.telec <= 1.0 ) { sprintf( chLine, " The optical depth to electron scattering is%10.2e", taucon.telec ); notein(chLine); } else if( taucon.telec > 1.0 && taucon.telec < 5. ) { sprintf( chLine, " C-The model is optically thick to electron scattering; tau=%10.2e", taucon.telec ); caunin(chLine); } else if( taucon.telec >= 5. ) { sprintf( chLine, " W-The model is optically thick to electron scattering; tau=%10.2e", taucon.telec ); warnin(chLine); } /* optical depth to 21 cm */ if( TauLines[ipH21cm].TauIn > 0.5 ) { sprintf( chLine, " !The optical depth in the HI 21 cm line is%10.2e",TauLines[ipH21cm].TauIn ); bangin(chLine); } /* heavy element molecules had negative poppulations */ if( comneg.lgComNeg && cmfrac.CarMolFrac > 0.5 ) { sprintf( chLine, " !Populations of heavy element molecules were negative because of large CO population." ); bangin(chLine); } else if( comneg.lgComNeg ) { sprintf( chLine, " C-Populations of heavy element molecules were negative; CO, H2O, etc populations may not be correct." ); caunin(chLine); } /* check on optical depth convergence of all hydrogenic lines */ for( ipZ=0; ipZ < LIMELM; ipZ++ ) { if( elmton.lgElmtOn[ipZ] ) { if( HydroLines[ipZ][3][IP2P].TauIn > 0.2 ) { differ = fabs(1.-HydroLines[ipZ][3][IP2P].TauIn* double_.DoubleTau/HydroLines[ipZ][3][IP2P].TauTot)*100.; /* check whether H-alpha optical depth changed by much on last iteration * no tolerance can be finer than autocv, the tolerance on the * iterate to convergence comamnd. It is 15% */ if( ((IterCnt.lgLastIt && HydroLines[ipZ][3][IP2P].TauIn > 0.8) && differ > 20.) && wind.windv == 0. ) { sprintf( chLine, " C-This is the last iteration and %2s%2ld Bal(a) optical depth" " changed by%6.1f%% (was %10.2e). Try another iteration.", elementnames.chElementSym[ipZ], ipZ+1, differ, HydroLines[ipZ][3][IP2P].TauTot ); caunin(chLine); itagan.lgIterAgain = TRUE; } /* only check on Lya convergence if Balmer lines are thick */ if( HydroLines[ipZ][IP2P][IP1S].TauIn > 0. ) { differ = fabs(1.-HydroLines[ipZ][IP2P][IP1S].TauIn* double_.DoubleTau/HydroLines[ipZ][IP2P][IP1S].TauTot)*100.; /* check whether Lya optical depth changed on last iteration * no tolerance can be finer than autocv, the tolerance on the * iterate to convergence comamnd. It is 15% */ if( ((IterCnt.lgLastIt && HydroLines[ipZ][IP2P][IP1S].TauIn > 0.8) && differ > 25.) && wind.windv == 0. ) { sprintf( chLine, " C-This is the last iteration and %2s%2ld Ly(a) optical depth" " changed by%7.0f%% (was %10.2e). Try another iteration.", elementnames.chElementSym[ipZ], ipZ+1,differ, HydroLines[ipZ][IP2P][IP1S].TauTot ); caunin(chLine); itagan.lgIterAgain = TRUE; } } } } } if( (he1tauCOM.he1tau[0][1] < 1e3) && (he1tauCOM.he1tau[0][1] - TAddLya.TAddHeI) > 0.5 ) { differ = fabs(1.-he1tauCOM.he1tau[0][1]*double_.DoubleTau/ he1tauCOM.he1lim[0][1])*100.; /* check whether he 1 h-alpha optical depth converged */ if( ((IterCnt.lgLastIt && he1tauCOM.he1tau[0][1] > 0.8) && differ > 10.) && wind.windv == 0. ) { sprintf( chLine, " C-This is the last iteration and HeI(Lya)" " optical depth changed by%6.1f%% (was %10.2e). Try another iteration", differ, he1tauCOM.he1lim[0][1] ); caunin(chLine); itagan.lgIterAgain = TRUE; } } /* check whether sphere was set if dr/r large */ if( radius.Radius/radius.rinner > 2. && !sphere.lgSphere ) { sprintf( chLine, " C-R(out)/R(in)=%10.2e and SPHERE was not set.", radius.Radius/radius.rinner ); caunin(chLine); } /* check if thin in Lyman continuum, and assumed thick */ if( ((opac.tauabs[0][nh.ipHn[0][IP1S]-1] < 2. && opac.tauabs[1][nh.ipHn[0][IP1S]-1] > 3.) && IterCnt.lgLastIt) && (!caseb.lgCaseb) ) { sprintf( chLine, " C-The Lyman continuum is thin, and I assumed" " that it was thick. Try another iteration." ); caunin(chLine); itagan.lgIterAgain = TRUE; } /* check if thin in helium continua, but assumed to be thick */ if( IterCnt.lgLastIt ) { if( (opac.tauabs[0][nh.ipHn[1][IP1S]-1] < 2. && opac.tauabs[0][nh.ipHn[1][IP1S]-1] > 2.) && rfield.nflux > he.nheii ) { sprintf( chLine, " C-The HeII continuum is thin and I assumed that it was thick." " Try another iteration." ); caunin(chLine); itagan.lgIterAgain = TRUE; } if( (opac.tauabs[0][nhe1Com.nhe1[0]-1] < 2. && opac.tauabs[0][nhe1Com.nhe1[0]-1] > 2.) && rfield.nflux > he.nhei1 ) { sprintf( chLine, " C-The HeI continuum is thin and I assumed that it was thick." " Try another iteration." ); caunin(chLine); itagan.lgIterAgain = TRUE; } } /* check whether column density changed by much on this iteration */ if( IterCnt.iter > 1 ) { if( coldenCom.ColDenSav[IterCnt.iter-2] <= 0. ) { fprintf( ioQQQ, " ColDenSav is insane in PrtComment.\n" ); ShowMe(); puts( "[Stop in PrtComment]" ); ShowMe(); exit(1); } differ = fabs(1.-coldenCom.colden[IPCOLUMNDENSITY-1]/ coldenCom.ColDenSav[IterCnt.iter-2]); if( differ > 0.1 && differ <= 0.3 ) { sprintf( chLine, " The column density changed by %10.2e%% between this and previous iteration.", differ*100. ); notein(chLine); } else if( differ > 0.3 ) { if( IterCnt.lgLastIt ) { sprintf( chLine, " !The column density changed by %10.2e%% and this is the last iteration. What happened?", differ*100. ); bangin(chLine); } else { sprintf( chLine, " !The column density changed by %10.2e%% What happened?", differ*100. ); bangin(chLine); } } } /* say if rad pressure caused by la and la optical depth changed too much */ differ = fabs(1.-HydroLines[0][IP2P][IP1S].TauIn/ HydroLines[0][IP2P][IP1S].TauTot)*100.; if( IterCnt.lgLastIt && (pressure.RadBetaMax > 0.1) && (differ > 50.) && (pressure.ipPradMax == 1) && (pressure.lgRadPresON) && wind.windv == 0. ) { sprintf( chLine, " C-This is the last iteration, radiation pressure was significant, and the L-a optical depth changed by %7.2f%% (was%10.2e)", differ, HydroLines[0][IP2P][IP1S].TauTot ); caunin(chLine); } /* say if la rad pressure capped by thermalization length */ if( pressure.lgPradDen ) { sprintf( chLine, " Line radiation pressure capped by thermalization length." ); notein(chLine); } /* print te failures */ nline = MIN2(conv.nTeFail,10); if( conv.nTeFail != 0 ) { if( conv.failmx < 0.1 ) { _o = sprintf( chLine, " There were%4ld minor temperature failures. zones:", conv.nTeFail ); for( i=0; i < nline; i++ ) { _o += sprintf( chLine+_o, "%4ld", failz.ifailz[i] ); } notein(chLine); } else { _o = sprintf( chLine, " !There were%4ld temperature failures, and some were large. The largest was%6.1f%%. What happened?", conv.nTeFail, conv.failmx*100. ); bangin(chLine); _o = sprintf( chLine , " !The zones were" ); for( i=0; i < nline; i++ ) { _o += sprintf( chLine+_o, "%4ld", failz.ifailz[i] ); } bangin(chLine); if( struc.testr[0] > 8e4 && phycon.te < 6e5 ) { sprintf( chLine, " !I think they may have been caused by the change from hot to nebular gas phase. The physics of this is unclear." ); bangin(chLine); } } } /* check for temperature jumps */ BigJump = 0.; j = 0; for( i=1; i < ZoneCnt.nzone; i++ ) { big = fabs(1.-struc.testr[i-1]/struc.testr[i]); if( big > BigJump ) { j = i; BigJump = big; } } if( BigJump > 0.2 ) { /* this is a sanity check, but only do it if jump detected */ if( j < 1 ) { fprintf( ioQQQ, " j too small bigjump check\n" ); puts( "[Stop in PrtComment]" ); ShowMe(); exit(1); } if( BigJump > 0.4 ) { sprintf( chLine, " C-A temperature discontinuity occurred at zone%4ld from%10.2eK to%10.2eK.", j, struc.testr[j-1], struc.testr[j] ); caunin(chLine); } else if( BigJump > 0.2 ) { sprintf( chLine, " !A temperature discontinuity occurred at zone%4ld from%10.2eK to%10.2eK.", j, struc.testr[j-1], struc.testr[j] ); bangin(chLine); } } /* check for largest error in local electron density */ if( fabs(conv.BigEdenError) > ElecDen.EdenError ) { /* this only produces a warning if not the very last zone */ if( fabs(conv.BigEdenError) > ElecDen.EdenError*10. && ElecDen.nzEdenBad != ZoneCnt.nzone ) { sprintf( chLine, " W-The local error in the electron density reached%5.1f%% at zone%5ld", conv.BigEdenError*100, ElecDen.nzEdenBad ); warnin(chLine); } else if( fabs(conv.BigEdenError) > ElecDen.EdenError*5. ) { sprintf( chLine, " C-The local error in the electron density reached%5.1f%% at zone%5ld", conv.BigEdenError*100, ElecDen.nzEdenBad ); caunin(chLine); } else { sprintf( chLine, " The local error in the electron density reached%5.1f%% at zone%5ld", conv.BigEdenError*100, ElecDen.nzEdenBad ); notein(chLine); } } /* check for temperature oscillations */ BigJump = 0.; j = 0; for( i=1; i < (ZoneCnt.nzone - 1); i++ ) { big = fabs( (struc.testr[i-1] - struc.testr[i])/struc.testr[i] ); bigm = fabs( (struc.testr[i] - struc.testr[i+1])/struc.testr[i] ); /* this is sign of change in temperatre, we are looking for change in sign */ rel = ( (struc.testr[i-1] - struc.testr[i])/struc.testr[i])* ( (struc.testr[i] - struc.testr[i+1])/struc.testr[i] ); if( rel < 0. && MIN2( bigm , big ) > BigJump ) { j = i; BigJump = MIN2( bigm , big ); } } if( BigJump > 0.1 ) { /* only do sanity check if jump detected */ if( j < 1 ) { fprintf( ioQQQ, " j too small bigjump2 check\n" ); puts( "[Stop in PrtComment]" ); ShowMe(); exit(1); } if( BigJump > 0.1 ) { sprintf( chLine, " C-A temperature oscillation occurred at zone%4ld by%5.0f%% from%10.2e to %10.2e to%10.2e", j, BigJump*100., struc.testr[j-1], struc.testr[j], struc.testr[j+1] ); caunin(chLine); } else if( BigJump > 0.3 ) { sprintf( chLine, " !A temperature oscillation occurred at zone%4ld by%5.0f%% from%10.2e to %10.2e to%10.2e", j, BigJump*100., struc.testr[j-1], struc.testr[j], struc.testr[j+1] ); bangin(chLine); } } /* check for eden oscillations */ if( strcmp(pressure.chCPres,"CDEN") == 0 ) { j = 0; BigJump = 0.; for( i=1; i < (ZoneCnt.nzone - 1); i++ ) { big = (struc.ednstr[i-1] - struc.ednstr[i])/struc.ednstr[i]; if( fabs(big) < ElecDen.EdenError ) big = 0.; bigm = (struc.ednstr[i] - struc.ednstr[i+1])/struc.ednstr[i]; if( fabs(bigm) < ElecDen.EdenError ) bigm = 0.; if( big*bigm < 0. && fabs(struc.ednstr[i-1]-struc.ednstr[i])/struc.ednstr[i] > BigJump ) { j = i; BigJump = fabs(struc.ednstr[i-1]-struc.ednstr[i])/ struc.ednstr[i]; } } if( j < 1 ) /* only check on j if there was a big jump detected, number must be * smallest jump */ if( BigJump > ElecDen.EdenError*3. ) { if( j < 1 ) { fprintf( ioQQQ, " j too small bigjump3 check\n" ); puts( "[Stop in PrtComment]" ); ShowMe(); exit(1); } if( BigJump > ElecDen.EdenError*10. ) { sprintf( chLine, " C-An electron density oscillation occurred at zone%4ld by%5.0f%% from%10.2e to %10.2e to%10.2e", j, BigJump*100., struc.ednstr[j-1], struc.ednstr[j], struc.ednstr[j+1] ); caunin(chLine); } else if( BigJump > ElecDen.EdenError*3. ) { sprintf( chLine, " !An electron density oscillation occurred at zone%4ld by%5.0f%% from%10.2e to %10.2e to%10.2e", j, BigJump*100., struc.ednstr[j-1], struc.ednstr[j], struc.ednstr[j+1] ); bangin(chLine); } } } if( the1loCom.nhe1lo != 0 ) { if( phycon.hden > 1e7 ) { sprintf( chLine, " C-Temp/ion too low and matrix not used for He singlets in some zones. Low density Case B assumed.%5ld THE1LO=%10.3e", the1loCom.nhe1lo, the1loCom.the1lo ); caunin(chLine); } else { sprintf( chLine, " Temp/ion too low and matrix not used for He singlets in some zones. Low density Case B assumed.%5ld THE1LO=%10.3e", the1loCom.nhe1lo, the1loCom.the1lo ); notein(chLine); } } /********************************************************** * check that the volumn integrates out ok * **********************************************************/ /* this was the number 1 fed through the line integrators, * the number 1e-10 is sent to linadd in lineset1 as follows:*/ /*linadd( 1.e-10 , 1 , "Unit" , 'i' );*/ i = cdLine( "Unit" , 1 , &rate , &absint ); assert( i> 0 ); /* this is now the linear vol, rel to inner radius */ VolComputed = LineSv[i].sumlin / 1e-10 ; /* set stored number to zero so it does not appear on the emission line list */ LineSv[i].sumlin = 0.; /* spherical or plane parallel case? */ if( radius.Radius/radius.rinner > 1.0001 ) { /* spherical case, * sphere.iEmissPower is usually 3, * and can be reset to 2 or 1 with slit option on sphere */ /*VolExpected = sphere.covgeo*filfac.FillFac*radius.rinner/3.* ( POW3( radius.Radius/radius.rinner) - 1. );*/ VolExpected = sphere.covgeo*filfac.FillFac*radius.rinner/sphere.iEmissPower* ( powi( radius.Radius/radius.rinner,sphere.iEmissPower ) - 1. ); } else { /* plane parallel case */ /* next can be zero for very thin model, depth is always positive */ /*VolExpected = radius.Radius - radius.rinner;*/ VolExpected = sphere.covgeo*filfac.FillFac*radius.depth; } /* now get the relative difference between computed and expected volumns */ error = fabs(VolComputed - VolExpected)/VolExpected; /* we need to ignore this test if filling factor changes with radius, or * cylinder geometry in place */ if( radius.lgCylnOn || filfac.filpow!=0. ) { error = 0.; } /* how large is relative error? */ if( error > 0.001 ) { sprintf( chLine, " W-PrtComment insanity - Line unit integration did not verify \n"); warnin(chLine); fprintf( ioQQQ, " PrtComment insanity - Line unit integration did not verify \n"); fprintf( ioQQQ, " expected, derived vols were %g %g \n", VolExpected , VolComputed ); fprintf( ioQQQ, " relative difference is %g, ratio is %g.\n",error,VolComputed/VolExpected); insane(); ShowMe(); } /* next do same thing for fake continuum point propagated in highest energy cell, plus 1 * = * the variable rfield.ThrowOut[rfield.nflux] and rfield.DiffLocal[rfield.nflux] * are set to * the number 1.e-10f * Dilution in RTDiffuse. this is the outward * local emissivity, per unit vol. It is then added to the outward beams * by the rest of the code, and then checked here. * * insanity will be detected if diffuse emission is thrown into the outward beam * in MadeDiffuse. this happens if the range of ionization encompasses the full * continuum array, up to nflux. */ ConComputed = (rfield.outcon[rfield.nflux]+rfield.ConOutNoInter[rfield.nflux]) / 1e-10; /* correct for fraction that went out, as set in StartZone, * this should now be the volumn of the emitting region */ ConComputed /= ( (1. + sphere.covrt)/2. ) ; /* we expect this to add up to the integral of unity over r^-2 */ if( radius.Radius/radius.rinner < 1.001 ) { /* plane parallel case, no dilution, use thickness */ ConExpected = radius.depth*filfac.FillFac; } else { /* spherical case */ ConExpected = radius.rinner*filfac.FillFac * (1. - radius.rinner/radius.Radius ); } /* this is impossible */ assert( ConExpected > 0. ); /* now get the relative difference between computed and expected volumns */ error = fabs(ConComputed - ConExpected)/ConExpected; /* we need to ignore this test if filling factor changes with radius, or * cylinder geometry in place */ if( radius.lgCylnOn || filfac.filpow!=0. ) { error = 0.; } /* how large is relative error? */ if( error > 0.001 ) { sprintf( chLine, " W-PrtComment insanity - Continuum unit integration did not verify \n"); warnin(chLine); fprintf( ioQQQ," PrtComment insanity - Continuum unit integration did not verify \n"); fprintf( ioQQQ," exact vol= %g, derived vol= %g relative difference is %g \n", ConExpected , ConComputed ,error); fprintf( ioQQQ," outcon= %g, ConOutNoInter= %g \n", rfield.outcon[rfield.nflux],rfield.ConOutNoInter[rfield.nflux] ); insane(); ShowMe(); } /* final prinout of warnings, cautions, notes, */ cdNwcns(&nw,&nc,&nn,&ns,&i,&j,&ipHi,&ipLo); if( nw > 0 ) { warngs.lgWarngs = TRUE; } else { warngs.lgWarngs = FALSE; } if( nc > 0 ) { warngs.lgCautns = TRUE; } else { warngs.lgCautns = FALSE; } if( called.lgTalk ) { /* print the title of the calculation */ fprintf( ioQQQ, " %s\n", input.chTitle ); /* say why the calculation stopped, and indicate the geometry*/ cdReasonGeo(ioQQQ); /* print all warnings */ cdWarnings(ioQQQ); /* all cautions */ cdCautions(ioQQQ); /* surprises, begining with a ! */ cdSurprises(ioQQQ); /* notes about the calculations */ cdNotes(ioQQQ); } /* option to print warnings on special io */ if( lgPrnErr ) { cdWarnings(ioPrnErr); } if( trace.lgTrace ) { fprintf( ioQQQ, " PrtComment returns.\n" ); } # ifdef DEBUG_FUN fputs( " <->PrtComment()\n", debug_fp ); # endif return; }