cloudy trunk
Loading...
Searching...
No Matches
cont_ffun.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/*ffun evaluate total flux for sum of all continuum sources */
4/*ffun1 derive flux at a specific energy, for one continuum */
5/*ReadTable called by TABLE READ to read in continuum from PUNCH TRANSMITTED CONTINUUM */
6#include "cddefines.h"
7#include "physconst.h"
8#include "rfield.h"
9#include "ipoint.h"
10#include "opacity.h"
11#include "continuum.h"
12
13double PlanckFunction( double Temp, double E_Ryd );
14
15/* evaluate sum of all individual continua at one energy, return is
16* continuum intensity */
17double ffun(
18 /* the energy in Rydbergs where the continuum will be evaluated */
19 double anu )
20{
21 double frac_beam_time;
22 /* fraction of beamed continuum that is constant */
23 double frac_beam_const;
24 /* fraction of continuum that is isotropic */
25 double frac_isotropic;
26 double a;
27
28 DEBUG_ENTRY( "ffun()" );
29
30 /* call real function */
31 a = ffun( anu , &frac_beam_time , &frac_beam_const , &frac_isotropic );
32 return a;
33}
34
35/* evaluate sum of all individual continua at one energy, return is
36 * continuum intensity */
37double ffun(
38 /* the energy in Rydbergs where the continuum will be evaluated */
39 double anu ,
40 /* fraction of beamed continuum that is varies with time */
41 double *frac_beam_time,
42 /* fraction of beamed continuum that is constant */
43 double *frac_beam_const,
44 /* fraction of continuum that is isotropic */
45 double *frac_isotropic )
46{
47 double ffun_v;
48 static bool lgWarn = false;
49 double flx_beam_time , flx_beam_const , flx_isotropic;
50
51 DEBUG_ENTRY( "ffun()" );
52
53 /* This routine, ffun, returns the sum of photons per unit time, area, energy,
54 * for all continua in the calculation. ffun1 is called and returns
55 * a single continuum. We loop over all nShape continuum sources -
56 * ipspec points to each */
57 ffun_v = 0.;
58 flx_beam_time = 0.;
59 flx_beam_const = 0.;
60 flx_isotropic = 0.;
61 for( rfield.ipSpec=0; rfield.ipSpec < rfield.nShape; rfield.ipSpec++ )
62 {
63 double one = ffun1(anu)*rfield.spfac[rfield.ipSpec];
64 ffun_v += one;
65
66 /* find fraction of total that is constant vs variable and
67 * isotropic vs beamed */
68 if( rfield.lgBeamed[rfield.ipSpec] )
69 {
70 if( rfield.lgTimeVary[rfield.ipSpec] )
71 flx_beam_time += one;
72 else
73 flx_beam_const += one;
74 }
75 else
76 flx_isotropic += one;
77 }
78
79 /* at this point rfield.flux is the sum of the following three continua
80 * now keep track of the three different types */
81 if( ffun_v < SMALLFLOAT )
82 {
83 *frac_beam_time = 0.;
84 *frac_beam_const = 1.;
85 *frac_isotropic = 0.;
86 }
87 else
88 {
89 /* fraction of beamed continuum that varies with time */
90 *frac_beam_time = flx_beam_time / ffun_v;
91 /* part of beamed continuum that is constant */
92 *frac_beam_const = flx_beam_const / ffun_v;
93 /* the constant isotropic continuum */
94 *frac_isotropic = flx_isotropic / ffun_v;
95 }
96 ASSERT( *frac_beam_time >=0. && *frac_beam_time<=1.+3.*DBL_EPSILON );
97 ASSERT( *frac_beam_const >=0.&& *frac_beam_const<=1.+3.*DBL_EPSILON );
98 ASSERT( *frac_isotropic >=0. && *frac_isotropic<=1.+3.*DBL_EPSILON );
99 ASSERT( fabs( 1.-*frac_beam_time-*frac_beam_const-*frac_isotropic)<
100 10.*DBL_EPSILON);
101
102 if( ffun_v > BIGFLOAT && !lgWarn )
103 {
104 lgWarn = true;
105 fprintf( ioQQQ, " FFUN: The net continuum is very intense.\n" );
106 fprintf( ioQQQ, " I will try to press on, but may have problems.\n" );
107 }
108 return( ffun_v );
109}
110
111/*ffun1 derive flux at a specific energy, for one continuum */
112double ffun1(double xnu)
113{
114 char chKey[6];
115 long int i;
116 double fac,
117 ffun1_v,
118 y;
119 static bool lgWarn = false;
120
121 DEBUG_ENTRY( "ffun1()" );
122
123
124 /* confirm that pointer is within range */
125 ASSERT( rfield.ipSpec >= 0);
126 ASSERT( rfield.ipSpec < rfield.nShape );
127
128 /* FFUN1 returns photons per unit time, area, energy, for one continuum
129 * ipspec is the pointer to the continuum source, in the order
130 * entered on the command lines */
131
132 /*begin sanity check */
133 ASSERT( xnu >= rfield.emm*0.99 );
134 ASSERT( xnu <= rfield.egamry*1.01 );
135 /*end sanity check */
136
137 strcpy( chKey, rfield.chSpType[rfield.ipSpec] );
138
139 if( strcmp(chKey,"AGN ") == 0 )
140 {
141 /* power law with cutoff at both ends
142 * nomenclature all screwed up - slope is cutoff energy in Ryd,
143 * cutoff[1][i] is ratio of two continua from alpha ox
144 * cutoff[2][i] is slope of bb temp */
145 ffun1_v = pow(xnu,-1. + rfield.cutoff[rfield.ipSpec][1])*
146 sexp(xnu/rfield.slope[rfield.ipSpec])*sexp(0.01/
147 xnu);
148 /* only add on x-ray for energies above 0.1 Ryd */
149 if( xnu > 0.1 )
150 {
151 if( xnu < 7350. )
152 {
153 /* cutoff is actually normalization constant
154 * below 100keV continuum is nu-1 */
155 ffun1_v += pow(xnu,rfield.cutoff[rfield.ipSpec][2] -
156 1.)*rfield.cutoff[rfield.ipSpec][0]*sexp(1./
157 xnu);
158 }
159 else
160 {
161 ffun1_v += pow(7350.,rfield.cutoff[rfield.ipSpec][2] -
162 1.)*rfield.cutoff[rfield.ipSpec][0]/
163 POW3(xnu/7350.);
164 }
165 }
166
167 }
168 else if( strcmp(chKey,"POWER") == 0 )
169 {
170 /* power law with cutoff at both ends */
171 ffun1_v = pow(xnu,-1. + rfield.slope[rfield.ipSpec])*
172 sexp(xnu/rfield.cutoff[rfield.ipSpec][0]+rfield.cutoff[rfield.ipSpec][1]/
173 xnu);
174
175 }
176 else if( strcmp(chKey,"DISKB") == 0 )
177 {
178 long numSteps = 100;
179 double TempHi, TempLo;
180 // Let temps take any values. If equal, just do blackbody...
181 if( fp_equal( rfield.slope[rfield.ipSpec], rfield.cutoff[rfield.ipSpec][0] ) )
182 {
183 ffun1_v = PlanckFunction( rfield.slope[rfield.ipSpec], xnu );
184 }
185 else
186 {
187 if( rfield.slope[rfield.ipSpec] > rfield.cutoff[rfield.ipSpec][0] )
188 {
189 TempHi = rfield.slope[rfield.ipSpec];
190 TempLo = rfield.cutoff[rfield.ipSpec][0];
191 }
192 else
193 {
194 TempLo = rfield.slope[rfield.ipSpec];
195 TempHi = rfield.cutoff[rfield.ipSpec][0];
196 }
197 ASSERT( TempLo < TempHi );
198 double LogDeltaT = (log10(TempHi) - log10(TempLo))/(numSteps-1.);
199 ffun1_v = 0.;
200 for( long i=0; i<numSteps; i++ )
201 {
202 double Temp = pow( 10., log10(TempHi) - i * LogDeltaT );
203 double relativeWeight = pow( TempHi/Temp, 2.6666 ) * pow( 10., LogDeltaT );
204 ffun1_v += PlanckFunction( Temp, xnu ) * relativeWeight;
205 }
206 }
207 }
208 else if( strcmp(chKey,"BLACK") == 0 )
209 {
210 ffun1_v = PlanckFunction( rfield.slope[rfield.ipSpec], xnu );
211 }
212 else if( strcmp(chKey,"INTER") == 0 )
213 {
214 /* interpolate on tabulated input spectrum, factor of 1.0001 to
215 * make sure that requested energy is within bounds of array */
216 if( xnu >= rfield.tNu[rfield.ipSpec][0].Ryd()*1.000001 )
217 {
218 /* loop starts at second array element, [1], since want to
219 * find next continuum energy greater than desired point */
220 i = 1;
221 /* up to NCELL tabulated points may be read in. Very fine
222 * continuum mesh such as that output by stellar atmospheres
223 * can have very large number of points */
224 while( i< NCELL-1 && rfield.tNu[rfield.ipSpec][i].Ryd()>0. )
225 {
226 if( xnu < rfield.tNu[rfield.ipSpec][i].Ryd() )
227 {
228 /* the energy xnu is between points rfield.tNuRyd[rfield.ipSpec][i-1]
229 * and rfield.tNuRyd[rfield.ipSpec][i] - do linear
230 * interpolation in log log space */
231 y = rfield.tFluxLog[rfield.ipSpec][i-1] +
232 rfield.tslop[rfield.ipSpec][i-1]*
233 log10(xnu/rfield.tNu[rfield.ipSpec][i-1].Ryd());
234
235 /* return value is photon density, div by energy */
236 ffun1_v = pow(10.,y);
237
238 /* this checks that overshoots did not occur - interpolated
239 * value must be between lowest and highest point */
240# ifndef NDEBUG
241 double ys1 = MIN2( rfield.tFluxLog[rfield.ipSpec][i-1],rfield.tFluxLog[rfield.ipSpec][i]);
242 double ys2 = MAX2( rfield.tFluxLog[rfield.ipSpec][i-1],rfield.tFluxLog[rfield.ipSpec][i]);
243 ys1 = pow( 10. , ys1 );
244 ys2 = pow( 10. , ys2 );
245 ASSERT( ffun1_v >= ys1/(1.+100.*FLT_EPSILON) );
246 ASSERT( ffun1_v <= ys2*(1.+100.*FLT_EPSILON) );
247# endif
248 /* return value is photon density, div by energy */
249 return( ffun1_v/xnu );
250 }
251 ++i;
252 }
253 /* energy above highest in table */
254 ffun1_v = 0.;
255 }
256 else
257 {
258 /* energy below lowest on table */
259 ffun1_v = 0.;
260 }
261 }
262
263 else if( strcmp(chKey,"BREMS") == 0 )
264 {
265 /* brems continuum, rough gaunt factor */
266 fac = TE1RYD*xnu/rfield.slope[rfield.ipSpec];
267 ffun1_v = sexp(fac)/pow(xnu,1.2);
268
269 }
270 else if( strcmp(chKey,"LASER") == 0 )
271 {
272 const double BIG = 1.e10;
273 const double SMALL = 1.e-25;
274 /* a laser, mostly one frequency */
275 /* >>chng 01 jul 01, was hard-wired 0.05 rel frac, change to optional
276 * second parameter, with default of 0.05 */
277 /*if( xnu > 0.95*rfield.slope[rfield.ipSpec] && xnu <
278 1.05*rfield.slope[rfield.ipSpec] )*/
279 if( xnu > (1.-rfield.cutoff[rfield.ipSpec][0])*rfield.slope[rfield.ipSpec] &&
280 xnu < (1.+rfield.cutoff[rfield.ipSpec][0])*rfield.slope[rfield.ipSpec] )
281 {
282 ffun1_v = BIG;
283 }
284 else
285 {
286 ffun1_v = SMALL;
287 }
288
289 }
290 else if( strcmp(chKey,"READ ") == 0 )
291 {
292 /* use array of values read in on TABLE READ command */
293 if( xnu >= rfield.egamry )
294 {
295 ffun1_v = 0.;
296 }
297 else
298 {
299 i = ipoint(xnu);
300 if( i > rfield.nupper || i < 1 )
301 {
302 ffun1_v = 0.;
303 }
304 else
305 {
306 // Fragile ASSERT to ensure consistency -- but should do something more to ensure
307 // consistency anyhow, using the INTERpolate option.
308 realnum tFluxLog = rfield.tFluxLog[rfield.ipSpec][i-1];
309 realnum tNu = (realnum)rfield.tNu[rfield.ipSpec][i-1].Ryd();
310 ASSERT( fp_equal(tFluxLog,(realnum)-70.) ||
311 fp_equal_tol(rfield.anu[i-1],(double)tNu,3e-3*rfield.anu[i-1]) );
312
313 ffun1_v = pow((realnum)10.,rfield.tFluxLog[rfield.ipSpec][i-1])/rfield.anu[i-1];
314 }
315 }
316 }
317
318 /* >>chng 06 jul 10, retired TABLE STARBURST command, PvH */
319 /* >>chng 05 nov 30, retired TABLE TLUSTY command, PvH */
320
321 else if( strcmp(chKey,"VOLK ") == 0 )
322 {
323 /* use array of values read in from Kevin Volk's rebinning of
324 * large atlas grids */
325 if( xnu >= rfield.egamry )
326 {
327 ffun1_v = 0.;
328 }
329 else
330 {
331 i = ipoint(xnu);
332 if( i > rfield.nupper )
333 {
334 fprintf( ioQQQ, " ffun1: Too many points - increase ncell\n" );
335 fprintf( ioQQQ, " cell needed=%4ld ncell=%4ld\n",
336 i, rfield.nupper );
338 }
339 if( i > rfield.nupper || i < 1 )
340 {
341 ffun1_v = 0.;
342 }
343 else
344 {
345 /* bug fixed Jul 9 93: FFUN1 = TSLOP(IPSPEC,I) / ANU(I) / ANU(I)
346 * i has value 939 */
347 ffun1_v = rfield.tslop[rfield.ipSpec][i-1]/ rfield.anu[i-1];
348 }
349 }
350 }
351 else
352 {
353 fprintf( ioQQQ, " ffun1: I do not understand continuum label \"%s\" for continuum %li.\n",
354 chKey , rfield.ipSpec);
356 }
357
358 if( ffun1_v > 1e35 && !lgWarn )
359 {
360 lgWarn = true;
361 fprintf( ioQQQ, " FFUN1: Continuum %ld is very intense.\n",
362 rfield.ipSpec );
363 fprintf( ioQQQ, " I will try to press on, but may have problems.\n" );
364 }
365 return( ffun1_v );
366}
367
368double PlanckFunction( double Temp, double E_Ryd )
369{
370 const double db_log = log(DBL_MAX);
371 double ffun1_v;
372 double fac;
373
374 /* black body */
375 fac = TE1RYD*E_Ryd/Temp;
376 /* >>>chng 00 apr 13 from 80 to log(dbl_max) */
377 if( fac > db_log )
378 {
379 ffun1_v = 0.;
380 }
381 else if( fac > 1.e-5 )
382 {
383 ffun1_v = E_Ryd*E_Ryd/(exp(fac) - 1.);
384 }
385 else
386 {
387 ffun1_v = E_Ryd*E_Ryd/(fac*(1. + fac/2.));
388 }
389
390 return ffun1_v;
391}
392/*outsum sum outward continuum beams */
393void outsum(double *outtot, double *outin, double *outout)
394{
395 long int i;
396
397 DEBUG_ENTRY( "outsum()" );
398
399 *outin = 0.;
400 *outout = 0.;
401 for( i=0; i < rfield.nflux; i++ )
402 {
403 /* N.B. in following en1ryd prevents overflow */
404 *outin += rfield.anu[i]*(rfield.flux[0][i]*EN1RYD);
405 *outout += rfield.anu[i]*(rfield.outlin[0][i] + rfield.outlin_noplot[i] +rfield.ConInterOut[i])*
406 EN1RYD;
407 }
408
409 *outtot = *outin + *outout;
410 return;
411}
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
sys_float sexp(sys_float x)
Definition service.cpp:914
#define POW3
Definition cddefines.h:936
#define MIN2
Definition cddefines.h:761
#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
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
bool fp_equal_tol(sys_float x, sys_float y, sys_float tol)
Definition cddefines.h:854
double ffun1(double xnu)
double ffun(double anu)
Definition cont_ffun.cpp:17
void outsum(double *outtot, double *outin, double *outout)
double PlanckFunction(double Temp, double E_Ryd)
long ipoint(double energy_ryd)
UNUSED const realnum BIGFLOAT
Definition cpu.h:189
const realnum SMALLFLOAT
Definition cpu.h:191
UNUSED const double EN1RYD
Definition physconst.h:179
UNUSED const double TE1RYD
Definition physconst.h:183
t_rfield rfield
Definition rfield.cpp:8
const int NCELL
Definition rfield.h:21
static const double BIG