cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
cont_ffun.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2017 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 "rfield.h"
8 #include "continuum.h"
9 
10 double PlanckFunction( double Temp, double E_Ryd );
11 
12 /* evaluate sum of all individual continua at one energy, return is
13 * continuum intensity */
14 double ffun(
15  /* the energy in Rydbergs where the continuum will be evaluated */
16  double anu )
17 {
18  double frac_beam_time;
19  /* fraction of beamed continuum that is constant */
20  double frac_beam_const;
21  /* fraction of continuum that is isotropic */
22  double frac_isotropic;
23  double a;
24 
25  DEBUG_ENTRY( "ffun()" );
26 
27  /* call real function */
28  a = ffun( anu , &frac_beam_time , &frac_beam_const , &frac_isotropic );
29  return a;
30 }
31 
32 /* evaluate sum of all individual continua at one energy, return is
33  * continuum intensity */
34 double ffun(
35  /* the energy in Rydbergs where the continuum will be evaluated */
36  double anu ,
37  /* fraction of beamed continuum that is varies with time */
38  double *frac_beam_time,
39  /* fraction of beamed continuum that is constant */
40  double *frac_beam_const,
41  /* fraction of continuum that is isotropic */
42  double *frac_isotropic )
43 {
44  double ffun_v;
45  static bool lgWarn = false;
46  double flx_beam_time , flx_beam_const , flx_isotropic;
47 
48  DEBUG_ENTRY( "ffun()" );
49 
50  /* This routine, ffun, returns the sum of photons per unit time, area, energy,
51  * for all continua in the calculation. ffun1 is called and returns
52  * a single continuum. We loop over all nShape continuum sources -
53  * ipspec points to each */
54  ffun_v = 0.;
55  flx_beam_time = 0.;
56  flx_beam_const = 0.;
57  flx_isotropic = 0.;
59  {
60  double one = ffun1(anu)*rfield.spfac[rfield.ipSpec];
61  ffun_v += one;
62 
63  /* find fraction of total that is constant vs variable and
64  * isotropic vs beamed */
66  {
68  flx_beam_time += one;
69  else
70  flx_beam_const += one;
71  }
72  else
73  flx_isotropic += one;
74  }
75 
76  /* at this point rfield.flux is the sum of the following three continua
77  * now keep track of the three different types */
78  if( ffun_v < SMALLDOUBLE )
79  {
80  *frac_beam_time = 0.;
81  *frac_beam_const = 1.;
82  *frac_isotropic = 0.;
83  ffun_v = 0.;
84  }
85  else
86  {
87  /* fraction of beamed continuum that varies with time */
88  *frac_beam_time = flx_beam_time / ffun_v;
89  /* part of beamed continuum that is constant */
90  *frac_beam_const = flx_beam_const / ffun_v;
91  /* the constant isotropic continuum */
92  *frac_isotropic = flx_isotropic / ffun_v;
93  }
94  ASSERT( *frac_beam_time >=0. && *frac_beam_time<=1.+3.*DBL_EPSILON );
95  ASSERT( *frac_beam_const >=0.&& *frac_beam_const<=1.+3.*DBL_EPSILON );
96  ASSERT( *frac_isotropic >=0. && *frac_isotropic<=1.+3.*DBL_EPSILON );
97  ASSERT( fabs( 1.-*frac_beam_time-*frac_beam_const-*frac_isotropic)<
98  10.*DBL_EPSILON);
99 
100  if( ffun_v > BIGFLOAT && !lgWarn )
101  {
102  lgWarn = true;
103  fprintf( ioQQQ, " FFUN: The net continuum is very intense.\n" );
104  fprintf( ioQQQ, " I will try to press on, but may have problems.\n" );
105  }
106  return ffun_v;
107 }
108 
109 /*ffun1 derive flux at a specific energy, for one continuum */
110 double ffun1(double xnu)
111 {
112  char chKey[6];
113  long int i;
114  double fac,
115  ffun1_v,
116  y;
117  static bool lgWarn = false;
118 
119  DEBUG_ENTRY( "ffun1()" );
120 
121 
122  /* confirm that pointer is within range */
123  ASSERT( rfield.ipSpec >= 0 );
125 
126  /* FFUN1 returns photons per unit time, area, energy, for one continuum
127  * ipspec is the pointer to the continuum source, in the order
128  * entered on the command lines */
129 
130  /*begin sanity check */
131  ASSERT( xnu > rfield.emm() );
132  ASSERT( xnu < rfield.egamry() );
133  /*end sanity check */
134 
135  strcpy( chKey, rfield.chSpType[rfield.ipSpec] );
136 
137  if( strcmp(chKey,"AGN ") == 0 )
138  {
139  /* power law with cutoff at both ends
140  * nomenclature all screwed up - slope is cutoff energy in Ryd,
141  * cutoff[1][i] is ratio of two continua from alpha ox
142  * cutoff[2][i] is slope of bb temp */
143  ffun1_v = pow(xnu,-1. + rfield.cutoff[rfield.ipSpec][1])*
144  sexp(xnu/rfield.slope[rfield.ipSpec])*sexp(0.01/
145  xnu);
146  /* only add on x-ray for energies above 0.1 Ryd */
147  if( xnu > 0.1 )
148  {
149  if( xnu < 7350. )
150  {
151  /* cutoff is actually normalization constant
152  * below 100keV continuum is nu-1 */
153  ffun1_v += pow(xnu,rfield.cutoff[rfield.ipSpec][2] -
154  1.)*rfield.cutoff[rfield.ipSpec][0]*sexp(1./
155  xnu);
156  }
157  else
158  {
159  ffun1_v += pow(7350.,rfield.cutoff[rfield.ipSpec][2] -
160  1.)*rfield.cutoff[rfield.ipSpec][0]/
161  POW3(xnu/7350.);
162  }
163  }
164 
165  }
166  else if( strcmp(chKey,"POWER") == 0 )
167  {
168  /* power law with cutoff at both ends */
169  ffun1_v = pow(xnu,-1. + rfield.slope[rfield.ipSpec])*
171  xnu);
172 
173  }
174  else if( strcmp(chKey,"DISKB") == 0 )
175  {
176  long numSteps = 100;
177  double TempHi, TempLo;
178  // Let temps take any values. If equal, just do blackbody...
180  {
181  ffun1_v = PlanckFunction( rfield.slope[rfield.ipSpec], xnu );
182  }
183  else
184  {
186  {
187  TempHi = rfield.slope[rfield.ipSpec];
188  TempLo = rfield.cutoff[rfield.ipSpec][0];
189  }
190  else
191  {
192  TempLo = rfield.slope[rfield.ipSpec];
193  TempHi = rfield.cutoff[rfield.ipSpec][0];
194  }
195  ASSERT( TempLo < TempHi );
196  double LogDeltaT = (log10(TempHi) - log10(TempLo))/(numSteps-1.);
197  ffun1_v = 0.;
198  for( long i=0; i<numSteps; i++ )
199  {
200  double Temp = exp10( log10(TempHi) - double(i) * LogDeltaT );
201  double relativeWeight = powpq( TempHi/Temp, 8, 3 ) * exp10( LogDeltaT );
202  ffun1_v += PlanckFunction( Temp, xnu ) * relativeWeight;
203  }
204  }
205  }
206  else if( strcmp(chKey,"BLACK") == 0 )
207  {
208  ffun1_v = PlanckFunction( rfield.slope[rfield.ipSpec], xnu );
209  }
210  else if( strcmp(chKey,"INTER") == 0 )
211  {
212  /* interpolate on tabulated input spectrum, factor of 1.0001 to
213  * make sure that requested energy is within bounds of array */
214  if( xnu >= rfield.tNu[rfield.ipSpec][0].Ryd()*1.000001 )
215  {
216  /* loop starts at second array element, [1], since want to
217  * find next continuum energy greater than desired point */
218  for( i=1; i < rfield.ncont[rfield.ipSpec]; ++i )
219  {
220  if( xnu < rfield.tNu[rfield.ipSpec][i].Ryd() )
221  {
222  /* the energy xnu is between points rfield.tNuRyd[rfield.ipSpec][i-1]
223  * and rfield.tNuRyd[rfield.ipSpec][i] - do linear
224  * interpolation in log log space */
225  y = rfield.tFluxLog[rfield.ipSpec][i-1] +
226  rfield.tslop[rfield.ipSpec][i-1]*
227  log10(xnu/rfield.tNu[rfield.ipSpec][i-1].Ryd());
228 
229  /* return value is photon density, div by energy */
230  ffun1_v = exp10(y);
231 
232  /* this checks that overshoots did not occur - interpolated
233  * value must be between lowest and highest point */
234 # ifndef NDEBUG
235  double ys1 = MIN2( rfield.tFluxLog[rfield.ipSpec][i-1],rfield.tFluxLog[rfield.ipSpec][i]);
236  double ys2 = MAX2( rfield.tFluxLog[rfield.ipSpec][i-1],rfield.tFluxLog[rfield.ipSpec][i]);
237  ys1 = exp10( ys1 );
238  ys2 = exp10( ys2 );
239  ASSERT( ffun1_v >= ys1/(1.+100.*FLT_EPSILON) );
240  ASSERT( ffun1_v <= ys2*(1.+100.*FLT_EPSILON) );
241 # endif
242  /* return value is photon density, div by energy */
243  return( ffun1_v/xnu );
244  }
245  }
246  /* energy above highest in table */
247  ffun1_v = 0.;
248  }
249  else
250  {
251  /* energy below lowest on table */
252  ffun1_v = 0.;
253  }
254  }
255  else if( strcmp(chKey,"BREMS") == 0 )
256  {
257  /* brems continuum, rough gaunt factor */
258  fac = TE1RYD*xnu/rfield.slope[rfield.ipSpec];
259  ffun1_v = sexp(fac)/pow(xnu,1.2);
260 
261  }
262  else if( strcmp(chKey,"LASER") == 0 )
263  {
264  const double BIG = 1.e10;
265  /* a laser, mostly one frequency */
266  /* >>chng 01 jul 01, was hard-wired 0.05 rel frac, change to optional
267  * second parameter, with default of 0.05 */
268  /*if( xnu > 0.95*rfield.slope[rfield.ipSpec] && xnu <
269  1.05*rfield.slope[rfield.ipSpec] )*/
270  if( xnu > (1.-rfield.cutoff[rfield.ipSpec][0])*rfield.slope[rfield.ipSpec] &&
271  xnu < (1.+rfield.cutoff[rfield.ipSpec][0])*rfield.slope[rfield.ipSpec] )
272  {
273  ffun1_v = BIG;
274  }
275  else
276  {
277  ffun1_v = 0.;
278  }
279 
280  }
281  else if( strcmp(chKey,"READ ") == 0 )
282  {
283  /* use array of values read in on TABLE READ command */
284  // this table is already on the Cloudy mesh, so no need to interpolate...
285  // this table is not defined for rfield.anu(rfield.nflux) (i.e. the cell for the unit test)
286  // we return 0 in that case since the flux value in that cell is irrelevant anyway
287  if( xnu >= rfield.anumax(rfield.nflux-1) )
288  {
289  ffun1_v = 0.;
290  }
291  else
292  {
293  i = rfield.ipointC(xnu);
295  ffun1_v = exp10((double)rfield.tFluxLog[rfield.ipSpec][i])/rfield.anu(i);
296  }
297  }
298  else if( strcmp(chKey,"VOLK ") == 0 )
299  {
300  /* use array of values from the TABLE STAR or TABLE HM12 command */
301  // this table is already on the Cloudy mesh, so no need to interpolate...
302  i = rfield.ipointC(xnu);
304  /* bug fixed Jul 9 93: FFUN1 = TSLOP(IPSPEC,I) / ANU(I) / ANU(I) */
305  ffun1_v = rfield.tslop[rfield.ipSpec][i]/rfield.anu(i);
306  }
307  else
308  {
309  fprintf( ioQQQ, " ffun1: I do not understand continuum label \"%s\" for continuum %li.\n",
310  chKey , rfield.ipSpec);
312  }
313 
314  if( ffun1_v > 1e35 && !lgWarn )
315  {
316  lgWarn = true;
317  fprintf( ioQQQ, " FFUN1: Continuum %ld is very intense.\n",
318  rfield.ipSpec );
319  fprintf( ioQQQ, " I will try to press on, but may have problems.\n" );
320  }
321  return ffun1_v;
322 }
323 
324 double PlanckFunction( double Temp, double E_Ryd )
325 {
326  const double db_log = log(DBL_MAX);
327  double ffun1_v;
328  double fac;
329 
330  /* black body */
331  fac = TE1RYD*E_Ryd/Temp;
332  /* >>>chng 00 apr 13 from 80 to log(dbl_max) */
333  if( fac > db_log )
334  {
335  ffun1_v = 0.;
336  }
337  else
338  {
339  ffun1_v = E_Ryd*E_Ryd/expm1(fac);
340  }
341 
342  return ffun1_v;
343 }
344 /*outsum sum outward continuum beams */
345 void outsum(double *outtot, double *outin, double *outout)
346 {
347  long int i;
348 
349  DEBUG_ENTRY( "outsum()" );
350 
351  *outin = 0.;
352  *outout = 0.;
353  for( i=0; i < rfield.nflux; i++ )
354  {
355  /* N.B. in following en1ryd prevents overflow */
356  *outin += rfield.anu(i)*(rfield.flux[0][i]*EN1RYD);
357  *outout += rfield.anu(i)*(rfield.outlin[0][i] + rfield.outlin_noplot[i] +rfield.ConInterOut[i])*
358  EN1RYD;
359  }
360 
361  *outtot = *outin + *outout;
362  return;
363 }
double emm() const
Definition: mesh.h:84
double PlanckFunction(double Temp, double E_Ryd)
Definition: cont_ffun.cpp:324
bool lgBeamed[LIMSPC]
Definition: rfield.h:296
double ffun(double anu)
Definition: cont_ffun.cpp:14
double exp10(double x)
Definition: cddefines.h:1383
realnum ** flux
Definition: rfield.h:70
realnum * outlin_noplot
Definition: rfield.h:191
const double SMALLDOUBLE
Definition: cpu.h:250
vector< Energy > tNu[LIMSPC]
Definition: rfield.h:316
double ffun1(double xnu)
Definition: cont_ffun.cpp:110
vector< realnum > tFluxLog[LIMSPC]
Definition: rfield.h:319
sys_float sexp(sys_float x)
Definition: service.cpp:1095
realnum ** outlin
Definition: rfield.h:191
FILE * ioQQQ
Definition: cddefines.cpp:7
long ncont[LIMSPC]
Definition: rfield.h:321
double spfac[LIMSPC]
Definition: rfield.h:286
vector< realnum > tslop[LIMSPC]
Definition: rfield.h:317
#define MIN2(a, b)
Definition: cddefines.h:807
double anu(size_t i) const
Definition: mesh.h:111
double cutoff[LIMSPC][3]
Definition: rfield.h:286
bool lgTimeVary[LIMSPC]
Definition: rfield.h:292
long int nflux_with_check
Definition: rfield.h:51
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:858
size_t ipointC(double anu) const
Definition: mesh.h:152
double slope[LIMSPC]
Definition: rfield.h:286
t_rfield rfield
Definition: rfield.cpp:9
realnum * ConInterOut
Definition: rfield.h:156
#define EXIT_FAILURE
Definition: cddefines.h:168
const realnum BIGFLOAT
Definition: cpu.h:244
#define cdEXIT(FAIL)
Definition: cddefines.h:484
static const double BIG
#define ASSERT(exp)
Definition: cddefines.h:617
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
double powpq(double x, int p, int q)
Definition: service.cpp:726
double egamry() const
Definition: mesh.h:88
#define MAX2(a, b)
Definition: cddefines.h:828
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1217
void outsum(double *outtot, double *outin, double *outout)
Definition: cont_ffun.cpp:345
double pow(double x, int i)
Definition: cddefines.h:782
long int nShape
Definition: rfield.h:308
long int ipSpec
Definition: rfield.h:308
double anumax(size_t i) const
Definition: mesh.h:143
#define POW3
Definition: cddefines.h:990
long int nflux
Definition: rfield.h:48
char chSpType[LIMSPC][6]
Definition: rfield.h:333