/*NextDR use adaptive logic to find next zone thickness */ #include #include #include "cddefines.h" #include "cdminmax.h"/* vfmax, vfmin */ #define DNGLOB 0.10 #define Z 1.0001 #include "nhe3lvl.h" #include "angllum.h" #include "showme.h" #include "opac.h" #include "globul.h" #include "drmin.h" #include "heating.h" #include "sdrmin.h" #include "converge.h" #include "hmi.h" #include "fluct.h" #include "zonecnt.h" #include "stopcalc.h" #include "colden.h" #include "h2opac.h" #include "ionfracs.h" #include "heat.h" #include "h.h" #include "phycon.h" #include "radius.h" #include "drminu.h" #include "trace.h" #include "filfac.h" #include "didz.h" #include "wind.h" #include "ipdr.h" #include "drneg.h" #include "taulines.h" #include "pressure.h" #include "densty.h" #include "itercnt.h" #include "sumoutinc.h" #include "tseton.h" #include "dtmase.h" #include "linelabl.h" #include "struc.h" #include "toler.h" #include "he3gams.h" #include "contrate.h" #include "chkrate.h" #include "fndlineht.h" #include "grainratedr.h" #include "insane.h" #include "tabden.h" #include "fabden.h" #include "chlinelbl.h" #include "nextdr.h" void NextDR(void) { char chDestAtom[9], chLbl[11]; int lgDoPun, lgOscil; long int icarstag, iironstag, initstag, ioxystag, ipStrong, istage, k, level, limit; double ThicknessSet , dThickness , drThickness , DepthToGo ; double DrGrainHeat, GlobDr, Out2Tot, SaveOHIIoHI, SpecDr, Strong, SvHe2oHe1, TauDTau, TauInwd, coleff, dDRCarbon, dDRIron, dDRNitrogen, dDROxygen, dDestRate, dEden, dHdStep, dRTaue, dTdStep, dnew, drConPres, drDepth, drDest, drEden, drFail, drFluc, drH2Band, drHMase, drHe1ion, drHion, drInter, drLineHeat, drOutFrac, drPressure, drSphere, drTab, drdHdStep, drdTdStep, drmax, dt, error, fac2, factor, freqm, grfreqm, gropacm, hdnew, opacm, OldDR , OldestEden, winddr, x; static double OHIIoHI, OHe2oHe1 = 0., OldEden, OldHeat = 0., OldTe = 0.; static double BigRadius = 1e30; # ifdef DEBUG_FUN fputs( "<+>NextDR()\n", debug_fp ); # endif /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> * * changes in logic * 95 oct 19, drSphere now 3% of radius, was 2%, fewer zone * 95 oct 19, subtracted grain opacity from total opacity used * to get thickness in routine ContRate * *>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> * * free statement labels >= 13 * *----------------------------------------------------------------------- * * this sub determines the thickness of the next zone * if is called one time for each zone * flag lgNxtDROff is true if this is initialization of NextDR, * is false if we are to use logic to find dr * *----------------------------------------------------------------------- */ if( trace.lgTrace ) { fprintf( ioQQQ, " NextDR called\n" ); } /* save current dr */ OldDR = radius.drad; /* 1 '' NextDR keys from change in H ionization'',e11.3)') * check on change in hydrogen ionizaiton */ if( ZoneCnt.nzone <= 1 ) { if( h.hii > 0. && h.hi > 0. ) { OHIIoHI = h.hii/h.hi; } else { OHIIoHI = 0.; } drHion = BigRadius; SaveOHIIoHI = h.hii/h.hi; /* else if(hii.gt.0. .and. hi.gt.0. .and. OHIIoHI.gt.0. ) then * >>chng 97 july 9, for deep in PDR vastly now ionz H slowed down works */ } else if( (h.hii > 0. && h.hi > 0.) && OHIIoHI > 1e-15 ) { /* ratio of current HII/HI to old value - < 1 when becoming more neutral */ x = (h.hii/h.hi)/OHIIoHI; if( x < 1. ) { x = 1. - x; /* >>chng 96 feb 23 from 0.7 to 0.3 because punched through H i-front * >>chng 97 may 5, from 0.3 to 0.2 since punched through H i-front */ drHion = radius.drad*MAX2( 0.2 , 0.2/MAX2(1e-10,x) ); } else { drHion = BigRadius; } SaveOHIIoHI = OHIIoHI; OHIIoHI = h.hii/h.hi; } else { SaveOHIIoHI = OHIIoHI; if( h.hii > 0. && h.hi > 0. ) { OHIIoHI = h.hii/h.hi; } else { OHIIoHI = 0.; } drHion = BigRadius; } /* "NextDR keys from H maser, dt, ij=" possible hydrogen maser action */ if( dtMase.dTauMase < -0.01 ) { /* maser so powerful, must limit inc in tay * >>chng 96 aug8, from 0.3 to 0.1 due to large maser crash */ drHMase = radius.drad*MAX2(0.1,-0.2/dtMase.dTauMase); } else { drHMase = BigRadius; } /* NextDR keys from change in He2/He1 ionization, old=%11.3e sv=%11.3e FrLya*/ /* check on change in HeI ionization */ SvHe2oHe1 = OHe2oHe1; if( ZoneCnt.nzone <= 1 ) { if( xIonFracs[2][ipHE] > 0. && xIonFracs[1][ipHE] > 0. ) { OHe2oHe1 = xIonFracs[2][ipHE]/xIonFracs[1][ipHE]; } else { OHe2oHe1 = 0.; } drHe1ion = BigRadius; } else if( (xIonFracs[2][ipHE] > 0. && xIonFracs[1][ipHE] > 0.) && OHe2oHe1 > 0. ) { /* >>chng 97 july 14, from 1e-2 to 1e-3 */ /* >>chng 99 feb 16, this logic included a test for whether * the ots fields of Lya were an important destruction mechanism * for the triplets. changes in the helium ionization were NOT * counted when ionization fraction was small and lya important. * changed to simpler logic, only count changes when there is * a significant amount of He+ */ /*if( (he3gams.He3FrLya > 1e-3) && (he.heii/he.hei <= 1e-3) )*/ if( xIonFracs[2][ipHE]/xIonFracs[1][ipHE] < 1e-3 ) { /* >>chng 97 jul 12, He mostly neutral but Lya photoionized, * oscillation could * take place - do not trigger dr from this oscillation */ drHe1ion = BigRadius; } else { /* this is the ratio of the current to old He2/He1 ratio */ x = (xIonFracs[2][ipHE]/xIonFracs[1][ipHE])/OHe2oHe1; if( x < 1. ) { x = 1. - x; /* >>chng 97 may 5, ran into i-fronts, * >>chng 99 feb 24, from 0.02 to 0.03 in below */ drHe1ion = radius.drad* MAX2(0.3,0.3/MAX2(1e-10,x)); } else { drHe1ion = BigRadius; } } OHe2oHe1 = xIonFracs[2][ipHE]/xIonFracs[1][ipHE]; } else { if( xIonFracs[2][ipHE] > 0. && xIonFracs[1][ipHE] > 0. ) { OHe2oHe1 = xIonFracs[2][ipHE]/xIonFracs[1][ipHE]; } else { OHe2oHe1 = 0.; } drHe1ion = BigRadius; } /* check on how much outward flux is being added relative to * attenuated incident beam * 1 '' NextDR keys from ratio out2 tot continuua'',e9.1)') */ if( (SumOutInc.SumIncCon > 0. && SumOutInc.SumOutCon > 0.) && !tseton.lgTSetOn ) { /* evaluated in metdif */ Out2Tot = SumOutInc.SumOutCon/SumOutInc.SumIncCon; /* Out2HTot = SumHOutCon / SumHIncCon * try to keep outward frac less than 5% of total */ drOutFrac = radius.drad*MAX2(0.7,0.05/MAX2(1e-10,Out2Tot)); } else { Out2Tot = FLT_MAX; drOutFrac = BigRadius; } /* check how heating is changing * '' NextDR keys from change in heating; current, delta='', */ if( ZoneCnt.nzone <= 1 || tseton.lgTSetOn ) { drdHdStep = BigRadius; dHdStep = FLT_MAX; } else { dHdStep = fabs(heat.htot-OldHeat)/heat.htot; if( dHdStep > 0. ) { if( phycon.hden >= 1e13 ) { /* drdHdStep = drad * MAX( 0.8 , 0.05/dHdStep ) */ drdHdStep = radius.drad*MAX2(0.8,0.075/dHdStep); } else if( phycon.hden >= 1e11 ) { /* drdHdStep = drad * MAX( 0.8 , 0.075/dHdStep ) */ drdHdStep = radius.drad*MAX2(0.8,0.100/dHdStep); } else { /* changed from .15 to .12 for outer edge of coolhii too steep dT * changed to .10 for symmetry, big change in some rates, 95nov14 * changed from .10 to .125 since parispn seemed to waste zones * >>chng 96 may 21, from .125 to .15 since pn's still waste zones */ drdHdStep = radius.drad*MAX2(0.8,0.15/dHdStep); } } else { drdHdStep = BigRadius; } } OldHeat = heat.htot; /* pressure due to incident continuum if in eos */ if( strcmp(pressure.chCPres,"CPRE") == 0 && pressure.lgConPres ) { if( ZoneCnt.nzone > 2. && pressure.pinzon > 0. ) { /* pinzon is pressrue from acceleration onto previos zone * want this less than some fraction of total pressure */ drConPres = 0.05*pressure.PressureInit/(wind.AccelTot* densty.xMassDensity*filfac.FillFac); } else { drConPres = BigRadius; } } else { drConPres = BigRadius; } /* check how temperature is changing * 1 '' NextDR keys from change in temperature; current, delta='', */ if( ZoneCnt.nzone <= 1 ) { drdTdStep = BigRadius; dTdStep = FLT_MAX; } else { dTdStep = fabs(phycon.te-OldTe)/phycon.te; if( phycon.te > 1e4 ) { /* >>chng 96 may 30, for high temps with flat cooling fcn, * do not want too fine a test */ x = 0.05; } else { /* >>chng 95 oct 13 from 0.05 to 0.03 for big te changes at * outer edges of pn's and hii regions */ x = 0.03; } if( dTdStep > 0. ) { /* >>chng 96 may 30, variable depending on temp * >>chng 96 may 31, allow 0.7 smaller, was 0.8 * >>chng 97 may 5, from 0.7 to 0.5 stop from punching through thermal front */ drdTdStep = radius.drad*MAX2(0.7,x/dTdStep); } else { drdTdStep = BigRadius; } } OldTe = phycon.te; /* find freq where opacity largest and interaction rate is fastest * 1 '' NextDR keys from cont inter nu='',e9.2, */ ContRate(&freqm,&opacm); /* drChange is optical depth at max interaction energy * >>chng 96 june 6 was drChange=0.15 changed to 0.3 for high Z models * taking too many zones * drInter = drChange / MAX(1e-30,opacm*FillFac) */ drInter = 0.3/MAX2(1e-30,opacm*filfac.FillFac*Angllum.AngleIllum); /* check on changes in destruction rates for various atoms * call ChkRate( 1 , dDRHydro ,istage ) * call ChkRate( 2 , dDRHelium ,istage ) */ ChkRate(6,&dDRCarbon,&icarstag); ChkRate(7,&dDRNitrogen,&initstag); ChkRate(8,&dDROxygen,&ioxystag); ChkRate(26,&dDRIron,&iironstag); dDestRate = vfmax(dDROxygen,dDRIron,dDRCarbon,dDRNitrogen,FEND); if( dDRCarbon == dDestRate ) { dDestRate = dDRCarbon; istage = icarstag; strcpy( chDestAtom, "Carbon " ); } else if( dDRNitrogen == dDestRate ) { dDestRate = dDRNitrogen; istage = initstag; strcpy( chDestAtom, "Nitrogen" ); } else if( dDROxygen == dDestRate ) { dDestRate = dDROxygen; strcpy( chDestAtom, "Oxygen " ); istage = ioxystag; } else if( dDRIron == dDestRate ) { dDestRate = dDRIron; istage = iironstag; strcpy( chDestAtom, "Iron " ); } else { fprintf( ioQQQ, " insanity following ChkRate\n" ); ShowMe(); puts( "[Stop in nextdr]" ); exit(1); } /* 1 '' NextDR keys from change in dest rates, atom='',a8,i3)') */ if( dDestRate > 0. ) { /* if( te.gt.40 000. ) then * added different accuracy for hot gas since tend to jump over * big te range for small change in heat (intrinsically unstable) * drDest = drad * MAX( 0.5 , 0.10/dDestRate ) * else * was drChange, changed to .15 for parishii go through HeII=HeI I front * drDest = drad * MAX( 0.5 , 0.15/dDestRate ) * >>chng 95 dec 18 from above to below to stop oscillations * >>chng 95 dec 27 from min of .5 to .75 to stop zone size from changing rapidly * drDest = drad * MAX( 0.8 , 0.20 /dDestRate ) * >>chng 96 jan 14 from .2 to .25 to use less zones * >>chng 96 may 30 from min of 0.8 to 0.5 to prevent crashing into He i-front */ drDest = radius.drad*MAX2(0.5,0.20/dDestRate); /* endif */ } else { drDest = BigRadius; } /* check whether change in wind velocity constrains DRAD, * this is usual hypersonic wind */ if( wind.windv > 0. ) { if( wind.dVeldRad != 0. ) { /* dVeldRad is D(vel)/vel / DRAD, computed in convpres */ winddr = 0.03/wind.dVeldRad; } else { winddr = 1.1*radius.drad; } } /* this is D-critical wind, dVdr is negative, and we want * small changes */ else if( wind.windv < 0. ) { if( wind.dVeldRad != 0. ) { /* this is change in velocity per zone, in cm/sec */ const double dVel=3e4; /* dVeldRad is D(vel)/vel / DRAD, computed in convpres */ winddr = dVel / fabs(wind.dVeldRad); } else { winddr = 1.1*radius.drad; } } else { winddr = BigRadius; } /* inside out globule */ if( strcmp(pressure.chCPres,"GLOB") == 0 ) { if( globul.glbdst < 0. ) { fprintf( ioQQQ, " Globule distance is negative, internal overflow has occured, sorry.\n" ); fprintf( ioQQQ, " This is routine NextDR, GLBDST is%10.2e\n", globul.glbdst ); puts( "[Stop in nextdr]" ); exit(1); } factor = globul.glbden*pow(globul.glbrad/globul.glbdst,globul.glbpow); fac2 = globul.glbden*pow(globul.glbrad/(globul.glbdst - radius.drad),globul.glbpow); if( fac2/factor > 1. + DNGLOB ) { /* DNGLOB is relative change in density allowed this zone, 0.10 */ GlobDr = radius.drad*DNGLOB/(fac2/factor - 1.); } else { GlobDr = BigRadius; } /* GlobDr = GLBDST * DNGLOB * (GLBRAD/GLBDST)**(-GLBPOW) / GLBPOW */ GlobDr = MIN2(GlobDr,globul.glbdst/20.); } else { GlobDr = BigRadius; } hdnew = 0.; if( strncmp( pressure.chCPres , "DLW" , 3) == 0 ) { /* one of the special density laws, first get density at possible next dr */ if( strcmp(pressure.chCPres,"DLW1") == 0 ) { hdnew = fabden(radius.Radius+radius.drad,radius.depth+ radius.drad); } else if( strcmp(pressure.chCPres,"DLW2") == 0 ) { hdnew = tabden(radius.Radius+radius.drad,radius.depth+ radius.drad); } else { fprintf( ioQQQ, " dlw insanity in NextDR\n" ); puts( "[Stop in nextdr]" ); exit(1); } drTab = fabs(hdnew-phycon.hden)/MAX2(hdnew,phycon.hden); drTab = radius.drad*MAX2(0.2,0.10/MAX2(0.01,drTab)); } else { drTab = BigRadius; } /* special density law */ if( strcmp(pressure.chCPres,"DLW1") == 0 ) { dnew = fabs(fabden(radius.Radius+radius.drad,radius.depth+ radius.drad)/phycon.hden-1.); /* DNGLOB is relative change in density allowed this zone, 0.10 */ if( dnew == 0. ) { SpecDr = radius.drad*3.; } else if( dnew/DNGLOB > 1.0 ) { SpecDr = radius.drad/(dnew/DNGLOB); } else { SpecDr = BigRadius; } } else { SpecDr = BigRadius; } /* check grain line heating dominates * this is important in PDR and HII region calculations * >>chng 97 july 3, added following check */ if( HeatingCom.heating[13][0]/heat.htot > 0.2 ) { GrainRateDr(&grfreqm,&gropacm); DrGrainHeat = 1.0/MAX2(1e-30,gropacm*filfac.FillFac*Angllum.AngleIllum); } else { gropacm = 0.; grfreqm = 0.; DrGrainHeat = BigRadius; } /* check if line heating dominates * this is important in high metallicity models */ if( HeatingCom.heating[22][0]/heat.htot > 0.2 ) { FndLineHt(&level,&ipStrong,&Strong); if( Strong/heat.htot > 0.1 ) { if( level == 1 ) { /* a level1 line was the heat source (usual case) * drLineHeat = MAX(1.0,TauLines(ipLnTauIn,ipStrong)*0.2) / * 1 TauLines(ipLnDTau,ipStrong) */ TauInwd = TauLines[ipStrong-1].TauIn; TauDTau = TauLines[ipStrong-1].dTau; } else if( level == 2 ) { /* a level2 line was the heat source * (bad case since not as well treated) * drLineHeat = MAX(1.0,WindLine(ipLnTauIn,ipStrong)*0.2) / * 1 WindLine(ipLnDTau,ipStrong) */ TauInwd = TauLine2[ipStrong-1].TauIn; TauDTau = TauLine2[ipStrong-1].dTau; } else { /* this is insane, since Strong was set, but not level */ insane(); fprintf( ioQQQ, " NextDR Strong line heating set, but not level.\n" ); drLineHeat = BigRadius; ipStrong = 1; Strong = 0.; puts( "[Stop in nextdr]" ); exit(1); } /* in following logic cannot use usual inverse opacity, * since line heating competes with escape probability, * so is effective at surprising optical depths */ if( TauDTau > 0. ) { drLineHeat = MAX2(1.,TauInwd)*0.4/TauDTau; } else { drLineHeat = BigRadius; } } else { TauInwd = 0.; drLineHeat = BigRadius; ipStrong = 1; Strong = 0.; } } else { TauInwd = 0.; drLineHeat = BigRadius; ipStrong = 1; level = 0; Strong = 0.; } /* Lyman band optical depth * second term total lyman band optical depth, do not check if thick * H2Opacity is opacity evaluated in hmole */ if( hmi.htwo*h2opac.H2Opacity > 1e-20 && coldenCom.colden[IPCOLH2-1]* h2opac.H2Opacity < 5. ) { if( coldenCom.colden[IPCOLH2-1]*h2opac.H2Opacity < 2. ) { /* special extra 2 since H2 will change by this amount early on */ drH2Band = (didz.drChange/2.)/(hmi.htwo*h2opac.H2Opacity* filfac.FillFac); } else { drH2Band = didz.drChange/(hmi.htwo*h2opac.H2Opacity*filfac.FillFac); } } else { drH2Band = BigRadius; } /* can't make drmax large deep in model since causes feedback * oscillations with changes in heating or destruction rates * >>chng 96 Oct 15, change from 5 to 10 */ if( ZoneCnt.nzone < 11 ) { lgOscil = FALSE; if( ZoneCnt.nzone < 5 ) { /* >>chng 96 nov 11, had been 4 * drad up to 11, change to following * to be similar to c90.01, max of 1.3 between 5 and 11 */ drmax = 4.*radius.drad/Angllum.AngleIllum; } else { drmax = 1.3*radius.drad/Angllum.AngleIllum; } } else { /* >>chng 96 Oct 15, do not let zones increase if oscillations present */ lgOscil = FALSE; /* this routine called before tauinc stores variables * >>chng 96 oct 31, error to declare oscillation propto toler, the *heating cooling tolerance */ error = tolerCom.toler*tolerCom.toler; limit = MIN2(NZLIM-2,ZoneCnt.nzone-2); for( k=ZoneCnt.nzone - 10; k < limit; k++ ) { /* small residiual is to allow 0.01 rel error */ if( (struc.testr[k-1] - struc.testr[k])/struc.testr[k]* (struc.testr[k] - struc.testr[k+1])/struc.testr[k] < -error ) { lgOscil = TRUE; } /* small residiual is to allow 0.01 rel error */ if( (struc.ednstr[k-1] - struc.ednstr[k])/struc.ednstr[k]* (struc.ednstr[k] - struc.ednstr[k+1])/struc.ednstr[k] < -error ) { lgOscil = TRUE; } } if( lgOscil ) { /* >>chng 96 Oct 15, do not let zones increase if oscillations present */ drmax = radius.drad; } else { /* >>chng 96 jan 15 was 1.15 - why so small? */ drmax = 1.3*radius.drad; } } /* is current pressure ok? */ if( conv.lgConvPres ) { drPressure = BigRadius; } else if( ZoneCnt.nzone > 0 ) { /* first zone usually has pressure fail due to bug in ptot */ drPressure = radius.drad/2.; } else { drPressure = BigRadius; } if( !conv.lgConvTemp ) { drFail = radius.drad/2.; } else { drFail = BigRadius; } /* change in electron density - NextDR keys from change in elec den, * remember old electron density during first call */ if( ZoneCnt.nzone == 0 ) { OldEden = phycon.eden; } /* this is low ionization solution */ if( phycon.eden/phycon.hden <= 1. ) { dEden = fabs(phycon.eden-OldEden)/phycon.eden; if( dEden > 0. ) { /* >>chng 96 may 15, down to .03 from 0.04 * >>chng 96 may 17, up to 0.1 from .03 */ /* >>chng 99 Nov 23, smallest change is 0.3, had been 0.7, * smaller value needed due to large initial changes in ionization * in 94 since code does not make attempt to identify drnext * and tay in middle of zone, as result some near-lte dense models * can have rapid changes in ionization when due to pumped lines, * see dalton.in for examples */ /* >>chng 99 nov 25, allow dr to become as much as 0.2 smaller, * for dense moels where initial change in eden is large at ill face. * had been 0.2 */ drEden = radius.drad*MAX2(0.2,0.10/dEden); } else { drEden = BigRadius; } } else { dEden = 0.; drEden = BigRadius; } OldestEden = OldEden; OldEden = phycon.eden; /* do not let thickness get too large * 1 '' NextDR keys from relative depth'',e11.3)') */ if( ZoneCnt.nzone > 20 ) { drDepth = radius.depth/5.; } else { drDepth = BigRadius; } /* case where stopping thickness or edge specified, need to approach slowly */ ThicknessSet = -1.; dThickness = -1.; DepthToGo = -1.; if( StopCalc.HColStop < 5e29 ) { ThicknessSet = StopCalc.HColStop; dThickness = phycon.hden*filfac.FillFac; DepthToGo = StopCalc.HColStop-coldenCom.colden[IPCOLUMNDENSITY-1]; } else if( StopCalc.colpls < 5e29 ) { ThicknessSet = StopCalc.colpls; dThickness = (double)(h.hii)*filfac.FillFac; DepthToGo = StopCalc.colpls-coldenCom.colden[IPCHII-1]; } else if( StopCalc.colnut < 5e29 ) { ThicknessSet = StopCalc.colnut; dThickness = (double)(h.hi)*filfac.FillFac; DepthToGo = StopCalc.colnut - coldenCom.colden[IPCHI-1] ; } /* this is case where outer radius is set */ if( radius.router[IterCnt.iter-1] < 5e29 ) { ThicknessSet = radius.router[IterCnt.iter-1] ; dThickness = 1.; DepthToGo = radius.router[IterCnt.iter-1] - radius.depth ; } /* set dr if one of above tests have triggered */ if( ThicknessSet > 0. && DepthToGo > 0.) { /* want to approach the outer edge slowly, * the need for this logic is most evident in brems.in - * HI fraction varies across coronal model */ DepthToGo = MIN2( ThicknessSet/10. , DepthToGo ); drThickness = DepthToGo / dThickness ; } else { drThickness = BigRadius ; } /*fprintf(ioQQQ," drThickness = %e %e %e \n", drThickness, DepthToGo , dThickness );*/ /* spherical models, do not want delta R/R big */ drSphere = radius.Radius*0.03; drSphere = radius.Radius*0.04; /* optical depth to electron scattering */ dRTaue = didz.drChange/(phycon.eden*6.65e-25); if( fluct.flong != 0. ) { /* FacAbunSav = (cfirst * COS( depth*flong ) + csecnd) * this is option for rapid density fluctuations * set DR to roughly 0.1 of fluctuation length * read flcsub first parameter (log of period) is used to get flong * flong = 6.2831853 / period in flcsub */ drFluc = 0.628/2./fluct.flong; } else { drFluc = BigRadius; } /* if density fluctuations in place then override change in heat * for dr set */ if( strcmp(pressure.chCPres,"SINE") == 0 && fluct.lgDenFluc ) { drdHdStep = BigRadius; } /*active dr sets */ /* we are deep into model, use logic since we have several zones * of old data */ radius.drNext = vfmin(drmax,drInter,drLineHeat,winddr,drFluc,GlobDr, DrGrainHeat,FEND); radius.drNext = vfmin(radius.drNext,SpecDr,drFail,drSphere,sdrminCom.sdrmax, drPressure,dRTaue,FEND); radius.drNext = vfmin(radius.drNext,drH2Band,drDest,drdTdStep,drdHdStep, drConPres,drTab,FEND); radius.drNext = vfmin(radius.drNext,drOutFrac,drHion,drHe1ion,drDepth, drEden,drHMase,drThickness,FEND); /* keep dr constant in first two zones, if it wants to increase, * but always allow it to decrease. * to guard against large increases in efrac for balmer cont photo dominated models, */ if( ZoneCnt.nzone <= 1 && radius.drNext > OldDR) { radius.drNext = OldDR; } /* option to force min drad */ if( radius.drNext < sdrminCom.sdrmin ) { radius.drNext = sdrminCom.sdrmin; } /* dr = zero is a logical mistake */ if( radius.drNext <= 0. ) { fprintf( ioQQQ, " NextDR finds insane drNext:%10.2e\n", radius.drNext ); fprintf( ioQQQ, " all drs follow:%10.2e%10.2e%10.2e%10.2e%10.2e\n all drs follow:%10.2e%10.2e%10.2e%10.2e%10.2e\n all drs follow:%10.2e%10.2e%10.2e%10.2e%10.2e\n", drmax, drInter, drLineHeat, winddr, drFluc, GlobDr, SpecDr, drFail, drSphere, sdrminCom.sdrmax, drPressure, dRTaue, drH2Band, drOutFrac, drDepth ); puts( "[Stop in nextdr]" ); exit(1); } /* all this is to only punch on last iteration * the punch dr command is not really a punch command, making this necessary * lgDRon is set true if "punch dr" entered */ if( ipdr.lgDROn ) { if( ipdr.lgDRPLst ) { /* lgDRPLst was set true if "punch" had "last" on it */ if( IterCnt.lgLastIt ) { lgDoPun = TRUE; } else { lgDoPun = FALSE; } } else { lgDoPun = TRUE; } } else { lgDoPun = FALSE; } if( (trace.lgTrace && trace.lgDrBug) || lgDoPun ) { if( !conv.lgConvTemp && ZoneCnt.nzone > 0 ) { fprintf( ipdr.ipDRout, " >>>> A temperature failure occured.\n" ); } if( !conv.lgConvPres && ZoneCnt.nzone > 0 ) { fprintf( ipdr.ipDRout, " >>>> A pressure failure occured.\n" ); } /*=======begin active dr sets */ if( radius.drNext == drLineHeat ) { if( level == 1 ) { strcpy( chLbl, chLineLbl(&TauLines[ipStrong-1]) ); fprintf( ipdr.ipDRout, " %4ld%11.3e NextDR keys from level 1 line heating,%10.10s TauIn%10.2e%10.2e%10.2e\n", ZoneCnt.nzone, radius.drNext, chLbl, TauInwd, TauLines[ipStrong-1].pump, TauLines[ipStrong-1].Pesc ); } else if( level == 2 ) { strcpy( chLbl, chLineLbl(&TauLine2[ipStrong-1])); fprintf( ipdr.ipDRout, " %4ld%11.3e NextDR keys from level 2 line heating,%10.10s TauIn%10.2e%10.2e%10.2e\n", ZoneCnt.nzone, radius.drNext, chLbl, TauInwd, TauLine2[ipStrong-1].pump, TauLine2[ipStrong-1].Pesc ); } else { fprintf( ioQQQ, " insanity pr line heat\n" ); puts( "[Stop in nextdr]" ); exit(1); } } else if( radius.drNext == drDepth ) { fprintf( ipdr.ipDRout, " %4ld%11.3e NextDR keys from relative depth\n", ZoneCnt.nzone, radius.drNext ); } else if( radius.drNext == drThickness ) { fprintf( ipdr.ipDRout, " %4ld%11.3e NextDR keys from depth to go\n", ZoneCnt.nzone, radius.drNext ); } else if( radius.drNext == drTab ) { fprintf( ipdr.ipDRout, " %4ld%11.3e NextDR keys from spec den law, new old den%10.2e%10.2e\n", ZoneCnt.nzone, radius.drNext, hdnew, phycon.hden ); } else if( radius.drNext == drHMase ) { fprintf( ipdr.ipDRout, " %4ld%11.3e NextDR keys from H maser, dt, ij=%10.2e%5ld\n", ZoneCnt.nzone, radius.drNext, dtMase.dTauMase, dtMase.ijMase ); } else if( radius.drNext == drHion ) { fprintf( ipdr.ipDRout, " %4ld%11.3e NextDR keys from change in H ionization fm to%11.3e%11.3e\n", ZoneCnt.nzone, radius.drNext, SaveOHIIoHI, OHIIoHI ); } else if( radius.drNext == drEden ) { fprintf( ipdr.ipDRout, " %4ld%11.3e NextDR keys from change in elec den, rel chng:%11.3e, cur %g old=%g\n", ZoneCnt.nzone, radius.drNext, dEden , OldEden , OldestEden ); } else if( radius.drNext == drHe1ion ) { fprintf( ipdr.ipDRout, " %4ld%11.3e NextDR keys from change in He2/He1 ionization, new=%11.3e old=%11.3e FrLya%11.3e\n", ZoneCnt.nzone, radius.drNext, OHe2oHe1, SvHe2oHe1, he3gams.He3FrLya ); } else if( radius.drNext == drDest ) { fprintf( ipdr.ipDRout, " %4ld%11.3e NextDR keys from change in dest rates, atom=%8.8s%3ld\n", ZoneCnt.nzone, radius.drNext, chDestAtom, istage ); } else if( radius.drNext == drOutFrac ) { fprintf( ipdr.ipDRout, " %4ld%11.3e NextDR keys from ratio out2 tot continuua%10.2e%10.2e%7ld %4.4s %4.4s\n", ZoneCnt.nzone, radius.drNext, Out2Tot, SumOutInc.SumOutMax, SumOutInc.ipSumOutMax, LineLabl.chLineLabel[SumOutInc.ipSumOutMax-1], LineLabl.chContLabel[SumOutInc.ipSumOutMax-1] ); /* in above SumOutMax is anu where peak occured */ } else if( radius.drNext == drdHdStep ) { fprintf( ipdr.ipDRout, " %4ld%11.3e NextDR keys from change in heating; current %10.3e delta=%10.3e\n", ZoneCnt.nzone, radius.drNext, heat.htot, dHdStep ); } else if( radius.drNext == drConPres ) { fprintf( ipdr.ipDRout, " %4ld%11.3e NextDR keys from change in con accel\n", ZoneCnt.nzone, radius.drNext ); } else if( radius.drNext == drdTdStep ) { fprintf( ipdr.ipDRout, " %4ld%11.3e NextDR keys from change in temperature; current= %10.3e, delta= %10.3e\n", ZoneCnt.nzone, radius.drNext, phycon.te, dTdStep ); } else if( radius.drNext == sdrminCom.sdrmin ) { fprintf( ipdr.ipDRout, " %4ld%11.3e NextDR keys from sdrmin\n", ZoneCnt.nzone, radius.drNext ); } else if( radius.drNext == drH2Band ) { fprintf( ipdr.ipDRout, " %4ld%11.3e NextDR keys from H2 Lyman Band - opacity=%10.2e\n", ZoneCnt.nzone, radius.drNext, hmi.htwo*h2opac.H2Opacity ); } else if( radius.drNext == sdrminCom.sdrmax ) { fprintf( ipdr.ipDRout, " %4ld%11.3e NextDR keys from sdrmax\n", ZoneCnt.nzone, radius.drNext ); } else if( radius.drNext == drSphere ) { fprintf( ipdr.ipDRout, " %4ld%11.3e NextDR keys from sphericity\n", ZoneCnt.nzone, radius.drNext ); } else if( radius.drNext == dRTaue ) { fprintf( ipdr.ipDRout, " %4ld%11.3e NextDR keys from optical depth to electron scattering\n", ZoneCnt.nzone, radius.drNext ); } else if( radius.drNext == drPressure ) { fprintf( ipdr.ipDRout, " %4ld%11.3e NextDR key from pressure failure\n", ZoneCnt.nzone, radius.drNext ); } else if( radius.drNext == drFail ) { fprintf( ipdr.ipDRout, " %4ld%11.3e NextDR key from temperature failure.\n", ZoneCnt.nzone, radius.drNext ); } else if( radius.drNext == drmax ) { fprintf( ipdr.ipDRout, " %4ld%11.3e NextDR key from DRMAX; nu opc dr=%10.2e%10.2e%10.2e oscil=%2c\n", ZoneCnt.nzone, radius.drNext, freqm, opacm, didz.drChange/ opacm, TorF(lgOscil) ); } else if( radius.drNext == drInter ) { fprintf( ipdr.ipDRout, " %4ld%11.3e NextDR keys from cont inter nu=%10.2e opac=%10.2e\n", ZoneCnt.nzone, radius.drNext, freqm, opacm ); } else if( radius.drNext == DrGrainHeat ) { fprintf( ipdr.ipDRout, " %4ld%11.3e NextDR keys from grain heating nu=%10.2e opac=%10.2e\n", ZoneCnt.nzone, radius.drNext, grfreqm, gropacm ); } else if( radius.drNext == winddr ) { fprintf( ipdr.ipDRout, " %4ld%11.3e NextDR key from Wind, dVeldRad=%10.3e\n", ZoneCnt.nzone, winddr, wind.dVeldRad ); } else if( radius.drNext == drFluc ) { fprintf( ipdr.ipDRout, " %4ld%11.3e NextDR key from density fluctuations\n", ZoneCnt.nzone, radius.drNext ); } else if( radius.drNext == GlobDr ) { fprintf( ipdr.ipDRout, " %4ld NextDR keys from GLOB law new dr=%10.2e HDEN=%10.2e\n", ZoneCnt.nzone, radius.drNext, phycon.hden ); } else if( radius.drNext == SpecDr ) { fprintf( ipdr.ipDRout, " %4ld NextDR keys from special law new dr=%10.2e HDEN=%10.2e\n", ZoneCnt.nzone, radius.drNext, phycon.hden ); } else if( radius.drNext == OldDR ) { fprintf( ipdr.ipDRout, " %4ld%11.3e NextDR uses old DR.\n", ZoneCnt.nzone, radius.drNext ); } else { fprintf( ipdr.ipDRout, " %4ld NextDR keys from insanity %10.2e\n", ZoneCnt.nzone, radius.drNext ); fprintf( ioQQQ, " %4ld NextDR keys from insanity %10.2e\n", ZoneCnt.nzone, radius.drNext ); puts( "[Stop in nextdr]" ); exit(1); } /*======end active dr sets */ } /* this is general code that prevents zone thickness drNext from * becoming too thin, something that can happen for various bad reasons * HOWEVER we do not want to do this test for certain density laws, * for which very small zone thicknesses are unavoidable * the special cases are: * special density law, * globule density law, * not at carbon +-0 i front * not flucuations command * drMinimum was set in FirstDR to either sdrmin (set drmin) or * some fraction of the initial radius - it is always set * to something non-trivial. * sdrmin is only set wih the "set dr" command */ if( ((strcmp(pressure.chCPres,"DLW1") != 0 && strcmp(pressure.chCPres ,"GLOB") != 0) && xIonFracs[1][5]/xIonFracs[0][5] < 0.05) && (fluct.flong == 0.) ) { /* don't let dr get smaller than drMinimum, if this resets drNext * then code stops with warning that zones got too thin */ if( radius.drNext < drmin.drMinimum ) { radius.drNext = drmin.drMinimum; /* leaving this at true will cause the model to stop with a warning * that the zone thickness is too small */ drminu.lgDrMinUsed = TRUE; } } /* following is to make thickness of model exact * n.b., on last zone, drNext can be NEGATIVE!! * DEPTH was incremented at start of zone calc in newrad, * has been outer edge of zone all throughout */ radius.drNext = MIN2(radius.drNext,(radius.router[IterCnt.iter-1]- radius.depth)*Z); /* this means outer limit exceeded */ if( radius.drNext < 0. ) { drneg.lgDrNeg = TRUE; } else { drneg.lgDrNeg = FALSE; } if( StopCalc.HColStop < 5e29 ) { coleff = phycon.hden*filfac.FillFac; radius.drNext = MIN2( radius.drNext, (StopCalc.HColStop-coldenCom.colden[IPCOLUMNDENSITY-1]-radius.drad*coleff)*Z/coleff); } else if( StopCalc.colpls < 5e29 ) { coleff = h.hii*filfac.FillFac; radius.drNext = MIN2( radius.drNext, (StopCalc.colpls-coldenCom.colden[IPCHII-1]- radius.drad*coleff)*Z/coleff); } else if( StopCalc.colnut < 5e29 ) { coleff = h.hi*filfac.FillFac; /* >>chng 97 oct 30, prevent overflow for very high U models */ if( (StopCalc.colnut - coldenCom.colden[IPCHI-1] - (float)(radius.drad)* coleff) < coleff*radius.drNext ) { radius.drNext = (StopCalc.colnut - coldenCom.colden[IPCHI-1] - radius.drad*coleff)*Z/coleff; } } if( StopCalc.iptnu != NCELL ) { /* end optical depth has been specified */ dt = opac.opac[StopCalc.iptnu-1]*filfac.FillFac; radius.drNext = MIN2(radius.drNext,(StopCalc.tauend-opac.tauabs[0][StopCalc.iptnu-1]- radius.drad*dt)/dt); } if( trace.lgTrace ) { fprintf( ioQQQ, " NextDR chooses next drad=%12.4e; this drad was%12.4e\n", radius.drNext, radius.drad ); } # ifdef DEBUG_FUN fputs( " <->NextDR()\n", debug_fp ); # endif return; }