cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
opacity_createall.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 /*OpacityCreateAll compute initial set of opacities for all species */
4 /*OpacityCreate1Element generate ionic subshell opacities by calling t_ADfA::Inst().phfit */
5 /*opacity_more_memory allocate more memory for opacity stack */
6 /*Opacity_iso_photo_cs returns photoionization cross section for isoelectronic sequences */
7 /*hmiopc derive total H- H minus opacity */
8 /*rayleh compute Rayleigh scattering cross section for Lya */
9 /*OpacityValenceRescale routine to rescale non-OP valence shell cross sections */
10 /*Yan_H2_CS - cross sections for the photoionization of H2 */
11 /******************************************************************************
12  *NB NB NB NB NB NB NB NB NB NB NB NB NB NB
13  * everything set here must be written to the opacity store files
14  *
15  ****************************************************************************** */
16 #include "cddefines.h"
17 #include "physconst.h"
18 #include "dense.h"
19 #include "continuum.h"
20 #include "iso.h"
21 #include "hydrogenic.h"
22 #include "oxy.h"
23 #include "trace.h"
24 #include "heavy.h"
25 #include "rfield.h"
26 #include "hmi.h"
27 #include "atmdat.h"
28 #include "punch.h"
29 #include "grains.h"
30 #include "thirdparty.h"
31 #include "hydro_bauman.h"
32 #include "opacity.h"
33 #include "helike_recom.h"
34 
35 static const int NCSH2P = 10;
36 
37 /* limit to number of opacity cells available in the opacity stack*/
38 static long int ndimOpacityStack = 2100000L;
39 
40 /*Yan_H2_CS - cross sections for the photoionization of H2
41  * may be represented analytically in the paper
42  *>>refer H2 photo cs Yan, M.; Sadeghpour, H. R.; Dalgarno, A. 1998, ApJ, 496, 1044
43  * This is revised version given in ERRATUM 2001, ApJ, 559, 1194
44  * return value is photo cs in cm-2 */
45 STATIC double Yan_H2_CS( double energy_ryd /* photon energy in ryd */)
46 {
47 
48  double energy_keV; /* keV */
49  double cross_section; /* barns */
50  double x; /* x = E/15.4 */
51  double xsqrt , x15 , x2;
52  double energy = energy_ryd * EVRYD;
53 
54  DEBUG_ENTRY( "Yan_H2_CS()" );
55 
56  /* energy relative to threshold */
57  x = energy / 15.4;
58  energy_keV = energy/1000.0;
59  xsqrt = sqrt(x);
60  x15 = x * xsqrt;
61  x2 = x*x;
62 
63  if( energy < 15.4 )
64  {
65  cross_section = 0.;
66  }
67 
68  else if(energy >= 15.4 && energy < 18.)
69  {
70  cross_section = 1e7 * (1 - 197.448/xsqrt + 438.823/x - 260.481/x15 + 17.915/x2);
71  /* this equation is obviously negative for x = 1 */
72  cross_section = MAX2( 0. , cross_section );
73  }
74 
75  else if(energy >= 18. && energy <= 30.)
76  {
77  cross_section = (-145.528 +351.394*xsqrt - 274.294*x + 74.320*x15)/pow(energy_keV,3.5);
78  }
79 
80  else if(energy > 30. && energy <= 85.)
81  {
82  cross_section = (65.304 - 91.762*xsqrt + 51.778*x - 9.364*x15)/pow(energy_keV,3.5);
83  }
84 
85  /* the high-energy tail */
86  else
87  {
88  cross_section = 45.57*(1 - 2.003/xsqrt - 4.806/x + 50.577/x15 - 171.044/x2 + 231.608/(xsqrt*x2) - 81.885/(x*x2))/pow(energy_keV,3.5);
89  }
90 
91  return( cross_section * 1e-24 );
92 }
93 
94 /*OpacityCreate1Element generate opacities for entire element by calling t_ADfA::Inst().phfit */
95 STATIC void OpacityCreate1Element(long int nelem);
96 
97 /*opacity_more_memory allocate more memory for opacity stack */
98 STATIC void opacity_more_memory(void);
99 
100 /*hmiopc derive total H- H minus opacity */
101 STATIC double hmiopc(double freq);
102 
103 /*rayleh compute Rayleigh scattering cross section for Lya */
104 STATIC double rayleh(double ener);
105 
106 /*Opacity_iso_photo_cs returns photoionization cross section for isoelectronic sequences */
107 STATIC double Opacity_iso_photo_cs( realnum energy , long ipISO , long nelem , long n );
108 
109 /*OpacityCreateReilMan generate photoionization cross sections from Reilman and Manson points */
110 STATIC void OpacityCreateReilMan(long int low,
111  long int ihi,
112  const realnum cross[],
113  long int ncross,
114  long int *ipop,
115  const char *chLabl);
116 
117 static bool lgRealloc = false;
118 
119 /*OpacityCreatePowerLaw generate array of cross sections using a simple power law fit */
121  /* lower energy limit on continuum mesh */
122  long int ilo,
123  /* upper energy limit on continuum mesh */
124  long int ihi,
125  /* threshold cross section */
126  double cross,
127  /* power law index */
128  double s,
129  /* pointer to opacity offset where this starts */
130  long int *ip);
131 
132 /*ofit compute cross sections for all shells of atomic oxygen */
133 STATIC double ofit(double e,
134  realnum opart[]);
135 
136 /*OpacityValenceRescale routine to rescale non-OP valence shell cross sections for atom */
138  /* element number on C scale */
139  long int nelem ,
140  /* scale factor, must be >= 0. */
141  double scale )
142 {
143 
144  long int ion , nshell , low , ihi , ipop , ip;
145 
146  DEBUG_ENTRY( "OpacityValenceRescale()" );
147 
148  /* return if element is not turned on
149  * >>chng 05 oct 19, this had not been done, so low in the opacity offset below was
150  * not set, and opacity index was negative - only problem when K turned off */
151  if( !dense.lgElmtOn[nelem] )
152  {
153  return;
154  }
155 
156  /* scale better be >= 0. */
157  ASSERT( scale >= 0. );
158 
159  ion = 0;
160  /* this is valence shell on C scale */
161  nshell = Heavy.nsShells[nelem][ion] - 1;
162 
163  /* set lower and upper limits to this range */
164  low = opac.ipElement[nelem][ion][nshell][0];
165  ihi = opac.ipElement[nelem][ion][nshell][1];
166  ipop = opac.ipElement[nelem][ion][nshell][2];
167 
168  /* loop over energy range of this shell */
169  for( ip=low-1; ip < ihi; ip++ )
170  {
171  opac.OpacStack[ip-low+ipop] *= scale;
172  }
173  return;
174 }
175 
177 {
178  long int i,
179  ipISO ,
180  n ,
181  need ,
182  nelem;
183 
184  realnum opart[7];
185 
186  double crs,
187  dx,
188  eps,
189  thom,
190  thres,
191  x;
192 
193  /* >>chng 02 may 29, change to lgOpacMalloced */
194  /* remember whether opacities have ever been evaluated in this coreload
195  static bool lgOpEval = false; */
196 
197  /* fits to cross section for photo dist of H_2^+ */
198  static const realnum csh2p[NCSH2P]={6.75f,0.24f,8.68f,2.5f,10.54f,7.1f,12.46f,
199  6.0f,14.28f,2.7f};
200 
201  DEBUG_ENTRY( "OpacityCreateAll()" );
202 
203  /* H2+ h2plus h2+ */
204 
205  /* make and print dust opacities
206  * fill in dstab and dstsc, totals, zero if no dust
207  * may be different if different grains are turned on */
208  GrainsInit();
209 
210  /* flag lgOpacMalloced says whether opacity stack has been generated
211  * only do this one time per core load */
212  /* >>chng 02 may 29, from lgOpEval to lgOpacMalloced */
213  if( lgOpacMalloced )
214  {
215  /* this is not the first time code called */
216  if( trace.lgTrace )
217  {
218  fprintf( ioQQQ, " OpacityCreateAll called but NOT evaluated since already done.\n" );
219  }
220  return;
221  }
222 
223  lgOpacMalloced = true;
224 
225  /* create the space for the opacity stack */
226  opac.OpacStack = (double*)MALLOC((size_t)ndimOpacityStack*sizeof(double));
227 
228  if( trace.lgTrace )
229  {
230  fprintf( ioQQQ, " OpacityCreateAll called, evaluating.\n" );
231  }
232 
233  /* zero out opac since this array sometimes addressed before OpacityAddTotal called */
234  for( i=0; i < rfield.nupper; i++ )
235  {
236  opac.opacity_abs[i] = 0.;
237  }
238 
239  /* nOpacTot is number of opacity cells in OpacStack filled so far by opacity generating routines */
240  opac.nOpacTot = 0;
241 
242  /* photoionization of h, he-like iso electronic sequences */
243  for( ipISO=ipH_LIKE; ipISO<=ipHE_LIKE; ++ipISO )
244  {
245  for( nelem=ipISO; nelem < LIMELM; nelem++ )
246  {
247  if( dense.lgElmtOn[nelem] )
248  {
249  long int nupper;
250 
251  /* this is the opacity offset in the general purpose pointer array */
252  /* indices are type, shell. ion, element, so this is the inner shell,
253  * NB - this works for H and He, but the second index should be 1 for Li */
254  opac.ipElement[nelem][nelem-ipISO][0][2] = opac.nOpacTot + 1;
255 
256  /* gound state goes to high-energy limit of code,
257  * but excited states only go up to edge for ground */
258  nupper = rfield.nupper;
259  for( n=0; n < iso.numLevels_max[ipISO][nelem]; n++ )
260  {
261  /* this is array index to the opacity offset */
262  iso.ipOpac[ipISO][nelem][n] = opac.nOpacTot + 1;
263 
264  /* first make sure that first energy point is at least near the limit */
265  /* >>chng 01 sep 23, increased factor form 0.98 to 0.94, needed since cells now go
266  * so far into radio, where resolution is poor */
267  ASSERT( rfield.AnuOrg[iso.ipIsoLevNIonCon[ipISO][nelem][n]-1] > 0.94f *
268  iso.xIsoLevNIonRyd[ipISO][nelem][n] );
269 
270  /* number of cells we will need to do this level */
271  need = nupper - iso.ipIsoLevNIonCon[ipISO][nelem][n] + 1;
272  ASSERT( need > 0 );
273 
274  if( opac.nOpacTot + need > ndimOpacityStack )
276 
277  for( i=iso.ipIsoLevNIonCon[ipISO][nelem][n]-1; i < nupper; i++ )
278  {
279  opac.OpacStack[i-iso.ipIsoLevNIonCon[ipISO][nelem][n]+iso.ipOpac[ipISO][nelem][n]] =
280  Opacity_iso_photo_cs( rfield.AnuOrg[i] , ipISO , nelem , n );
281  }
282 
283  opac.nOpacTot += need;
284  /* for all excited levels high-energy limit is edge for ground state */
285  nupper = iso.ipIsoLevNIonCon[ipISO][nelem][0];
286  }
287  }
288  }
289  }
290 
291  /* This check will get us through Klein-Nishina below. */
292  /* >>chng 02 may 08, by Ryan. Added this and other checks for allotted memory. */
295 
296  /* Lyman alpha damping wings - Rayleigh scattering */
297  opac.ipRayScat = opac.nOpacTot + 1;
298  for( i=0; i < iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]; i++ )
299  {
301  }
303 
304  /* ==============================================================
305  * this block of code defines the electron scattering cross section
306  * for all energies */
307 
308  /* assume Thomson scattering up to ipCKshell, 20.6 Ryd=0.3 keV */
309  thom = 6.65e-25;
310  opac.iopcom = opac.nOpacTot + 1;
311  for( i=0; i < opac.ipCKshell; i++ )
312  {
313  opac.OpacStack[i-1+opac.iopcom] = thom;
314  /*fprintf(ioQQQ,"%.3e\t%.3e\n",
315  rfield.AnuOrg[i]*EVRYD , opac.OpacStack[i-1+opac.iopcom] );*/
316  }
317 
318  /* Klein-Nishina from eqn 7.5,
319  * >>refer Klein-Nishina cs Rybicki and Lightman */
320  for( i=opac.ipCKshell; i < rfield.nupper; i++ )
321  {
322  dx = rfield.AnuOrg[i]/3.7573e4;
323 
324  opac.OpacStack[i-1+opac.iopcom] = thom*3.e0/4.e0*((1.e0 +
325  dx)/POW3(dx)*(2.e0*dx*(1.e0 + dx)/(1.e0 + 2.e0*dx) - log(1.e0+
326  2.e0*dx)) + 1.e0/2.e0/dx*log(1.e0+2.e0*dx) - (1.e0 + 3.e0*
327  dx)/POW3(1.e0 + 2.e0*dx));
328  /*fprintf(ioQQQ,"%.3e\t%.3e\n",
329  rfield.AnuOrg[i]*EVRYD , opac.OpacStack[i-1+opac.iopcom] );*/
330  }
331  opac.nOpacTot += rfield.nupper - 1 + 1;
332 
333  /* ============================================================== */
334 
335  /* This check will get us through "H- hminus H minus bound-free opacity" below. */
336  /* >>chng 02 may 08, by Ryan. Added this and other checks for allotted memory. */
339 
340  /* pair production */
341  opac.ioppr = opac.nOpacTot + 1;
342  for( i=opac.ippr-1; i < rfield.nupper; i++ )
343  {
344  /* pair production heating rate for unscreened H + He
345  * fit to figure 41 of Experimental Nuclear Physics,
346  * Vol1, E.Segre, ed */
347 
348  x = rfield.AnuOrg[i]/7.512e4*2.;
349 
350  opac.OpacStack[i-opac.ippr+opac.ioppr] = 5.793e-28*
351  POW2((-0.46737118 + x*(0.349255416 + x*0.002179893))/(1. +
352  x*(0.130471301 + x*0.000524906)));
353  /*fprintf(ioQQQ,"%.3e\t%.3e\n",
354  rfield.AnuOrg[i]*EVRYD , opac.OpacStack[i-opac.ippr+opac.ioppr] );*/
355  }
356  opac.nOpacTot += rfield.nupper - opac.ippr + 1;
357 
358  /* brems (free-free) opacity */
359  opac.ipBrems = opac.nOpacTot + 1;
360 
361  for( i=0; i < rfield.nupper; i++ )
362  {
363  /* missing factor of 1E-20 to avoid underflow
364  * free free opacity needs g(ff)*(1-exp(hn/kT))/SQRT(T)*1E-20 */
365  opac.OpacStack[i-1+opac.ipBrems] =
366  /*(realnum)(1.03680e-18/POW3(rfield.AnuOrg[i]));*/
367  /* >>chng 00 jun 05, descale by 1e10 so that underflow at high-energy
368  * end does not occur */
369  1.03680e-8/POW3(rfield.AnuOrg[i]);
370  }
371  opac.nOpacTot += rfield.nupper - 1 + 1;
372 
373  opac.iphmra = opac.nOpacTot + 1;
374  for( i=0; i < rfield.nupper; i++ )
375  {
376  /* following is ratio of h minus to neut h bremss opacity */
377  opac.OpacStack[i-1+opac.iphmra] = 0.1175*rfield.anusqr[i];
378  }
379  opac.nOpacTot += rfield.nupper - 1 + 1;
380 
381  opac.iphmop = opac.nOpacTot + 1;
382  for( i=hmi.iphmin-1; i < iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0]; i++ )
383  {
384  /* H- hminus H minus bound-free opacity */
386  hmiopc(rfield.AnuOrg[i]);
387  }
389 
390  /* ============================================================== */
391 
392  /* This check will get us through "H2 photoionization cross section" below. */
393  /* >>chng 07 oct 10, by Ryan. Added this check for allotted memory. */
396 
398  for( i=opac.ipH2_photo_thresh-1; i < rfield.nupper; i++ )
399  {
400  /* >>chng 05 nov 24, add H2 photoionization cross section */
402  Yan_H2_CS(rfield.AnuOrg[i]);
403  /*fprintf(ioQQQ,"DEBUG H2 photo\t%.3e\t%.3e\t%.3e\n",
404  rfield.AnuOrg[i],
405  opac.OpacStack[i-opac.ipH2_photo_thresh + opac.ipH2_photo_opac_offset] ,
406  Opacity_iso_photo_cs( rfield.AnuOrg[i] , 0 , 0 , 0 )*2. );*/
407  }
409 
410  /* H2+ H2P h2plus photoabsorption
411  * cs from
412  * >>refer H2+ photodissoc Buckingham, R.A., Reid, S., & Spence, R. 1952, MNRAS 112, 382, 0 K temp */
414  "H2+ ");
415 
416  /* This check will get us through "HeI singlets neutral helium ground" below. */
417  /* >>chng 02 may 08, by Ryan. Added this and other checks for allotted memory. */
420 
421  /* HeI singlets neutral helium ground */
422  opac.iophe1[0] = opac.nOpacTot + 1;
423  opac.ipElement[ipHELIUM][0][0][2] = opac.iophe1[0];
424  for( i=iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0]-1; i < rfield.nupper; i++ )
425  {
426  crs = t_ADfA::Inst().phfit(2,2,1,rfield.AnuOrg[i]*EVRYD);
428  crs*1e-18;
429  }
431 
432  /* these are opacity offset points that would be defined in OpacityCreate1Element,
433  * but this routine will not be called for H and He
434  * generate all heavy element opacities, everything heavier than He,
435  * nelem is on the C scale, so Li is 2 */
436  /*>>chng 99 jan 27, do not reevaluate hydrogenic opacity below */
437  for( nelem=2; nelem < LIMELM; nelem++ )
438  {
439  if( dense.lgElmtOn[nelem] )
440  {
441  OpacityCreate1Element(nelem);
442  }
443  }
444 
445  /* option to rescale atoms of some elements that were not done by opacity project
446  * the valence shell - two arguments - element number on C scale, and scale factor */
447  /*>>chng 05 sep 26, fudge factor to get atomic K fraction along well defined line of sight
448  * to be observed value - this is ratio of cross sections, actual value is very uncertain since
449  * differences betweeen Verner & opacity project are huge */
451 
452  /* now add on some special cases, where exicted states, etc, come in */
453  /* Nitrogen
454  * >>refer n1 photo Henry, R., ApJ 161, 1153.
455  * photoionization of excited level of N+ */
456  OpacityCreatePowerLaw(opac.in1[0],opac.in1[1],9e-18,1.75,&opac.in1[2]);
457 
458  /* atomic Oxygen
459  * only do this if 1996 Verner results are used */
460  if( dense.lgElmtOn[ipOXYGEN] && t_ADfA::Inst().get_version() == PHFIT96 )
461  {
462  /* This check will get us through this loop. */
463  /* >>chng 02 may 08, by Ryan. Added this and other checks for allotted memory. */
464  if( opac.nOpacTot + opac.ipElement[ipOXYGEN][0][2][1] - opac.ipElement[ipOXYGEN][0][2][0] + 1 > ndimOpacityStack )
466 
467  /* integrate over energy range of the valence shell of atomic oxygen*/
468  for( i=opac.ipElement[ipOXYGEN][0][2][0]-1; i < opac.ipElement[ipOXYGEN][0][2][1]; i++ )
469  {
470  /* call special routine to evaluate partial cross section for OI shells */
471  eps = rfield.AnuOrg[i]*EVRYD;
472  crs = ofit(eps,opart);
473 
474  /* this will be total cs of all processes leaving shell 3 */
475  crs = opart[0];
476  for( n=1; n < 6; n++ )
477  {
478  /* add up table of cross sections */
479  crs += opart[n];
480  }
481  /* convert to cgs and overwrite cross sections set by OpacityCreate1Element */
482  crs *= 1e-18;
483  opac.OpacStack[i-opac.ipElement[ipOXYGEN][0][2][0]+opac.ipElement[ipOXYGEN][0][2][2]] = crs;
484  }
485  /* >>chng 02 may 09 - this was a significant error */
486  /* >>chng 02 may 08, by Ryan. This loop did not update total slots filled. */
487  opac.nOpacTot += opac.ipElement[ipOXYGEN][0][2][1] - opac.ipElement[ipOXYGEN][0][2][0] + 1;
488  }
489 
490  /* Henry nubmers for 1S excit state of OI, OP data very sparse */
491  OpacityCreatePowerLaw(opac.ipo1exc[0],opac.ipo1exc[1],4.64e-18,0.,&opac.ipo1exc[2]);
492 
493  /* photoionization of excited level of O2+ 1D making 5007
494  * fit to TopBase Opacity Project cs */
495  OpacityCreatePowerLaw(opac.ipo3exc[0],opac.ipo3exc[1],3.8e-18,0.,&opac.ipo3exc[2]);
496 
497  /* photoionization of excited level of O2+ 1S making 4363 */
498  OpacityCreatePowerLaw(opac.ipo3exc3[0],opac.ipo3exc3[1],5.5e-18,0.01,
499  &opac.ipo3exc3[2]);
500 
501  /* This check will get us through the next two steps. */
502  /* >>chng 02 may 08, by Ryan. Added this and other checks for allotted memory. */
506 
507  /* photoionization to excited states of O+ */
508  opac.iopo2d = opac.nOpacTot + 1;
509  thres = rfield.AnuOrg[oxy.i2d-1];
510  for( i=oxy.i2d-1; i < iso.ipIsoLevNIonCon[ipH_LIKE][ipHELIUM][ipH1s]; i++ )
511  {
512  crs = 3.85e-18*(4.4*pow(rfield.AnuOrg[i]/thres,-1.5) - 3.38*
513  pow(rfield.AnuOrg[i]/thres,-2.5));
514 
515  opac.OpacStack[i-oxy.i2d+opac.iopo2d] = crs;
516  }
518 
519  /* magnesium
520  * photoionization of excited level of Mg+
521  * fit to opacity project data Dima got */
522  opac.ipOpMgEx = opac.nOpacTot + 1;
523  for( i=opac.ipmgex-1; i < iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]; i++ )
524  {
526  (0.2602325880970085 +
527  445.8558249365131*exp(-rfield.AnuOrg[i]/0.1009243952792674))*
528  1e-18;
529  }
531 
533 
534  /* Calcium
535  * excited states of Ca+ */
537 
538  if( trace.lgTrace )
539  {
540  fprintf( ioQQQ,
541  " OpacityCreateAll return OK, number of opacity cells used in OPSC= %ld and OPSV is dimensioned %ld\n",
543  }
544 
545  /* option to compile opacities into file for later use
546  * this is executed if the 'compile opacities' command is entered */
547  if( opac.lgCompileOpac )
548  {
549  fprintf( ioQQQ, "The COMPILE OPACITIES command is currently not supported\n" );
550  cdEXIT(EXIT_FAILURE);
551  }
552 
553  if( lgRealloc )
554  fprintf(ioQQQ,
555  " Please consider revising ndimOpacityStack in OpacityCreateAll, a total of %li cells were needed.\n\n" , opac.nOpacTot);
556  return;
557 }
558 /*OpacityCreatePowerLaw generate array of cross sections using a simple power law fit */
560  /* lower energy limit on continuum mesh */
561  long int ilo,
562  /* upper energy limit on continuum mesh */
563  long int ihi,
564  /* threshold cross section */
565  double cross,
566  /* power law index */
567  double s,
568  /* pointer to opacity offset where this starts */
569  long int *ip)
570 {
571  long int i;
572  double thres;
573 
574  DEBUG_ENTRY( "OpacityCreatePowerLaw()" );
575 
576  /* non-positive cross section is unphysical */
577  ASSERT( cross > 0. );
578 
579  /* place in the opacity stack where we will stuff cross sections */
580  *ip = opac.nOpacTot + 1;
581  ASSERT( *ip > 0 );
582  ASSERT( ilo > 0 );
583  thres = rfield.anu[ilo-1];
584 
585  if( opac.nOpacTot + ihi - ilo + 1 > ndimOpacityStack )
587 
588  for( i=ilo-1; i < ihi; i++ )
589  {
590  opac.OpacStack[i-ilo+*ip] = cross*pow(rfield.anu[i]/thres,-s);
591  }
592 
593  opac.nOpacTot += ihi - ilo + 1;
594  return;
595 }
596 
597 /*OpacityCreateReilMan generate photoionization cross sections from Reilman and Manson points */
598 STATIC void OpacityCreateReilMan(long int low,
599  long int ihi,
600  const realnum cross[],
601  long int ncross,
602  long int *ipop,
603  const char *chLabl)
604 {
605  long int i,
606  ics,
607  j,
608  ncr;
609 
610  const int NOP = 100;
611  realnum cs[NOP],
612  en[NOP],
613  slope;
614 
615  DEBUG_ENTRY( "OpacityCreateReilMan()" );
616 
617  /* this is the opacity entering routine designed for
618  * the Reilman and Manson tables. It works with incident
619  * photon energy (entered in eV) and cross sections in megabarns
620  * */
621  *ipop = opac.nOpacTot + 1;
622  ASSERT( *ipop > 0 );
623 
624  ncr = ncross/2;
625  if( ncr > NOP )
626  {
627  fprintf( ioQQQ, " Too many opacities were entered into OpacityCreateReilMan. Increase the value of NOP.\n" );
628  fprintf( ioQQQ, " chLabl was %4.4s\n", chLabl );
629  cdEXIT(EXIT_FAILURE);
630  }
631 
632  /* the array CROSS has ordered pairs of elements.
633  * the first is the energy in eV (not Ryd)
634  * and the second is the cross section in megabarns */
635  for( i=0; i < ncr; i++ )
636  {
637  en[i] = cross[i*2]/13.6f;
638  cs[i] = cross[(i+1)*2-1]*1e-18f;
639  }
640 
641  ASSERT( low>0 );
642  if( en[0] > rfield.anu[low-1] )
643  {
644  fprintf( ioQQQ,
645  " OpacityCreateReilMan: The entered opacity energy bandwidth is not large enough (low fail).\n" );
646  fprintf( ioQQQ,
647  " The desired energy (Ryd) was%12.5eeV and the lowest entered in the array was%12.5e eV\n",
648  rfield.anu[low-1]*EVRYD, en[0]*EVRYD );
649 
650  fprintf( ioQQQ, " chLabl was %4.4s\n", chLabl );
651  fprintf( ioQQQ, " The original energy (eV) and cross section (mb) arrays follow:\n" );
652  fprintf( ioQQQ, " " );
653 
654  for( i=0; i < ncross; i++ )
655  {
656  fprintf( ioQQQ, "%11.4e", cross[i] );
657  }
658 
659  fprintf( ioQQQ, "\n" );
660  cdEXIT(EXIT_FAILURE);
661  }
662 
663  slope = (cs[1] - cs[0])/(en[1] - en[0]);
664  ics = 1;
665 
666  if( opac.nOpacTot + ihi - low + 1 > ndimOpacityStack )
668 
669  /* now fill in the opacities using linear interpolation */
670  for( i=low-1; i < ihi; i++ )
671  {
672  if( rfield.anu[i] > en[ics-1] && rfield.anu[i] <= en[ics] )
673  {
674  opac.OpacStack[i-low+*ipop] = cs[ics-1] + slope*(rfield.anu[i] -
675  en[ics-1]);
676  }
677 
678  else
679  {
680  ics += 1;
681  if( ics + 1 > ncr )
682  {
683  fprintf( ioQQQ, " OpacityCreateReilMan: The entered opacity energy bandwidth is not large enough (high fail).\n" );
684  fprintf( ioQQQ, " The entered energy was %10.2eeV and the highest in the array was %10.2eeV\n",
685  rfield.anu[i]*13.6, en[ncr-1]*13.6 );
686  fprintf( ioQQQ, " chLabl was %4.4s\n", chLabl
687  );
688  fprintf( ioQQQ, " The lowest energy enterd in the array was%10.2e eV\n",
689  en[0]*13.65 );
690  fprintf( ioQQQ, " The highest energy ever needed would be%10.2eeV\n",
691  rfield.anu[ihi-1]*13.6 );
692  fprintf( ioQQQ, " The lowest energy needed was%10.2eeV\n",
693  rfield.anu[low-1]*13.6 );
694  cdEXIT(EXIT_FAILURE);
695  }
696 
697  slope = (cs[ics] - cs[ics-1])/(en[ics] - en[ics-1]);
698  if( rfield.anu[i] > en[ics-1] && rfield.anu[i] <= en[ics] )
699  {
700  opac.OpacStack[i-low+*ipop] = cs[ics-1] + slope*(rfield.anu[i] -
701  en[ics-1]);
702  }
703  else
704  {
705  ASSERT( i > 0);
706  fprintf( ioQQQ, " Internal logical error in OpacityCreateReilMan.\n" );
707  fprintf( ioQQQ, " The desired energy (%10.2eeV), I=%5ld, is not within the next energy bound%10.2e%10.2e\n",
708  rfield.anu[i]*13.6, i, en[ics-1], en[ics] );
709 
710  fprintf( ioQQQ, " The previous energy (eV) was%10.2e\n",
711  rfield.anu[i-1]*13.6 );
712 
713  fprintf( ioQQQ, " Here comes the energy array. ICS=%4ld\n",
714  ics );
715 
716  for( j=0; j < ncr; j++ )
717  {
718  fprintf( ioQQQ, "%10.2e", en[j] );
719  }
720  fprintf( ioQQQ, "\n" );
721 
722  fprintf( ioQQQ, " chLabl was %4.4s\n", chLabl );
723  cdEXIT(EXIT_FAILURE);
724  }
725  }
726  }
727  /* >>chng 02 may 09, this was a significant logcal error */
728  /* >>chng 02 may 08, by Ryan. This routine did not update the total slots filled. */
729  opac.nOpacTot += ihi - low + 1;
730  return;
731 }
732 
733 
734 /*ofit compute cross sections for all shells of atomic oxygen */
735 STATIC double ofit(double e,
736  realnum opart[])
737 {
738  double otot,
739  q,
740  x;
741 
742  static const double y[7][5] = {
743  {8.915,3995.,3.242,10.44,0.0},
744  {11.31,1498.,5.27,7.319,0.0},
745  {10.5,1.059e05,1.263,13.04,0.0},
746  {19.49,48.47,8.806,5.983,0.0},
747  {50.,4.244e04,0.1913,7.012,4.454e-02},
748  {110.5,0.1588,148.3,-3.38,3.589e-02},
749  {177.4,32.37,381.2,1.083,0.0}
750  };
751  static const double eth[7]={13.62,16.94,18.79,28.48,50.,110.5,538.};
752  static const long l[7]={1,1,1,0,1,1,0};
753 
754  DEBUG_ENTRY( "ofit()" );
755  /*compute cross sections for all shells of atomic oxygen
756  * Photoionization of OI
757  * Input parameter: e - photon energy, eV
758  * Output parameters: otot - total photoionization cross section, Mb
759  * opart(1) - 2p-shell photoionization, goes to 4So
760  * opart(2) - 2p-shell photoionization, goes to 2Do
761  * opart(3) - 2p-shell photoionization, goes to 2Po
762  * opart(4) - 2s-shell photoionization
763  * opart(5) - double photoionization, goes to O++
764  * opart(6) - triple photoionization, goes to O+++
765  * opart(7) - 1s-shell photoionization */
766 
767  otot = 0.0;
768 
769  for( int i=0; i < 7; i++ )
770  {
771  opart[i] = 0.0;
772  }
773 
774  for( int i=0; i < 7; i++ )
775  {
776  if( e >= eth[i] )
777  {
778  q = 5.5 - 0.5*y[i][3] + l[i];
779 
780  x = e/y[i][0];
781 
782  opart[i] = (realnum)(y[i][1]*(POW2(x - 1.0) + POW2(y[i][4]))/
783  pow(x,q)/pow(1.0 + sqrt(x/y[i][2]),y[i][3]));
784 
785  otot += opart[i];
786  }
787  }
788  return(otot);
789 }
790 
791 /******************************************************************************/
792 
793 /*OpacityCreate1Element generate ionic subshell opacities by calling t_ADfA::Inst().phfit */
795  /* atomic number on the C scale, lowest ever called will be Li=2 */
796  long int nelem)
797 {
798  long int ihi,
799  ip,
800  ipop,
801  limit,
802  low,
803  need,
804  nelec,
805  ion,
806  nshell;
807  double cs;
808  double energy;
809 
810  DEBUG_ENTRY( "OpacityCreate1Element()" );
811 
812  /* confirm range of validity of atomic number, Li=2 should be the lightest */
813  ASSERT( nelem >= 2 );
814  ASSERT( nelem < LIMELM );
815 
816  /*>>chng 99 jan 27, no longer redo hydrogenic opacity here */
817  /* for( ion=0; ion <= nelem; ion++ )*/
818  for( ion=0; ion < nelem; ion++ )
819  {
820 
821  /* will be used for a sanity check on number of hits in a cell*/
822  for( ip=0; ip < rfield.nupper; ip++ )
823  {
824  opac.opacity_abs[ip] = 0.;
825  }
826 
827  /* number of bound electrons */
828  nelec = nelem+1 - ion;
829 
830  /* loop over all shells, from innermost K shell to valence */
831  for( nshell=0; nshell < Heavy.nsShells[nelem][ion]; nshell++ )
832  {
833  /* this is array index for start of this shell within large opacity stack */
834  opac.ipElement[nelem][ion][nshell][2] = opac.nOpacTot + 1;
835 
836  /* this is continuum upper limit to array index for energy range of this shell */
837  limit = opac.ipElement[nelem][ion][nshell][1];
838 
839  /* this is number of cells in continuum needed to store opacity */
840  need = limit - opac.ipElement[nelem][ion][nshell][0] + 1;
841 
842  /* check that opac will have enough frequeny cells */
843  if( opac.nOpacTot + need > ndimOpacityStack )
845 
846  /* set lower and upper limits to this range */
847  low = opac.ipElement[nelem][ion][nshell][0];
848  ihi = opac.ipElement[nelem][ion][nshell][1];
849  ipop = opac.ipElement[nelem][ion][nshell][2];
850 
851  /* make sure indices are within correct bounds,
852  * mainly check on logic for detecting missing shells */
853  ASSERT( low <= ihi || low<5 );
854 
855  /* loop over energy range of this shell */
856  for( ip=low-1; ip < ihi; ip++ )
857  {
858  /* photo energy MAX so that we never eval below threshold */
859  energy = MAX2(rfield.AnuOrg[ip]*EVRYD ,
860  t_ADfA::Inst().ph1(nshell,nelec-1,nelem,0));
861 
862  /* the cross section in mega barns */
863  cs = t_ADfA::Inst().phfit(nelem+1,nelec,nshell+1,energy);
864  /* cannot assert that cs is positive since, at edge of shell,
865  * energy might be slightly below threshold and hence zero,
866  * due to finite size of continuum bins */
867  opac.OpacStack[ip-low+ipop] = cs*1e-18;
868 
869  /* add this to total opacity, which we will confirm to be greater than zero below */
870  opac.opacity_abs[ip] += cs;
871  }
872 
873  opac.nOpacTot += ihi - low + 1;
874 
875  /* punch pointers option */
876  if( punch.lgPunPoint )
877  {
878  fprintf( punch.ipPoint, "%3ld%3ld%3ld%10.2e%10.2e%10.2e%10.2e\n",
879  nelem, ion, nshell, rfield.anu[low-1], rfield.anu[ihi-1],
880  opac.OpacStack[ipop-1], opac.OpacStack[ihi-low+ipop-1] );
881  }
882  }
883 
884  ASSERT( Heavy.nsShells[nelem][ion] >= 1 );
885  /*confirm that total opacity is greater than zero */
886  for( ip=opac.ipElement[nelem][ion][Heavy.nsShells[nelem][ion]-1][0]-1;
887  ip < continuum.KshellLimit; ip++ )
888  {
889  ASSERT( opac.opacity_abs[ip] > 0. );
890  }
891 
892  }
893  return;
894 }
895 
896 /*opacity_more_memory allocate more memory for opacity stack */
898 {
899 
900  DEBUG_ENTRY( "opacity_more_memory()" );
901 
902  /* double size */
903  ndimOpacityStack *= 2;
904  opac.OpacStack = (double *)REALLOC( opac.OpacStack , (size_t)ndimOpacityStack*sizeof(double) );
905  fprintf( ioQQQ, " NOTE OpacityCreate1Element needed more opacity cells than ndimOpacityStack, please consider increasing it.\n" );
906  fprintf( ioQQQ, " NOTE OpacityCreate1Element doubled memory allocation to %li.\n",ndimOpacityStack );
907  lgRealloc = true;
908  return;
909 }
910 
911 /*Opacity_iso_photo_cs returns photoionization cross section for isoelectronic sequences */
913  /* photon energy ryd */
914  realnum energy ,
915  /* iso sequence */
916  long ipISO ,
917  /* charge, 0 for H */
918  long nelem ,
919  /* n - meaning depends on iso */
920  long n )
921 {
922  /* >>chng 01 dec 23, from float to double */
923  double thres/*,ejected_electron_energy*/;
924  double crs=-DBL_MAX;
925 
926  DEBUG_ENTRY( "Opacity_iso_photo_cs()" );
927 
928  if( ipISO==ipH_LIKE )
929  {
930  double photon;
931 
932  if( n==0 )
933  {
934  /* this is the ground state, use Dima's routine, which works in eV
935  * and returns megabarns */
936  thres = MAX2(energy*(realnum)EVRYD , t_ADfA::Inst().ph1(0,0,nelem,0));
937  crs = t_ADfA::Inst().phfit(nelem+1,1,1,thres)* 1e-18;
938  /* make sure cross section is reasonable */
939  ASSERT( crs > 0. && crs < 1e-10 );
940  }
941  else if( n < iso.numLevels_max[ipISO][nelem] - iso.nCollapsed_max[ipISO][nelem] )
942  {
943  /* photon energy relative to threshold */
944  photon = MAX2( energy/iso.xIsoLevNIonRyd[ipISO][nelem][n], 1. + FLT_EPSILON*2. );
945 
946  crs = H_photo_cs( photon , N_(n), L_(n), nelem+1 );
947  /* make sure cross section is reasonable */
948  ASSERT( crs > 0. && crs < 1e-10 );
949  }
950  else if( N_(n) <= NHYDRO_MAX_LEVEL )
951  {
952  /* for first cell, depending on the current resolution of the energy mesh,
953  * the center of the first cell can be below the ionization limit of the
954  * level. do not let the energy fall below this limit */
955  /* This will make sure that we don't call epsilon below threshold,
956  * the factor 1.001 was chosen so that t_ADfA::Inst().hpfit, which works
957  * in terms of Dima's Rydberg constant, is not tripped below threshold */
958  thres = MAX2( energy , iso.xIsoLevNIonRyd[ipISO][nelem][n]*1.001f );
959 
960  crs = t_ADfA::Inst().hpfit(nelem+1,N_(n),thres*EVRYD);
961  /* make sure cross section is reasonable */
962  ASSERT( crs > 0. && crs < 1e-10 );
963  }
964  else
965  {
966  /* photon energy relative to threshold */
967  photon = MAX2( energy/iso.xIsoLevNIonRyd[ipISO][nelem][n], 1. + FLT_EPSILON*2. );
968 
969  /* cross section for collapsed level should be
970  * roughly equal to cross-section for yrast level,
971  * so third parameter is n - 1. */
972  crs = H_photo_cs( photon , N_(n), N_(n)-1, nelem+1 );
973 
974  /* make sure cross section is reasonable */
975  ASSERT( crs > 0. && crs < 1e-10 );
976  }
977  }
978  else if( ipISO==ipHE_LIKE )
979  {
980  thres = MAX2( energy , iso.xIsoLevNIonRyd[ipISO][nelem][n]);
981  /* this would be a collapsed level */
982  if( n >= iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem] )
983  {
984  long int nup = iso.n_HighestResolved_max[ipHE_LIKE][nelem] + n + 1 -
986 
987  /* this is a collapsed level - this is hydrogenic routine and
988  * first he-like energy may not agree exactly with threshold for H */
989  crs = t_ADfA::Inst().hpfit(nelem,nup ,thres*EVRYD);
990  /* make sure cross section is reasonable if away from threshold */
991  ASSERT(
992  (energy < iso.xIsoLevNIonRyd[ipISO][nelem][n]*1.02) ||
993  (crs > 0. && crs < 1e-10) );
994  }
995  else
996  {
997 
998  /* He_cross_section returns cross section (cm^-2),
999  * given EgammaRyd, the photon energy in Ryd,
1000  * ipLevel, the index of the level, 0 is ground, 3 within 2 3P,
1001  * nelem is charge, equal to 1 for Helium,
1002  * this is a wrapper for cross_section */
1003  crs = He_cross_section( thres , n , nelem );
1004 
1005  /* make sure cross section is reasonable */
1006  ASSERT( crs > 0. && crs < 1e-10 );
1007  }
1008  }
1009  else
1010  TotalInsanity();
1011  return(crs);
1012 }
1013 
1014 /*hmiopc derive total H- H minus opacity */
1015 static const int NCRS = 33;
1016 
1017 STATIC double hmiopc(double freq)
1018 {
1019  double energy,
1020  hmiopc_v,
1021  x,
1022  y;
1023  static double y2[NCRS];
1024  static double crs[NCRS]={0.,0.124,0.398,0.708,1.054,1.437,1.805,
1025  2.176,2.518,2.842,3.126,3.377,3.580,3.741,3.851,3.913,3.925,
1026  3.887,3.805,3.676,3.511,3.306,3.071,2.810,2.523,2.219,1.898,
1027  1.567,1.233,.912,.629,.39,.19};
1028  static double ener[NCRS]={0.,0.001459,0.003296,0.005256,0.007351,
1029  0.009595,0.01201,0.01460,0.01741,0.02044,0.02375,0.02735,0.03129,
1030  0.03563,0.04043,0.04576,0.05171,0.05841,0.06601,0.07469,0.08470,
1031  0.09638,0.1102,0.1268,0.1470,0.1723,0.2049,0.2483,0.3090,0.4001,
1032  0.5520,0.8557,1.7669};
1033  static bool lgFirst = true;
1034 
1035  DEBUG_ENTRY( "hmiopc()" );
1036 
1037  /* bound free cross section (x10**-17 cm^2) from Doughty et al
1038  * 1966, MNRAS 132, 255; good agreement with Wishart MNRAS 187, 59p. */
1039 
1040  /* photoelectron energy, add 0.05552 to get incoming energy (Ryd) */
1041 
1042 
1043  if( lgFirst )
1044  {
1045  /* set up coefficients for spline */
1046  spline(ener,crs,NCRS,2e31,2e31,y2);
1047  lgFirst = false;
1048  }
1049 
1050  energy = freq - 0.05552;
1051  if( energy < ener[0] || energy > ener[NCRS-1] )
1052  {
1053  hmiopc_v = 0.;
1054  }
1055  else
1056  {
1057  x = energy;
1058  splint(ener,crs,y2,NCRS,x,&y);
1059  hmiopc_v = y*1e-17;
1060  }
1061  return( hmiopc_v );
1062 }
1063 
1064 /*rayleh compute Rayleigh scattering cross section for Lya */
1065 STATIC double rayleh(double ener)
1066 {
1067  double rayleh_v;
1068 
1069  DEBUG_ENTRY( "rayleh()" );
1070 
1072  /* do hydrogen Rayleigh scattering cross sections;
1073  * fits to
1074  *>>refer Ly scattering Gavrila, M., 1967, Physical Review 163, 147
1075  * and Mihalas radiative damping
1076  *
1077  * >>chng 96 aug 15, changed logic to do more terms for each part of
1078  * rayleigh scattering
1079  * if( ener.lt.0.05 ) then
1080  * rayleh = 8.41e-25 * ener**4 * DampOnFac
1081  * */
1082  if( ener < 0.05 )
1083  {
1084  rayleh_v = (8.41e-25*powi(ener,4) + 3.37e-24*powi(ener,6))*
1085  hydro.DampOnFac;
1086  }
1087 
1088  else if( ener < 0.646 )
1089  {
1090  rayleh_v = (8.41e-25*powi(ener,4) + 3.37e-24*powi(ener,6) +
1091  4.71e-22*powi(ener,14))*hydro.DampOnFac;
1092  }
1093 
1094  else if( ener >= 0.646 && ener < 1.0 )
1095  {
1096  rayleh_v = fabs(0.74959-ener);
1097  rayleh_v = 1.788e5/POW2(FR1RYD*MAX2(0.001,rayleh_v));
1098  /* typical energy between Ly-a and Ly-beta */
1099  rayleh_v = MAX2(rayleh_v,1e-24)*hydro.DampOnFac;
1100  }
1101 
1102  else
1103  {
1104  rayleh_v = 0.;
1105  }
1106  return( rayleh_v );
1107 }

Generated for cloudy by doxygen 1.8.3.1