/home66/gary/public_html/cloudy/c08_branch/source/atmdat_readin.cpp

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*atmdat_readin read in some data files, but only if this is very first call, 
 * called by Cloudy */
#include "cddefines.h"
#include "physconst.h"
#include "taulines.h"
#include "mewecoef.h"
#include "iterations.h"
#include "heavy.h"
#include "mole.h"
#include "yield.h"
#include "trace.h"
#include "lines.h"
#include "lines_service.h"
#include "ionbal.h"
#include "struc.h"
#include "geometry.h"
#include "dynamics.h"
#include "elementnames.h"
#include "hyperfine.h"
#include "atmdat.h"
/* */
/* this was needed to get array to crash out of bounds if not set.
 * std limits on limits.h did not work with visual studio! */
#define INTBIG 2000000000

/* these are the individual pointers to the level 1 lines, they are set to
 * very large negative.  
 * NB NB NB!!
 * these occur two times in the code!!
 * They are declared in taulines.h, and defined here  */
long 
   ipT1656=INTBIG,    ipT9830=INTBIG,  /*ipH21cm=INTBIG, ipHe3cm=INTBIG,*/
   ipT8727=INTBIG,    ipT1335=INTBIG,  ipT1909=INTBIG,
        ipT977=INTBIG,    ipT1550=INTBIG,  ipT1548=INTBIG, ipT386=INTBIG, 
        ipT310=INTBIG,    ipc31175=INTBIG, ipT291=INTBIG,  ipT280=INTBIG,
        ipT274=INTBIG,    ipT270=INTBIG,   ipT312=INTBIG,  ipT610=INTBIG, 
        ipT370=INTBIG,    ipT157=INTBIG,   ipT1085=INTBIG, 
        ipT990=INTBIG,    ipT1486=INTBIG,  ipT765=INTBIG,  ipT1243=INTBIG, 
        ipT1239=INTBIG,   ipT374g=INTBIG,  ipT374x=INTBIG, ipT1200=INTBIG,
        ipT2140=INTBIG,   ipT671=INTBIG,   ipT315=INTBIG,  ipT324=INTBIG, 
        ipT333=INTBIG,    ipT209=INTBIG,   ipT122=INTBIG,  ipT205=INTBIG,
        ipT57=INTBIG,     ipT6300=INTBIG,  ipT6363=INTBIG, ipT5577=INTBIG, 
        ipT834=INTBIG,    ipT1661=INTBIG,  ipT1666=INTBIG, ipT835=INTBIG,
        ipT789=INTBIG,    ipT630=INTBIG,   ipT1304=INTBIG, ipSi10_606=INTBIG,
        ipT1039=INTBIG,   ipT8446=INTBIG,  ipT4368=INTBIG, ipTOI13=INTBIG,
        ipTOI11=INTBIG,   ipTOI29=INTBIG,  ipTOI46=INTBIG, ipTO1025=INTBIG, 
        ipT304=INTBIG,    ipT1214=INTBIG,  ipT150=INTBIG,  ipT146=INTBIG,
        ipT63=INTBIG,     ipTO88=INTBIG,   ipT52=INTBIG,   ipT26=INTBIG, 
        ipT1032=INTBIG,   ipT1037=INTBIG,  ipF0229=INTBIG, ipF0267=INTBIG,
        ipF444=INTBIG,    ipF425=INTBIG,   ipT770=INTBIG, ipT780=INTBIG, 
        ipxNe0676=INTBIG, ipT895=INTBIG,   ipT88=INTBIG, ipTNe13=INTBIG,
        ipTNe36=INTBIG,   ipTNe16=INTBIG,  ipTNe14=INTBIG, ipTNe24=INTBIG, 
        ipT5895=INTBIG,   ipfsNa373=INTBIG,ipfsNa490=INTBIG, ipfsNa421=INTBIG,
        ipxNa6143=INTBIG, ipxNa6862=INTBIG,ipxNa0746=INTBIG, ipMgI2853=INTBIG, 
        ipMgI2026=INTBIG, ipT2796=INTBIG,  ipT2804=INTBIG,
        ipT705=INTBIG,    ipT4561=INTBIG,  ipxMg51325=INTBIG, ipxMg52417=INTBIG, 
        ipxMg52855=INTBIG,ipxMg71190=INTBIG, ipxMg72261=INTBIG,
        ipxMg72569=INTBIG,ipxMg08303=INTBIG, ipTMg610=INTBIG, ipTMg625=INTBIG, 
        ipT58=INTBIG,     ipTMg4=INTBIG, ipTMg14=INTBIG, ipTMg6=INTBIG,
        ipfsMg790=INTBIG, ipfsMg755=INTBIG, ipAlI3957=INTBIG, ipAlI3090=INTBIG, 
        ipT1855=INTBIG,   ipT1863=INTBIG, ipT2670=INTBIG,
        ipAl529=INTBIG,   ipAl6366=INTBIG, ipAl6912=INTBIG, ipAl8575=INTBIG, 
        ipAl8370=INTBIG,  ipAl09204=INTBIG, ipT639=INTBIG,
        ipTAl550=INTBIG,  ipTAl568=INTBIG, ipTAl48=INTBIG, ipSii2518=INTBIG, 
        ipSii2215=INTBIG, ipT1808=INTBIG,
        ipT1207=INTBIG,   ipT1895=INTBIG, ipT1394=INTBIG, ipT1403=INTBIG, 
        ipT1527=INTBIG,   ipT1305=INTBIG, ipT1260=INTBIG, ipSi619=INTBIG,
        ipSi10143=INTBIG, ipTSi499=INTBIG, ipTSi521=INTBIG, ipTSi41=INTBIG, 
        ipTSi35=INTBIG,   ipTSi25=INTBIG, ipTSi65=INTBIG,
        ipTSi3=INTBIG,    ipTSi4=INTBIG, ipP0260=INTBIG, ipP0233=INTBIG, 
        ipP0318=INTBIG,   ipP713=INTBIG, ipP848=INTBIG, ipP817=INTBIG,
        ipP1027=INTBIG,   ipP1018=INTBIG, ipT1256=INTBIG, ipT1194=INTBIG, 
        ipTS1720=INTBIG,  ipT1198=INTBIG, ipT786=INTBIG,
        ipT933=INTBIG,    ipT944=INTBIG, ipfsS810=INTBIG, ipfsS912=INTBIG, 
        ipfsS938=INTBIG,  ipfsS1119=INTBIG, ipfsS1114=INTBIG, ipfsS1207=INTBIG,
        ipTSu418=INTBIG,  ipTSu446=INTBIG, ipTSu30=INTBIG, ipTS19=INTBIG, 
        ipTS34=INTBIG,    ipTS11=INTBIG, ipfsCl214=INTBIG, ipfsCl233=INTBIG,
        ipCl04203=INTBIG, ipCl04117=INTBIG, ipCl973=INTBIG, ipCl1030=INTBIG, 
        ipCl1092=INTBIG,  ipT354=INTBIG, ipT389=INTBIG, ipT25=INTBIG,
        ipTAr7=INTBIG,    ipTAr9=INTBIG, ipTAr22=INTBIG, ipTAr13=INTBIG, 
        ipTAr8=INTBIG,    ipAr06453=INTBIG, ipAr1055=INTBIG, ipAr1126=INTBIG,
        ipAr1178=INTBIG,  ipKI7745=INTBIG, ipxK03462=INTBIG, ipxK04598=INTBIG, 
        ipxK04154=INTBIG, ipxK06882=INTBIG, ipxK06557=INTBIG,
        ipxK07319=INTBIG, ipxK11425=INTBIG, ipCaI4228=INTBIG, ipT3934=INTBIG, 
        ipT3969=INTBIG,   ipT8498=INTBIG, ipT8542=INTBIG,
        ipT8662=INTBIG,   ipT7291=INTBIG, ipT7324=INTBIG, ipTCa302=INTBIG, 
        ipTCa345=INTBIG,  ipTCa19=INTBIG, ipTCa3=INTBIG, ipTCa12=INTBIG,
        ipTCa4=INTBIG,    ipCa0741=INTBIG, ipCa0761=INTBIG, ipCa08232=INTBIG, 
        ipCa12333=INTBIG, ipSc05231=INTBIG, ipSc13264=INTBIG,
        ipTi06172=INTBIG, ipTi14212=INTBIG, ipVa07130=INTBIG, ipVa15172=INTBIG, 
        ipCr08101=INTBIG, ipCr16141=INTBIG, ipxMn0979=INTBIG,
        ipxMn1712=INTBIG, ipFeI3884=INTBIG, ipFeI3729=INTBIG, ipFeI3457=INTBIG, 
        ipFeI3021=INTBIG, ipFeI2966=INTBIG, ipTuv3=INTBIG,
        ipTr48=INTBIG,    ipTFe16=INTBIG, ipTFe26=INTBIG, ipTFe34=INTBIG, 
        ipTFe35=INTBIG,   ipTFe46=INTBIG, ipTFe56=INTBIG, ipT1122=INTBIG,
        ipFe0795=INTBIG,  ipFe0778=INTBIG, ipT245=INTBIG, ipT352=INTBIG, 
        ipFe106375=INTBIG,ipT353=INTBIG, /*ipFe1310=INTBIG, ipFe1311=INTBIG,*/
        ipT347=INTBIG,    ipT192=INTBIG, ipT255=INTBIG, ipT11=INTBIG, 
        ipT191=INTBIG,    /*ipTFe07=INTBIG, ipTFe61=INTBIG,*/ ipFe18975=INTBIG, ipTFe23=INTBIG,
        ipTFe13=INTBIG,   ipCo11527=INTBIG, ipxNi1242=INTBIG;
