cloudy  trunk
 All Data Structures 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-2008 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 
13 /*ReadTable called by TABLE READ to read in continuum from PUNCH TRANSMITTED CONTINUUM */
14 STATIC void ReadTable(const string& fnam);
15 
16 /* evaluate sum of all individual continua at one energy, return is
17 * continuum intensity */
18 double ffun(
19  /* the energy in Rydbergs where the continuum will be evaluated */
20  double anu )
21 {
22  double frac_beam_time;
23  /* fraction of beamed continuum that is constant */
24  double frac_beam_const;
25  /* fraction of continuum that is isotropic */
26  double frac_isotropic;
27  double a;
28 
29  DEBUG_ENTRY( "ffun()" );
30 
31  /* call real function */
32  a = ffun( anu , &frac_beam_time , &frac_beam_const , &frac_isotropic );
33  return a;
34 }
35 
36 /* evaluate sum of all individual continua at one energy, return is
37  * continuum intensity */
38 double ffun(
39  /* the energy in Rydbergs where the continuum will be evaluated */
40  double anu ,
41  /* fraction of beamed continuum that is varies with time */
42  double *frac_beam_time,
43  /* fraction of beamed continuum that is constant */
44  double *frac_beam_const,
45  /* fraction of continuum that is isotropic */
46  double *frac_isotropic )
47 {
48  double ffun_v;
49  static bool lgWarn = false;
50  double flx_beam_time , flx_beam_const , flx_isotropic;
51 
52  DEBUG_ENTRY( "ffun()" );
53 
54  /* This routine, ffun, returns the sum of photons per unit time, area, energy,
55  * for all continua in the calculation. ffun1 is called and returns
56  * a single continuum. We loop over all nspec continuum sources -
57  * ipspec points to each */
58  ffun_v = 0.;
59  flx_beam_time = 0.;
60  flx_beam_const = 0.;
61  flx_isotropic = 0.;
63  {
64  double one = ffun1(anu)*rfield.spfac[rfield.ipspec];
65  ffun_v += one;
66 
67  /* find fraction of total that is constant vs variable and
68  * isotropic vs beamed */
70  {
72  flx_beam_time += one;
73  else
74  flx_beam_const += one;
75  }
76  else
77  flx_isotropic += one;
78  }
79 
80  /* at this point rfield.flux is the sum of the following three continua
81  * now keep track of the three different types */
82  if( ffun_v < SMALLFLOAT )
83  {
84  *frac_beam_time = 0.;
85  *frac_beam_const = 1.;
86  *frac_isotropic = 0.;
87  }
88  else
89  {
90  /* fraction of beamed continuum that varies with time */
91  *frac_beam_time = flx_beam_time / ffun_v;
92  /* part of beamed continuum that is constant */
93  *frac_beam_const = flx_beam_const / ffun_v;
94  /* the constant isotropic continuum */
95  *frac_isotropic = flx_isotropic / ffun_v;
96  }
97  ASSERT( *frac_beam_time >=0. && *frac_beam_time<=1.+3.*DBL_EPSILON );
98  ASSERT( *frac_beam_const >=0.&& *frac_beam_const<=1.+3.*DBL_EPSILON );
99  ASSERT( *frac_isotropic >=0. && *frac_isotropic<=1.+3.*DBL_EPSILON );
100  ASSERT( fabs( 1.-*frac_beam_time-*frac_beam_const-*frac_isotropic)<
101  10.*DBL_EPSILON);
102 
103  if( ffun_v > BIGFLOAT && !lgWarn )
104  {
105  lgWarn = true;
106  fprintf( ioQQQ, " FFUN: The net continuum is very intense.\n" );
107  fprintf( ioQQQ, " I will try to press on, but may have problems.\n" );
108  }
109  return( ffun_v );
110 }
111 
112 /*ffun1 derive flux at a specific energy, for one continuum */
113 double ffun1(double xnu)
114 {
115  char chKey[6];
116  long int i;
117  double fac,
118  ffun1_v,
119  y;
120  static bool lgWarn = false;
121 
122  DEBUG_ENTRY( "ffun1()" );
123 
124 
125  /* confirm that pointer is within range */
126  ASSERT( rfield.ipspec >= 0);
128 
129  /* FFUN1 returns photons per unit time, area, energy, for one continuum
130  * ipspec is the pointer to the continuum source, in the order
131  * entered on the command lines */
132 
133  /*begin sanity check */
134  ASSERT( xnu >= rfield.emm*0.99 );
135  ASSERT( xnu <= rfield.egamry*1.01 );
136  /*end sanity check */
137 
138  strcpy( chKey, rfield.chSpType[rfield.ipspec] );
139 
140  if( strcmp(chKey,"AGN ") == 0 )
141  {
142  /* power law with cutoff at both ends
143  * nomenclature all screwed up - slope is cutoff energy in Ryd,
144  * cutoff[1][i] is ratio of two continua from alpha ox
145  * cutoff[2][i] is slope of bb temp */
146  ffun1_v = pow(xnu,-1. + rfield.cutoff[rfield.ipspec][1])*
147  sexp(xnu/rfield.slope[rfield.ipspec])*sexp(0.01/
148  xnu);
149  /* only add on x-ray for energies above 0.1 Ryd */
150  if( xnu > 0.1 )
151  {
152  if( xnu < 7350. )
153  {
154  /* cutoff is actually normalization constant
155  * below 100keV continuum is nu-1 */
156  ffun1_v += pow(xnu,rfield.cutoff[rfield.ipspec][2] -
157  1.)*rfield.cutoff[rfield.ipspec][0]*sexp(1./
158  xnu);
159  }
160  else
161  {
162  ffun1_v += pow(7350.,rfield.cutoff[rfield.ipspec][2] -
163  1.)*rfield.cutoff[rfield.ipspec][0]/
164  POW3(xnu/7350.);
165  }
166  }
167 
168  }
169  else if( strcmp(chKey,"POWER") == 0 )
170  {
171  /* power law with cutoff at both ends */
172  ffun1_v = pow(xnu,-1. + rfield.slope[rfield.ipspec])*
174  xnu);
175 
176  }
177  else if( strcmp(chKey,"BLACK") == 0 )
178  {
179  const double db_log = log(DBL_MAX);
180 
181  /* black body */
182  fac = TE1RYD*xnu/rfield.slope[rfield.ipspec];
183  /* >>>chng 00 apr 13 from 80 to log(dbl_max) */
184  if( fac > db_log )
185  {
186  ffun1_v = 0.;
187  }
188  else if( fac > 1.e-5 )
189  {
190  ffun1_v = xnu*xnu/(exp(fac) - 1.);
191  }
192  else
193  {
194  ffun1_v = xnu*xnu/(fac*(1. + fac/2.));
195  }
196 
197  }
198  else if( strcmp(chKey,"INTER") == 0 )
199  {
200  /* interpolate on tabulated input spectrum, factor of 1.0001 to
201  * make sure that requested energy is within bounds of array */
202  if( xnu >= rfield.tNuRyd[rfield.ipspec][0]*1.000001 )
203  {
204  /* loop starts at second array element, [1], since want to
205  * find next continuum energy greater than desired point */
206  i = 1;
207  /* up to NCELL tabulated points may be read in. Very fine
208  * continuum mesh such as that output by stellar atmospheres
209  * can have very large number of points */
210  while( i< NCELL-1 && rfield.tNuRyd[rfield.ipspec][i]>0. )
211  {
212  if( xnu < rfield.tNuRyd[rfield.ipspec][i] )
213  {
214  /* the energy xnu is between points rfield.tNuRyd[rfield.ipspec][i-1]
215  * and rfield.tNuRyd[rfield.ipspec][i] - do linear
216  * interpolation in log log space */
217  y = rfield.tFluxLog[rfield.ipspec][i-1] +
218  rfield.tslop[rfield.ipspec][i-1]*
219  log10(xnu/rfield.tNuRyd[rfield.ipspec][i-1]);
220 
221  /* return value is photon density, div by energy */
222  ffun1_v = pow(10.,y);
223 
224  /* this checks that overshoots did not occur - interpolated
225  * value must be between lowest and highest point */
226 # ifndef NDEBUG
227  double ys1 = MIN2( rfield.tFluxLog[rfield.ipspec][i-1],rfield.tFluxLog[rfield.ipspec][i]);
228  double ys2 = MAX2( rfield.tFluxLog[rfield.ipspec][i-1],rfield.tFluxLog[rfield.ipspec][i]);
229  ys1 = pow( 10. , ys1 );
230  ys2 = pow( 10. , ys2 );
231  ASSERT( ffun1_v >= ys1/(1.+100.*FLT_EPSILON) );
232  ASSERT( ffun1_v <= ys2*(1.+100.*FLT_EPSILON) );
233 # endif
234  /* return value is photon density, div by energy */
235  return( ffun1_v/xnu );
236  }
237  ++i;
238  }
239  /* energy above highest in table */
240  ffun1_v = 0.;
241  }
242  else
243  {
244  /* energy below lowest on table */
245  ffun1_v = 0.;
246  }
247  }
248 
249  else if( strcmp(chKey,"BREMS") == 0 )
250  {
251  /* brems continuum, rough gaunt factor */
252  fac = TE1RYD*xnu/rfield.slope[rfield.ipspec];
253  ffun1_v = sexp(fac)/pow(xnu,1.2);
254 
255  }
256  else if( strcmp(chKey,"LASER") == 0 )
257  {
258  const double BIG = 1.e10;
259  const double SMALL = 1.e-25;
260  /* a laser, mostly one frequency */
261  /* >>chng 01 jul 01, was hard-wired 0.05 rel frac, change to optional
262  * second parameter, with default of 0.05 */
263  /*if( xnu > 0.95*rfield.slope[rfield.ipspec] && xnu <
264  1.05*rfield.slope[rfield.ipspec] )*/
265  if( xnu > (1.-rfield.cutoff[rfield.ipspec][0])*rfield.slope[rfield.ipspec] &&
266  xnu < (1.+rfield.cutoff[rfield.ipspec][0])*rfield.slope[rfield.ipspec] )
267  {
268  ffun1_v = BIG;
269  }
270  else
271  {
272  ffun1_v = SMALL;
273  }
274 
275  }
276  else if( strcmp(chKey,"READ ") == 0 )
277  {
278  /* implement the table read command
279  * if this is first time called then we need to read in file */
280  if( rfield.tslop[rfield.ipspec][0] == 0. )
281  {
282  /* >>chng 01 nov 01, array index for below was bogus, only worked for a single continuum */
284  rfield.tslop[rfield.ipspec][0] = 1.;
285  }
286  /* use array of values read in on TABLE READ command */
287  if( xnu >= rfield.egamry )
288  {
289  ffun1_v = 0.;
290  }
291  else
292  {
293  i = ipoint(xnu);
294  if( i > rfield.nupper || i < 1 )
295  {
296  ffun1_v = 0.;
297  }
298  else
299  {
300  ffun1_v = rfield.ConTabRead[i-1]/rfield.anu[i-1]/rfield.anu[i-1];
301  }
302  }
303  }
304 
305  /* >>chng 06 jul 10, retired TABLE STARBURST command, PvH */
306  /* >>chng 05 nov 30, retired TABLE TLUSTY command, PvH */
307 
308  else if( strcmp(chKey,"VOLK ") == 0 )
309  {
310  /* use array of values read in from Kevin Volk's rebinning of
311  * large atlas grids */
312  if( xnu >= rfield.egamry )
313  {
314  ffun1_v = 0.;
315  }
316  else
317  {
318  i = ipoint(xnu);
319  if( i > rfield.nupper )
320  {
321  fprintf( ioQQQ, " ffun1: Too many points - increase ncell\n" );
322  fprintf( ioQQQ, " cell needed=%4ld ncell=%4ld\n",
323  i, rfield.nupper );
324  cdEXIT(EXIT_FAILURE);
325  }
326  if( i > rfield.nupper || i < 1 )
327  {
328  ffun1_v = 0.;
329  }
330  else
331  {
332  /* bug fixed Jul 9 93: FFUN1 = TSLOP(IPSPEC,I) / ANU(I) / ANU(I)
333  * i has value 939 */
334  ffun1_v = rfield.tslop[rfield.ipspec][i-1]/ rfield.anu[i-1];
335  }
336  }
337  }
338  else
339  {
340  fprintf( ioQQQ, " ffun1: I do not understand continuum label \"%s\" for continuum %li.\n",
341  chKey , rfield.ipspec);
342  cdEXIT(EXIT_FAILURE);
343  }
344 
345  if( ffun1_v > 1e35 && !lgWarn )
346  {
347  lgWarn = true;
348  fprintf( ioQQQ, " FFUN1: Continuum %ld is very intense.\n",
349  rfield.ipspec );
350  fprintf( ioQQQ, " I will try to press on, but may have problems.\n" );
351  }
352  return( ffun1_v );
353 }
354 
355 /*ReadTable called by TABLE READ to read in continuum from PUNCH TRANSMITTED CONTINUUM */
356 STATIC void ReadTable(const string& fnam)
357 {
358  char chLine[INPUT_LINE_LENGTH];
359  long int i,
360  nPoints;
361  double Differ,
362  EnerLast;
363  FILE *io;
364 
365  DEBUG_ENTRY( "ReadTable()" );
366 
367  io = open_data( fnam.c_str(), "r", AS_LOCAL_ONLY );
368 
369  /* read in first line of header */
370  if( read_whole_line( chLine , (int)sizeof(chLine) , io )==NULL )
371  {
372  fprintf( ioQQQ, " error 1 reading input continuum file.\n" );
373  cdEXIT(EXIT_FAILURE);
374  }
375 
376  /* now read in the file of numbers */
377  i = 0;
378  /* keep reading until we hit eol or run out of room in the array */
379  while( (read_whole_line(chLine, (int)sizeof(chLine),io)!=NULL) && (i<rfield.nupper) )
380  {
381  double help[2];
382  sscanf( chLine, "%lf%lf ", &help[0], &help[1] );
383  opac.tmn[i] = (realnum)help[0];
384  rfield.ConTabRead[i] = (realnum)help[1];
385  ++i;
386  }
387  /* put pointer at last good value */
388  nPoints = i - 1;
389 
390  /* check that sane number of values entered */
391  if( nPoints < 1 )
392  {
393  fprintf( ioQQQ, " ReadTable, file for TABLE READ has too few points, =%5ld\n",
394  nPoints );
395  cdEXIT(EXIT_FAILURE);
396  }
397 
398  /* check on units of energy scale, convert to Rydbergs if necessary
399  * tmn is scale read in from punch file, anu is correct scale
400  * EnerLast is energy of last point in Rydbergs */
401  EnerLast = opac.tmn[nPoints];
402  if( fabs(opac.tmn[0]-rfield.anu[0])/rfield.anu[0] > 1e-3 )
403  {
404  /* first energy does not appear to have been in Rydbergs */
405  if( opac.tmn[0] <= 0. )
406  {
407  fprintf( ioQQQ,
408  " DISASTER ReadTable: the first energy in table read file is not positive. Something is wrong.\n" );
409  cdEXIT(EXIT_FAILURE);
410  }
411  else if( fabs(opac.tmn[0]-RYDLAM/rfield.anu[0]*1e-4)/opac.tmn[0] <
412  1e-3 )
413  {
414  /* wavelength in microns */
415  EnerLast = RYDLAM/opac.tmn[nPoints]/1e4;
416  }
417  else if( fabs(opac.tmn[0]-RYDLAM/rfield.anu[0])/opac.tmn[0] <
418  1e-3 )
419  {
420  /* wavelength in Angstroms */
421  EnerLast = RYDLAM/opac.tmn[nPoints];
422  }
423  else if( fabs(opac.tmn[0]-rfield.anu[0]*EVRYD*1e-3)/opac.tmn[0] <
424  1e-3 )
425  {
426  /* wavelength in keV */
427  EnerLast = opac.tmn[nPoints]/EVRYD/1e-3;
428  }
429  else if( fabs(opac.tmn[0]-rfield.anu[0]*EVRYD)/opac.tmn[0] <
430  1e-3 )
431  {
432  /* wavelength in eV */
433  EnerLast = opac.tmn[nPoints]/EVRYD;
434  }
435  else
436  {
437  fprintf( ioQQQ, " DISASTER ReadTable: the energy scale in the table read file makes no sense to me.\n" );
438  cdEXIT(EXIT_FAILURE);
439  }
440  }
441 
442  /* now check now the energies of the highest points agree */
443  Differ = fabs(EnerLast-rfield.anu[nPoints])/rfield.anu[nPoints];
444  if( Differ > 0.001 )
445  {
446  fprintf( ioQQQ, " DISASTER ReadTable: The energy mesh of the file read in by the TABLE READ command does not agree with this version of Cloudy.\n" );
447  fprintf( ioQQQ, " DISASTER ReadTable: Was the file generated by an older version of the code?\n" );
448  fprintf( ioQQQ, " DISASTER ReadTable: Did the first model do more than one iteration, but the LAST option was missing on PUNCH LAST TRANSMITTED CONTINUUM?\n");
449  fprintf( ioQQQ, " DISASTER ReadTable: Number of points read in=%5ld\n",
450  nPoints );
451  fprintf( ioQQQ, " ReadTable: input, internal energies=%12.4e%12.4e\n",
452  opac.tmn[nPoints], rfield.anu[nPoints] );
453  cdEXIT(EXIT_FAILURE);
454  }
455 
456  /* make sure rest of array has valid zeros */
457  for( i=nPoints + 1; i < rfield.nupper; i++ )
458  {
459  rfield.ConTabRead[i] = 0.;
460  }
461 
462  fclose(io);
463  return;
464 }
465 

Generated for cloudy by doxygen 1.8.1.2