cloudy trunk
Loading...
Searching...
No Matches
rt_continuum.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/*RT_continuum attenuation of diffuse and beamed continua */
4/*pnegopc save negative opacities on io unit, iff 'set negopc' command was given */
5#include "cddefines.h"
6#include "rfield.h"
7#include "opacity.h"
8#include "dense.h"
9#include "colden.h"
10#include "opacity.h"
11#include "geometry.h"
12#include "trace.h"
13#include "radius.h"
14#include "iso.h"
15#include "hextra.h"
16
17/*cmshft - so Compton scattering shift of spectrum
18 * this code is a placeholder */
19STATIC void cmshft(void)
20{
21 DEBUG_ENTRY( "cmshft()" );
22
23 /* first check whether Compton scattering is in as heat/cool */
24 if( !rfield.lgComptonOn )
25 {
26 return;
27 }
28
29 if( rfield.lgComptonOn )
30 {
31 return;
32 }
33
34 /* do reshuffle */
35 for( long i=0; i < rfield.nflux; i++ )
36 {
37 continue;
38 }
39 return;
40}
41
42
43#if !defined(NDEBUG)
44/*pnegopc save negative opacities on io unit, iff 'set negopc' command was given */
45STATIC void pnegopc(void)
46{
47 FILE *ioFile;
48
49 DEBUG_ENTRY( "pnegopc()" );
50
51 if( opac.lgNegOpacIO )
52 {
53 /* option to save negative opacities */
54 ioFile = open_data( "negopc.txt", "w", AS_LOCAL_ONLY );
55 for( long i=0; i < rfield.nflux; i++ )
56 {
57 fprintf( ioFile, "%10.2e %10.2e \n", rfield.anu[i],
58 opac.opacity_abs[i] );
59 }
60 fclose( ioFile);
61 }
62 return;
63}
64#endif
65
66/*RT_continuum attenuation of diffuse and beamed continua */
67void RT_continuum(void)
68{
69
70 DEBUG_ENTRY( "RT_continuum()" );
71
72 if( trace.lgTrace && trace.lgConBug )
73 {
74 fprintf( ioQQQ, " Energy, flux, OTS:\n" );
75 for( long i=0; i < rfield.nflux; i++ )
76 {
77 fprintf( ioQQQ, "%6ld%10.2e%10.2e%10.2e", i, rfield.anu[i],
78 rfield.flux[0][i] + rfield.outlin[0][i] + rfield.ConInterOut[i],
79 rfield.otscon[i] + rfield.otslin[i] + rfield.outlin_noplot[i]);
80 }
81 fprintf( ioQQQ, "\n" );
82 }
83
84 /* begin sanity check in debug mode */
85# if !defined(NDEBUG)
86 bool lgFlxNeg = false;
87 for( long i=0; i < rfield.nflux; i++ )
88 {
89 if( rfield.flux[0][i] < 0. )
90 {
91 fprintf( ioQQQ, " radius_increment finds negative intensity in flux.\n" );
92 fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld\n",
93 rfield.flux[0][i], rfield.anu[i], i );
94 lgFlxNeg = true;
95 }
96 if( rfield.otscon[i] < 0. )
97 {
98 fprintf( ioQQQ, " radius_increment finds negative intensity in otscon.\n" );
99 fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld\n",
100 rfield.otscon[i], rfield.anu[i], i );
101 lgFlxNeg = true;
102 }
103 if( opac.tmn[i] < 0. )
104 {
105 fprintf( ioQQQ, " radius_increment finds negative tmn.\n" );
106 fprintf( ioQQQ, " value, frequency, pointer=%11.3e%11.3e%6ld %4.4s\n",
107 opac.tmn[i], rfield.anu[i], i, rfield.chLineLabel[i] );
108 lgFlxNeg = true;
109 }
110 if( rfield.otslin[i] < 0. )
111 {
112 fprintf( ioQQQ, " radius_increment finds negative intensity in otslin.\n" );
113 fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld %4.4s\n",
114 rfield.otslin[i], rfield.anu[i], i, rfield.chLineLabel[i] );
115 lgFlxNeg = true;
116 }
117 if( rfield.outlin[0][i] < 0. )
118 {
119 fprintf( ioQQQ, " radius_increment finds negative intensity in outlin.\n" );
120 fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld %4.4s\n",
121 rfield.outlin[0][i], rfield.anu[i], i, rfield.chLineLabel[i] );
122 lgFlxNeg = true;
123 }
124 if( rfield.ConInterOut[i] < 0. )
125 {
126 fprintf( ioQQQ, " radius_increment finds negative intensity in ConInterOut.\n" );
127 fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld %4.4s\n",
128 rfield.ConInterOut[i], rfield.anu[i], i, rfield.chContLabel[i] );
129 lgFlxNeg = true;
130 }
131 if( opac.opacity_abs[i] < 0. )
132 {
133 opac.lgOpacNeg = true;
134 /* this sub will save negative opacities on io unit,
135 * iff 'set negopc' command was given */
136 pnegopc();
137 }
138 }
139 if( lgFlxNeg )
140 {
141 fprintf( ioQQQ, " Insanity has occurred, this is zone%4ld\n",
142 nzone );
143 ShowMe();
145 }
146 /*end sanity check*/
147# endif
148
149 // ratio of inner to outer radii, at this point
150 // radius is the outer radius of this zone
151 double DilutionHere = POW2((radius.Radius - radius.drad*radius.dRadSign)/
152 radius.Radius);
153
154 rfield.EnergyIncidCont = 0.;
155 rfield.EnergyDiffCont = 0.;
156
157 // this loop should not be to <= nflux since highest cell is for
158 // continuum unit integration
159 // scattering opacities are included in energy exchange here in the
160 // sphere case, since photons diffuse out of the closed sphere.
161 // scattering opacities are not included as extinction sources in the
162 // sphere case
163 for( long i=0; i < rfield.nflux; i++ )
164 {
165 double dTau_abs = opac.opacity_abs[i]*radius.drad_x_fillfac;
166 double dTau_sct = opac.opacity_sct[i]*radius.drad_x_fillfac;
167
168 // sum total continuous optical depths
169 opac.TauAbsGeo[0][i] += (realnum)(dTau_abs);
170 opac.TauScatGeo[0][i] += (realnum)(dTau_sct);
171
172 // following only optical depth to illuminated face
173 opac.TauAbsFace[i] += (realnum)(dTau_abs);
174 opac.TauScatFace[i] += (realnum)(dTau_sct);
175
176 // these are total in inward direction, large if spherical
177 opac.TauTotalGeo[0][i] = opac.TauAbsGeo[0][i] + opac.TauScatGeo[0][i];
178
179 // attenuation of flux by optical depths IN THIS ZONE
180 // DirectionalCosin is 1/COS(theta), is usually 1, reset with illuminate command,
181 // option for illumination of slab at an angle
182 opac.ExpZone[i] = sexp(dTau_abs*geometry.DirectionalCosin);
183
184 // e(-tau) in inward direction, up to illuminated face
185 opac.ExpmTau[i] *= (realnum)opac.ExpZone[i];
186
187 // e2(tau) in inward direction, up to illuminated face
188 opac.E2TauAbsFace[i] = (realnum)e2(opac.TauAbsFace[i]);
189 ASSERT( opac.E2TauAbsFace[i] <= 1. && opac.E2TauAbsFace[i] >= 0. );
190
191 // on second and later iterations define outward E2
192 if( iteration > 1 )
193 {
194 // e2 from current position to outer edge of shell
195 realnum tau = MAX2(SMALLFLOAT , opac.TauAbsTotal[i] - opac.TauAbsFace[i] );
196 opac.E2TauAbsOut[i] = (realnum)e2( tau );
197 ASSERT( opac.E2TauAbsOut[i]>=0. && opac.E2TauAbsOut[i]<=1. );
198 }
199
200 // DilutionHere is square of ratio of inner to outer radius
201 double AttenuationDilutionFactor = opac.ExpZone[i]*DilutionHere;
202 ASSERT( AttenuationDilutionFactor <= 1.0 );
203
204 // continuum has three parts
205 rfield.flux_beam_const[i] *= (realnum)AttenuationDilutionFactor;
206 rfield.flux_beam_time[i] *= (realnum)AttenuationDilutionFactor;
207 rfield.flux_isotropic[i] *= (realnum)AttenuationDilutionFactor;
208 rfield.flux[0][i] = rfield.flux_beam_const[i] + rfield.flux_beam_time[i] +
209 rfield.flux_isotropic[i];
210
211 // update SummedCon here since flux changed
212 rfield.SummedCon[i] = rfield.flux[0][i] + rfield.SummedDif[i];
213
214 // outward lines and diffuse continua
215 rfield.outlin[0][i] *= (realnum)AttenuationDilutionFactor;
216 rfield.outlin_noplot[i] *= (realnum)AttenuationDilutionFactor;
217
218 // interactive outward diffuse continuum
219 TestCode();// move this over from rt_diffuse - this preserves
220 // the original order and is incorrect
221 rfield.ConInterOut[i] += rfield.DiffuseEscape[i]*(realnum)radius.dVolOutwrd;
222 rfield.ConInterOut[i] *= (realnum)AttenuationDilutionFactor;
223
224 // this is not the interacting continuum
225 rfield.ConEmitOut[0][i] *= (realnum)AttenuationDilutionFactor;
226 rfield.ConEmitOut[0][i] += rfield.ConEmitLocal[nzone][i]*(realnum)radius.dVolOutwrd*opac.tmn[i]/**/;
227
228 // set occupation numbers, first attenuated incident continuum
229 rfield.OccNumbIncidCont[i] = rfield.flux[0][i]*rfield.convoc[i];
230
231 // local diffuse continua
232 rfield.OccNumbDiffCont[i] = rfield.ConEmitLocal[nzone][i]*rfield.convoc[i];
233
234 // outward diffuse continuum
235 rfield.OccNumbContEmitOut[i] = rfield.ConEmitOut[0][i]*rfield.convoc[i];
236
237 // integrated energy flux, ergs s^-1 cm^-2
238 rfield.EnergyIncidCont += rfield.flux[0][i]*rfield.anu[i];
239 rfield.EnergyDiffCont += (rfield.outlin[0][i] + rfield.outlin_noplot[i] +
240 rfield.ConInterOut[i])* rfield.anu[i];
241 }
242
243 // convert Ryd to erg
244 rfield.EnergyIncidCont *= (realnum)EN1RYD;
245 rfield.EnergyDiffCont *= (realnum)EN1RYD;
246
247 // sanity check, compare total Lyman continuum optical depth
248 // with amount of extinction there
249 // this is amount continuum attenuated to illuminated face,
250 // but only do test if flux positive, not counting scattering opacity,
251 // and correction for spherical dilution not important
252 if( rfield.flux[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1]>SMALLFLOAT &&
253 (rfield.flux[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1]/
254 SDIV(rfield.flux_total_incident[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1]) ) > SMALLFLOAT &&
255 !opac.lgScatON &&
256 radius.depth/radius.Radius < 1e-4 )
257 {
258 // ratio of current to incident continuum, converted to optical depth
259 /* >>chng 99 apr 23, this crashed on alpha due to underflow to zero in argy
260 * to log. defended two ways - above checks that ratio of fluxes is large enough,
261 * and here convert to double.
262 * error found by Peter van Hoof */
263 double tau_effec = -log((double)rfield.flux[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1]/
264 (double)opac.tmn[iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1]/
265 (double)rfield.flux_total_incident[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1]);
266
267 // this is computed absorption optical depth to illuminated face
268 double tau_true = opac.TauAbsFace[iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1]*geometry.DirectionalCosin;
269
270 // first test is relative error, second is to absolute error and comes
271 // in for very small tau, where differences are in the round-off
272 if( fabs( tau_effec - tau_true ) / MAX2(tau_effec , tau_true) > 0.01 &&
273 // for very small inner optical depths, the tmn correction is major,
274 // and this test is not precise
275 fabs(tau_effec-tau_true)>MAX2(0.001,1.-opac.tmn[iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1]) )
276 {
277 // in print below must add extra HI col den since this will later be
278 // incremented in RT_tau_inc
279 fprintf( ioQQQ,
280 " PROBLEM radius_increment Lyman continuum insanity zone %li, effective tau=%g, atomic tau=%g simple tau=%g\n",
281 nzone,
282 tau_effec ,
283 tau_true ,
284 6.327e-18*(dense.xIonDense[ipHYDROGEN][0]*radius.drad_x_fillfac+colden.colden[ipCOL_H0]) );
286 }
287 }
288
289 // do scattering opacity, not included when sphere is set
290 if( opac.lgScatON )
291 {
292 for( long i=0; i < rfield.nflux; i++ )
293 {
294 // Lightman and White equation 11 in small epsilon limit,
295 // >>refer continuum RT Lightman, A.P., & White, T.R. 1988, ApJ, 335, 57 */
296 double AttenuationScatteringFactor = 1./(1. + radius.drad_x_fillfac*opac.opacity_sct[i]);
297 ASSERT( AttenuationScatteringFactor <= 1.0 );
298 rfield.flux_beam_const[i] *= (realnum)AttenuationScatteringFactor;
299 rfield.flux_beam_time[i] *= (realnum)AttenuationScatteringFactor;
300 rfield.flux_isotropic[i] *= (realnum)AttenuationScatteringFactor;
301 rfield.flux[0][i] = rfield.flux_beam_const[i] + rfield.flux_beam_time[i] +
302 rfield.flux_isotropic[i];
303
304 rfield.ConInterOut[i] *= (realnum)AttenuationScatteringFactor;
305 rfield.ConEmitOut[0][i] *= (realnum)AttenuationScatteringFactor;
306 rfield.outlin[0][i] *= (realnum)AttenuationScatteringFactor;
307 rfield.outlin_noplot[i] *= (realnum)AttenuationScatteringFactor;
308 }
309 }
310
311 // this dilution is needed to conserve volume in spherical models. tests such
312 // as parispn.in will fault if this is removed
313 realnum Dilution = (realnum)POW2( radius.rinner / (radius.Radius-radius.drad/2.) );
314
315 // this is a unit of energy that will be passed through the code as a test
316 // that all integrations are carried out. A similar test is set in lineset1
317 // and verified in PrtFinal. The opacity at this cell is zero so only
318 // geometrical dilution will affect the integral
319 // Radius is currently outer edge of zone, so radius-drad/2 is radius
320 // of center of zone
321 rfield.ConEmitLocal[nzone][rfield.nflux] = 1.e-10f * Dilution;
322 rfield.DiffuseEscape[rfield.nflux] = 1.e-10f * Dilution;
323 // must do unit integration somewhere
324 rfield.ConInterOut[rfield.nflux] +=
325 rfield.DiffuseEscape[rfield.nflux]*(realnum)radius.dVolOutwrd;
326
327 // opacity should be zero at this energy so J not changed elsewhere
328 ASSERT( opac.opacity_abs[rfield.nflux] == 0. );
329
330 // placeholder code for Compton scattering
331 cmshft();
332
333 // attenuate neutrons if they are present
334 if( hextra.lgNeutrnHeatOn )
335 {
336 // correct for optical depth effects
337 hextra.totneu *= (realnum)sexp(hextra.CrsSecNeutron*
338 dense.gas_phase[ipHYDROGEN]*radius.drad_x_fillfac*geometry.DirectionalCosin);
339 // correct for spherical effects
340 hextra.totneu *= (realnum)DilutionHere;
341 }
342
343 // following radiation factors are extinguished by 1/r**2ilution, electron
344 // scattering by free and bound electrons
345
346 // do all emergent spectrum from illuminated face if model is NOT spherical
347 if( !geometry.lgSphere )
348 {
349 double Reflec_Diffuse_Cont;
350
351 // emission starting at the the plasma frequency
352 for( long i=rfield.ipPlasma-1; i < rfield.nflux; i++ )
353 {
354 if( opac.TauAbsGeo[0][i] < 30. )
355 {
356 // ConEmitLocal is diffuse emission per unit vol, fill factor
357 // the 1/2 comes from isotropic emission
358 Reflec_Diffuse_Cont = rfield.ConEmitLocal[nzone][i]/2.*
359 radius.drad_x_fillfac * opac.E2TauAbsFace[i]*radius.r1r0sq;
360
361 // ConEmitReflec - reflected diffuse continuum
362 rfield.ConEmitReflec[0][i] += (realnum)(Reflec_Diffuse_Cont);
363
364 // the reflected part of the incident continuum
365 rfield.ConRefIncid[0][i] += (realnum)(rfield.flux[0][i]*opac.opacity_sct[i]*
366 radius.drad_x_fillfac*opac.E2TauAbsFace[i]*radius.r1r0sq);
367
368 // reflected line emission
369 rfield.reflin[0][i] += (realnum)(rfield.outlin[0][i]*opac.opacity_sct[i]*
370 radius.drad_x_fillfac*opac.E2TauAbsFace[i]*radius.r1r0sq);
371 }
372 }
373 }
374 return;
375}
long int nzone
Definition cddefines.cpp:14
FILE * ioQQQ
Definition cddefines.cpp:7
long int iteration
Definition cddefines.cpp:16
#define ASSERT(exp)
Definition cddefines.h:578
sys_float sexp(sys_float x)
Definition service.cpp:914
#define STATIC
Definition cddefines.h:97
void TestCode(void)
Definition service.cpp:972
#define POW2
Definition cddefines.h:929
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
double e2(double t)
Definition service.cpp:299
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 ShowMe(void)
Definition service.cpp:181
t_colden colden
Definition colden.cpp:5
#define ipCOL_H0
Definition colden.h:22
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition cpu.cpp:625
@ AS_LOCAL_ONLY
Definition cpu.h:208
const realnum SMALLFLOAT
Definition cpu.h:191
t_dense dense
Definition dense.cpp:24
t_geometry geometry
Definition geometry.cpp:5
t_hextra hextra
Definition hextra.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
const int ipH_LIKE
Definition iso.h:62
t_opac opac
Definition opacity.cpp:5
UNUSED const double EN1RYD
Definition physconst.h:179
t_radius radius
Definition radius.cpp:5
t_rfield rfield
Definition rfield.cpp:8
void RT_continuum(void)
STATIC void pnegopc(void)
STATIC void cmshft(void)
t_trace trace
Definition trace.cpp:5