cloudy
trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*radius_next use adaptive logic to find next zone thickness */ 00004 /*ContRate called by radius_next to find energy of maximum continuum-gas interaction */ 00005 /*GrainRateDr called by radius_next to find grain heating rate dr */ 00014 #include "cddefines.h" 00015 #include "lines_service.h" 00016 #include "iso.h" 00017 #include "geometry.h" 00018 #include "h2.h" 00019 #include "mole.h" 00020 #include "hyperfine.h" 00021 #include "opacity.h" 00022 #include "dense.h" 00023 #include "heavy.h" 00024 #include "grainvar.h" 00025 #include "elementnames.h" 00026 #include "conv.h" 00027 #include "rfield.h" 00028 #include "dynamics.h" 00029 #include "thermal.h" 00030 #include "hmi.h" 00031 #include "coolheavy.h" 00032 #include "timesc.h" 00033 #include "doppvel.h" 00034 #include "stopcalc.h" 00035 #include "colden.h" 00036 #include "phycon.h" 00037 #include "rt.h" 00038 #include "trace.h" 00039 #include "wind.h" 00040 #include "punch.h" 00041 #include "taulines.h" 00042 #include "pressure.h" 00043 #include "iterations.h" 00044 #include "struc.h" 00045 #include "radius.h" 00046 00047 #if 0 00048 /*ChkRate called by radius_next to check how rates of destruction of various species changes */ 00049 STATIC void ChkRate( 00050 /* element number on physical scale */ 00051 long int nelem, 00052 /* change in destruction rate */ 00053 double *dDestRate, 00054 /* old and new destruction rates */ 00055 double *DestRateOld, 00056 double *DestRateNew, 00057 /* stage of ionization on the physical scale */ 00058 long int *istage); 00059 #endif 00060 00061 /*ContRate called by radius_next to find energy of maximum continuum-gas interaction */ 00062 STATIC void ContRate(double *freqm, 00063 double *opacm); 00064 00065 /*GrainRateDr called by radius_next to find grain heating rate dr */ 00066 STATIC void GrainRateDr(double *grfreqm, 00067 double *gropacm); 00068 00069 /*radius_next use adaptive logic to find next zone thickness 00070 * return 0 if ok, 1 for abort */ 00071 int radius_next(void) 00072 { 00073 char chLbl[11]; 00074 bool lgDoPun; 00075 00076 long int level , ipStrong; 00077 00078 double thickness_total , drThickness , DepthToGo , AV_to_go; 00079 int mole_dr_change; 00080 00081 double DrGrainHeat, 00082 GlobDr, 00083 SaveOHIIoHI, 00084 SpecDr, 00085 Strong, 00086 TauDTau, 00087 TauInwd, 00088 drSolomon_BigH2 , 00089 dEfrac, 00090 dHdStep, 00091 dRTaue, 00092 dTdStep, 00093 dnew, 00094 drConPres, 00095 drH2_heat_cool = 0. , 00096 dH2_heat_cool = 0., 00097 drH2_abund = 0. , 00098 dr_mole_abund = 0., 00099 dH2_abund=0., 00100 dCO_abund=0., 00101 drDepth, 00102 drDest, 00103 drEfrac, 00104 drFail, 00105 drFluc, 00106 drHMase, 00107 drHe1ovHe2, 00108 drHion, 00109 drInter, 00110 drLeiden_hack , 00111 drLineHeat, 00112 drSphere, 00113 drTab, 00114 drdHdStep, 00115 drdTdStep, 00116 drThermalFront , 00117 drmax, 00118 dVelRelative, 00119 fac2, 00120 factor, 00121 freqm, 00122 grfreqm=0., 00123 gropacm=0., 00124 hdnew, 00125 opacm, 00126 recom_dr_last_iter , 00127 OldDR , 00128 winddr, 00129 WindAccelDR, 00130 x; 00131 00132 double change_heavy_frac=-1. , change_heavy_frac_big , dr_change_heavy , 00133 frac_change_heavy_frac_big, Efrac_old , Efrac_new; 00134 long int ichange_heavy_nelem=-1 , nelem , ion , ichange_heavy_ion=-1; 00135 00136 static double OHIIoHI, 00137 OldHeat = 0., 00138 OldTe = 0., 00139 OlddTdStep = 0., 00140 OldH2Abund=0., 00141 OldWindVelocity=0., 00142 Old_He_atom_ov_ion = 0, 00143 Old_H2_heat_cool; 00144 static long int iteration_last=-1; 00145 00146 static double BigRadius = 1e30; 00147 bool lgFirstCall; 00148 00149 DEBUG_ENTRY( "radius_next()" ); 00150 00151 00152 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 00153 * 00154 * changes in logic 00155 * 95 oct 19, drSphere now 3% of radius, was 2%, fewer zone 00156 * 95 oct 19, subtracted grain opacity from total opacity used 00157 * to get thickness in routine ContRate 00158 * 00159 *>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 00160 * 00161 * free statement labels >= 13 00162 * 00163 *----------------------------------------------------------------------- 00164 * 00165 * this sub determines the thickness of the next zone 00166 * if is called one time for each zone 00167 * flag lgNxtDROff is true if this is initialization of radius_next, 00168 * is false if we are to use logic to find dr 00169 * 00170 *----------------------------------------------------------------------- */ 00171 00172 /* >>chng 03 sep 21 - decide whether this is the first call */ 00173 if( (iteration != iteration_last) && (nzone==0) ) 00174 { 00175 /* this is the first call in this iteration */ 00176 iteration_last = iteration; 00177 lgFirstCall = true; 00178 } 00179 else 00180 lgFirstCall = false; 00181 00182 if( trace.lgTrace ) 00183 { 00184 fprintf( ioQQQ, " radius_next called\n" ); 00185 } 00186 00187 /* save current dr */ 00188 OldDR = radius.drad; 00189 00190 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/ 00191 /* 1 '' radius_next keys from change in H ionization'',e11.3)') 00192 * check on change in hydrogen ionizaiton */ 00193 if( lgFirstCall ) 00194 { 00195 if( dense.xIonDense[ipHYDROGEN][1] > 0. && dense.xIonDense[ipHYDROGEN][0] > 0. ) 00196 { 00197 OHIIoHI = dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0]; 00198 } 00199 else 00200 { 00201 OHIIoHI = 0.; 00202 } 00203 SaveOHIIoHI = OHIIoHI; 00204 drHion = BigRadius; 00205 /* else if(hii.gt.0. .and. hi.gt.0. .and. OHIIoHI.gt.0. ) then 00206 * >>chng 97 jul 9, for deep in PDR vastly now ionz H slowed down works */ 00207 } 00208 00209 else if( (dense.xIonDense[ipHYDROGEN][1] > 0. && dense.xIonDense[ipHYDROGEN][0] > 0.) && OHIIoHI > 1e-15 ) 00210 { 00211 double atomic_frac = (dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0]); 00212 /* ratio of current HII/HI to old value - < 1 when becoming more neutral */ 00213 /* this is now change in atomic fraction */ 00214 x = 1. - atomic_frac /OHIIoHI; 00215 if( atomic_frac > 0.05 && atomic_frac < 0.9 ) 00216 { 00217 /* >>chng 96 feb 23 from 0.7 to 0.3 because punched through H i-front 00218 * >>chng 97 may 5, from 0.3 to 0.2 since punched through H i-front */ 00219 /* >>chng 02 jun 05 from 0.2 to 0.05 poorly resolved i-front, also added two-branch logic*/ 00220 drHion = radius.drad*MAX2( 0.2 , 0.05/MAX2(1e-10,x) ); 00221 } 00222 else if( x > 0. ) 00223 { 00224 /* >>chng 96 feb 23 from 0.7 to 0.3 because punched through H i-front 00225 * >>chng 97 may 5, from 0.3 to 0.2 since punched through H i-front */ 00226 drHion = radius.drad*MAX2( 0.2 , 0.2/MAX2(1e-10,x) ); 00227 } 00228 else 00229 { 00230 drHion = BigRadius; 00231 } 00232 SaveOHIIoHI = OHIIoHI; 00233 OHIIoHI = dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0]; 00234 } 00235 00236 else 00237 { 00238 SaveOHIIoHI = OHIIoHI; 00239 if( dense.xIonDense[ipHYDROGEN][1] > 0. && dense.xIonDense[ipHYDROGEN][0] > 0. ) 00240 { 00241 OHIIoHI = dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0]; 00242 } 00243 else 00244 { 00245 OHIIoHI = 0.; 00246 } 00247 drHion = BigRadius; 00248 } 00249 00250 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/ 00251 /* "radius_next keys from H maser, dt, ij=" possible hydrogen maser action */ 00252 if( rt.dTauMase < -0.01 ) 00253 { 00254 /* maser so powerful, must limit inc in tay 00255 * >>chng 96 aug 08, from 0.3 to 0.1 due to large maser crash */ 00256 drHMase = radius.drad*MAX2(0.1,-0.2/rt.dTauMase); 00257 } 00258 else 00259 { 00260 drHMase = BigRadius; 00261 } 00262 00263 /* >>chng 03 nov 09, try doing he in following, not above */ 00264 Old_He_atom_ov_ion = 0.; 00265 drHe1ovHe2 = BigRadius; 00266 00267 /* check on N0 - 1 ionization changes, 00268 * >>chng 03 jun 06, add this test due to smashing into H ifront in blr89.in 00269 */ 00270 dr_change_heavy = BigRadius; 00271 change_heavy_frac_big = -1.; 00272 frac_change_heavy_frac_big = -1.; 00273 /* >chng 03 jun 09, from 0.05 to 0.1, initial tests with zoning */ 00274 # define CHANGE_ION_HEAV 0.2f 00275 # define CHANGE_ION_HHE 0.15f 00276 if( nzone > 4 ) 00277 { 00278 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem ) 00279 { 00280 if( dense.lgElmtOn[nelem] ) 00281 { 00282 realnum change; 00283 /* this is the limit to the ionization we will check - 00284 * also used in prt_comment to check on whether oscillations occurred */ 00285 realnum frac_limit; 00286 if( nelem<=ipHELIUM || nelem==ipCARBON ) 00287 { 00288 frac_limit = 1e-4f; 00289 change = CHANGE_ION_HHE; 00290 } 00291 else 00292 { 00293 /* this var is used to print warnings, 00294 * make sure we converge below it */ 00295 frac_limit = struc.dr_ionfrac_limit/2.f; 00296 change = CHANGE_ION_HEAV; 00297 } 00298 /* >>chng 03 dec 09, go up to full range of ion, not just =2 */ 00299 for( ion=0; ion<=nelem+1; ++ion ) 00300 { 00301 realnum abundnew = dense.xIonDense[nelem][ion]/SDIV( dense.gas_phase[nelem]); 00302 realnum abundold = struc.xIonDense[nelem][ion][nzone-3]/SDIV( struc.gas_phase[nelem][nzone-3]); 00303 if( abundnew > frac_limit && abundold > frac_limit ) 00304 { 00305 /* NB must make sure this test is not done when nzone-x <0 */ 00306 /* -2 because previous zone, and nzone is off by one (on physical, not C, scale) */ 00307 /*realnum abundold = struc.xIonDense[nelem][ion][nzone-3]/SDIV( struc.gas_phase[nelem][nzone-3]);*/ 00308 realnum abundolder = struc.xIonDense[nelem][ion][nzone-4]/SDIV( struc.gas_phase[nelem][nzone-4]); 00309 realnum abundoldest = struc.xIonDense[nelem][ion][nzone-5]/SDIV( struc.gas_phase[nelem][nzone-5]); 00310 /* this is fractional change */ 00311 /* >>chng 04 feb 28, take min of old and new abund, to try to pick up 00312 * rapid changing Ca+ sooner */ 00313 change_heavy_frac = fabs(abundnew-abundold)/MIN2(abundold,abundnew); 00314 /* want fractional change to be less than this factor */ 00315 if( (change_heavy_frac > change) && (change_heavy_frac > change_heavy_frac_big) && 00316 /* >>chng 03 dec 07, add test that abund is not oscillating */ 00317 /* also test that abundance is increasing - we are headed into a front */ 00318 ((abundnew-abundold)>0.) && 00319 ((abundold-abundolder)>0.) && 00320 ((abundolder-abundoldest)>0.) ) 00321 { 00322 ichange_heavy_nelem = nelem; 00323 ichange_heavy_ion = ion; 00324 change_heavy_frac_big = change_heavy_frac; 00325 frac_change_heavy_frac_big = abundnew; 00326 /* >>chng 03 dec 06, from min of 0.5 to min of 0.25, crash into He i-front 00327 * in hizqso.in */ 00328 /* >>chng 04 mar 03, min had become 0.1, forced oscillations in nova.in 00329 * in silicon, zoning changed greatly, causing change in diffuse lin 00330 * pumping. put back to 0.25 */ 00331 dr_change_heavy = radius.drad * MAX2(0.25, change / change_heavy_frac ); 00332 } 00333 } 00334 } 00335 } 00336 } 00337 } 00338 00339 /* if Leiden hacks are on then use increase in dust extinction as control 00340 * >>chng 05 aug 13, add this */ 00341 if(!co.lgUMISTrates) 00342 { 00343 /* Draine field is only defined over narrow range in FUV - must not let change 00344 * in extinction become too large - 00345 * prefactor is change in optical depth */ 00346 drLeiden_hack = MAX2( 0.02 , 0.05*rfield.extin_mag_V_point) / SDIV( geometry.FillFac * rfield.opac_mag_V_point ); 00347 } 00348 else 00349 { 00350 drLeiden_hack = BigRadius; 00351 } 00352 /* >>chng 04 feb 15, kill this block - not used */ 00353 00354 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/ 00355 /* check how heating is changing 00356 * '' radius_next keys from change in heating; current, delta='', */ 00357 if( nzone <= 1 || thermal.lgTemperatureConstant ) 00358 { 00359 drdHdStep = BigRadius; 00360 dHdStep = FLT_MAX; 00361 } 00362 else 00363 { 00364 dHdStep = fabs(thermal.htot-OldHeat)/thermal.htot; 00365 if( dHdStep > 0. ) 00366 { 00367 if( dense.gas_phase[ipHYDROGEN] >= 1e13 ) 00368 { 00369 /* drdHdStep = drad * MAX( 0.8 , 0.05/dHdStep ) */ 00370 drdHdStep = radius.drad*MAX2(0.8,0.075/dHdStep); 00371 } 00372 else if( dense.gas_phase[ipHYDROGEN] >= 1e11 ) 00373 { 00374 /* drdHdStep = drad * MAX( 0.8 , 0.075/dHdStep ) */ 00375 drdHdStep = radius.drad*MAX2(0.8,0.100/dHdStep); 00376 } 00377 else 00378 { 00379 /* changed from .15 to .12 for outer edge of coolhii too steep dT 00380 * changed to .10 for symmetry, big change in some rates, 95nov14 00381 * changed from .10 to .125 since parispn seemed to waste zones 00382 * >>chng 96 may 21, from .125 to .15 since pn's still waste zones */ 00383 drdHdStep = radius.drad*MAX2(0.8,0.15/dHdStep); 00384 } 00385 } 00386 else 00387 { 00388 drdHdStep = BigRadius; 00389 } 00390 } 00391 OldHeat = thermal.htot; 00392 00393 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/ 00394 /* pressure due to incident continuum if in eos */ 00395 if( strcmp(dense.chDenseLaw,"CPRE") == 0 && pressure.lgContRadPresOn ) 00396 { 00397 if( nzone > 2 && pressure.pinzon > 0. ) 00398 { 00399 /* pinzon is pressrue from acceleration onto previos zone 00400 * want this less than some fraction of total pressure */ 00401 /* >>chng 06 feb 01, change from init pres to current total pressure 00402 * in const press high U ulirgs current pressure may be quite larger 00403 * than init pressure due to continuum absorption */ 00404 drConPres = 0.05*pressure.PresTotlCurr/(wind.AccelTot* 00405 dense.xMassDensity*geometry.FillFac); 00406 } 00407 else 00408 { 00409 drConPres = BigRadius; 00410 } 00411 } 00412 else 00413 { 00414 drConPres = BigRadius; 00415 } 00416 00417 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/ 00418 /* check how temperature is changing 00419 * 1 '' radius_next keys from change in temperature; current, delta='', */ 00420 if( nzone <= 1 ) 00421 { 00422 drdTdStep = BigRadius; 00423 dTdStep = FLT_MAX; 00424 OlddTdStep = dTdStep; 00425 } 00426 else 00427 { 00428 /* change in temperature; current= */ 00429 dTdStep = (phycon.te-OldTe)/phycon.te; 00430 /* >>chng 02 dec 08, desired change in temperature per zone must not 00431 * be smaller than allower error in temperature. For now use relative error 00432 * in heating - - cooling balance. Better would be to also use c-h deriv wrt t 00433 * to get slope */ 00434 x = conv.HeatCoolRelErrorAllowed*2.; 00435 x = MAX2( 0.01 , x ); 00436 x = MIN2( 0.05 , x ); 00437 /* >>chng 02 dec 11 rjrw change back to 0.03 -- improve dynamics.dRad criterion instead */ 00438 x = 0.03; 00439 /* >>chng 02 dec 07, do not do this if there is mild te jitter, 00440 * so that dT is changing sign - this happens 00441 * in ism.ini, very stable temperature with slight noise up and down */ 00442 if( dTdStep*OlddTdStep > 0. ) 00443 { 00444 /* >>chng 96 may 30, variable depending on temp 00445 * >>chng 96 may 31, allow 0.7 smaller, was 0.8 00446 * >>chng 97 may 05, from 0.7 to 0.5 stop from punching through thermal front 00447 * >>chng 04 mar 30, from 0.7 to 0.5 stop from punching through thermal front, 00448 * for some reason factor was 0.7, not 0.5 as per previous change */ 00449 double absdt = fabs(dTdStep); 00450 drdTdStep = radius.drad*MAX2(0.5,x/absdt); 00451 } 00452 else 00453 { 00454 drdTdStep = BigRadius; 00455 } 00456 } 00457 OlddTdStep = dTdStep; 00458 OldTe = phycon.te; 00459 00460 /* >>chng 02 oct 06, only check on opacity - interaction if not 00461 * constant temperature - there were constant temperature models that 00462 * extended to infinity but hung with last few photons and this test. 00463 * better to ignore this test which is really for thermal models */ 00464 /* >>chng 06 mar 20, do not call if recombination logic in place */ 00465 if( !thermal.lgTemperatureConstant && !dynamics.lgRecom ) 00466 { 00467 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/ 00468 /* find freq where opacity largest and interaction rate is fastest 00469 * "cont inter nu=%10.2e opac=%10.2e\n" */ 00470 ContRate(&freqm,&opacm); 00471 00472 /* use optical depth at max interaction energy 00473 * >>chng 96 jun 06 was drChange=0.15 changed to 0.3 for high Z models 00474 * taking too many zones 00475 * drInter = drChange / MAX(1e-30,opacm*FillFac) */ 00476 00477 drInter = 0.3/MAX2(1e-30,opacm*geometry.FillFac*geometry.DirectionalCosin); 00478 } 00479 else 00480 { 00481 drInter = BigRadius; 00482 freqm = 0.; 00483 opacm = 1.; 00484 } 00485 00486 # if 0 00487 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/ 00488 /* check on changes in destruction rates for various atoms */ 00489 ChkRate(ipCARBON,&dDRCarbon,&DestOldCarb, &DestNewCarb,&icarstag); 00490 ChkRate(ipNITROGEN,&dDRNitrogen,&DestOldNit, &DestNewNit,&initstag); 00491 ChkRate(ipOXYGEN,&dDROxygen,&DestOldOxy, &DestNewOxy,&ioxystag); 00492 ChkRate(ipIRON,&dDRIron,&DestOldIron, &DestNewIron,&iironstag); 00493 00494 dDestRate = MAX4(dDROxygen,dDRIron,dDRCarbon,dDRNitrogen); 00495 00496 if( fp_equal( dDRCarbon, dDestRate ) ) 00497 { 00498 dDestRate = dDRCarbon; 00499 DestOld = DestOldCarb; 00500 DestNew = DestNewCarb; 00501 istage = icarstag; 00502 strcpy( chDestAtom, "Carbon " ); 00503 } 00504 00505 else if( fp_equal( dDRNitrogen, dDestRate ) ) 00506 { 00507 dDestRate = dDRNitrogen; 00508 DestOld = DestOldNit; 00509 DestNew = DestNewNit; 00510 istage = initstag; 00511 strcpy( chDestAtom, "Nitrogen" ); 00512 } 00513 00514 else if( fp_equal( dDROxygen, dDestRate ) ) 00515 { 00516 dDestRate = dDROxygen; 00517 DestOld = DestOldOxy; 00518 DestNew = DestNewOxy; 00519 strcpy( chDestAtom, "Oxygen " ); 00520 istage = ioxystag; 00521 } 00522 00523 else if( fp_equal( dDRIron, dDestRate ) ) 00524 { 00525 dDestRate = dDRIron; 00526 DestOld = DestOldIron; 00527 DestNew = DestNewIron; 00528 istage = iironstag; 00529 strcpy( chDestAtom, "Iron " ); 00530 } 00531 00532 else 00533 { 00534 fprintf( ioQQQ, " insanity following ChkRate\n" ); 00535 ShowMe(); 00536 cdEXIT(EXIT_FAILURE); 00537 } 00538 00539 /* radius_next keys from change in dest rates, atom= */ 00540 if( dDestRate > 0. ) 00541 { 00542 /* if( te.gt.40 000. ) then 00543 * added different accuracy for hot gas since tend to jump over 00544 * big te range for small change in heat (intrinsically unstable) 00545 * drDest = drad * MAX( 0.5 , 0.10/dDestRate ) 00546 * else 00547 * was drChange, changed to .15 for parishii go through HeII=HeI I front 00548 * drDest = drad * MAX( 0.5 , 0.15/dDestRate ) 00549 * >>chng 95 dec 18 from above to below to stop oscillations 00550 * >>chng 95 dec 27 from min of .5 to .75 to stop zone size from changing rapidly 00551 * drDest = drad * MAX( 0.8 , 0.20 /dDestRate ) 00552 * >>chng 96 jan 14 from .2 to .25 to use less zones 00553 * >>chng 96 may 30 from min of 0.8 to 0.5 to prevent crashing into He i-front */ 00554 drDest = radius.drad*MAX2(0.5,0.20/dDestRate); 00555 /* endif */ 00556 } 00557 else 00558 { 00559 drDest = BigRadius; 00560 } 00561 # endif 00562 drDest = BigRadius;/*03 dec 12 is this needed? */ 00563 00564 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/ 00565 /* check whether change in wind velocity constrains DRAD 00566 * WJH 22 May 2004: disable when we are near the sonic point since 00567 * the velocity may be jumping all over the place but we just want 00568 * to push through it as quickly as we can */ 00569 if( wind.windv!=0. && !pressure.lgSonicPoint && !pressure.lgStrongDLimbo ) 00570 { 00571 double v = fabs(wind.windv); 00572 /* this is relative change in velocity */ 00573 dVelRelative = fabs(wind.windv-OldWindVelocity)/ 00574 MAX2(v,0.1*timesc.sound_speed_isothermal); 00575 00576 # define WIND_CHNG_VELOCITY_RELATIVE 0.01 00577 /*fprintf(ioQQQ,"DEBUG rad %.3f vel %.2e dVelRelative %.2e", 00578 log10(radius.Radius) , wind.windv , dVelRelative );*/ 00579 00580 /* do not use this on first zone since do not have old velocity */ 00581 if( dVelRelative > WIND_CHNG_VELOCITY_RELATIVE && nzone>1 ) 00582 { 00583 /* factor less than one if change larger than WIND_CHNG_VELOCITY_RELATIVE */ 00584 double factor = WIND_CHNG_VELOCITY_RELATIVE / dVelRelative; 00585 /*fprintf(ioQQQ," DEBUG factor %.2e", factor );*/ 00586 factor = MAX2(0.8 , factor ); 00587 winddr = radius.drad * factor; 00588 } 00589 else 00590 { 00591 winddr = BigRadius; 00592 } 00593 /*fprintf(ioQQQ," DEBUG \n" );*/ 00594 00595 WindAccelDR = BigRadius; 00596 00597 /* set dr from advective term, 00598 * the 1/20 came from looking at one set of structure plots */ 00599 if( dynamics.lgAdvection && iteration > dynamics.n_initial_relax+1) 00600 { 00601 /* >>chng 02 dec 11, from 5 to 20 */ 00602 winddr = MIN2( winddr , dynamics.dRad / 20. ); 00603 /*>>chng 04 oct 06, set dVelRelative to dynamics.dRad since dVelRelative is printed as part 00604 * of reason for choosing this criteria, want to reflect valid reason. */ 00605 dVelRelative = dynamics.dRad; 00606 } 00607 } 00608 else 00609 { 00610 winddr = BigRadius; 00611 WindAccelDR = BigRadius; 00612 dVelRelative = 0.; 00613 } 00614 /* remember old velocity */ 00615 OldWindVelocity = wind.windv; 00616 00617 /* inside out globule */ 00618 if( strcmp(dense.chDenseLaw,"GLOB") == 0 ) 00619 { 00620 # define DNGLOB 0.10 00621 if( radius.glbdst < 0. ) 00622 { 00623 fprintf( ioQQQ, " Globule distance is negative, internal overflow has occured, sorry.\n" ); 00624 fprintf( ioQQQ, " This is routine radius_next, GLBDST is%10.2e\n", 00625 radius.glbdst ); 00626 cdEXIT(EXIT_FAILURE); 00627 } 00628 factor = radius.glbden*pow(radius.glbrad/radius.glbdst,radius.glbpow); 00629 fac2 = radius.glbden*pow(radius.glbrad/(radius.glbdst - (realnum)radius.drad),radius.glbpow); 00630 if( fac2/factor > 1. + DNGLOB ) 00631 { 00632 /* DNGLOB is relative change in density allowed this zone, 0.10 */ 00633 GlobDr = radius.drad*DNGLOB/(fac2/factor - 1.); 00634 } 00635 else 00636 { 00637 GlobDr = BigRadius; 00638 } 00639 /* GlobDr = GLBDST * DNGLOB * (GLBRAD/GLBDST)**(-GLBPOW) / GLBPOW */ 00640 GlobDr = MIN2(GlobDr,radius.glbdst/20.); 00641 } 00642 else 00643 { 00644 GlobDr = BigRadius; 00645 } 00646 00647 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/ 00648 hdnew = 0.; 00649 if( strncmp( dense.chDenseLaw , "DLW" , 3) == 0 ) 00650 { 00651 /* one of the special density laws, first get density at possible next dr */ 00652 if( strcmp(dense.chDenseLaw,"DLW1") == 0 ) 00653 { 00654 hdnew = dense_fabden(radius.Radius+radius.drad,radius.depth+ 00655 radius.drad); 00656 } 00657 else if( strcmp(dense.chDenseLaw,"DLW2") == 0 ) 00658 { 00659 hdnew = dense_tabden(radius.Radius+radius.drad,radius.depth+ 00660 radius.drad); 00661 } 00662 else 00663 { 00664 fprintf( ioQQQ, " dlw insanity in radius_next\n" ); 00665 cdEXIT(EXIT_FAILURE); 00666 } 00667 drTab = fabs(hdnew-dense.gas_phase[ipHYDROGEN])/MAX2(hdnew,dense.gas_phase[ipHYDROGEN]); 00668 drTab = radius.drad*MAX2(0.2,0.10/MAX2(0.01,drTab)); 00669 } 00670 else 00671 { 00672 drTab = BigRadius; 00673 } 00674 00675 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/ 00676 /* special density law */ 00677 if( strcmp(dense.chDenseLaw,"DLW1") == 0 ) 00678 { 00679 dnew = fabs(dense_fabden(radius.Radius+radius.drad,radius.depth+ 00680 radius.drad)/dense.gas_phase[ipHYDROGEN]-1.); 00681 /* DNGLOB is relative change in density allowed this zone, 0.10 */ 00682 if( dnew == 0. ) 00683 { 00684 SpecDr = radius.drad*3.; 00685 } 00686 else if( dnew/DNGLOB > 1.0 ) 00687 { 00688 SpecDr = radius.drad/(dnew/DNGLOB); 00689 } 00690 else 00691 { 00692 SpecDr = BigRadius; 00693 } 00694 } 00695 else 00696 { 00697 SpecDr = BigRadius; 00698 } 00699 00700 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/ 00701 /* check grain line heating dominates 00702 * this is important in PDR and HII region calculations 00703 * >>chng 97 jul 03, added following check */ 00704 if( thermal.heating[0][13]/thermal.htot > 0.2 ) 00705 { 00706 /* >>chng 01 jan 03, following returns 0 when NO light at all, 00707 * had failed with botched assert */ 00708 GrainRateDr(&grfreqm,&gropacm); 00709 DrGrainHeat = 1.0/MAX2(1e-30,gropacm*geometry.FillFac*geometry.DirectionalCosin); 00710 } 00711 else 00712 { 00713 gropacm = 0.; 00714 grfreqm = 0.; 00715 DrGrainHeat = BigRadius; 00716 } 00717 00718 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/ 00719 /* check if line heating dominates 00720 * this is important in high metallicity models */ 00721 if( thermal.heating[0][22]/thermal.htot > 0.2 ) 00722 { 00723 FndLineHt(&level,&ipStrong,&Strong); 00724 if( Strong/thermal.htot > 0.1 ) 00725 { 00726 if( level == 1 ) 00727 { 00728 /* a level1 line was the heat source (usual case) 00729 * drLineHeat = MAX(1.0,TauLines(ipLnTauIn,ipStrong)*0.2) / 00730 * 1 TauLines(ipLnDTau,ipStrong) */ 00731 TauInwd = TauLines[ipStrong].Emis->TauIn; 00732 TauDTau = TauLines[ipStrong].Emis->PopOpc * TauLines[ipStrong].Emis->opacity / 00733 DoppVel.doppler[TauLines[ipStrong].Hi->nelem-1]; 00734 } 00735 else if( level == 2 ) 00736 { 00737 /* a atom_level2 line was the heat source 00738 * (bad case since not as well treated) 00739 * drLineHeat = MAX(1.0,WindLine(ipLnTauIn,ipStrong)*0.2) / 00740 * 1 WindLine(ipLnDTau,ipStrong) */ 00741 TauInwd = TauLine2[ipStrong].Emis->TauIn; 00742 TauDTau = TauLine2[ipStrong].Emis->PopOpc * TauLine2[ipStrong].Emis->opacity / 00743 DoppVel.doppler[TauLine2[ipStrong].Hi->nelem-1]; 00744 } 00745 else if( level == 3 ) 00746 { 00747 /* a 12CO line was the heat source 00748 * (bad case since not as well treated) 00749 * drLineHeat = MAX(1.0,WindLine(ipLnTauIn,ipStrong)*0.2) / 00750 * 1 WindLine(ipLnDTau,ipStrong) */ 00751 TauInwd = C12O16Rotate[ipStrong].Emis->TauIn; 00752 TauDTau = C12O16Rotate[ipStrong].Emis->PopOpc * C12O16Rotate[ipStrong].Emis->opacity / 00753 DoppVel.doppler[C12O16Rotate[ipStrong].Hi->nelem-1]; 00754 } 00755 else if( level == 4 ) 00756 { 00757 /* a 13CO line was the heat source 00758 * (bad case since not as well treated) 00759 * drLineHeat = MAX(1.0,WindLine(ipLnTauIn,ipStrong)*0.2) / 00760 * 1 WindLine(ipLnDTau,ipStrong) */ 00761 TauInwd = C13O16Rotate[ipStrong].Emis->TauIn; 00762 TauDTau = C13O16Rotate[ipStrong].Emis->PopOpc * C13O16Rotate[ipStrong].Emis->opacity / 00763 DoppVel.doppler[C13O16Rotate[ipStrong].Hi->nelem-1]; 00764 } 00765 else if( level == 5 ) 00766 { 00767 /* >>chng 03 dec 07, this branch had been left off, caught by Hiroaki Oyaizu */ 00768 /* a hyperfine transition */ 00769 TauInwd = HFLines[ipStrong].Emis->TauIn; 00770 TauDTau = HFLines[ipStrong].Emis->PopOpc * HFLines[ipStrong].Emis->opacity / 00771 DoppVel.doppler[HFLines[ipStrong].Hi->nelem-1]; 00772 } 00773 else 00774 { 00775 /* this is insane, since Strong was set, but not level */ 00776 fprintf( ioQQQ, " PROBLEM radius_next Strong line heating set, but not level.\n" ); 00777 TotalInsanity(); 00778 } 00779 /* in following logic cannot use usual inverse opacity, 00780 * since line heating competes with escape probability, 00781 * so is effective at surprising optical depths */ 00782 if( TauDTau > 0. ) 00783 { 00784 drLineHeat = MAX2(1.,TauInwd)*0.4/TauDTau; 00785 } 00786 else 00787 { 00788 drLineHeat = BigRadius; 00789 } 00790 } 00791 else 00792 { 00793 TauInwd = 0.; 00794 drLineHeat = BigRadius; 00795 ipStrong = 0; 00796 Strong = 0.; 00797 } 00798 } 00799 else 00800 { 00801 TauInwd = 0.; 00802 drLineHeat = BigRadius; 00803 ipStrong = 0; 00804 level = 0; 00805 Strong = 0.; 00806 } 00807 00808 /* >>chng 03 mar 03, add this logic */ 00809 /* do not let change in cooling/heating due to H2 become too large */ 00810 drH2_heat_cool = BigRadius; 00811 if( lgFirstCall ) 00812 { 00813 Old_H2_heat_cool = 0.; 00814 } 00815 else if( !thermal.lgTemperatureConstant ) 00816 { 00817 /* this is case where temperature has not been set */ 00818 /* compare total heating - cooling due to h2 with total due to everything */ 00819 double H2_heat_cool = (fabs(hmi.HeatH2Dexc_used)+fabs(hmi.HeatH2Dish_used)) / thermal.htot; 00820 if( H2_heat_cool > 0.1 ) 00821 { 00822 dH2_heat_cool = fabs( H2_heat_cool - Old_H2_heat_cool ); 00823 dH2_heat_cool = SDIV(dH2_heat_cool); 00824 /* >>chng 05 oct 27, had been taking 20% of original radius - this caused zoning 00825 * to become very fine and may not have been needed - change from 0.2 to 0.5 */ 00826 /*drH2_heat_cool = radius.drad*MAX2(0.2,0.05/dH2_heat_cool);*/ 00827 drH2_heat_cool = radius.drad*MAX2(0.5,0.05/dH2_heat_cool); 00828 } 00829 else 00830 { 00831 drH2_heat_cool = BigRadius; 00832 } 00833 } 00834 Old_H2_heat_cool = (fabs(hmi.HeatH2Dexc_used)+fabs(hmi.HeatH2Dish_used)) / thermal.htot; 00835 00836 /* >>chng 03 mar 04, add this logic */ 00837 /* do not let change in H2 and heavy element molecular abundances become too large 00838 * "change in heav ele mole abundance, d(mole)/elem" */ 00839 drH2_abund = BigRadius; 00840 dr_mole_abund = BigRadius; 00841 mole_dr_change = -1; 00842 if( nzone>=4 ) 00843 { 00844 /* first do H2 abundance */ 00845 double Old_H2_abund; 00846 /* >>chng 04 dec 15, do not use special logic when large h2 is turned on */ 00847 int i; 00848 Old_H2_abund = struc.H2_molec[ipMH2g][nzone-3] + struc.H2_molec[ipMH2s][nzone-3]; 00849 /* >>chng 03 jun 16, limit from 0.01 to 0.001, some models fell over H2 front due to 00850 * large zoning, when large H2 was just this caused oscillations in solomon process */ 00851 /* >>chng 03 nov 22, from > 0.001 to > 3e-4, models that start almost in H2 have 00852 * rapid increase in H2 at shallow depths, start sensing this sooner */ 00853 /* >>chng 03 dec 10, from 3e-4 to 1e-4, to get smaller chagned in leiden1.in test */ 00854 /* radius_next keys from change in H2 abundance, d(H2)/H */ 00855 /* >>chng 04 mar 11, start sensing H2 earlier since otherwise step size 00856 * needs to become way too small tooo quickly. limit changed from 1e-4 to 1e-6 */ 00857 /* >>chng 04 jun 29, fromo > 1e-6 to >1e-8, some pdr's had too large a change in H2 */ 00858 00859 if( 2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN] > 1e-8 ) 00860 { 00861 double fac = 0.2; 00862 /* this is percentage change in H2 density - "change in H2 abundance" */ 00863 dH2_abund = 2.*fabs( hmi.H2_total - Old_H2_abund ) / hmi.H2_total; 00864 /* in testing th85ism the dH2_abund did come out zero */ 00865 /* >>chng 03 jun 16, change d(H2) from 0.05 to 0.1, fine resolution of H2/H exposed 00866 * small numerical oscillations in Solomon process */ 00867 /* >>chng 03 nov 22, smallest possible ratio of dr(next)/dr changed from 00868 * 0.2 to 0.05, models that started almost in H2 front need to be able to sense it */ 00869 /*drH2_abund = radius.drad*MAX2(0.05,fac/SDIV(dH2_abund) );*/ 00870 /* >>chng 04 mar 15, with such small possible changes in dr, only 0.05, a thermal front 00871 * can easily cause large changes in H2 and T that are not due to zoning, but to the 00872 * discontinuity. make smallest change larger to prevent hald due to dr */ 00873 /* >>chng 05 aug 23, thermal front allowed dr to become much too small 00874 * chng from 0.02 to .6 */ 00875 dH2_abund = SDIV(dH2_abund); 00876 drH2_abund = radius.drad*MAX2(0.2,fac/dH2_abund ); 00877 } 00878 else 00879 { 00880 drH2_abund = BigRadius; 00881 } 00882 00883 /* check how molecular fractions of all heavy elements are chaning relative 00884 * to total gas phase abundance */ 00885 dr_mole_abund = BigRadius; 00886 /* >>chng 04 jun 02, upper limit had been all species but now limit to real 00887 * molecules since do not want this logic to work with the ions */ 00888 for(i=0; i<mole.num_comole_calc; ++i) 00889 { 00890 realnum abund, 00891 abund_limit; 00892 if( !dense.lgElmtOn[COmole[i]->nelem_hevmol] || COmole[i]->n_nuclei <= 1 ) 00893 continue; 00894 /* >>chng 03 sep 21, add CO logic */ 00895 /* >>chng 04 mar 30, generalize to any molecule at all */ 00896 /* >>chng 04 mar 31 lower limit to abund had been 0.01, change 00897 * to 0.001 to pick up approach to molecular fronts */ 00898 abund = COmole[i]->hevmol/dense.gas_phase[COmole[i]->nelem_hevmol]; 00899 /* is this an ice? need special logic for ice since density increases 00900 * exponentially when grain temp falls below sublimation temperature 00901 * >>chng 05 dec 06 - detect changes for smaller abundances for ices 00902 * due to large changes in ice abundances */ 00903 if( COmole[i]->lgGas_Phase ) 00904 { 00905 abund_limit = 1e-3f; 00906 } 00907 else 00908 { 00909 /* this is an ice - track its abundance at lower values so that 00910 * we resolve the sublimation transition region */ 00911 abund_limit = 1e-5f; 00912 } 00913 00914 if( abund > abund_limit ) 00915 { 00916 double drmole_one, relative_change, relative_change_denom; 00917 /* >>chng 05 dec 08, use smaller abundance for the denominator since just taking 00918 * current abundance will overlook case where current density is vastly large 00919 * than old density */ 00920 if( struc.CO_molec[i][nzone-3]>SMALLFLOAT ) 00921 { 00922 relative_change_denom = MIN2( COmole[i]->hevmol , struc.CO_molec[i][nzone-3] ); 00923 } 00924 else 00925 { 00926 relative_change_denom = COmole[i]->hevmol; 00927 } 00928 /* the relative change in the abundance */ 00929 relative_change = 00930 /*fabs( COmole[i]->hevmol - struc.CO_molec[i][nzone-3] ) / COmole[i]->hevmol;*/ 00931 /* >>chng 05 dec 08, use smaller abundance */ 00932 fabs( COmole[i]->hevmol - struc.CO_molec[i][nzone-3] ) / relative_change_denom; 00933 /* >>chng 03 jun 16, change from 0.05 to 0.1, fine resolution of H2/H exposed 00934 * small numerical oscillations in Solomon process */ 00935 /* >>chng 04 jun 02, from 0.1 back to 0.05, more extensive CO etc network 00936 * caused oscillations in SiO abundance and Si Si+ density. */ 00937 /* >>chng 04 aug 03, from 0.05 to 0.035, leiden pdr model v2 had major 00938 * jump in eden */ 00939 /* >>chng 04 oct 18, from 0.035 back to 0.05, leiden pdr v2 actually due to having 00940 * PAHs in fully molecular limit (??), this caused cool flow pdr grid to trip on 00941 * too small dr */ 00942 relative_change = SDIV(relative_change); 00943 /*>>chng 05 dec 08, relative_change must be less than one - with 00944 * revised logic above can be bigger than one */ 00945 if( relative_change > 1. ) 00946 relative_change = 1./relative_change; 00947 /*drmole_one = radius.drad*MAX2(0.2,0.035/relative_change );*/ 00948 /* >>chng 05 aug 23, thermal front allowed dr to become much too small 00949 * chng from 0.02 to .6 */ 00950 /*drmole_one = radius.drad*MAX2(0.2,0.05/relative_change );*/ 00951 drmole_one = radius.drad*MAX2(0.6,0.05/relative_change ); 00952 /* final dr will be the smallest we encounter */ 00953 if( drmole_one < dr_mole_abund ) 00954 { 00955 /* this is the dr used to set next dr - keep track of which moe was changing */ 00956 dr_mole_abund = drmole_one; 00957 mole_dr_change = i; 00958 dCO_abund = relative_change; 00959 } 00960 } 00961 } 00962 } 00963 00964 /* some consideration due to big H2 molecule */ 00965 drSolomon_BigH2 = H2_DR(); 00966 00967 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/ 00968 /* can't make drmax large deep in model since causes feedback 00969 * oscillations with changes in heating or destruction rates 00970 * >>chng 96 oct 15, change from 5 to 10 */ 00971 if( nzone < 5 ) 00972 { 00973 /* >>chng 96 nov 11, had been 4 * drad up to 11, change to following 00974 * to be similar to c90.01, max of 1.3 between 5 and 11 00975 * >>chng 04 oct 29 geometry.DirectionalCosin was ioncorrect applied to this factor */ 00976 drmax = 4.*radius.drad; 00977 } 00978 else 00979 { 00980 drmax = 1.3*radius.drad; 00981 } 00982 00983 /* >>chng 05 apr 05, do not sense temp oscillation, so that we can move past this 00984 * point if it occurs */ 00985 # if 0 00986 /* look for oscillations in electron density or tempeature - freeze dr if these occur */ 00987 dr_ne_oscil = BigRadius; 00988 dr_te_oscil = BigRadius; 00989 if( nzone >= 11 ) 00990 { 00991 /* >>chng 96 oct 15, do not let zones increase if oscillations present */ 00992 /* >>chng 96 oct 31, error to declare oscillation propto toler, the 00993 *heating cooling tolerance */ 00994 realnum errorHC = POW2(conv.HeatCoolRelErrorAllowed); 00995 realnum errorNe = (realnum)POW2(conv.EdenErrorAllowed ); 00996 limit = nzone -2; 00997 ASSERT( limit < struc.nzlim ); 00998 for( k=nzone - 10; k < limit; k++ ) 00999 { 01000 /* check that square of change both chng sign and is 01001 * greater than square of heat-tool error */ 01002 if( (struc.testr[k-1] - struc.testr[k])/struc.testr[k]* 01003 (struc.testr[k] - struc.testr[k+1])/struc.testr[k] < 01004 -errorHC ) 01005 { 01006 dr_te_oscil = radius.drad; 01007 } 01008 /* small residiual is to allow 0.01 rel error */ 01009 if( (struc.ednstr[k-1] - struc.ednstr[k])/struc.ednstr[k]* 01010 (struc.ednstr[k] - struc.ednstr[k+1])/struc.ednstr[k] < 01011 -errorNe ) 01012 { 01013 /* >>chng 96 oct 15, do not let zones increase if oscillations present */ 01014 /* radius_next keys from electron density oscillation*/ 01015 dr_ne_oscil = radius.drad; 01016 } 01017 } 01018 } 01019 /* >>chng 05 apr 05, do not sense temp oscillation, so that we can move past this 01020 * point if it occurs */ 01021 dr_ne_oscil = BigRadius; 01022 dr_te_oscil = BigRadius; 01023 # endif 01024 01025 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/ 01026 /* check on several convergence criteria */ 01027 01028 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/ 01029 if( !conv.lgConvTemp ) 01030 { 01031 drFail = radius.drad/2.; 01032 } 01033 else 01034 { 01035 drFail = BigRadius; 01036 } 01037 01038 /* time dependent recombination - conditions become very homogeneous and 01039 * crash into I fronts - use last iteration's radius as guess of current case */ 01040 recom_dr_last_iter = BigRadius; 01041 if( dynamics.lgStatic && dynamics.lgRecom ) 01042 { 01043 /* first time through nzone == 1 */ 01044 static long int nzone_recom = -1; 01045 if( nzone<=1 ) 01046 { 01047 /* initialization case */ 01048 nzone_recom = 0; 01049 } 01050 else if( radius.depth < struc.depth_last[struc.nzone_last-1] && 01051 nzone_recom < struc.nzone_last ) 01052 { 01053 ASSERT( nzone_recom>=0 && nzone_recom<struc.nzone_last ); 01054 /* case where we are within previous computed structure 01055 * first possibly increase nzone_recom so that it points deeper 01056 * than this zone */ 01057 while( struc.depth_last[nzone_recom] < radius.depth && 01058 nzone_recom < struc.nzone_last-1 ) 01059 ++nzone_recom; 01060 recom_dr_last_iter = struc.drad_last[nzone_recom]*3.; 01061 } 01062 else 01063 { 01064 /* case where we have overrun the previous iteration's saved structure */ 01065 recom_dr_last_iter = BigRadius; 01066 } 01067 } 01068 01069 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/ 01070 /* change in electron density - radius_next keys from change in elec den, 01071 * remember old electron density during first call */ 01072 /* this is low ionization solution */ 01073 if( nzone > 2 ) 01074 { 01075 /* next is-2 since nzone is on physics not c scale, and we want previous zone */ 01076 Efrac_old = struc.ednstr[nzone-3]/struc.gas_phase[ipHYDROGEN][nzone-3]; 01077 Efrac_new = dense.eden/dense.gas_phase[ipHYDROGEN]; 01078 dEfrac = fabs(Efrac_old-Efrac_new) / Efrac_new; 01079 01080 if( dEfrac > SMALLFLOAT ) 01081 { 01082 double fac = 0.04; 01083 /* >>chng 03 dec 25, use smaller rel change in elec frac when most elec in ipMole or grains */ 01084 /* >>chng 04 sep 14, change to from from metals but comment out */ 01085 /* >>chng 04 sep 17, change to from from metals - uncomment */ 01086 if( dense.eden_from_metals > 0.5 ) 01087 { 01088 /* >>chng 04 sep 18, change 0.02 from 0.01 */ 01089 /* >>chng 04 sep 18, change 0.02 from 0.04 */ 01090 fac = 0.04; 01091 } 01092 /* >>chng 04 feb 23, add test on hydrogen being predomintly 01093 * recombined due to three-body recom, which is very sensitive 01094 * to the electron density - but only do this in partially ionized medium */ 01095 else if( iso.RecomCollisFrac[ipH_LIKE][ipHYDROGEN] > 0.8 && 01096 dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN]>0.1 && 01097 dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN]<0.8 ) 01098 01099 { 01100 fac = 0.02; 01101 } 01102 /* >>chng 04 sep 17, change to 0.1 from 0.2 */ 01103 drEfrac = radius.drad*MAX2(0.1,fac/dEfrac); 01104 } 01105 else 01106 { 01107 drEfrac = BigRadius; 01108 } 01109 } 01110 else 01111 { 01112 dEfrac = 0.; 01113 drEfrac = BigRadius; 01114 Efrac_old = 0.; 01115 Efrac_new = 0.; 01116 } 01117 01118 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/ 01119 /* do not let thickness get too large 01120 * 1 '' radius_next keys from relative depth'',e11.3)') */ 01121 if( nzone > 20 ) 01122 { 01123 /*drDepth = radius.depth/20.;*/ 01124 /* >>chng 02 nov 05, change from 1/20 to 1/10 wasted zones early on */ 01125 drDepth = radius.depth/10.; 01126 } 01127 else 01128 { 01129 drDepth = BigRadius; 01130 } 01131 01132 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/ 01133 /* case where stopping thickness or edge specified, need to approach slowly */ 01134 thickness_total = BigRadius; 01135 DepthToGo = BigRadius; 01136 if( StopCalc.HColStop < 5e29 ) 01137 { 01138 double coleff = SDIV( dense.gas_phase[ipHYDROGEN]*geometry.FillFac); 01139 DepthToGo = MIN2(DepthToGo , 01140 (StopCalc.HColStop-colden.colden[ipCOL_HTOT]) / coleff ); 01141 /* >>chng 03 oct 28, forgot to div col den above by eff density */ 01142 thickness_total = MIN2(thickness_total , StopCalc.HColStop / coleff ); 01143 } 01144 01145 if( StopCalc.colpls < 5e29 ) 01146 { 01147 double coleff = (double)SDIV( (dense.xIonDense[ipHYDROGEN][1])*geometry.FillFac); 01148 DepthToGo = MIN2(DepthToGo , 01149 (StopCalc.colpls-colden.colden[ipCOL_Hp]) / coleff ); 01150 thickness_total = MIN2(thickness_total , StopCalc.colpls / coleff ); 01151 } 01152 01153 if( StopCalc.col_h2 < 5e29 ) 01154 { 01155 /* >>chng 03 apr 15, add this molecular hydrogen */ 01156 double coleff = (double)SDIV( hmi.H2_total*geometry.FillFac); 01157 DepthToGo = MIN2(DepthToGo , 01158 (StopCalc.col_h2-colden.colden[ipCOL_H2g]-colden.colden[ipCOL_H2s]) / coleff ); 01159 thickness_total = MIN2(thickness_total , StopCalc.col_h2 / coleff ); 01160 } 01161 01162 if( StopCalc.col_h2_nut < 5e29 ) 01163 { 01164 /* >>chng 03 apr 15, add this molecular hydrogen */ 01165 double coleff = (double)SDIV( (2*hmi.H2_total+dense.xIonDense[ipHYDROGEN][1])*geometry.FillFac); 01166 DepthToGo = MIN2(DepthToGo , 01167 (StopCalc.col_h2_nut-(2*(colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s])+dense.xIonDense[ipHYDROGEN][1])) / coleff ); 01168 thickness_total = MIN2(thickness_total , StopCalc.col_h2_nut / coleff ); 01169 } 01170 01171 if( StopCalc.col_H0_ov_Tspin < 5e29 ) 01172 { 01173 /* >>chng 05 jan 09, add n(H0)/Tspin */ 01174 double coleff = (double)SDIV( dense.xIonDense[ipHYDROGEN][0] / hyperfine.Tspin21cm*geometry.FillFac ); 01175 DepthToGo = MIN2(DepthToGo , 01176 (StopCalc.col_H0_ov_Tspin - colden.H0_ov_Tspin) / coleff ); 01177 thickness_total = MIN2(thickness_total , StopCalc.col_H0_ov_Tspin / coleff ); 01178 } 01179 01180 if( StopCalc.col_monoxco < 5e29 ) 01181 { 01182 /* >>chng 03 apr 15, add this, CO */ 01183 double coleff = (double)SDIV( (findspecies("CO")->hevmol)*geometry.FillFac); 01184 DepthToGo = MIN2(DepthToGo , 01185 (StopCalc.col_monoxco-findspecies("CO")->hevcol) / coleff ); 01186 thickness_total = MIN2(thickness_total , StopCalc.col_monoxco / coleff ); 01187 } 01188 01189 if( StopCalc.colnut < 5e29 ) 01190 { 01191 double coleff = (double)SDIV( (dense.xIonDense[ipHYDROGEN][0])*geometry.FillFac); 01192 DepthToGo = MIN2(DepthToGo , 01193 (StopCalc.colnut - colden.colden[ipCOL_H0]) / coleff ); 01194 thickness_total = MIN2(thickness_total , StopCalc.colnut / coleff ); 01195 } 01196 01197 /* this is case where outer radius is set */ 01198 if( radius.router[iteration-1] < 5e29 ) 01199 { 01200 thickness_total = MIN2(thickness_total , radius.router[iteration-1] ); 01201 DepthToGo = MIN2(DepthToGo , 01202 radius.router[iteration-1] - radius.depth ); 01203 } 01204 01205 /* this is case where stopping optical depth was specified */ 01206 if( StopCalc.iptnu != rfield.nupper ) 01207 { 01208 /* end optical depth has been specified */ 01209 double dt = SDIV(opac.opacity_abs[StopCalc.iptnu-1]*geometry.FillFac); 01210 DepthToGo = MIN2(DepthToGo , 01211 (StopCalc.tauend-opac.TauAbsGeo[0][StopCalc.iptnu-1])/dt); 01212 } 01213 /* stop AV - usually this is dust, but we consider all opacity sources, 01214 * so always include this */ 01215 /* compute some average grain properties */ 01216 AV_to_go = BigRadius; 01217 if( rfield.opac_mag_V_extended > SMALLFLOAT && rfield.opac_mag_V_point>SMALLFLOAT ) 01218 { 01219 /* by default stop av is very large, and opacity can be very small, so ratio 01220 * goes to inf - work with logs to see how big the number is */ 01221 double ave = log10(StopCalc.AV_extended - rfield.extin_mag_V_extended) - 01222 log10(rfield.opac_mag_V_extended); 01223 double avp = log10(StopCalc.AV_point - rfield.extin_mag_V_point) - 01224 log10(rfield.opac_mag_V_point); 01225 AV_to_go = MIN2( ave , avp ); 01226 if( AV_to_go < 37. ) 01227 { 01228 AV_to_go = pow(10., AV_to_go ); 01229 /* this is to make sure that we go slightly deeper than AV so that 01230 * we ttrigger this stop */ 01231 AV_to_go *= 1.0001; 01232 } 01233 else 01234 AV_to_go = BigRadius; 01235 /*fprintf(ioQQQ,"DEBUG next dr %.3e %.3e %.3e\n", AV_to_go , ave, avp );*/ 01236 } 01237 01238 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/ 01239 /* set dr if one of above tests have triggered */ 01240 if( DepthToGo <= 0. ) 01241 { 01242 TotalInsanity(); 01243 } 01244 else if( DepthToGo < BigRadius ) 01245 { 01246 /* want to approach the outer edge slowly, 01247 * the need for this logic is most evident in brems.in - 01248 * HI fraction varies across coronal model */ 01249 drThickness = MIN2( thickness_total/10. , DepthToGo ); 01250 } 01251 else 01252 { 01253 drThickness = BigRadius; 01254 } 01255 01256 /*fprintf(ioQQQ, 01257 "DEBUG depth2go z%li drThickness2 = %e %e %e\n", 01258 nzone , drThickness , thickness_total , DepthToGo );*/ 01259 01260 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/ 01261 /* spherical models, do not want delta R/R big */ 01262 drSphere = radius.Radius*0.04; 01263 01264 /* optical depth to electron scattering */ 01265 /* >>chng 04 jun 16, add filling factor, was missing */ 01266 dRTaue = radius.drChange/(dense.eden*6.65e-25*geometry.FillFac); 01267 /* >>chng 02 oct 06, increase dr when constant temperature, 01268 * to prevent some ct models from taking too many cells */ 01269 if( thermal.lgTemperatureConstant ) 01270 dRTaue *= 3.; 01271 01272 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/ 01273 if( dense.flong != 0. ) 01274 { 01275 drFluc = 0.628/2./dense.flong; 01276 /* >>chng 04 sep 18, caused cautions that ionization jumps occurred. 01277 * set to have the value */ 01278 drFluc /= 2.; 01279 } 01280 else 01281 { 01282 drFluc = BigRadius; 01283 } 01284 01285 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/ 01286 /* if density fluctuations in place then override change in heat 01287 * for dr set */ 01288 if( strcmp(dense.chDenseLaw,"SINE") == 0 && dense.lgDenFlucOn ) 01289 { 01290 drdHdStep = BigRadius; 01291 } 01292 01293 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/ 01294 /*active dr sets */ 01295 /* we are deep into model, use logic since we have several zones 01296 * of old data */ 01297 radius.drNext = MIN2( MIN4( drmax, drInter, drLineHeat, winddr ), 01298 MIN4( drFluc, GlobDr, DrGrainHeat, dr_change_heavy ) ); 01299 radius.drNext = MIN2( MIN4( radius.drNext, SpecDr, drFail, WindAccelDR ), 01300 MIN3( drSphere, radius.sdrmax, dRTaue ) ); 01301 radius.drNext = MIN3( MIN3( radius.drNext, drDest, drdTdStep ), 01302 MIN3( drdHdStep, drConPres, drTab ), 01303 MIN3( drSolomon_BigH2, drLeiden_hack, recom_dr_last_iter ) ); 01304 radius.drNext = MIN3( MIN4( radius.drNext, drHion, drDepth, dr_mole_abund ), 01305 MIN4( AV_to_go, drEfrac, drHMase, drThickness ), 01306 MIN3( drHe1ovHe2, drH2_heat_cool, drH2_abund ) ); 01307 01308 /*fprintf(ioQQQ, 01309 "DEBUG depth2go drNext = %e \n", radius.drNext );*/ 01310 01311 /* keep dr constant in first two zones, if it wants to increase, 01312 * but always allow it to decrease. 01313 * to guard against large increases in efrac for balmer cont photo dominated models, 01314 */ 01315 if( nzone <= 1 && radius.drNext > OldDR ) 01316 { 01317 radius.drNext = OldDR; 01318 } 01319 01320 /* option to force min drad */ 01321 if( radius.drNext < radius.sdrmin ) 01322 { 01323 radius.drNext = radius.sdrmin; 01324 } 01325 01328 /* a pressure failure has occurred - keep zone the same time, hoping to pass through 01329 * troubled region */ 01330 if( !conv.lgConvPres || !conv.lgConvTemp ) 01331 { 01332 radius.drNext = radius.drad; 01333 } 01334 01335 /* >>chng 05 aug 05, in case of thermal front, where temp is falling quickly and 01336 * conditions change very fast, the zone thickness does not really affect the change 01337 * in conditions and can cause zoning to become very very thin, which causes 01338 * an abort. this occurs between 200 and 1000K. if we are doing temp soln, 01339 * temp is between these values, and temp is changing rapidly, do not make sone 01340 * thickness much smaller */ 01341 drThermalFront = BigRadius; 01342 /* do not use thermal front logic in case of recombination time dep cloud */ 01343 if( nzone >=5 && !dynamics.lgStatic ) 01344 { 01345 /* >>chng 05 aug 23, upper bound of thermal from from 1000K to 4000K */ 01346 /*if( phycon.te > 200. && phycon.te < 1000. && */ 01347 if( phycon.te > 200. && phycon.te < 3000. && 01348 /* >>chng 05 aug 23, from > 10% in zone to to two zones > 5%, 01349 * to fix leiden v3 with large H2 */ 01350 (struc.testr[nzone-3] - phycon.te)/phycon.te > 0.02 && 01351 (struc.testr[nzone-4] - phycon.te)/phycon.te > 0.02 && 01352 (struc.testr[nzone-5] - phycon.te)/phycon.te > 0.02 ) 01353 { 01354 /* the 0.91 is to make dr unique, so that print statement that 01355 * follows will identify this reason */ 01356 drThermalFront = radius.drad * 0.91; 01357 radius.drNext = drThermalFront; 01358 } 01359 } 01360 01361 /* dr = zero is a logical mistake */ 01362 if( radius.drNext <= 0. ) 01363 { 01364 fprintf( ioQQQ, " radius_next finds insane drNext:%10.2e\n", 01365 radius.drNext ); 01366 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", 01367 drmax, drInter, drLineHeat, winddr, drFluc, GlobDr, SpecDr, 01368 drFail, drSphere, radius.sdrmax, dRTaue, 01369 OldH2Abund, drDepth ); 01370 cdEXIT(EXIT_FAILURE); 01371 } 01372 01373 /* set flag if dr set by maser */ 01374 if( fp_equal( radius.drNext, drHMase ) ) 01375 { 01376 rt.lgMaserSetDR = true; 01377 } 01378 01379 /* all this is to only punch on last iteration 01380 * the punch dr command is not really a punch command, making this necessary 01381 * lgDRon is set true if "punch dr" entered */ 01382 if( punch.lgDROn ) 01383 { 01384 if( punch.lgDRPLst ) 01385 { 01386 /* lgDRPLst was set true if "punch" had "last" on it */ 01387 if( iterations.lgLastIt ) 01388 { 01389 lgDoPun = true; 01390 } 01391 else 01392 { 01393 lgDoPun = false; 01394 } 01395 } 01396 else 01397 { 01398 lgDoPun = true; 01399 } 01400 } 01401 else 01402 { 01403 lgDoPun = false; 01404 } 01405 if( (trace.lgTrace && trace.lgDrBug) || lgDoPun ) 01406 { 01407 if( !conv.lgConvTemp && nzone > 0 ) 01408 { 01409 fprintf( punch.ipDRout, "#>>>> A temperature failure occured.\n" ); 01410 } 01411 if( !conv.lgConvPres && nzone > 0 ) 01412 { 01413 fprintf( punch.ipDRout, "#>>>> A pressure failure occured.\n" ); 01414 } 01415 01416 /* this is common part of each line, the zone count, depth, chosen dr, and depth2go */ 01417 fprintf( punch.ipDRout , "%ld\t%.5e\t%.3e\t%.3e\t", nzone, radius.depth, radius.drNext, radius.Depth2Go ); 01418 01419 /*=======begin active dr sets */ 01420 if( fp_equal( radius.drNext, drLineHeat ) ) 01421 { 01422 if( level == 1 ) 01423 { 01424 strcpy( chLbl, chLineLbl(&TauLines[ipStrong]) ); 01425 fprintf( punch.ipDRout, "level 1 line heating,%10.10s TauIn%10.2e%10.2e%10.2e\n", 01426 chLbl, TauInwd, TauLines[ipStrong].Emis->pump, 01427 TauLines[ipStrong].Emis->Pesc ); 01428 } 01429 else if( level == 2 ) 01430 { 01431 strcpy( chLbl, chLineLbl(&TauLine2[ipStrong])); 01432 fprintf( punch.ipDRout, "level 2 line heating,%10.10s TauIn%10.2e%10.2e%10.2e\n", 01433 chLbl, TauInwd, TauLine2[ipStrong].Emis->pump, 01434 TauLine2[ipStrong].Emis->Pesc ); 01435 } 01436 else 01437 { 01438 fprintf( ioQQQ, " insanity pr line heat\n" ); 01439 cdEXIT(EXIT_FAILURE); 01440 } 01441 } 01442 01443 else if( fp_equal( radius.drNext, drDepth ) ) 01444 { 01445 fprintf( punch.ipDRout, "relative depth\n"); 01446 } 01447 01448 else if( fp_equal( radius.drNext, drThermalFront ) ) 01449 { 01450 fprintf( punch.ipDRout, "thermal front logic\n"); 01451 } 01452 01453 else if( fp_equal( radius.drNext, dr_change_heavy ) ) 01454 { 01455 ASSERT( ichange_heavy_nelem < LIMELM+3 ); 01456 fprintf( punch.ipDRout, 01457 "change in ionization, element %s ion (C scale) %li rel change %.2e ion frac %.2e\n", 01458 elementnames.chElementName[ichange_heavy_nelem], 01459 ichange_heavy_ion , 01460 change_heavy_frac_big , 01461 frac_change_heavy_frac_big); 01462 } 01463 01464 else if( fp_equal( radius.drNext, drThickness ) ) 01465 { 01466 fprintf( punch.ipDRout, "depth to go\n"); 01467 } 01468 01469 else if( fp_equal( radius.drNext, AV_to_go ) ) 01470 { 01471 fprintf( punch.ipDRout, "A_V to go\n"); 01472 } 01473 01474 else if( fp_equal( radius.drNext, drTab ) ) 01475 { 01476 fprintf( punch.ipDRout, "spec den law, new old den%10.2e%10.2e\n", 01477 hdnew, dense.gas_phase[ipHYDROGEN] ); 01478 } 01479 01480 else if( fp_equal( radius.drNext, drHMase ) ) 01481 { 01482 fprintf( punch.ipDRout, "H maser dTauMase=%10.2e %li %li %li %li\n", 01483 rt.dTauMase, 01484 rt.mas_species, 01485 rt.mas_ion, 01486 rt.mas_hi, 01487 rt.mas_lo ); 01488 } 01489 01490 else if( fp_equal( radius.drNext, drHe1ovHe2 ) ) 01491 { 01492 /* radius_next keys from change in He2/He1 ionization, old=%11.3e sv=%11.3e FrLya*/ 01493 fprintf( punch.ipDRout, "change in He0/He+ ionization, ratio %.2e\n", 01494 Old_He_atom_ov_ion ); 01495 } 01496 01497 else if( fp_equal( radius.drNext, drH2_heat_cool ) ) 01498 { 01499 /* radius_next keys from change in H2 heating, old=%11.3e sv=%11.3e FrLya*/ 01500 fprintf( punch.ipDRout, "change in H2 heating/cooling, d(c,h)/H %.2e\n", 01501 dH2_heat_cool ); 01502 } 01503 01504 else if( fp_equal( radius.drNext, drH2_abund ) ) 01505 { 01506 /* radius_next keys from change in H2 abundance, old=%11.3e sv=%11.3e FrLya*/ 01507 fprintf( punch.ipDRout, "change in H2 abundance, d(H2)/H %.2e\n", 01508 dH2_abund ); 01509 } 01510 01511 else if( fp_equal( radius.drNext, dr_mole_abund ) ) 01512 { 01513 /* radius_next keys from change in CO mole abundance */ 01514 fprintf( punch.ipDRout, "change in heav ele mole abundance, d(mole)/elem %.2e mole=%i=%s\n", 01515 dCO_abund , mole_dr_change , COmole[mole_dr_change]->label); 01516 } 01517 01518 else if( fp_equal( radius.drNext, drSolomon_BigH2 ) ) 01519 { 01520 /* radius_next keys from change in H2 abundance, old=%11.3e sv=%11.3e FrLya*/ 01521 fprintf( punch.ipDRout, "change in big H2 Solomon rate line opt depth\n"); 01522 } 01523 01524 else if( fp_equal( radius.drNext, recom_dr_last_iter ) ) 01525 { 01526 /* radius_next keys from previous iteration recom logic */ 01527 fprintf( punch.ipDRout, "previous iteration recom logic\n"); 01528 } 01529 01530 else if( fp_equal( radius.drNext, drHion ) ) 01531 { 01532 fprintf( punch.ipDRout, "change in H ionization fm to%11.3e%11.3e\n", 01533 SaveOHIIoHI, OHIIoHI ); 01534 } 01535 01536 else if( fp_equal( radius.drNext, drEfrac ) ) 01537 { 01538 fprintf( punch.ipDRout, 01539 "change in elec den, rel chng:%11.3e, cur %11.3e old=%11.3e\n", 01540 dEfrac , Efrac_old , Efrac_new ); 01541 } 01542 01543 else if( fp_equal( radius.drNext, drdHdStep ) ) 01544 { 01545 fprintf( punch.ipDRout, 01546 "change in heating; current %10.3e delta=%10.3e\n", 01547 thermal.htot, dHdStep ); 01548 } 01549 01550 else if( fp_equal( radius.drNext, drLeiden_hack ) ) 01551 { 01552 fprintf( punch.ipDRout, 01553 "Leiden hack\n" ); 01554 } 01555 01556 else if( fp_equal( radius.drNext, drConPres ) ) 01557 { 01558 fprintf( punch.ipDRout, "change in con accel\n" ); 01559 } 01560 01561 else if( fp_equal( radius.drNext, drdTdStep ) ) 01562 { 01563 fprintf( punch.ipDRout, 01564 "change in temperature; current= %10.3e, dT/T= %.3f\n", 01565 phycon.te, dTdStep ); 01566 } 01567 01568 else if( fp_equal( radius.drNext, radius.sdrmin ) ) 01569 { 01570 fprintf( punch.ipDRout, "sdrmin\n" ); 01571 } 01572 01573 else if( fp_equal( radius.drNext, radius.sdrmax ) ) 01574 { 01575 fprintf( punch.ipDRout, "sdrmax\n" ); 01576 } 01577 01578 else if( fp_equal( radius.drNext, drSphere ) ) 01579 { 01580 fprintf( punch.ipDRout, "sphericity\n" ); 01581 } 01582 01583 else if( fp_equal( radius.drNext, dRTaue ) ) 01584 { 01585 fprintf( punch.ipDRout, 01586 "optical depth to electron scattering\n" ); 01587 } 01588 01589 else if( fp_equal( radius.drNext, drFail ) ) 01590 { 01591 fprintf( punch.ipDRout, 01592 "temperature failure.\n" ); 01593 } 01594 01595 else if( fp_equal( radius.drNext, drmax ) ) 01596 { 01597 fprintf( punch.ipDRout, 01598 "DRMAX; nu opc dr=%10.2e%10.2e%10.2e\n", 01599 freqm, opacm, radius.drChange/ 01600 SDIV(opacm) ); 01601 } 01602 01603 else if( fp_equal( radius.drNext, drInter ) ) 01604 { 01605 fprintf( punch.ipDRout, 01606 "cont inter nu=%10.2e opac=%10.2e\n", 01607 freqm, opacm ); 01608 } 01609 01610 else if( fp_equal( radius.drNext, DrGrainHeat ) ) 01611 { 01612 01613 fprintf( punch.ipDRout, 01614 "grain heating nu=%10.2e opac=%10.2e\n", 01615 grfreqm, gropacm ); 01616 } 01617 01618 else if( fp_equal( radius.drNext, winddr ) ) 01619 { 01620 fprintf( punch.ipDRout, 01621 "Wind, dVelRelative=%.3e windv=%.3e\n", 01622 dVelRelative , 01623 wind.windv); 01624 } 01625 01626 else if( fp_equal( radius.drNext, drFluc ) ) 01627 { 01628 fprintf( punch.ipDRout, 01629 "density fluctuations\n" ); 01630 } 01631 01632 else if( fp_equal( radius.drNext, GlobDr ) ) 01633 { 01634 fprintf( punch.ipDRout, 01635 "GLOB law new dr=%10.2e HDEN=%10.2e\n", 01636 GlobDr, 01637 dense.gas_phase[ipHYDROGEN] ); 01638 } 01639 01640 else if( fp_equal( radius.drNext, SpecDr ) ) 01641 { 01642 fprintf( punch.ipDRout, 01643 "special law new dr=%10.2e HDEN=%10.2e\n", 01644 SpecDr, 01645 dense.gas_phase[ipHYDROGEN] ); 01646 } 01647 01648 else if( fp_equal( radius.drNext, OldDR ) ) 01649 { 01650 fprintf( punch.ipDRout, "old DR.\n" ); 01651 } 01652 01653 else 01654 { 01655 fprintf( punch.ipDRout, 01656 " %4ld radius_next keys from insanity %10.2e\n", 01657 nzone, radius.drNext ); 01658 01659 fprintf( ioQQQ, 01660 " %4ld radius_next keys from insanity %10.2e\n", 01661 nzone, radius.drNext ); 01662 cdEXIT(EXIT_FAILURE); 01663 } 01664 01665 /*======end active dr sets */ 01666 } 01667 01668 /* this is general code that prevents zone thickness drNext from 01669 * becoming too thin, something that can happen for various bad reasons 01670 * HOWEVER we do not want to do this test for certain density laws, 01671 * for which very small zone thicknesses are unavoidable 01672 * the special cases are: 01673 * special density law, 01674 * globule density law, 01675 * carbon +-0 i front 01676 * the fluctuations command 01677 * drMinimum was set in radius_first to either sdrmin (set drmin) or 01678 * some fraction of the initial radius - it is always set 01679 * to something non-trivial. 01680 * sdrmin is set with the "set dr" command - is SMALLFLOAT by default */ 01681 if( ((strcmp(dense.chDenseLaw,"DLW1") != 0 && 01682 strcmp(dense.chDenseLaw ,"GLOB") != 0) )&& 01683 /* >>chng 04 feb 19, do not use this test - errors can still happen 01684 * when all C is atomic! */ 01685 /*dense.xIonDense[ipCARBON][0]/dense.gas_phase[ipCARBON] < 0.05) && */ 01686 (dense.flong == 0.) && 01687 /* >>chng 01 aug 11, add check against stopping on depth to go */ 01688 radius.drNext != DepthToGo ) 01689 { 01690 /* don't let dr get smaller than drMinimum, if this resets drNext 01691 * then code stops with warning that zones got too thin 01692 * this can happen due to numerical oscillations, which the code 01693 * tries to damp out by making the zones thinner. 01694 * scale by radius.Radius/radius.rinner to stop very spherical 01695 * simulations from false trigger on dr too small */ 01696 if( radius.drNext* radius.Radius/radius.rinner * 01697 /* drMinimum is drad * hden, to make it proportional to optical depth. 01698 * This avoids false trigger across thermal fronts. */ 01699 dense.gas_phase[ipHYDROGEN] < 01700 radius.drMinimum ) 01701 { 01702 radius.drNext = radius.drMinimum/ dense.gas_phase[ipHYDROGEN]; 01703 /* leaving this at true will cause the model to stop with a warning 01704 * that the zone thickness is too small */ 01705 radius.lgDrMinUsed = true; 01706 /* set abort handler */ 01707 lgAbort = true; 01708 /* must decrement nzone, since we will not complete this zone, and will not have 01709 * valid structure data for it */ 01710 --nzone; 01711 fprintf( ioQQQ, 01712 "\n DISASTER PROBLEM radius_next finds dr too small and aborts. " 01713 "This is zone %ld iteration %ld. The thickness was %.2e\n", 01714 nzone, 01715 iteration, 01716 radius.drNext); 01717 fprintf( ioQQQ, 01718 " If this simulation seems reasonable you can override this limit " 01719 "with the command SET DRMIN %.2e\n\n", 01720 radius.drNext*2); 01721 /* set abort flag */ 01722 lgAbort = true; 01723 return(1); 01724 } 01725 } 01726 01727 /* factor to allow for slop in floating numbers */ 01728 # define Z 1.0001 01729 01730 /* following is to make thickness of model exact 01731 * n.b., on last zone, drNext can be NEGATIVE!! 01732 * DEPTH was incremented at start of zone calc in newrad, 01733 * has been outer edge of zone all throughout */ 01734 radius.drNext = MIN2(radius.drNext,(radius.router[iteration-1]- 01735 radius.depth)*Z); 01736 01737 /* this means outer limit exceeded */ 01738 if( radius.drNext < 0. ) 01739 { 01740 radius.lgDrNeg = true; 01741 } 01742 else 01743 { 01744 radius.lgDrNeg = false; 01745 } 01746 01747 ASSERT( radius.drNext > 0. ); 01748 01749 if( trace.lgTrace ) 01750 { 01751 fprintf( ioQQQ, " radius_next chooses next drad drNext=%.4e; this drad was%12.4e\n", 01752 radius.drNext, radius.drad ); 01753 } 01754 01755 return( 0 ); 01756 } 01757 01758 /*ContRate called by radius_next to find energy of maximum continuum-gas interaction */ 01759 STATIC void ContRate(double *freqm, 01760 double *opacm) 01761 { 01762 long int i, 01763 ipHi, 01764 iplow, 01765 limit; 01766 double FreqH, 01767 Freq_nonIonizing, 01768 Opac_Hion, 01769 Opac_nonIonizing, 01770 Rate_max_Hion, 01771 Rate_max_nonIonizing; 01772 01773 DEBUG_ENTRY( "ContRate()" ); 01774 01775 /* 01776 * find maximum continuum interaction rate, 01777 * these should be reset in following logic without exception, 01778 * if they are still zero at the end we have a logical error 01779 */ 01780 Rate_max_nonIonizing = 0.; 01781 Freq_nonIonizing = 0.; 01782 Opac_nonIonizing = 0.; 01783 01784 /* this must be reset to val >= 0 */ 01785 *opacm = -1.; 01786 *freqm = -1.; 01787 01788 /* do up to carbon photo edge if carbon is turned on */ 01789 /* >>>chng 00 apr 07, add test for whether element is turned on */ 01790 if( dense.lgElmtOn[ipCARBON] ) 01791 { 01792 /* carbon is turned on, use carbon 1 edge */ 01793 ipHi = Heavy.ipHeavy[ipCARBON][0] - 1; 01794 } 01795 else 01796 { 01797 /* carbon turned off, use hydrogen Balmer continuum */ 01798 ipHi = iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2s]-1; 01799 } 01800 01801 /* start at plasma frequency */ 01802 for( i=rfield.ipPlasma; i < ipHi; i++ ) 01803 { 01804 /* this does not have grain opacity since grains totally passive 01805 * at energies smaller than CI edge */ 01806 if( rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*(opac.opacity_abs[i] - 01807 gv.dstab[i]*dense.gas_phase[ipHYDROGEN]) > Rate_max_nonIonizing ) 01808 { 01809 Rate_max_nonIonizing = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]* 01810 (opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN]); 01811 Freq_nonIonizing = rfield.anu[i]; 01812 Opac_nonIonizing = opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN]; 01813 } 01814 } 01815 01816 /* not every continuum extends beyond C edge-this whole loop can add to zero 01817 * use total opacity here 01818 * test is to put in fir continuum if free free heating is important */ 01819 if( CoolHeavy.brems_heat_total/thermal.htot > 0.05 ) 01820 { 01821 /* this is index for energy where cloud free free optical depth is unity, 01822 * is zero if no freq are opt thin */ 01823 iplow = MAX2(1 , rfield.ipEnergyBremsThin ); 01824 } 01825 else 01826 { 01827 /* >>>chng 00 apr 07, from Heavy.ipHeavy[0][5] to ipHi defined above, since 01828 * would crash if element not defined */ 01829 iplow = ipHi; 01830 } 01831 /* do not go below plasma frequency */ 01832 iplow = MAX2( rfield.ipPlasma , iplow ); 01833 01834 /* this energy range from carbon edge to hydrogen edge */ 01835 limit = MIN2(rfield.nflux,iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1); 01836 for( i=iplow; i < limit; i++ ) 01837 { 01838 if( rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*(opac.opacity_abs[i] - 01839 gv.dstab[i]*dense.gas_phase[ipHYDROGEN]) > Rate_max_nonIonizing ) 01840 { 01841 Rate_max_nonIonizing = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]* 01842 (opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN]); 01843 Freq_nonIonizing = rfield.anu[i]; 01844 Opac_nonIonizing = opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN]; 01845 } 01846 } 01847 01848 /* variables to check continuum interactions over Lyman continuum */ 01849 Rate_max_Hion = 0.; 01850 Opac_Hion = 0.; 01851 FreqH = 0.; 01852 01853 /* not every continuum extends beyond 1 Ryd-this whole loop can add to zero */ 01854 for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1; i < rfield.nflux; i++ ) 01855 { 01856 if( rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*(opac.opacity_abs[i] - 01857 gv.dstab[i]*dense.gas_phase[ipHYDROGEN]) > Rate_max_Hion ) 01858 { 01859 /* Rate_max_Hion = anu(i)*flux(i)/widflx(i)*opac(i) */ 01860 Rate_max_Hion = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]* 01861 (opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN]); 01862 FreqH = rfield.anu[i]; 01863 Opac_Hion = opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN]; 01864 } 01865 } 01866 01867 /* use Lyman continuum if its opacity is larger than non-h ion */ 01868 if( Rate_max_nonIonizing < 1e-30 && Opac_Hion > SMALLFLOAT ) 01869 { 01870 /* this happens for laser source - use Lyman continuum */ 01871 *opacm = Opac_Hion; 01872 *freqm = FreqH; 01873 } 01874 /* >>chng 05 aug 03, add last test on Opac_Hion for case where we go very 01875 * deep and very little radiation is left */ 01876 else if( Opac_Hion > Opac_nonIonizing && Rate_max_Hion/Rate_max_nonIonizing > 1e-10 && Opac_Hion > SMALLFLOAT ) 01877 { 01878 /* use Lyman continuum */ 01879 *opacm = Opac_Hion; 01880 *freqm = FreqH; 01881 } 01882 else 01883 { 01884 /* not much rate in the Lyman continuum, stick with low energy */ 01885 *opacm = Opac_nonIonizing; 01886 *freqm = Freq_nonIonizing; 01887 } 01888 01889 # if 0 01890 /*>>chng 06 aug 12, do not let zones become optically thick to 01891 * peak of dust emission - one of NA's vastly optically thick dust clouds 01892 * did not conserve total energy - very opticallly thick to ir dust emission 01893 * so ir was absorbed and reemitted many times 01894 * this makes sure the cells remain optically thin to dust emission */ 01895 if( gv.lgDustOn && gv.lgGrainPhysicsOn ) 01896 { 01897 double fluxGrainPeak = -1.; 01898 long int ipGrainPeak = -1; 01899 long int ipGrainPeak2 = -1; 01900 01901 i = 0; 01902 while( rfield.anu[i] < 0.0912 && i < (rfield.nflux-2) ) 01903 { 01904 /* >>chng 06 aug 23. Only want to do this test for the IR dust continuum, therefore only 01905 * check on optical depth for wavelengths greater than 1 micron */ 01906 if( gv.GrainEmission[i]*rfield.anu[i]*opac.opacity_abs[i] > fluxGrainPeak ) 01907 { 01908 ipGrainPeak = i; 01909 fluxGrainPeak = gv.GrainEmission[i]*rfield.anu[i]*opac.opacity_abs[i]; 01910 } 01911 ++i; 01912 } 01913 ASSERT( fluxGrainPeak>=0. && ipGrainPeak >=0 ); 01914 01915 /* >>chng 06 aug 23. We have just found the wavelength and flux at the peak of the IR emission. 01916 * Now we want to find the wavelength, short of the peak, which corresponds to 5% of the 01917 * peak emission. This wavelength will be where the code checks to make sure the zone has 01918 * not become optically thick. Since dust opacity generally decreases with wavelength, this assures that 01919 * the optical depths are small for all wavelengths where the flux is appreciable. Tests show that 01920 * this allows for flux/luminosity conservation to within ~5% for an AV of 1e4 mag and at least 2 iterations*/ 01921 i = ipGrainPeak; 01922 /* >>chng 06 aug 23. Only want to do this test for the IR dust continuum, therefore only 01923 * check on opacities for wavelengths shortward of the peak and greater than 1 micron 01924 * this routine can be called with the dust emission totally zero - it is only evaluated 01925 * late to save time. don't do the following tests if peak dust emission is zero */ 01926 if( fluxGrainPeak > 0. ) 01927 { 01928 while( rfield.anu[i] < 0.0912 && i < (rfield.nflux-2) ) 01929 { 01930 /* find wavelength where flux = 5% of peak flux, shortward of the peak */ 01931 if( gv.GrainEmission[i]*rfield.anu[i]*opac.opacity_abs[i] > 0.05*fluxGrainPeak ) 01932 { 01933 /* This will be the array number and flux at 10% of the peak */ 01934 ipGrainPeak2 = i; 01935 } 01936 ++i; 01937 } 01938 ASSERT( ipGrainPeak2>=0 ); 01939 01940 /* use this as a limit to the dr if opacity is larger */ 01941 if( opac.opacity_abs[ipGrainPeak2] > *opacm ) 01942 { 01943 /* scale opacity by factor of 10, which further decreases the zone thickness*/ 01944 *opacm = opac.opacity_abs[ipGrainPeak2]*10.; 01945 *freqm = rfield.anu[ipGrainPeak2]; 01946 /*fprintf(ioQQQ,"DEBUT opac peak %.2e %.2e \n", 01947 *opacm , *freqm );*/ 01948 } 01949 } 01950 } 01951 # endif 01952 01953 { 01954 /* following should be set true to print contributors */ 01955 enum {DEBUG_LOC=false}; 01956 if( DEBUG_LOC ) 01957 { 01958 fprintf(ioQQQ,"conratedebug \t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", 01959 Rate_max_nonIonizing,Freq_nonIonizing,Opac_nonIonizing, 01960 Rate_max_Hion,FreqH ,Opac_Hion,*freqm,*opacm 01961 ); 01962 01963 } 01964 } 01965 01966 /* these were set to -1 at start, and must have been reset in one of the 01967 * two loops. Logic error if still <0. */ 01968 /* >>chng 05 aug 03, change logic to -1 on entry and check at least zero 01969 * here - will be zero if NO radiation field exists, perhaps since totally 01970 * attenuated */ 01971 ASSERT( *opacm >= 0. && *freqm >= 0. ); 01972 return; 01973 } 01974 01975 /*GrainRateDr called by radius_next to find grain heating rate dr */ 01976 STATIC void GrainRateDr(double *grfreqm, 01977 double *gropacm) 01978 { 01979 long int i, 01980 iplow; 01981 double xMax; 01982 01983 DEBUG_ENTRY( "GrainRateDr()" ); 01984 01985 /* in all following changed from anu2 to anu july 25 95 01986 * 01987 * find maximum continuum interaction rate */ 01988 01989 /* not every continuum extends beyond C edge-this whole loop can add to zero 01990 * use total opacity here 01991 * test is to put in fir continuum if free free heating is important */ 01992 if( CoolHeavy.brems_heat_total/thermal.htot > 0.05 ) 01993 { 01994 /* this is pointer to energy where cloud free free optical depth is unity, 01995 * is zero if no freq are opt thin */ 01996 iplow = MAX2(1 , rfield.ipEnergyBremsThin ); 01997 } 01998 else 01999 { 02000 /* do up to carbon photo edge if carbon is turned on */ 02001 /* >>>chng 00 apr 07, add test for whether element is turned on */ 02002 if( dense.lgElmtOn[ipCARBON] ) 02003 { 02004 /* carbon is turned on, use carbon 1 edge */ 02005 iplow = Heavy.ipHeavy[ipCARBON][0]; 02006 } 02007 else 02008 { 02009 /* carbon truned off, use hydrogen balmer continuum */ 02010 iplow = iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2s]; 02011 } 02012 02013 } 02014 02015 xMax = -1.; 02016 /* integrate up to H edge */ 02017 for( i=iplow-1; i < Heavy.ipHeavy[ipHYDROGEN][0]; i++ ) 02018 { 02019 if( rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*opac.opacity_abs[i] > xMax ) 02020 { 02021 xMax = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]* 02022 opac.opacity_abs[i]; 02023 *grfreqm = rfield.anu[i]; 02024 *gropacm = opac.opacity_abs[i]; 02025 } 02026 } 02027 /* integrate up to heii edge if he is turned on, 02028 * this logic will not make sense if grains on but he off, which in itself makes no sense*/ 02029 if( dense.lgElmtOn[ipHELIUM] ) 02030 { 02031 for( i=Heavy.ipHeavy[ipHYDROGEN][0]-1; i < Heavy.ipHeavy[ipHELIUM][1]; i++ ) 02032 { 02033 if( rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*opac.opacity_abs[i] > xMax ) 02034 { 02035 xMax = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]* 02036 opac.opacity_abs[i]; 02037 *grfreqm = rfield.anu[i]; 02038 *gropacm = opac.opacity_abs[i]; 02039 } 02040 } 02041 } 02042 02043 /* possible that there is NO ionizing radiation, in extreme cases, 02044 * if so then xMax will still be negative */ 02045 if( xMax <= 0. ) 02046 { 02047 *gropacm = 0.; 02048 *grfreqm = 0.; 02049 } 02050 return; 02051 } 02052 02053 #undef WIND_CHNG_VELOCITY_RELATIVE 02054 02055 #if 0 02056 /*ChkRate called by radius_next to check how rates of destruction of various species changes */ 02057 STATIC void ChkRate( 02058 /* element number on C scale */ 02059 long int nelem, 02060 /* change in destruction rate */ 02061 double *dDestRate, 02062 /* old and new destruction rates */ 02063 double *DestRateOld, 02064 double *DestRateNew, 02065 /* stage of ionization on the physical scale */ 02066 long int *istage) 02067 { 02068 long int i; 02069 02070 double average, 02071 dDest; 02072 02073 static double OldDest[LIMELM][LIMELM]; 02074 02075 DEBUG_ENTRY( "ChkRate()" ); 02076 02077 /* return if this element is not turned on */ 02078 if( !dense.lgElmtOn[nelem] ) 02079 { 02080 *dDestRate = 1e-3; 02081 *DestRateOld = 1e-3; 02082 *DestRateNew = 1e-3; 02083 *istage = 0; 02084 return; 02085 } 02086 02087 /* for first zone, and during search, we will do nothing but 02088 * still must return finite numbers */ 02089 *istage = 1; 02090 *dDestRate = 0.; 02091 *DestRateOld = 0.; 02092 *DestRateNew = 0.; 02093 *dDestRate = 0.; 02094 02095 if( nzone <= 1 ) 02096 { 02097 for( i=0; i < nelem+1; i++ ) 02098 { 02099 OldDest[nelem][i] = ionbal.RateIonizTot[nelem][i]; 02100 } 02101 } 02102 02103 else if( dense.xIonDense[nelem][0]/dense.gas_phase[nelem] < 0.9 ) 02104 { 02105 /* do not use this method if everything is atomic */ 02106 for( i=0; i < (nelem); i++ ) 02107 { 02108 /* last check below, .5 chosen so that do not key off 02109 * predominantly neutral species where self-opacity 02110 * could cause oscillations */ 02111 if( ((dense.xIonDense[nelem][i]/dense.gas_phase[nelem] > 0.01 && 02112 dense.xIonDense[nelem][i]/dense.gas_phase[nelem] < 0.9) && 02113 dense.xIonDense[nelem][i+1]/dense.gas_phase[nelem] > .05) && 02114 OldDest[nelem][i] > 0. ) 02115 { 02116 /* last check on old dest in case just lowered ionization 02117 * stage, so no history */ 02118 /* check that rate is positive */ 02119 if( ionbal.RateIonizTot[nelem][i] <= 0. ) 02120 { 02121 fprintf( ioQQQ, " ChkRate gets insane destruction rate for ion%4ld%4ld%10.2e\n", 02122 nelem+1, i, ionbal.RateIonizTot[nelem][i] ); 02123 cdEXIT(EXIT_FAILURE); 02124 } 02125 02126 /* do not consider unless of middling ionization, and 02127 * rate is going down (to prevent dr osciallating) 02128 * no absolute value in following since do not want to 02129 * consider cases where ionization rate increases */ 02130 average = (OldDest[nelem][i] + ionbal.RateIonizTot[nelem][i])* 0.5; 02131 02132 dDest = (OldDest[nelem][i] - ionbal.RateIonizTot[nelem][i])/ average; 02133 /* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ term + if rate going down */ 02134 02135 if( *dDestRate < dDest ) 02136 { 02137 /* biggest rate so far, remember change in rates and ionization stage */ 02138 *dDestRate = dDest; 02139 *istage = i+1; 02140 *DestRateOld = OldDest[nelem][i]; 02141 *DestRateNew = ionbal.RateIonizTot[nelem][i]; 02142 } 02143 02144 } 02145 OldDest[nelem][i] = ionbal.RateIonizTot[nelem][i]; 02146 } 02147 } 02148 return; 02149 } 02150 #endif