long    ipS4_1405=INTBIG,ipS4_1398=INTBIG,ipS4_1424=INTBIG,ipS4_1417=INTBIG,ipS4_1407=INTBIG,
                ipO4_1400=INTBIG,ipO4_1397=INTBIG,ipO4_1407=INTBIG,ipO4_1405=INTBIG,ipO4_1401=INTBIG,
                ipN3_1749=INTBIG,ipN3_1747=INTBIG,ipN3_1754=INTBIG,ipN3_1752=INTBIG,ipN3_1751=INTBIG,
                ipC2_2325=INTBIG,ipC2_2324=INTBIG,ipC2_2329=INTBIG,ipC2_2328=INTBIG,ipC2_2327=INTBIG,
                ipSi2_2334=INTBIG,ipSi2_2329=INTBIG,ipSi2_2350=INTBIG,ipSi2_2344=INTBIG,ipSi2_2336=INTBIG,
                ipFe22_247=INTBIG,ipFe22_217=INTBIG,ipFe22_348=INTBIG,ipFe22_292=INTBIG,
                ipFe22_253=INTBIG,ipFe22_846=INTBIG,ipTFe20_721=INTBIG,ipTFe20_578=INTBIG , ipZn04363=INTBIG,
                ipS12_520=INTBIG,ipS1_25m=INTBIG,ipS1_56m=INTBIG,ipCl1_11m=INTBIG,
                ipFe1_24m=INTBIG,ipFe1_35m=INTBIG,ipFe1_54m=INTBIG,ipFe1_111m=INTBIG,
                ipNi1_7m=INTBIG ,ipNi1_11m=INTBIG,ipSi1_130m=INTBIG,ipSi1_68m=INTBIG;

/* above are the individual pointers to the level 1 lines, they are set to
 * very large negative.  
 * NB NB NB!!
 * these occur two times in the code!!
 * They are declared in taulines.h, and defined here  */

/* definition for whether level 2 lines are enabled, will be set to -1 
 * with no level2 command */
/*long nWindLine = NWINDDIM;*/
/*realnum TauLine2[NWINDDIM][NTA];*/
/*realnum **TauLine2;*/
#include "atomfeii.h"

