cloudy trunk
Loading...
Searching...
No Matches
radius_next.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3/*radius_next use adaptive logic to find next zone thickness */
4/*ContRate called by radius_next to find energy of maximum continuum-gas interaction */
5/*GrainRateDr called by radius_next to find grain heating rate dr */
14#include "cddefines.h"
15#include "lines_service.h"
16#include "iso.h"
17#include "geometry.h"
18#include "h2.h"
19#include "mole.h"
20#include "hyperfine.h"
21#include "opacity.h"
22#include "dense.h"
23#include "heavy.h"
24#include "grainvar.h"
25#include "elementnames.h"
26#include "rfield.h"
27#include "dynamics.h"
28#include "thermal.h"
29#include "hmi.h"
30#include "coolheavy.h"
31#include "timesc.h"
32#include "doppvel.h"
33#include "stopcalc.h"
34#include "colden.h"
35#include "phycon.h"
36#include "rt.h"
37#include "trace.h"
38#include "wind.h"
39#include "save.h"
40#include "taulines.h"
41#include "pressure.h"
42#include "iterations.h"
43#include "struc.h"
44#include "radius.h"
45#include "dark_matter.h"
46
47/*ContRate called by radius_next to find energy of maximum continuum-gas interaction */
48STATIC void ContRate(double *freqm,
49 double *opacm);
50
51/*GrainRateDr called by radius_next to find grain heating rate dr */
52STATIC void GrainRateDr(double *grfreqm,
53 double *gropacm);
54
55/*radius_next use adaptive logic to find next zone thickness
56 * return 0 if ok, 1 for abort */
57int radius_next(void)
58{
59 static double OHIIoHI,
60 OldHeat = 0.,
61 OldTe = 0.,
62 OlddTdStep = 0.,
63 OldWindVelocity = 0.,
64 Old_H2_heat_cool;
65
66 const double BigRadius = 1e30;
67
68 DEBUG_ENTRY( "radius_next()" );
69
70 /*--------------------------------------------------------------
71 *
72 * this sub determines the thickness of the next zone
73 * if is called one time for each zone
74 *
75 *-------------------------------------------------------------- */
76
77 if( trace.lgTrace )
78 fprintf( ioQQQ, " radius_next called\n" );
79
80 /* >>chng 03 sep 21 - decide whether this is the first call */
81 bool lgFirstCall = ( nzone == 0 );
82
83 multimap<double,string> drChoice;
84
85 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
86 * check on change in hydrogen ionization */
87 if( !lgFirstCall && OHIIoHI > 1e-15 &&
88 (dense.xIonDense[ipHYDROGEN][1] > 0. && dense.xIonDense[ipHYDROGEN][0] > 0.) )
89 {
90 double atomic_frac = dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0];
91 double drHion;
92 /* ratio of current HII/HI to old value - < 1 when becoming more neutral */
93 /* this is now change in atomic fraction */
94 double x = 1. - atomic_frac/OHIIoHI;
95 if( atomic_frac > 0.05 && atomic_frac < 0.9 )
96 {
97 /* >>chng 96 feb 23 from 0.7 to 0.3 because punched through H i-front
98 * >>chng 97 may 5, from 0.3 to 0.2 since punched through H i-front */
99 /* >>chng 02 jun 05 from 0.2 to 0.05 poorly resolved i-front, also added two-branch logic*/
100 drHion = radius.drad*MAX2( 0.2 , 0.05/MAX2(1e-10,x) );
101 }
102 else if( x > 0. )
103 {
104 /* >>chng 96 feb 23 from 0.7 to 0.3 because punched through H i-front
105 * >>chng 97 may 5, from 0.3 to 0.2 since punched through H i-front */
106 drHion = radius.drad*MAX2( 0.2 , 0.2/MAX2(1e-10,x) );
107 }
108 else
109 {
110 drHion = BigRadius;
111 }
112 double SaveOHIIoHI = OHIIoHI;
113 OHIIoHI = dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0];
114 ostringstream oss;
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() ) );
118 }
119 else
120 {
121 if( dense.xIonDense[ipHYDROGEN][1] > 0. && dense.xIonDense[ipHYDROGEN][0] > 0. )
122 OHIIoHI = dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0];
123 else
124 OHIIoHI = 0.;
125 }
126
127 if( rt.dTauMase < -0.01 )
128 {
129 /* maser so powerful, must limit inc in tay
130 * >>chng 96 aug 08, from 0.3 to 0.1 due to large maser crash */
131 double drHMase = radius.drad*MAX2(0.1,-0.2/rt.dTauMase);
132 ostringstream oss;
133 // NB NB -- DON'T ALTER THIS STRING, setting rt.lgMaserSetDR below keys from this!
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() ) );
137 }
138
139 /* check change in relative ionization - this is to insure a gradual approach to
140 * ionization fronts - were dr allowed to remain constant we would crash through
141 * a front in a single zone */
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;
146 if( nzone > 4 )
147 {
148 long ichange_heavy_nelem = -1L;
149 long ichange_heavy_ion = -1L;
150 double dr_change_heavy = BigRadius;
151 for( int nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
152 {
153 if( dense.lgElmtOn[nelem] )
154 {
155 realnum change;
156 /* will only track ions with fractional abundances greater than this
157 * He and C are special cases because of special roles in PDRs */
158 realnum frac_limit;
159 if( nelem<=ipHELIUM || nelem==ipCARBON )
160 {
161 frac_limit = 1e-4f;
162 change = CHANGE_ION_HHE;
163 }
164 else
165 {
166 /* struc.dr_ionfrac_limit set in init_coreload, is limit to fractional
167 * abundances of ions that are used in prt_comment to indicate oscillations
168 * or rapid change in ionization */
169 frac_limit = struc.dr_ionfrac_limit/2.f;
170 change = CHANGE_ION_HEAV;
171 }
172
173 for( int ion=0; ion<=nelem+1; ++ion )
174 {
175 realnum abundnew = dense.xIonDense[nelem][ion]/SDIV( dense.gas_phase[nelem]);
176 realnum abundold = struc.xIonDense[nelem][ion][nzone-3]/SDIV( struc.gas_phase[nelem][nzone-3]);
177 if( abundnew > frac_limit && abundold > frac_limit )
178 {
179 /* this test is not done when [nzone-n] <0 */
180 realnum abundolder = struc.xIonDense[nelem][ion][nzone-4]/SDIV( struc.gas_phase[nelem][nzone-4]);
181 realnum abundoldest = struc.xIonDense[nelem][ion][nzone-5]/SDIV( struc.gas_phase[nelem][nzone-5]);
182 /* this is fractional change, use min of old and new abundances, to pick up
183 * rapid changing Ca+ sooner */
184 double change_heavy_frac = fabs(abundnew-abundold)/MIN2(abundold,abundnew);
185 /* want fractional change to be less than this factor */
186 if( (change_heavy_frac > change) && (change_heavy_frac > change_heavy_frac_big) &&
187 /* >>chng 03 dec 07, add test that abund is not oscillating */
188 /* also test that abundance is increasing - we are headed into a front */
189 ((abundnew-abundold)>0.) &&
190 ((abundold-abundolder)>0.) &&
191 ((abundolder-abundoldest)>0.) )
192 {
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;
197 /* factor in max has been adjusted to prevent code from running
198 * through various ionization fronts.
199 * >>chng 03 dec 06, from min of 0.5 to min of 0.25, crash into He i-front
200 * in hizqso.in */
201 /* >>chng 04 mar 03, min had become 0.1, forced oscillations in nova.in
202 * in silicon, zoning changed greatly, causing change in diffuse lin
203 * pumping. put back to 0.25 */
204 dr_change_heavy = radius.drad * MAX2(0.25, change / change_heavy_frac );
205 }
206 }
207 }
208 }
209 }
210 /* set to -1 before loop over ions, if no significant changes in ionization occurred
211 * then may still be -1
212 */
213 ostringstream oss;
214 if(ichange_heavy_nelem>=0)
215 {
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;
219 }
220 else
221 {
222 oss << "none";
223 dr_change_heavy = BigRadius;
224 }
225 drChoice.insert( pair<const double,string>( dr_change_heavy, oss.str() ) );
226 }
227
228 /* if Leiden hacks are on then use increase in dust extinction as control
229 * >>chng 05 aug 13, add this */
230 if( mole_global.lgLeidenHack )
231 {
232 /* Draine field is only defined over narrow range in FUV - must not let change
233 * in extinction become too large -
234 * prefactor is change in optical depth */
235 double drLeiden_hack = MAX2( 0.02 , 0.05*rfield.extin_mag_V_point ) /
236 SDIV( geometry.FillFac * rfield.opac_mag_V_point );
237 drChoice.insert( pair<const double,string>( drLeiden_hack, "Leiden hack" ) );
238 }
239
240 // limit to dust optical depth per zone
241 static double extin_mag_V_point_old = -1.;
242 if( nzone>1 )
243 {
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;
246 if( dExtin > SMALLFLOAT )
247 {
248 double drExtintion = radius.drad / dExtin;
249 ostringstream oss;
250 oss << "change in V mag extinction; current=" << scientific << setprecision(3) <<
251 rfield.extin_mag_V_point;
252 oss << fixed << " delta=" << dExtin;
253 drChoice.insert( pair<const double,string>( drExtintion, oss.str() ) );
254 }
255
256 }
257 extin_mag_V_point_old = rfield.extin_mag_V_point;
258
259 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
260 /* check how heating is changing */
261 if( nzone > 1 && !thermal.lgTemperatureConstant &&
262 /* if density fluctuations in place then override change in heat for dr set */
263 !( strcmp(dense.chDenseLaw,"SINE") == 0 && dense.lgDenFlucOn ) )
264 {
265 double dHdStep = fabs(thermal.htot-OldHeat)/thermal.htot;
266 if( dHdStep > 0. )
267 {
268 double drdHdStep;
269 if( dense.gas_phase[ipHYDROGEN] >= 1e13 )
270 {
271 drdHdStep = radius.drad*MAX2(0.8,0.075/dHdStep);
272 }
273 else if( dense.gas_phase[ipHYDROGEN] >= 1e11 )
274 {
275 drdHdStep = radius.drad*MAX2(0.8,0.100/dHdStep);
276 }
277 else
278 {
279 /* changed from .15 to .12 for outer edge of coolhii too steep dT
280 * changed to .10 for symmetry, big change in some rates, 95nov14
281 * changed from .10 to .125 since parispn seemed to waste zones
282 * >>chng 96 may 21, from .125 to .15 since pn's still waste zones */
283 drdHdStep = radius.drad*MAX2(0.8,0.15/dHdStep);
284 }
285 ostringstream oss;
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() ) );
289 }
290 }
291 OldHeat = thermal.htot;
292
293 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
294 /* pressure due to incident continuum if in eos */
295 if( strcmp(dense.chDenseLaw,"CPRE") == 0 && pressure.lgContRadPresOn )
296 {
297 if( nzone > 2 && pressure.pinzon > 0. )
298 {
299 /* pinzon is pressrue from acceleration onto previos zone
300 * want this less than some fraction of total pressure */
301 /* >>chng 06 feb 01, change from init pres to current total pressure
302 * in const press high U ulirgs current pressure may be quite larger
303 * than init pressure due to continuum absorption */
304 double drConPres = 0.05*pressure.PresTotlCurr/
305 (wind.AccelTotalOutward*dense.xMassDensity*geometry.FillFac);
306 drChoice.insert( pair<const double,string>( drConPres, "change in con accel" ) );
307 }
308 }
309
310 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
311 /* gravitational pressure */
312 if( strcmp(dense.chDenseLaw,"CPRE") == 0 && (dark.lgNFW_Set || pressure.gravity_symmetry>=0) )
313 {
314 if( fabs( pressure.RhoGravity ) > 0. )
315 {
316 double drGravPres = 0.05 * pressure.PresTotlCurr * radius.drad / fabs( pressure.RhoGravity );
317 ASSERT( drGravPres > 0. );
318 drChoice.insert( pair<const double,string>( drGravPres, "change in gravitational pressure" ) );
319 }
320 }
321
322 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
323 /* check how temperature is changing */
324 double dTdStep = FLT_MAX;
325 if( nzone > 1 )
326 {
327 /* change in temperature; current= */
328 dTdStep = (phycon.te-OldTe)/phycon.te;
329 /* >>chng 02 dec 08, desired change in temperature per zone must not
330 * be smaller than allower error in temperature. For now use relative error
331 * in heating - - cooling balance. Better would be to also use c-h deriv wrt t
332 * to get slope */
333 /* >>chng 02 dec 11 rjrw change back to 0.03 -- improve dynamics.dRad criterion instead */
334 double x = 0.03;
335 /* >>chng 02 dec 07, do not do this if there is mild te jitter,
336 * so that dT is changing sign - this happens
337 * in ism.ini, very stable temperature with slight noise up and down */
338 if( dTdStep*OlddTdStep > 0. )
339 {
340 /* >>chng 96 may 30, variable depending on temp
341 * >>chng 96 may 31, allow 0.7 smaller, was 0.8
342 * >>chng 97 may 05, from 0.7 to 0.5 stop from punching through thermal front
343 * >>chng 04 mar 30, from 0.7 to 0.5 stop from punching through thermal front,
344 * for some reason factor was 0.7, not 0.5 as per previous change */
345 double absdt = fabs(dTdStep);
346 double drdTdStep = radius.drad*MAX2(0.5,x/absdt);
347 ostringstream oss;
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() ) );
351 }
352 }
353 OlddTdStep = dTdStep;
354 OldTe = phycon.te;
355
356 /* >>chng 02 oct 06, only check on opacity - interaction if not
357 * constant temperature - there were constant temperature models that
358 * extended to infinity but hung with last few photons and this test.
359 * better to ignore this test which is really for thermal models */
360 /* >>chng 06 mar 20, do not call if recombination logic in place */
361 if( !thermal.lgTemperatureConstant && !dynamics.lgRecom )
362 {
363 double freqm = 0., opacm = 0.;
364 /* find freq where opacity largest and interaction rate is fastest */
365 ContRate(&freqm,&opacm);
366
367 /* use optical depth at max interaction energy
368 * >>chng 96 jun 06 was drChange=0.15 changed to 0.3 for high Z models
369 * taking too many zones
370 * drInter = drChange / MAX(1e-30,opacm*FillFac) */
371
372 double drInter = 0.3/MAX2(1e-30,opacm*geometry.FillFac*geometry.DirectionalCosin);
373 ostringstream oss;
374 oss << "cont inter nu=" << scientific << setprecision(2) << freqm << " opac=" << opacm;
375 drChoice.insert( pair<const double,string>( drInter, oss.str() ) );
376 }
377
378 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
379 /* check whether change in wind velocity constrains DRAD
380 * WJH 22 May 2004: disable when we are near the sonic point since
381 * the velocity may be jumping all over the place but we just want
382 * to push through it as quickly as we can */
383 if( !wind.lgStatic() && !pressure.lgSonicPoint && !pressure.lgStrongDLimbo )
384 {
385 double v = fabs(wind.windv);
386 /* this is relative change in velocity */
387 double dVelRelative = fabs(wind.windv-OldWindVelocity)/
388 MAX2(v,0.1*timesc.sound_speed_isothermal);
389
390 const double WIND_CHNG_VELOCITY_RELATIVE = 0.01;
391 /*fprintf(ioQQQ,"DEBUG rad %.3f vel %.2e dVelRelative %.2e",
392 log10(radius.Radius) , wind.windv , dVelRelative );*/
393
394 /* do not use this on first zone since do not have old velocity */
395 double winddr;
396 if( dVelRelative > WIND_CHNG_VELOCITY_RELATIVE && nzone>1 )
397 {
398 /* factor less than one if change larger than WIND_CHNG_VELOCITY_RELATIVE */
399 double factor = WIND_CHNG_VELOCITY_RELATIVE / dVelRelative;
400 /*fprintf(ioQQQ," DEBUG factor %.2e", factor );*/
401 factor = MAX2(0.8 , factor );
402 winddr = radius.drad * factor;
403 }
404 else
405 {
406 winddr = BigRadius;
407 }
408
409 /* set dr from advective term */
410 if( dynamics.lgAdvection && iteration > dynamics.n_initial_relax)
411 {
412 winddr = MIN2( winddr , dynamics.dRad );
413 /*>>chng 04 oct 06, set dVelRelative to dynamics.dRad since dVelRelative is printed as part
414 * of reason for choosing this criteria, want to reflect valid reason. */
415 dVelRelative = dynamics.dRad;
416 }
417 ostringstream oss;
418 oss << "Wind, dVelRelative=" << scientific << setprecision(3);
419 oss << dVelRelative << " windv=" << wind.windv;
420 drChoice.insert( pair<const double,string>( winddr, oss.str() ) );
421 }
422 /* remember old velocity */
423 OldWindVelocity = wind.windv;
424
425 const double DNGLOB = 0.10;
426
427 /* inside out globule */
428 if( strcmp(dense.chDenseLaw,"GLOB") == 0 )
429 {
430 if( radius.glbdst < 0. )
431 {
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",
434 radius.glbdst );
436 }
437 double factor = radius.glbden*pow(radius.glbrad/radius.glbdst,radius.glbpow);
438 double fac2 = radius.glbden*pow(radius.glbrad/(radius.glbdst - (realnum)radius.drad),radius.glbpow);
439 /* DNGLOB is relative change in density allowed this zone, 0.10 */
440 double GlobDr = ( fac2/factor > 1. + DNGLOB ) ? radius.drad*DNGLOB/(fac2/factor - 1.) : BigRadius;
441 GlobDr = MIN2(GlobDr,radius.glbdst/20.);
442 ostringstream oss;
443 oss << "GLOB law, HDEN=" << scientific << setprecision(2) << dense.gas_phase[ipHYDROGEN];
444 drChoice.insert( pair<const double,string>( GlobDr, oss.str() ) );
445 }
446
447 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
448 double hdnew = 0.;
449 if( strncmp( dense.chDenseLaw , "DLW" , 3) == 0 )
450 {
451 /* one of the special density laws, first get density at possible next dr */
452 if( strcmp(dense.chDenseLaw,"DLW1") == 0 )
453 {
454 hdnew = dense_fabden(radius.Radius+radius.drad,radius.depth+
455 radius.drad);
456 }
457 else if( strcmp(dense.chDenseLaw,"DLW2") == 0 )
458 {
459 hdnew = dense_tabden(radius.Radius+radius.drad,radius.depth+
460 radius.drad);
461 }
462 else if( strcmp(dense.chDenseLaw,"DLW3") == 0 )
463 {
464 hdnew = dense_parametric_wind(radius.Radius+radius.drad);
465 }
466 else
467 {
468 fprintf( ioQQQ, " dlw insanity in radius_next\n" );
470 }
471 double drTab = fabs(hdnew-dense.gas_phase[ipHYDROGEN])/MAX2(hdnew,dense.gas_phase[ipHYDROGEN]);
472 drTab = radius.drad*MAX2(0.2,0.10/MAX2(0.01,drTab));
473 ostringstream oss;
474 oss << "spec den law, new old den " << scientific << setprecision(2);
475 oss << hdnew << " " << dense.gas_phase[ipHYDROGEN];
476 drChoice.insert( pair<const double,string>( drTab, oss.str() ) );
477 }
478
479 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
480 /* special density law */
481 if( strcmp(dense.chDenseLaw,"DLW1") == 0 )
482 {
483 double dnew = fabs(dense_fabden(radius.Radius+radius.drad,radius.depth+
484 radius.drad)/dense.gas_phase[ipHYDROGEN]-1.);
485 /* DNGLOB is relative change in density allowed this zone, 0.10 */
486 double SpecDr;
487 if( dnew == 0. )
488 {
489 SpecDr = radius.drad*3.;
490 }
491 else if( dnew/DNGLOB > 1.0 )
492 {
493 SpecDr = radius.drad/(dnew/DNGLOB);
494 }
495 else
496 {
497 SpecDr = BigRadius;
498 }
499 ostringstream oss;
500 oss << "special law, HDEN=" << scientific << setprecision(2) << dense.gas_phase[ipHYDROGEN];
501 drChoice.insert( pair<const double,string>( SpecDr, oss.str() ) );
502 }
503
504 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
505 /* check grain line heating dominates
506 * this is important in PDR and HII region calculations
507 * >>chng 97 jul 03, added following check */
508 if( thermal.heating[0][13]/thermal.htot > 0.2 )
509 {
510 double grfreqm = 0., gropacm = 0.;
511 /* >>chng 01 jan 03, following returns 0 when NO light at all,
512 * had failed with botched monitor */
513 GrainRateDr(&grfreqm,&gropacm);
514 double DrGrainHeat = 1.0/MAX2(1e-30,gropacm*geometry.FillFac*geometry.DirectionalCosin);
515 ostringstream oss;
516 oss << "grain heating nu=" << scientific << setprecision(2) << grfreqm << " opac=" << gropacm;
517 drChoice.insert( pair<const double,string>( DrGrainHeat, oss.str() ) );
518 }
519
520 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
521 /* check if line heating dominates
522 * this is important in high metallicity models */
523 if( thermal.heating[0][22]/thermal.htot > 0.2 )
524 {
525 long level = -1;
526 TransitionProxy t = FndLineHt(&level);
527
528 if( t.Coll().heat()/thermal.htot > 0.1 )
529 {
530 ASSERT( t.associated() );
531
532 double TauInwd = t.Emis().TauIn();
533 double DopplerWidth = t.Emis().dampXvel()/t.Emis().damp();
534 double TauDTau = t.Emis().PopOpc()*t.Emis().opacity()/DopplerWidth;
535
536 fixit(); // all other line stacks need to be included here.
537 // can we just sweep over line stack? Is that ready yet?
538
539 /* in following logic cannot use usual inverse opacity,
540 * since line heating competes with escape probability,
541 * so is effective at surprising optical depths */
542 double drLineHeat = ( TauDTau > 0. ) ? MAX2(1.,TauInwd)*0.4/TauDTau : BigRadius;
543 ostringstream oss;
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() ) );
547 }
548 }
549
550 /* do not let change in cooling/heating due to H2 become too large */
551 if( lgFirstCall )
552 Old_H2_heat_cool = 0.;
553 else if( !thermal.lgTemperatureConstant )
554 {
555 /* this is case where temperature has not been set */
556 /* compare total heating - cooling due to h2 with total due to everything */
557 double H2_heat_cool = (fabs(hmi.HeatH2Dexc_used)+fabs(hmi.HeatH2Dish_used)) / thermal.htot;
558 // Leiden hack PDR sims can have H2 heating set by sims equation, does not change
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 )
561 {
562 /* >>chng 05 oct 27, had been taking 20% of original radius - this caused zoning
563 * to become very fine and may not have been needed - change from 0.2 to 0.5 */
564 double drH2_heat_cool = radius.drad*MAX2(0.5,0.05/dH2_heat_cool);
565 ostringstream oss;
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() ) );
569 }
570 }
571 Old_H2_heat_cool = (fabs(hmi.HeatH2Dexc_used)+fabs(hmi.HeatH2Dish_used)) / thermal.htot;
572
573 /* >>chng 03 mar 04, add this logic */
574 /* do not let change in H2 and heavy element molecular abundances become too large
575 * "change in heav ele mole abundance, d(mole)/elem" */
576 int mole_dr_change = -1;
577 if( nzone >= 4 )
578 {
579 /* first do H2 abundance */
580 double Old_H2_abund;
581 /* >>chng 04 dec 15, do not use special logic when large h2 is turned on */
582 int i;
583 Old_H2_abund = struc.H2_abund[nzone-3];
584 /* >>chng 03 jun 16, limit from 0.01 to 0.001, some models fell over H2 front due to
585 * large zoning, when large H2 was just this caused oscillations in solomon process */
586 /* >>chng 03 nov 22, from > 0.001 to > 3e-4, models that start almost in H2 have
587 * rapid increase in H2 at shallow depths, start sensing this sooner */
588 /* >>chng 03 dec 10, from 3e-4 to 1e-4, to get smaller chagned in leiden1.in test */
589 /* radius_next keys from change in H2 abundance, d(H2)/H */
590 /* >>chng 04 mar 11, start sensing H2 earlier since otherwise step size
591 * needs to become way too small tooo quickly. limit changed from 1e-4 to 1e-6 */
592 /* >>chng 04 jun 29, fromo > 1e-6 to >1e-8, some pdr's had too large a change in H2 */
593
594 if( 2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN] > 1e-8 )
595 {
596 double fac = 0.2;
597 /* this is percentage change in H2 density - "change in H2 abundance" */
598 double dH2_abund = 2.*fabs( hmi.H2_total - Old_H2_abund ) / hmi.H2_total;
599 /* in testing th85ism the dH2_abund did come out zero */
600 /* >>chng 03 jun 16, change d(H2) from 0.05 to 0.1, fine resolution of H2/H exposed
601 * small numerical oscillations in Solomon process */
602 /* >>chng 03 nov 22, smallest possible ratio of dr(next)/dr changed from
603 * 0.2 to 0.05, models that started almost in H2 front need to be able to sense it */
604 /* >>chng 04 mar 15, with such small possible changes in dr, only 0.05, a thermal front
605 * can easily cause large changes in H2 and T that are not due to zoning, but to the
606 * discontinuity. make smallest change larger to prevent hald due to dr */
607 /* >>chng 05 aug 23, thermal front allowed dr to become much too small
608 * chng from 0.02 to .6 */
609 dH2_abund = SDIV(dH2_abund);
610 double drH2_abund = radius.drad*MAX2(0.2,fac/dH2_abund );
611 ostringstream oss;
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() ) );
614 }
615
616 /* check how molecular fractions of all heavy elements are chaning relative
617 * to total gas phase abundance */
618 double max_change = 0.0;
619 /* >>chng 04 jun 02, upper limit had been all species but now limit to real
620 * molecules since do not want this logic to work with the ions */
621 for( i=0; i < mole_global.num_calc; ++i )
622 {
623 realnum abund,
624 abund1,
625 abund_limit;
626 if( !mole_global.list[i]->isMonatomic() && mole_global.list[i]->parentLabel.empty() )
627 {
628 /* >>chng 03 sep 21, add CO logic */
629 /* >>chng 04 mar 30, generalize to any molecule at all */
630 /* >>chng 04 mar 31 lower limit to abund had been 0.01, change
631 * to 0.001 to pick up approach to molecular fronts */
632 abund = 0.;
633 for( molecule::nAtomsMap::iterator atom=mole_global.list[i]->nAtom.begin();
634 atom != mole_global.list[i]->nAtom.end(); ++atom)
635 {
636 long int nelem = atom->first->el->Z-1;
637 if (dense.lgElmtOn[nelem])
638 {
639 abund1 = (realnum)mole.species[i].den*atom->second/SDIV(dense.gas_phase[nelem]);
640 if (abund1 > abund)
641 {
642 abund = abund1;
643 }
644 }
645 }
646 /* is this an ice? need special logic for ice since density increases
647 * exponentially when grain temp falls below sublimation temperature
648 * >>chng 05 dec 06 - detect changes for smaller abundances for ices
649 * due to large changes in ice abundances */
650 if( mole_global.list[i]->lgGas_Phase )
651 {
652 abund_limit = 1e-3f;
653 }
654 else
655 {
656 /* this is an ice - track its abundance at lower values so that
657 * we resolve the sublimation transition region */
658 abund_limit = 1e-5f;
659 }
660
661 if( abund > abund_limit )
662 {
663 /* the relative change in the abundance */
664 double relative_change =
665 2.0 * fabs( mole.species[i].den - struc.molecules[i][nzone-3] ) /
666 SDIV ( mole.species[i].den + struc.molecules[i][nzone-3] ) ;
667 if (relative_change > max_change)
668 {
669 max_change = relative_change;
670 mole_dr_change = i;
671 }
672 }
673 }
674 }
675 if( mole_dr_change >= 0 )
676 {
677 double dr_mole_abund = radius.drad * MAX2(0.6, 0.05/SDIV(max_change));
678 ostringstream oss;
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() ) );
683 }
684 }
685
686 /* some consideration due to big H2 molecule */
687 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
688 drChoice.insert( pair<const double,string>( (*diatom)->H2_DR(), "change in big H2 Solomon rate line opt depth" ) );
689
690 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
691 /* can't make drmax large deep in model since causes feedback
692 * oscillations with changes in heating or destruction rates
693 * >>chng 96 oct 15, change from 5 to 10 */
694 /* >>chng 96 nov 11, had been 4 * drad up to 11, change to following
695 * to be similar to c90.01, max of 1.3 between 5 and 11
696 * >>chng 04 oct 29 geometry.DirectionalCosin was ioncorrect applied to this factor */
697 double drmax = ( nzone < 5 ) ? 4.*radius.drad : 1.3*radius.drad;
698 drChoice.insert( pair<const double,string>( drmax, "DRMAX" ) );
699
700 /* time dependent recombination - conditions become very homogeneous and
701 * crash into I fronts - use last iteration's radius as guess of current case */
702 double recom_dr_last_iter = BigRadius;
703 if( dynamics.lgTimeDependentStatic && dynamics.lgRecom )
704 {
705 /* first time through nzone == 1 */
706 static long int nzone_recom = -1;
707 if( nzone <= 1 )
708 {
709 /* initialization case */
710 nzone_recom = 0;
711 }
712 else if( radius.depth < struc.depth_last[struc.nzonePreviousIteration-1] &&
713 nzone_recom < struc.nzonePreviousIteration )
714 {
715 ASSERT( nzone_recom>=0 && nzone_recom<struc.nzonePreviousIteration );
716 /* case where we are within previous computed structure
717 * first possibly increase nzone_recom so that it points deeper
718 * than this zone */
719 while( struc.depth_last[nzone_recom] < radius.depth &&
720 nzone_recom < struc.nzonePreviousIteration-1 )
721 ++nzone_recom;
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" ) );
725 }
726 }
727
728 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
729 /* change in electron density - radius_next keys from change in elec den,
730 * remember old electron density during first call */
731 /* this is low ionization solution */
732 if( nzone > 2 )
733 {
734 /* next is-2 since nzone is on physics not c scale, and we want previous zone */
735 double Efrac_old = struc.ednstr[nzone-3]/struc.gas_phase[ipHYDROGEN][nzone-3];
736 double Efrac_new = dense.eden/dense.gas_phase[ipHYDROGEN];
737 double dEfrac = fabs(Efrac_old-Efrac_new) / Efrac_new;
738 double drEfrac;
739 if( dEfrac > SMALLFLOAT )
740 {
741 double fac = 0.04;
742 /* >>chng 03 dec 25, use smaller rel change in elec frac when most elec in ipMole or grains */
743 /* >>chng 04 sep 14, change to from from metals but comment out */
744 /* >>chng 04 sep 17, change to from from metals - uncomment */
745 if( dense.eden_from_metals > 0.5 )
746 {
747 fac = 0.04;
748 }
749 /* >>chng 04 feb 23, add test on hydrogen being predomintly
750 * recombined due to three-body recom, which is very sensitive
751 * to the electron density - but only do this in partially ionized medium */
752 else if( iso_sp[ipH_LIKE][ipHYDROGEN].RecomCollisFrac > 0.8 &&
753 dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN] > 0.1 &&
754 dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN] < 0.8 )
755
756 {
757 fac = 0.02;
758 }
759 /* >>chng 04 sep 17, change to 0.1 from 0.2 */
760 drEfrac = radius.drad*MAX2(0.1,fac/dEfrac);
761 }
762 else
763 {
764 drEfrac = BigRadius;
765 }
766 ostringstream oss;
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() ) );
770 }
771
772 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
773 /* do not let thickness get too large */
774 if( nzone > 20 )
775 {
776 /* >>chng 02 nov 05, change from 1/20 to 1/10 wasted zones early on */
777 drChoice.insert( pair<const double,string>( radius.depth/10., "relative depth" ) );
778 }
779
780 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
781 /* case where stopping thickness or edge specified, need to approach slowly */
782 double thickness_total = BigRadius;
783 double DepthToGo = BigRadius;
784 if( StopCalc.HColStop < 5e29 )
785 {
786 double coleff = SDIV( dense.gas_phase[ipHYDROGEN]*geometry.FillFac);
787 DepthToGo = MIN2(DepthToGo ,
788 (StopCalc.HColStop-colden.colden[ipCOL_HTOT]) / coleff );
789 /* >>chng 03 oct 28, forgot to div col den above by eff density */
790 thickness_total = MIN2(thickness_total , StopCalc.HColStop / coleff );
791 }
792
793 if( StopCalc.colpls < 5e29 )
794 {
795 double coleff = (double)SDIV( (dense.xIonDense[ipHYDROGEN][1])*geometry.FillFac);
796 DepthToGo = MIN2(DepthToGo ,
797 (StopCalc.colpls-colden.colden[ipCOL_Hp]) / coleff );
798 thickness_total = MIN2(thickness_total , StopCalc.colpls / coleff );
799 }
800
801 if( StopCalc.col_h2 < 5e29 )
802 {
803 /* >>chng 03 apr 15, add this molecular hydrogen */
804 double coleff = (double)SDIV( hmi.H2_total*geometry.FillFac);
805 DepthToGo = MIN2(DepthToGo ,
806 (StopCalc.col_h2-colden.colden[ipCOL_H2g]-colden.colden[ipCOL_H2s]) / coleff );
807 thickness_total = MIN2(thickness_total , StopCalc.col_h2 / coleff );
808 }
809
810 if( StopCalc.col_h2_nut < 5e29 )
811 {
812 /* >>chng 03 apr 15, add this molecular hydrogen */
813 double coleff = (double)SDIV( (2*hmi.H2_total+dense.xIonDense[ipHYDROGEN][0])*geometry.FillFac);
814 DepthToGo = MIN2(DepthToGo ,
815 (StopCalc.col_h2_nut-(2*(colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s])+colden.colden[ipCOL_H0])) / coleff );
816 thickness_total = MIN2(thickness_total , StopCalc.col_h2_nut / coleff );
817 }
818
819 if( StopCalc.col_H0_ov_Tspin < 5e29 )
820 {
821 /* >>chng 05 jan 09, add n(H0)/Tspin */
822 double coleff = (double)SDIV( dense.xIonDense[ipHYDROGEN][0] / hyperfine.Tspin21cm*geometry.FillFac );
823 DepthToGo = MIN2(DepthToGo ,
824 (StopCalc.col_H0_ov_Tspin - colden.H0_ov_Tspin) / coleff );
825 thickness_total = MIN2(thickness_total , StopCalc.col_H0_ov_Tspin / coleff );
826 }
827
828 if( StopCalc.col_monoxco < 5e29 )
829 {
830 /* >>chng 03 apr 15, add this, CO */
831 double coleff = (double)SDIV( (findspecieslocal("CO")->den)*geometry.FillFac);
832 DepthToGo = MIN2(DepthToGo ,
833 (StopCalc.col_monoxco-findspecieslocal("CO")->column) / coleff );
834 thickness_total = MIN2(thickness_total , StopCalc.col_monoxco / coleff );
835 }
836
837 if( StopCalc.colnut < 5e29 )
838 {
839 double coleff = (double)SDIV( (dense.xIonDense[ipHYDROGEN][0])*geometry.FillFac);
840 DepthToGo = MIN2(DepthToGo ,
841 (StopCalc.colnut - colden.colden[ipCOL_H0]) / coleff );
842 thickness_total = MIN2(thickness_total , StopCalc.colnut / coleff );
843 }
844
845 /* this is case where outer radius is set */
846 if( radius.StopThickness[iteration-1] < 5e29 )
847 {
848 thickness_total = MIN2(thickness_total , radius.StopThickness[iteration-1] );
849 DepthToGo = MIN2(DepthToGo ,
850 radius.StopThickness[iteration-1] - radius.depth );
851 }
852
853 /* this is case where stopping optical depth was specified */
854 if( StopCalc.iptnu != rfield.nupper )
855 {
856 /* end optical depth has been specified */
857 double dt = SDIV(opac.opacity_abs[StopCalc.iptnu-1]*geometry.FillFac);
858 DepthToGo = MIN2(DepthToGo ,
859 (StopCalc.tauend-opac.TauAbsGeo[0][StopCalc.iptnu-1])/dt);
860 }
861
862 if( DepthToGo <= 0. )
864
865 /* set dr if one of above tests have triggered */
866 if( DepthToGo < BigRadius )
867 {
868 double depth_min = MIN2( DepthToGo , thickness_total/100. );
869 DepthToGo = MAX2( DepthToGo / 10. , depth_min );
870 /* want to
871
872 approach the outer edge slowly,
873 * the need for this logic is most evident in brems.in -
874 * HI fraction varies across coronal model */
875 double drThickness = MIN2( thickness_total/10. , DepthToGo );
876 drChoice.insert( pair<const double,string>( drThickness, "depth to go" ) );
877 }
878
879 /* stop AV - usually this is dust, but we consider all opacity sources,
880 * so always include this */
881 /* compute some average grain properties */
882 double AV_to_go = BigRadius;
883 if( rfield.opac_mag_V_extended > SMALLFLOAT && rfield.opac_mag_V_point > SMALLFLOAT )
884 {
885 double SAFETY = 1. + 10.*DBL_EPSILON;
886 /* by default stop av is very large, and opacity can be very small, so ratio
887 * goes to inf - work with logs to see how big the number is */
888 /* apply safety margin to avoid updates to the total extinction getting lost
889 * in the numerical precision */
890 double ave = log10(SAFETY*StopCalc.AV_extended - rfield.extin_mag_V_extended) -
891 log10(rfield.opac_mag_V_extended);
892 double avp = log10(SAFETY*StopCalc.AV_point - rfield.extin_mag_V_point) -
893 log10(rfield.opac_mag_V_point);
894 AV_to_go = MIN2( ave , avp );
895 if( AV_to_go < 37. )
896 {
897 AV_to_go = pow(10., AV_to_go)/geometry.FillFac;
898 /* this is to make sure that we go slightly deeper than AV so that
899 * we trigger this stop */
900 AV_to_go *= 1.0001;
901 }
902 else
903 AV_to_go = BigRadius;
904 /*fprintf(ioQQQ,"DEBUG next dr %.3e %.3e %.3e\n", AV_to_go , ave, avp );*/
905 }
906 drChoice.insert( pair<const double,string>( AV_to_go, "A_V to go" ) );
907
908 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
909 /* spherical models, do not want delta R/R big */
910 double drSphere = radius.Radius*0.04;
911 drChoice.insert( pair<const double,string>( drSphere, "sphericity" ) );
912
913 /* optical depth to electron scattering */
914 double dRTaue = radius.drChange/(dense.eden*6.65e-25*geometry.FillFac);
915 /* allow larger dr when constant temperature,
916 * to prevent some ct models from taking too many cells */
917 if( thermal.lgTemperatureConstant )
918 dRTaue *= 3.;
919 drChoice.insert( pair<const double,string>( dRTaue, "optical depth to electron scattering" ) );
920
921 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
922 if( dense.flong != 0. )
923 {
924 double drFluc = 0.628/2./dense.flong;
925 /* >>chng 04 sep 18, caused cautions that ionization jumps occurred.
926 * set to have the value */
927 drFluc /= 2.;
928 drChoice.insert( pair<const double,string>( drFluc, "density fluctuations" ) );
929 }
930
931 /* keep dr constant in first two zones, if it wants to increase,
932 * but always allow it to decrease.
933 * to guard against large increases in efrac for balmer cont photo dominated models,
934 */
935 double drStart = ( nzone <= 1 ) ? radius.drad : BigRadius;
936 drChoice.insert( pair<const double,string>( drStart, "capped to old DR in first zone" ) );
937
938 // lgSdrmaxRel true if sdrmax is relative to current radius, false if limit in cm
939 double rfacmax = radius.lgSdrmaxRel ? radius.Radius : 1.;
940 drChoice.insert( pair<const double,string>( rfacmax*radius.sdrmax, "sdrmax" ) );
941
942 /* >>chng 05 aug 05, in case of thermal front, where temp is falling quickly and
943 * conditions change very fast, the zone thickness does not really affect the change
944 * in conditions and can cause zoning to become very very thin, which causes
945 * an abort. this occurs between 200 and 1000K. if we are doing temp soln,
946 * temp is between these values, and temp is changing rapidly, do not make sone
947 * thickness much smaller */
948 /* do not use thermal front logic in case of recombination time dep cloud */
949 if( nzone >= 5 && !dynamics.lgTimeDependentStatic )
950 {
951 /* >>chng 05 aug 23, upper bound of thermal from from 1000K to 4000K */
952 /*if( phycon.te > 200. && phycon.te < 1000. && */
953 if( phycon.te > 200. && phycon.te < 3000. &&
954 /* >>chng 05 aug 23, from > 10% in zone to to two zones > 5%,
955 * to fix leiden v3 with large H2 */
956 (struc.testr[nzone-3] - phycon.te)/phycon.te > 0.02 &&
957 (struc.testr[nzone-4] - phycon.te)/phycon.te > 0.02 &&
958 (struc.testr[nzone-5] - phycon.te)/phycon.te > 0.02 )
959 {
960 /* the 0.91 is to make dr unique, so that print statement that
961 * follows will identify this reason */
962 double drThermalFront = radius.drad * 0.91;
963 drChoice.clear();
964 drChoice.insert( pair<const double,string>( drThermalFront, "thermal front logic" ) );
965 }
966 }
967
968 ASSERT( drChoice.size() > 0 );
969
970 /* dr = zero is a logical mistake */
971 if( drChoice.begin()->first <= 0. )
972 {
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() )
977 {
978 fprintf( ioQQQ, " %.2e %s\n", ptr->first, ptr->second.c_str() );
979 ++ptr;
980 }
982 }
983
984 /* option to force min drad relative to depth */
985 if( !radius.lgFixed && drChoice.begin()->first < radius.depth*radius.sdrmin_rel_depth )
986 {
987 drChoice.clear();
988 drChoice.insert( pair<const double,string>( radius.depth*radius.sdrmin_rel_depth, "sdrmin_rel_depth" ) );
989 }
990
991 /* option to force min drad */
992 double rfacmin = radius.lgSdrminRel ? radius.Radius : 1.;
993 if( drChoice.begin()->first < rfacmin*radius.sdrmin )
994 {
995 drChoice.clear();
996 drChoice.insert( pair<const double,string>( rfacmin*radius.sdrmin, "sdrmin" ) );
997 }
998
999 /* this is general code that prevents zone thickness drNext from
1000 * becoming too thin, something that can happen for various bad reasons
1001 * HOWEVER we do not want to do this test for certain density laws,
1002 * for which very small zone thicknesses are unavoidable
1003 * the special cases are:
1004 * special density law,
1005 * globule density law,
1006 * carbon +-0 i front
1007 * the fluctuations command
1008 * drMinimum was set in radius_first to either sdrmin (set drmin) or
1009 * some fraction of the initial radius - it is always set
1010 * to something non-trivial.
1011 * sdrmin is set with the "set dr" command - is SMALLFLOAT by default */
1012 if( strcmp(dense.chDenseLaw,"DLW1") != 0 &&
1013 strcmp(dense.chDenseLaw ,"GLOB") != 0 &&
1014 dense.flong == 0. &&
1015 // stopping on depth to go is not a fault
1016 drChoice.begin()->first != DepthToGo &&
1017 // stopping on Av to go is not a fault
1018 drChoice.begin()->first != AV_to_go )
1019 {
1020 /* don't let dr get smaller than drMinimum, if this resets drNext
1021 * then code stops with warning that zones got too thin
1022 * this can happen due to numerical oscillations, which the code
1023 * tries to damp out by making the zones thinner.
1024 * scale by radius.Radius/radius.rinner to stop very spherical
1025 * simulations from false trigger on dr too small */
1026 /* drMinimum is drad * hden, to make it proportional to optical depth.
1027 * This avoids false trigger across thermal fronts. */
1028 if( drChoice.begin()->first*radius.Radius/radius.rinner*dense.gas_phase[ipHYDROGEN] < radius.drMinimum )
1029 {
1030 radius.drNext = radius.drMinimum/dense.gas_phase[ipHYDROGEN];
1031 /* leaving this at true will cause the model to stop with a warning
1032 * that the zone thickness is too small */
1033 radius.lgDrMinUsed = true;
1034 /* set abort handler */
1035 lgAbort = true;
1036 /* must decrement nzone, since we will not complete this zone, and will not have
1037 * valid structure data for it */
1038 --nzone;
1039 fprintf( ioQQQ,
1040 "\n DISASTER PROBLEM radius_next finds dr too small and aborts. "
1041 "This is zone %ld iteration %ld. The thickness was %.2e\n",
1042 nzone,
1043 iteration,
1044 radius.drNext);
1045 fprintf( ioQQQ,
1046 " If this simulation seems reasonable you can override this limit "
1047 "with the command SET DRMIN %.2e\n\n",
1048 radius.drNext*2);
1049 /* set abort flag */
1050 lgAbort = true;
1051 return 1;
1052 }
1053 }
1054
1055 /* factor to allow for slop in floating numbers */
1056 const double Z = 1.0001;
1057
1058 /* following is to make thickness of model exact
1059 * n.b., on last zone, drNext can be NEGATIVE!!
1060 * DEPTH was incremented at start of zone calc in newrad,
1061 * has been outer edge of zone all throughout */
1062 double drOuterRadius = (radius.StopThickness[iteration-1]-radius.depth)*Z;
1063 drChoice.insert( pair<const double,string>( drOuterRadius, "outer radius" ) );
1064
1065 // choose the smallest dR as the next choice
1066 radius.drNext = drChoice.begin()->first;
1067
1068 /* set flag if dr set by maser */
1069 if( drChoice.begin()->second.find( "H maser" ) != string::npos )
1070 {
1071 rt.lgMaserSetDR = true;
1072 }
1073
1074 /* all this is to only save on last iteration
1075 * the save dr command is not really a save command, making this necessary
1076 * lgDRon is set true if "save dr" entered */
1077 /* lgDRPLst was set true if "save" had "last" on it */
1078 bool lgDoPun = ( save.lgDROn && !( save.lgDRPLst && !iterations.lgLastIt ) );
1079 if( (trace.lgTrace && trace.lgDrBug) || lgDoPun )
1080 {
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() );
1083 }
1084
1085 ASSERT( radius.drNext > 0. );
1086
1087 if( trace.lgTrace )
1088 {
1089 fprintf( ioQQQ, " radius_next chooses next drad drNext=%.4e; this drad was %.4e\n",
1090 radius.drNext, radius.drad );
1091 }
1092
1093 return 0;
1094}
1095
1096/*ContRate called by radius_next to find energy of maximum continuum-gas interaction */
1097STATIC void ContRate(double *freqm,
1098 double *opacm)
1099{
1100 long int i,
1101 ipHi,
1102 iplow,
1103 limit;
1104 double FreqH,
1105 Freq_nonIonizing,
1106 Opac_Hion,
1107 Opac_nonIonizing,
1108 Rate_max_Hion,
1109 Rate_max_nonIonizing;
1110
1111 DEBUG_ENTRY( "ContRate()" );
1112
1113 /*
1114 * find maximum continuum interaction rate,
1115 * these should be reset in following logic without exception,
1116 * if they are still zero at the end we have a logical error
1117 */
1118 Rate_max_nonIonizing = 0.;
1119 Freq_nonIonizing = 0.;
1120 Opac_nonIonizing = 0.;
1121
1122 /* this must be reset to val >= 0 */
1123 *opacm = -1.;
1124 *freqm = -1.;
1125
1126 /* do up to carbon photo edge if carbon is turned on */
1127 /* >>>chng 00 apr 07, add test for whether element is turned on */
1128 if( dense.lgElmtOn[ipCARBON] )
1129 {
1130 /* carbon is turned on, use carbon 1 edge */
1131 ipHi = Heavy.ipHeavy[ipCARBON][0] - 1;
1132 }
1133 else
1134 {
1135 /* carbon turned off, use hydrogen Balmer continuum */
1136 ipHi = iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH2s].ipIsoLevNIonCon-1;
1137 }
1138
1139 /* start at plasma frequency */
1140 for( i=rfield.ipPlasma; i < ipHi; i++ )
1141 {
1142 /* this does not have grain opacity since grains totally passive
1143 * at energies smaller than CI edge */
1144 if( rfield.anu[i]*rfield.flux[0][i]/rfield.widflx[i]*(opac.opacity_abs[i] -
1145 gv.dstab[i]*dense.gas_phase[ipHYDROGEN]) > Rate_max_nonIonizing )
1146 {
1147 Rate_max_nonIonizing = rfield.anu[i]*rfield.flux[0][i]/rfield.widflx[i]*
1148 (opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN]);
1149 Freq_nonIonizing = rfield.anu[i];
1150 Opac_nonIonizing = opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN];
1151 }
1152 }
1153
1154 /* not every continuum extends beyond C edge-this whole loop can add to zero
1155 * use total opacity here
1156 * test is to put in fir continuum if free free heating is important */
1157 if( CoolHeavy.brems_heat_total/thermal.htot > 0.05 )
1158 {
1159 /* this is index for energy where cloud free free optical depth is unity,
1160 * is zero if no freq are opt thin */
1161 iplow = MAX2(1 , rfield.ipEnergyBremsThin );
1162 }
1163 else
1164 {
1165 /* >>>chng 00 apr 07, from Heavy.ipHeavy[0][5] to ipHi defined above, since
1166 * would crash if element not defined */
1167 iplow = ipHi;
1168 }
1169 /* do not go below plasma frequency */
1170 iplow = MAX2( rfield.ipPlasma , iplow );
1171
1172 /* this energy range from carbon edge to hydrogen edge */
1173 limit = MIN2(rfield.nflux,iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1);
1174 for( i=iplow; i < limit; i++ )
1175 {
1176 if( rfield.anu[i]*rfield.flux[0][i]/rfield.widflx[i]*(opac.opacity_abs[i] -
1177 gv.dstab[i]*dense.gas_phase[ipHYDROGEN]) > Rate_max_nonIonizing )
1178 {
1179 Rate_max_nonIonizing = rfield.anu[i]*rfield.flux[0][i]/rfield.widflx[i]*
1180 (opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN]);
1181 Freq_nonIonizing = rfield.anu[i];
1182 Opac_nonIonizing = opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN];
1183 }
1184 }
1185
1186 /* variables to check continuum interactions over Lyman continuum */
1187 Rate_max_Hion = 0.;
1188 Opac_Hion = 0.;
1189 FreqH = 0.;
1190
1191 /* not every continuum extends beyond 1 Ryd-this whole loop can add to zero */
1192 for( i=iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1; i < rfield.nflux; i++ )
1193 {
1194 if( rfield.anu[i]*rfield.flux[0][i]/rfield.widflx[i]*(opac.opacity_abs[i] -
1195 gv.dstab[i]*dense.gas_phase[ipHYDROGEN]) > Rate_max_Hion )
1196 {
1197 /* Rate_max_Hion = anu(i)*flux(i)/widflx(i)*opac(i) */
1198 Rate_max_Hion = rfield.anu[i]*rfield.flux[0][i]/rfield.widflx[i]*
1199 (opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN]);
1200 FreqH = rfield.anu[i];
1201 Opac_Hion = opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN];
1202 }
1203 }
1204
1205 /* use Lyman continuum if its opacity is larger than non-h ion */
1206 if( Rate_max_nonIonizing < 1e-30 && Opac_Hion > SMALLFLOAT )
1207 {
1208 /* this happens for laser source - use Lyman continuum */
1209 *opacm = Opac_Hion;
1210 *freqm = FreqH;
1211 }
1212 /* >>chng 05 aug 03, add last test on Opac_Hion for case where we go very
1213 * deep and very little radiation is left */
1214 else if( Opac_Hion > Opac_nonIonizing && Rate_max_Hion/SDIV(Rate_max_nonIonizing) > 1e-10 && Opac_Hion > SMALLFLOAT )
1215 {
1216 /* use Lyman continuum */
1217 *opacm = Opac_Hion;
1218 *freqm = FreqH;
1219 }
1220 else
1221 {
1222 /* not much rate in the Lyman continuum, stick with low energy */
1223 *opacm = Opac_nonIonizing;
1224 *freqm = Freq_nonIonizing;
1225 }
1226
1227# if 0
1228 /*>>chng 06 aug 12, do not let zones become optically thick to
1229 * peak of dust emission - one of NA's vastly optically thick dust clouds
1230 * did not conserve total energy - very opticallly thick to ir dust emission
1231 * so ir was absorbed and reemitted many times
1232 * this makes sure the cells remain optically thin to dust emission */
1233 if( gv.lgDustOn() && gv.lgGrainPhysicsOn )
1234 {
1235 double fluxGrainPeak = -1.;
1236 long int ipGrainPeak = -1;
1237 long int ipGrainPeak2 = -1;
1238
1239 i = 0;
1240 while( rfield.anu[i] < 0.0912 && i < (rfield.nflux-2) )
1241 {
1242 /* >>chng 06 aug 23. Only want to do this test for the IR dust continuum, therefore only
1243 * check on optical depth for wavelengths greater than 1 micron */
1244 if( gv.GrainEmission[i]*rfield.anu[i]*opac.opacity_abs[i] > fluxGrainPeak )
1245 {
1246 ipGrainPeak = i;
1247 fluxGrainPeak = gv.GrainEmission[i]*rfield.anu[i]*opac.opacity_abs[i];
1248 }
1249 ++i;
1250 }
1251 ASSERT( fluxGrainPeak>=0. && ipGrainPeak >=0 );
1252
1253 /* >>chng 06 aug 23. We have just found the wavelength and flux at the peak of the IR emission.
1254 * Now we want to find the wavelength, short of the peak, which corresponds to 5% of the
1255 * peak emission. This wavelength will be where the code checks to make sure the zone has
1256 * not become optically thick. Since dust opacity generally decreases with wavelength, this assures that
1257 * the optical depths are small for all wavelengths where the flux is appreciable. Tests show that
1258 * this allows for flux/luminosity conservation to within ~5% for an AV of 1e4 mag and at least 2 iterations*/
1259 i = ipGrainPeak;
1260 /* >>chng 06 aug 23. Only want to do this test for the IR dust continuum, therefore only
1261 * check on opacities for wavelengths shortward of the peak and greater than 1 micron
1262 * this routine can be called with the dust emission totally zero - it is only evaluated
1263 * late to save time. don't do the following tests if peak dust emission is zero */
1264 if( fluxGrainPeak > 0. )
1265 {
1266 while( rfield.anu[i] < 0.0912 && i < (rfield.nflux-2) )
1267 {
1268 /* find wavelength where flux = 5% of peak flux, shortward of the peak */
1269 if( gv.GrainEmission[i]*rfield.anu[i]*opac.opacity_abs[i] > 0.05*fluxGrainPeak )
1270 {
1271 /* This will be the array number and flux at 10% of the peak */
1272 ipGrainPeak2 = i;
1273 }
1274 ++i;
1275 }
1276 ASSERT( ipGrainPeak2>=0 );
1277
1278 /* use this as a limit to the dr if opacity is larger */
1279 if( opac.opacity_abs[ipGrainPeak2] > *opacm )
1280 {
1281 /* scale opacity by factor of 10, which further decreases the zone thickness*/
1282 *opacm = opac.opacity_abs[ipGrainPeak2]*10.;
1283 *freqm = rfield.anu[ipGrainPeak2];
1284 /*fprintf(ioQQQ,"DEBUT opac peak %.2e %.2e \n",
1285 *opacm , *freqm );*/
1286 }
1287 }
1288 }
1289# endif
1290
1291 {
1292 /* following should be set true to print contributors */
1293 enum {DEBUG_LOC=false};
1294 if( DEBUG_LOC )
1295 {
1296 fprintf(ioQQQ,"conratedebug \t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
1297 Rate_max_nonIonizing,Freq_nonIonizing,Opac_nonIonizing,
1298 Rate_max_Hion,FreqH ,Opac_Hion,*freqm,*opacm
1299 );
1300
1301 }
1302 }
1303
1304 /* these were set to -1 at start, and must have been reset in one of the
1305 * two loops. Logic error if still <0. */
1306 /* >>chng 05 aug 03, change logic to -1 on entry and check at least zero
1307 * here - will be zero if NO radiation field exists, perhaps since totally
1308 * attenuated */
1309 ASSERT( *opacm >= 0. && *freqm >= 0. );
1310 return;
1311}
1312
1313/*GrainRateDr called by radius_next to find grain heating rate dr */
1314STATIC void GrainRateDr(double *grfreqm,
1315 double *gropacm)
1316{
1317 long int i,
1318 iplow;
1319 double xMax;
1320
1321 DEBUG_ENTRY( "GrainRateDr()" );
1322
1323 /* in all following changed from anu2 to anu july 25 95
1324 *
1325 * find maximum continuum interaction rate */
1326
1327 /* not every continuum extends beyond C edge-this whole loop can add to zero
1328 * use total opacity here
1329 * test is to put in fir continuum if free free heating is important */
1330 if( CoolHeavy.brems_heat_total/thermal.htot > 0.05 )
1331 {
1332 /* this is pointer to energy where cloud free free optical depth is unity,
1333 * is zero if no freq are opt thin */
1334 iplow = MAX2(1 , rfield.ipEnergyBremsThin );
1335 }
1336 else
1337 {
1338 /* do up to carbon photo edge if carbon is turned on */
1339 /* >>>chng 00 apr 07, add test for whether element is turned on */
1340 if( dense.lgElmtOn[ipCARBON] )
1341 {
1342 /* carbon is turned on, use carbon 1 edge */
1343 iplow = Heavy.ipHeavy[ipCARBON][0];
1344 }
1345 else
1346 {
1347 /* carbon truned off, use hydrogen balmer continuum */
1348 iplow = iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH2s].ipIsoLevNIonCon;
1349 }
1350 }
1351
1352 xMax = -1.;
1353 /* integrate up to H edge */
1354 for( i=iplow-1; i < Heavy.ipHeavy[ipHYDROGEN][0]; i++ )
1355 {
1356 if( rfield.anu[i]*rfield.flux[0][i]/rfield.widflx[i]*opac.opacity_abs[i] > xMax )
1357 {
1358 xMax = rfield.anu[i]*rfield.flux[0][i]/rfield.widflx[i]*
1359 opac.opacity_abs[i];
1360 *grfreqm = rfield.anu[i];
1361 *gropacm = opac.opacity_abs[i];
1362 }
1363 }
1364 /* integrate up to heii edge if he is turned on,
1365 * this logic will not make sense if grains on but he off, which in itself makes no sense*/
1366 if( dense.lgElmtOn[ipHELIUM] )
1367 {
1368 for( i=Heavy.ipHeavy[ipHYDROGEN][0]-1; i < Heavy.ipHeavy[ipHELIUM][1]; i++ )
1369 {
1370 if( rfield.anu[i]*rfield.flux[0][i]/rfield.widflx[i]*opac.opacity_abs[i] > xMax )
1371 {
1372 xMax = rfield.anu[i]*rfield.flux[0][i]/rfield.widflx[i]*
1373 opac.opacity_abs[i];
1374 *grfreqm = rfield.anu[i];
1375 *gropacm = opac.opacity_abs[i];
1376 }
1377 }
1378 }
1379
1380 /* possible that there is NO ionizing radiation, in extreme cases,
1381 * if so then xMax will still be negative */
1382 if( xMax <= 0. )
1383 {
1384 *gropacm = 0.;
1385 *grfreqm = 0.;
1386 }
1387 return;
1388}
t_abund abund
Definition abund.cpp:5
long int nzone
Definition cddefines.cpp:14
FILE * ioQQQ
Definition cddefines.cpp:7
bool lgAbort
Definition cddefines.cpp:10
long int iteration
Definition cddefines.cpp:16
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
const int ipCARBON
Definition cddefines.h:310
#define MIN2
Definition cddefines.h:761
const int LIMELM
Definition cddefines.h:258
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
const int ipHELIUM
Definition cddefines.h:306
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
NORETURN void TotalInsanity(void)
Definition service.cpp:886
sys_float SDIV(sys_float x)
Definition cddefines.h:952
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void fixit(void)
Definition service.cpp:991
double & heat() const
Definition collision.h:194
realnum & damp() const
Definition emission.h:563
double & PopOpc() const
Definition emission.h:603
realnum & Pesc() const
Definition emission.h:523
realnum & TauIn() const
Definition emission.h:423
realnum & opacity() const
Definition emission.h:593
double & pump() const
Definition emission.h:473
realnum & dampXvel() const
Definition emission.h:553
CollisionProxy Coll() const
Definition transition.h:424
bool associated() const
Definition transition.h:50
EmissionList::reference Emis() const
Definition transition.h:408
t_colden colden
Definition colden.cpp:5
#define ipCOL_HTOT
Definition colden.h:12
#define ipCOL_H2s
Definition colden.h:18
#define ipCOL_H2g
Definition colden.h:16
#define ipCOL_Hp
Definition colden.h:26
#define ipCOL_H0
Definition colden.h:22
t_CoolHeavy CoolHeavy
Definition coolheavy.cpp:5
const realnum SMALLFLOAT
Definition cpu.h:191
t_dark_matter dark
t_dense dense
Definition dense.cpp:24
double dense_tabden(double r0, double depth)
double dense_fabden(double radius, double depth)
double dense_parametric_wind(double rad)
t_dynamics dynamics
Definition dynamics.cpp:44
t_elementnames elementnames
t_geometry geometry
Definition geometry.cpp:5
static const double SAFETY
GrainVar gv
Definition grainvar.cpp:5
vector< diatomics * > diatoms
Definition h2.cpp:8
vector< diatomics * >::iterator diatom_iter
Definition h2.h:13
t_Heavy Heavy
Definition heavy.cpp:5
t_hmi hmi
Definition hmi.cpp:5
t_hyperfine hyperfine
Definition hyperfine.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
const int ipH1s
Definition iso.h:27
const int ipH2s
Definition iso.h:28
const int ipH_LIKE
Definition iso.h:62
t_iterations iterations
Definition iterations.cpp:5
const TransitionProxy FndLineHt(long int *level)
t_mole_global mole_global
Definition mole.cpp:6
t_mole_local mole
Definition mole.cpp:7
molezone * findspecieslocal(const char buf[])
t_opac opac
Definition opacity.cpp:5
t_phycon phycon
Definition phycon.cpp:6
t_pressure pressure
Definition pressure.cpp:5
t_radius radius
Definition radius.cpp:5
int radius_next(void)
STATIC void GrainRateDr(double *grfreqm, double *gropacm)
STATIC void ContRate(double *freqm, double *opacm)
t_rfield rfield
Definition rfield.cpp:8
t_rt rt
Definition rt.cpp:5
t_save save
Definition save.cpp:5
t_StopCalc StopCalc
Definition stopcalc.cpp:5
t_struc struc
Definition struc.cpp:6
t_thermal thermal
Definition thermal.cpp:5
t_timesc timesc
Definition timesc.cpp:5
t_trace trace
Definition trace.cpp:5
char * chLineLbl(const TransitionProxy &t)
Wind wind
Definition wind.cpp:5