cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
radius_increment.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*radius_increment do work associated with geometry increments of this zone, called before RT_tau_inc */
4 /*pnegopc punch negative opacities on io unit, iff 'set negopc' command was given */
5 #include "cddefines.h"
6 #include "physconst.h"
7 #include "iso.h"
8 #include "hydrogenic.h"
9 #include "dynamics.h"
10 #include "rfield.h"
11 #include "hextra.h"
12 #include "colden.h"
13 #include "geometry.h"
14 #include "opacity.h"
15 #include "thermal.h"
16 #include "doppvel.h"
17 #include "dense.h"
18 #include "mole.h"
19 #include "h2.h"
20 #include "timesc.h"
21 #include "hmi.h"
22 #include "lines_service.h"
23 #include "taulines.h"
24 #include "trace.h"
25 #include "wind.h"
26 #include "phycon.h"
27 #include "pressure.h"
28 #include "grainvar.h"
29 #include "molcol.h"
30 #include "conv.h"
31 #include "hyperfine.h"
32 #include "mean.h"
33 #include "struc.h"
34 #include "radius.h"
35 /*cmshft - so Compton scattering shift of spectrum */
36 STATIC void cmshft(void);
37 
38 #if !defined(NDEBUG)
39 /*pnegopc punch negative opacities on io unit, iff 'set negopc' command was
40  * given. This function is only used in debug mode */
41 STATIC void pnegopc(void);
42 #endif
43 
44 void radius_increment(void)
45 {
46 # if !defined(NDEBUG)
47  bool lgFlxNeg;
48 # endif
49 
50  long int nelem,
51  i,
52  j ,
53  ion,
54  nzone_minus_1 ,
55  mol;
56  double
57  avWeight,
58  DilutionHere ,
59  escatc,
60  fmol,
61  opfac,
62  relec,
63  rforce,
64  t;
65 
66  double ajmass,
67  Error,
68  rjeans;
69 
70  realnum dradfac;
71 
72  DEBUG_ENTRY( "radius_increment()" );
73 
74  /* when this sub is called radius is the outer edge of zone */
75 
76  if( lgAbort )
77  {
78  return;
79  }
80 
81  /* following block only set of asserts to check that iso populations are sane */
82 # if !defined(NDEBUG)
83  if( !dynamics.lgAdvection )
84  {
85  long int ipISO;
86  for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
87  {
88  for( nelem=ipISO; nelem<LIMELM; ++nelem )
89  {
90  /* >>chng 05 feb 05, rm div by gas_phase for element, to eden,
91  * see explain for this date below */
92  if( dense.lgElmtOn[nelem] &&
93  dense.xIonDense[nelem][nelem-ipISO]/dense.eden>conv.EdenErrorAllowed / 10. &&
94  !(iso.chTypeAtomUsed[ipISO][nelem][0]=='L' &&
95  iso.xIonSimple[ipISO][nelem]<1e-10) )
96  {
97  /* we will check that ground pops are within this factor of the xIonDense value
98  * these are only converged down to some extent - probably this should be related
99  * to the ionization convergence error */
100  /* >>chng 05 feb 05, for very highly ionized sims, where N or O is He-like,
101  * the error can approach 1%, and depends on the convergence criteria.
102  * the code does not try to converge these quantities in an absolute sense,
103  * only the quantities relative to the electron or hydrogen density. so it is not
104  * safe to do the assert for small values
105  * change is to only do this branch if abundance is large enough to be detected by
106  * the ionization convergence solvers */
107 # define ERR_CHK 1.001
108  double err_chk = ERR_CHK;
109  /* >>chng 05 sep 02, when low-T solver used solution is approximate,
110  * and it must not matter (lot-T solver should not be used if it
111  * does matter - so use more liberal expectation */
112  if( iso.chTypeAtomUsed[ipISO][nelem][0]=='L' )
113  err_chk *= 10.;
114  /* make sure that populations are valid, first check Pop2Ion */
115  if( StatesElem[ipISO][nelem][0].Pop >
116  dense.xIonDense[nelem][nelem-ipISO]/(double)SDIV(dense.xIonDense[nelem][nelem+1-ipISO])*err_chk &&
117  /*>>chng 06 jul 22, only do check if pops are within tolerance of double */
118  iso.pop_ion_ov_neut[ipISO][nelem] > 1e-15 )
119  {
120  fprintf(ioQQQ,
121  " DISASTER radius_increment finds inconsistent populations, Pop2Ion zn %.2f ipISO %li nelem %li %.4e %.4e solver:%s\n",
122  fnzone,
123  ipISO , nelem ,
124  StatesElem[ipISO][nelem][0].Pop*dense.xIonDense[nelem][nelem+1-ipISO] ,
125  dense.xIonDense[nelem][nelem-ipISO] ,
126  iso.chTypeAtomUsed[ipISO][nelem]);
127  fprintf(ioQQQ," level solver %s found pop_ion_ov_neut of %.5e",
128  iso.chTypeAtomUsed[ipISO][nelem] ,
129  iso.pop_ion_ov_neut[ipISO][nelem] );
130  fprintf(ioQQQ," ion_solver found pop_ion_ov_neut of %.5e\n",
131  dense.xIonDense[nelem][nelem+1-ipISO]/SDIV( dense.xIonDense[nelem][nelem-ipISO] ) );
132  fprintf(ioQQQ,
133  "simple %.3e Test shows Pop2Ion (%.6e) > xIonDense[nelem]/[nelem+1] (%.6e) \n",
134  iso.xIonSimple[ipISO][nelem],
135  StatesElem[ipISO][nelem][0].Pop ,
136  dense.xIonDense[nelem][nelem-ipISO]/(double)SDIV(dense.xIonDense[nelem][nelem+1-ipISO]) );
137  }
138 
139  /* >>chng 06 jul 24, add assert that this is not search
140  * phase to confirm that removing other asserts on this
141  * are ok - we cannot be in search phase in this routine
142  * if ok, rm this and all ref to lgSearch, perhaps conv.h header */
143  ASSERT( !conv.lgSearch );
144  ASSERT( StatesElem[ipISO][nelem][0].Pop <=
145  dense.xIonDense[nelem][nelem-ipISO]/
146  (double)SDIV(dense.xIonDense[nelem][nelem+1-ipISO])*err_chk ||
147  /*>>chng 06 jul 22, only do check if pops are within tolerance of double */
148  iso.pop_ion_ov_neut[ipISO][nelem] < 1e-15);
149 
150  /* make sure that populations are valid, first check Pop2Ion */
151  if( StatesElem[ipISO][nelem][0].Pop*
152  dense.xIonDense[nelem][nelem+1-ipISO] >
153  dense.xIonDense[nelem][nelem-ipISO]*err_chk &&
154  /*>>chng 06 jul 22, only do check if pops are within tolerance of double */
155  iso.pop_ion_ov_neut[ipISO][nelem] > 1e-15)
156  {
157  fprintf(ioQQQ," DISASTER PopLo fnz %.2f ipISO %li nelem %li pop_ion_ov_neut %.2e 2pops %.6e %.6e solver:%s\n",
158  fnzone,
159  ipISO , nelem ,
160  iso.pop_ion_ov_neut[ipISO][nelem] ,
161  StatesElem[ipISO][nelem][0].Pop*dense.xIonDense[nelem][nelem+1-ipISO],
162  dense.xIonDense[nelem][nelem-ipISO]*err_chk ,
163  iso.chTypeAtomUsed[ipISO][nelem]);
164  }
165  ASSERT(
166  /*(conv.lgSearch || !conv.nPres2Ioniz) || */
167  (StatesElem[ipISO][nelem][0].Pop*dense.xIonDense[nelem][nelem+1-ipISO]<=
168  dense.xIonDense[nelem][nelem-ipISO]*err_chk) ||
169  /*>>chng 06 jul 22, only do check if pops are within tolerance of double */
170  iso.pop_ion_ov_neut[ipISO][nelem] < 1e-15);
171 # undef ERR_CHK
172  }
173  }
174  }
175  }
176 # endif
177 
178  if( trace.lgTrace )
179  {
180  fprintf( ioQQQ,
181  " radius_increment called; radius=%10.3e rinner=%10.3e DRAD=%10.3e drNext=%10.3e ROUTER=%10.3e DEPTH=%10.3e\n",
184  }
185 
186  /* remember mean and largest errors on electron density */
187  Error = fabs(dense.eden - dense.EdenTrue)/SDIV(dense.EdenTrue);
188  if( Error > conv.BigEdenError )
189  {
190  conv.BigEdenError = (realnum)Error;
192  }
193  conv.AverEdenError += (realnum)Error;
194 
195  /* remember mean and largest errors between heating and cooling */
196  Error = fabs(thermal.ctot - thermal.htot) / thermal.ctot;
198  conv.AverHeatCoolError += (realnum)Error;
199 
200  /* remember mean and largest pressure errors */
203  conv.AverPressError += (realnum)Error;
204 
205  /* integrate total mass over model */
207 
208  /* check cooling time for this zone, remember longest */
210  thermal.ctot);
211 
212  /* fmol is fraction of hydrogen in form of any molecule or ion */
213  /* if N_H_MOLEC != 8, we probably need to change this line... */
214  ASSERT( N_H_MOLEC == 8);
216 
217  /* remember the largest H2 fraction that occurs in the model */
219 
220  /* largest fraction of atoms in molecules */
221  for( i=0; i<mole.num_comole_calc; ++i )
222  {
223  if( dense.lgElmtOn[COmole[i]->nelem_hevmol] )
224  {
226  frac *= (realnum) COmole[i]->nElem[COmole[i]->nelem_hevmol];
227  COmole[i]->xMoleFracMax = MAX2(frac,COmole[i]->xMoleFracMax);
228  }
229  }
230 
231  /* H 21 cm equilibrium timescale, H21cm returns H (not e) collisional
232  * deexcitation rate (not cs) */
234  /* >>chng 02 feb 14, add electron term as per discussion in */
235  /* >>refer H1 21cm Liszt, H., 2001, A&A, 371, 698 */
237 
238  /* only update time scale if t is significant */
239  if( t > SMALLFLOAT )
240  timesc.TimeH21cm = MAX2( 1./t, timesc.TimeH21cm );
241 
242  /* remember longest CO timescale */
244  {
245  double timemin;
246  double product = SDIV(dense.xIonDense[ipCARBON][0]*dense.xIonDense[ipOXYGEN][0]);
247  struct molecule *co_molecule;
248  co_molecule = findspecies("CO");
249  timemin = MIN2(timesc.AgeCOMoleDest[co_molecule->index] ,
250  timesc.AgeCOMoleDest[co_molecule->index] * co_molecule->hevmol/ product );
251  /* this is rate CO is destroyed, equal to formation rate in equilibrium */
253  }
254 
255  /* remember longest H2 destruction timescale timescale */
257 
258  /* remember longest H2 formation timescale timescale */
260 
261  /* increment counter if this zone possibly thermally unstable
262  * this flag was set in conv_temp_eden_ioniz.cpp,
263  * derivative of heating and cooling negative */
264  if( thermal.lgUnstable )
265  thermal.nUnstable += 1;
266 
267  /* remember Stromgren radius - where hydrogen ionization falls below half */
269  {
271  rfield.lgUSphON = true;
272  }
273 
274  /* ratio of inner to outer radii, at this point
275  * radius is the outer radius of this zone */
276  DilutionHere = POW2((radius.Radius - radius.drad*radius.dRadSign)/
277  radius.Radius);
278 
279  if( trace.lgTrace && trace.lgConBug )
280  {
281  fprintf( ioQQQ, " Energy, flux, OTS:\n" );
282  for( i=0; i < rfield.nflux; i++ )
283  {
284  fprintf( ioQQQ, "%6ld%10.2e%10.2e%10.2e", i, rfield.anu[i],
285  rfield.flux[i] + rfield.outlin[i] + rfield.ConInterOut[i],
287  }
288  fprintf( ioQQQ, "\n" );
289  }
290 
291  /* diffuse continua factor
292  * if lgSphere then all diffuse radiation is outward (COVRT=1)
293  * lgSphere not set then COVRT=0.0 */
294 
295  /* begin sanity check */
296  /* only do this test in debug mode */
297 # if !defined(NDEBUG)
298  lgFlxNeg = false;
299  for( i=0; i < rfield.nflux; i++ )
300  {
301  if( rfield.flux[i] < 0. )
302  {
303  fprintf( ioQQQ, " radius_increment finds negative intensity in flux.\n" );
304  fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld\n",
305  rfield.flux[i], rfield.anu[i], i );
306  lgFlxNeg = true;
307  }
308  if( rfield.otscon[i] < 0. )
309  {
310  fprintf( ioQQQ, " radius_increment finds negative intensity in otscon.\n" );
311  fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld\n",
312  rfield.otscon[i], rfield.anu[i], i );
313  lgFlxNeg = true;
314  }
315  if( opac.tmn[i] < 0. )
316  {
317  fprintf( ioQQQ, " radius_increment finds negative tmn.\n" );
318  fprintf( ioQQQ, " value, frequency, pointer=%11.3e%11.3e%6ld %4.4s\n",
319  opac.tmn[i], rfield.anu[i], i, rfield.chLineLabel[i] );
320  lgFlxNeg = true;
321  }
322  if( rfield.otslin[i] < 0. )
323  {
324  fprintf( ioQQQ, " radius_increment finds negative intensity in otslin.\n" );
325  fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld %4.4s\n",
326  rfield.otslin[i], rfield.anu[i], i, rfield.chLineLabel[i] );
327  lgFlxNeg = true;
328  }
329  if( rfield.outlin[i] < 0. )
330  {
331  fprintf( ioQQQ, " radius_increment finds negative intensity in outlin.\n" );
332  fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld %4.4s\n",
333  rfield.outlin[i], rfield.anu[i], i, rfield.chLineLabel[i] );
334  lgFlxNeg = true;
335  }
336  if( rfield.ConInterOut[i] < 0. )
337  {
338  fprintf( ioQQQ, " radius_increment finds negative intensity in ConInterOut.\n" );
339  fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld %4.4s\n",
340  rfield.ConInterOut[i], rfield.anu[i], i, rfield.chContLabel[i] );
341  lgFlxNeg = true;
342  }
343  if( opac.opacity_abs[i] < 0. )
344  {
345  opac.lgOpacNeg = true;
346  /* this sub will punch negative opacities on io unit,
347  * iff 'set negopc' command was given */
348  pnegopc();
349  }
350  }
351  if( lgFlxNeg )
352  {
353  fprintf( ioQQQ, " Insanity has occurred, this is zone%4ld\n",
354  nzone );
355  ShowMe();
356  cdEXIT(EXIT_FAILURE);
357  }
358  /*end sanity check*/
359 # endif
360 
361  /* rfield.lgOpacityFine flag set false with no fine opacities command
362  * tests show that always evaluating this changes fast run of
363  * pn_paris from 26.7 sec to 35.1 sec
364  * but must always update fine opacities since used for transmission */
365  if( rfield.lgOpacityFine )
366  {
367  dradfac = (realnum)radius.drad_x_fillfac;
368  /* increment the fine opacity array */
369  for( i=0; i<rfield.nfine; ++i )
370  {
371  rfield.fine_opt_depth[i] +=
372  rfield.fine_opac_zone[i]*dradfac;
373  }
374 
375  /* evaluate fine opacity transmission coefficient for transmission
376  * through this zone. This is assumed to be a scattering opacity,
377  * only do this if including scattering */
378  if( opac.lgScatON )
379  {
380  /* sum over coarse continuum */
381  for( i=0; i < rfield.nflux-1; i++ )
382  {
383  /* find transmission coefficient if lower and upper bounds of coarse continuum is within
384  * boundaries of fine continuum
385  * unity is default */
387  {
388  rfield.trans_coef_zone[i] = 0.;
389  for( j=rfield.ipnt_coarse_2_fine[i]; j<rfield.ipnt_coarse_2_fine[i+1]; ++j )
390  {
391  /* find transmission coefficient for this zone */
392  rfield.trans_coef_zone[i] += 1.f / (1.f + rfield.fine_opac_zone[j]*dradfac );
393  }
394  /* first branch is normal case, where fine continuum is finer than
395  * coarse continuum. But, when end temp is very high, fine continuum is
396  * very coarse, so may be just one cell, and following will not pass */
398  {
400  }
401  else
402  {
403  /* in case where fine is coarser than coarse, just use firs cell */
404  rfield.trans_coef_zone[i] =
405  1.f / (1.f + rfield.fine_opac_zone[rfield.ipnt_coarse_2_fine[i]]*dradfac );
406  }
407  }
408  else
409  {
410  /* this branch, not within find continuum energy range - transmission unity */
411  rfield.trans_coef_zone[i] = 1.f;
412  }
413 
414  /* the total is just the current total times the transmission coefficient in this zone. */
416  }
417  }
418  }
419 
420  /* electron scattering opacity */
421  escatc = 6.65e-25*dense.eden;
422 
423  /* this loop should not be to <= nflux since we not want to clobber the
424  * continuum unit integration */
425  relec = 0.;
426  rforce = 0.;
427  for( i=0; i < rfield.nflux; i++ )
428  {
429  /* sum total continuous optical depths */
432 
433  /* following only optical depth to illuminated face */
436 
437  /* these are total in inward direction, large if spherical */
438  opac.TauTotalGeo[0][i] = opac.TauAbsGeo[0][i] + opac.TauScatGeo[0][i];
439 
440  /* TMN is array of scale factors which account for attenuation
441  * of continuum across the zone (necessary to conserve energy
442  * at the 1% - 5% level.) sphere effects in
443  * drNext was set by NEXTDR and will be next dr */
444  /* compute both total and Thomson scat rad acceleration */
445  rforce += (rfield.flux[i] + rfield.outlin[i] + rfield.outlin_noplot[i]+
447  opac.opacity_sct[i]);
448 
449  relec += ((rfield.flux[i] + rfield.outlin[i] + rfield.outlin_noplot[i]+
450  rfield.ConInterOut[i])* escatc)*rfield.anu[i];
451 
452  /* attenuation of flux by optical depths IN THIS ZONE
453  * AngleIllumRadian is 1/COS(theta), is usually 1, reset with illuminate command,
454  * option for illumination of slab at an angle */
457 
458  /* e(-tau) in inward direction, up to illuminated face */
459  opac.ExpmTau[i] *= (realnum)opac.ExpZone[i];
460 
461  /* e2(tau) in inward direction, up to illuminated face */
463  ASSERT( opac.E2TauAbsFace[i] <= 1. && opac.E2TauAbsFace[i] >= 0. );
464 
465  /* on second and later iterations define outward E2 */
466  if( iteration > 1 )
467  {
468  /* e2 from current position to outer edge of shell */
470  opac.E2TauAbsOut[i] = (realnum)e2( tau );
471  ASSERT( opac.E2TauAbsOut[i]>=0. && opac.E2TauAbsOut[i]<=1. );
472  }
473 
474  /* DilutionHere is square of ratio of inner to outer radius */
475  opfac = opac.ExpZone[i]*DilutionHere;
476  ASSERT( opfac <= 1.0 );
477 
478  /* >>chng 07 may 22, break continuum into three parts */
479  rfield.flux_beam_const[i] *= (realnum)opfac;
480  rfield.flux_beam_time[i] *= (realnum)opfac;
481  rfield.flux_isotropic[i] *= (realnum)opfac;
484 
485  /* update SummedCon here since flux changed */
487 
488  /* outward lines and continua */
489  rfield.ConInterOut[i] *= (realnum)opfac;
490  rfield.outlin[i] *= (realnum)opfac;
491  rfield.outlin_noplot[i] *= (realnum)opfac;
492 
493  rfield.ConEmitOut[i] *= (realnum)opfac;
495 
496  /* set occupation numbers, first attenuated incident continuum */
498 
499  /* >>chng 00 oct 03, add diffuse continua */
500  /* local diffuse continua */
502 
503  /* >>chng 05 feb 16, occupation number of outward diffuse continuum */
505  }
506 
507  /* begin sanity check, compare total Lyman continuum optical depth
508  * with amount of extinction there */
509 
510  /* this is amount continuum attenuated to illuminated face,
511  * but only do test if flux positive, not counting scattering opacity,
512  * and correction for spherical dilution not important */
513  /* >>chng 02 dec 05, add test for small float, protect against models
514  * where we have gone below smallfloat, and so float is ragged */
518  !opac.lgScatON &&
519  radius.depth/radius.Radius < 1e-4 )
520  {
521  /* ratio of current to incident continuum, converted to optical depth */
522  /* >>chng 99 apr 23, this crashed on alpha due to underflow to zero in argy
523  * to log. defended two ways - above checks that ratio of fluxes is large enough,
524  * and here convert to double.
525  * error found by Peter van Hoof */
526  double tau_effec = -log((double)rfield.flux[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]/
527  (double)opac.tmn[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]/
529 
530  /* this is computed absorption optical depth to illuminated face */
531  /* >>chng 06 jul 11, add geometry factor for illumination at an angle - bug fix - should have
532  * always been there */
534 
535  /* first test is relative error, second is to absolute error and comes
536  * in for very small tau, where differences are in the round-off */
537  if( fabs( tau_effec - tau_true ) / t > 0.01 &&
538  /* for very small inner optical depths, the tmn correction is major,
539  * and this test is not precise */
540  fabs(tau_effec-tau_true)>MAX2(0.001,1.-opac.tmn[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]) )
541  {
542  /* in print below must add extra HI col den since this will later be
543  * incremented in RT_tau_inc */
544  fprintf( ioQQQ,
545  " PROBLEM radius_increment Lyman continuum insanity zone %li, effective tau=%g, atomic tau=%g simple tau=%g\n",
546  nzone,
547  tau_effec ,
548  tau_true ,
550  TotalInsanity();
551  }
552  }
553  /* end sanity check */
554 
555  /* do scattering opacity, intensity goes as 1/(1+tau)
556  * scattering opacity is turned off when sphere is set */
557  if( opac.lgScatON )
558  {
559  for( i=0; i < rfield.nflux; i++ )
560  {
561  /* assume half forward scattering
562  * opfac = 1. / ( 1. + dReff*0.5 * scatop(i) )
563  * >>chng 97 apr 25, remove .5 to get agreement with
564  * Lightman and White equation 11 in small epsilon limit,
565  * >>refer continuum RT Lightman, A.P., & White, T.R. 1988, ApJ, 335, 57 */
566  opfac = 1./(1. + radius.drad_x_fillfac*opac.opacity_sct[i]);
567  ASSERT( opfac <= 1.0 );
568  rfield.flux_beam_const[i] *= (realnum)opfac;
569  rfield.flux_beam_time[i] *= (realnum)opfac;
570  rfield.flux_isotropic[i] *= (realnum)opfac;
573 
574  rfield.ConInterOut[i] *= (realnum)opfac;
575  rfield.ConEmitOut[i] *= (realnum)opfac;
576  rfield.outlin[i] *= (realnum)opfac;
577  rfield.outlin_noplot[i] *= (realnum)opfac;
578  }
579  }
580 
581  /* now do slight reshuffling of energy due to Compton scattering */
582  cmshft();
583 
584  /* remember the largest value */
586 
587  /* force multiplier; relec can be zero for very low densities - so use double
588  * form of safe_div - fmul is of order unity - wind.AccelLine and wind.AccelCont
589  * are both floats to will underflow long before relec will - fmul is only used
590  * in output, not in any physics */
591  relec *= (EN1RYD/SPEEDLIGHT/dense.xMassDensity);
592  if( relec > SMALLFLOAT )
593  wind.fmul = (realnum)( (wind.AccelLine + wind.AccelCont) / relec);
594  else
595  wind.fmul = 0.;
596 
597  /* keep track of average acceleration */
600 
601  /* following is integral of radiative force */
603  /*fprintf(ioQQQ," debuggg pinzon %.2f %.2e %.2e %.2e\n",
604  fnzone,pressure.pinzon,dense.xMassDensity,wind.AccelTot);*/
606 
607  /* sound is sound travel time, sqrt term is sound speed */
609  /* adiabatic sound speed assuming mono-atomic gas - gamma is 5/3*/
612 
613  /* attenuate neutrons if they are present */
614  if( hextra.lgNeutrnHeatOn )
615  {
616  /* correct for optical depth effects */
618  /* correct for spherical effects */
619  hextra.totneu *= (realnum)DilutionHere;
620  }
621 
622  /* following radiation factors are extinguished by 1/r**2 and e- opacity
623  * also bound electrons */
624 
625  /* do all emergent spectrum from illuminated face if model is NOT spherical */
626  if( !geometry.lgSphere )
627  {
628  double Reflec_Diffuse_Cont;
629 
630  /* >>chng 01 jul 14, from lower limit of 0 to plasma frequency -
631  * never should have added diffuse emission from below the plasma frequency */
632  for( i=rfield.ipPlasma-1; i < rfield.nflux; i++ )
633  {
634  if( opac.TauAbsGeo[0][i] < 30. )
635  {
636  /* ConEmitLocal is diffuse emission per unit vol, fill factor
637  * the 1/2 comes from isotropic emission */
638  Reflec_Diffuse_Cont = rfield.ConEmitLocal[nzone][i]/2.*
640 
641  /* ConEmitReflec - reflected diffuse continuum */
642  rfield.ConEmitReflec[i] += (realnum)(Reflec_Diffuse_Cont);
643 
644  /* the reflected part of the incident continuum */
647 
648  /* reflected line emission */
651  }
652  }
653  }
654 
655  /* following is general method to find means weighted by various functions
656  * called in IterStart to initialize to zero, called here to put in numbers
657  * results will be weighted by radius and volume
658  * this is the only place things must be entered to create averages
659  * code is in mean.c */
660  aver("zone",1.,1.," ");
661  aver("doit",phycon.te,1.," Te ");
662  aver("doit",phycon.te,dense.eden," Te(Ne) ");
663  aver("doit",phycon.te,dense.eden*dense.xIonDense[ipHYDROGEN][1]," Te(NeNp) ");
664  aver("doit",phycon.te,dense.eden*dense.xIonDense[ipHELIUM][1] ," Te(NeHe+)");
665  aver("doit",phycon.te,dense.eden*dense.xIonDense[ipHELIUM][2] ,"Te(NeHe2+)");
666  aver("doit",phycon.te,dense.eden*dense.xIonDense[ipOXYGEN][1] ," Te(NeO+) " );
667  aver("doit",phycon.te,dense.eden*dense.xIonDense[ipOXYGEN][2] ," Te(NeO2+)");
668  /*aver("doit",phycon.te,hmi.Hmolec[ipMH2g]," Te(H2) ");*/
669  aver("doit",phycon.te,hmi.H2_total," Te(H2) ");
670  aver("doit",dense.gas_phase[ipHYDROGEN],1.," N(H) ");
671  aver("doit",dense.eden,dense.xIonDense[ipOXYGEN][2]," Ne(O2+) ");
672  aver("doit",dense.eden,dense.xIonDense[ipHYDROGEN][1]," Ne(Np) ");
673 
674  /* save information about structure of model, now used to get t^2 */
675  /* max because if program aborts during search phase, will get to here
676  * with nzone = -1 */
677  nzone_minus_1 = MAX2( 0, nzone-1 );
678  /* this is number of struc.xx zones with valid data */
679  struc.nzone = nzone_minus_1+1;
680  ASSERT(nzone_minus_1>=0 && nzone_minus_1 < struc.nzlim );
681  struc.testr[nzone_minus_1] = (realnum)phycon.te;
682  /* number of particles per unit vol */
683  struc.DenParticles[nzone_minus_1] = dense.pden;
684  /* >>chng 02 May 2001 rjrw: add hden for dilution */
685  struc.hden[nzone_minus_1] = (realnum)dense.gas_phase[ipHYDROGEN];
686  /* total grams per unit vol */
687  struc.DenMass[nzone_minus_1] = dense.xMassDensity;
688  struc.heatstr[nzone_minus_1] = thermal.htot;
689  struc.coolstr[nzone_minus_1] = thermal.ctot;
690  struc.volstr[nzone_minus_1] = (realnum)radius.dVeff;
691  struc.drad[nzone_minus_1] = (realnum)radius.drad;
693  struc.histr[nzone_minus_1] = dense.xIonDense[ipHYDROGEN][0];
694  struc.hiistr[nzone_minus_1] = dense.xIonDense[ipHYDROGEN][1];
695  struc.ednstr[nzone_minus_1] = (realnum)dense.eden;
696  struc.o3str[nzone_minus_1] = dense.xIonDense[ipOXYGEN][2];
697  struc.pressure[nzone_minus_1] = (realnum)pressure.PresTotlCurr;
698  struc.windv[nzone_minus_1] = (realnum)wind.windv;
699  struc.AccelTot[nzone_minus_1] = wind.AccelTot;
700  struc.AccelGravity[nzone_minus_1] = wind.AccelGravity;
702  struc.GasPressure[nzone_minus_1] = (realnum)pressure.PresGasCurr;
703  struc.depth[nzone_minus_1] = (realnum)radius.depth;
704  /* save absorption optical depth from illuminated face to current position */
706  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
707  {
708  struc.gas_phase[nelem][nzone_minus_1] = dense.gas_phase[nelem];
709  for( ion=0; ion<nelem+2; ++ion )
710  {
711  struc.xIonDense[nelem][ion][nzone_minus_1] = dense.xIonDense[nelem][ion];
712  }
713  }
714 
715  /* the hydrogen molecules */
716  for(mol=0;mol<N_H_MOLEC;mol++)
717  {
718  struc.H2_molec[mol][nzone_minus_1] = hmi.Hmolec[mol];
719  }
720 
721  /* the heavy element molecules */
722  for(mol=0;mol<mole.num_comole_calc;mol++)
723  {
724  struc.CO_molec[mol][nzone_minus_1] = COmole[mol]->hevmol;
725  }
726 
731 
732  /* >>chng 05 jan 09, add integral of n(H0) / Tspin */
734 
735  /* >>chng 05 Mar 07, add integral of n(OH) / Tspin */
737 
738  /* this is Lya excitation temperature */
740 
741  /* count number of times Lya excitation temp hotter than gas */
742  if( hydro.TexcLya > phycon.te )
743  {
744  hydro.nLyaHot += 1;
745  if( hydro.TexcLya > hydro.TLyaMax )
746  {
749  hydro.nZTLaMax = nzone;
750  }
751  }
752 
753  /* column densities in various species */
757  /* >>chng 02 sep 20, from htwo to H2_total */
758  /* >>chng 05 mar 14, rather than H2_total, give H2g and H2s */
761  /* this is a special form of column density - should be proportional to total shielding */
768 
770  /* He I t2S column density*/
773 
774  /* the ortho and para column densities */
777  ASSERT( fabs( h2.ortho_density + h2.para_density - hmi.H2_total ) / SDIV( hmi.H2_total ) < 1e-4 );
778  /*fprintf(ioQQQ,"DEBUG ortho para\t%.3e\t%.3e\ttot\t%.3e\t or pa colden\t%.3e\t%.3e\n",
779  h2.ortho_density, h2.para_density,hmi.H2_total,
780  h2.ortho_colden , h2.para_colden);*/
781 
782  /*>>chng 27mar, GS, Column density of F=0 and F=1 levels of H0*/
785  /*fprintf(ioQQQ,"DEBUG pophi-poplo\t%.3e\t%.3e\radius\t%.3e\t col_hi\t%.3e\t%.3e\n",
786  HFLines[0].PopHi, HFLines[0].PopLo, radius.drad_x_fillfac,
787  HFLines[0].PopHi*radius.drad_x_fillfac,colden.H0_21cm_upper );*/
788 
789  /* the CII and SiII atoms are resolved */
790  for( i=0; i < 5; i++ )
791  {
792  /* pops and column density for C II atom */
794  /* pops and column density for SiII atom */
796  }
797  for( i=0; i < 3; i++ )
798  {
799  /* pops and column density for CI atom */
801  /* pops and column density for OI atom */
803  }
804  for( i=0; i < 4; i++ )
805  {
806  /* pops and column density for CIII atom */
808  }
809 
810  /* now add total molecular column densities */
811  molcol("ADD ",ioQQQ);
812 
813  /* increment forming the mean ionization and temperature */
814  MeanInc();
815 
816  /*-----------------------------------------------------------------------*/
817 
818  /* calculate average atomic weight per hydrogen of the plasma */
819  avWeight = 0.;
820  for( nelem=0; nelem < LIMELM; nelem++ )
821  {
822  avWeight += dense.gas_phase[nelem]*dense.AtomicWeight[nelem];
823  }
824  avWeight /= dense.gas_phase[ipHYDROGEN];
825 
826  /* compute some average grain properties */
831  long int nd;
832  for( nd=0; nd < gv.nBin; nd++ )
833  {
834 
835  /* this is total extinction in magnitudes at V and B, for a point source
836  * total absorption and scattering,
837  * does not discount forward scattering to be similar to stellar extinction
838  * measurements made within ism */
840  gv.bin[nd]->pure_sc1[rfield.ipB_filter-1])*gv.bin[nd]->dstAbund*
843  gv.bin[nd]->pure_sc1[rfield.ipB_filter-1])*gv.bin[nd]->dstAbund*
845 
847  gv.bin[nd]->pure_sc1[rfield.ipV_filter-1])*gv.bin[nd]->dstAbund*
849 
851  gv.bin[nd]->pure_sc1[rfield.ipV_filter-1])*gv.bin[nd]->dstAbund*
853 
854  /* this is total extinction in magnitudes at V and B, for an extended source
855  * total absorption and scattering,
856  * DOES discount forward scattering to apply for extended source like Orion */
863 
870 
874  gv.bin[nd]->avDGRatio += (realnum)(gv.bin[nd]->dustp[1]*gv.bin[nd]->dustp[2]*
875  gv.bin[nd]->dustp[3]*gv.bin[nd]->dustp[4]*gv.bin[nd]->dstAbund/avWeight*
877  }
878 
879  /* there are some quantities needed to calculation the Jeans mass and radius */
883 
884  /* now find minimum Jeans length and mass; length in cm */
885  rjeans = 7.79637 + (phycon.alogte - log10(dense.wmole) - log10(dense.xMassDensity*
886  geometry.FillFac))/2.;
887 
888  /* minimum Jeans mass in gm */
889  ajmass = 3.*(rjeans - 0.30103) + log10(4.*PI/3.*dense.xMassDensity*
890  geometry.FillFac);
891 
892  /* now remember smallest */
893  colden.rjnmin = MIN2(colden.rjnmin,(realnum)rjeans);
894  colden.ajmmin = MIN2(colden.ajmmin,(realnum)ajmass);
895 
896  if( trace.lgTrace )
897  {
898  fprintf( ioQQQ, " radius_increment returns\n" );
899  }
900  return;
901 }
902 
903 #if !defined(NDEBUG)
904 /*pnegopc punch negative opacities on io unit, iff 'set negopc' command was given */
905 STATIC void pnegopc(void)
906 {
907  long int i;
908  FILE *ioFile;
909 
910  DEBUG_ENTRY( "pnegopc()" );
911 
912  if( opac.lgNegOpacIO )
913  {
914  /* option to punch negative opacities */
915  ioFile = open_data( "negopc.txt", "w", AS_LOCAL_ONLY );
916  for( i=0; i < rfield.nflux; i++ )
917  {
918  fprintf( ioFile, "%10.2e %10.2e \n", rfield.anu[i],
919  opac.opacity_abs[i] );
920  }
921  fclose( ioFile);
922  }
923  return;
924 }
925 #endif
926 /* Note - when this file is compiled in fast optimized mode,
927  * a good compiler will complain that the routine pnegopc is a
928  * static, local routine, but that it has not been used. It is
929  * indeed used , but only when NDEBUG is not set, so it is only
930  * used when the code is compiled in debug mode. So this is
931  * not a problem. */
932 
933 
934 /*cmshft - so Compton scattering shift of spectrum */
935 STATIC void cmshft(void)
936 {
937  long int i;
938 
939  DEBUG_ENTRY( "cmshft()" );
940 
941  /* first check whether Compton scattering is in as heat/cool */
942  if( rfield.comoff == 0. )
943  {
944  return;
945  }
946 
947  if( rfield.comoff != 0. )
948  {
949  return;
950  }
951 
952  /* do reshuffle */
953  for( i=0; i < rfield.nflux; i++ )
954  {
955  continue;
956  /* watch this space for some really great code!!!!
957  * COMUP needs factor of TE to be Compton cooling */
958  }
959  return;
960 }

Generated for cloudy by doxygen 1.8.4