void atmdat_readin(void)
{
        long int i, 
          iCase ,
          ion, 
          ipDens ,
          ipISO ,
          ipTemp ,
          j, 
          J,
          ig0, 
          ig1, 
          imax, 
          nelem, 
          nelec, 
          magic1, 
          magic2,
          mol ,

          nUTA_Romas,
          nFeIonRomas[ipIRON],

          nUTA_Behar ,
          nFeIonBehar[ipIRON],

          nUTA_Gu06,
          nFeIonGu[ipIRON];
        double
          WLloRomas[ipIRON],
          WLhiRomas[ipIRON],
          WLloBehar[ipIRON],
          WLhiBehar[ipIRON],
          WLloGu[ipIRON],
          WLhiGu[ipIRON];

        char cha;
        char chS2[3];

        long ipZ;
        bool lgEOL;

        FILE *ioDATA, *ioROMAS , *ioBEHAR=NULL;
        FILE *ioGU06=NULL;
        char chLine[FILENAME_PATH_LENGTH_2] , 
                chFilename[FILENAME_PATH_LENGTH_2];

        static bool lgFirstCall = true;

        DEBUG_ENTRY( "atmdat_readin()" );

        /* do nothing if not first call */
        if( !lgFirstCall )
        { 
                /* do not do anything, but make sure that number of zones has not increased */
                bool lgTooBig = false;
                for( j=0; j < iterations.iter_malloc; j++ )
                {
                        if( geometry.nend[j]>=struc.nzlim )
                                lgTooBig = true;
                }
                if( lgTooBig )
                {
                        fprintf(ioQQQ," This is the second or later calculation in a grid.\n");
                        fprintf(ioQQQ," The number of zones has been increased beyond what it was on the first calculation.\n");
                        fprintf(ioQQQ," This can\'t be done since space has already been allocated.\n");
                        fprintf(ioQQQ," Have the first calculation do the largest number of zones so that an increase is not needed.\n");
                        fprintf(ioQQQ," Sorry.\n");
                        cdEXIT(EXIT_FAILURE);
                }
                return;
        }

        lgFirstCall = false; /* do not reevaluate again */

        /* make sure that molecules have been initialized - this will fail
         * if this routine is called before size of molecular network is known */
        if( !mole.num_comole_calc )
        {
                /* mole.num_comole_calc can't be zero */
                TotalInsanity();
        }

        /* create space for the structure variables */
        /* nzlim will be limit, and number allocated */
        /* >>chng 01 jul 28, define this var, do all following mallocs */
        for( j=0; j < iterations.iter_malloc; j++ )
        {
                struc.nzlim = MAX2( struc.nzlim , geometry.nend[j] );
        }

        /* sloppy, but add one extra for safely */
        ++struc.nzlim;

        struc.coolstr = ((double*)MALLOC( (size_t)(struc.nzlim)*sizeof(double )));

        struc.heatstr = ((double*)MALLOC( (size_t)(struc.nzlim)*sizeof(double )));

        struc.testr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));

        struc.volstr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));

        struc.drad_x_fillfac = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));

        struc.histr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));

        struc.hiistr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));

        struc.ednstr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));

        struc.o3str = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));

        struc.pressure = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));

        struc.windv = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));

        struc.AccelTot = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));

        struc.AccelGravity = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));

        struc.pres_radiation_lines_curr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));

        struc.GasPressure = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));

        struc.hden = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));

        struc.DenParticles = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));

        struc.DenMass = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));

        struc.drad = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));

        struc.depth = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));

        struc.depth_last = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));

        struc.drad_last = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));

        struc.xLyman_depth = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));

        struc.H2_molec = ((realnum**)MALLOC( (size_t)(N_H_MOLEC)*sizeof(realnum* )));

        struc.CO_molec = ((realnum**)MALLOC( (size_t)(mole.num_comole_calc)*sizeof(realnum* )));

        for(mol=0;mol<N_H_MOLEC;mol++)
        {
                struc.H2_molec[mol] = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
        }

        for(mol=0;mol<mole.num_comole_calc;mol++)
        {
                struc.CO_molec[mol] = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
        }

        /* create space for gas phase abundances array, first create space for the elements */
        struc.gas_phase = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)LIMELM );
        /* last the zones */
        for( ipZ=0; ipZ< LIMELM;++ipZ )
        {
                struc.gas_phase[ipZ] = 
                        (realnum*)MALLOC(sizeof(realnum)*(unsigned)(struc.nzlim) );
        }

        /* create space for struc.xIonDense array, first create space for the zones */
        struc.xIonDense = (realnum ***)MALLOC(sizeof(realnum **)*(unsigned)(LIMELM+3) );

        for( ipZ=0; ipZ< (LIMELM+3);++ipZ )
        {
                /* space for the ionization stage */
                struc.xIonDense[ipZ] = 
                        (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(LIMELM+1) );
                /* now create diagonal of space for zones */
                for( ion=0; ion < (LIMELM+1); ++ion )
                {
                        struc.xIonDense[ipZ][ion] = 
                                (realnum*)MALLOC(sizeof(realnum)*(unsigned)(struc.nzlim) );
                }
        }

        /* some structure variables */
        for( i=0; i < struc.nzlim; i++ )
        {
                struc.testr[i] = 0.;
                struc.volstr[i] = 0.;
                struc.drad_x_fillfac[i] = 0.;
                struc.histr[i] = 0.;
                struc.hiistr[i] = 0.;
                struc.ednstr[i] = 0.;
                struc.o3str[i] = 0.;
                struc.heatstr[i] = 0.;
                struc.coolstr[i] = 0.;
                struc.pressure[i] = 0.;
                struc.pres_radiation_lines_curr[i] = 0.;
                struc.GasPressure[i] = 0.;
                struc.DenParticles[i] = 0.;
                struc.depth[i] = 0.;
                for(mol=0;mol<N_H_MOLEC;mol++)
                {
                        struc.H2_molec[mol][i] = 0.;
                }
                for(mol=0;mol<mole.num_comole_calc;mol++)
                {
                        struc.CO_molec[mol][i] = 0.;
                }
        }

        /* allocate space for some arrays used by dynamics routines, and zero out vars */
        DynaCreateArrays( );

        /*************************************************************
         *                                                           *
         * get the level 1 lines, with real collision data set       *
         *                                                           *
         *************************************************************/

        if( trace.lgTrace )
                fprintf( ioQQQ," atmdat_readin opening level1.dat:");

        ioDATA = open_data( "level1.dat", "r" );

        /* first line is a version number and does not count */
        if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
        {
                fprintf( ioQQQ, " atmdat_readin could not read first line of level1.dat.\n");
                cdEXIT(EXIT_FAILURE);
        }

        /* count how many lines are in the file, ignoring all lines
         * starting with '#' */
        nLevel1 = 0;
        while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
        {
                /* we want to count the lines that do not start with #
                 * since these contain data */
                if( chLine[0] != '#')
                        ++nLevel1;
        }

        /* now malloc the TauLines array */
        TauLines = (transition *)MALLOC( (size_t)(nLevel1+1)*sizeof(transition ) );

        for( i=0; i<nLevel1+1; i++ )
        {
                TransitionJunk( &TauLines[i] );
                TauLines[i].Hi = AddState2Stack();
                TauLines[i].Lo = AddState2Stack();
                TauLines[i].Emis = AddLine2Stack( true );
        }

        /* now rewind the file so we can read it a second time*/
        if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
        {
                fprintf( ioQQQ, " atmdat_readin could not rewind level1.dat.\n");
                cdEXIT(EXIT_FAILURE);
        }

        /* check that magic number is ok */
        if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
        {
                fprintf( ioQQQ, " atmdat_readin could not read first line of level1.dat.\n");
                cdEXIT(EXIT_FAILURE);
        }
        i = 1;
        /* level 1 magic number */
        nelem = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
        nelec = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
        ion = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);

        /* magic number
         * the following is the set of numbers that appear at the start of level1.dat  */
        if( ( nelem != 8 ) || ( nelec != 12 ) || ( ion != 20 ) )
        {
                fprintf( ioQQQ, 
                        " atmdat_readin: the version of level1.dat is not the current version.\n" );
                fprintf( ioQQQ, 
                        " Please obtain the current version from the Cloudy web site.\n" );
                fprintf( ioQQQ, 
                        " I expected to find the number 04 11 5 and got %li %li %li instead.\n" ,
                        nelem , nelec , ion );
                fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
                cdEXIT(EXIT_FAILURE);
        }

        /* this starts at 1 not 0 since zero is reserved for the
         * dummy line - we counted number of lines above */
        i = 1;
        while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
        {
                if( i >= (nLevel1+1) )
                {
                        fprintf(ioQQQ," version conflict.\n");
                }
                /* only look at lines without '#' in first col */
                if( chLine[0] != '#')
                {
                        /* >>chng 00 apr 24, prevents crash in PresTotCurrent on call from ContSetIntensity, by PvH */


                        /*sscanf( chLine ,
                        "%2s%2li%9li %11f %2li %2li %10g %10g %2li\n",
                        chS2,&i2,&i3,&f1,&i4,&i5,&f2,&f3,&i6 );*/

                        /* first two cols of line has chem element symbol, 
                         * try to match it with the list of names 
                         * there are no H lines in the level 1 set so zero
                         * is an invalid result */

                        /* copy first two char into chS2 and null terminate */
                        strncpy( chS2 , chLine , 2);
                        chS2[2] = 0;

                        ipZ = 0;
                        for( j=0; j<LIMELM; ++j)
                        {
                                if( strcmp( elementnames.chElementSym[j], chS2 ) == 0 )
                                {
                                        /* ipZ is the actual atomic number starting from 1 */
                                        ipZ = j + 1;
                                        break;
                                }
                        }

                        /* this happens if no valid chemical element in first two cols on this line */
                        if( ipZ == 0 )
                        {
                                fprintf( ioQQQ, 
                                        " atmdat_readin could not identify chem symbol on this level 1line:\n");
                                fprintf( ioQQQ, "%s\n", chLine );
                                fprintf( ioQQQ, "looking for this string==%2s==\n",chS2 );
                                cdEXIT(EXIT_FAILURE);
                        }

                        /* now stuff them into the array, will convert over to proper units later */
                        TauLines[i].Hi->nelem = (int)ipZ;
                        /* start the scan early on the line -- the element label will not be
                         * picked up, first number will be the ion stage after the label, as in c  4 */
                        j = 1;
                        TauLines[i].Hi->IonStg = (int)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                fprintf( ioQQQ, " There should have been a number on this level1 line 1.   Sorry.\n" );
                                fprintf( ioQQQ, "string==%s==\n" ,chLine );
                                cdEXIT(EXIT_FAILURE);
                        }

                        TauLines[i].Lo->nelem = TauLines[i].Hi->nelem;
                        TauLines[i].Lo->IonStg = TauLines[i].Hi->IonStg;

                        TauLines[i].WLAng = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                fprintf( ioQQQ, " There should have been a number on this level1 line 2.   Sorry.\n" );
                                fprintf( ioQQQ, "string==%s==\n" ,chLine );
                                cdEXIT(EXIT_FAILURE);
                        }
                        /*TauLines[i].WLAng = (realnum)i3;*/
                        TauLines[i].EnergyWN = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                fprintf( ioQQQ, " There should have been a number on this level1 line 3.   Sorry.\n" );
                                fprintf( ioQQQ, "string==%s==\n" ,chLine );
                                cdEXIT(EXIT_FAILURE);
                        }
                        /*TauLines[i].EnergyWN = (realnum)f1;*/
                        TauLines[i].Lo->g = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                fprintf( ioQQQ, " There should have been a number on this level1 line 4.   Sorry.\n" );
                                fprintf( ioQQQ, "string==%s==\n" ,chLine );
                                cdEXIT(EXIT_FAILURE);
                        }
                        /*TauLines[i].Lo->g = (realnum)i4;*/
                        TauLines[i].Hi->g = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                fprintf( ioQQQ, " There should have been a number on this level1 line 5.   Sorry.\n" );
                                fprintf( ioQQQ, "string==%s==\n" ,chLine );
                                cdEXIT(EXIT_FAILURE);
                        }
                        /*TauLines[i].Hi->g = (realnum)i5;*/

                        /* gf is log if negative */
                        TauLines[i].Emis->gf = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                fprintf( ioQQQ, " There should have been a number on this level1 line 6.   Sorry.\n" );
                                fprintf( ioQQQ, "string==%s==\n" ,chLine );
                                cdEXIT(EXIT_FAILURE);
                        }
                        /*TauLines[i].Emis->gf = (realnum)f2;*/
                        if( TauLines[i].Emis->gf < 0. )
                                TauLines[i].Emis->gf = (realnum)pow((realnum)10.f,TauLines[i].Emis->gf);

                        /* Emis.Aul is log if negative */
                        TauLines[i].Emis->Aul = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                fprintf( ioQQQ, " There should have been a number on this level1 line 7.   Sorry.\n" );
                                fprintf( ioQQQ, "string==%s==\n" ,chLine );
                                cdEXIT(EXIT_FAILURE);
                        }
                        if( TauLines[i].Emis->Aul < 0. )
                                TauLines[i].Emis->Aul = (realnum)pow((realnum)10.f,TauLines[i].Emis->Aul);
                        /*TauLines[i].Emis->Aul = f3;*/

                        TauLines[i].Emis->iRedisFun = (int)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                fprintf( ioQQQ, " There should have been a number on this level1 line 8.   Sorry.\n" );
                                fprintf( ioQQQ, "string==%s==\n" ,chLine );
                                cdEXIT(EXIT_FAILURE);
                        }
                        /*TauLines[i].Emis->iRedisFun = i6;*/

                        /* this is new in c94.01 and returns nothing (0.) for most lines */
                        TauLines[i].Coll.cs = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);

                        /* finally increment i */
                        ++i;
                }
        }

        /* now close the file */
        fclose(ioDATA);

        /*************************************************************
         *                                                           *
         * get the level 2 line, opacity project, data set           *
         *                                                           *
         *************************************************************/

        /* nWindLine is initialized to the dimension of the vector when it is
         * initialized in the definition at the start of this file.  
         * it is set to -1 with the "no level2" command, which
         * stops us from trying to establish this vector */
        if( nWindLine > 0 )
        {
                /* begin level2 level 2 wind block */
                /* open file with level 2 line data */

                /* create the TauLine2 emline array */
                TauLine2 = ((transition *)MALLOC( (size_t)nWindLine*sizeof(transition )));
                cs1_flag_lev2 = ((realnum *)MALLOC( (size_t)nWindLine*sizeof(realnum )));

                /* first initialize entire array to dangerously large negative numbers */
                for( i=0; i< nWindLine; ++i )
                {
                        /* >>chng 99 jul 16, from setting all t[] to flt_max, to call for
                         * following, each member of structure set to own type of impossible value */
                        TransitionJunk( &TauLine2[i] );

                        TauLine2[i].Hi = AddState2Stack();
                        TauLine2[i].Lo = AddState2Stack();
                        TauLine2[i].Emis = AddLine2Stack( true );
                }

                if( trace.lgTrace )
                        fprintf( ioQQQ," atmdat_readin opening level2.dat:");

                ioDATA = open_data( "level2.dat", "r" );

                /* get magic number off first line */
                if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
                {
                        fprintf( ioQQQ, " level2.dat error getting magic number\n" );
                        cdEXIT(EXIT_FAILURE);
                }

                /* scan off three numbers, should be the yr mn dy of creation date */
                i = 1;
                nelem = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
                nelec = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
                ion = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);

                /* level2.dat magic number
                 * the following is the set of numbers that appear at the start of level2.dat 01 05 29 
                 * >>chng 06 aug 10, chng to 06, 08, 10 */
                if( lgEOL || ( nelem != 6 ) || ( nelec != 8 ) || ( ion != 10 ) )
                {
                        fprintf( ioQQQ, 
                                " atmdat_readin: the version of level2.dat is not the current version.\n" );
                        fprintf( ioQQQ, 
                                " I expected to find the number 06 08 10 and got %li %li %li instead.\n" ,
                                nelem , nelec , ion );
                        fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
                        fprintf( ioQQQ, "Please obtain the correct version.\n");
                        cdEXIT(EXIT_FAILURE);
                }

                /* now get the actual data */
                i = 0;
                while( i < nWindLine )
                {
                        /* this must be double for sscanf to work below */
                        double tt[7];
                        if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
                        {
                                fprintf( ioQQQ, " level2.dat error getting line  %li\n", i );
                                cdEXIT(EXIT_FAILURE);
                        }
                        /* option to skip any line starting with # */
                        if( chLine[0]!='#' )
                        {
                                /*printf("string is %s\n",chLine );*/
                                sscanf( chLine , "%lf %lf %lf %lf %lf %lf %lf " , 
                                        &tt[0] ,
                                        &tt[1] ,
                                        &tt[2] ,
                                        &tt[3] ,
                                        &tt[4] ,
                                        &tt[5] ,
                                        &tt[6] );
                                /* these are readjusted into their final form in the structure 
                                 * in routine lines_setup*/
                                TauLine2[i].Hi->nelem = (int)tt[0];
                                TauLine2[i].Hi->IonStg = (int)tt[1];
                                TauLine2[i].Lo->g = (realnum)tt[2];
                                TauLine2[i].Hi->g = (realnum)tt[3];
                                TauLine2[i].Emis->gf = (realnum)tt[4];
                                TauLine2[i].EnergyWN = (realnum)tt[5];
                                cs1_flag_lev2[i] = (realnum)tt[6];
                                ++i;
                        }
                }
                /* get magic number off last line */
                if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
                {
                        fprintf( ioQQQ, " level2.dat error getting last magic number\n" );
                        cdEXIT(EXIT_FAILURE);
                }
                sscanf( chLine , "%ld" , &magic2 );
                if( 999 != magic2 )
                {
                        fprintf( ioQQQ, " level2.dat ends will wrong magic number=%ld \n", 
                          magic2 );
                        cdEXIT(EXIT_FAILURE);
                }
                fclose( ioDATA );
                if( trace.lgTrace )
                        fprintf( ioQQQ," OK \n");

                /* end of level2 block*/
        }

        /*************************
         *  
         * get the UTA line sets   
         * >>chng 06 nov 26, use Gu et al. 2006 UTA data by default
         * with option to use older Behar
         *  
         *************************/

        /* these will count number of UTA lines of each type */
        nUTA = 0;
        nUTA_Gu06 = 0;
        nUTA_Behar = 0;
        nUTA_Romas = 0;
        /* this is option to use older Behar data rather than new Gu data
         * variable lgInnerShell_Gu06 is true by default, set false with 
         * SET UTA BEHAR command */
        if( ionbal.lgInnerShell_Gu06 )
        {
                if( trace.lgTrace )
                        fprintf( ioQQQ," atmdat_readin opening UTA_Gu06.dat:");

                ioGU06 = open_data( "UTA_Gu06.dat", "r" );

                /* count how many non-comment lines there are in Gu06 file */
                while( read_whole_line( chLine , (int)sizeof(chLine) , ioGU06 ) != NULL )
                {
                        /* we want to count the lines that do not start with #
                        * since these contain data */
                        if( chLine[0] != '#')
                                ++nUTA_Gu06;
                }
                /* we counted a magic number as a real line, back off by one */
                ASSERT( nUTA_Gu06 > 0 );
                --nUTA_Gu06;

                /* now rewind the file so we can read it a second time*/
                if( fseek( ioGU06 , 0 , SEEK_SET ) != 0 )
                {
                        fprintf( ioQQQ, " atmdat_readin could not rewind UTA_Gu06.dat.\n");
                        cdEXIT(EXIT_FAILURE);
                }

                /* get up to magic number in Gu 06 file - there is a large header of 
                 * comments with the first data the magic number */
                /* count how many non-comment lines there are in Gu06 file */
                while( read_whole_line( chLine , (int)sizeof(chLine) , ioGU06 ) != NULL )
                {
                        /* we want to break on first line that does not start with #
                        * since that contains data */
                        if( chLine[0] != '#')
                                break;
                }
                /* now get Gu magic number */
                j = 1;
                nelem = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
                nelec = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
                ion = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);

                /* is magic number correct?
                 * the following is the set of numbers that appear at the start of UTA_Gu06.dat 2006 11 17 */
                if( ( nelem != 2007 ) || ( nelec != 1 ) || ( ion != 23 ) )
                {
                        fprintf( ioQQQ, 
                                " atmdat_readin: the version of UTA_Gu06.dat is not the current version.\n" );
                        fprintf( ioQQQ, 
                                " I expected to find the number 2007 1 23 and got %li %li %li instead.\n" ,
                                nelem , nelec , ion );
                        fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
                        cdEXIT(EXIT_FAILURE);
                }
        }
        else
        {
                /* this branch, get the Behar 01 data */
                /* >>refer      Fe      inner shell UTA Behar, E., Sako, M, Kahn, S.M., 
                 * >>refercon   2001, ApJ, 563, 497-504
                 * the fits were based on the paper
                 * >>refer      Fe      inner shell UTA Behar, E., & Netzer, H., 2002, ApJ, 570, 165-170 */

                if( trace.lgTrace )
                        fprintf( ioQQQ," atmdat_readin opening UTA_Behar.dat:");

                ioBEHAR = open_data( "UTA_Behar.dat", "r" );

                /* count number of lines in Behar file 
                 * first line is a version number and does not count */
                if( read_whole_line( chLine , (int)sizeof(chLine) , ioBEHAR ) == NULL )
                {
                        fprintf( ioQQQ, " atmdat_readin could not read first line of UTA_Behar.dat.\n");
                        cdEXIT(EXIT_FAILURE);
                }

                j = 1;
                /* UTA_Behar.dat magic number */
                nelem = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
                nelec = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
                ion = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);

                /* magic number
                * the following is the set of numbers that appear at the start of UTA_Behar.dat 2002 8 19 */
                /* >>chng 02 oct 22, had incorrect stat weights for all stages of ionization except atom.
                * no longer used stored values when stat weight is zero */
                if( ( nelem != 2002 ) || ( nelec != 10 ) || ( ion != 22 ) )
                {
                        fprintf( ioQQQ, 
                                " atmdat_readin: the version of UTA_Behar.dat is not the current version.\n" );
                        fprintf( ioQQQ, 
                                " I expected to find the number 2002 10 22 and got %li %li %li instead.\n" ,
                                nelem , nelec , ion );
                        fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
                        cdEXIT(EXIT_FAILURE);
                }
        }

        if( trace.lgTrace )
                fprintf( ioQQQ," atmdat_readin UTA data files open OK\n");

        /* count how many lines are in the file, ignoring all lines
         * starting with '#' */
        nUTA_Behar = 0;

        /* first fill in Behar & Netzer interpolation - first pass just count number of lines */
        for( ipISO=ipLI_LIKE; ipISO<=ipF_LIKE; ++ipISO )
        {
                for( nelem=ipISO; nelem<LIMELM; ++nelem )
                {
                        ++nUTA_Behar;
                }
        }

        if( !ionbal.lgInnerShell_Gu06 )
        {
                /* count number of lines in Behar file if we are going to use it
                 * wee always use above inner shell UTA from B&N */
                while( read_whole_line( chLine , (int)sizeof(chLine) , ioBEHAR ) != NULL )
                {
                        /* we want to count the lines that do not start with #
                         * since these contain data */
                        if( chLine[0] != '#')
                                ++nUTA_Behar;
                }
                /* now rewind the file so we can read it a second time*/
                if( fseek( ioBEHAR , 0 , SEEK_SET ) != 0 )
                {
                        fprintf( ioQQQ, " atmdat_readin could not rewind UTA_Behar.dat.\n");
                        cdEXIT(EXIT_FAILURE);
                }
                /* first line is a version number and does not count */
                if( read_whole_line( chLine , (int)sizeof(chLine) , ioBEHAR ) == NULL )
                {
                        fprintf( ioQQQ, " atmdat_readin could not read first line of UTA_Behar.dat.\n");
                        cdEXIT(EXIT_FAILURE);
                }
        }

        /* now read Romas Kisielius data, taken from
        * >>refer       Fe      UTA     Kisielius, R., Hibbert, A., Ferland, G.J., Foord, M.E., Rose, S.J.,
        * >>refercon    van Hoof, P.A.M. & Keenan, F.P. 2003, MNRAS, 344, 696 */

        if( trace.lgTrace )
                fprintf( ioQQQ," atmdat_readin opening UTA_Kisielius.dat:");

        ioROMAS = open_data( "UTA_Kisielius.dat", "r" );

        nUTA_Romas = 0;
        while( read_whole_line( chLine , (int)sizeof(chLine) , ioROMAS ) != NULL )
        {
                /* we want to count the lines that do not start with #
                 * since these contain data */
                if( chLine[0] != '#')
                {
                        ++nUTA_Romas;
                }
        }

        /* now malloc the UTA lines array with enough space for both
         * data sets */
        UTALines = (transition *)MALLOC( (size_t)(nUTA_Behar+nUTA_Romas+nUTA_Gu06)*sizeof(transition) );

        for( i=0; i<nUTA_Behar+nUTA_Romas+nUTA_Gu06; i++ )
        {
                TransitionJunk( &UTALines[i] );
                UTALines[i].Hi = AddState2Stack();
                UTALines[i].Lo = AddState2Stack();
                UTALines[i].Emis = AddLine2Stack( true );
        }

        /* these will be used to keep track of what is in each data set */
        for( ion=0; ion<ipIRON; ++ion )
        {
                nFeIonRomas[ion] = 0;
                WLloRomas[ion] = 1e10;
                WLhiRomas[ion] = 0.;

                nFeIonBehar[ion] = 0;
                WLloBehar[ion] = 1e10;
                WLhiBehar[ion] = 0.;

                nFeIonGu[ion] = 0;
                WLloGu[ion] = 1e10;
                WLhiGu[ion] = 0.;
        }

        /* this will count number of lines, must equal sum of UTAs at end */
        i = 0;
        /* now rewind the file so we can read it a second time*/
        if( fseek( ioROMAS , 0 , SEEK_SET ) != 0 )
        {
                fprintf( ioQQQ, " atmdat_readin could not rewind UTA_Kisielius.dat.\n");
                cdEXIT(EXIT_FAILURE);
        }

        i = 0;
        /* first fill in Behar & Netzer interpolation - now do real thing,
         * these are Nick Abel's fits to the BN data */
        for( ipISO=ipLI_LIKE; ipISO<=ipF_LIKE; ++ipISO )
        {
                for( nelem=ipISO; nelem<LIMELM; ++nelem )
                {
                        double f1;
                        double sum_spon_auto;
#                       define N 7

                        double alam[N]={ 1.859035 , 4.692138 ,  10.324543 , 19.294364 ,
                                40.401292  ,  44.442754 , 72.959719 };
                        double blam[N]={-9.7855968,-11.159739, -12.914931 , -14.987272,
                                -18.853446 , -19.271781 ,-23.828388 };
                        double clam[N]={ 8.2874628, 8.3043002,  8.3164038 ,  8.3269937,
                                8.4312895  , 8.3730893  , 8.493802 };
                        /* >>chng 05 mar 07, by NA, update fits to better form */

                        double aA[N] = {1.9067 , 1.8715 ,  2.1033, 2.2815 ,
                                1.9511 , 1.9966 ,  2.0852};
                        double bA[N] = {8.4538 , 8.5528, 7.616, 6.9247 ,
                                7.9479, 7.9777 , 7.8349 };

                        double aDep[N] = {29.54 , 12.07 , 24.23 , 7.35 , 7.37 , 11.14 , 11.99 };
                        double bDep[N] = {-14.853 , -7.606 , -7.844 ,  9.987 , 10.503 , 8.541 , 8.865 };

                        /* NB - these are plus one is because 1 will be subtracted when used */
                        /* derive element and ion */
                        UTALines[i].Hi->nelem = nelem+1;
                        ion = nelem + 1 - ipISO;
                        UTALines[i].Hi->IonStg = ion;

                        /* these were the statistical weights */
                        UTALines[i].Hi->g = 4.f;
                        UTALines[i].Lo->g = 2.f;

                        f1 = alam[ipISO-2]*1e-4 + blam[ipISO-2]*1e-4*(nelem+1) + 
                                clam[ipISO-2]*1e-4*POW2(nelem+1.);
                        f1 = (1./f1);
                        UTALines[i].WLAng = (realnum)f1;
                        UTALines[i].EnergyWN = (realnum)(1e8/f1);
                        UTALines[i].EnergyErg = (realnum)(ERG1CM)*UTALines[i].EnergyWN;

                        f1 = aA[ipISO-2]*log((double)(nelem+1)) + bA[ipISO-2];  

                        UTALines[i].Emis->Aul = (realnum)pow(10., f1 );
                        /* generate Emis.gf from A */
                        UTALines[i].Emis->gf = 
                                (realnum)(UTALines[i].Emis->Aul*UTALines[i].Hi->g*
                                1.4992e-8*POW2(1e4/UTALines[i].EnergyWN));

                        UTALines[i].Emis->iRedisFun = 1;

                        /* find sum of spontaneous relax and autoionization As */
                        f1 = log((double)(nelem+1));
                        f1 = POW2( f1 );
                        sum_spon_auto = aDep[ipISO-2]*1e-2*f1 + bDep[ipISO-2]*1e-1;

                        if( ipISO == ipBE_LIKE )
                        {
                                sum_spon_auto = exp( sum_spon_auto );
                                sum_spon_auto = pow(10., sum_spon_auto )*1e13;
                        }
                        else if( ipISO <= ipF_LIKE )
                        {
                                sum_spon_auto = pow(10., sum_spon_auto )*1e13;
                        }
                        else
                                TotalInsanity();

                        /* last number is negative branching ratio for autoionization, 
                         * which we will store as negative of cs so not misinterpreted as cs */
                        /* >>chng 03 feb 18, change to negative of br */

                        /*ASSERT( UTALines[i].Emis->Aul / sum_spon_auto <= 1. && UTALines[i].Aul / sum_spon_auto >= 0. );*/
                        UTALines[i].Coll.cs = (realnum)MIN2(0., UTALines[i].Emis->Aul / sum_spon_auto - 1. );

                        /* option to print UTAs */
#                       if 0
                        fprintf(ioQQQ," %2s%2li\t%.4f\t%.3e\t%.3e\t%.3e\n",
                                elementnames.chElementSym[nelem], 
                                ion, 
                                UTALines[i].WLAng ,
                                UTALines[i].Emis->Aul ,
                                sum_spon_auto , 
                                -UTALines[i].cs);
#                       endif

                        /* finally increment i */
                        ++i;
#                       undef N
                }
        }

        /* next read in the Gu file */
        if( ionbal.lgInnerShell_Gu06 )
        {
                ioDATA = ioGU06;
                while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
                {
                        /* only look at lines without '#' in first col */
                        if( chLine[0] != '#')
                        {
                                long int i2;
                                double WLAng, Aul, oscill , Aauto;

                                sscanf( chLine ,"%4li%5li%8lf%13le%12le",
                                        &ion,&i2,&WLAng,&Aul,&Aauto);
                                sscanf( &chLine[54] ,"%13le", &oscill);

                                /* all these are iron, first number was ion stage with 0 the atom */
                                UTALines[i].Hi->nelem = ipIRON+1;

                                /* now do stage of ionization */
                                ASSERT( ion >= 0 && ion < ipIRON );
                                /* the ion stage - 1 is atom - this data file has different
                                 * format from other two - others are 0 for atom */
                                UTALines[i].Hi->IonStg = ion;

                                /* these were the statistical weights 
                                 * not included in the original data files, and not needed,
                                 * so use unity */
                                UTALines[i].Hi->g = 1.f;
                                UTALines[i].Lo->g = 1.f;

                                /* wavelength in Angstroms */
                                UTALines[i].WLAng = (realnum)WLAng;

                                /* keep track of what we have */
                                ++nFeIonGu[ion-1];
                                WLloGu[ion-1] = MIN2( WLAng , WLloGu[ion-1] );
                                WLhiGu[ion-1] = MAX2( WLAng , WLloGu[ion-1] );

                                /* energy in wavenumbers */
                                UTALines[i].EnergyWN = (realnum)(1e8/WLAng);
                                UTALines[i].EnergyErg = (realnum)(ERG1CM)*UTALines[i].EnergyWN;

                                /* this is true spontaneous rate for doubly excited state to inner 
                                 * shell UTA, and is used for pumping, and also relaxing back to inner shell */
                                UTALines[i].Emis->Aul = (realnum)Aul;

                                /* cs is branching ratio for autoionization, 
                                 * which we will store as negative cs so not misinterpreted as real cs
                                 * this multiplies the direct pump rate to get ionization rate 
                                 * rate evaluated in conv_base.cpp*/
                                double frac_ioniz = Aauto/(Aul + Aauto);
                                ASSERT( frac_ioniz >=0. &&  frac_ioniz <=1. );
                                UTALines[i].Coll.cs = -(realnum)frac_ioniz;

                                /* save gf scanned from line */
                                UTALines[i].Emis->gf = UTALines[i].Lo->g * (realnum)oscill;

                                /* statistical weights are not given in the data file -
                                 * assume unity and convert absorption oscillator strength
                                 * into an effective Aul, which is correct for unit statistical
                                 * weight */
                                UTALines[i].Emis->Aul = (realnum)eina(UTALines[i].Emis->gf,
                                        UTALines[i].EnergyWN,UTALines[i].Hi->g);

                                UTALines[i].Emis->iRedisFun = 1;
                                {
                                        /* DEBUG OTS rates - turn this on to debug line, continuum or both rates */
                                        /*@-redef@*/
                                        enum {DEBUG_LOC=false};
                                        /*@+redef@*/
                                        if( DEBUG_LOC && UTALines[i].Hi->nelem==ipIRON+1 && UTALines[i].Hi->IonStg==9)
                                        {
                                                fprintf(ioQQQ,"DEBUG Gu UTA\t%li\t%.3f\t%.3e\t%.3e\t%.3e\n",
                                                        ion , WLAng,Aul , 
                                                        eina(UTALines[i].Emis->gf,
                                                        UTALines[i].EnergyWN,UTALines[i].Hi->g),
                                                        Aauto );
                                        }
                                }

                                /* finally increment i */
                                ++i;
                        }
                }
        }
        else
        {
                /* ionbal.lgInnerShell_Gu06 is false, do not use this data */
                nUTA_Gu06 = 0;

                /* the Behar 01 data */
                /* read in the Behar data */
                ioDATA = ioBEHAR;
                while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
                {
                        /* only look at lines without '#' in first col */
                        if( chLine[0] != '#')
                        {
                                long int i1, i2, i3;
                                double f1, f2, oscill;
                                double frac_relax;

                                sscanf( chLine ,"%li\t%li\t%li\t%lf\t%le\t%le\t%le",
                                        &i1,&i2,&i3,&f1,&f2,&frac_relax,&oscill );

                                /* all these are iron, first number was ion stage with 0 the atom */
                                UTALines[i].Hi->nelem = ipIRON+1;

                                /* now do stage of ionization */
                                ASSERT( i1 >= 0 && i1 < ipIRON );
                                /* NB - the plus one is because 1 will be subtracted when used,
                                * in the original data file i1 is 0 for the atom */
                                UTALines[i].Hi->IonStg = i1 + 1;

                                /* these were the statistical weights */
                                if( i2>0 && i3>0 )
                                {
                                        UTALines[i].Hi->g = (realnum)i2;
                                        UTALines[i].Lo->g = (realnum)i3;
                                }
                                else
                                {
                                        /* nearly all stages of ionization do not include the stat weights,
                                        * so use unity */
                                        UTALines[i].Hi->g = 1.f;
                                        UTALines[i].Lo->g = 1.f;
                                }

                                /* wavelength in Angstroms */
                                UTALines[i].WLAng = (realnum)f1;

                                /* keep track of what we have */
                                ++nFeIonBehar[i1];
                                WLloBehar[i1] = MIN2( f1 , WLloBehar[i1] );
                                WLhiBehar[i1] = MAX2( f1 , WLhiBehar[i1] );

                                /* energy in wavenumbers */
                                UTALines[i].EnergyWN = (realnum)(1e8/f1);
                                UTALines[i].EnergyErg = (realnum)(ERG1CM)*UTALines[i].EnergyWN;

                                /* this is true spontaneous rate for doubly excited state to inner shell UTA,
                                * and is used for pumping, and also relaxing back to inner shell */
                                UTALines[i].Emis->Aul = (realnum)f2;

                                /* >>chng 02 oct 22, use absorption oscillator strength in data file
                                * since don't have stat weights for most lines */
                                UTALines[i].Emis->gf = UTALines[i].Lo->g * (realnum)oscill;

                                UTALines[i].Emis->iRedisFun = 1;

                                /* last number is branching ratio for autoionizing, 
                                * which we will store in negative in cs so not misinterpreted as real cs*/
                                /* >>chng 03 feb 18, cs to negative of branching ratio */
                                ASSERT( frac_relax >=0. &&  frac_relax <=1. );
                                UTALines[i].Coll.cs = (realnum)(frac_relax-1.);
#                               if 0
                                fprintf(ioQQQ,"data set %li line %li %2s%2li\t%.4f\t%.3e\t%.3e\n",
                                        nDataSet,
                                        i,
                                        elementnames.chElementSym[UTALines[i].Hi->nelem-1], 
                                        UTALines[i].Hi->IonStg, 
                                        UTALines[i].WLAng ,
                                        UTALines[i].Emis->Aul ,
                                        -UTALines[i].cs);
#                               endif

                                /* finally increment the counter i */
                                ++i;
                        }
                }
        }

        /* last read in the Romas Kisielius data last since will use it in preference to other
         * two sets 
         *>>refer       Fe      UTA             Kisielius, R.; Hibbert, A.; Ferland, G. J.; Foord, M. E.; Rose, S. J.; 
         *>>refercon    van Hoof, P. A. M.; Keenan, F. P. 2003, MNRAS, 344, 696 
         */
        ioDATA = ioROMAS;
        while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
        {
                /* only look at lines without '#' in first col */
                if( chLine[0] != '#')
                {
                        long int i1, i2, i3;
                        double f1, f2, oscill;
                        double frac_relax;

                        sscanf( chLine ,"%li\t%li\t%li\t%lf\t%le\t%le\t%le",
                                &i1,&i2,&i3,&f1,&f2,&frac_relax,&oscill );

                        /* all these are iron, first number was ion stage with 0 the atom */
                        UTALines[i].Hi->nelem = ipIRON+1;

                        /* now do stage of ionization */
                        ASSERT( i1 >= 0 && i1 < ipIRON );
                        /* NB - the plus one is because 1 will be subtracted when used,
                        * in the original data file i1 is 0 for the atom */
                        UTALines[i].Hi->IonStg = i1 + 1;

                        /* these were the statistical weights */
                        if( i2>0 && i3>0 )
                        {
                                UTALines[i].Hi->g = (realnum)i2;
                                UTALines[i].Lo->g = (realnum)i3;
                        }
                        else
                        {
                                /* nearly all stages of ionization do not include the stat weights,
                                * so use unity */
                                UTALines[i].Hi->g = 1.f;
                                UTALines[i].Lo->g = 1.f;
                        }

                        /*TauLines[i].Hi->IonStg = i2;*/
                        UTALines[i].WLAng = (realnum)f1;

                        /* keep track of what we have */
                        ++nFeIonRomas[i1];
                        WLloRomas[i1] = MIN2( f1 , WLloRomas[i1] );
                        WLhiRomas[i1] = MAX2( f1 , WLhiRomas[i1] );

                        UTALines[i].EnergyWN = (realnum)(1e8/f1);
                        UTALines[i].EnergyErg = (realnum)(ERG1CM)*UTALines[i].EnergyWN;

                        /* this is true spontaneous rate for doubly excited state to inner shell,
                        * and is used for pumping, and also relaxing back to inner shell */
                        UTALines[i].Emis->Aul = (realnum)f2;
                        /* generate gf from A */
                        /* >>chng 02 oct 22, use absorption oscillator strength in data file
                        * since don't have stat weights for most lines */
                        UTALines[i].Emis->gf = UTALines[i].Lo->g * (realnum)oscill;

                        UTALines[i].Emis->iRedisFun = 1;

                        /* last number is branching ratio for autoionizing, 
                        * which we will store in negative in cs so not misinterpreted as real cs*/
                        /* >>chng 03 feb 18, cs to negative of br */
                        ASSERT( frac_relax >=0. &&  frac_relax <=1. );
                        UTALines[i].Coll.cs = (realnum)(frac_relax-1.);

                        /* finally increment i */
                        ++i;
                }
        }

        /* these have to agree */
        ASSERT( i == (nUTA_Behar+nUTA_Romas+nUTA_Gu06) );

        if( trace.lgTrace )
        {
                fprintf(ioQQQ , "summary of UTA data sets read in\nion numb     WLlo    WLhi\n");
                fprintf(ioQQQ , "Behar 01 total lines %li\n" , nUTA_Behar );
                for( ion=0; ion<ipIRON;++ion )
                {
                        if( nFeIonBehar[ion] )
                        {
                                fprintf(ioQQQ, "%3li %6li %7.3f %7.3f \n", 
                                        ion,
                                        nFeIonBehar[ion],
                                        WLloBehar[ion],
                                        WLhiBehar[ion]);
                        }
                }

                fprintf(ioQQQ , "Gu 06 total lines %li\n" , nUTA_Gu06 );
                for( ion=0; ion<ipIRON;++ion )
                {
                        if( nFeIonGu[ion] )
                        {
                                fprintf(ioQQQ, "%3li %6li %7.3f %7.3f \n", 
                                        ion,
                                        nFeIonGu[ion],
                                        WLloGu[ion],
                                        WLhiGu[ion]);
                        }
                }

                fprintf(ioQQQ , "Romas Kisielius 03 total lines %li\n", nUTA_Romas );
                for( ion=0; ion<ipIRON;++ion )
                {
                        if( nFeIonRomas[ion] )
                        {
                                fprintf(ioQQQ, "%3li %6li %7.3f %7.3f \n", 
                                        ion,
                                        nFeIonRomas[ion],
                                        WLloRomas[ion],
                                        WLhiRomas[ion] );
                        }
                }
        }

        /* this is default option to use the data Romas Kisielius computed
         * to replace any existing data .  */

        /* here there are two options - we can add the Romas Kisielius data to the
        * first block of UTA data, but zero out any overlapping data,
        * or simply ignore the Romas Kisielius data */
        if( ionbal.lgInnerShell_Kisielius )
        {
                for( i=0; i<nUTA_Behar+nUTA_Gu06; ++i )
                {
                        if( UTALines[i].Hi->nelem-1 == ipIRON )
                        {
                                /* nFeIonRomas counts how many of the Romas Kisielius lines
                                 * fall in this stage of ionization - if zero then
                                 * none, so keep Behar data, if positive, then
                                 * there are new data, and zero out Behar data */
                                ASSERT( UTALines[i].Hi->IonStg-1< ipIRON );
                                if( nFeIonRomas[UTALines[i].Hi->IonStg-1] )
                                {
                                        UTALines[i].Emis->Aul = -1.;
                                }
                        }
                }
                /* the total number of UTAs, a global variable */
                nUTA = nUTA_Behar + nUTA_Gu06 + nUTA_Romas;
        }
        else
        {
                /* the Kisielius data come last, and nUTA is the global variable 
                 * that counts the number of UTA lines.  Setting nUTA to only the
                 * first two amounts to ignoring the Kisielius data */
                nUTA = nUTA_Behar + nUTA_Gu06;
        }

        /* now close the files */
        fclose(ioROMAS);
        if( ionbal.lgInnerShell_Gu06 )
                fclose(ioGU06);
        else
                fclose(ioBEHAR);

        /* end UTA */

        /* allocate space for the carbon monoxide CO rotation levels */
        ASSERT( nCORotate > 0);
        C12O16Rotate = (transition *)MALLOC( (size_t)nCORotate*sizeof(transition) );
        C13O16Rotate = (transition *)MALLOC( (size_t)nCORotate*sizeof(transition) );

        /* this says we cannot change the number of levels again, since space is allocated */
        lgCORotateMalloc = true;

        /* next initialize entire array to dangerously large negative numbers */
        for( J=0; J< nCORotate; ++J )
        {
                TransitionJunk( &C12O16Rotate[J] );
                C12O16Rotate[J].Hi = AddState2Stack();
                C12O16Rotate[J].Lo = AddState2Stack();
                C12O16Rotate[J].Emis = AddLine2Stack( true );

                TransitionJunk( &C13O16Rotate[J] );
                C13O16Rotate[J].Hi = AddState2Stack();
                C13O16Rotate[J].Lo = AddState2Stack();
                C13O16Rotate[J].Emis = AddLine2Stack( true );
        }

        /* read in data for the set of hyperfine structure lines, and allocate
         * space for the transition HFLines[nHFLines] structure */
        HyperfineCreate();

        /* this is Humeshkar's work on generic atoms and molecules. */
        //Nemala_Start();

        /* initialize the large block of level 1 real lines, and OP level 2 lines */
        lines_setup();

        /* mewe_gbar.dat mewe_gbar.dat mewe_gbar.dat mewe_gbar.dat mewe_gbar.dat mewe_gbar.dat ========*/
        /* read in g-bar data taken from
         *>>refer       all     cs-gbar Mewe, R.; Gronenschild, E. H. B. M.; van den Oord, G. H. J., 1985,
         *>>refercon    A&AS, 62, 197 */
        /* open file with Mewe coefficients */

        if( trace.lgTrace )
                fprintf( ioQQQ," atmdat_readin opening mewe_gbar.dat:");

        ioDATA = open_data( "mewe_gbar.dat", "r" );

        /* get magic number off first line */
        if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
        {
                fprintf( ioQQQ, " mewe_gbar.dat error getting magic number\n" );
                cdEXIT(EXIT_FAILURE);
        }
        /* check the number */
        sscanf( chLine , "%ld" , &magic1 );
        if( magic1 != 9101 )
        {
                fprintf( ioQQQ, " mewe_gbar.dat starts with wrong magic number=%ld \n", 
                  magic1 );
                cdEXIT(EXIT_FAILURE);
        }

        /* now get the actual data, indices are correct for c, in Fort went from 2 to 210 */
        for( i=1; i < 210; i++ )
        {
                if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
                {
                        fprintf( ioQQQ, " mewe_gbar.dat error getting line  %li\n", i );
                        cdEXIT(EXIT_FAILURE);
                }
                /*printf("%s\n",chLine);*/
                double help[4];
                sscanf( chLine, "%lf %lf %lf %lf ", &help[0], &help[1], &help[2], &help[3] );
                for( int l=0; l < 4; ++l )
                        MeweCoef.g[i][l] = (realnum)help[l];
        }

        /* get magic number off last line */
        if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
        {
                fprintf( ioQQQ, " mewe_gbar.dat error getting last magic number\n" );
                cdEXIT(EXIT_FAILURE);
        }

        sscanf( chLine , "%ld" , &magic2 );

        if( magic1 != magic2 )
        {
                fprintf( ioQQQ, " mewe_gbar.dat ends will wrong magic number=%ld \n", 
                  magic2 );
                cdEXIT(EXIT_FAILURE);
        }

        fclose( ioDATA );

        if( trace.lgTrace )
                fprintf( ioQQQ," OK \n");

         /* This is what remains of the t_yield initialization
          * this should not be in the constructor of t_yield
          * since it initializes a different struct */
        for( nelem=0; nelem < LIMELM; nelem++ )
                for( ion=0; ion < LIMELM; ion++ )
                        Heavy.nsShells[nelem][ion] = LONG_MAX;

        /* now read in all auger yields
         * will do elements from li on up,
         * skip any line starting with #
         * this loop goes from lithium to Zn */
        for( nelem=2; nelem < LIMELM; nelem++ )
        {
                /* nelem is on the shifted C scale, so 2 is Li */
                for( ion=0; ion <= nelem; ion++ )
                {
                        /* number of bound electrons, = atomic number for neutral */
                        nelec = nelem - ion + 1;
                        /* one of dima's routines to determine the number of electrons
                         * for this species, nelem +1 to shift to physical number */
                        /* subroutine atmdat_outer_shell(iz,in,imax,ig0,ig1)
                         * iz - atomic number from 1 to 30 (integer) 
                         * in - number of electrons from 1 to iz (integer)
                         * Output: imax - number of the outer shell
                         */
                        atmdat_outer_shell(nelem+1,nelec,&imax,&ig0,&ig1);

                        ASSERT( imax > 0 && imax <= 10 );

                        /* nsShells[nelem][ion] is outer shell number for ion with nelec electrons
                         * on physics scale, with K shell being 1 */
                        Heavy.nsShells[nelem][ion] = imax;
                }
        }

        /*************************************************************
         *                                                           *
         * get the Auger electron yield data set                     *
         *                                                           *
         *************************************************************/

        t_yield::Inst().init_yield();

        /****************************************************************
         *                                                              *
         * get the Hummer and Storey model case A, B results, these are *
         * the two data files e1b.dat and e2b.dat, for H and He         *
         *                                                              *
         ****************************************************************/

        /* now get Hummer and Storey case b data, loop over H, He */
        /* >>chng 01 aug 08, generalized this to both case A and B, and all elements
         * up to HS_NZ, now 8 */
        for( ipZ=0; ipZ<HS_NZ; ++ipZ )
        {
                /* don't do the minor elements, Li-B */
                if( ipZ>1 && ipZ<5 ) continue;

                for( iCase=0; iCase<2; ++iCase )
                {
                        /* open file with case B data */
                        /* >>chng 01 aug 08, add HS_ to start of file names to indicate Hummer Storey origin */
                        /* create file name for this charge
                         * first character after e is charge number for this element,
                         * then follows whether this is case A or case B data */
                        sprintf( chFilename, "HS_e%ld%c.dat", ipZ+1, ( iCase == 0 ) ? 'a' : 'b' );

                        if( trace.lgTrace )
                                fprintf( ioQQQ," atmdat_readin opening Hummer Storey emission file %s:",chFilename);

                        /* open the file */
                        ioDATA = open_data( chFilename, "r" );

                        if( trace.lgTrace )
                        {
                                fprintf( ioQQQ," OK\n");
                        }

                        /* read in the number of temperatures and densities*/
                        i = fscanf( ioDATA, "%li %li ", 
                                &atmdat.ntemp[iCase][ipZ], &atmdat.nDensity[iCase][ipZ] );

                        /* check that ntemp and nDensity are below NHSDIM, 
                         * set to 15 in atmdat_HS_caseB.h */
                        assert (atmdat.ntemp[iCase][ipZ] >0 && atmdat.ntemp[iCase][ipZ] <= NHSDIM );
                        assert (atmdat.nDensity[iCase][ipZ] > 0 && atmdat.nDensity[iCase][ipZ] <= NHSDIM);

                        /* loop reading in line emissivities for all temperatures*/
                        for( ipTemp=0; ipTemp < atmdat.ntemp[iCase][ipZ]; ipTemp++ )
                        {
                                for( ipDens=0; ipDens < atmdat.nDensity[iCase][ipZ]; ipDens++ )
                                {
                                        long int junk, junk2 , ne;
                                        fscanf( ioDATA, " %lf %li %lf %c %li %ld ", 
                                                &atmdat.Density[iCase][ipZ][ipDens], &junk , 
                                                &atmdat.ElecTemp[iCase][ipZ][ipTemp], &cha , &junk2 , 
                                                &atmdat.ncut[iCase][ipZ] );

                                        ne = atmdat.ncut[iCase][ipZ]*(atmdat.ncut[iCase][ipZ] - 1)/2;
                                        ASSERT( ne<=NLINEHS );
                                        for( j=0; j < ne; j++ )
                                        {
                                                fscanf( ioDATA, "%lf ", &atmdat.Emiss[iCase][ipZ][ipTemp][ipDens][j] );
                                        }
                                }
                        }

                        /*this is end of read-in loop */
                        fclose(ioDATA); 
#                       if 0
                        /* print list of densities and temperatures */
                        for( ipDens=0; ipDens<atmdat.nDensity[iCase][ipZ]; ipDens++ )
                        {
                                fprintf(ioQQQ," %e,", atmdat.Density[iCase][ipZ][ipDens]);
                        }
                        fprintf(ioQQQ,"\n");
                        for( ipTemp=0; ipTemp<atmdat.ntemp[iCase][ipZ]; ipTemp++ )
                        {
                                fprintf(ioQQQ," %e,", atmdat.ElecTemp[iCase][ipZ][ipTemp]);
                        }
                        fprintf(ioQQQ,"\n");
#                       endif
                }
        }

        /* Verner's model FeII atom
         * do not call if Netzer model used, or it Fe(2) is zero
         * exception is when code is searching for first soln */
        /* read the atomic data from files */
        /* convert line arrays to internal form needed by code */
        FeIICreate();
        return;
}

