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