root/branches/newmole/source/atmdat_readin.cpp

Revision 2582, 57.2 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/*atmdat_readin read in some data files, but only if this is very first call,
4 * called by Cloudy */
5#include "cddefines.h"
6#include "physconst.h"
7#include "taulines.h"
8#include "mewecoef.h"
9#include "iterations.h"
10#include "heavy.h"
11#include "mole.h"
12#include "yield.h"
13#include "trace.h"
14#include "lines.h"
15#include "lines_service.h"
16#include "ionbal.h"
17#include "struc.h"
18#include "geometry.h"
19#include "dynamics.h"
20#include "elementnames.h"
21#include "hyperfine.h"
22#include "atmdat.h"
23/* */
24/* this was needed to get array to crash out of bounds if not set.
25 * std limits on limits.h did not work with visual studio! */
26const long INTBIG = 2000000000;
27
28/* these are the individual pointers to the level 1 lines, they are set to
29 * very large negative. 
30 * NB NB NB!!
31 * these occur two times in the code!!
32 * They are declared in taulines.h, and defined here  */
33long    ipT1656=INTBIG,    ipT9830=INTBIG,  /*ipH21cm=INTBIG, ipHe3cm=INTBIG,*/
34        ipT8727=INTBIG,    ipT1335=INTBIG,  ipT1909=INTBIG,
35        ipT977=INTBIG,    ipT1550=INTBIG,  ipT1548=INTBIG, ipT386=INTBIG,
36        ipT310=INTBIG,    ipc31175=INTBIG, ipT291=INTBIG,  ipT280=INTBIG,
37        ipT274=INTBIG,    ipT270=INTBIG,   ipT312=INTBIG,  ipT610=INTBIG,
38        ipT370=INTBIG,    ipT157=INTBIG,   ipT1085=INTBIG,
39        ipT990=INTBIG,    ipT1486=INTBIG,  ipT765=INTBIG,  ipT1243=INTBIG,
40        ipT1239=INTBIG,   ipT374g=INTBIG,  ipT374x=INTBIG, ipT1200=INTBIG,
41        ipT2140=INTBIG,   ipT671=INTBIG,   ipT315=INTBIG,  ipT324=INTBIG,
42        ipT333=INTBIG,    ipT209=INTBIG,   ipT122=INTBIG,  ipT205=INTBIG,
43        ipT57=INTBIG,     ipT6300=INTBIG,  ipT6363=INTBIG, ipT5577=INTBIG,
44        ipT834=INTBIG,    ipT1661=INTBIG,  ipT1666=INTBIG, ipT835=INTBIG,
45        ipT789=INTBIG,    ipT630=INTBIG,   ipT1304=INTBIG, ipSi10_606=INTBIG,
46        ipT1039=INTBIG,   ipT8446=INTBIG,  ipT4368=INTBIG, ipTOI13=INTBIG,
47        ipTOI11=INTBIG,   ipTOI29=INTBIG,  ipTOI46=INTBIG, ipTO1025=INTBIG,
48        ipT304=INTBIG,    ipT1214=INTBIG,  ipT150=INTBIG,  ipT146=INTBIG,
49        ipT63=INTBIG,     ipTO88=INTBIG,   ipT52=INTBIG,   ipT26=INTBIG,
50        ipT1032=INTBIG,   ipT1037=INTBIG,  ipF0229=INTBIG, ipF0267=INTBIG,
51        ipF444=INTBIG,    ipF425=INTBIG,   ipT770=INTBIG, ipT780=INTBIG,
52        ipxNe0676=INTBIG, ipT895=INTBIG,   ipT88=INTBIG, ipTNe13=INTBIG,
53        ipTNe36=INTBIG,   ipTNe16=INTBIG,  ipTNe14=INTBIG, ipTNe24=INTBIG,
54        ipT5895=INTBIG,   ipfsNa373=INTBIG,ipfsNa490=INTBIG, ipfsNa421=INTBIG,
55        ipxNa6143=INTBIG, ipxNa6862=INTBIG,ipxNa0746=INTBIG, ipMgI2853=INTBIG,
56        ipMgI2026=INTBIG, ipT2796=INTBIG,  ipT2804=INTBIG,
57        ipT705=INTBIG,    ipT4561=INTBIG,  ipxMg51325=INTBIG, ipxMg52417=INTBIG,
58        ipxMg52855=INTBIG,ipxMg71190=INTBIG, ipxMg72261=INTBIG,
59        ipxMg72569=INTBIG,ipxMg08303=INTBIG, ipTMg610=INTBIG, ipTMg625=INTBIG,
60        ipT58=INTBIG,     ipTMg4=INTBIG, ipTMg14=INTBIG, ipTMg6=INTBIG,
61        ipfsMg790=INTBIG, ipfsMg755=INTBIG, ipAlI3957=INTBIG, ipAlI3090=INTBIG,
62        ipT1855=INTBIG,   ipT1863=INTBIG, ipT2670=INTBIG,
63        ipAl529=INTBIG,   ipAl6366=INTBIG, ipAl6912=INTBIG, ipAl8575=INTBIG,
64        ipAl8370=INTBIG,  ipAl09204=INTBIG, ipT639=INTBIG,
65        ipTAl550=INTBIG,  ipTAl568=INTBIG, ipTAl48=INTBIG, ipSii2518=INTBIG,
66        ipSii2215=INTBIG, ipT1808=INTBIG,
67        ipT1207=INTBIG,   ipT1895=INTBIG, ipT1394=INTBIG, ipT1403=INTBIG,
68        ipT1527=INTBIG,   ipT1305=INTBIG, ipT1260=INTBIG, ipSi619=INTBIG,
69        ipSi10143=INTBIG, ipTSi499=INTBIG, ipTSi521=INTBIG, ipTSi41=INTBIG,
70        ipTSi35=INTBIG,   ipTSi25=INTBIG, ipTSi65=INTBIG,
71        ipTSi3=INTBIG,    ipTSi4=INTBIG, ipP0260=INTBIG, ipP0233=INTBIG,
72        ipP0318=INTBIG,   ipP713=INTBIG, ipP848=INTBIG, ipP817=INTBIG,
73        ipP1027=INTBIG,   ipP1018=INTBIG, ipT1256=INTBIG, ipT1194=INTBIG,
74        ipTS1720=INTBIG,  ipT1198=INTBIG, ipT786=INTBIG,
75        ipT933=INTBIG,    ipT944=INTBIG, ipfsS810=INTBIG, ipfsS912=INTBIG,
76        ipfsS938=INTBIG,  ipfsS1119=INTBIG, ipfsS1114=INTBIG, ipfsS1207=INTBIG,
77        ipTSu418=INTBIG,  ipTSu446=INTBIG, ipTSu30=INTBIG, ipTS19=INTBIG,
78        ipTS34=INTBIG,    ipTS11=INTBIG, ipfsCl214=INTBIG, ipfsCl233=INTBIG,
79        ipCl04203=INTBIG, ipCl04117=INTBIG, ipCl973=INTBIG, ipCl1030=INTBIG,
80        ipCl1092=INTBIG,  ipT354=INTBIG, ipT389=INTBIG, ipT25=INTBIG,
81        ipTAr7=INTBIG,    ipTAr9=INTBIG, ipTAr22=INTBIG, ipTAr13=INTBIG,
82        ipTAr8=INTBIG,    ipAr06453=INTBIG, ipAr1055=INTBIG, ipAr1126=INTBIG,
83        ipAr1178=INTBIG,  ipKI7745=INTBIG, ipxK03462=INTBIG, ipxK04598=INTBIG,
84        ipxK04154=INTBIG, ipxK06882=INTBIG, ipxK06557=INTBIG,
85        ipxK07319=INTBIG, ipxK11425=INTBIG, ipCaI4228=INTBIG, ipT3934=INTBIG,
86        ipT3969=INTBIG,   ipT8498=INTBIG, ipT8542=INTBIG,
87        ipT8662=INTBIG,   ipT7291=INTBIG, ipT7324=INTBIG, ipTCa302=INTBIG,
88        ipTCa345=INTBIG,  ipTCa19=INTBIG, ipTCa3=INTBIG, ipTCa12=INTBIG,
89        ipTCa4=INTBIG,    ipCa0741=INTBIG, ipCa0761=INTBIG, ipCa08232=INTBIG,
90        ipCa12333=INTBIG, ipSc05231=INTBIG, ipSc13264=INTBIG,
91        ipTi06172=INTBIG, ipTi14212=INTBIG, ipVa07130=INTBIG, ipVa15172=INTBIG,
92        ipCr08101=INTBIG, ipCr16141=INTBIG, ipxMn0979=INTBIG,
93        ipxMn1712=INTBIG, ipFeI3884=INTBIG, ipFeI3729=INTBIG, ipFeI3457=INTBIG,
94        ipFeI3021=INTBIG, ipFeI2966=INTBIG, ipTuv3=INTBIG,
95        ipTr48=INTBIG,    ipTFe16=INTBIG, ipTFe26=INTBIG, ipTFe34=INTBIG,
96        ipTFe35=INTBIG,   ipTFe46=INTBIG, ipTFe56=INTBIG, ipT1122=INTBIG,
97        ipFe0795=INTBIG,  ipFe0778=INTBIG, ipT245=INTBIG, ipT352=INTBIG,
98        ipFe106375=INTBIG,ipT353=INTBIG, /*ipFe1310=INTBIG, ipFe1311=INTBIG,*/
99        ipT347=INTBIG,    ipT192=INTBIG, ipT255=INTBIG, ipT11=INTBIG,
100        ipT191=INTBIG,    /*ipTFe07=INTBIG, ipTFe61=INTBIG,*/ ipFe18975=INTBIG, ipTFe23=INTBIG,
101        ipTFe13=INTBIG,   ipCo11527=INTBIG, ipxNi1242=INTBIG;
102long    ipS4_1405=INTBIG,ipS4_1398=INTBIG,ipS4_1424=INTBIG,ipS4_1417=INTBIG,ipS4_1407=INTBIG,
103        ipO4_1400=INTBIG,ipO4_1397=INTBIG,ipO4_1407=INTBIG,ipO4_1405=INTBIG,ipO4_1401=INTBIG,
104        ipN3_1749=INTBIG,ipN3_1747=INTBIG,ipN3_1754=INTBIG,ipN3_1752=INTBIG,ipN3_1751=INTBIG,
105        ipC2_2325=INTBIG,ipC2_2324=INTBIG,ipC2_2329=INTBIG,ipC2_2328=INTBIG,ipC2_2327=INTBIG,
106        ipSi2_2334=INTBIG,ipSi2_2329=INTBIG,ipSi2_2350=INTBIG,ipSi2_2344=INTBIG,ipSi2_2336=INTBIG,
107        ipFe22_247=INTBIG,ipFe22_217=INTBIG,ipFe22_348=INTBIG,ipFe22_292=INTBIG,
108        ipFe22_253=INTBIG,ipFe22_846=INTBIG,ipTFe20_721=INTBIG,ipTFe20_578=INTBIG , ipZn04363=INTBIG,
109        ipS12_520=INTBIG,ipS1_25m=INTBIG,ipS1_56m=INTBIG,ipCl1_11m=INTBIG,
110        ipFe1_24m=INTBIG,ipFe1_35m=INTBIG,ipFe1_54m=INTBIG,ipFe1_111m=INTBIG,
111        ipNi1_7m=INTBIG ,ipNi1_11m=INTBIG,ipSi1_130m=INTBIG,ipSi1_68m=INTBIG;
112
113/* above are the individual pointers to the level 1 lines, they are set to
114 * very large negative. 
115 * NB NB NB!!
116 * these occur two times in the code!!
117 * They are declared in taulines.h, and defined here  */
118
119/* definition for whether level 2 lines are enabled, will be set to -1
120 * with no level2 command */
121/*long nWindLine = NWINDDIM;*/
122/*realnum TauLine2[NWINDDIM][NTA];*/
123/*realnum **TauLine2;*/
124#include "atomfeii.h"
125
126// NB NB - IS_TOP should always be the last entry of this enum!
127typedef enum { IS_NONE, IS_K_SHELL, IS_L1_SHELL, IS_L2_SHELL, IS_TOP } exc_type;
128
129struct t_BadnellLevel
130{
131        string config;
132        int irsl;
133        int S;
134        int L;
135        int g; // 2*J+1
136        realnum energy; // in cm^-1
137        bool lgAutoIonizing;
138        exc_type WhichShell;
139};
140
141// read autoionization data from Badnell data file
142STATIC void ReadBadnellAIData(const string& fnam,      // filename containing the Badnell data
143                              long nelem,              // nelem is on C scale
144                              long ion,                // ion is on C scale
145                              vector<transition>& UTA, // UTA lines will be pushed on this stack
146                              bitset<IS_TOP> Skip);    // option to skip transitions from a particular shell
147// simple helper functions for ReadBadnellAIData
148inline transition& InitTransition(transition& t);
149inline int irsl2ind(vector<t_BadnellLevel>& level, int irsl);
150
151// UTA lines below this gf value will be ignored
152// const realnum gf_cutoff = 5.5e-4f;
153const realnum gf_cutoff = 1.e-5f;
154
155void atmdat_readin(void)
156{
157        long int i,
158          iCase ,
159          ion,
160          ipDens ,
161          ipISO ,
162          ipTemp ,
163          j,
164          J,
165          ig0,
166          ig1,
167          imax,
168          nelem,
169          nelec,
170          magic1,
171          magic2,
172          mol;
173
174        char cha;
175        char chS2[3];
176
177        long ipZ;
178        bool lgEOL;
179
180        FILE *ioDATA;
181        char chLine[FILENAME_PATH_LENGTH_2] ,
182                chFilename[FILENAME_PATH_LENGTH_2];
183
184        static bool lgFirstCall = true;
185
186        DEBUG_ENTRY( "atmdat_readin()" );
187
188        /* do nothing if not first call */
189        if( !lgFirstCall )
190        {
191                /* do not do anything, but make sure that number of zones has not increased */
192                bool lgTooBig = false;
193                for( j=0; j < iterations.iter_malloc; j++ )
194                {
195                        if( geometry.nend[j]>=struc.nzlim )
196                                lgTooBig = true;
197                }
198                if( lgTooBig )
199                {
200                        fprintf(ioQQQ," This is the second or later calculation in a grid.\n");
201                        fprintf(ioQQQ," The number of zones has been increased beyond what it was on the first calculation.\n");
202                        fprintf(ioQQQ," This can\'t be done since space has already been allocated.\n");
203                        fprintf(ioQQQ," Have the first calculation do the largest number of zones so that an increase is not needed.\n");
204                        fprintf(ioQQQ," Sorry.\n");
205                        cdEXIT(EXIT_FAILURE);
206                }
207                return;
208        }
209
210        lgFirstCall = false; /* do not reevaluate again */
211
212        /* make sure that molecules have been initialized - this will fail
213         * if this routine is called before size of molecular network is known */
214        if( !mole.num_total )
215        {
216                /* mole.num_comole_calc can't be zero */
217                TotalInsanity();
218        }
219
220        /* create space for the structure variables */
221        /* nzlim will be limit, and number allocated */
222        /* >>chng 01 jul 28, define this var, do all following mallocs */
223        for( j=0; j < iterations.iter_malloc; j++ )
224        {
225                struc.nzlim = MAX2( struc.nzlim , geometry.nend[j] );
226        }
227
228        /* sloppy, but add one extra for safely */
229        ++struc.nzlim;
230
231        struc.coolstr = ((double*)MALLOC( (size_t)(struc.nzlim)*sizeof(double )));
232
233        struc.heatstr = ((double*)MALLOC( (size_t)(struc.nzlim)*sizeof(double )));
234
235        struc.testr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
236
237        struc.volstr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
238
239        struc.drad_x_fillfac = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
240
241        struc.histr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
242
243        struc.hiistr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
244
245        struc.ednstr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
246
247        struc.o3str = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
248
249        struc.pressure = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
250
251        struc.windv = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
252
253        struc.AccelTotalOutward = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
254
255        struc.AccelGravity = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
256
257        struc.pres_radiation_lines_curr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
258
259        struc.GasPressure = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
260
261        struc.hden = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
262
263        struc.DenParticles = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
264
265        struc.DenMass = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
266
267        struc.drad = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
268
269        struc.depth = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
270
271        struc.depth_last = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
272
273        struc.drad_last = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
274
275        struc.xLyman_depth = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
276
277        if (mole.num_calc != 0)
278        {
279                struc.molecules = ((realnum**)MALLOC( (size_t)(mole.num_calc)*sizeof(realnum* )));
280        }
281        else
282        {
283                struc.molecules = NULL;
284        }
285
286        for(mol=0;mol<mole.num_calc;mol++)
287        {
288                struc.molecules[mol] = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
289        }
290
291        struc.H2_abund = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
292
293        /* create space for gas phase abundances array, first create space for the elements */
294        struc.gas_phase = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)LIMELM );
295        /* last the zones */
296        for( ipZ=0; ipZ< LIMELM;++ipZ )
297        {
298                struc.gas_phase[ipZ] =
299                        (realnum*)MALLOC(sizeof(realnum)*(unsigned)(struc.nzlim) );
300        }
301
302        /* create space for struc.xIonDense array, first create space for the zones */
303        struc.xIonDense = (realnum ***)MALLOC(sizeof(realnum **)*(unsigned)(LIMELM+3) );
304
305        for( ipZ=0; ipZ< (LIMELM+3);++ipZ )
306        {
307                /* space for the ionization stage */
308                struc.xIonDense[ipZ] =
309                        (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(LIMELM+1) );
310                /* now create diagonal of space for zones */
311                for( ion=0; ion < (LIMELM+1); ++ion )
312                {
313                        struc.xIonDense[ipZ][ion] =
314                                (realnum*)MALLOC(sizeof(realnum)*(unsigned)(struc.nzlim) );
315                }
316        }
317
318        /* some structure variables */
319        for( i=0; i < struc.nzlim; i++ )
320        {
321                struc.testr[i] = 0.;
322                struc.volstr[i] = 0.;
323                struc.drad_x_fillfac[i] = 0.;
324                struc.histr[i] = 0.;
325                struc.hiistr[i] = 0.;
326                struc.ednstr[i] = 0.;
327                struc.o3str[i] = 0.;
328                struc.heatstr[i] = 0.;
329                struc.coolstr[i] = 0.;
330                struc.pressure[i] = 0.;
331                struc.pres_radiation_lines_curr[i] = 0.;
332                struc.GasPressure[i] = 0.;
333                struc.DenParticles[i] = 0.;
334                struc.depth[i] = 0.;
335                for(mol=0;mol<mole.num_calc;mol++)
336                {
337                        struc.molecules[mol][i] = 0.;
338                }
339                struc.H2_abund[i] = 0.;
340        }
341
342        /* allocate space for some arrays used by dynamics routines, and zero out vars */
343        DynaCreateArrays( );
344
345        /*************************************************************
346         *                                                           *
347         * get the level 1 lines, with real collision data set       *
348         *                                                           *
349         *************************************************************/
350
351        if( trace.lgTrace )
352                fprintf( ioQQQ," atmdat_readin reading level1.dat\n");
353
354        ioDATA = open_data( "level1.dat", "r" );
355
356        /* first line is a version number and does not count */
357        if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
358        {
359                fprintf( ioQQQ, " atmdat_readin could not read first line of level1.dat.\n");
360                cdEXIT(EXIT_FAILURE);
361        }
362
363        /* count how many lines are in the file, ignoring all lines
364         * starting with '#' */
365        nLevel1 = 0;
366        while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
367        {
368                /* we want to count the lines that do not start with #
369                 * since these contain data */
370                if( chLine[0] != '#')
371                        ++nLevel1;
372        }
373
374        /* now malloc the TauLines array */
375        TauLines = (transition *)MALLOC( (size_t)(nLevel1+1)*sizeof(transition ) );
376
377        for( i=0; i<nLevel1+1; i++ )
378        {
379                TransitionJunk( &TauLines[i] );
380                TauLines[i].Hi = AddState2Stack();
381                TauLines[i].Lo = AddState2Stack();
382                TauLines[i].Emis = AddLine2Stack( true );
383        }
384
385        /* now rewind the file so we can read it a second time*/
386        if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
387        {
388                fprintf( ioQQQ, " atmdat_readin could not rewind level1.dat.\n");
389                cdEXIT(EXIT_FAILURE);
390        }
391
392        /* check that magic number is ok */
393        if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
394        {
395                fprintf( ioQQQ, " atmdat_readin could not read first line of level1.dat.\n");
396                cdEXIT(EXIT_FAILURE);
397        }
398        i = 1;
399        /* level 1 magic number */
400        nelem = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
401        nelec = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
402        ion = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
403
404        /* magic number
405         * the following is the set of numbers that appear at the start of level1.dat  */
406        if( ( nelem != 8 ) || ( nelec != 12 ) || ( ion != 20 ) )
407        {
408                fprintf( ioQQQ,
409                        " atmdat_readin: the version of level1.dat is not the current version.\n" );
410                fprintf( ioQQQ,
411                        " Please obtain the current version from the Cloudy web site.\n" );
412                fprintf( ioQQQ,
413                        " I expected to find the number 04 11 5 and got %li %li %li instead.\n" ,
414                        nelem , nelec , ion );
415                fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
416                cdEXIT(EXIT_FAILURE);
417        }
418
419        /* this starts at 1 not 0 since zero is reserved for the
420         * dummy line - we counted number of lines above */
421        i = 1;
422        while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
423        {
424                if( i >= (nLevel1+1) )
425                {
426                        fprintf(ioQQQ," version conflict.\n");
427                }
428                /* only look at lines without '#' in first col */
429                if( chLine[0] != '#')
430                {
431                        /* >>chng 00 apr 24, prevents crash in PresTotCurrent on call from ContSetIntensity, by PvH */
432
433
434                        /*sscanf( chLine ,
435                        "%2s%2li%9li %11f %2li %2li %10g %10g %2li\n",
436                        chS2,&i2,&i3,&f1,&i4,&i5,&f2,&f3,&i6 );*/
437
438                        /* first two cols of line has chem element symbol,
439                         * try to match it with the list of names
440                         * there are no H lines in the level 1 set so zero
441                         * is an invalid result */
442
443                        /* copy first two char into chS2 and null terminate */
444                        strncpy( chS2 , chLine , 2);
445                        chS2[2] = 0;
446
447                        ipZ = 0;
448                        for( j=0; j<LIMELM; ++j)
449                        {
450                                if( strcmp( elementnames.chElementSym[j], chS2 ) == 0 )
451                                {
452                                        /* ipZ is the actual atomic number starting from 1 */
453                                        ipZ = j + 1;
454                                        break;
455                                }
456                        }
457
458                        /* this happens if no valid chemical element in first two cols on this line */
459                        if( ipZ == 0 )
460                        {
461                                fprintf( ioQQQ,
462                                        " atmdat_readin could not identify chem symbol on this level 1line:\n");
463                                fprintf( ioQQQ, "%s\n", chLine );
464                                fprintf( ioQQQ, "looking for this string==%2s==\n",chS2 );
465                                cdEXIT(EXIT_FAILURE);
466                        }
467
468                        /* now stuff them into the array, will convert over to proper units later */
469                        TauLines[i].Hi->nelem = (int)ipZ;
470                        /* start the scan early on the line -- the element label will not be
471                         * picked up, first number will be the ion stage after the label, as in c  4 */
472                        j = 1;
473                        TauLines[i].Hi->IonStg = (int)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
474                        if( lgEOL )
475                        {
476                                fprintf( ioQQQ, " There should have been a number on this level1 line 1.   Sorry.\n" );
477                                fprintf( ioQQQ, "string==%s==\n" ,chLine );
478                                cdEXIT(EXIT_FAILURE);
479                        }
480
481                        TauLines[i].Lo->nelem = TauLines[i].Hi->nelem;
482                        TauLines[i].Lo->IonStg = TauLines[i].Hi->IonStg;
483
484                        TauLines[i].WLAng = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
485                        if( lgEOL )
486                        {
487                                fprintf( ioQQQ, " There should have been a number on this level1 line 2.   Sorry.\n" );
488                                fprintf( ioQQQ, "string==%s==\n" ,chLine );
489                                cdEXIT(EXIT_FAILURE);
490                        }
491                        /*TauLines[i].WLAng = (realnum)i3;*/
492                        TauLines[i].EnergyWN = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
493                        if( lgEOL )
494                        {
495                                fprintf( ioQQQ, " There should have been a number on this level1 line 3.   Sorry.\n" );
496                                fprintf( ioQQQ, "string==%s==\n" ,chLine );
497                                cdEXIT(EXIT_FAILURE);
498                        }
499                        /*TauLines[i].EnergyWN = (realnum)f1;*/
500                        TauLines[i].Lo->g = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
501                        if( lgEOL )
502                        {
503                                fprintf( ioQQQ, " There should have been a number on this level1 line 4.   Sorry.\n" );
504                                fprintf( ioQQQ, "string==%s==\n" ,chLine );
505                                cdEXIT(EXIT_FAILURE);
506                        }
507                        /*TauLines[i].Lo->g = (realnum)i4;*/
508                        TauLines[i].Hi->g = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
509                        if( lgEOL )
510                        {
511                                fprintf( ioQQQ, " There should have been a number on this level1 line 5.   Sorry.\n" );
512                                fprintf( ioQQQ, "string==%s==\n" ,chLine );
513                                cdEXIT(EXIT_FAILURE);
514                        }
515                        /*TauLines[i].Hi->g = (realnum)i5;*/
516
517                        /* gf is log if negative */
518                        TauLines[i].Emis->gf = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
519                        if( lgEOL )
520                        {
521                                fprintf( ioQQQ, " There should have been a number on this level1 line 6.   Sorry.\n" );
522                                fprintf( ioQQQ, "string==%s==\n" ,chLine );
523                                cdEXIT(EXIT_FAILURE);
524                        }
525                        /*TauLines[i].Emis->gf = (realnum)f2;*/
526                        if( TauLines[i].Emis->gf < 0. )
527                                TauLines[i].Emis->gf = (realnum)pow((realnum)10.f,TauLines[i].Emis->gf);
528
529                        /* Emis.Aul is log if negative */
530                        TauLines[i].Emis->Aul = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
531                        if( lgEOL )
532                        {
533                                fprintf( ioQQQ, " There should have been a number on this level1 line 7.   Sorry.\n" );
534                                fprintf( ioQQQ, "string==%s==\n" ,chLine );
535                                cdEXIT(EXIT_FAILURE);
536                        }
537                        if( TauLines[i].Emis->Aul < 0. )
538                                TauLines[i].Emis->Aul = (realnum)pow((realnum)10.f,TauLines[i].Emis->Aul);
539                        /*TauLines[i].Emis->Aul = f3;*/
540
541                        TauLines[i].Emis->iRedisFun = (int)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
542                        if( lgEOL )
543                        {
544                                fprintf( ioQQQ, " There should have been a number on this level1 line 8.   Sorry.\n" );
545                                fprintf( ioQQQ, "string==%s==\n" ,chLine );
546                                cdEXIT(EXIT_FAILURE);
547                        }
548                        /*TauLines[i].Emis->iRedisFun = i6;*/
549
550                        /* this is new in c94.01 and returns nothing (0.) for most lines */
551                        TauLines[i].Coll.col_str = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
552
553                        /* finally increment i */
554                        ++i;
555                }
556        }
557
558        /* now close the file */
559        fclose(ioDATA);
560        if( trace.lgTrace )
561                fprintf( ioQQQ, " reading level1.dat OK\n" );
562
563
564        /*************************************************************
565         *                                                           *
566         * get the level 2 line, opacity project, data set           *
567         *                                                           *
568         *************************************************************/
569
570        /* nWindLine is initialized to the dimension of the vector when it is
571         * initialized in the definition at the start of this file. 
572         * it is set to -1 with the "no level2" command, which
573         * stops us from trying to establish this vector */
574        if( nWindLine > 0 )
575        {
576                /* begin level2 level 2 wind block */
577                /* open file with level 2 line data */
578
579                /* create the TauLine2 emline array */
580                TauLine2 = ((transition *)MALLOC( (size_t)nWindLine*sizeof(transition )));
581                cs1_flag_lev2 = ((realnum *)MALLOC( (size_t)nWindLine*sizeof(realnum )));
582
583                /* first initialize entire array to dangerously large negative numbers */
584                for( i=0; i< nWindLine; ++i )
585                {
586                        /* >>chng 99 jul 16, from setting all t[] to flt_max, to call for
587                         * following, each member of structure set to own type of impossible value */
588                        TransitionJunk( &TauLine2[i] );
589
590                        TauLine2[i].Hi = AddState2Stack();
591                        TauLine2[i].Lo = AddState2Stack();
592                        TauLine2[i].Emis = AddLine2Stack( true );
593                }
594
595                if( trace.lgTrace )
596                        fprintf( ioQQQ," atmdat_readin reading level2.dat\n");
597
598                ioDATA = open_data( "level2.dat", "r" );
599
600                /* get magic number off first line */
601                if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
602                {
603                        fprintf( ioQQQ, " level2.dat error getting magic number\n" );
604                        cdEXIT(EXIT_FAILURE);
605                }
606
607                /* scan off three numbers, should be the yr mn dy of creation date */
608                i = 1;
609                nelem = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
610                nelec = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
611                ion = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
612
613                /* level2.dat magic number
614                 * the following is the set of numbers that appear at the start of level2.dat 01 05 29
615                 * >>chng 06 aug 10, chng to 06, 08, 10 */
616                if( lgEOL || ( nelem != 6 ) || ( nelec != 8 ) || ( ion != 10 ) )
617                {
618                        fprintf( ioQQQ,
619                                " atmdat_readin: the version of level2.dat is not the current version.\n" );
620                        fprintf( ioQQQ,
621                                " I expected to find the number 06 08 10 and got %li %li %li instead.\n" ,
622                                nelem , nelec , ion );
623                        fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
624                        fprintf( ioQQQ, "Please obtain the correct version.\n");
625                        cdEXIT(EXIT_FAILURE);
626                }
627
628                /* now get the actual data */
629                i = 0;
630                while( i < nWindLine )
631                {
632                        /* this must be double for sscanf to work below */
633                        double tt[7];
634                        if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
635                        {
636                                fprintf( ioQQQ, " level2.dat error getting line  %li\n", i );
637                                cdEXIT(EXIT_FAILURE);
638                        }
639                        /* option to skip any line starting with # */
640                        if( chLine[0]!='#' )
641                        {
642                                /*printf("string is %s\n",chLine );*/
643                                sscanf( chLine , "%lf %lf %lf %lf %lf %lf %lf " ,
644                                        &tt[0] ,
645                                        &tt[1] ,
646                                        &tt[2] ,
647                                        &tt[3] ,
648                                        &tt[4] ,
649                                        &tt[5] ,
650                                        &tt[6] );
651                                /* these are readjusted into their final form in the structure
652                                 * in routine lines_setup*/
653                                TauLine2[i].Hi->nelem = (int)tt[0];
654                                TauLine2[i].Hi->IonStg = (int)tt[1];
655                                TauLine2[i].Lo->g = (realnum)tt[2];
656                                TauLine2[i].Hi->g = (realnum)tt[3];
657                                TauLine2[i].Emis->gf = (realnum)tt[4];
658                                TauLine2[i].EnergyWN = (realnum)tt[5];
659                                cs1_flag_lev2[i] = (realnum)tt[6];
660                                ++i;
661                        }
662                }
663                /* get magic number off last line */
664                if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
665                {
666                        fprintf( ioQQQ, " level2.dat error getting last magic number\n" );
667                        cdEXIT(EXIT_FAILURE);
668                }
669                sscanf( chLine , "%ld" , &magic2 );
670                if( 999 != magic2 )
671                {
672                        fprintf( ioQQQ, " level2.dat ends will wrong magic number=%ld \n",
673                          magic2 );
674                        cdEXIT(EXIT_FAILURE);
675                }
676                fclose( ioDATA );
677                if( trace.lgTrace )
678                        fprintf( ioQQQ," reading level2.dat OK\n");
679
680                /* end of level2 block*/
681        }
682
683        /*************************
684         * 
685         * get the UTA line sets   
686         * >>chng 06 nov 26, use Gu et al. 2006 UTA data by default
687         * with option to use older Behar
688         * 
689         *************************/
690
691        /* reserve space for all data sets, no worries if too small though... */
692        UTALines.reserve( 4000 );
693
694        // version of element symbols in lower case and without spaces....
695        const char* chElmSymLC[] =
696                { "h", "he", "li", "be", "b", "c", "n", "o", "f", "ne", "na", "mg", "al", "si", "p",
697                  "s", "cl", "ar", "k", "ca", "sc", "ti", "v", "cr", "mn", "fe", "co", "ni", "cu", "zn" };
698
699        /* first read in the Badnell data */
700        for( ipISO=ipLI_LIKE; ipISO < ipAL_LIKE; ++ipISO )
701        {
702                for( nelem=ipISO; nelem < LIMELM; ++nelem )
703                {
704                        // ion = 0 for neutral atom
705                        long ion = nelem - ipISO;
706                        bitset<IS_TOP> Skip;
707                        if( ipISO < ipNA_LIKE )
708                        {
709                                // construct file name
710                                ostringstream oss;
711                                // Badnell calls Li-like series He-like, etc...
712                                oss << "nrb00#" << chElmSymLC[ipISO-1] << "_";
713                                // Badnell uses ion = 1 for neutral atom
714                                oss << chElmSymLC[nelem] << ion+1 << "ic1-2.dat";
715                                // now read the data...
716                                ReadBadnellAIData( oss.str(), nelem, ion, UTALines, Skip );
717                        }
718                        else
719                        {
720                                // from Na-like onwards both K-shell (ic1-3) and L-shell (ic2-3)
721                                // excitation are treated in two separate files
722                                ostringstream oss;
723                                oss << "nrb00#" << chElmSymLC[ipISO-1] << "_";
724                                oss << chElmSymLC[nelem] << ion+1 << "ic1-3.dat";
725                                ReadBadnellAIData( oss.str(), nelem, ion, UTALines, Skip );
726
727                                // Kisielius L2-shell data set takes precedence for Na, Mg, and Al-like iron
728                                if( ionbal.lgInnerShell_Kisielius && nelem == ipIRON &&
729                                    ipISO >= ipNA_LIKE && ipISO <= ipAL_LIKE )
730                                        Skip[IS_L2_SHELL] = 1;
731
732                                ostringstream oss2;
733                                oss2 << "nrb00#" << chElmSymLC[ipISO-1] << "_";
734                                oss2 << chElmSymLC[nelem] << ion+1 << "ic2-3.dat";
735                                ReadBadnellAIData( oss2.str(), nelem, ion, UTALines, Skip );
736                        }
737                }
738        }
739
740        /* these are the statistical weights for the ground levels of all ions of iron */
741        const realnum StatWeightGroundLevelIron[] =
742                { 9.f, 10.f, 9.f, 6.f, 1.f, 4.f, 5.f, 4.f, 1.f, 4.f, 5.f, 4.f, 1.f,
743                  2.f, 1.f, 2.f, 1.f, 4.f, 5.f, 4.f, 1.f, 2.f, 1.f, 2.f, 1.f, 2.f };
744
745        // blank line that will be pushed on the UTA line stack
746        transition BlankLine;
747        TransitionJunk( &BlankLine );
748
749        /* next read in the Gu file */
750        if( ionbal.lgInnerShell_Gu06 )
751        {
752                /* read the Gu et al. (2006) data
753                 * >>refer      Fe      inner shell UTA Gu, M. F.; Holczer T.; Behar E.; Kahn S. M.;
754                 * >>refercon   2006, ApJ 641, 1227-1232 */
755                if( trace.lgTrace )
756                        fprintf( ioQQQ," atmdat_readin reading UTA_Gu06.dat\n");
757
758                FILE *ioGU06 = open_data( "UTA_Gu06.dat", "r" );
759
760                nelem = 0;
761                /* get up to magic number in Gu 06 file - there is a large header of
762                 * comments with the first data the magic number */
763                while( read_whole_line( chLine , (int)sizeof(chLine) , ioGU06 ) != NULL )
764                {
765                        /* we want to break on first line that does not start with #
766                         * since that contains data */
767                        if( chLine[0] != '#')
768                                break;
769                }
770                /* now get Gu magic number */
771                sscanf( chLine, "%li %li %li", &nelem, &nelec, &ion );
772
773                /* is magic number correct?
774                 * the following is the set of numbers that appear at the start of UTA_Gu06.dat 2006 11 17 */
775                if( nelem != 2007 || nelec != 1 || ion != 23 )
776                {
777                        fprintf( ioQQQ,
778                                " atmdat_readin: the version of UTA_Gu06.dat is not the current version.\n" );
779                        fprintf( ioQQQ,
780                                " I expected to find the number 2007 1 23 and got %li %li %li instead.\n" ,
781                                nelem , nelec , ion );
782                        fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
783                        cdEXIT(EXIT_FAILURE);
784                }
785
786                while( read_whole_line( chLine, (int)sizeof(chLine), ioGU06 ) != NULL )
787                {
788                        if( chLine[0] != '#' )
789                        {
790                                long int i2;
791                                double WLAng, Aul, oscill, Aauto;
792
793                                sscanf( chLine, "%4li%5li%8lf%13lf%12lf",
794                                        &ion, &i2, &WLAng, &Aul, &Aauto );
795                                sscanf( &chLine[54], "%13lf", &oscill );
796
797                                // avoid duplication of ions: anything upto and including the Mg-like
798                                // series is covered by the Badnell data set, Al-like is covered by
799                                // the Kisielius data set if it is turned on.
800                                ipISO = ipIRON - ion + 1;
801                                int ipThres = ionbal.lgInnerShell_Kisielius ? ipAL_LIKE : ipMG_LIKE;
802                                if( ipISO <= ipThres )
803                                        continue;
804
805                                UTALines.push_back( InitTransition( BlankLine ) );
806
807                                /* all these are iron, first number was ion stage with 0 the atom */
808                                UTALines.back().Hi->nelem = ipIRON+1;
809
810                                /* now do stage of ionization */
811                                ASSERT( ion > 0 && ion <= ipIRON );
812                                /* the ion stage - 1 is atom - this data file has different
813                                 * format from other two - others are 0 for atom */
814                                UTALines.back().Hi->IonStg = ion;
815
816                                /* these are the statistical weights
817                                 * lower levels are not included in the original data file */
818                                if( strstr( chLine, "(J=1/2)" ) != NULL )
819                                        UTALines.back().Hi->g = 2.f;
820                                else if( strstr( chLine, "(J=1)" ) != NULL )
821                                        UTALines.back().Hi->g = 3.f;
822                                else if( strstr( chLine, "(J=3/2)" ) != NULL )
823                                        UTALines.back().Hi->g = 4.f;
824                                else if( strstr( chLine, "(J=2)" ) != NULL )
825                                        UTALines.back().Hi->g = 5.f;
826                                else if( strstr( chLine, "(J=5/2)" ) != NULL )
827                                        UTALines.back().Hi->g = 6.f;
828                                else if( strstr( chLine, "(J=3)" ) != NULL )
829                                        UTALines.back().Hi->g = 7.f;
830                                else if( strstr( chLine, "(J=7/2)" ) != NULL )
831                                        UTALines.back().Hi->g = 8.f;
832                                else if( strstr( chLine, "(J=4)" ) != NULL )
833                                        UTALines.back().Hi->g = 9.f;
834                                else if( strstr( chLine, "(J=9/2)" ) != NULL )
835                                        UTALines.back().Hi->g = 10.f;
836                                else if( strstr( chLine, "(J=5)" ) != NULL )
837                                        UTALines.back().Hi->g = 11.f;
838                                else if( strstr( chLine, "(J=11/2)" ) != NULL )
839                                        UTALines.back().Hi->g = 12.f;
840                                else
841                                        TotalInsanity();
842                                UTALines.back().Lo->g = StatWeightGroundLevelIron[ion-1];
843
844                                /* wavelength in Angstroms */
845                                UTALines.back().WLAng = (realnum)WLAng;
846                                UTALines.back().EnergyWN = (realnum)(1e8/WLAng);
847                                UTALines.back().EnergyErg = (realnum)(ERG1CM)*UTALines.back().EnergyWN;
848
849                                /* store branching ratio for autoionization */
850                                double frac_ioniz = Aauto/(Aul + Aauto);
851                                ASSERT( frac_ioniz >= 0. &&  frac_ioniz <= 1. );
852                                UTALines.back().Emis->AutoIonizFrac = (realnum)frac_ioniz;
853
854                                /* save gf scanned from line */
855                                UTALines.back().Emis->gf = UTALines.back().Lo->g * (realnum)oscill;
856
857                                /* this is true spontaneous rate for doubly excited state to inner
858                                 * shell UTA, and is used for pumping, and also relaxing back to inner shell */
859                                UTALines.back().Emis->Aul =
860                                        (realnum)eina( UTALines.back().Emis->gf,
861                                                       UTALines.back().EnergyWN,
862                                                       UTALines.back().Hi->g );
863
864                                ASSERT( fp_equal_tol( (realnum)Aul, UTALines.back().Emis->Aul, 1.e-3f*(realnum)Aul ) );
865
866                                UTALines.back().Emis->iRedisFun = ipPRD;