t_yield::t_yield()
{
         /* preset two arrays that will hold auger data 
          * set to very large values so that code will blow
          * if not set properly below */
        for( int nelem=0; nelem < LIMELM; nelem++ )
        {
                for( int ion=0; ion < LIMELM; ion++ )
                {
                        for( int ns=0; ns < 7; ns++ )
                        {
                                n_elec_eject[nelem][ion][ns] = LONG_MAX;
                                for( int nelec=0; nelec < 10; nelec++ )
                                {
                                        frac_elec_eject[nelem][ion][ns][nelec] = FLT_MAX;
                                }
                        }
                }
        }

        lgKillAuger = false;
}

void t_yield::init_yield()
{
        char chLine[FILENAME_PATH_LENGTH_2];
        const char* chFilename;

        /* following is double for sscanf to work */
        double temp[15];

        DEBUG_ENTRY( "init_yield()" );

        /*************************************************************
         *                                                           *
         * get the Auger electron yield data set                     *
         *                                                           *
         *************************************************************/

        /* NB NB -- This test of Heavy.nsShells remains here to assure
         * that it contains meaningful values since needed below, once
         * t_Heavy is a Singleton, it can be removed !!! */
        ASSERT( Heavy.nsShells[2][0] > 0 );

        /* hydrogen and helium will not be done below, so set yields here*/
        n_elec_eject[ipHYDROGEN][0][0] = 1;
        n_elec_eject[ipHELIUM][0][0] = 1;
        n_elec_eject[ipHELIUM][1][0] = 1;

        frac_elec_eject[ipHYDROGEN][0][0][0] = 1;
        frac_elec_eject[ipHELIUM][0][0][0] = 1;
        frac_elec_eject[ipHELIUM][1][0][0] = 1;

        /* open file auger.dat, yield data file that came from
         * >>refer      all     auger   Kaastra, J.S., & Mewe, R. 1993, A&AS, 97, 443-482 */
        chFilename = "mewe_nelectron.dat";

        if( trace.lgTrace )
                fprintf( ioQQQ, " init_yield opening %s:", chFilename );

        FILE *ioDATA;

        /* open the file */
        ioDATA = open_data( chFilename, "r" );

        /* now read in all auger yields
         * will do elements from li on up,
         * skip any line starting with #
         * this loop goes from lithium to Zn */
        for( int nelem=2; nelem < LIMELM; nelem++ )
        {
                /* nelem is on the shifted C scale, so 2 is Li */
                for( int ion=0; ion <= nelem; ion++ )
                {
                        for( int ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
                        {
                                char ch1 = '#';
                                /* the * is the old comment char, accept it, but really want # */
                                while( ch1 == '#' || ch1 == '*' )
                                {
                                        if( read_whole_line( chLine, (int)sizeof(chLine), ioDATA ) == NULL )
                                        {
                                                fprintf( ioQQQ, " %s error getting line %i\n", chFilename, ns );
                                                cdEXIT(EXIT_FAILURE);
                                        }
                                        ch1 = chLine[0];
                                }
                                /*printf("string is %s\n",chLine );*/
                                sscanf( chLine, "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf", 
                                        &temp[0], &temp[1], &temp[2], &temp[3], &temp[4],
                                        &temp[5], &temp[6], &temp[7], &temp[8], &temp[9],
                                        &temp[10],&temp[11],&temp[12],&temp[13],&temp[14] );
                                n_elec_eject[nelem][ion][ns] = (long int)temp[3];

                                ASSERT( n_elec_eject[nelem][ion][ns] >= 0 && n_elec_eject[nelem][ion][ns] < 11 );

                                /* can pop off up to 10 auger electrons, these are the probabilities*/
                                for( int j=0; j < 10; j++ )
                                {
                                        frac_elec_eject[nelem][ion][ns][j] = (realnum)temp[j+4];
                                        ASSERT( frac_elec_eject[nelem][ion][ns][j] >= 0. );
                                }
                        }
                }
                /* activate this print statement to get yield for k-shell of all atoms */
                /*fprintf(ioQQQ,"yyyield\t%li\t%.4e\n", nelem+1, frac_elec_eject[nelem][0][0][0] );*/
        }

        fclose( ioDATA );

        if( trace.lgTrace )
        {
                /* this is set with no auger command */
                if( lgKillAuger )
                        fprintf( ioQQQ, " Auger yields will be killed.\n");
                fprintf( ioQQQ, " OK\n");
        }

        /* open file mewe_fluor.dat, yield data file that came from
         * >>refer      all     auger   Kaastra, J.S., & Mewe, R. 1993, 97, 443-482 */
        chFilename = "mewe_fluor.dat";

        if( trace.lgTrace )
                fprintf( ioQQQ, " init_yield opening %s:", chFilename );

        /* open the file */
        ioDATA = open_data( chFilename, "r" );

        /* now read in all auger yields
         * will do elements from li on up,
         * skip any line starting with # */
        do 
        {
                if( read_whole_line( chLine, (int)sizeof(chLine), ioDATA ) == NULL )
                {
                        fprintf( ioQQQ, " %s error getting line %i\n", chFilename, 0 );
                        cdEXIT(EXIT_FAILURE);
                }
        }
        while( chLine[0] == '#' );

        bool lgEOL = false;

        nfl_lines = 0;
        do
        {
                const int NKM = 10;
                int nDima[NKM] = { 0, 1, 2, 2, 3, 4, 4, 5, 5, 6 };
                int nAuger;

                if( nfl_lines >= MEWE_FLUOR )
                        TotalInsanity();

                /*printf("string is %s\n",chLine );*/
                sscanf( chLine, "%lf %lf %lf %lf %lf %lf %lf", 
                        &temp[0], &temp[1], &temp[2], &temp[3], &temp[4], 
                        &temp[5], &temp[6] );

                /* the atomic number, C is 5 */
                nfl_nelem[nfl_lines] = (int)temp[0]-1;
                ASSERT( nfl_nelem[nfl_lines] >= 0 && nfl_nelem[nfl_lines] < LIMELM );

                /* the ion stage for target, atom is 0 */
                nfl_ion[nfl_lines] = (int)temp[1]-1;
                ASSERT( nfl_ion[nfl_lines] >= 0 && nfl_ion[nfl_lines] <= nfl_nelem[nfl_lines]+1 );

                /* the target's shell */
                nfl_nshell[nfl_lines] = nDima[(long)temp[2]-1];
                ASSERT( nfl_nshell[nfl_lines] >= 0 && 
                        /* nsShells is shell number, where K is 1 */
                        nfl_nshell[nfl_lines] < Heavy.nsShells[nfl_nelem[nfl_lines]][nfl_ion[nfl_lines]]-1 );

                /* this is the number of Auger electrons ejected */
                nAuger = (int)temp[3];
                /* so this is the spectrum of the photons that are emitted */
                nfl_ion_emit[nfl_lines] = nfl_ion[nfl_lines] + nAuger + 1;
                /* must be gt 0 since at least photoelectron comes off */
                ASSERT( nfl_ion_emit[nfl_lines] > 0 && nfl_ion_emit[nfl_lines] <= nfl_nelem[nfl_lines]+1);

                /* this is the type of line as defined in their paper */
                nfl_nLine[nfl_lines] = (int)temp[4];

                /* energy in Ryd */
                fl_energy[nfl_lines] = (realnum)temp[5] / (realnum)EVRYD;
                ASSERT( fl_energy[nfl_lines] > 0. );

                /* fluor yield */
                fl_yield[nfl_lines] = (realnum)temp[6];
                /* NB cannot assert <=1 since data file has yields around 1.3 - 1.4 */
                ASSERT( fl_yield[nfl_lines] >= 0 );

                ++nfl_lines;

                do
                {
                        if( read_whole_line( chLine, (int)sizeof(chLine), ioDATA ) == NULL )
                                lgEOL = true;
                } 
                while( chLine[0]=='#' && !lgEOL );
        } 
        while( !lgEOL );

        fclose( ioDATA );
}

---