76 long int level , ipStrong;
78 double thickness_total , drThickness , DepthToGo , AV_to_go;
132 double change_heavy_frac=-1. , change_heavy_frac_big , dr_change_heavy ,
133 frac_change_heavy_frac_big, Efrac_old , Efrac_new;
134 long int ichange_heavy_nelem=-1 , nelem , ion , ichange_heavy_ion=-1;
136 static double OHIIoHI,
142 Old_He_atom_ov_ion = 0,
144 static long int iteration_last=-1;
146 static double BigRadius = 1e30;
184 fprintf(
ioQQQ,
" radius_next called\n" );
203 SaveOHIIoHI = OHIIoHI;
214 x = 1. - atomic_frac /OHIIoHI;
215 if( atomic_frac > 0.05 && atomic_frac < 0.9 )
232 SaveOHIIoHI = OHIIoHI;
238 SaveOHIIoHI = OHIIoHI;
264 Old_He_atom_ov_ion = 0.;
265 drHe1ovHe2 = BigRadius;
270 dr_change_heavy = BigRadius;
271 change_heavy_frac_big = -1.;
272 frac_change_heavy_frac_big = -1.;
274 # define CHANGE_ION_HEAV 0.2f
275 # define CHANGE_ION_HHE 0.15f
299 for( ion=0; ion<=nelem+1; ++ion )
303 if( abundnew > frac_limit && abundold > frac_limit )
313 change_heavy_frac = fabs(abundnew-abundold)/
MIN2(abundold,abundnew);
315 if( (change_heavy_frac > change) && (change_heavy_frac > change_heavy_frac_big) &&
318 ((abundnew-abundold)>0.) &&
319 ((abundold-abundolder)>0.) &&
320 ((abundolder-abundoldest)>0.) )
322 ichange_heavy_nelem = nelem;
323 ichange_heavy_ion = ion;
324 change_heavy_frac_big = change_heavy_frac;
325 frac_change_heavy_frac_big = abundnew;
331 dr_change_heavy =
radius.
drad *
MAX2(0.25, change / change_heavy_frac );
350 drLeiden_hack = BigRadius;
359 drdHdStep = BigRadius;
388 drdHdStep = BigRadius;
409 drConPres = BigRadius;
414 drConPres = BigRadius;
422 drdTdStep = BigRadius;
424 OlddTdStep = dTdStep;
435 x =
MAX2( 0.01 , x );
436 x =
MIN2( 0.05 , x );
442 if( dTdStep*OlddTdStep > 0. )
449 double absdt = fabs(dTdStep);
454 drdTdStep = BigRadius;
457 OlddTdStep = dTdStep;
489 ChkRate(
ipCARBON,&dDRCarbon,&DestOldCarb, &DestNewCarb,&icarstag);
490 ChkRate(
ipNITROGEN,&dDRNitrogen,&DestOldNit, &DestNewNit,&initstag);
491 ChkRate(
ipOXYGEN,&dDROxygen,&DestOldOxy, &DestNewOxy,&ioxystag);
492 ChkRate(
ipIRON,&dDRIron,&DestOldIron, &DestNewIron,&iironstag);
494 dDestRate =
MAX4(dDROxygen,dDRIron,dDRCarbon,dDRNitrogen);
496 if(
fp_equal( dDRCarbon, dDestRate ) )
498 dDestRate = dDRCarbon;
499 DestOld = DestOldCarb;
500 DestNew = DestNewCarb;
502 strcpy( chDestAtom,
"Carbon " );
505 else if(
fp_equal( dDRNitrogen, dDestRate ) )
507 dDestRate = dDRNitrogen;
508 DestOld = DestOldNit;
509 DestNew = DestNewNit;
511 strcpy( chDestAtom,
"Nitrogen" );
514 else if(
fp_equal( dDROxygen, dDestRate ) )
516 dDestRate = dDROxygen;
517 DestOld = DestOldOxy;
518 DestNew = DestNewOxy;
519 strcpy( chDestAtom,
"Oxygen " );
523 else if(
fp_equal( dDRIron, dDestRate ) )
526 DestOld = DestOldIron;
527 DestNew = DestNewIron;
529 strcpy( chDestAtom,
"Iron " );
534 fprintf(
ioQQQ,
" insanity following ChkRate\n" );
573 dVelRelative = fabs(
wind.
windv-OldWindVelocity)/
576 # define WIND_CHNG_VELOCITY_RELATIVE 0.01
586 factor =
MAX2(0.8 , factor );
595 WindAccelDR = BigRadius;
611 WindAccelDR = BigRadius;
623 fprintf(
ioQQQ,
" Globule distance is negative, internal overflow has occured, sorry.\n" );
624 fprintf(
ioQQQ,
" This is routine radius_next, GLBDST is%10.2e\n",
630 if( fac2/factor > 1. +
DNGLOB )
664 fprintf(
ioQQQ,
" dlw insanity in radius_next\n" );
686 else if( dnew/
DNGLOB > 1.0 )
715 DrGrainHeat = BigRadius;
735 else if( level == 2 )
745 else if( level == 3 )
755 else if( level == 4 )
765 else if( level == 5 )
776 fprintf(
ioQQQ,
" PROBLEM radius_next Strong line heating set, but not level.\n" );
784 drLineHeat =
MAX2(1.,TauInwd)*0.4/TauDTau;
788 drLineHeat = BigRadius;
794 drLineHeat = BigRadius;
802 drLineHeat = BigRadius;
810 drH2_heat_cool = BigRadius;
813 Old_H2_heat_cool = 0.;
820 if( H2_heat_cool > 0.1 )
822 dH2_heat_cool = fabs( H2_heat_cool - Old_H2_heat_cool );
823 dH2_heat_cool =
SDIV(dH2_heat_cool);
831 drH2_heat_cool = BigRadius;
839 drH2_abund = BigRadius;
840 dr_mole_abund = BigRadius;
875 dH2_abund =
SDIV(dH2_abund);
880 drH2_abund = BigRadius;
885 dr_mole_abund = BigRadius;
914 if( abund > abund_limit )
916 double drmole_one, relative_change, relative_change_denom;
942 relative_change =
SDIV(relative_change);
945 if( relative_change > 1. )
946 relative_change = 1./relative_change;
953 if( drmole_one < dr_mole_abund )
956 dr_mole_abund = drmole_one;
958 dCO_abund = relative_change;
965 drSolomon_BigH2 =
H2_DR();
987 dr_ne_oscil = BigRadius;
988 dr_te_oscil = BigRadius;
998 for( k=
nzone - 10; k < limit; k++ )
1021 dr_ne_oscil = BigRadius;
1022 dr_te_oscil = BigRadius;
1040 recom_dr_last_iter = BigRadius;
1044 static long int nzone_recom = -1;
1065 recom_dr_last_iter = BigRadius;
1078 dEfrac = fabs(Efrac_old-Efrac_new) / Efrac_new;
1107 drEfrac = BigRadius;
1113 drEfrac = BigRadius;
1129 drDepth = BigRadius;
1134 thickness_total = BigRadius;
1135 DepthToGo = BigRadius;
1139 DepthToGo =
MIN2(DepthToGo ,
1148 DepthToGo =
MIN2(DepthToGo ,
1157 DepthToGo =
MIN2(DepthToGo ,
1166 DepthToGo =
MIN2(DepthToGo ,
1175 DepthToGo =
MIN2(DepthToGo ,
1184 DepthToGo =
MIN2(DepthToGo ,
1192 DepthToGo =
MIN2(DepthToGo ,
1201 DepthToGo =
MIN2(DepthToGo ,
1210 DepthToGo =
MIN2(DepthToGo ,
1216 AV_to_go = BigRadius;
1225 AV_to_go =
MIN2( ave , avp );
1226 if( AV_to_go < 37. )
1228 AV_to_go = pow(10., AV_to_go );
1234 AV_to_go = BigRadius;
1240 if( DepthToGo <= 0. )
1244 else if( DepthToGo < BigRadius )
1249 drThickness =
MIN2( thickness_total/10. , DepthToGo );
1253 drThickness = BigRadius;
1290 drdHdStep = BigRadius;
1298 MIN4( drFluc, GlobDr, DrGrainHeat, dr_change_heavy ) );
1302 MIN3( drdHdStep, drConPres, drTab ),
1303 MIN3( drSolomon_BigH2, drLeiden_hack, recom_dr_last_iter ) );
1305 MIN4( AV_to_go, drEfrac, drHMase, drThickness ),
1306 MIN3( drHe1ovHe2, drH2_heat_cool, drH2_abund ) );
1315 if( nzone <= 1 && radius.drNext > OldDR )
1341 drThermalFront = BigRadius;
1364 fprintf(
ioQQQ,
" radius_next finds insane drNext:%10.2e\n",
1366 fprintf(
ioQQQ,
" all drs follow:%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\n",
1367 drmax, drInter, drLineHeat, winddr, drFluc, GlobDr, SpecDr,
1369 OldH2Abund, drDepth );
1409 fprintf(
punch.
ipDRout,
"#>>>> A temperature failure occured.\n" );
1413 fprintf(
punch.
ipDRout,
"#>>>> A pressure failure occured.\n" );
1425 fprintf(
punch.
ipDRout,
"level 1 line heating,%10.10s TauIn%10.2e%10.2e%10.2e\n",
1429 else if( level == 2 )
1432 fprintf(
punch.
ipDRout,
"level 2 line heating,%10.10s TauIn%10.2e%10.2e%10.2e\n",
1438 fprintf(
ioQQQ,
" insanity pr line heat\n" );
1457 "change in ionization, element %s ion (C scale) %li rel change %.2e ion frac %.2e\n",
1460 change_heavy_frac_big ,
1461 frac_change_heavy_frac_big);
1476 fprintf(
punch.
ipDRout,
"spec den law, new old den%10.2e%10.2e\n",
1482 fprintf(
punch.
ipDRout,
"H maser dTauMase=%10.2e %li %li %li %li\n",
1493 fprintf(
punch.
ipDRout,
"change in He0/He+ ionization, ratio %.2e\n",
1494 Old_He_atom_ov_ion );
1500 fprintf(
punch.
ipDRout,
"change in H2 heating/cooling, d(c,h)/H %.2e\n",
1507 fprintf(
punch.
ipDRout,
"change in H2 abundance, d(H2)/H %.2e\n",
1514 fprintf(
punch.
ipDRout,
"change in heav ele mole abundance, d(mole)/elem %.2e mole=%i=%s\n",
1515 dCO_abund , mole_dr_change ,
COmole[mole_dr_change]->
label);
1521 fprintf(
punch.
ipDRout,
"change in big H2 Solomon rate line opt depth\n");
1527 fprintf(
punch.
ipDRout,
"previous iteration recom logic\n");
1532 fprintf(
punch.
ipDRout,
"change in H ionization fm to%11.3e%11.3e\n",
1533 SaveOHIIoHI, OHIIoHI );
1539 "change in elec den, rel chng:%11.3e, cur %11.3e old=%11.3e\n",
1540 dEfrac , Efrac_old , Efrac_new );
1546 "change in heating; current %10.3e delta=%10.3e\n",
1564 "change in temperature; current= %10.3e, dT/T= %.3f\n",
1586 "optical depth to electron scattering\n" );
1592 "temperature failure.\n" );
1598 "DRMAX; nu opc dr=%10.2e%10.2e%10.2e\n",
1606 "cont inter nu=%10.2e opac=%10.2e\n",
1614 "grain heating nu=%10.2e opac=%10.2e\n",
1621 "Wind, dVelRelative=%.3e windv=%.3e\n",
1629 "density fluctuations\n" );
1635 "GLOB law new dr=%10.2e HDEN=%10.2e\n",
1643 "special law new dr=%10.2e HDEN=%10.2e\n",
1656 " %4ld radius_next keys from insanity %10.2e\n",
1660 " %4ld radius_next keys from insanity %10.2e\n",
1712 "\n DISASTER PROBLEM radius_next finds dr too small and aborts. "
1713 "This is zone %ld iteration %ld. The thickness was %.2e\n",
1718 " If this simulation seems reasonable you can override this limit "
1719 "with the command SET DRMIN %.2e\n\n",
1751 fprintf(
ioQQQ,
" radius_next chooses next drad drNext=%.4e; this drad was%12.4e\n",
1771 Rate_max_nonIonizing;
1780 Rate_max_nonIonizing = 0.;
1781 Freq_nonIonizing = 0.;
1782 Opac_nonIonizing = 0.;
1836 for( i=iplow; i < limit; i++ )
1868 if( Rate_max_nonIonizing < 1e-30 && Opac_Hion >
SMALLFLOAT )
1876 else if( Opac_Hion > Opac_nonIonizing && Rate_max_Hion/Rate_max_nonIonizing > 1e-10 && Opac_Hion > SMALLFLOAT )
1885 *opacm = Opac_nonIonizing;
1886 *freqm = Freq_nonIonizing;
1897 double fluxGrainPeak = -1.;
1898 long int ipGrainPeak = -1;
1899 long int ipGrainPeak2 = -1;
1913 ASSERT( fluxGrainPeak>=0. && ipGrainPeak >=0 );
1926 if( fluxGrainPeak > 0. )
1938 ASSERT( ipGrainPeak2>=0 );
1955 enum {DEBUG_LOC=
false};
1958 fprintf(
ioQQQ,
"conratedebug \t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
1959 Rate_max_nonIonizing,Freq_nonIonizing,Opac_nonIonizing,
1960 Rate_max_Hion,FreqH ,Opac_Hion,*freqm,*opacm
1971 ASSERT( *opacm >= 0. && *freqm >= 0. );
2053 #undef WIND_CHNG_VELOCITY_RELATIVE
2063 double *DestRateOld,
2064 double *DestRateNew,
2081 *DestRateOld = 1e-3;
2082 *DestRateNew = 1e-3;
2097 for( i=0; i < nelem+1; i++ )
2106 for( i=0; i < (nelem); i++ )
2114 OldDest[nelem][i] > 0. )
2121 fprintf(
ioQQQ,
" ChkRate gets insane destruction rate for ion%4ld%4ld%10.2e\n",
2135 if( *dDestRate < dDest )
2140 *DestRateOld = OldDest[nelem][i];