cloudy trunk
Loading...
Searching...
No Matches
cont_createpointers.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/*ContCreatePointers set up pointers for lines and continua called by cloudy after input read in
4 * and continuum mesh has been set */
5/*fiddle adjust energy bounds of certain cells so that match ionization edges exactly */
6/*ipShells assign continuum energy pointers to shells for all atoms,
7 * called by ContCreatePointers */
8/*LimitSh sets upper energy limit to subshell integrations */
9/*ContBandsCreate - read set of continuum bands to enter total emission into line stack*/
10#include "cddefines.h"
11#include "physconst.h"
12#include "lines_service.h"
13#include "iso.h"
14#include "secondaries.h"
15#include "taulines.h"
16#include "elementnames.h"
17#include "ionbal.h"
18#include "rt.h"
19#include "opacity.h"
20#include "yield.h"
21#include "dense.h"
22#include "he.h"
23#include "fe.h"
24#include "rfield.h"
25#include "oxy.h"
26#include "atomfeii.h"
27#include "atoms.h"
28#include "trace.h"
29#include "hmi.h"
30#include "heavy.h"
31#include "atmdat.h"
32#include "ipoint.h"
33#include "h2.h"
34#include "continuum.h"
35
36/*LimitSh sets upper energy limit to subshell integrations */
37STATIC long LimitSh(long int ion,
38 long int nshell,
39 long int nelem);
40
41STATIC void ipShells(
42 /* nelem is the atomic number on the C scale, Li is 2 */
43 long int nelem);
44
45/*ContBandsCreate - read set of continuum bands to enter total emission into line*/
47 /* chFile is optional filename, if void then use default bands,
48 * if not void then use file specified,
49 * return value is 0 for success, 1 for failure */
50 const char chFile[] );
51
52/*fiddle adjust energy bounds of certain cells so that match ionization edges exactly */
53STATIC void fiddle(long int ipnt,
54 double exact);
55
57{
58 char chLab[5];
59 /* counter to say whether pointers have ever been evaluated */
60 static int nCalled = 0;
61
62 DEBUG_ENTRY( "ContCreatePointers()" );
63
64 /* create the hydrogen atom for this core load, routine creates space then zeros it out
65 * on first call, on second and later calls it only zeros things out */
66 iso_create();
67
68 /* create internal static variables needed to do the H2 molecule */
69 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
70 (*diatom)->init();
71
72 /* nCalled is local static variable defined 0 when defined.
73 * it is incremented below, so that space only allocated one time per coreload. */
74 if( nCalled > 0 )
75 {
76 if( trace.lgTrace )
77 fprintf( ioQQQ, " ContCreatePointers called, not evaluating.\n" );
78 return;
79 }
80 else
81 {
82 if( trace.lgTrace )
83 fprintf( ioQQQ, " ContCreatePointers called first time.\n" );
84 ++nCalled;
85 }
86
87 for( long i=0; i < rfield.nupper; i++ )
88 {
89 /* this is array of labels for lines and continua, set to blanks at first */
90 strcpy( rfield.chContLabel[i], " ");
91 strcpy( rfield.chLineLabel[i], " ");
92 }
93
94 /* we will generate a set of array indices to ionization edges for
95 * the first thirty elements. First set all array indices to
96 * totally bogus values so we will crash if misused */
97 for( long nelem=0; nelem<LIMELM; ++nelem )
98 {
99 if( dense.lgElmtOn[nelem] )
100 {
101 for( long ion=0; ion<LIMELM; ++ion )
102 {
103 for( long nshells=0; nshells<7; ++nshells )
104 {
105 for( long j=0; j<3; ++j )
106 {
107 opac.ipElement[nelem][ion][nshells][j] = INT_MIN;
108 }
109 }
110 }
111 }
112 }
113
114 /* pointer to excited state of O+2 */
115 opac.ipo3exc[0] = ipContEnergy(3.855,"O3ex");
116
117 /* main hydrogenic arrays - THIS OCCURS TWICE!! H and He here, then the
118 * remaining hydrogenic species near the bottom. This is so that H and HeII get
119 * their labels stuffed into the arrays, and the rest of the hydrogenic series
120 * get whatever is left over after the level 1 lines.
121 * to find second block, search on "ipZ=2" */
122 /* NB note that no test for H or He being on exists here - we will always
123 * define the H and He arrays even when He is off, since we need to
124 * know where the he edges are for the bookkeeping that occurs in continuum
125 * binning routines */
126 /* this loop is over H, He-like only - DO NOT CHANGE TO NISO */
127 for( long ipISO=ipH_LIKE; ipISO<=ipHE_LIKE; ++ipISO )
128 {
129 /* this will be over HI, HeII, then HeI only */
130 for( long nelem=ipISO; nelem < 2; nelem++ )
131 {
132 /* generate label for this ion */
133 sprintf( chLab, "%2s%2ld",elementnames.chElementSym[nelem], nelem+1-ipISO );
134
135 /* array index for continuum edges for ground */
136 iso_sp[ipISO][nelem].fb[0].ipIsoLevNIonCon = ipContEnergy(iso_sp[ipISO][nelem].fb[0].xIsoLevNIonRyd,chLab);
137 for( long ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
138 {
139 /* array index for continuum edges for excited levels */
140 iso_sp[ipISO][nelem].fb[ipHi].ipIsoLevNIonCon = ipContEnergy(iso_sp[ipISO][nelem].fb[ipHi].xIsoLevNIonRyd,chLab);
141
142 /* define all line array indices */
143 for( long ipLo=0; ipLo < ipHi; ipLo++ )
144 {
145 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
146 continue;
147
148 /* some lines have negative or zero energy */
149 /* >>chng 03 apr 22, this check was if less than or equal to zero,
150 * changed to lowest energy point so that ultra low energy transitions are
151 * not considered. */
152 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyRyd() < continuum.filbnd[0] )
153 continue;
154
155 /* some energies are negative for inverted levels */
156 iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont() =
157 ipLineEnergy(iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyRyd() , chLab ,
158 iso_sp[ipISO][nelem].fb[ipLo].ipIsoLevNIonCon);
159 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().ipFine() =
160 ipFineCont(iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyRyd() );
161 /* check that energy scales are the same, to within energy resolution of arrays */
162# ifndef NDEBUG
163 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().ipFine() > 0 )
164 {
165 realnum anuCoarse = rfield.anu[iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont()-1];
166 realnum anuFine = rfield.fine_anu[iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().ipFine()];
167 realnum widCoarse = rfield.widflx[iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont()-1];
168 realnum widFine = anuFine - rfield.fine_anu[iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().ipFine()-1];
169 realnum width = MAX2( widFine , widCoarse );
170 /* NB - can only assert more than width of coarse cell
171 * since close to ionization limit, coarse index is
172 * kept below next ionization edge
173 * >>chng 05 mar 02, pre factor below had been 1.5, chng to 2
174 * tripped near H grnd photo threshold */
175 ASSERT( fabs(anuCoarse - anuFine) / anuCoarse <
176 2.*width/anuCoarse);
177 }
178# endif
179 }
180 }/* ipHi loop */
181 }/* nelem loop */
182 }/* ipISO */
183
184 {
185 /* need to increase the cell for the HeII balmer continuum by one, so that
186 * hydrogen Lyman continuum will pick it up. */
187 long nelem = ipHELIUM;
188 /* copy label over to next higher cell */
189 if( strcmp( rfield.chContLabel[iso_sp[ipH_LIKE][nelem].fb[ipH2s].ipIsoLevNIonCon-1] , "He 2" ) == 0)
190 {
191 strcpy( rfield.chContLabel[iso_sp[ipH_LIKE][nelem].fb[ipH2s].ipIsoLevNIonCon],
192 rfield.chContLabel[iso_sp[ipH_LIKE][nelem].fb[ipH2s].ipIsoLevNIonCon-1] );
193 /* set previous spot to blank */
194 strcpy( rfield.chContLabel[iso_sp[ipH_LIKE][nelem].fb[ipH2s].ipIsoLevNIonCon-1] , " ");
195 }
196 else if( strcmp( rfield.chContLabel[iso_sp[ipH_LIKE][nelem].fb[ipH2s].ipIsoLevNIonCon-1] , "H 1" ) == 0)
197 {
198 /* set previous spot to blank */
199 strcpy( rfield.chContLabel[iso_sp[ipH_LIKE][nelem].fb[ipH2s].ipIsoLevNIonCon] , "He 2");
200 }
201 else
202 {
203 fprintf(ioQQQ," insanity heii pointer fix contcreatepointers\n");
204 }
205 /* finally increment the two HeII pointers so that they are above the Lyman continuum */
206 ++iso_sp[ipH_LIKE][nelem].fb[ipH2s].ipIsoLevNIonCon;
207 ++iso_sp[ipH_LIKE][nelem].fb[ipH2p].ipIsoLevNIonCon;
208 }
209
210 /* we will define an array of either 1 or 0 to show whether photooelectron
211 * energy is large enough to produce secondary ionizations
212 * below 100eV no secondary ionization - this is the threshold*/
213 secondaries.ipSecIon = ipoint(7.353);
214
215 /* this is highest energy where k-shell opacities are counted
216 * can be adjusted with "set kshell" command */
217 continuum.KshellLimit = ipoint(continuum.EnergyKshell);
218
219 /* pointers for molecules
220 * H2+ dissociation energy 2.647 eV but cs small below 0.638 Ryd */
221 opac.ih2pnt[0] = ipContEnergy(0.502,"H2+ ");
222 opac.ih2pnt[1] = ipoint(1.03);
223
224 //pointers for most prominent PAH features
225 {
226 /* energies given to ipContEnergy are only to put lave in the right place
227 * wavelengths are rough observed values of blends
228 * 7.6 microns */
229 (void) ipContEnergy(0.0117, "PAH " );
230
231 /* feature near 6.2 microns */
232 (void) ipContEnergy(0.0147, "PAH " );
233
234 /* 3.3 microns */
235 (void) ipContEnergy(0.028, "PAH " );
236
237 /* 11.2 microns */
238 (void) ipContEnergy(0.0080, "PAH " );
239
240 /* 12.3 microns */
241 (void) ipContEnergy(0.0077, "PAH " );
242
243 /* 13.5 microns */
244 (void) ipContEnergy(0.0069, "PAH " );
245 }
246
247 /* fix pointers for hydrogen and helium */
248
249 /* pointer to Bowen 374A resonance line */
250 he.ip374 = ipLineEnergy(1.92,"He 2",0);
251 he.ip660 = ipLineEnergy(1.38,"He 2",0);
252
253 /* pointer to energy defining effect x-ray column */
254 rt.ipxry = ipoint(73.5);
255
256 /* pointer to Hminus edge at 0.754eV */
257 hmi.iphmin = ipContEnergy(0.05544,"H- ");
258
259 /* pointer to threshold for H2 photoionization at 15.4 eV */
260 fixit(); // need to generalize ionization energy and label!
261 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
262 (*diatom)->ip_photo_opac_thresh = ipContEnergy( 15.4/EVRYD , "H2 ");
263
264 hmi.iheh1 = ipoint(1.6);
265 hmi.iheh2 = ipoint(2.3);
266
267 /* pointer to carbon k-shell ionization */
268 opac.ipCKshell = ipoint(20.6);
269
270 /* pointer to threshold for pair production */
271 opac.ippr = ipoint(7.51155e4) + 1;
272
273 /* pointer to x-ray - gamma ray bound; 100 kev */
274 rfield.ipEnerGammaRay = ipoint(rfield.EnerGammaRay);
275
277
278 /* make low energy edges of some cells exact,
279 * index is on c scale
280 * 0.99946 is correct energy of hydrogen in inf mass ryd */
281 fiddle(iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1,0.99946);
282 fiddle(iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH2p].ipIsoLevNIonCon-1,0.24986);
283 /* confirm that labels are in correct location */
284 ASSERT( strcmp( rfield.chContLabel[iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1], "H 1" ) ==0 );
285 ASSERT( strcmp( rfield.chContLabel[iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH2p].ipIsoLevNIonCon-1], "H 1" ) ==0 );
286
287 fiddle(iso_sp[ipH_LIKE][ipHELIUM].fb[ipH1s].ipIsoLevNIonCon-1,4.00000);
288 ASSERT( strcmp( rfield.chContLabel[iso_sp[ipH_LIKE][ipHELIUM].fb[ipH1s].ipIsoLevNIonCon-1], "He 2" ) ==0 );
289
290 /* pointer to excited state of O+2 */
291 fiddle(opac.ipo3exc[0]-1,3.855);
292
293 /* now fix widflx array so that it is correct */
294 for( long i=1; i<rfield.nupper-1; ++i )
295 {
296 rfield.widflx[i] = ((rfield.anu[i+1] - rfield.anu[i]) + (rfield.anu[i] - rfield.anu[i-1]))/2.f;
297 }
298
299 /* these are indices for centers of B and V filters,
300 * taken from table on page 202 of Allen, AQ, 3rd ed */
301 /* the B filter array offset */
302 rfield.ipB_filter = ipoint( RYDLAM / WL_B_FILT );
303 /* the V filter array offset */
304 rfield.ipV_filter = ipoint( RYDLAM / WL_V_FILT );
305
306 /* these are the lower and upper bounds for the G0 radiation field
307 * used by Tielens & Hollenbach in their PDR work */
308 rfield.ipG0_TH85_lo = ipoint( 6.0 / EVRYD );
309 rfield.ipG0_TH85_hi = ipoint( 13.6 / EVRYD );
310
311 /* this is the limits for Draine & Bertoldi Habing field */
312 rfield.ipG0_DB96_lo = ipoint( RYDLAM / 1110. );
313 rfield.ipG0_DB96_hi = ipoint( RYDLAM / 912. );
314
315 /* this is special form of G0 that could be used in some future version, for now,
316 * use default, TH85 */
317 rfield.ipG0_spec_lo = ipoint( 6.0 / EVRYD );
318 rfield.ipG0_spec_hi = ipoint( RYDLAM / 912. );
319
320 /* this index is to 1000A to obtain the extinction at 1000A */
321 rfield.ip1000A = ipoint( RYDLAM / 1000. );
322
323 /* now save current form of array */
324 for( long i=0; i < rfield.nupper; i++ )
325 {
326 rfield.AnuOrg[i] = rfield.anu[i];
327 rfield.anusqr[i] = (realnum)sqrt(rfield.AnuOrg[i]);
328 }
329
330 /* following order of elements is in roughly decreasing abundance
331 * when ipShells gets a cell for the valence IP it does so through
332 * ipContEnergy, which makes sure that no two ions get the same cell
333 * earliest elements have most precise ip mapping */
334
335 /* set up shell pointers for hydrogen */
336 {
337 long nelem = ipHYDROGEN;
338 long ion = 0;
339
340 /* the number of shells */
341 Heavy.nsShells[nelem][0] = 1;
342
343 /*pointer to ionization threshold in energy array*/
344 Heavy.ipHeavy[nelem][ion] = iso_sp[ipH_LIKE][nelem].fb[ipH1s].ipIsoLevNIonCon;
345 opac.ipElement[nelem][ion][0][0] = iso_sp[ipH_LIKE][nelem].fb[ipH1s].ipIsoLevNIonCon;
346
347 /* upper limit to energy integration */
348 opac.ipElement[nelem][ion][0][1] = rfield.nupper;
349 }
350
351 if( dense.lgElmtOn[ipHELIUM] )
352 {
353 /* set up shell pointers for helium */
354 long nelem = ipHELIUM;
355 long ion = 0;
356
357 /* the number of shells */
358 Heavy.nsShells[nelem][0] = 1;
359
360 /*pointer to ionization threshold in energy array*/
361 Heavy.ipHeavy[nelem][ion] = iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon;
362 opac.ipElement[nelem][ion][0][0] = iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon;
363
364 /* upper limit to energy integration */
365 opac.ipElement[nelem][ion][0][1] = rfield.nupper;
366
367 /* (hydrogenic) helium ion */
368 ion = 1;
369 /* the number of shells */
370 Heavy.nsShells[nelem][1] = 1;
371
372 /*pointer to ionization threshold in energy array*/
373 Heavy.ipHeavy[nelem][ion] = iso_sp[ipH_LIKE][nelem].fb[ipH1s].ipIsoLevNIonCon;
374 opac.ipElement[nelem][ion][0][0] = iso_sp[ipH_LIKE][nelem].fb[ipH1s].ipIsoLevNIonCon;
375
376 /* upper limit to energy integration */
377 opac.ipElement[nelem][ion][0][1] = rfield.nupper;
378 }
379
380 /* check that ionization potential of neutral carbon valence shell is
381 * positive */
382 ASSERT( t_ADfA::Inst().ph1(2,5,5,0) > 0. );
383
384 /* now fill in all sub-shell ionization array indices for elements heavier than He,
385 * this must be done after previous loop on iso.ipIsoLevNIonCon[ipH_LIKE] since hydrogenic species use
386 * iso.ipIsoLevNIonCon[ipH_LIKE] rather than ipoint in getting array index within continuum array */
387 for( long i=NISO; i<LIMELM; ++i )
388 {
389 /* i is the atomic number on the c scale, 2 for Li */
390 ipShells(i);
391 }
392
393 /* most of these are set in ipShells, but not h-like or he-like, so do these here */
394 Heavy.Valence_IP_Ryd[ipHYDROGEN][0] = t_ADfA::Inst().ph1(0,0,ipHYDROGEN,0)/EVRYD* 0.9998787;
395 Heavy.Valence_IP_Ryd[ipHELIUM][0] = t_ADfA::Inst().ph1(0,1,ipHELIUM,0)/EVRYD* 0.9998787;
396 Heavy.Valence_IP_Ryd[ipHELIUM][1] = t_ADfA::Inst().ph1(0,0,ipHELIUM,0)/EVRYD* 0.9998787;
397 for( long nelem=2; nelem<LIMELM; ++nelem )
398 {
399 Heavy.Valence_IP_Ryd[nelem][nelem-1] = t_ADfA::Inst().ph1(0,1,nelem,0)/EVRYD* 0.9998787;
400 Heavy.Valence_IP_Ryd[nelem][nelem] = t_ADfA::Inst().ph1(0,0,nelem,0)/EVRYD* 0.9998787;
401 if( dense.lgElmtOn[nelem])
402 {
403 /* now confirm that all are properly set */
404 for( long j=0; j<=nelem; ++j )
405 {
406 ASSERT( Heavy.Valence_IP_Ryd[nelem][j] > 0.05 );
407 }
408 for( long j=0; j<nelem; ++j )
409 {
410 ASSERT( Heavy.Valence_IP_Ryd[nelem][j] < Heavy.Valence_IP_Ryd[nelem][j+1]);
411 }
412 }
413 }
414
415 /* array indices to bound Compton electron recoil ionization of all elements */
416 for( long nelem=0; nelem<LIMELM; ++nelem )
417 {
418 if( dense.lgElmtOn[nelem])
419 {
420 for( long ion=0; ion<nelem+1; ++ion )
421 {
422 /* this is the threshold energy to Compton ionize valence shell electrons */
423 double energy = sqrt( Heavy.Valence_IP_Ryd[nelem][ion] * EN1RYD * ELECTRON_MASS * SPEEDLIGHT * SPEEDLIGHT ) / EN1RYD;
424 /* the array index for this energy */
425 ionbal.ipCompRecoil[nelem][ion] = ipoint( energy );
426 }
427 }
428 }
429
430 /* oxygen pointers for excited states
431 * IO3 is pointer to O++ exc state, is set above */
432 oxy.i2d = ipoint(1.242);
433 oxy.i2p = ipoint(1.367);
434 opac.ipo1exc[0] = ipContEnergy(0.856,"O1ex");//2s^2 2p^4, ^1D level, J=2, energy relative to ground level is ~0.1446 Ry
435 opac.ipo1exc[1] = ipoint(2.0);// upper limit to range of energies where opacity for 1D absorption is defined
436
437 /* upper limit for excited state photoionization
438 * do not use ipContEnergy since only upper limit */
439 opac.ipo3exc[1] = ipoint(5.0);
440
441 /* upper level of 4363 */
442 opac.ipo3exc3[0] = ipContEnergy(3.646,"O3ex");
443 opac.ipo3exc3[1] = ipoint(5.0);
444
445 /* following are various pointers for OI - Ly-beta pump problem
446 * these are delta energies for Boltzmann factors */
451 atoms.ipoiex[0] = ipoint(0.7005);
452 atoms.ipoiex[1] = ipoint(0.10791);
453 atoms.ipoiex[2] = ipoint(0.06925);
454 atoms.ipoiex[3] = ipoint(0.01151);
455 atoms.ipoiex[4] = ipoint(0.01999);
456
457 /* >>chng 97 jan 27, move nitrogen after oxygen so that O gets the
458 * most accurate pointers
459 * Nitrogen
460 * in1(1) is thresh for photo from excited state */
461 opac.in1[0] = ipContEnergy(0.893,"N1ex");
462
463 /* upper limit */
464 opac.in1[1] = ipoint(2.);
465
466 if( (trace.lgTrace && trace.lgConBug) || (trace.lgTrace && trace.lgPointBug) )
467 {
468 fprintf( ioQQQ, " ContCreatePointers:%ld energy cells used. N(1R):%4ld N(1.8):%4ld N(4Ryd):%4ld N(O3)%4ld N(x-ray):%5ld N(rcoil)%5ld\n",
469 rfield.nupper,
470 iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon,
471 iso_sp[ipHE_LIKE][ipHELIUM].fb[ipH1s].ipIsoLevNIonCon,
472 iso_sp[ipH_LIKE][ipHELIUM].fb[ipH1s].ipIsoLevNIonCon,
473 opac.ipo3exc[0],
474 opac.ipCKshell,
475 ionbal.ipCompRecoil[ipHYDROGEN][0] );
476
477
478 fprintf( ioQQQ, " ContCreatePointers: ipEnerGammaRay: %5ld IPPRpari produc%5ld\n",
479 rfield.ipEnerGammaRay, opac.ippr );
480
481 fprintf( ioQQQ, " ContCreatePointers: H pointers;" );
482 for( long i=0; i <= 6; i++ )
483 {
484 fprintf( ioQQQ, "%4ld%4ld", i, iso_sp[ipH_LIKE][ipHYDROGEN].fb[i].ipIsoLevNIonCon );
485 }
486 fprintf( ioQQQ, "\n" );
487
488 fprintf( ioQQQ, " ContCreatePointers: Oxy pnters;" );
489
490 for( long i=1; i <= 8; i++ )
491 {
492 fprintf( ioQQQ, "%4ld%4ld", i, Heavy.ipHeavy[ipOXYGEN][i-1] );
493 }
494 fprintf( ioQQQ, "\n" );
495
496 }
497
498 /* Magnesium
499 * following is energy for phot of MG+ from exc state producing 2798 */
500 opac.ipmgex = ipoint(0.779);
501
502 /* lower, upper edges of Ca+ excited term photoionizaiton */
503 opac.ica2ex[0] = ipContEnergy(0.72,"Ca2x");
504 opac.ica2ex[1] = ipoint(1.);
505
506 /* set up factors and pointers for Fe continuum fluorescence */
507 fe.ipfe10 = ipoint(2.605);
508
509 /* following is WL(CM)**2/(8PI) * conv fac for RYD to NU *A21 */
510 fe.pfe10 = (realnum)(2.00e-18/rfield.widflx[fe.ipfe10-1]);
511
512 /* this is 353 pump, f=0.032 */
513 fe.pfe11a = (realnum)(4.54e-19/rfield.widflx[fe.ipfe10-1]);
514
515 /* this is 341.1 f=0.012 */
516 fe.pfe11b = (realnum)(2.75e-19/rfield.widflx[fe.ipfe10-1]);
517 fe.pfe14 = (realnum)(1.15e-18/rfield.widflx[fe.ipfe10-1]);
518
519 /* set up energy pointers for line optical depth arrays
520 * this also increments flux, sets other parameters for lines */
521
522 /* level 1 heavy elements line array */
523 for( long i=1; i <= nLevel1; i++ )
524 {
525 /* put null terminated line label into chLab using line array*/
526 chIonLbl(chLab, TauLines[i]);
527
528 TauLines[i].ipCont() =
529 ipLineEnergy(TauLines[i].EnergyRyd() , chLab ,0);
530 TauLines[i].Emis().ipFine() =
531 ipFineCont(TauLines[i].EnergyRyd());
532 /* for debugging pointer - index on f not c scale,
533 * this will find all lines that entered a given cell
534 if( TauLines[i].ipCont()==603 )
535 fprintf(ioQQQ,"( level1 %s\n", chLab);*/
536
537 if( TauLines[i].Emis().gf() > 0. )
538 {
539 /* the gf value was entered
540 * derive the A, call to function is gf, wl (A), g_up */
541 TauLines[i].Emis().Aul() =
542 (realnum)(eina(TauLines[i].Emis().gf(),
543 TauLines[i].EnergyWN(),(*TauLines[i].Hi()).g()));
544 }
545 else if( TauLines[i].Emis().Aul() > 0. )
546 {
547 /* the Einstein A value was entered
548 * derive the gf, call to function is A, wl (A), g_up */
549 TauLines[i].Emis().gf() =
550 (realnum)(GetGF(TauLines[i].Emis().Aul(),
551 TauLines[i].EnergyWN(),(*TauLines[i].Hi()).g()));
552 }
553 else
554 {
555 fprintf( ioQQQ, " level 1 line does not have valid gf or A\n" );
556 fprintf( ioQQQ, " This is ContCreatePointers\n" );
558 }
559
560 /* used to get damping constant */
561 TauLines[i].Emis().dampXvel() =
562 (realnum)(TauLines[i].Emis().Aul() / TauLines[i].EnergyWN()/PI4);
563
564 /* derive the abs coefficient, call to function is gf, wl (A), g_low */
565 TauLines[i].Emis().opacity() =
566 (realnum)(abscf(TauLines[i].Emis().gf(),
567 TauLines[i].EnergyWN(),(*TauLines[i].Lo()).g()));
568 /*fprintf(ioQQQ,"TauLinesss\t%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
569 i,TauLines[i].Emis().opacity() , TauLines[i].Emis().gf() , TauLines[i].Emis().Aul() ,TauLines[i].EnergyWN(),(*TauLines[i].Lo()).g());*/
570
571 /*line wavelength must be gt 0 */
572 ASSERT( TauLines[i].WLAng() > 0 );
573 }
574
575 /*Beginning of the dBaseLines*/
576 for (int ipSpecies=0; ipSpecies < nSpecies; ++ipSpecies)
577 {
578 for( EmissionList::iterator em=dBaseTrans[ipSpecies].Emis().begin();
579 em != dBaseTrans[ipSpecies].Emis().end(); ++em)
580 {
581 /* upper level lifetime to calculate the damping parameter.*/
582 (*em).dampXvel() = (realnum)(1./
583 dBaseStates[ipSpecies][em->Tran().ipHi()].lifetime()/em->Tran().EnergyWN()/PI4);
584 (*em).damp() = -1000.0;
585
586 /*Put null terminated line label into chLab*/
587 strncpy(chLab,(*(*em).Tran().Hi()).chLabel(),4);
588 chLab[4]='\0';
589
590 static const double minAul = 1e-29;
591 if( (*em).Aul() > minAul )
592 {
593 (*em).Tran().ipCont() = ipLineEnergy((*em).Tran().EnergyRyd(), chLab ,0);
594 (*em).ipFine() = ipFineCont((*em).Tran().EnergyRyd() );
595 }
596 else
597 {
598 (*em).Tran().ipCont() = -1;
599 (*em).ipFine() = -1;
600 }
601 /* derive the abs coefficient, call to function is gf, wl (A), g_low */
602 (*em).opacity() =
603 (realnum)(abscf((*em).gf(),(*em).Tran().EnergyWN(),
604 (*(*em).Tran().Lo()).g()));
605 }
606 }
607 /*end of the dBaseLines*/
608
609 /* set the ipCont struc element for the H2 molecule */
610 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
611 (*diatom)->H2_ContPoint();
612
613 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
614 {
615 /* do remaining part of the iso sequences */
616 for( long nelem=2; nelem < LIMELM; nelem++ )
617 {
618 if( dense.lgElmtOn[nelem])
619 {
620 /* generate label for this ion */
621 sprintf( chLab, "%2s%2ld",elementnames.chElementSym[nelem], nelem+1-ipISO );
622 /* array index for continuum edges */
623 iso_sp[ipISO][nelem].fb[0].ipIsoLevNIonCon =
624 ipContEnergy(iso_sp[ipISO][nelem].fb[0].xIsoLevNIonRyd,chLab);
625
626 for( long ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
627 {
628 /* array index for continuum edges */
629 iso_sp[ipISO][nelem].fb[ipHi].ipIsoLevNIonCon = ipContEnergy(iso_sp[ipISO][nelem].fb[ipHi].xIsoLevNIonRyd,chLab);
630
631 /* define all line pointers */
632 for( long ipLo=0; ipLo < ipHi; ipLo++ )
633 {
634
635 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
636 continue;
637
638 /* some lines have negative or zero energy */
639 /* >>chng 03 apr 22, this check was if less than or equal to zero,
640 * changed to lowest energy point so that ultra low energy transitions are
641 * not considered. */
642 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyRyd() < continuum.filbnd[0] )
643 continue;
644
645 iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont() =
646 ipLineEnergy(iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyRyd() , chLab ,
647 iso_sp[ipISO][nelem].fb[ipLo].ipIsoLevNIonCon);
648 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().ipFine() =
649 ipFineCont(iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyRyd() );
650 }
651 }
652 iso_sp[ipISO][nelem].fb[0].ipIsoLevNIonCon = ipContEnergy(iso_sp[ipISO][nelem].fb[0].xIsoLevNIonRyd,chLab);
653 }
654 }
655 }
656 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
657 {
658 /* this will be over HI, HeII, then HeI only */
659 for( long nelem=ipISO; nelem < LIMELM; nelem++ )
660 {
661 if( dense.lgElmtOn[nelem])
662 {
663 /* these are the extra Lyman lines */
664 for( long ipHi=2; ipHi < iso_ctrl.nLyman_malloc[ipISO]; ipHi++ )
665 {
666 long ipLo = 0;
667 /* some energies are negative for inverted levels */
668 /*lint -e772 chLab not initialized */
669 TransitionList::iterator tr = ExtraLymanLines[ipISO][nelem].begin()+ipExtraLymanLines[ipISO][nelem][ipHi];
670 (*tr).ipCont() =
671 ipLineEnergy((*tr).EnergyRyd() , chLab ,
672 iso_sp[ipISO][nelem].fb[ipLo].ipIsoLevNIonCon);
673
674 (*tr).Emis().ipFine() =
675 ipFineCont((*tr).EnergyRyd() );
676 /*lint +e772 chLab not initialized */
677 }
678
679 if( iso_ctrl.lgDielRecom[ipISO] )
680 {
681 ASSERT( ipISO>ipH_LIKE );
682 for( long ipLo=0; ipLo<iso_sp[ipISO][nelem].numLevels_max; ipLo++ )
683 {
684
685 SatelliteLines[ipISO][nelem][ipSatelliteLines[ipISO][nelem][ipLo]].ipCont() = ipLineEnergy(
686 SatelliteLines[ipISO][nelem][ipSatelliteLines[ipISO][nelem][ipLo]].EnergyRyd() , chLab ,
687 0);
688
689 SatelliteLines[ipISO][nelem][ipSatelliteLines[ipISO][nelem][ipLo]].Emis().ipFine() =
690 ipFineCont(SatelliteLines[ipISO][nelem][ipSatelliteLines[ipISO][nelem][ipLo]].EnergyRyd() );
691 }
692 }
693 }
694 }
695 }
696
697 fixit(); /* is this redundant? */
698 /* for He-like sequence the majority of the transitions are bogus - A set to special value in this case */
699 {
700 long ipISO = ipHE_LIKE;
701 /* do remaining part of the he iso sequence */
702 for( long nelem=ipISO; nelem < LIMELM; nelem++ )
703 {
704 if( dense.lgElmtOn[nelem])
705 {
706 for( long ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
707 {
708 for( long ipLo=0; ipLo < ipHi; ipLo++ )
709 {
710 TransitionProxy tr = iso_sp[ipISO][nelem].trans(ipHi,ipLo);
711 if( tr.ipCont() <= 0 )
712 continue;
713
714 if( fabs(tr.Emis().Aul() - iso_ctrl.SmallA) < SMALLFLOAT )
715 {
716 /* iso_ctrl.SmallA is value assigned to bogus transitions */
717 tr.ipCont() = -1;
718 tr.Emis().ipFine() = -1;
719 }
720 }
721 }
722 }
723 }
724 }
725
726 /* inner shell transitions */
727 for( long i=0; i<nUTA; ++i )
728 {
729 if( UTALines[i].Emis().Aul() > 0. )
730 {
731
732 // dampXvel is derived in atmdat_readin because autoionization rates
733 // must be included in the total rate for the damping constant
734 ASSERT( UTALines[i].Emis().dampXvel() >0. );
735
736 /* derive the abs coefficient, call to function is gf, wl (A), g_low */
737 UTALines[i].Emis().opacity() =
738 (realnum)(abscf( UTALines[i].Emis().gf(), UTALines[i].EnergyWN(), (*UTALines[i].Lo()).g()));
739
740 /* put null terminated line label into chLab using line array*/
741 chIonLbl(chLab, UTALines[i]);
742
743 /* get pointer to energy in continuum mesh */
744 UTALines[i].ipCont() = ipLineEnergy(UTALines[i].EnergyRyd(), chLab,0 );
745 UTALines[i].Emis().ipFine() = ipFineCont(UTALines[i].EnergyRyd() );
746
747 /* find heating per absorption,
748 * first find threshold for this shell in ergs */
749 /* ionization threshold in erg */
750 double thresh = Heavy.Valence_IP_Ryd[(*UTALines[i].Hi()).nelem()-1][(*UTALines[i].Hi()).IonStg()-1] *EN1RYD;
751 UTALines[i].Coll().heat() = (UTALines[i].EnergyErg()-thresh);
752 ASSERT( UTALines[i].Coll().heat()> 0. );
753 }
754 }
755
756 /* level 2 heavy element lines */
757 for( long i=0; i < nWindLine; i++ )
758 {
759 /* derive the A, call to function is gf, wl (A), g_up */
760 TauLine2[i].Emis().Aul() =
761 (realnum)(eina(TauLine2[i].Emis().gf(),
762 TauLine2[i].EnergyWN(),(*TauLine2[i].Hi()).g()));
763
764 /* coefficient needed for damping constant - units cm s-1 */
765 TauLine2[i].Emis().dampXvel() =
766 (realnum)(TauLine2[i].Emis().Aul()/
767 TauLine2[i].EnergyWN()/PI4);
768
769 /* derive the abs coefficient, call to function is gf, wl (A), g_low */
770 TauLine2[i].Emis().opacity() =
771 (realnum)(abscf(TauLine2[i].Emis().gf(),
772 TauLine2[i].EnergyWN(),(*TauLine2[i].Lo()).g()));
773
774 /* put null terminated line label into chLab using line array*/
775 chIonLbl(chLab, TauLine2[i]);
776
777 /* get pointer to energy in continuum mesh */
778 TauLine2[i].ipCont() = ipLineEnergy(TauLine2[i].EnergyRyd(), chLab,0 );
779 TauLine2[i].Emis().ipFine() = ipFineCont(TauLine2[i].EnergyRyd() );
780 /*if( TauLine2[i].ipCont()==751 )
781 fprintf(ioQQQ,"( atom_level2 %s\n", chLab);*/
782 }
783
784 /* hyperfine structure lines */
785 for( long i=0; i < nHFLines; i++ )
786 {
787 ASSERT( HFLines[i].Emis().Aul() > 0. );
788 /* coefficient needed for damping constant */
789 HFLines[i].Emis().dampXvel() =
790 (realnum)(HFLines[i].Emis().Aul()/
791 HFLines[i].EnergyWN()/PI4);
792 HFLines[i].Emis().damp() = 1e-20f;
793
794 /* derive the abs coefficient, call to function is gf, wl (A), g_low */
795 HFLines[i].Emis().opacity() =
796 (realnum)(abscf(HFLines[i].Emis().gf(),
797 HFLines[i].EnergyWN(),
798 (*HFLines[i].Lo()).g()));
799 /* gf from this and 21 cm do not agree, A for HFS is 10x larger than level1 dat */
800 /*fprintf(ioQQQ,"HFLinesss\t%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
801 i,HFLines[i].Emis().opacity() , HFLines[i].Emis().gf() , HFLines[i].Emis().Aul() , HFLines[i].EnergyWN(),(*HFLines[i].Lo()).g());*/
802
803 /* put null terminated line label into chLab using line array*/
804 chIonLbl(chLab, HFLines[i]);
805
806 /* get pointer to energy in continuum mesh */
807 HFLines[i].ipCont() = ipLineEnergy(HFLines[i].EnergyRyd() , chLab,0 );
808 HFLines[i].Emis().ipFine() = ipFineCont(HFLines[i].EnergyRyd() );
809 }
810
811 /* Verner's FeII lines - do first so following labels will over write this
812 * only call if large FeII atom is turned on */
813 FeIIPoint();
814
815 /* the group of inner shell fluorescent lines */
816 for( long i=0; i < t_yield::Inst().nlines(); ++i )
817 {
818 strcpy( chLab , elementnames.chElementSym[t_yield::Inst().nelem(i)] );
819 strcat( chLab , elementnames.chIonStage[t_yield::Inst().ion_emit(i)] );
820
821 t_yield::Inst().set_ipoint( i, ipLineEnergy( t_yield::Inst().energy(i) , chLab , 0 ) );
822 }
823
824 /* ================================================================================== */
825 /* two photon two-photon 2-nu 2 nu 2 photon 2-photon */
826
827 /* now loop over the two iso-sequences with two photon continua */
828 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
829 {
830 /* set up two photon emission */
831 for( long nelem=ipISO; nelem<LIMELM; ++nelem )
832 {
833 if( dense.lgElmtOn[nelem] )
834 {
835 // upper level for two-photon emission in H and He iso sequences
836 // The 1+ipISO is not rigorous, it just works for H- and He-like
837 const int TwoS = (1+ipISO);
838 double Aul;
839 /* 2s two photon */
840 if( ipISO==ipH_LIKE )
841 Aul = 8.226*pow((double)(nelem+1.),6.);
842 else
843 {
844 ASSERT( ipISO==ipHE_LIKE );
845 /* >>refer Helike 2pho Derevianko, A., & Johnson, W. R. 1997, Phys. Rev. A, 56, 1288
846 * numbers are not explicitly given in this paper for Z=21-24,26,27,and 29.
847 * So numbers given here are interpolated. */
848 fixit(); // where is 51.02 from? Value is 51.3 from the Derevianko & Johnson paper cited above.
849 const double As2nuFrom1S[29] = {51.02,1940.,1.82E+04,9.21E+04,3.30E+05,9.44E+05,2.31E+06,5.03E+06,1.00E+07,
850 1.86E+07,3.25E+07,5.42E+07,8.69E+07,1.34E+08,2.02E+08,2.96E+08,4.23E+08,5.93E+08,8.16E+08,
851 1.08E+09,1.43E+09,1.88E+09,2.43E+09,3.25E+09,3.95E+09,4.96E+09,6.52E+09,7.62E+09,9.94E+09};
852 Aul = As2nuFrom1S[nelem-1];
853 }
854
855 TwoPhotonSetup( iso_sp[ipISO][nelem].TwoNu, TwoS, 0,
856 Aul,
857 iso_sp[ipISO][nelem].trans(TwoS,0),
858 ipISO, nelem );
859 }
860 }
861 }
862
863 // add He-like 2nu 2^3S - 1^1S
864 {
865 const long ipISO = ipHE_LIKE;
866 for( long nelem=ipISO; nelem<LIMELM; ++nelem )
867 {
868 if( dense.lgElmtOn[nelem] )
869 {
870 /* Important clarification, according to Derevianko & Johnson (see ref above), 2^3S can decay
871 * to ground in one of two ways: through a two-photon process, or through a single-photon M1 decay,
872 * but the M1 rates are about 10^4 greater that the two-photon decays throughout the entire
873 * sequence. Thus these numbers, are much weaker than the effective decay rate, but should probably
874 * be treated in as a two-photon decay at some point */
875 // >> refer He As Drake, G. W. F., Victor, G. A., & Dalgarno, A. 1969, Physical Review, 180, 25
876 const double As2nuFrom3S[29] = {4.09e-9,1.25E-06,5.53E-05,8.93E-04,8.05E-03,4.95E-02,2.33E-01,8.94E-01,2.95E+00,
877 8.59E+00,2.26E+01,5.49E+01,1.24E+02,2.64E+02,5.33E+02,1.03E+03,1.91E+03,3.41E+03,5.91E+03,
878 9.20E+03,1.50E+04,2.39E+04,3.72E+04,6.27E+04,8.57E+04,1.27E+05,2.04E+05,2.66E+05,4.17E+05};
879
880 TwoPhotonSetup( iso_sp[ipISO][nelem].TwoNu, ipHe2s3S, ipHe1s1S,
881 As2nuFrom3S[nelem-1],
882 iso_sp[ipISO][nelem].trans(ipHe2s3S,ipHe1s1S),
883 ipISO, nelem );
884 }
885 }
886 }
887
888 {
889 /* this is an option to print out one of the two photon continua */
890 enum {DEBUG_LOC=false};
891 if( DEBUG_LOC )
892 {
893 const int nCRS = 21;
894 double ener[nCRS]={
895 0., 0.03738, 0.07506, 0.1124, 0.1498, 0.1875,
896 0.225, 0.263, 0.300, 0.3373, 0.375, 0.4127,
897 0.4500, 0.487, 0.525, 0.5625, 0.6002, 0.6376,
898 0.6749, 0.7126, 0.75};
899
900 long nelem = ipHYDROGEN;
901 long ipISO = ipHYDROGEN;
902 two_photon& tnu = iso_sp[ipISO][nelem].TwoNu[0];
903
904 double limit = tnu.ipTwoPhoE;
905
906 for( long i=0; i < nCRS; i++ )
907 {
908 fprintf(ioQQQ,"%.3e\t%.3e\n", ener[i] ,
909 atmdat_2phot_shapefunction( ener[i]/0.75, ipISO, nelem ) );
910 }
911
912 double xnew = 0.;
914 for( long i=0; i < limit; i++ )
915 {
916 double fach = tnu.As2nu[i]/2.*rfield.anu2[i]/rfield.widflx[i]*EN1RYD;
917 fprintf(ioQQQ,"%.3e\t%.3e\t%.3e\n",
918 rfield.anu[i] ,
919 tnu.As2nu[i] / rfield.widflx[i] ,
920 fach );
921 xnew += tnu.As2nu[i];
922 }
923 fprintf(ioQQQ," sum is %.3e\n", xnew );
925 }
926 }
927
928 {
929 enum {DEBUG_LOC=false};
930 if( DEBUG_LOC )
931 {
932 for( long i=0; i<11; ++i )
933 {
934 char chLsav[10];
935 (*TauDummy).WLAng() = (realnum)(PI * pow(10.,(double)i));
936 strcpy( chLsav, chLineLbl(*TauDummy) );
937 fprintf(ioQQQ,"%.2f\t%s\n", (*TauDummy).WLAng() , chLsav );
938 }
940 }
941 }
942
943 /* option to print out whole thing with "trace lines" command */
944 if( trace.lgTrLine )
945 {
946 fprintf( ioQQQ, " WL(Ang) E(RYD) IP gl gu gf A damp abs K\n" );
947 for( long i=1; i <= nLevel1; i++ )
948 {
949 strcpy( chLab, chLineLbl(TauLines[i]) );
950 long iWL_Ang = (long)TauLines[i].WLAng();
951 if( iWL_Ang > 1000000 )
952 {
953 iWL_Ang /= 10000;
954 }
955 else if( iWL_Ang > 10000 )
956 {
957 iWL_Ang /= 1000;
958 }
959
960 fprintf( ioQQQ, " %10.10s%5ld%10.3e %4li%4ld%4ld%10.2e%10.2e%10.2e%10.2e\n",
961 chLab, iWL_Ang, RYDLAM/TauLines[i].WLAng(),
962 TauLines[i].ipCont(), (long)((*TauLines[i].Lo()).g()),
963 (long)((*TauLines[i].Hi()).g()), TauLines[i].Emis().gf(),
964 TauLines[i].Emis().Aul(), TauLines[i].Emis().dampXvel(),
965 TauLines[i].Emis().opacity() );
966 }
967
968 /*Atomic Or Molecular lines*/
969 for (int ipSpecies=0; ipSpecies < nSpecies; ++ipSpecies)
970 {
971 for( EmissionList::iterator em=dBaseTrans[ipSpecies].Emis().begin();
972 em != dBaseTrans[ipSpecies].Emis().end(); ++em)
973 {
974 strcpy( chLab, chLineLbl((*em).Tran()));
975
976 long iWL_Ang = (long)(*em).Tran().WLAng();
977
978 if( iWL_Ang > 1000000 )
979 {
980 iWL_Ang /= 10000;
981 }
982 else if( iWL_Ang > 10000 )
983 {
984 iWL_Ang /= 1000;
985 }
986 fprintf( ioQQQ, " %10.10s%5ld%10.3e %4li%4ld%4ld%10.2e%10.2e%10.2e%10.2e\n",
987 chLab, iWL_Ang, RYDLAM/(*em).Tran().WLAng(),
988 (*em).Tran().ipCont(), (long)((*(*em).Tran().Lo()).g()),
989 (long)((*(*em).Tran().Hi()).g()),(*em).gf(),
990 (*em).Aul(),(*em).dampXvel(),
991 (*em).opacity());
992 }
993 }
994
995 for( long i=0; i < nWindLine; i++ )
996 {
997 strcpy( chLab, chLineLbl(TauLine2[i]) );
998
999 long iWL_Ang = (long)TauLine2[i].WLAng();
1000
1001 if( iWL_Ang > 1000000 )
1002 {
1003 iWL_Ang /= 10000;
1004 }
1005 else if( iWL_Ang > 10000 )
1006 {
1007 iWL_Ang /= 1000;
1008 }
1009 fprintf( ioQQQ, " %10.10s%5ld%10.3e %4li%4ld%4ld%10.2e%10.2e%10.2e%10.2e\n",
1010 chLab, iWL_Ang, RYDLAM/TauLine2[i].WLAng(),
1011 TauLine2[i].ipCont(), (long)((*TauLine2[i].Lo()).g()),
1012 (long)((*TauLine2[i].Hi()).g()), TauLine2[i].Emis().gf(),
1013 TauLine2[i].Emis().Aul(), TauLine2[i].Emis().dampXvel(),
1014 TauLine2[i].Emis().opacity() );
1015 }
1016 for( long i=0; i < nHFLines; i++ )
1017 {
1018 strcpy( chLab, chLineLbl(HFLines[i]) );
1019
1020 long iWL_Ang = (long)HFLines[i].WLAng();
1021
1022 if( iWL_Ang > 1000000 )
1023 {
1024 iWL_Ang /= 10000;
1025 }
1026 else if( iWL_Ang > 10000 )
1027 {
1028 iWL_Ang /= 1000;
1029 }
1030 fprintf( ioQQQ, " %10.10s%5ld%10.3e %4li%4ld%4ld%10.2e%10.2e%10.2e%10.2e\n",
1031 chLab, iWL_Ang, RYDLAM/HFLines[i].WLAng(),
1032 HFLines[i].ipCont(), (long)((*HFLines[i].Lo()).g()),
1033 (long)((*HFLines[i].Hi()).g()), HFLines[i].Emis().gf(),
1034 HFLines[i].Emis().Aul(), HFLines[i].Emis().dampXvel(),
1035 HFLines[i].Emis().opacity() );
1036 }
1037 }
1038
1039 /* this is an option to kill fine structure line optical depths */
1040 if( !rt.lgFstOn )
1041 {
1042 for( long i=1; i <= nLevel1; i++ )
1043 {
1044 if( TauLines[i].EnergyWN() < 10000. )
1045 {
1046 TauLines[i].Emis().opacity() = 0.;
1047 }
1048 }
1049
1050 /*Atomic or Molecular Lines-Humeshkar Nemala*/
1051 for (int ipSpecies=0; ipSpecies < nSpecies; ++ipSpecies)
1052 {
1053 for( EmissionList::iterator em=dBaseTrans[ipSpecies].Emis().begin();
1054 em != dBaseTrans[ipSpecies].Emis().end(); ++em)
1055 {
1056 if((*em).Tran().EnergyWN() < 10000. )
1057 {
1058 (*em).opacity() = 0.;
1059 }
1060 }
1061 }
1062 }
1063
1064 /* read in continuum bands data set */
1065 ContBandsCreate( "" );
1066
1067 /* we're done adding lines and states to the stacks.
1068 * This flag is used to make sure we never add them again in this coreload. */
1069 lgLinesAdded = true;
1070 lgStatesAdded = true;
1071
1073
1074 return;
1075}
1076
1077/*fiddle adjust energy bounds of cell with index ipnt so that lower energy
1078 * matches ionization edges exactly, called by createpoint */
1079/* ipnt is index on c scale */
1080STATIC void fiddle(long int ipnt,
1081 double exact)
1082{
1083 realnum Ehi,
1084 Elo,
1085 OldEner;
1086
1087 DEBUG_ENTRY( "fiddle()" );
1088
1089 /* make low edge of cell exact for photo integrals */
1090 ASSERT( ipnt >= 0 );
1091 ASSERT( ipnt < rfield.nupper-1 );
1092
1093 /* upper edge of higher cell*/
1094 Ehi = rfield.anu[ipnt] + rfield.widflx[ipnt]*0.5f;
1095 /* lower edge of lower cell */
1096 Elo = rfield.anu[ipnt-1] - rfield.widflx[ipnt-1]*0.5f;
1097
1098 /* >>chng 02 nov 11, do nothing if already very close,
1099 * because VERY high res models had negative widths for some very close edges */
1100 if( fabs( Elo/exact - 1. ) < 0.001 )
1101 return;
1102
1103 ASSERT( Elo <= exact );
1104
1105 OldEner = rfield.anu[ipnt];
1106
1107 /* centroid of desired lower energy and upper edge */
1108 rfield.anu[ipnt] = (realnum)((Ehi + exact)/2.);
1109 /* same for lower */
1110 rfield.anu[ipnt-1] = (realnum)((exact + Elo)/2.);
1111
1112 rfield.widflx[ipnt] = (realnum)(Ehi - exact);
1113 rfield.widflx[ipnt-1] = (realnum)(exact - Elo);
1114
1115 /* bring upper cell down a bit too, since we dont want large change in widflx */
1116 rfield.anu[ipnt+1] -= (OldEner - rfield.anu[ipnt])/2.f;
1117
1118 /* sanity check */
1119 ASSERT( rfield.widflx[ipnt-1] > 0. );
1120 ASSERT( rfield.widflx[ipnt] > 0. );
1121 return;
1122}
1123
1124/*ipShells assign continuum energy pointers to shells for all atoms,
1125 * called by ContCreatePointers */
1127 /* nelem is the atomic number on the C scale, Li is 2 */
1128 long int nelem)
1129{
1130 char chLab[5];
1131 long int
1132 imax,
1133 ion,
1134 nelec,
1135 ns,
1136 nshell;
1137 /* following value cannot be used - will be set to proper threshold */
1138 double thresh=-DBL_MAX;
1139
1140 DEBUG_ENTRY( "ipShells()" );
1141
1142 ASSERT( nelem >= NISO);
1143 ASSERT( nelem < LIMELM );
1144
1145 /* fills in pointers to valence shell ionization threshold
1146 * PH1(a,b,c,d)
1147 * a=1 => thresh, others fitting parameters
1148 * b atomic number
1149 * c number of electrons
1150 * d shell number 7-1 */
1151
1152 /* threshold in Ryd
1153 * ion=0 for atom, up to nelem-1 for helium like, hydrogenic is elsewhere */
1154 for( ion=0; ion < nelem; ion++ )
1155 {
1156 /* generate label for ionization stage */
1157 /* this is short form of element name */
1158 strcpy( chLab, elementnames.chElementSym[nelem] );
1159
1160 /* this is a number between 1 and 31 */
1161 strcat( chLab, elementnames.chIonStage[ion] );
1162
1163 /* this is the iso sequence - must not redo sequence if done as iso */
1164 long int ipISO = nelem-ion;
1165
1166 /* number of bound electrons */
1167 nelec = ipISO+1;
1168
1169 /* nsShells(nelem,ion) is the number of shells for ion with nelec electrons,
1170 * physical not c scale */
1171 imax = Heavy.nsShells[nelem][ion];
1172
1173 /* loop on all inner shells, valence shell */
1174 for( nshell=0; nshell < imax; nshell++ )
1175 {
1176 /* ionization potential of this shell in rydbergs */
1177 thresh = (double)(t_ADfA::Inst().ph1(nshell,nelec-1,nelem,0)/EVRYD* 0.9998787);
1178 if( thresh <= 0.1 )
1179 {
1180 /* negative ip shell does not exist, set upper limit
1181 * to less than lower limit so this never looped upon
1182 * these are used as flags by LimitSh to check whether
1183 * this is a real shell - if 1 or 2 is changed - change LimitSh!! */
1184 opac.ipElement[nelem][ion][nshell][0] = 2;
1185 opac.ipElement[nelem][ion][nshell][1] = 1;
1186 }
1187 else
1188 {
1189 /* this is lower bound to energy range for this shell */
1190 /* >>chng 02 may 27, change to version of ip with label, so that
1191 * inner shell edges will appear */
1192 /*opac.ipElement[nelem][ion][nshell][0] = ipoint(thresh);*/
1193 opac.ipElement[nelem][ion][nshell][0] =
1194 ipContEnergy( thresh , chLab );
1195
1196 /* this is upper bound to energy range for this shell
1197 * LimitSh is an integer function, returns pointer
1198 * to threshold of next major shell. For k-shell it
1199 * returns the values KshellLimit, default=7.35e4
1200 * >>chng 96 sep 26, had been below, result zero cross sec at
1201 * many energies where opacity project did not produce state specific
1202 * cross section */
1203 opac.ipElement[nelem][ion][nshell][1] =
1204 LimitSh(ion+1, nshell+1,nelem+1);
1205 ASSERT( opac.ipElement[nelem][ion][nshell][1] > 0);
1206 }
1207 }
1208
1209 ASSERT( imax > 0 && imax <= 7 );
1210
1211 /* this will be index pointing to valence edge */
1212 /* [0] is pointer to threshold in energy array */
1213 opac.ipElement[nelem][ion][imax-1][0] =
1214 ipContEnergy(thresh, chLab);
1215
1216 /* pointer to valence electron ionization potential */
1217 Heavy.ipHeavy[nelem][ion] = opac.ipElement[nelem][ion][imax-1][0];
1218 ASSERT( Heavy.ipHeavy[nelem][ion]>0 );
1219
1220 /* ionization potential of valence shell in Ryd
1221 * thresh was evaluated above, now has last value, the valence shell */
1222 Heavy.Valence_IP_Ryd[nelem][ion] = thresh;
1223
1224 Heavy.xLyaHeavy[nelem][ion] = 0.;
1225 if( ipISO >= NISO )
1226 {
1227 /* this is set of 3/4 of valence shell IP, this is important
1228 * source of ots deep in cloud */
1229 Heavy.ipLyHeavy[nelem][ion] =
1230 ipLineEnergy(thresh*0.75,chLab , 0);
1231
1232 Heavy.ipBalHeavy[nelem][ion] =
1233 ipLineEnergy(thresh*0.25,chLab , 0);
1234 }
1235 else
1236 {
1237 /* do not treat this simple way since done exactly with iso
1238 * sequences */
1239 Heavy.ipLyHeavy[nelem][ion] = -1;
1240 Heavy.ipBalHeavy[nelem][ion] = -1;
1241 }
1242 }
1243
1244 /* above loop did up to hydrogenic, now do hydrogenic -
1245 * hydrogenic is special since arrays already set up */
1246 Heavy.nsShells[nelem][nelem] = 1;
1247
1248 /* this is lower limit to range */
1249 /* hydrogenic photoionization set to special hydro array
1250 * this is pointer to threshold energy */
1251 /* this statement is in ContCreatePointers but has not been done when this routine called */
1252 /*iso_sp[ipH_LIKE][ipZ].fb[ipLo].ipIsoLevNIonCon = ipContEnergy(iso_sp[ipH_LIKE][ipZ].fb[ipLo].xIsoLevNIonRyd,chLab);*/
1253 /*opac.ipElement[nelem][nelem][0][0] = iso_sp[ipH_LIKE][nelem].fb[ipH1s].ipIsoLevNIonCon;*/
1254 opac.ipElement[nelem][nelem][0][0] = ipoint( t_ADfA::Inst().ph1(0,0,nelem,0)/EVRYD* 0.9998787 );
1255 ASSERT( opac.ipElement[nelem][nelem][0][0] > 0 );
1256
1257 /* this is the high-energy limit */
1258 opac.ipElement[nelem][nelem][0][1] = continuum.KshellLimit;
1259
1260 Heavy.ipHeavy[nelem][nelem] = opac.ipElement[nelem][nelem][0][0];
1261
1262 /* this is for backwards computability with Cambridge code */
1263 if( trace.lgTrace && trace.lgPointBug )
1264 {
1265 for( ion=0; ion < (nelem+1); ion++ )
1266 {
1267 fprintf( ioQQQ, "Ion:%3ld%3ld %2.2s%2.2s total shells:%3ld\n",
1268 nelem, ion+1, elementnames.chElementSym[nelem], elementnames.chIonStage[ion]
1269 , Heavy.nsShells[nelem][ion] );
1270 for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
1271 {
1272 fprintf( ioQQQ, " shell%3ld %2.2s range eV%10.2e-%8.2e\n",
1273 ns+1, Heavy.chShell[ns], rfield.anu[opac.ipElement[nelem][ion][ns][0]-1]*
1274 EVRYD, rfield.anu[opac.ipElement[nelem][ion][ns][1]-1]*EVRYD );
1275 }
1276 }
1277 }
1278 return;
1279}
1280
1281/*LimitSh sets upper energy limit to subshell integrations */
1282STATIC long LimitSh(long int ion,
1283 long int nshell,
1284 long int nelem)
1285{
1286 long int LimitSh_v;
1287
1288 DEBUG_ENTRY( "LimitSh()" );
1289
1290 /* this routine returns the high-energy limit to the energy range
1291 * for photoionization of a given shell
1292 * */
1293 if( nshell == 1 )
1294 {
1295 /* this limit is high-energy limit to code unless changed with set kshell */
1296 LimitSh_v = continuum.KshellLimit;
1297
1298 }
1299 else if( nshell == 2 )
1300 {
1301 /* this is 2s shell, upper limit is 1s
1302 * >>chng 96 oct 08, up to high-energy limit
1303 * LimitSh = ipElement(nelem,ion , 1,1)-1 */
1304 LimitSh_v = continuum.KshellLimit;
1305
1306 }
1307 else if( nshell == 3 )
1308 {
1309 /* this is 2p shell, upper limit is 1s
1310 * >>chng 96 oct 08, up to high-energy limit
1311 * LimitSh = ipElement(nelem,ion , 1,1)-1 */
1312 LimitSh_v = continuum.KshellLimit;
1313
1314 }
1315 else if( nshell == 4 )
1316 {
1317 /* this is 3s shell, upper limit is 2p
1318 * >>chng 96 oct 08, up to K-shell edge
1319 * LimitSh = ipElement(nelem,ion , 3,1)-1 */
1320 LimitSh_v = opac.ipElement[nelem-1][ion-1][0][0] - 1;
1321
1322 }
1323 else if( nshell == 5 )
1324 {
1325 /* this is 3p shell, upper limit is 2p
1326 * >>chng 96 oct 08, up to K-shell edge
1327 * LimitSh = ipElement(nelem,ion , 3,1)-1 */
1328 LimitSh_v = opac.ipElement[nelem-1][ion-1][0][0] - 1;
1329
1330 }
1331 else if( nshell == 6 )
1332 {
1333 /* this is 3d shell, upper limit is 2p
1334 * >>chng 96 oct 08, up to K-shell edge
1335 * LimitSh = ipElement(nelem,ion , 3,1)-1 */
1336 LimitSh_v = opac.ipElement[nelem-1][ion-1][0][0] - 1;
1337
1338 }
1339 else if( nshell == 7 )
1340 {
1341 /* this is 4s shell, upper limit is 3d */
1342 if( opac.ipElement[nelem-1][ion-1][5][0] < 3 )
1343 {
1344 /* this is check for empty shell 6, 3d
1345 * if so then set to 3p instead */
1346 LimitSh_v = opac.ipElement[nelem-1][ion-1][4][0] -
1347 1;
1348 }
1349 else
1350 {
1351 LimitSh_v = opac.ipElement[nelem-1][ion-1][5][0] -
1352 1;
1353 }
1354 /* >>chng 96 sep 26, set upper limit down to 2s */
1355 LimitSh_v = opac.ipElement[nelem-1][ion-1][2][0] - 1;
1356
1357 }
1358 else
1359 {
1360 fprintf( ioQQQ, " LimitSh cannot handle nshell as large as%4ld\n",
1361 nshell );
1363 }
1364 return( LimitSh_v );
1365}
1366
1367/*ContBandsCreate - read set of continuum bands to enter total emission into line*/
1369 /* chFile is optional filename, if void then use default bands,
1370 * if not void then use file specified,
1371 * return value is 0 for success, 1 for failure */
1372 const char chFile[] )
1373{
1374 char chLine[FILENAME_PATH_LENGTH_2];
1375 const char* chFilename;
1376 FILE *ioDATA;
1377
1378 bool lgEOL;
1379 long int i,k;
1380
1381 /* keep track of whether we have been called - want to be
1382 * called a total of one time */
1383 static bool lgCalled=false;
1384
1385 DEBUG_ENTRY( "ContBandsCreate()" );
1386
1387 /* do nothing if second or later call*/
1388 if( lgCalled )
1389 {
1390 /* success */
1391 return;
1392 }
1393 lgCalled = true;
1394
1395 /* use default filename if void string, else use file specified */
1396 chFilename = ( strlen(chFile) == 0 ) ? "continuum_bands.ini" : chFile;
1397
1398 /* get continuum band data */
1399 if( trace.lgTrace )
1400 {
1401 fprintf( ioQQQ, " ContBandsCreate opening %s:", chFilename );
1402 }
1403
1404 ioDATA = open_data( chFilename, "r" );
1405
1406 /* now count how many bands are in the file */
1407 continuum.nContBand = 0;
1408
1409 /* first line is a magic number and does not count as a band*/
1410 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1411 {
1412 fprintf( ioQQQ, " ContBandsCreate could not read first line of %s.\n", chFilename );
1414 }
1415 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
1416 {
1417 /* we want to count the lines that do not start with #
1418 * since these contain data */
1419 if( chLine[0] != '#')
1420 ++continuum.nContBand;
1421 }
1422
1423 /* now rewind the file so we can read it a second time*/
1424 if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
1425 {
1426 fprintf( ioQQQ, " ContBandsCreate could not rewind %s.\n", chFilename );
1428 }
1429
1430 continuum.ContBandWavelength = (realnum *)MALLOC(sizeof(realnum)*(unsigned)(continuum.nContBand) );
1431 continuum.chContBandLabels = (char **)MALLOC(sizeof(char *)*(unsigned)(continuum.nContBand) );
1432 continuum.ipContBandLow = (long int *)MALLOC(sizeof(long int)*(unsigned)(continuum.nContBand) );
1433 continuum.ipContBandHi = (long int *)MALLOC(sizeof(long int)*(unsigned)(continuum.nContBand) );
1434 continuum.BandEdgeCorrLow = (realnum *)MALLOC(sizeof(realnum)*(unsigned)(continuum.nContBand) );
1435 continuum.BandEdgeCorrHi = (realnum *)MALLOC(sizeof(realnum)*(unsigned)(continuum.nContBand) );
1436
1437 /* now make second dim, id wavelength, and lower and upper bounds */
1438 for( i=0; i<continuum.nContBand; ++i )
1439 {
1440 /* array of labels, each 4 long plus 0 at [4] */
1441 continuum.chContBandLabels[i] = (char *)MALLOC(sizeof(char)*(unsigned)(5) );
1442 }
1443
1444 /* first line is a versioning magic number - now confirm that it is valid */
1445 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1446 {
1447 fprintf( ioQQQ, " ContBandsCreate could not read first line of %s.\n", chFilename );
1449 }
1450 /* bands_continuum magic number here <- this string is in band_continuum.dat
1451 * with comment to search for this to find magic number */
1452 {
1453 long int m1 , m2 , m3,
1454 // the magic number at the start of the data file
1455 myr = 11, mmo = 9, mdy = 10;
1456
1457 i = 1;
1458 m1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
1459 m2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
1460 m3 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
1461 if( ( m1 != myr ) ||
1462 ( m2 != mmo ) ||
1463 ( m3 != mdy ) )
1464 {
1465 fprintf( ioQQQ,
1466 " ContBandsCreate: the version of the data file %s I found (%li %li %li)is not the current version (%li %li %li).\n",
1467 chFilename ,
1468 m1 , m2 , m3 ,
1469 myr , mmo , mdy );
1470 fprintf( ioQQQ,
1471 " ContBandsCreate: you need to update this file.\n");
1473 }
1474 }
1475
1476 /* now read in data again, but save it this time */
1477 k = 0;
1478 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
1479 {
1480 /* we want to count the lines that do not start with #
1481 * since these contain data */
1482 if( chLine[0] != '#')
1483 {
1484 /* copy 4 char label plus termination */
1485 strncpy( continuum.chContBandLabels[k] , chLine , 4 );
1486 continuum.chContBandLabels[k][4] = 0;
1487
1488 /* now get central band wavelength
1489 * >>chng 06 aug 11 from 4 to 6, the first 4 char are labels and
1490 * these can contain numbers, next comes a space, then the number */
1491 i = 6;
1492 continuum.ContBandWavelength[k] = (realnum)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
1493 /* >>chng 06 feb 21, multiply by 1e4 to convert micron wavelength into Angstroms,
1494 * which is assumed by the code. before this correction the band centroid
1495 * wavelength was given in the output incorrectly listed as Angstroms.
1496 * results were correct just label was wrong */
1497 continuum.ContBandWavelength[k] *= 1e4f;
1498
1499 /* these are short and long wave limits, which are high and
1500 * low energy limits - these are now wl in microns but are
1501 * converted to Angstroms */
1502 double xHi = FFmtRead(chLine,&i,sizeof(chLine),&lgEOL)*1e4;
1503 double xLow = FFmtRead(chLine,&i,sizeof(chLine),&lgEOL)*1e4;
1504 if( lgEOL )
1505 {
1506 fprintf( ioQQQ, " There should have been 3 numbers on this band line. Sorry.\n" );
1507 fprintf( ioQQQ, " string==%s==\n" ,chLine );
1509 }
1510
1511 {
1512 enum {DEBUG_LOC=false};
1513 if( DEBUG_LOC )
1514 {
1515 fprintf(ioQQQ, "READ:%s\n", chLine );
1516 fprintf(ioQQQ, "GOT: %s %g %g %g\n",continuum.chContBandLabels[k],
1517 continuum.ContBandWavelength[k] , xHi , xLow );
1518 }
1519 }
1520
1521 /* make sure bands bounds are in correct order, shorter - longer wavelength*/
1522 if( xHi >= xLow )
1523 {
1524 fprintf( ioQQQ, " ContBandWavelength band %li "
1525 "edges are in improper order.\n" ,k);
1526 fprintf(ioQQQ,"band: %s %.3e %.3e %.3e \n",
1527 continuum.chContBandLabels[k],
1528 continuum.ContBandWavelength[k],
1529 xHi ,
1530 xLow);
1532 }
1533
1534 // check that central wavelength is indeed between the limits
1535 // xHi & xLow are hi and low energy limits to band so logic reversed
1536 if( continuum.ContBandWavelength[k] < xHi ||
1537 continuum.ContBandWavelength[k] > xLow )
1538 {
1539 fprintf( ioQQQ, " ContBandWavelength band %li central "
1540 "wavelength not within band.\n" ,k);
1541 fprintf(ioQQQ,"band: %s %.3e %li %li \n",
1542 continuum.chContBandLabels[k],
1543 continuum.ContBandWavelength[k],
1544 continuum.ipContBandHi[k] ,
1545 continuum.ipContBandLow[k]);
1547 }
1548
1549 /* get continuum index - RYDLAM is 911.6A = 1 Ryd so 1e4 converts
1550 * micron to Angstrom - xHi is high energy (not wavelength)
1551 * edge of the band */
1552 continuum.ipContBandHi[k] = ipoint( RYDLAM / xHi );
1553 continuum.ipContBandLow[k] = ipoint( RYDLAM / xLow );
1554
1555 /* fraction of first and last bin to include */
1556 continuum.BandEdgeCorrLow[k] =
1557 ((rfield.anu[continuum.ipContBandLow[k]-1]+
1558 rfield.widflx[continuum.ipContBandLow[k]-1]/2.f)-
1559 (realnum)(RYDLAM/xLow)) /
1560 rfield.widflx[continuum.ipContBandLow[k]-1];
1561 ASSERT( continuum.BandEdgeCorrLow[k]>=0. && continuum.BandEdgeCorrLow[k]<=1.);
1562 continuum.BandEdgeCorrHi[k] = ((realnum)(RYDLAM/xHi) -
1563 (rfield.anu[continuum.ipContBandHi[k]-1]-
1564 rfield.widflx[continuum.ipContBandHi[k]-1]/2.f)) /
1565 rfield.widflx[continuum.ipContBandHi[k]-1];
1566 ASSERT( continuum.BandEdgeCorrHi[k]>=0. && continuum.BandEdgeCorrHi[k]<=1.);
1567 /*fprintf(ioQQQ,"DEBUG bands_continuum %s %.3e %li %li \n",
1568 continuum.chContBandLabels[k],
1569 continuum.ContBandWavelength[k],
1570 continuum.ipContBandHi[k] ,
1571 continuum.ipContBandLow[k]);*/
1572
1573 if( trace.lgTrace && trace.lgConBug )
1574 {
1575 if( k==0 )
1576 fprintf( ioQQQ, " ContCreatePointer trace bands\n");
1577 fprintf( ioQQQ,
1578 " band %ld label %s low wl= %.3e low ipnt= %li "
1579 " hi wl= %.3e hi ipnt= %li \n",
1580 k,
1581 continuum.chContBandLabels[k] ,
1582 xLow,
1583 continuum.ipContBandLow[k] ,
1584 xHi,
1585 continuum.ipContBandHi[k] );
1586 }
1587# if 0
1588 // hazy table giving band properties
1589# include "prt.h"
1590 fprintf(ioQQQ,
1591 "DEBUG %s & ",
1592 continuum.chContBandLabels[k] );
1593 prt_wl( ioQQQ , continuum.ContBandWavelength[k] );
1594 fprintf(ioQQQ," & ");
1595 prt_wl( ioQQQ , xHi );
1596 fprintf(ioQQQ," -- ");
1597 prt_wl( ioQQQ , xLow );
1598 fprintf(ioQQQ,"\\\\ \n");
1599# endif
1600 ++k;
1601 }
1602 }
1603 /* now validate this incoming data */
1604 for( i=0; i<continuum.nContBand; ++i )
1605 {
1606 /* make sure all are positive */
1607 if( continuum.ContBandWavelength[i] <=0. )
1608 {
1609 fprintf( ioQQQ, " ContBandWavelength band %li has non-positive entry.\n",i );
1611 }
1612 }
1613
1614 fclose(ioDATA);
1615 return;
1616}
double atmdat_2phot_shapefunction(double EbyE2nu, long ipISO, long nelem)
void FeIIPoint(void)
t_atoms atoms
Definition atoms.cpp:5
FILE * ioQQQ
Definition cddefines.cpp:7
const int FILENAME_PATH_LENGTH_2
Definition cddefines.h:249
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
#define MALLOC(exp)
Definition cddefines.h:501
const int ipOXYGEN
Definition cddefines.h:312
const int LIMELM
Definition cddefines.h:258
#define EXIT_FAILURE
Definition cddefines.h:140
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition service.cpp:381
#define cdEXIT(FAIL)
Definition cddefines.h:434
const int NISO
Definition cddefines.h:261
const int ipHELIUM
Definition cddefines.h:306
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition service.cpp:70
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void fixit(void)
Definition service.cpp:991
static bool lgCalled
Definition cddrive.cpp:425
long nWindLine
Definition cdinit.cpp:19
EmissionProxy::iterator iterator
Definition emission.h:317
long int & ipFine() const
Definition emission.h:413
realnum & Aul() const
Definition emission.h:613
static t_fe2ovr_la & Inst()
Definition cddefines.h:175
TransitionProxy::iterator iterator
Definition transition.h:280
long & ipCont() const
Definition transition.h:450
EmissionList::reference Emis() const
Definition transition.h:408
realnum ph1(int i, int j, int k, int l) const
Definition atmdat.h:329
void init_pointers()
void set_ipoint(long n, long val)
Definition yield.h:73
int nlines() const
Definition yield.h:76
vector< realnum > As2nu
Definition two_photon.h:46
long ipTwoPhoE
Definition two_photon.h:41
STATIC void fiddle(long int ipnt, double exact)
STATIC void ipShells(long int nelem)
STATIC long LimitSh(long int ion, long int nshell, long int nelem)
STATIC void ContBandsCreate(const char chFile[])
void ContCreatePointers(void)
long ipFineCont(double energy_ryd)
long ipoint(double energy_ryd)
long ipLineEnergy(double energy, const char *chLabel, long ipIonEnergy)
long ipContEnergy(double energy, const char *chLabel)
t_continuum continuum
Definition continuum.cpp:5
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition cpu.cpp:625
const realnum SMALLFLOAT
Definition cpu.h:191
t_dense dense
Definition dense.cpp:24
t_elementnames elementnames
t_fe fe
Definition fe.cpp:5
vector< diatomics * > diatoms
Definition h2.cpp:8
vector< diatomics * >::iterator diatom_iter
Definition h2.h:13
t_he he
Definition he.cpp:5
t_Heavy Heavy
Definition heavy.cpp:5
t_hmi hmi
Definition hmi.cpp:5
t_ionbal ionbal
Definition ionbal.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
t_isoCTRL iso_ctrl
Definition iso.cpp:6
const int ipHe1s1S
Definition iso.h:41
const int ipH1s
Definition iso.h:27
const int ipHe2s3S
Definition iso.h:44
const int ipHE_LIKE
Definition iso.h:63
const int ipH2p
Definition iso.h:29
const int ipH2s
Definition iso.h:28
const int ipH_LIKE
Definition iso.h:62
void iso_create(void)
double eina(double gf, double enercm, double gup)
double abscf(double gf, double enercm, double gl)
double GetGF(double trans_prob, double enercm, double gup)
t_opac opac
Definition opacity.cpp:5
t_oxy oxy
Definition oxy.cpp:5
UNUSED const double PI
Definition physconst.h:29
UNUSED const double PI4
Definition physconst.h:35
UNUSED const double SPEEDLIGHT
Definition physconst.h:100
UNUSED const double ELECTRON_MASS
Definition physconst.h:91
UNUSED const double EN1RYD
Definition physconst.h:179
UNUSED const double EVRYD
Definition physconst.h:189
UNUSED const double RYDLAM
Definition physconst.h:176
void prt_wl(FILE *ioOUT, realnum wl)
Definition prt.cpp:13
t_rfield rfield
Definition rfield.cpp:8
const double WL_B_FILT
Definition rfield.h:14
const double WL_V_FILT
Definition rfield.h:11
t_rt rt
Definition rt.cpp:5
t_secondaries secondaries
static double * g
Definition species2.cpp:28
long int nSpecies
Definition taulines.cpp:21
vector< vector< TransitionList > > SatelliteLines
Definition taulines.cpp:38
vector< qList > dBaseStates
Definition taulines.cpp:15
vector< vector< TransitionList > > ExtraLymanLines
Definition taulines.cpp:25
TransitionList UTALines("UTALines", &AnonStates)
TransitionList TauLine2("TauLine2", &AnonStates)
vector< TransitionList > dBaseTrans
Definition taulines.cpp:17
TransitionProxy::iterator TauDummy
Definition taulines.cpp:60
TransitionList HFLines("HFLines", &AnonStates)
bool lgStatesAdded
Definition taulines.cpp:10
bool lgLinesAdded
Definition taulines.cpp:11
multi_arr< int, 3 > ipSatelliteLines
Definition taulines.cpp:37
long int nUTA
Definition taulines.cpp:26
long int nLevel1
Definition taulines.cpp:28
long int nHFLines
Definition taulines.cpp:31
multi_arr< int, 3 > ipExtraLymanLines
Definition taulines.cpp:24
TransitionList TauLines("TauLines", &AnonStates)
void checkTransitionListOfLists(vector< TransitionList > &list)
Definition taulines.cpp:42
vector< TransitionList > AllTransitions
Definition taulines.cpp:8
t_trace trace
Definition trace.cpp:5
void chIonLbl(char *chIonLbl_v, const TransitionProxy &t)
char * chLineLbl(const TransitionProxy &t)
void TwoPhotonSetup(vector< two_photon > &tnu_vec, const long &ipHi, const long &ipLo, const double &Aul, const TransitionProxy &tr, const long ipISO, const long nelem)
Definition two_photon.cpp:9