root/branches/newmole/source/cont_createpointers.cpp

Revision 2582, 57.7 kB (checked in by rporter, 2 weeks ago)

branches/newmole - bring in trunk r2541:2581

  • Property svn:eol-style set to native
Line 
1/* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
2 * others.  For conditions of distribution and use see copyright notice in license.txt */
3/*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 "predcont.h"
32#include "atmdat.h"
33#include "ipoint.h"
34#include "h2.h"
35#include "continuum.h"
36
37/*LimitSh sets upper energy limit to subshell integrations */
38STATIC long LimitSh(long int ion,
39  long int nshell,
40  long int nelem);
41
42STATIC void ipShells(
43        /* nelem is the atomic number on the C scale, Li is 2 */
44        long int nelem);
45
46/*ContBandsCreate - read set of continuum bands to enter total emission into line*/
47STATIC void ContBandsCreate(
48        /* chFile is optional filename, if void then use default bands,
49         * if not void then use file specified,
50         * return value is 0 for success, 1 for failure */
51         const char chFile[] );
52
53/* upper level for two-photon emission in H and He iso sequences */
54#define TwoS    (1+ipISO)
55
56/*fiddle adjust energy bounds of certain cells so that match ionization edges exactly */
57STATIC void fiddle(long int ipnt,
58  double exact);
59
60void ContCreatePointers(void)
61{
62        char chLab[5];
63        long int 
64          i,
65          ion,
66          ipHi,
67          ipLo,
68          ipISO,
69          iWL_Ang,
70          j,
71          nelem,
72          nshells;
73        double energy,
74                xnew;
75        /* counter to say whether pointers have ever been evaluated */
76        static int nCalled = 0;
77
78        DEBUG_ENTRY( "ContCreatePointers()" );
79
80        /* create the hydrogen atom for this core load, routine creates space then zeros it out
81         * on first call, on second and later calls it only zeros things out */
82        iso_create();
83
84        /* create internal static variables needed to do the H2 molecule */
85        H2_Create();
86
87        /* nCalled is local static variable defined 0 when defined.
88         * it is incremented below, so that space only allocated one time per coreload. */
89        if( nCalled > 0 )
90        {
91                if( trace.lgTrace )
92                        fprintf( ioQQQ, " ContCreatePointers called, not evaluating.\n" );
93                return;
94        }
95        else
96        {
97                if( trace.lgTrace )
98                        fprintf( ioQQQ, " ContCreatePointers called first time.\n" );
99                ++nCalled;
100        }
101
102        for( i=0; i < rfield.nupper; i++ )
103        {
104                /* this is array of labels for lines and continua, set to blanks at first */
105                strcpy( rfield.chContLabel[i], "    ");
106                strcpy( rfield.chLineLabel[i], "    ");
107        }
108
109        /* we will generate a set of array indices to ionization edges for
110         * the first thirty elements.  First set all array indices to
111         * totally bogus values so we will crash if misused */
112        for( nelem=0; nelem<LIMELM; ++nelem )
113        {
114                if( dense.lgElmtOn[nelem] )
115                {
116                        for( ion=0; ion<LIMELM; ++ion )
117                        {
118                                for( nshells=0; nshells<7; ++nshells )
119                                {
120                                        for( j=0; j<3; ++j )
121                                        {
122                                                opac.ipElement[nelem][ion][nshells][j] = INT_MIN;
123                                        }
124                                }
125                        }
126                }
127        }
128
129        /* pointer to excited state of O+2 */
130        opac.ipo3exc[0] = ipContEnergy(3.855,"O3ex");
131
132        /* main hydrogenic arrays - THIS OCCURS TWICE!! H and He here, then the
133         * remaining hydrogenic species near the bottom.  This is so that H and HeII get
134         * their labels stuffed into the arrays, and the rest of the hydrogenic series
135         * get whatever is left over after the level 1 lines.
136         * to find second block, search on "ipZ=2" */
137        /* NB note that no test for H or He being on exists here - we will always
138         * define the H and He arrays even when He is off, since we need to
139         * know where the he edges are for the bookkeeping that occurs in continuum
140         * binning routines */
141        /* this loop is over H, He-like only - DO NOT CHANGE TO NISO */
142        for( ipISO=ipH_LIKE; ipISO<=ipHE_LIKE; ++ipISO )
143        {
144                /* this will be over HI, HeII, then HeI only */
145                for( nelem=ipISO; nelem < 2; nelem++ )
146                {
147                        /* generate label for this ion */
148                        sprintf( chLab, "%2s%2ld",elementnames.chElementSym[nelem], nelem+1-ipISO );
149
150                        /* array index for continuum edges for ground */
151                        iso.ipIsoLevNIonCon[ipISO][nelem][0] = ipContEnergy(iso.xIsoLevNIonRyd[ipISO][nelem][0],chLab);
152                        for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
153                        {
154                                /* array index for continuum edges for excited levels */
155                                iso.ipIsoLevNIonCon[ipISO][nelem][ipHi] = ipContEnergy(iso.xIsoLevNIonRyd[ipISO][nelem][ipHi],chLab);
156
157                                /* define all line array indices */
158                                for( ipLo=0; ipLo < ipHi; ipLo++ )
159                                {
160                                        if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA )
161                                                continue;
162
163                                        /* some lines have negative or zero energy */
164                                        /* >>chng 03 apr 22, this check was if less than or equal to zero,
165                                         * changed to lowest energy point so that ultra low energy transitions are
166                                         * not considered.      */
167                                        if( Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN * WAVNRYD < continuum.filbnd[0] )
168                                                continue;
169
170                                        /* some energies are negative for inverted levels */
171                                        Transitions[ipISO][nelem][ipHi][ipLo].ipCont =
172                                                ipLineEnergy(Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN * WAVNRYD , chLab ,
173                                                iso.ipIsoLevNIonCon[ipISO][nelem][ipLo]);
174                                        Transitions[ipISO][nelem][ipHi][ipLo].Emis->ipFine =
175                                                ipFineCont(Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN * WAVNRYD );
176                                        /* check that energy scales are the same, to within energy resolution of arrays */
177#                                       ifndef NDEBUG
178                                        if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->ipFine > 0 )
179                                        {
180                                                realnum anuCoarse = rfield.anu[Transitions[ipISO][nelem][ipHi][ipLo].ipCont-1];
181                                                realnum anuFine = rfield.fine_anu[Transitions[ipISO][nelem][ipHi][ipLo].Emis->ipFine];
182                                                realnum widCoarse = rfield.widflx[Transitions[ipISO][nelem][ipHi][ipLo].ipCont-1];
183                                                realnum widFine = anuFine - rfield.fine_anu[Transitions[ipISO][nelem][ipHi][ipLo].Emis->ipFine-1];
184                                                realnum width = MAX2( widFine , widCoarse );
185                                                /* NB - can only assert more than width of coarse cell
186                                                 * since close to ionization limit, coarse index is
187                                                 * kept below next ionization edge
188                                                 * >>chng 05 mar 02, pre factor below had been 1.5, chng to 2
189                                                 * tripped near H grnd photo threshold */
190                                                ASSERT( fabs(anuCoarse - anuFine) / anuCoarse <
191                                                        2.*width/anuCoarse);
192                                        }
193#                                       endif
194                                }
195                        }/* ipHi loop */
196                }/* nelem loop */
197        }/* ipISO */
198
199        /* need to increase the cell for the HeII balmer continuum by one, so that
200         * hydrogen Lyman continuum will pick it up. */
201        nelem = ipHELIUM;
202        /* copy label over to next higher cell */
203        if( strcmp( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s]-1] , "He 2" ) == 0)
204        {
205                strcpy( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s]],
206                                 rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s]-1] );
207                /* set previous spot to blank */
208                strcpy( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s]-1] , "    ");
209        }
210        else if( strcmp( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s]-1] , "H  1" ) == 0)
211        {
212                /* set previous spot to blank */
213                strcpy( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s]] , "He 2");
214        }
215        else
216        {
217                fprintf(ioQQQ," insanity heii pointer fix contcreatepointers\n");
218        }
219        /* finally increment the two HeII pointers so that they are above the Lyman continuum */
220        ++iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s];
221        ++iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2p];
222
223        /* we will define an array of either 1 or 0 to show whether photooelectron
224         * energy is large enough to produce secondary ionizations
225         * below 100eV no secondary ionization - this is the threshold*/
226        secondaries.ipSecIon = ipoint(7.353);
227
228        /* this is highest energy where k-shell opacities are counted
229         * can be adjusted with "set kshell" command */
230        continuum.KshellLimit = ipoint(continuum.EnergyKshell);
231
232        /* pointers for molecules
233         * H2+ dissociation energy 2.647 eV but cs small below 0.638 Ryd */
234        opac.ih2pnt[0] = ipContEnergy(0.502,"H2+ ");
235        opac.ih2pnt[1] = ipoint(1.03);
236
237        /* pointers for most prominent PAH features
238         * energies given to ipContEnergy are only to put lave in the right place
239         * wavelengths are rough observed values of blends
240         * 7.6 microns */
241        i = ipContEnergy(0.0117, "PAH " );
242
243        /* feature near 6.2 microns */
244        i = ipContEnergy(0.0147, "PAH " );
245
246        /* 3.3 microns */
247        i = ipContEnergy(0.028, "PAH " );
248
249        /* 11.2 microns */
250        i = ipContEnergy(0.0080, "PAH " );
251
252        /* 12.3 microns */
253        i = ipContEnergy(0.0077, "PAH " );
254
255        /* 13.5 microns */
256        i = ipContEnergy(0.0069, "PAH " );
257
258
259        /* fix pointers for hydrogen and helium */
260
261        /* pointer to Bowen 374A resonance line */
262        he.ip374 = ipLineEnergy(1.92,"He 2",0);
263        he.ip660 = ipLineEnergy(1.38,"He 2",0);
264
265        /* set up series of continuum pointers to be used in routine lines to
266         * enter continuum points into line array */
267        for( i=0; i < NPREDCONT; i++ )
268        {
269                /* EnrPredCont contains array of energy points, as set in zerologic.c */
270                ipPredCont[i] = ipoint(EnrPredCont[i]) - 1;
271        }
272
273        /* pointer to energy defining effect x-ray column */
274        rt.ipxry = ipoint(73.5);
275
276        /* pointer to Hminus edge at 0.754eV */
277        hmi.iphmin = ipContEnergy(0.05544,"H-  ");
278
279        /* pointer to threshold for H2 photoionization at 15.4 eV */
280        opac.ipH2_photo_thresh = ipContEnergy( 15.4/EVRYD , "H2  ");
281
282        hmi.iheh1 = ipoint(1.6);
283        hmi.iheh2 = ipoint(2.3);
284
285        /* pointer to carbon k-shell ionization */
286        opac.ipCKshell = ipoint(20.6);
287
288        /* pointer to threshold for pair production */
289        opac.ippr = ipoint(7.51155e4) + 1;
290
291        /* pointer to x-ray - gamma ray bound; 100 kev */
292        rfield.ipEnerGammaRay = ipoint(rfield.EnerGammaRay);
293
294        t_fe2ovr_la::Inst().init_pointers();
295
296        /* make low energy edges of some cells exact,
297         * index is on c scale
298         * 0.99946 is correct energy of hydrogen in inf mass ryd */
299        fiddle(iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1,0.99946);
300        fiddle(iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2p]-1,0.24986);
301        /* confirm that labels are in correct location */
302        ASSERT( strcmp( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1], "H  1" ) ==0 );
303        ASSERT( strcmp( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2p]-1], "H  1" ) ==0 );
304
305        fiddle(iso.ipIsoLevNIonCon[ipH_LIKE][ipHELIUM][ipH1s]-1,4.00000);
306        ASSERT( strcmp( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][ipHELIUM][ipH1s]-1], "He 2" ) ==0 );
307
308        /* pointer to excited state of O+2 */
309        fiddle(opac.ipo3exc[0]-1,3.855);
310
311        /* now fix widflx array so that it is correct */
312        for( i=1; i<rfield.nupper-1; ++i )
313        {
314                rfield.widflx[i] = ((rfield.anu[i+1] - rfield.anu[i]) + (rfield.anu[i] - rfield.anu[i-1]))/2.f;
315        }
316
317        /* these are indices for centers of B and V filters,
318         * taken from table on page 202 of Allen, AQ, 3rd ed */
319        /* the B filter array offset */
320        rfield.ipB_filter = ipoint( RYDLAM / WL_B_FILT );
321        /* the V filter array offset */
322        rfield.ipV_filter = ipoint( RYDLAM / WL_V_FILT );
323
324        /* these are the lower and upper bounds for the G0 radiation field
325         * used by Tielens & Hollenbach in their PDR work */
326        rfield.ipG0_TH85_lo =  ipoint(  6.0 / EVRYD );
327        rfield.ipG0_TH85_hi =  ipoint( 13.6 / EVRYD );
328
329        /* this is the limits for Draine & Bertoldi Habing field */
330        rfield.ipG0_DB96_lo =  ipoint(  RYDLAM / 1110. );
331        rfield.ipG0_DB96_hi =  ipoint( RYDLAM / 912. );
332
333        /* this is special form of G0 that could be used in some future version, for now,
334         * use default, TH85 */
335        rfield.ipG0_spec_lo = ipoint(  6.0 / EVRYD );
336        rfield.ipG0_spec_hi = ipoint( RYDLAM / 912. );
337
338        /* this index is to 1000A to obtain the extinction at 1000A */
339        rfield.ip1000A = ipoint( RYDLAM / 1000. );
340
341        /* now save current form of array */
342        for( i=0; i < rfield.nupper; i++ )
343        {
344                rfield.AnuOrg[i] = rfield.anu[i];
345                rfield.anusqr[i] = (realnum)sqrt(rfield.AnuOrg[i]);
346        }
347
348        /* following order of elements is in roughly decreasing abundance
349         * when ipShells gets a cell for the valence IP it does so through
350         * ipContEnergy, which makes sure that no two ions get the same cell
351         * earliest elements have most precise ip mapping */
352
353        /* set up shell pointers for hydrogen */
354        nelem = ipHYDROGEN;
355        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.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s];
362        opac.ipElement[nelem][ion][0][0] = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s];
363
364        /* upper limit to energy integration */
365        opac.ipElement[nelem][ion][0][1] = rfield.nupper;
366
367        if( dense.lgElmtOn[ipHELIUM] )
368        {
369                /* set up shell pointers for helium */
370                nelem = ipHELIUM;
371                ion = 0;
372
373                /* the number of shells */
374                Heavy.nsShells[nelem][0] = 1;
375
376                /*pointer to ionization threshold in energy array*/
377                Heavy.ipHeavy[nelem][ion] = iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0];
378                opac.ipElement[nelem][ion][0][0] = iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0];
379
380                /* upper limit to energy integration */
381                opac.ipElement[nelem][ion][0][1] = rfield.nupper;
382
383                /* (hydrogenic) helium ion */
384                ion = 1;
385                /* the number of shells */
386                Heavy.nsShells[nelem][1] = 1;
387
388                /*pointer to ionization threshold in energy array*/
389                Heavy.ipHeavy[nelem][ion] = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s];
390                opac.ipElement[nelem][ion][0][0] = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s];
391
392                /* upper limit to energy integration */
393                opac.ipElement[nelem][ion][0][1] = rfield.nupper;
394        }
395
396        /* check that ionization potential of neutral carbon valence shell is
397         * positive */
398        ASSERT( t_ADfA::Inst().ph1(2,5,5,0) > 0. );
399
400        /* now fill in all sub-shell ionization array indices for elements heavier than He,
401         * this must be done after previous loop on iso.ipIsoLevNIonCon[ipH_LIKE] since hydrogenic species use
402         * iso.ipIsoLevNIonCon[ipH_LIKE] rather than ipoint in getting array index within continuum array */
403        for( i=NISO; i<LIMELM; ++i )
404        {
405                if( dense.lgElmtOn[i])
406                {
407                        /* i is the atomic number on the c scale, 2 for Li */
408                        ipShells(i);
409                }
410        }
411
412        /* most of these are set in ipShells, but not h-like or he-like, so do these here */
413        Heavy.Valence_IP_Ryd[ipHYDROGEN][0] = t_ADfA::Inst().ph1(0,0,ipHYDROGEN,0)/EVRYD* 0.9998787;
414        Heavy.Valence_IP_Ryd[ipHELIUM][0] = t_ADfA::Inst().ph1(0,1,ipHELIUM,0)/EVRYD* 0.9998787;
415        Heavy.Valence_IP_Ryd[ipHELIUM][1] = t_ADfA::Inst().ph1(0,0,ipHELIUM,0)/EVRYD* 0.9998787;
416        for( nelem=2; nelem<LIMELM; ++nelem )
417        {
418                Heavy.Valence_IP_Ryd[nelem][nelem-1] = t_ADfA::Inst().ph1(0,1,nelem,0)/EVRYD* 0.9998787;
419                Heavy.Valence_IP_Ryd[nelem][nelem] = t_ADfA::Inst().ph1(0,0,nelem,0)/EVRYD* 0.9998787;
420                if( dense.lgElmtOn[nelem])
421                {
422                        /* now confirm that all are properly set */
423                        for( j=0; j<=nelem; ++j )
424                        {
425                                ASSERT( Heavy.Valence_IP_Ryd[nelem][j] > 0.05 );
426                        }
427                        for( j=0; j<nelem; ++j )
428                        {
429                                ASSERT( Heavy.Valence_IP_Ryd[nelem][j] < Heavy.Valence_IP_Ryd[nelem][j+1]);
430                        }
431                }
432        }
433
434        /* array indices to bound Compton electron recoil ionization of all elements */
435        for( nelem=0; nelem<LIMELM; ++nelem )
436        {
437                if( dense.lgElmtOn[nelem])
438                {
439                        for( ion=0; ion<nelem+1; ++ion )
440                        {
441                                /* this is the threshold energy to Compton ionize valence shell electrons */
442                                energy = sqrt( Heavy.Valence_IP_Ryd[nelem][ion] * EN1RYD * ELECTRON_MASS * SPEEDLIGHT * SPEEDLIGHT ) / EN1RYD;
443                                /* the array index for this energy */
444                                ionbal.ipCompRecoil[nelem][ion] = ipoint( energy );
445                        }
446                }
447        }
448
449        /* oxygen pointers for excited states
450         * IO3 is pointer to O++ exc state, is set above */
451        oxy.i2d = ipoint(1.242);
452        oxy.i2p = ipoint(1.367);
453        opac.ipo1exc[0] = ipContEnergy(0.856,"O1ex");
454        opac.ipo1exc[1] = ipoint(2.0);
455
456        /* upper limit for excited state photoionization
457         * do not use ipContEnergy since only upper limit */
458        opac.ipo3exc[1] = ipoint(5.0);
459
460        /* upper level of 4363 */
461        opac.ipo3exc3[0] = ipContEnergy(3.646,"O3ex");
462        opac.ipo3exc3[1] = ipoint(5.0);
463
464        /* following are various pointers for OI - Ly-beta pump problem
465         * these are delta energies for Boltzmann factors */
466        /** \todo       2       this is redundant with contents of oxygen line arrays
467         * use them instead
468         * when removing this, make sure all line intensity predictions
469         * also go into oi line arrays */
470        atoms.ipoiex[0] = ipoint(0.7005);
471        atoms.ipoiex[1] = ipoint(0.10791);
472        atoms.ipoiex[2] = ipoint(0.06925);
473        atoms.ipoiex[3] = ipoint(0.01151);
474        atoms.ipoiex[4] = ipoint(0.01999);
475
476        /* >>chng 97 jan 27, move nitrogen after oxygen so that O gets the
477         * most accurate pointers
478         * Nitrogen
479         * in1(1) is thresh for photo from excited state */
480        opac.in1[0] = ipContEnergy(0.893,"N1ex");
481
482        /* upper limit */
483        opac.in1[1] = ipoint(2.);
484
485        if( (trace.lgTrace && trace.lgConBug) || (trace.lgTrace && trace.lgPointBug) )
486        {
487                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",
488                  rfield.nupper,
489                  iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s],
490                  iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][ipH1s],
491                  iso.ipIsoLevNIonCon[ipH_LIKE][ipHELIUM][ipH1s],
492                  opac.ipo3exc[0],
493                  opac.ipCKshell,
494                  ionbal.ipCompRecoil[ipHYDROGEN][0] );
495
496
497                fprintf( ioQQQ, "   ContCreatePointers: ipEnerGammaRay: %5ld IPPRpari produc%5ld\n",
498                  rfield.ipEnerGammaRay, opac.ippr );
499
500                fprintf( ioQQQ, "   ContCreatePointers: H pointers;" );
501                for( i=0; i <= 6; i++ )
502                {
503                        fprintf( ioQQQ, "%4ld%4ld", i, iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][i] );
504                }
505                fprintf( ioQQQ, "\n" );
506
507                fprintf( ioQQQ, "   ContCreatePointers: Oxy pnters;" );
508
509                for( i=1; i <= 8; i++ )
510                {
511                        fprintf( ioQQQ, "%4ld%4ld", i, Heavy.ipHeavy[ipOXYGEN][i-1] );
512                }
513                fprintf( ioQQQ, "\n" );
514        }
515
516        /* Magnesium
517         * following is energy for phot of MG+ from exc state producing 2798 */
518        opac.ipmgex = ipoint(0.779);
519
520        /* lower, upper edges of Ca+ excited term photoionizaiton */
521        opac.ica2ex[0] = ipContEnergy(0.72,"Ca2x");
522        opac.ica2ex[1] = ipoint(1.);
523
524        /* set up factors and pointers for Fe continuum fluorescence */
525        fe.ipfe10 = ipoint(2.605);
526
527        /* following is WL(CM)**2/(8PI) * conv fac for RYD to NU *A21 */
528        fe.pfe10 = (realnum)(2.00e-18/rfield.widflx[fe.ipfe10-1]);
529
530        /* this is 353 pump, f=0.032 */
531        fe.pfe11a = (realnum)(4.54e-19/rfield.widflx[fe.ipfe10-1]);
532
533        /* this is 341.1 f=0.012 */
534        fe.pfe11b = (realnum)(2.75e-19/rfield.widflx[fe.ipfe10-1]);
535        fe.pfe14 = (realnum)(1.15e-18/rfield.widflx[fe.ipfe10-1]);
536
537        /* set up energy pointers for line optical depth arrays
538         * this also increments flux, sets other parameters for lines */
539
540        /* level 1 heavy elements line array */
541        for( i=1; i <= nLevel1; i++ )
542        {
543                /* put null terminated line label into chLab using line array*/
544                chIonLbl(chLab, &TauLines[i]);
545
546                TauLines[i].ipCont =
547                        ipLineEnergy(TauLines[i].EnergyWN * WAVNRYD, chLab ,0);
548                TauLines[i].Emis->ipFine =
549                        ipFineCont(TauLines[i].EnergyWN * WAVNRYD );
550                /* for debugging pointer - index on f not c scale,
551                 * this will find all lines that entered a given cell
552                   if( TauLines[i].ipCont==603 )
553                        fprintf(ioQQQ,"( level1 %s\n", chLab);*/
554
555                if( TauLines[i].Emis->gf > 0. )
556                {
557                        /* the gf value was entered
558                         * derive the A, call to function is gf, wl (A), g_up */
559                        TauLines[i].Emis->Aul =
560                                (realnum)(eina(TauLines[i].Emis->gf,
561                          TauLines[i].EnergyWN,TauLines[i].Hi->g));
562                }
563                else if( TauLines[i].Emis->Aul > 0. )
564                {
565                        /* the Einstein A value was entered
566                         * derive the gf, call to function is A, wl (A), g_up */
567                        TauLines[i].Emis->gf =
568                                (realnum)(GetGF(TauLines[i].Emis->Aul,
569                          TauLines[i].EnergyWN,TauLines[i].Hi->g));
570                }
571                else
572                {
573                        fprintf( ioQQQ, " level 1 line does not have valid gf or A\n" );
574                        fprintf( ioQQQ, " This is ContCreatePointers\n" );
575                        cdEXIT(EXIT_FAILURE);
576                }
577
578                /* used to get damping constant */
579                TauLines[i].Emis->dampXvel =
580                        (realnum)(TauLines[i].Emis->Aul / TauLines[i].EnergyWN/PI4);
581
582                /* derive the abs coefficient, call to function is gf, wl (A), g_low */
583                TauLines[i].Emis->opacity =
584                        (realnum)(abscf(TauLines[i].Emis->gf,
585                  TauLines[i].EnergyWN,TauLines[i].Lo->g));
586                /*fprintf(ioQQQ,"TauLinesss\t%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
587                        i,TauLines[i].Emis->opacity , TauLines[i].Emis->gf , TauLines[i].Emis->Aul ,TauLines[i].EnergyWN,TauLines[i].Lo->g);*/
588
589                /* excitation energy of transition in degrees kelvin */
590                TauLines[i].EnergyK =
591                        (realnum)(T1CM)*TauLines[i].EnergyWN;
592
593                /* energy of photon in ergs */
594                TauLines[i].EnergyErg =
595                        (realnum)(ERG1CM)*TauLines[i].EnergyWN;
596
597                /*line wavelength must be gt 0 */
598                ASSERT( TauLines[i].WLAng > 0 );
599        }
600
601        /*Beginning of the atmolEmis*/
602        for( i=0; i <linesAdded2; i++)
603        {
604                atmolEmis[i].gf = (realnum)GetGF(atmolEmis[i].Aul,
605                        atmolEmis[i].tran->EnergyWN,atmolEmis[i].tran->Hi->g);
606                atmolEmis[i].dampXvel = (realnum)(atmolEmis[i].Aul/
607                                        atmolEmis[i].tran->EnergyWN/PI4);
608                atmolEmis[i].damp = -1000.0;
609                /*Put null terminated line label into chLab*/
610                strncpy(chLab,atmolEmis[i].tran->Hi->chLabel,4);
611                chLab[4]='\0';
612
613                atmolEmis[i].tran->ipCont =
614                        ipLineEnergy(atmolEmis[i].tran->EnergyWN * WAVNRYD, chLab ,0);
615                atmolEmis[i].ipFine = ipFineCont(atmolEmis[i].tran->EnergyWN * WAVNRYD );
616                /* derive the abs coefficient, call to function is gf, wl (A), g_low */
617                atmolEmis[i].opacity =
618                        (realnum)(abscf(atmolEmis[i].gf,atmolEmis[i].tran->EnergyWN,
619                                atmolEmis[i].tran->Lo->g));
620                /* excitation energy of in degrees kelvin */
621                atmolEmis[i].tran->EnergyK =
622                        (realnum)(T1CM)*atmolEmis[i].tran->EnergyWN;
623                /* energy of photon in ergs */
624                atmolEmis[i].tran->EnergyErg =
625                        (realnum)(ERG1CM)*atmolEmis[i].tran->EnergyWN ;                                                           
626        }
627        /*end of the atmolEmis*/
628
629        /* set the ipCont struc element for the H2 molecule */
630        H2_ContPoint();
631
632        for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
633        {
634                /* do remaining part of the iso sequences */
635                for( nelem=2; nelem < LIMELM; nelem++ )
636                {
637                        if( dense.lgElmtOn[nelem])
638                        {
639                                /* generate label for this ion */
640                                sprintf( chLab, "%2s%2ld",elementnames.chElementSym[nelem], nelem+1-ipISO );
641                                /* array index for continuum edges */
642                                iso.ipIsoLevNIonCon[ipISO][nelem][0] =
643                                        ipContEnergy(iso.xIsoLevNIonRyd[ipISO][nelem][0],chLab);
644
645                                for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
646                                {
647                                        /* array index for continuum edges */
648                                        iso.ipIsoLevNIonCon[ipISO][nelem][ipHi] = ipContEnergy(iso.xIsoLevNIonRyd[ipISO][nelem][ipHi],chLab);
649
650                                        /* define all line pointers */
651                                        for( ipLo=0; ipLo < ipHi; ipLo++ )
652                                        {
653
654                                                if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA )
655                                                        continue;
656
657                                                /* some lines have negative or zero energy */
658                                                /* >>chng 03 apr 22, this check was if less than or equal to zero,
659                                                 * changed to lowest energy point so that ultra low energy transitions are
660                                                 * not considered.      */
661                                                if( Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN * WAVNRYD < continuum.filbnd[0] )
662                                                        continue;
663
664                                                Transitions[ipISO][nelem][ipHi][ipLo].ipCont =
665                                                        ipLineEnergy(Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN * WAVNRYD , chLab ,
666                                                        iso.ipIsoLevNIonCon[ipISO][nelem][ipLo]);
667                                                Transitions[ipISO][nelem][ipHi][ipLo].Emis->ipFine =
668                                                        ipFineCont(Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN * WAVNRYD );
669                                        }
670                                }
671                                iso.ipIsoLevNIonCon[ipISO][nelem][0] = ipContEnergy(iso.xIsoLevNIonRyd[ipISO][nelem][0],chLab);
672                        }
673                }
674        }
675        for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
676        {
677                /* this will be over HI, HeII, then HeI only */
678                for( nelem=ipISO; nelem < LIMELM; nelem++ )
679                {
680                        if( dense.lgElmtOn[nelem])
681                        {
682                                ipLo = 0;
683                                /* these are the extra Lyman lines */
684                                for( ipHi=2; ipHi < iso.nLyman[ipISO]; ipHi++ )
685                                {
686                                        /* some energies are negative for inverted levels */
687                                        /*lint -e772 chLab not initialized */
688                                        ExtraLymanLines[ipISO][nelem][ipHi].ipCont =
689                                                ipLineEnergy(ExtraLymanLines[ipISO][nelem][ipHi].EnergyWN * WAVNRYD , chLab ,
690                                                iso.ipIsoLevNIonCon[ipISO][nelem][ipLo]);
691
692                                        ExtraLymanLines[ipISO][nelem][ipHi].Emis->ipFine =
693                                                ipFineCont(ExtraLymanLines[ipISO][nelem][ipHi].EnergyWN * WAVNRYD );
694                                        /*lint +e772 chLab not initialized */
695                                }
696
697                                if( iso.lgDielRecom[ipISO] )
698                                {
699                                        ASSERT( ipISO>ipH_LIKE );
700                                        for( ipHi=0; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ )
701                                        {
702                                               
703                                                SatelliteLines[ipISO][nelem][ipHi].ipCont = ipLineEnergy(
704                                                        SatelliteLines[ipISO][nelem][ipHi].EnergyWN * WAVNRYD , chLab ,
705                                                        0);
706
707                                                SatelliteLines[ipISO][nelem][ipHi].Emis->ipFine = 
708                                                        ipFineCont(SatelliteLines[ipISO][nelem][ipHi].EnergyWN * WAVNRYD );
709                                        }
710                                }
711                        }
712                }
713        }
714
715        /* some lines do not exist, the flag indicating this is ipCont == -1 */
716        /* for H-like sequence, these are 2p-2s (energies degenerate) and 2s-1s, two photon */
717        ipISO = ipHYDROGEN;
718        /* do remaining part of the he iso sequence */
719        for( nelem=ipISO; nelem < LIMELM; nelem++ )
720        {
721                if( dense.lgElmtOn[nelem])
722                {
723                        /* for H-like sequence want 2p-2s (energies degenerate) and 2s-1s, two photon */
724                        Transitions[ipISO][nelem][ipH2p][ipH2s].ipCont = -1;
725                        //Transitions[ipISO][nelem][ipH2p][ipH2s].Emis->ipFine = -1;
726                        Transitions[ipISO][nelem][ipH2s][ipH1s].ipCont = -1;
727                        //Transitions[ipISO][nelem][ipH2s][ipH1s].Emis->ipFine = -1;
728                }
729        }
730
731        fixit(); /* is this redundant? */
732        /* for He-like sequence the majority of the transitions are bogus - A set to special value in this case */
733        ipISO = ipHELIUM;
734        /* do remaining part of the he iso sequence */
735        for( nelem=ipISO; nelem < LIMELM; nelem++ )
736        {
737                if( dense.lgElmtOn[nelem])
738                {
739                        /* this is the two photon transition in the singlets */
740                        Transitions[ipISO][nelem][ipHe2s1S][ipHe1s1S].ipCont = -1;
741                        Transitions[ipISO][nelem][ipHe2s1S][ipHe1s1S].Emis->ipFine = -1;
742
743                        for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
744                        {
745                                for( ipLo=0; ipLo < ipHi; ipLo++ )
746                                {
747                                        if( Transitions[ipISO][nelem][ipHi][ipLo].ipCont <= 0 )
748                                                continue;
749
750                                        if( fabs(Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul - iso.SmallA) < SMALLFLOAT )
751                                        {
752                                                /* iso.SmallA is value assigned to bogus transitions */
753                                                Transitions[ipISO][nelem][ipHi][ipLo].ipCont = -1;
754                                                Transitions[ipISO][nelem][ipHi][ipLo].Emis->ipFine = -1;
755                                        }
756                                }
757                        }
758                }
759        }
760
761        /* co carbon monoxide line array */
762        for( i=0; i < nCORotate; i++ )
763        {
764                /* excitation energy of transition in degrees kelvin */
765                C12O16Rotate[i].EnergyK =
766                        (realnum)(T1CM)*C12O16Rotate[i].EnergyWN;
767                C13O16Rotate[i].EnergyK =
768                        (realnum)(T1CM)*C13O16Rotate[i].EnergyWN;
769
770                /* energy of photon in ergs */
771                C12O16Rotate[i].EnergyErg =
772                        (realnum)(ERG1CM)*C12O16Rotate[i].EnergyWN;
773                C13O16Rotate[i].EnergyErg =
774                        (realnum)(ERG1CM)*C13O16Rotate[i].EnergyWN;
775
776                /* put null terminated line label into chLab using line array*/
777                chIonLbl(chLab, &C12O16Rotate[i]);
778                chIonLbl(chLab, &C13O16Rotate[i]);
779
780                C12O16Rotate[i].ipCont =
781                        ipLineEnergy(C12O16Rotate[i].EnergyWN * WAVNRYD, "12CO" ,0);
782                C12O16Rotate[i].Emis->ipFine =
783                        ipFineCont(C12O16Rotate[i].EnergyWN * WAVNRYD );
784                C13O16Rotate[i].ipCont =
785                        ipLineEnergy(C13O16Rotate[i].EnergyWN * WAVNRYD, "13CO" ,0);
786                C13O16Rotate[i].Emis->ipFine =
787                        ipFineCont(C13O16Rotate[i].EnergyWN * WAVNRYD );
788
789                if( C12O16Rotate[i].Emis->gf > 0. )
790                {
791                        /* the gf value was entered
792                         * derive the A, call to function is gf, wl (A), g_up */
793                        C12O16Rotate[i].Emis->Aul =
794                                (realnum)(eina(C12O16Rotate[i].Emis->gf,
795                          C12O16Rotate[i].EnergyWN,C12O16Rotate[i].Hi->g));
796                }
797                else if( C12O16Rotate[i].Emis->Aul > 0. )
798                {
799                        /* the Einstein A value was entered
800                         * derive the gf, call to function is A, wl (A), g_up */
801                        C12O16Rotate[i].Emis->gf =
802                                (realnum)(GetGF(C12O16Rotate[i].Emis->Aul,
803                          C12O16Rotate[i].EnergyWN,C12O16Rotate[i].Hi->g));
804                }
805                else
806                {
807                        fprintf( ioQQQ, " 12CO line does not have valid gf or A\n" );
808                        fprintf( ioQQQ, " This is ContCreatePointers\n" );
809                        cdEXIT(EXIT_FAILURE);
810                }
811                if( C13O16Rotate[i].Emis->gf > 0. )
812                {
813                        /* the gf value was entered
814                         * derive the A, call to function is gf, wl (A), g_up */
815                        C13O16Rotate[i].Emis->Aul =
816                                (realnum)(eina(C13O16Rotate[i].Emis->gf,
817                          C13O16Rotate[i].EnergyWN,C13O16Rotate[i].Hi->g));
818                }
819                else if( C13O16Rotate[i].Emis->Aul > 0. )
820                {
821                        /* the Einstein A value was entered
822                         * derive the gf, call to function is A, wl (A), g_up */
823                        C13O16Rotate[i].Emis->gf =
824                                (realnum)(GetGF(C13O16Rotate[i].Emis->Aul,
825                          C13O16Rotate[i].EnergyWN,C13O16Rotate[i].Hi->g));
826                }
827                else
828                {
829                        fprintf( ioQQQ, " 13CO line does not have valid gf or A\n" );
830                        fprintf( ioQQQ, " This is ContCreatePointers\n" );
831                        cdEXIT(EXIT_FAILURE);
832                }
833
834                /*line wavelength must be gt 0 */
835                ASSERT( C12O16Rotate[i].WLAng > 0 );
836                ASSERT( C13O16Rotate[i].WLAng > 0 );
837
838                C12O16Rotate[i].Emis->dampXvel =
839                        (realnum)(C12O16Rotate[i].Emis->Aul/
840                  C12O16Rotate[i].EnergyWN/PI4);
841
842                C13O16Rotate[i].Emis->dampXvel =
843                        (realnum)(C13O16Rotate[i].Emis->Aul/
844                  C13O16Rotate[i].EnergyWN/PI4);
845
846                /* derive the abs coefficient, call to function is gf, wl (A), g_low */
847                C12O16Rotate[i].Emis->opacity =
848                        (realnum)(abscf(C12O16Rotate[i].Emis->gf,
849                  C12O16Rotate[i].EnergyWN,C12O16Rotate[i].Lo->g));
850                C13O16Rotate[i].Emis->opacity =
851                        (realnum)(abscf(C13O16Rotate[i].Emis->gf,
852                  C13O16Rotate[i].EnergyWN,C13O16Rotate[i].Lo->g));
853        }
854
855        /* inner shell transitions */
856        for( i=0; i<nUTA; ++i )
857        {
858                if( UTALines[i].Emis->Aul > 0. )
859                {
860                        UTALines[i].Emis->dampXvel =
861                                (realnum)(UTALines[i].Emis->Aul/ UTALines[i].EnergyWN/PI4);
862
863                        /* derive the abs coefficient, call to function is gf, wl (A), g_low */
864                        UTALines[i].Emis->opacity =
865                                (realnum)(abscf( UTALines[i].Emis->gf, UTALines[i].EnergyWN, UTALines[i].Lo->g));
866
867                        /* excitation energy of transition in degrees kelvin */
868                        UTALines[i].EnergyK =
869                                (realnum)(T1CM*UTALines[i].EnergyWN);
870
871                        /* energy of photon in ergs */
872                        UTALines[i].EnergyErg =
873                                (realnum)(ERG1CM*UTALines[i].EnergyWN);
874
875                        /* put null terminated line label into chLab using line array*/
876                        chIonLbl(chLab, &UTALines[i]);
877
878                        /* get pointer to energy in continuum mesh */
879                        UTALines[i].ipCont = ipLineEnergy(UTALines[i].EnergyWN * WAVNRYD , chLab,0 );
880                        UTALines[i].Emis->ipFine = ipFineCont(UTALines[i].EnergyWN * WAVNRYD  );
881
882                        /* find heating per absorption,
883                         * first find threshold for this shell in ergs */
884                        /* ionization threshold in erg */
885                        double thresh = Heavy.Valence_IP_Ryd[UTALines[i].Hi->nelem-1][UTALines[i].Hi->IonStg-1] *EN1RYD;
886                        UTALines[i].Coll.heat = (UTALines[i].EnergyErg-thresh);
887                        ASSERT( UTALines[i].Coll.heat> 0. );
888                }
889        }
890
891        /* level 2 heavy element lines */
892        for( i=0; i < nWindLine; i++ )
893        {
894
895                /* derive the A, call to function is gf, wl (A), g_up */
896                TauLine2[i].Emis->Aul =
897                        (realnum)(eina(TauLine2[i].Emis->gf,
898                  TauLine2[i].EnergyWN,TauLine2[i].Hi->g));
899
900                /* coefficient needed for damping constant - units cm s-1 */
901                TauLine2[i].Emis->dampXvel =
902                        (realnum)(TauLine2[i].Emis->Aul/
903                  TauLine2[i].EnergyWN/PI4);
904
905                /* derive the abs coefficient, call to function is gf, wl (A), g_low */
906                TauLine2[i].Emis->opacity =
907                        (realnum)(abscf(TauLine2[i].Emis->gf,
908                  TauLine2[i].EnergyWN,TauLine2[i].Lo->g));
909
910                /* excitation energy of transition in degrees kelvin */
911                TauLine2[i].EnergyK =
912                        (realnum)(T1CM*TauLine2[i].EnergyWN);
913
914                /* energy of photon in ergs */
915                TauLine2[i].EnergyErg =
916                        (realnum)(ERG1CM*TauLine2[i].EnergyWN);