59 static double OHIIoHI,
66 const double BigRadius = 1e30;
78 fprintf(
ioQQQ,
" radius_next called\n" );
81 bool lgFirstCall = (
nzone == 0 );
83 multimap<double,string> drChoice;
87 if( !lgFirstCall && OHIIoHI > 1e-15 &&
94 double x = 1. - atomic_frac/OHIIoHI;
95 if( atomic_frac > 0.05 && atomic_frac < 0.9 )
112 double SaveOHIIoHI = OHIIoHI;
115 oss <<
"change in H ionization old=" << scientific << setprecision(3);
116 oss << SaveOHIIoHI <<
" new=" << OHIIoHI;
117 drChoice.insert( pair<const double,string>( drHion, oss.str() ) );
127 if(
rt.dTauMase < -0.01 )
131 double drHMase =
radius.drad*
MAX2(0.1,-0.2/
rt.dTauMase);
134 oss <<
"H maser dTauMase=" << scientific << setprecision(2) <<
rt.dTauMase <<
" ";
135 oss <<
rt.mas_species <<
" " <<
rt.mas_ion <<
" " <<
rt.mas_hi <<
" " <<
rt.mas_lo;
136 drChoice.insert( pair<const double,string>( drHMase, oss.str() ) );
142 double change_heavy_frac_big = -1.;
143 double frac_change_heavy_frac_big = -1.;
144 const realnum CHANGE_ION_HEAV = 0.2f;
145 const realnum CHANGE_ION_HHE = 0.15f;
148 long ichange_heavy_nelem = -1L;
149 long ichange_heavy_ion = -1L;
150 double dr_change_heavy = BigRadius;
153 if(
dense.lgElmtOn[nelem] )
162 change = CHANGE_ION_HHE;
169 frac_limit =
struc.dr_ionfrac_limit/2.f;
170 change = CHANGE_ION_HEAV;
173 for(
int ion=0; ion<=nelem+1; ++ion )
177 if( abundnew > frac_limit && abundold > frac_limit )
184 double change_heavy_frac = fabs(abundnew-abundold)/
MIN2(abundold,abundnew);
186 if( (change_heavy_frac > change) && (change_heavy_frac > change_heavy_frac_big) &&
189 ((abundnew-abundold)>0.) &&
190 ((abundold-abundolder)>0.) &&
191 ((abundolder-abundoldest)>0.) )
193 ichange_heavy_nelem = nelem;
194 ichange_heavy_ion = ion;
195 change_heavy_frac_big = change_heavy_frac;
196 frac_change_heavy_frac_big = abundnew;
204 dr_change_heavy =
radius.drad *
MAX2(0.25, change / change_heavy_frac );
214 if(ichange_heavy_nelem>=0)
216 oss <<
"change in ionization, element " <<
elementnames.chElementName[ichange_heavy_nelem];
217 oss <<
" ion (C scale) " << ichange_heavy_ion <<
" rel change " << scientific << setprecision(2);
218 oss << change_heavy_frac_big <<
" ion frac " << frac_change_heavy_frac_big;
223 dr_change_heavy = BigRadius;
225 drChoice.insert( pair<const double,string>( dr_change_heavy, oss.str() ) );
235 double drLeiden_hack =
MAX2( 0.02 , 0.05*
rfield.extin_mag_V_point ) /
237 drChoice.insert( pair<const double,string>( drLeiden_hack,
"Leiden hack" ) );
241 static double extin_mag_V_point_old = -1.;
244 const double extin_mag_V_limit = 20.+0.01*extin_mag_V_point_old;
245 double dExtin = (
rfield.extin_mag_V_point - extin_mag_V_point_old)/extin_mag_V_limit;
248 double drExtintion =
radius.drad / dExtin;
250 oss <<
"change in V mag extinction; current=" << scientific << setprecision(3) <<
252 oss << fixed <<
" delta=" << dExtin;
253 drChoice.insert( pair<const double,string>( drExtintion, oss.str() ) );
257 extin_mag_V_point_old =
rfield.extin_mag_V_point;
263 !( strcmp(
dense.chDenseLaw,
"SINE") == 0 &&
dense.lgDenFlucOn ) )
271 drdHdStep =
radius.drad*
MAX2(0.8,0.075/dHdStep);
275 drdHdStep =
radius.drad*
MAX2(0.8,0.100/dHdStep);
283 drdHdStep =
radius.drad*
MAX2(0.8,0.15/dHdStep);
286 oss <<
"change in heating; current=" << scientific << setprecision(3) <<
thermal.htot;
287 oss << fixed <<
" delta=" << dHdStep;
288 drChoice.insert( pair<const double,string>( drdHdStep, oss.str() ) );
295 if( strcmp(
dense.chDenseLaw,
"CPRE") == 0 &&
pressure.lgContRadPresOn )
304 double drConPres = 0.05*
pressure.PresTotlCurr/
306 drChoice.insert( pair<const double,string>( drConPres,
"change in con accel" ) );
312 if( strcmp(
dense.chDenseLaw,
"CPRE") == 0 && (
dark.lgNFW_Set ||
pressure.gravity_symmetry>=0) )
314 if( fabs(
pressure.RhoGravity ) > 0. )
317 ASSERT( drGravPres > 0. );
318 drChoice.insert( pair<const double,string>( drGravPres,
"change in gravitational pressure" ) );
324 double dTdStep = FLT_MAX;
338 if( dTdStep*OlddTdStep > 0. )
345 double absdt = fabs(dTdStep);
346 double drdTdStep =
radius.drad*
MAX2(0.5,x/absdt);
348 oss <<
"change in temperature; current=" << scientific << setprecision(3) <<
phycon.te;
349 oss << fixed <<
" dT/T=" << dTdStep;
350 drChoice.insert( pair<const double,string>( drdTdStep, oss.str() ) );
353 OlddTdStep = dTdStep;
363 double freqm = 0., opacm = 0.;
374 oss <<
"cont inter nu=" << scientific << setprecision(2) << freqm <<
" opac=" << opacm;
375 drChoice.insert( pair<const double,string>( drInter, oss.str() ) );
385 double v = fabs(
wind.windv);
387 double dVelRelative = fabs(
wind.windv-OldWindVelocity)/
390 const double WIND_CHNG_VELOCITY_RELATIVE = 0.01;
396 if( dVelRelative > WIND_CHNG_VELOCITY_RELATIVE &&
nzone>1 )
399 double factor = WIND_CHNG_VELOCITY_RELATIVE / dVelRelative;
401 factor =
MAX2(0.8 , factor );
402 winddr =
radius.drad * factor;
418 oss <<
"Wind, dVelRelative=" << scientific << setprecision(3);
419 oss << dVelRelative <<
" windv=" <<
wind.windv;
420 drChoice.insert( pair<const double,string>( winddr, oss.str() ) );
423 OldWindVelocity =
wind.windv;
425 const double DNGLOB = 0.10;
428 if( strcmp(
dense.chDenseLaw,
"GLOB") == 0 )
432 fprintf(
ioQQQ,
" Globule distance is negative, internal overflow has occured, sorry.\n" );
433 fprintf(
ioQQQ,
" This is routine radius_next, GLBDST is%10.2e\n",
440 double GlobDr = ( fac2/factor > 1. + DNGLOB ) ?
radius.drad*DNGLOB/(fac2/factor - 1.) : BigRadius;
443 oss <<
"GLOB law, HDEN=" << scientific << setprecision(2) <<
dense.gas_phase[
ipHYDROGEN];
444 drChoice.insert( pair<const double,string>( GlobDr, oss.str() ) );
449 if( strncmp(
dense.chDenseLaw ,
"DLW" , 3) == 0 )
452 if( strcmp(
dense.chDenseLaw,
"DLW1") == 0 )
457 else if( strcmp(
dense.chDenseLaw,
"DLW2") == 0 )
462 else if( strcmp(
dense.chDenseLaw,
"DLW3") == 0 )
468 fprintf(
ioQQQ,
" dlw insanity in radius_next\n" );
474 oss <<
"spec den law, new old den " << scientific << setprecision(2);
476 drChoice.insert( pair<const double,string>( drTab, oss.str() ) );
481 if( strcmp(
dense.chDenseLaw,
"DLW1") == 0 )
491 else if( dnew/DNGLOB > 1.0 )
493 SpecDr =
radius.drad/(dnew/DNGLOB);
500 oss <<
"special law, HDEN=" << scientific << setprecision(2) <<
dense.gas_phase[
ipHYDROGEN];
501 drChoice.insert( pair<const double,string>( SpecDr, oss.str() ) );
510 double grfreqm = 0., gropacm = 0.;
516 oss <<
"grain heating nu=" << scientific << setprecision(2) << grfreqm <<
" opac=" << gropacm;
517 drChoice.insert( pair<const double,string>( DrGrainHeat, oss.str() ) );
542 double drLineHeat = ( TauDTau > 0. ) ?
MAX2(1.,TauInwd)*0.4/TauDTau : BigRadius;
544 oss <<
"level " << level <<
" line heating, " <<
chLineLbl(t) <<
" TauIn " << scientific;
545 oss << setprecision(2) << TauInwd <<
" " << t.
Emis().
pump() <<
" " << t.
Emis().
Pesc();
546 drChoice.insert( pair<const double,string>( drLineHeat, oss.str() ) );
552 Old_H2_heat_cool = 0.;
553 else if( !
thermal.lgTemperatureConstant )
557 double H2_heat_cool = (fabs(
hmi.HeatH2Dexc_used)+fabs(
hmi.HeatH2Dish_used)) /
thermal.htot;
559 double dH2_heat_cool = fabs( H2_heat_cool - Old_H2_heat_cool );
560 if( H2_heat_cool > 0.1 && Old_H2_heat_cool>0. && dH2_heat_cool>
SMALLFLOAT )
564 double drH2_heat_cool =
radius.drad*
MAX2(0.5,0.05/dH2_heat_cool);
566 oss <<
"change in H2 heating/cooling, d(c,h)/H " << scientific << setprecision(2);
567 oss << dH2_heat_cool;
568 drChoice.insert( pair<const double,string>( drH2_heat_cool, oss.str() ) );
571 Old_H2_heat_cool = (fabs(
hmi.HeatH2Dexc_used)+fabs(
hmi.HeatH2Dish_used)) /
thermal.htot;
576 int mole_dr_change = -1;
598 double dH2_abund = 2.*fabs(
hmi.H2_total - Old_H2_abund ) /
hmi.H2_total;
609 dH2_abund =
SDIV(dH2_abund);
610 double drH2_abund =
radius.drad*
MAX2(0.2,fac/dH2_abund );
612 oss <<
"change in H2 abundance, d(H2)/H " << scientific << setprecision(2) << dH2_abund;
613 drChoice.insert( pair<const double,string>( drH2_abund, oss.str() ) );
618 double max_change = 0.0;
633 for( molecule::nAtomsMap::iterator atom=
mole_global.list[i]->nAtom.begin();
636 long int nelem = atom->first->el->Z-1;
637 if (
dense.lgElmtOn[nelem])
661 if(
abund > abund_limit )
664 double relative_change =
667 if (relative_change > max_change)
669 max_change = relative_change;
675 if( mole_dr_change >= 0 )
677 double dr_mole_abund =
radius.drad *
MAX2(0.6, 0.05/
SDIV(max_change));
679 oss <<
"change in molecular abundance, old=" << scientific << setprecision(2);
680 oss <<
struc.molecules[mole_dr_change][
nzone-3] <<
" new=" <<
mole.species[mole_dr_change].den;
681 oss <<
" mole=" << mole_dr_change <<
"=" <<
mole_global.list[mole_dr_change]->label;
682 drChoice.insert( pair<const double,string>( dr_mole_abund, oss.str() ) );
688 drChoice.insert( pair<const double,string>( (*diatom)->H2_DR(),
"change in big H2 Solomon rate line opt depth" ) );
698 drChoice.insert( pair<const double,string>( drmax,
"DRMAX" ) );
702 double recom_dr_last_iter = BigRadius;
706 static long int nzone_recom = -1;
713 nzone_recom <
struc.nzonePreviousIteration )
715 ASSERT( nzone_recom>=0 && nzone_recom<
struc.nzonePreviousIteration );
719 while(
struc.depth_last[nzone_recom] <
radius.depth &&
720 nzone_recom <
struc.nzonePreviousIteration-1 )
722 recom_dr_last_iter =
struc.drad_last[nzone_recom]*3.;
723 drChoice.insert( pair<const double,string>( recom_dr_last_iter,
724 "previous iteration recom logic" ) );
737 double dEfrac = fabs(Efrac_old-Efrac_new) / Efrac_new;
745 if(
dense.eden_from_metals > 0.5 )
767 oss <<
"change in elec den, rel chng: " << scientific << setprecision(3) << dEfrac;
768 oss <<
" cur=" << Efrac_new <<
" old=" << Efrac_old;
769 drChoice.insert( pair<const double,string>( drEfrac, oss.str() ) );
777 drChoice.insert( pair<const double,string>(
radius.depth/10.,
"relative depth" ) );
782 double thickness_total = BigRadius;
783 double DepthToGo = BigRadius;
787 DepthToGo =
MIN2(DepthToGo ,
790 thickness_total =
MIN2(thickness_total ,
StopCalc.HColStop / coleff );
796 DepthToGo =
MIN2(DepthToGo ,
798 thickness_total =
MIN2(thickness_total ,
StopCalc.colpls / coleff );
805 DepthToGo =
MIN2(DepthToGo ,
807 thickness_total =
MIN2(thickness_total ,
StopCalc.col_h2 / coleff );
814 DepthToGo =
MIN2(DepthToGo ,
816 thickness_total =
MIN2(thickness_total ,
StopCalc.col_h2_nut / coleff );
819 if(
StopCalc.col_H0_ov_Tspin < 5e29 )
823 DepthToGo =
MIN2(DepthToGo ,
825 thickness_total =
MIN2(thickness_total ,
StopCalc.col_H0_ov_Tspin / coleff );
832 DepthToGo =
MIN2(DepthToGo ,
834 thickness_total =
MIN2(thickness_total ,
StopCalc.col_monoxco / coleff );
840 DepthToGo =
MIN2(DepthToGo ,
842 thickness_total =
MIN2(thickness_total ,
StopCalc.colnut / coleff );
849 DepthToGo =
MIN2(DepthToGo ,
858 DepthToGo =
MIN2(DepthToGo ,
862 if( DepthToGo <= 0. )
866 if( DepthToGo < BigRadius )
868 double depth_min =
MIN2( DepthToGo , thickness_total/100. );
869 DepthToGo =
MAX2( DepthToGo / 10. , depth_min );
875 double drThickness =
MIN2( thickness_total/10. , DepthToGo );
876 drChoice.insert( pair<const double,string>( drThickness,
"depth to go" ) );
882 double AV_to_go = BigRadius;
885 double SAFETY = 1. + 10.*DBL_EPSILON;
891 log10(
rfield.opac_mag_V_extended);
893 log10(
rfield.opac_mag_V_point);
894 AV_to_go =
MIN2( ave , avp );
897 AV_to_go = pow(10., AV_to_go)/
geometry.FillFac;
903 AV_to_go = BigRadius;
906 drChoice.insert( pair<const double,string>( AV_to_go,
"A_V to go" ) );
910 double drSphere =
radius.Radius*0.04;
911 drChoice.insert( pair<const double,string>( drSphere,
"sphericity" ) );
917 if(
thermal.lgTemperatureConstant )
919 drChoice.insert( pair<const double,string>( dRTaue,
"optical depth to electron scattering" ) );
922 if(
dense.flong != 0. )
924 double drFluc = 0.628/2./
dense.flong;
928 drChoice.insert( pair<const double,string>( drFluc,
"density fluctuations" ) );
935 double drStart = (
nzone <= 1 ) ?
radius.drad : BigRadius;
936 drChoice.insert( pair<const double,string>( drStart,
"capped to old DR in first zone" ) );
939 double rfacmax =
radius.lgSdrmaxRel ?
radius.Radius : 1.;
940 drChoice.insert( pair<const double,string>( rfacmax*
radius.sdrmax,
"sdrmax" ) );
962 double drThermalFront =
radius.drad * 0.91;
964 drChoice.insert( pair<const double,string>( drThermalFront,
"thermal front logic" ) );
968 ASSERT( drChoice.size() > 0 );
971 if( drChoice.begin()->first <= 0. )
973 multimap<double,string>::const_iterator ptr = drChoice.begin();
974 fprintf(
ioQQQ,
" radius_next finds insane drNext: %.2e\n", ptr->first );
975 fprintf(
ioQQQ,
" all drs follow:\n" );
976 while( ptr != drChoice.end() )
978 fprintf(
ioQQQ,
" %.2e %s\n", ptr->first, ptr->second.c_str() );
985 if( !
radius.lgFixed && drChoice.begin()->first <
radius.depth*
radius.sdrmin_rel_depth )
988 drChoice.insert( pair<const double,string>(
radius.depth*
radius.sdrmin_rel_depth,
"sdrmin_rel_depth" ) );
992 double rfacmin =
radius.lgSdrminRel ?
radius.Radius : 1.;
993 if( drChoice.begin()->first < rfacmin*
radius.sdrmin )
996 drChoice.insert( pair<const double,string>( rfacmin*
radius.sdrmin,
"sdrmin" ) );
1012 if( strcmp(
dense.chDenseLaw,
"DLW1") != 0 &&
1013 strcmp(
dense.chDenseLaw ,
"GLOB") != 0 &&
1014 dense.flong == 0. &&
1016 drChoice.begin()->first != DepthToGo &&
1018 drChoice.begin()->first != AV_to_go )
1033 radius.lgDrMinUsed =
true;
1040 "\n DISASTER PROBLEM radius_next finds dr too small and aborts. "
1041 "This is zone %ld iteration %ld. The thickness was %.2e\n",
1046 " If this simulation seems reasonable you can override this limit "
1047 "with the command SET DRMIN %.2e\n\n",
1056 const double Z = 1.0001;
1063 drChoice.insert( pair<const double,string>( drOuterRadius,
"outer radius" ) );
1066 radius.drNext = drChoice.begin()->first;
1069 if( drChoice.begin()->second.find(
"H maser" ) != string::npos )
1071 rt.lgMaserSetDR =
true;
1079 if( (
trace.lgTrace &&
trace.lgDrBug) || lgDoPun )
1081 fprintf(
save.ipDRout ,
"%ld\t%.5e\t%.3e\t%.3e\t%s\n",
nzone,
radius.depth,
1082 radius.drNext,
radius.Depth2Go, drChoice.begin()->second.c_str() );
1089 fprintf(
ioQQQ,
" radius_next chooses next drad drNext=%.4e; this drad was %.4e\n",