/home66/gary/public_html/cloudy/c13_branch/source/atmdat_chianti.cpp

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
* others.  For conditions of distribution and use see copyright notice in license.txt */
#include "cddefines.h"
#include "lines_service.h"
#include "physconst.h"
#include "taulines.h"
#include "trace.h"
#include "string.h"
#include "thirdparty.h"
#include "rfield.h"
#include "atmdat.h"

static const bool DEBUGSTATE = false;
// minimum energy difference (wavenumbers) we will accept
const double ENERGY_MIN_WN = 1e-10;

void atmdat_STOUT_readin( long intNS, char *chPrefix )
{
        DEBUG_ENTRY( "atmdat_STOUT_readin()" );

        long int nMolLevs;
        char chLine[FILENAME_PATH_LENGTH_2],
                chNRGFilename[FILENAME_PATH_LENGTH_2],
                chTPFilename[FILENAME_PATH_LENGTH_2],
                chCOLLFilename[FILENAME_PATH_LENGTH_2];

        static const int MAX_NUM_LEVELS = 999;

        dBaseSpecies[intNS].lgMolecular = false;
        dBaseSpecies[intNS].lgLTE = false;
        dBaseSpecies[intNS].lgLAMDA = false;

        strcpy( chNRGFilename , chPrefix );
        strcpy( chTPFilename , chPrefix );
        strcpy( chCOLLFilename , chPrefix );

        /*Open the energy levels file*/
        strcat( chNRGFilename , ".nrg");
        uncaps( chNRGFilename );
        if(DEBUGSTATE)
                printf("The data file is %s \n",chNRGFilename);

        /*Open the files*/
        if( trace.lgTrace )
                fprintf( ioQQQ," atmdat_STOUT_readin opening %s:",chNRGFilename);

        FILE *ioDATA;
        ioDATA = open_data( chNRGFilename, "r" );
        long int i = 0;
        bool lgEOL = false;
        long index = 0;
        double nrg = 0.0;
        double stwt = 0.0;

        /* first line is a version number - now confirm that it is valid */
        if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
        {
                fprintf( ioQQQ, " atmdat_STOUT_readin could not read first line of %s.\n", chNRGFilename );
                cdEXIT(EXIT_FAILURE);
        }
        i = 1;
        const long int iyr = 11, imo=10 , idy = 14;
        long iyrread, imoread , idyread;
        iyrread = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
        imoread = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
        idyread = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);

        if(( iyrread != iyr ) ||
          (  imoread != imo ) ||
          (  idyread != idy ) )
        {
                fprintf( ioQQQ,
                        " PROBLEM atmdat_STOUT_readin: the version of %s is not the current version.\n", chNRGFilename );
                fprintf( ioQQQ,
                        " atmdat_STOUT_readin: I expected the magic numbers %li %li %li but found %li %li %li.\n",
                        iyr, imo , idy ,iyrread, imoread , idyread  );
                cdEXIT(EXIT_FAILURE);
        }

        nMolLevs = 0;
        //Count number of levels
        while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
        {
                /* # - comment, *** ends data */
                if( chLine[0] != '#' && chLine[0] != '\n' && chLine[0] != '*' )
                {
                        i = 1;
                        long n = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
                        if( n < 0 )
                                break;
                        nMolLevs++;
                }
        }

        /* now rewind the file so we can read it a second time*/
        if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
        {
                fprintf( ioQQQ, " atmdat_STOUT_readin could not rewind %s.\n", chNRGFilename );
                cdEXIT(EXIT_FAILURE);
        }
        //Skip the magic numbers this time
        read_whole_line( chLine , (int)sizeof(chLine) , ioDATA );

        long HighestIndexInFile = nMolLevs;

        nMolLevs = MIN3(atmdat.nStoutMaxLevels,nMolLevs, MAX_NUM_LEVELS );

        if( atmdat.lgStoutPrint == true)
        {
                fprintf( ioQQQ,"     Using STOUT model %s with %li levels of %li available.\n",
                                dBaseSpecies[intNS].chLabel , nMolLevs , HighestIndexInFile );
        }

        dBaseSpecies[intNS].numLevels_max = nMolLevs;
        dBaseSpecies[intNS].numLevels_local = dBaseSpecies[intNS].numLevels_max;

        /*Resize the States array*/
        dBaseStates[intNS].resize(nMolLevs);
        /*Zero out the maximum wavenumber for each species */
        dBaseSpecies[intNS].maxWN = 0.;

        /* allocate the Transition array*/
        ipdBaseTrans[intNS].reserve(nMolLevs);
        for( long ipHi = 1; ipHi < nMolLevs; ipHi++)
                ipdBaseTrans[intNS].reserve(ipHi,ipHi);
        ipdBaseTrans[intNS].alloc();
        dBaseTrans[intNS].resize(ipdBaseTrans[intNS].size());
        dBaseTrans[intNS].states() = &dBaseStates[intNS];
        dBaseTrans[intNS].chLabel() = "Stout";

        //This is creating transitions that we don't have collision data for
        //Maybe use gbar or keep all of the Fe 2 even if it was assumed (1e-5)
        int ndBase = 0;
        for( long ipHi = 1; ipHi < nMolLevs; ipHi++)
        {
                for( long ipLo = 0; ipLo < ipHi; ipLo++)
                {
                        ipdBaseTrans[intNS][ipHi][ipLo] = ndBase;
                        dBaseTrans[intNS][ndBase].Junk();
                        dBaseTrans[intNS][ndBase].setLo(ipLo);
                        dBaseTrans[intNS][ndBase].setHi(ipHi);
                        ++ndBase;
                }
        }


        /* Check for an end of file sentinel */
        bool lgSentinelReached = false;

        //Read the first line of data
        if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
        {
                fprintf( ioQQQ, " atmdat_STOUT_readin could not read first line of the energy level file.\n");
                cdEXIT(EXIT_FAILURE);
        }
        //Read the remaining lines of the energy level file
        do
        {
                i = 1;

                /* Stop on *** */
                if( chLine[0] == '*' )
                {
                        lgSentinelReached = true;
                        break;
                }
                //Comments start with #, skip them
                if( chLine[0] != '#' )
                {
                        /* Skip blank lines */
                        if( chLine[0] == '\n')
                                continue;

                        index = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL) -1 ;

                        if( index < 0 )
                        {
                                fprintf( ioQQQ, " PROBLEM An energy level index was less than 1 in file %s\n",chNRGFilename);
                                fprintf( ioQQQ, " The line being read is between the braces {%.*s}\n",int(strlen(chLine)-1),chLine);
                                cdEXIT(EXIT_FAILURE);
                        }

                        nrg = (double)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
                        stwt = (double)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);

                        if( lgEOL )
                        {
                                fprintf( ioQQQ, " PROBLEM End of line reached prematurely in file %s\n",chNRGFilename);
                                fprintf( ioQQQ, " The line being read is between the braces {%.*s}\n",int(strlen(chLine)-1),chLine);
                                cdEXIT(EXIT_FAILURE);
                        }

                        if( index >= nMolLevs )
                        {
                                //Skip unused levels
                                continue;
                        }

                        dBaseStates[intNS][index].energy().set(nrg,"cm^-1");
                        dBaseStates[intNS][index].g() = stwt;

                        /*information needed for label*/
                        strcpy(dBaseStates[intNS][index].chLabel(), "    ");
                        strncpy(dBaseStates[intNS][index].chLabel(),dBaseSpecies[intNS].chLabel, 4);
                        dBaseStates[intNS][index].chLabel()[4] = '\0';
                }
        }
        while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL );
        if( !lgSentinelReached )
        {
                fprintf( ioQQQ, " PROBLEM End of data sentinel was not reached in file %s\n",chNRGFilename);
                fprintf( ioQQQ, " Check that stars (*****) appear after the last line of data and start in the first column of that line.\n");
                cdEXIT(EXIT_FAILURE);
        }
        fclose(ioDATA);

        if( DEBUGSTATE )
        {
                /*Test whether energy levels are read in properly*/
                printf("DEBUG STOUT ENERGY LEVELS:\n");
                for( int i = 0; i< nMolLevs; i++ )
                {
                        printf("DEBUG\t%i\t%11.4e\t%3.1f\n",i+1,dBaseStates[intNS][i].energy().WN(),dBaseStates[intNS][i].g());
                }
        }

        /* fill in all transition energies, can later overwrite for specific radiative transitions */
        double fenergyWN = 0.;
        for(TransitionList::iterator tr=dBaseTrans[intNS].begin();
                 tr!= dBaseTrans[intNS].end(); ++tr)
        {
                int ipHi = (*tr).ipHi();
                int ipLo = (*tr).ipLo();
                ASSERT(ipHi > ipLo );
                fenergyWN = dBaseStates[intNS][ipHi].energy().WN() - dBaseStates[intNS][ipLo].energy().WN();
                if( fenergyWN <= 0.)
                {
                        printf("\nWARNING: The %s transition between level %i and %i has zero or negative energy.\n",
                                        dBaseStates[intNS][ipHi].chLabel(),ipLo+1,ipHi+1);
                        printf("Check the Stout energy level data file (*.nrg) of this species for missing or duplicate energies.\n");
                        //cdEXIT(EXIT_FAILURE);
                }
                (*tr).EnergyWN() = (realnum)MAX2(ENERGY_MIN_WN,fenergyWN);
                (*tr).WLAng() = (realnum)(1e+8/(*tr).EnergyWN()/RefIndex((*tr).EnergyWN()));
                dBaseSpecies[intNS].maxWN = MAX2(dBaseSpecies[intNS].maxWN,(*tr).EnergyWN());
        }

        /******************************************************
         ************* Transition Probability File ************
         ******************************************************/
        strcat( chTPFilename , ".tp");
        uncaps( chTPFilename );

        ioDATA = open_data( chTPFilename, "r" );

        /* first line is a version number - now confirm that it is valid */
        if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
        {
                fprintf( ioQQQ, " atmdat_STOUT_readin could not read first line of %s.\n", chTPFilename );
                cdEXIT(EXIT_FAILURE);
        }
        i = 1;
        iyrread = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
        imoread = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
        idyread = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);

        if(( iyrread != iyr ) ||
          (  imoread != imo ) ||
          (  idyread != idy ) )
        {
                fprintf( ioQQQ,
                        " PROBLEM atmdat_STOUT_readin: the version of %s is not the current version.\n", chTPFilename );
                fprintf( ioQQQ,
                        " atmdat_STOUT_readin: I expected the magic numbers %li %li %li but found %li %li %li.\n",
                        iyr, imo , idy ,iyrread, imoread , idyread  );
                cdEXIT(EXIT_FAILURE);
        }

        long numtrans = 0;
        //Count number of transitions
        while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
        {
                /* # is comment, *** is end of data */
                if( chLine[0] != '#' && chLine[0] != '\n' && chLine[0] != '*' )
                {
                        i = 1;
                        long n = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
                        if( n < 0 )
                                break;
                        numtrans++;
                }
        }
        /* now rewind the file so we can read it a second time*/
        if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
        {
                fprintf( ioQQQ, " atmdat_STOUT_readin could not rewind %s.\n", chTPFilename );
                cdEXIT(EXIT_FAILURE);
        }
        //Skip the magic numbers this time
        read_whole_line( chLine , (int)sizeof(chLine) , ioDATA );


        //Read the first line of data
        if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
        {
                fprintf( ioQQQ, " atmdat_STOUT_readin could not read first line of the transition probability file.\n");
                cdEXIT(EXIT_FAILURE);
        }

        long ipLo = 0;
        long ipHi = 0;
        double Aul = 0.0;
        lgSentinelReached = false;
        //Read the remaining lines of the transition probability file
        do
        {
                if( chLine[0] == '*' )
                {
                        lgSentinelReached = true;
                        break;
                }

                //Comments start with #, skip them
                if( chLine[0] != '#' )
                {
                        i = 1;

                        /* skip null lines */
                        if( chLine[0] == '\n')
                                continue;

                        //This means last data column has Aul.
                        if( nMatch("A",chLine) )
                        {
                                /* reset read pointer */
                                i = 1;

                                ipLo= (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL) - 1;
                                ipHi = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL) - 1;
                                Aul = (double)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
                                if( lgEOL )
                                {
                                        fprintf( ioQQQ, " PROBLEM End of line reached prematurely in file %s\n",chTPFilename);
                                        fprintf( ioQQQ, " The line being read is between the braces {%.*s}\n",int(strlen(chLine)-1),chLine);
                                        cdEXIT(EXIT_FAILURE);
                                }

                                if( ipLo >= nMolLevs || ipHi >= nMolLevs )
                                {
                                        // skip these lines
                                        continue;
                                }

                                TransitionList::iterator tr = dBaseTrans[intNS].begin()+ipdBaseTrans[intNS][ipHi][ipLo];
                                if( (*tr).hasEmis() )
                                {
                                        fprintf(ioQQQ," PROBLEM duplicate transition found by atmdat_STOUT_readin in %s, "
                                                        "wavelength=%f\n", chTPFilename,dBaseStates[intNS][ipHi].energy().WN()
                                                        - dBaseStates[intNS][ipLo].energy().WN());
                                        cdEXIT(EXIT_FAILURE);
                                }

                                if( (*tr).EnergyWN() > ENERGY_MIN_WN )
                                {
                                        (*tr).AddLine2Stack();
                                        (*tr).Emis().Aul() = Aul;
                                        (*tr).Emis().gf() = (realnum)GetGF((*tr).Emis().Aul(), (*tr).EnergyWN(), (*(*tr).Hi()).g());
                                }
                        }
                        else
                        {
                                fprintf( ioQQQ, " PROBLEM File %s contains an invalid line.\n",chTPFilename);
                                fprintf( ioQQQ, " The line being read is between the braces {%.*s}\n",int(strlen(chLine)-1),chLine);
                                cdEXIT(EXIT_FAILURE);
                        }
                }
        }
        while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL );
        if( !lgSentinelReached )
        {
                fprintf( ioQQQ, " PROBLEM End of data sentinel was not reached in file %s\n",chTPFilename);
                fprintf( ioQQQ, " Check that stars (*****) appear after the last line of data and start in the first column of that line.");
                cdEXIT(EXIT_FAILURE);
        }
        fclose(ioDATA);

        if( DEBUGSTATE )
        {
                /*Test whether stout A's are read in properly */
                printf("DEBUG STOUT A's:\n");
                for(TransitionList::iterator tr=dBaseTrans[intNS].begin();
                                 tr!= dBaseTrans[intNS].end(); ++tr)
                {
                        printf("DEBUG.STOUT.A:\t%i\t%i\t%8.2e\n",(*tr).ipLo()+1,(*tr).ipHi()+1,(*tr).Emis().Aul());
                }
        }



        /******************************************************
         ************* Collision Data File ********************
         ******************************************************/

        strcat( chCOLLFilename , ".coll");
        uncaps( chCOLLFilename );

        ioDATA = open_data( chCOLLFilename, "r" );

        /* first line is a version number - now confirm that it is valid */
        if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
        {
                fprintf( ioQQQ, " atmdat_STOUT_readin could not read first line of %s.\n", chCOLLFilename );
                cdEXIT(EXIT_FAILURE);
        }
        i = 1;
        iyrread = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
        imoread = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
        idyread = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);

        if(( iyrread != iyr ) ||
          (  imoread != imo ) ||
          (  idyread != idy ) )
        {
                fprintf( ioQQQ,
                        " PROBLEM atmdat_STOUT_readin: the version of %s is not the current version.\n", chCOLLFilename );
                fprintf( ioQQQ,
                        " atmdat_STOUT_readin: I expected the magic numbers %li %li %li but found %li %li %li.\n",
                        iyr, imo , idy ,iyrread, imoread , idyread  );
                cdEXIT(EXIT_FAILURE);
        }

        /****** Could add ability to count number of temperature changes, electron CS, and proton CS ****/


        //Read the first line of data
        if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
        {
                fprintf( ioQQQ, " atmdat_STOUT_readin could not read first line of the collision data file.\n");
                cdEXIT(EXIT_FAILURE);
        }

        /* Malloc space for collision strengths */
        StoutCollData[intNS] = (StoutColls***)MALLOC(nMolLevs *sizeof(StoutColls**));
        for( long ipHi=0; ipHi<nMolLevs; ipHi++ )
        {
                StoutCollData[intNS][ipHi] = (StoutColls **)MALLOC((unsigned long)(nMolLevs)*sizeof(StoutColls *));
                for( long ipLo=0; ipLo<nMolLevs; ipLo++ )
                {
                        StoutCollData[intNS][ipHi][ipLo] =
                                (StoutColls *)MALLOC((unsigned long)(ipNCOLLIDER)*sizeof(StoutColls ));

                        for( long k=0; k<ipNCOLLIDER; k++ )
                        {
                                /* initialize all spline variables */
                                StoutCollData[intNS][ipHi][ipLo][k].ntemps = -1;
                                StoutCollData[intNS][ipHi][ipLo][k].temps = NULL;
                                StoutCollData[intNS][ipHi][ipLo][k].collstrs = NULL;
                                StoutCollData[intNS][ipHi][ipLo][k].lgIsRate = false;
                        }
                }
        }

        ipLo = 0;
        ipHi = 0;
        int numpoints = 0;
        double *temps = NULL;
        long ipCollider = -1;
        lgSentinelReached = false;

        //Read the remaining lines of the collision data file
        do
        {
                /* Stop on *** */
                if( chLine[0] == '*' )
                {
                        lgSentinelReached = true;
                        break;
                }

                //Comments start with #, skip them
                if( chLine[0] != '#' )
                {
                        i = 1;

                        /* Skip blank lines */
                        if( chLine[0] == '\n')
                                continue;

                        //This is a temperature line
                        if( nMatch("TEMP",chLine) )
                        {
                                if( temps != NULL)
                                        free( temps );

                                i = 1;
                                numpoints = (int)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
                                ASSERT( numpoints > 0 );

                                temps = (double*)MALLOC(numpoints*sizeof(double));
                                memset( temps, 0, (unsigned long)(numpoints)*sizeof(double) );
                                for( int j = 0; j < numpoints; j++ )
                                {
                                        temps[j] = (double)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
                                        ASSERT( temps[j] > 0 );
                                }
                        }
                        else if( nMatch("CS",chLine) || nMatch("RATE",chLine) )
                        {

                                bool isRate = false;
                                if( nMatch("RATE", chLine) )
                                        isRate = true;

                                if( nMatch( "ELECTRON",chLine ) )
                                {
                                        ipCollider = ipELECTRON;
                                }
                                else if( nMatch( "PROTON",chLine ) )
                                {
                                        ipCollider = ipPROTON;
                                }
                                else
                                {
                                        fprintf( ioQQQ, " PROBLEM atmdat_STOUT_readin: Each line of the collision data"
                                                        "file must specify the collider.\n");
                                        fprintf( ioQQQ, " Possible colliders are: ELECTRON, PROTON\n");
                                        cdEXIT(EXIT_FAILURE);
                                }

                                if( temps == NULL )
                                {
                                        fprintf( ioQQQ, " PROBLEM atmdat_STOUT_readin: The collision "
                                                        "file must specify temperatures before the collision strengths");
                                        cdEXIT(EXIT_FAILURE);
                                }

                                i = 1;
                                ipLo= (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL) - 1;
                                ipHi = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL) - 1;

                                if( ipLo >= nMolLevs || ipHi >= nMolLevs )
                                {
                                        // skip these lines
                                        continue;
                                }

                                /* Set this as a collision strength not a collision rate */
                                StoutCollData[intNS][ipHi][ipLo][ipCollider].lgIsRate = isRate;

                                StoutCollData[intNS][ipHi][ipLo][ipCollider].ntemps = numpoints;
                                ASSERT( numpoints > 0 );

                                /*malloc the space here for the temps and collision strengths*/
                                StoutCollData[intNS][ipHi][ipLo][ipCollider].temps =
                                        (double *)MALLOC((unsigned long)(numpoints)*sizeof(double));
                                StoutCollData[intNS][ipHi][ipLo][ipCollider].collstrs =
                                        (double *)MALLOC((unsigned long)(numpoints)*sizeof(double));
                                /* Loop over all but one CS value */
                                for( int j = 0; j < numpoints; j++ )
                                {
                                        StoutCollData[intNS][ipHi][ipLo][ipCollider].temps[j] = temps[j];
                                        ASSERT( StoutCollData[intNS][ipHi][ipLo][ipCollider].temps[j] > 0 );
                                        StoutCollData[intNS][ipHi][ipLo][ipCollider].collstrs[j] = (double)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
                                        ASSERT( StoutCollData[intNS][ipHi][ipLo][ipCollider].collstrs[j] > 0 );
                                }
                        }
                        else
                        {
                                fprintf( ioQQQ, " PROBLEM File %s contains an invalid line.\n",chCOLLFilename);
                                fprintf( ioQQQ, " The line being read is between the braces {%.*s}\n",int(strlen(chLine)-1),chLine);
                                cdEXIT(EXIT_FAILURE);
                        }

                        if( lgEOL )
                        {
                                fprintf( ioQQQ, " PROBLEM End of line reached prematurely in file %s\n",chCOLLFilename);
                                fprintf( ioQQQ, " The line being read is between the braces {%.*s}\n",int(strlen(chLine)-1),chLine);
                                cdEXIT(EXIT_FAILURE);
                        }
                }
        }
        while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL );
        if( !lgSentinelReached )
        {
                fprintf( ioQQQ, " PROBLEM End of data sentinel was not reached in file %s\n",chCOLLFilename);
                fprintf( ioQQQ, " Check that stars (*****) appear after the last line of data and start in the first column.");
                cdEXIT(EXIT_FAILURE);
        }
        free(temps);
        fclose(ioDATA);

        return;
}

void atmdat_CHIANTI_readin( long intNS, char *chPrefix )
{
        DEBUG_ENTRY( "atmdat_CHIANTI_readin()" );

        int intsplinepts,intTranType,intxs;
        long int nMolLevs,nMolExpLevs,nElvlcLines;// number of experimental and total levels
        FILE *ioElecCollData=NULL, *ioProtCollData=NULL;
        realnum  fstatwt,fenergyWN,fWLAng,fenergy,feinsteina;
        double fScalingParam,fEnergyDiff,*xs,*spl,*spl2;

        char chLine[FILENAME_PATH_LENGTH_2] ,
                chEnFilename[FILENAME_PATH_LENGTH_2],
                chTraFilename[FILENAME_PATH_LENGTH_2],
                chEleColFilename[FILENAME_PATH_LENGTH_2],
                chProColFilename[FILENAME_PATH_LENGTH_2];

        realnum *dBaseStatesEnergy;
        long int *intNewIndex=NULL,*intOldIndex=NULL, *intExperIndex=NULL;
        int interror;
        bool lgProtonData=false,lgEneLevOrder;

        // this is the largest number of levels allowed by the chianti format, I3
        static const int MAX_NUM_LEVELS = 999;

        dBaseSpecies[intNS].lgMolecular = false;
        dBaseSpecies[intNS].lgLAMDA = false;
        dBaseSpecies[intNS].lgLTE = false;

        strcpy( chEnFilename , chPrefix );
        strcpy( chTraFilename , chPrefix );     
        strcpy( chEleColFilename , chPrefix );          
        strcpy( chProColFilename , chPrefix );                  

        /*For the CHIANTI DATABASE*/
        /*Open the energy levels file*/
        strcat( chEnFilename , ".elvlc");
        uncaps( chEnFilename );
        if(DEBUGSTATE)
                printf("The data file is %s \n",chEnFilename);

        /*Open the files*/
        if( trace.lgTrace )
                fprintf( ioQQQ," atmdat_CHIANTI_readin opening %s:",chEnFilename);

        fstream elvlcstream;
        open_data(elvlcstream, chEnFilename,mode_r);

        /*Open the transition probabilities file*/
        strcat( chTraFilename , ".wgfa");
        uncaps( chTraFilename );
        if(DEBUGSTATE)
                printf("The data file is %s \n",chTraFilename);

        /*Open the files*/
        if( trace.lgTrace )
                fprintf( ioQQQ," atmdat_CHIANTI_readin opening %s:",chTraFilename);

        fstream wgfastream;
        open_data(wgfastream, chTraFilename,mode_r);

        /*Open the electron collision data*/
        strcat( chEleColFilename , ".splups");
        uncaps( chEleColFilename );
        if(DEBUGSTATE)
                printf("The data file is %s \n",chEleColFilename);
        /*Open the files*/
        if( trace.lgTrace )
                fprintf( ioQQQ," atmdat_CHIANTI_readin opening %s:",chEleColFilename);

        ioElecCollData = open_data( chEleColFilename, "r" );

        /*Open the proton collision data*/
        strcat( chProColFilename , ".psplups");
        uncaps( chProColFilename );
        if(DEBUGSTATE)
                printf("The data file is %s \n",chProColFilename);
        /*Open the files*/
        if( trace.lgTrace )
                fprintf( ioQQQ," atmdat_CHIANTI_readin opening %s:",chProColFilename);

        /*We will set a flag here to indicate if the proton collision strengths are available */
        if( ( ioProtCollData = fopen( chProColFilename , "r" ) ) != NULL )
        {
                lgProtonData = true;
                //fclose( ioProtCollData );
                //ioProtCollData = NULL;
        }
        else
        {
                lgProtonData = false;
        }

        /*Loop finds how many theoretical and experimental levels are in the elvlc file */
        //eof_col is used get the first 4 charcters per line to find end of file
        const int eof_col = 5;
        //length (+1) of the nrg in the elvlc file
        const int lvl_nrg_col=16;
        //# of columns skipped from the left to get to nrg start
        const int lvl_skipto_nrg = 40;
        /* # of columns to skip from eof check to nrg start */
        const int lvl_eof_to_nrg = lvl_skipto_nrg - eof_col + 1;
        nElvlcLines = 0;
        nMolExpLevs = 1;
        lgEneLevOrder = true;
        if (elvlcstream.is_open())
        {
                int nj = 0;
                char otemp[eof_col];
                char exptemp[lvl_nrg_col];
                double tempexpenergy = 0.;
                /*This loop counts the number of valid rows within the elvlc file
                  as well as the number of experimental energy levels.*/
                while(nj != -1)
                {
                        elvlcstream.get(otemp,eof_col);
                        nj = atoi(otemp);
                        if( nj == -1)
                                break;
                        nElvlcLines++;

                        elvlcstream.seekg(lvl_eof_to_nrg,ios::cur);
                        elvlcstream.get(exptemp,lvl_nrg_col);
                        tempexpenergy = (realnum) atof(exptemp);
                        if( tempexpenergy != 0.)
                                nMolExpLevs++;

                        elvlcstream.ignore(INT_MAX,'\n');

                }
                elvlcstream.seekg(0,ios::beg);
        }

        if(DEBUGSTATE)
        {
                printf("DEBUG: The file %s contains %li experimental energy levels of the %li total levels. \n",chEnFilename,nMolExpLevs,nElvlcLines);
        }

        /* The total number of levels depends on the experimental Chianti switch */
        if( atmdat.lgChiantiExp )
        {
                if( tolower(dBaseSpecies[intNS].chLabel[0]) == 'f' && tolower(dBaseSpecies[intNS].chLabel[1]) == 'e')
                {
                        // Fe is special case with more levels
                        nMolLevs = MIN3(nMolExpLevs, atmdat.nChiantiMaxLevelsFe, MAX_NUM_LEVELS );
                }
                else
                {
                        nMolLevs = MIN3(nMolExpLevs, atmdat.nChiantiMaxLevels, MAX_NUM_LEVELS );
                }
        }
        else
        {
                if( tolower(dBaseSpecies[intNS].chLabel[0]) == 'f' && tolower(dBaseSpecies[intNS].chLabel[1]) == 'e')
                {
                        // Fe is special case with more levels
                        nMolLevs = MIN3(nElvlcLines, atmdat.nChiantiMaxLevelsFe,MAX_NUM_LEVELS );
                }
                else
                {
                        nMolLevs = MIN3(nElvlcLines, atmdat.nChiantiMaxLevels,MAX_NUM_LEVELS );
                }
        }

        if( nMolLevs <= 0 )
        {
                fprintf( ioQQQ, "The number of energy levels is non-positive in datafile %s.\n", chEnFilename );
                fprintf( ioQQQ, "The file must be corrupted.\n" );
                fclose( ioProtCollData );
                cdEXIT( EXIT_FAILURE );
        }

        dBaseSpecies[intNS].numLevels_max = nMolLevs;
        dBaseSpecies[intNS].numLevels_local = dBaseSpecies[intNS].numLevels_max;

        if( atmdat.lgChiantiPrint == true)
        {
                if( atmdat.lgChiantiExp )
                {
                        fprintf( ioQQQ,"     Using CHIANTI model %s with %li experimental energy levels of %li available.\n",
                                dBaseSpecies[intNS].chLabel , nMolLevs , nMolExpLevs );
                }
                else
                {
                        fprintf( ioQQQ,"     Using CHIANTI model %s with %li theoretical energy levels of %li available.\n",
                                dBaseSpecies[intNS].chLabel , nMolLevs , nElvlcLines );
                }
        }

        /*Malloc the energy array to check if the energy levels are in order*/
        dBaseStatesEnergy = (realnum*)MALLOC((unsigned long)(nMolLevs)*sizeof(realnum));

        /*malloc the States array*/
        dBaseStates[intNS].resize(nMolLevs);

        /* allocate the Transition array*/
        ipdBaseTrans[intNS].reserve(nMolLevs);
        for( long ipHi = 1; ipHi < nMolLevs; ipHi++)
                ipdBaseTrans[intNS].reserve(ipHi,ipHi);
        ipdBaseTrans[intNS].alloc();
        dBaseTrans[intNS].resize(ipdBaseTrans[intNS].size());
        dBaseTrans[intNS].states() = &dBaseStates[intNS];
        dBaseTrans[intNS].chLabel() = "Chianti";

        int ndBase = 0;
        for( long ipHi = 1; ipHi < nMolLevs; ipHi++)
        {
                for( long ipLo = 0; ipLo < ipHi; ipLo++)
                {
                        ipdBaseTrans[intNS][ipHi][ipLo] = ndBase;
                        dBaseTrans[intNS][ndBase].Junk();
                        dBaseTrans[intNS][ndBase].setLo(ipLo);
                        dBaseTrans[intNS][ndBase].setHi(ipHi);
                        ++ndBase;
                }
        }

        /*Reading in the energy and checking that they are in order*/
        // C++ io works in terms of cursor movement rather than position on line
        //# of columns to skip over the rydberg energy, we don't use it
        const int lvl_skip_ryd = 15;

        /*Keep track of which levels have experimental data and then create a vector
        which relates their indices to the default chianti energy indices.
         */
        long ncounter = 0;

        if( atmdat.lgChiantiExp )
        {
                //Relate Chianti level indices to a set that only include experimental levels
                intExperIndex = (long int*)MALLOC((unsigned long)(nElvlcLines)* sizeof(long int));

                //Fill relational vector with bad values
                for(int i = 0; i < nElvlcLines; i++)
                {
                        intExperIndex[i] = -1;
                }
        }

        //Read in nrg levels to see if they are in order
        for( long ipLev=0; ipLev<nElvlcLines; ipLev++ )
        {
                if(elvlcstream.is_open())
                {
                        char thtemp[lvl_nrg_col],obtemp[lvl_nrg_col];
                        elvlcstream.seekg(lvl_skipto_nrg,ios::cur);
                        elvlcstream.get(thtemp,lvl_nrg_col);
                        fenergy = (realnum) atof(thtemp);

                        if( atmdat.lgChiantiExp )
                        {
                                /* Go through the entire level list selectively choosing only experimental level energies.
                                 * Store them, not zeroes, in order using ncounter to count the index.
                                 * Any row on the level list where there is no experimental energy, put a -1 in the relational vector.
                                 * If it is a valid experimental energy level store the new ncounter index.
                                 */

                                //Once we collect enough experimental levels stop looking for more.
                                if( ncounter == nMolLevs)
                                        break;
                                ASSERT( ncounter < nMolLevs );

                                if(fenergy != 0. || ipLev == 0 )
                                {
                                        dBaseStatesEnergy[ncounter] = fenergy;
                                        intExperIndex[ipLev] = ncounter;
                                        ncounter++;
                                }
                                else
                                {
                                        intExperIndex[ipLev] = -1;
                                }

                        }
                        else
                        {
                                //Stop looking for levels when the array is full
                                if( ipLev == nMolLevs )
                                        break;

                                elvlcstream.seekg(lvl_skip_ryd,ios::cur);
                                elvlcstream.get(obtemp,lvl_nrg_col);
                                if(atof(obtemp) != 0.)
                                {
                                        fenergy = (realnum) atof(obtemp);
                                }
                                //else
                                        //fprintf( ioQQQ," WARNING: Switched to theoretical energies for species %s, level %3li\n", dBaseSpecies[intNS].chLabel, ipLev );
                                dBaseStatesEnergy[ipLev] = fenergy;
                        }

                        elvlcstream.ignore(INT_MAX,'\n');

                }
                else
                {
                        fprintf( ioQQQ, " The data file %s is corrupted .\n",chEnFilename);
                        fclose( ioProtCollData );
                        cdEXIT(EXIT_FAILURE);
                }
        }

        for( long ipLev=1; ipLev<nMolLevs; ipLev++ )
        {
                /*To check if energy levels are in order*/
                /*If the energy levels are not in order a flag is set*/
                if (dBaseStatesEnergy[ipLev] < dBaseStatesEnergy[ipLev-1])
                {
                        lgEneLevOrder = false;
                }
        }

        // malloc vector for new energy-sorted indices
        intNewIndex = (long int*)MALLOC((unsigned long)(nElvlcLines)* sizeof(long int));

        if(lgEneLevOrder == false)
        {
                /*If the energy levels are not in order and the flag is set(lgEneLevOrder=FALSE)
                then we sort the energy levels to find the relation matrix between the old and new indices*/
                intOldIndex = (long int*)MALLOC((unsigned long)(nMolLevs)* sizeof(long int));
                /*We now do a modified quick sort*/
                spsort(dBaseStatesEnergy,nMolLevs,intOldIndex,2,&interror);
                /*intNewIndex has the key to map*/
                for( long i=0; i<nMolLevs; i++ )
                {
                        ASSERT( intOldIndex[i] < nMolLevs );
                        intNewIndex[intOldIndex[i]] = i;        
                }

                if( nMolLevs < nElvlcLines )
                {
                        /* If any chianti energy levels do not have experimental values
                         * the vector that relates experimental levels to the default will not be the same size
                         * as the stored experimental energy levels.
                         * Therefore this code is required to pad the rest of the sorted vector.
                         * Any index larger than nMolLevs set to -1 */
                        for( long i=nMolLevs; i<nElvlcLines; i++)
                        {
                                intNewIndex[i] = -1;
                        }
                }

                if(DEBUGSTATE)
                {
                        for( long i=0; i<nMolLevs; i++ )
                        {
                                printf("The %ld value of the relation matrix is %ld \n",i,intNewIndex[i]);
                        }
                        for( long i=0; i<nMolLevs; i++ )
                        {
                                printf("The %ld value of the sorted energy array is %f \n",i,dBaseStatesEnergy[i]);
                        }
                }
                free( intOldIndex );
        }
        else
        {
                /* if energies already in correct order, new index is same as original. */
                for( long i=0; i<nMolLevs; i++ )
                {
                        intNewIndex[i] = i;     
                }

                if( nMolLevs < nElvlcLines )
                {
                        /*Pad the sorted vector so that it has the correct number of elements
                         * Any index larger than nMolLevs set to -1.*/
                        for( long i=nMolLevs; i<nElvlcLines; i++)
                        {
                                intNewIndex[i] = -1;
                        }
                }
        }

        /* insist that there the intNewIndex values are unique */
        for( long i=0; i<nMolLevs; i++ )
        {
                for( long j=i+1; j<nMolLevs; j++ )
                {
                        ASSERT( intNewIndex[i] != intNewIndex[j] );
                }
        }

        /*This will print out the relational matrix
         * which includes the original index,the new experimental index,
         * and the sorted experimental index.
         * The numbers NOT be on the C scale, so 1 is the lowest level. */
        if( DEBUGSTATE && atmdat.lgChiantiExp )
        {
                printf("\n\n%s Energy Level Matrix\n",dBaseSpecies[intNS].chLabel);
                printf("(Original, Experimental, Sorted)\n");
                for( long ipLevOld=0; ipLevOld<nElvlcLines; ipLevOld++ )
                {
                        if( intExperIndex[ipLevOld] != -1)
                                printf("%li\t%li\t%li\n",ipLevOld+1,intExperIndex[ipLevOld]+1,intNewIndex[intExperIndex[ipLevOld]]+1);
                }
                printf("\n");
        }

        //After we read in the energies we rewind the energy levels file
        elvlcstream.seekg(0,ios::beg);

        //lvl_skipto_statwt is the # of columns to skip to statwt from left
        const int lvl_skipto_statwt = 37;
        //lvl_statwt_col is the length (+1) of statwt
        const int lvl_statwt_col = 4;
        //Read in stat weight and energy
        for( long ipLevOld=0; ipLevOld<nElvlcLines; ipLevOld++ )
        {
                long ipLevNew = 0;
                if( atmdat.lgChiantiExp )
                {
                        /* We know that non-experimental levels are stored as -1
                         * in the observed/experimental vector.
                         * Use that to skip over those values. */

                        if( intExperIndex[ipLevOld] == -1 )
                        {
                                elvlcstream.ignore(INT_MAX,'\n');
                                continue;
                        }
                        else
                        {
                                /* Store values based on the sorted experimental indices */
                                ipLevNew = intNewIndex[intExperIndex[ipLevOld]];
                        }
                }
                else
                {
                        /* With level trimming on it is possible that there can be rows that
                         * have to be skipped when using theoretical
                         * since the levels no longer exist */
                        if( intNewIndex[ipLevOld] == -1 )
                        {
                                elvlcstream.ignore(INT_MAX,'\n');
                                continue;
                        }
                        else
                        {
                                ipLevNew = intNewIndex[ipLevOld];
                        }
                }

                char gtemp[lvl_statwt_col],thtemp[lvl_nrg_col],obtemp[lvl_nrg_col];
        
                /*information needed for label*/
                strcpy(dBaseStates[intNS][ipLevNew].chLabel(), "    ");
                strncpy(dBaseStates[intNS][ipLevNew].chLabel(),dBaseSpecies[intNS].chLabel, 4);
                dBaseStates[intNS][ipLevNew].chLabel()[4] = '\0';
                
                //Move cursor to and read statwt
                elvlcstream.seekg(lvl_skipto_statwt,ios::cur);
                elvlcstream.get(gtemp,lvl_statwt_col);
                fstatwt = (realnum)atof(gtemp);

                if(fstatwt <= 0.)
                {
                        fprintf( ioQQQ, " WARNING: A positive non zero value is expected for the "
                                         "statistical weight but was not found in %s"
                                        " level %li\n", chEnFilename,ipLevNew);
                        fstatwt = 1.f;
                        //cdEXIT(EXIT_FAILURE);
                }
                dBaseStates[intNS][ipLevNew].g() = fstatwt;

                //Read nrg
                elvlcstream.get(thtemp,lvl_nrg_col);

                fenergy = (realnum) atof(thtemp);

                /* If we are looking for theoretical energies
                 * move over a couple columns before reading in the energy*/
                if( !atmdat.lgChiantiExp )
                {
                        elvlcstream.seekg(lvl_skip_ryd,ios::cur);
                        elvlcstream.get(obtemp,lvl_nrg_col);
                        if(atof(obtemp) != 0.)
                        {
                                fenergy = (realnum) atof(obtemp);
                        }
                }               
                elvlcstream.ignore(INT_MAX,'\n');
                dBaseStates[intNS][ipLevNew].energy().set(fenergy,"cm^-1");
        }
        //Close the level stream
        elvlcstream.close();

        // highest energy transition in chianti
        dBaseSpecies[intNS].maxWN = 0.;
        /* fill in all transition energies, can later overwrite for specific radiative transitions */
        for(TransitionList::iterator tr=dBaseTrans[intNS].begin();
                 tr!= dBaseTrans[intNS].end(); ++tr)
        {
                int ipHi = (*tr).ipHi();
                int ipLo = (*tr).ipLo();
                fenergyWN = (realnum)MAX2( ENERGY_MIN_WN , dBaseStates[intNS][ipHi].energy().WN() - dBaseStates[intNS][ipLo].energy().WN() );

                (*tr).EnergyWN() = fenergyWN;
                (*tr).WLAng() = (realnum)(1e+8/fenergyWN/RefIndex(fenergyWN));
                dBaseSpecies[intNS].maxWN = MAX2(dBaseSpecies[intNS].maxWN,fenergyWN);
        }

        if(DEBUGSTATE)
        {
                /* Print out the stored data for each level.
                 * Level index is NOT on C scale. */
                printf("\n%s Stored Energy Levels\n",dBaseSpecies[intNS].chLabel);
                printf("(Index,Energy in WN, Statistical Weight)\n");
                for( long ipLo = 0; ipLo < nMolLevs; ipLo++ )
                {
                        printf("%li\t%e\t%.1f\n",ipLo+1,dBaseStates[intNS][ipLo].energy().WN(),dBaseStates[intNS][ipLo].g());
                }
        }

        /************************************************************************/
        /*** Read in the level numbers, Einstein A and transition wavelength  ***/
        /************************************************************************/

        //Count the number of rows first
        long wgfarows = -1;
        if (wgfastream.is_open())
        {
                int nj = 0;
                char otemp[eof_col];
                while(nj != -1)
                {
                        wgfastream.get(otemp,eof_col);
                        wgfastream.ignore(INT_MAX,'\n');
                        nj = atoi(otemp);
                        wgfarows++;
                }
                wgfastream.seekg(0,ios::beg);
        }
        else 
                fprintf( ioQQQ, " The data file %s is corrupted .\n",chTraFilename);

        //line_index_col is the length(+1) of the level indexes in the WGFA file
        const int line_index_col = 6;
        //line_nrg_to_eina is the # of columns to skip from wavelength to eina in WGFA file
        const int line_nrg_to_eina =15;
        //line_eina_col is the length(+1) of einsteinA in WGFA
        const int line_eina_col = 16;
        char lvltemp[line_index_col];
        //Start reading WGFA file
        if (wgfastream.is_open())
        {
                for (long i = 0;i<wgfarows;i++)
                {
                        wgfastream.get(lvltemp,line_index_col);
                        long ipLo = atoi(lvltemp)-1;
                        wgfastream.get(lvltemp,line_index_col);
                        long ipHi = atoi(lvltemp)-1;

                        if( atmdat.lgChiantiExp )
                        {
                                /* If either upper or lower index is -1 in the relational vector,
                                 * skip that line in the wgfa file.
                                 * Otherwise translate the level indices.*/
                                if( intExperIndex[ipLo] == - 1 || intExperIndex[ipHi] == -1 )
                                {
                                        wgfastream.ignore(INT_MAX,'\n');
                                        continue;
                                }
                                else
                                {
                                        ipLo = intNewIndex[intExperIndex[ipLo]];
                                        ipHi = intNewIndex[intExperIndex[ipHi]];
                                }
                        }
                        else
                        {
                                /* With level trimming on it is possible that there can be rows that
                                 * have to be skipped when using theoretical
                                * since the levels no longer exist */
                                if( intNewIndex[ipLo] == - 1 || intNewIndex[ipHi] == -1 )
                                {
                                        wgfastream.ignore(INT_MAX,'\n');
                                        continue;
                                }
                                else
                                {
                                        ipLo = intNewIndex[ipLo];
                                        ipHi = intNewIndex[ipHi];
                                }
                        }

                        if( ipLo >= nMolLevs || ipHi >= nMolLevs )
                        {
                                // skip these lines
                                wgfastream.ignore(INT_MAX,'\n');
                                continue;
                        }
        
                        if( ipHi == ipLo )
                        {
                                fprintf( ioQQQ," WARNING: Upper level = lower for a radiative transition in %s. Ignoring.\n", chTraFilename );
                                wgfastream.ignore(INT_MAX,'\n');
                                continue;
                        }

        
                        ASSERT( ipHi != ipLo );
                        ASSERT(ipHi >= 0);
                        ASSERT(ipLo >= 0);

                        // sometimes the CHIANTI datafiles list the highest index first as in the middle of these five lines in ne_10.wgfa:
                        //    ...
                        //    8   10       187.5299      0.000e+00      4.127e+05                 3d 2D1.5 -                  4s 2S0.5           E2
                        //    9   10       187.6573      0.000e+00      6.197e+05                 3d 2D2.5 -                  4s 2S0.5           E2
                        //   11   10   4842624.0000      1.499e-05      9.423e-06                 4p 2P0.5 -                  4s 2S0.5           E1
                        //    1   11         9.7085      1.892e-02      6.695e+11                 1s 2S0.5 -                  4p 2P0.5           E1
                        //    2   11        48.5157      6.787e-02      9.618e+10                 2s 2S0.5 -                  4p 2P0.5           E1
                        //    ...
                        // so, just set ipHi (ipLo) equal to the max (min) of the two indices.
                        // NB NB NB it looks like this may depend upon whether one uses observed or theoretical energies.

                        //Read in wavelengths
                        char trantemp[lvl_nrg_col];
                        wgfastream.get(trantemp,lvl_nrg_col);
                        fWLAng = (realnum)atof(trantemp);

                        /* \todo 2 CHIANTI labels the H 1 2-photon transition as z wavelength of zero.
                         * Should we just ignore all of the wavelengths in this file and use the
                         * difference of level energies instead. */

                        if( atmdat.lgChiantiExp )
                        {
                                /* Sometimes the wgfa file lists the levels to where ipLo > ipHi.
                                 * This seems to be related to the theoretical energies being out of order for those transitions.
                                 * When this is true it seems that the associated wavelength is negative which means they are using theoretical energies,
                                 * even though they have experimental energies.
                                 * For now, print that the levels are switched and ignore the lines if the wavelength is negative. */

                                if( ipHi < ipLo )
                                {
                                        if( strcmp(dBaseSpecies[intNS].chLabel,"Fe 3") == 0)
                                        {
                                                long ipLoTemp = ipLo;
                                                long ipHiTemp = ipHi;
                                                ipHi = max( ipHiTemp, ipLoTemp );
                                                ipLo = min( ipHiTemp, ipLoTemp );
                                                if( atmdat.lgChiantiPrint )
                                                {
                                                        fprintf( ioQQQ," WARNING: Swapped level indices for species %6s, indices %3li %3li with wavelength %e \n",
                                                                dBaseSpecies[intNS].chLabel, ipLoTemp, ipHiTemp,fWLAng );
                                                }
                                        }
                                        else if( atmdat.lgChiantiPrint )
                                        {
                                                fprintf( ioQQQ," WARNING: Upper and Lower Indices are reversed.Species: %6s, Indices: %3li %3li Wavelength: %e \n",
                                                        dBaseSpecies[intNS].chLabel, ipLo+1, ipHi+1,fWLAng );
                                        }
                                }

                                if( fWLAng < 0.)
                                {
                                        // skip these lines
                                        wgfastream.ignore(INT_MAX,'\n');
                                        if( atmdat.lgChiantiPrint )
                                        {
                                                fprintf(ioQQQ,"WARNING: Skipping Transition %li to %li of %s.\n",ipLo+1,ipHi+1,dBaseSpecies[intNS].chLabel);
                                        }
                                        continue;
                                }

                        }
                        else
                        {
                                /* If Cloudy is using theoretical energies, then just flip the levels. */
                                if( ipHi < ipLo )
                                {
                                        long ipLoTemp = ipLo;
                                        long ipHiTemp = ipHi;
                                        ipHi = max( ipHiTemp, ipLoTemp );
                                        ipLo = min( ipHiTemp, ipLoTemp );
                                        fprintf( ioQQQ," WARNING: Swapped level indices for species %6s, indices %3li %3li with wavelength %e \n",
                                                        dBaseSpecies[intNS].chLabel, ipLoTemp, ipHiTemp,fWLAng );
                                }
                        }

                        /* If the given wavelength is negative, then theoretical enegies are being used.
                         * Take the difference in stored theoretical energies.
                         * It should equal the absolute value of the wavelength in the wgfa file. */
                        if( !atmdat.lgChiantiExp && fWLAng <= 0. )
                        {
                                //if( fWLAng < 0.)
                                        //fprintf( ioQQQ," WARNING: Negative wavelength for species %6s, indices %3li %3li \n", dBaseSpecies[intNS].chLabel, ipLo, ipHi);
                                fWLAng = (realnum)(1e8/
                                        (dBaseStates[intNS][ipHi].energy().WN() - dBaseStates[intNS][ipLo].energy().WN()));
                        }
                        //Skip from end of Wavelength to Einstein A and read in
                        wgfastream.seekg(line_nrg_to_eina,ios::cur);
                        wgfastream.get(trantemp,line_eina_col);
                        feinsteina = (realnum)atof(trantemp);
                        if( feinsteina < 1e-20 )
                        {
                                static bool notPrintedYet = true;
                                if( notPrintedYet && atmdat.lgChiantiPrint)
                                {
                                        fprintf( ioQQQ," WARNING: Radiative rate(s) equal to zero in %s.\n", chTraFilename );
                                        notPrintedYet = false;
                                }
                                wgfastream.ignore(INT_MAX,'\n');
                                continue;
                        }

                        fixit(); // may need to do something with these rates
                        //Read in the rest of the line and look for auto
                        wgfastream.getline(chLine,INT_MAX);
                        TransitionList::iterator tr = dBaseTrans[intNS].begin()+ipdBaseTrans[intNS][ipHi][ipLo];
                        if( nMatch("auto", chLine) )
                        {
                                if( (*tr).hasEmis() )
                                {
                                        (*tr).Emis().AutoIonizFrac() =
                                                feinsteina/((*tr).Emis().Aul() + feinsteina);
                                        ASSERT( (*tr).Emis().AutoIonizFrac() >= 0. );
                                        ASSERT( (*tr).Emis().AutoIonizFrac() <= 1. );
                                }
                                continue;
                        }

                        if( (*tr).hasEmis() )
                        {
                                fprintf(ioQQQ," PROBLEM duplicate transition found by atmdat_chianti in %s, "
                                                "wavelength=%f\n", chTraFilename,fWLAng);
                                fclose( ioProtCollData );
                                cdEXIT(EXIT_FAILURE);
                        }
                        (*tr).AddLine2Stack();
                        (*tr).Emis().Aul() = feinsteina;
                        (*tr).WLAng() = fWLAng;
                 
                        fenergyWN = (realnum)(1e+8/fWLAng);

                        // TODO::Check the wavelength in the file with the difference in energy levels

                        (*tr).EnergyWN() = fenergyWN;
                        (*tr).WLAng() = (realnum)(1e+8/fenergyWN/RefIndex(fenergyWN));
                        (*tr).Emis().gf() = (realnum)GetGF((*tr).Emis().Aul(), (*tr).EnergyWN(), (*(*tr).Hi()).g());

                        if(DEBUGSTATE)
                        {
                                fprintf( ioQQQ, "The lower level is %ld \n",ipLo);
                                fprintf( ioQQQ, "The upper level is %ld \n",ipHi);
                                fprintf( ioQQQ, "The Einstein A is %f \n",(*tr).Emis().Aul());
                                fprintf( ioQQQ, "The wavelength of the transition is %f \n",(*tr).WLAng() );
                        }
                }
        }
        else fprintf( ioQQQ, " The data file %s is corrupted .\n",chTraFilename);
        wgfastream.close();

        /* Malloc space for splines */
        AtmolCollSplines[intNS] = (CollSplinesArray***)MALLOC(nMolLevs *sizeof(CollSplinesArray**));
        for( long ipHi=0; ipHi<nMolLevs; ipHi++ )
        {
                AtmolCollSplines[intNS][ipHi] = (CollSplinesArray **)MALLOC((unsigned long)(nMolLevs)*sizeof(CollSplinesArray *));
                for( long ipLo=0; ipLo<nMolLevs; ipLo++ )
                {
                        AtmolCollSplines[intNS][ipHi][ipLo] = 
                                (CollSplinesArray *)MALLOC((unsigned long)(ipNCOLLIDER)*sizeof(CollSplinesArray ));

                        for( long k=0; k<ipNCOLLIDER; k++ )
                        {
                                /* initialize all spline variables */
                                AtmolCollSplines[intNS][ipHi][ipLo][k].collspline = NULL;
                                AtmolCollSplines[intNS][ipHi][ipLo][k].SplineSecDer = NULL;
                                AtmolCollSplines[intNS][ipHi][ipLo][k].nSplinePts = -1; 
                                AtmolCollSplines[intNS][ipHi][ipLo][k].intTranType = -1;
                                AtmolCollSplines[intNS][ipHi][ipLo][k].EnergyDiff = BIGDOUBLE;
                                AtmolCollSplines[intNS][ipHi][ipLo][k].ScalingParam = BIGDOUBLE;
                        }
                }
        }

        /************************************/
        /*** Read in the collisional data ***/
        /************************************/

        // ipCollider 0 is electrons, 1 is protons
        for( long ipCollider=0; ipCollider<=1; ipCollider++ )
        {
                char chFilename[FILENAME_PATH_LENGTH_2];

                if( ipCollider == ipELECTRON )
                {
                        strcpy( chFilename, chEleColFilename );
                }
                else if( ipCollider == ipPROTON )
                {
                        if( !lgProtonData )
                                break;
                        fprintf( ioQQQ," WARNING: Chianti proton collision data have different format than electron data files!\n" );
                        strcpy( chFilename, chProColFilename );
                }
                else
                        TotalInsanity();

                /*Dummy string used for convenience*/
                strcpy(chLine,"A");

                fstream splupsstream;
                open_data(splupsstream, chFilename,mode_r);

                //cs_eof_col is the length(+1) of the first column used for finding the end of file
                const int cs_eof_col = 4;
                //cs_index_col is the length(+1) of the indexes in the CS file
                const int cs_index_col = 4;
                //cs_trantype_col is the length(+1) of the transition type in the CS file
                const int cs_trantype_col = 4;
                //cs_values_col is the length(+1) of the other values in the CS file
                //including: GF, nrg diff, scaling parameter, and spline points
                const int cs_values_col = 11;
                //Determine the number of rows in the CS file
                if (splupsstream.is_open())
                {
                        int nj = 0;
                        //splupslines is -1 since the loop runs 1 extra time
                        long splupslines = -1;
                        char otemp[cs_eof_col];
                        while(nj != -1)
                        {
                                splupsstream.get(otemp,cs_eof_col);
                                splupsstream.ignore(INT_MAX,'\n');
                                nj = atoi(otemp);
                                splupslines++;
                        }
                        splupsstream.seekg(0,ios::beg);
        
                        for (int m = 0;m<splupslines;m++)
                        {
                                if( ipCollider == ipELECTRON )
                                {
                                        splupsstream.seekg(6,ios::cur);
                                }

                                /* level indices */
                                splupsstream.get(otemp,cs_index_col);
                                long ipLo = atoi(otemp)-1;
                                splupsstream.get(otemp,cs_index_col);
                                long ipHi = atoi(otemp)-1;

                                /* If either upper or lower index is -1 in the relational vector,
                                * skip that line in the splups file.
                                * Otherwise translate the level indices.*/
                                if( atmdat.lgChiantiExp )
                                {
                                        if( intExperIndex[ipLo] == - 1 || intExperIndex[ipHi] == -1 )
                                        {
                                                splupsstream.ignore(INT_MAX,'\n');
                                                continue;
                                        }
                                        else
                                        {
                                                ipLo = intNewIndex[intExperIndex[ipLo]];
                                                ipHi = intNewIndex[intExperIndex[ipHi]];
                                        }

                                }
                                else
                                {
                                        /* With level trimming on it is possible that there can be rows that
                                         * have to be skipped when using theoretical
                                         * since the levels no longer exist */
                                        if( intNewIndex[ipLo] == - 1 || intNewIndex[ipHi] == -1 )
                                        {
                                                splupsstream.ignore(INT_MAX,'\n');
                                                continue;
                                        }
                                        else
                                        {
                                                ipLo = intNewIndex[ipLo];
                                                ipHi = intNewIndex[ipHi];
                                        }
                                }

                                if( ipLo >= nMolLevs || ipHi >= nMolLevs )
                                {
                                        // skip these transitions
                                        splupsstream.ignore(INT_MAX,'\n');
                                        continue;
                                }

                                ASSERT( ipLo < nMolLevs );
                                ASSERT( ipHi < nMolLevs );
                                /*Transition Type*/
                                splupsstream.get(otemp,cs_trantype_col);
                                intTranType = atoi(otemp);
                                char qtemp[cs_values_col];
                                splupsstream.get(qtemp,cs_values_col);
                                /*Energy Difference*/
                                splupsstream.get(qtemp,cs_values_col);
                                fEnergyDiff = atof(qtemp);
                                /*Scaling Parameter*/
                                splupsstream.get(qtemp,cs_values_col);
                                fScalingParam = atof(qtemp);

                                ASSERT( ipLo < nMolLevs );
                                ASSERT( ipHi < nMolLevs );
                                ASSERT( AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline == NULL );
                                ASSERT( AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].SplineSecDer == NULL );

                                const int CHIANTI_SPLINE_MAX=9, CHIANTI_SPLINE_MIN=5;
                                STATIC_ASSERT(CHIANTI_SPLINE_MAX > CHIANTI_SPLINE_MIN);

                                /*We malloc the space here*/
                                AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline =
                                        (double *)MALLOC((unsigned long)(CHIANTI_SPLINE_MAX)*sizeof(double));
                                AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].SplineSecDer =
                                        (double *)MALLOC((unsigned long)(CHIANTI_SPLINE_MAX)*sizeof(double));

                                /* always read at least CHIANTI_SPLINE_MIN */
                                for( intsplinepts=0; intsplinepts<=CHIANTI_SPLINE_MAX; intsplinepts++ )
                                {
                                        //Look at the next character to see if it is the end of line.
                                        char p = splupsstream.peek();
                                        if( p == '\n' )
                                        {
                                                /* there should be 5 or 9 spline points.  If we got EOL,
                                                 * insist that we were trying to read the 6th or 10th. */
                                                if( (intsplinepts != CHIANTI_SPLINE_MAX) && (intsplinepts != CHIANTI_SPLINE_MIN) )
                                                {
                                                        static bool notPrintedYet = true;
                                                        if( notPrintedYet )
                                                        {
                                                                fprintf( ioQQQ, " WARNING: Wrong number of spline points in %s.\n", chFilename);
                                                                notPrintedYet = false;
                                                        }
                                                        for( long j=intsplinepts-1; j<CHIANTI_SPLINE_MAX; j++ )
                                                                AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline[j] = 0.;
                                                }
                                                else
                                                        break;
                                        }
                                        else
                                        {
                                                if( intsplinepts >= CHIANTI_SPLINE_MAX )
                                                {
                                                        fprintf( ioQQQ, " WARNING: More spline points than expected in %s, indices %3li %3li.  Ignoring extras.\n", chFilename, ipHi, ipLo );
                                                        break;
                                                }
                                                ASSERT( intsplinepts < CHIANTI_SPLINE_MAX );
                                                double temp;
                                                //Store a single spline point then look for more
                                                splupsstream.get(qtemp,cs_values_col);
                                                temp = atof(qtemp);
                                                if(temp < 0)
                                                {
                                                        temp = 0.;
                                                }
                                                AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline[intsplinepts] = temp;
                                        }

                                }

                                ASSERT( (intsplinepts == CHIANTI_SPLINE_MAX) || (intsplinepts == CHIANTI_SPLINE_MIN) );

                                /*The zeroth element contains the number of spline points*/
                                AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].nSplinePts = intsplinepts;
                                /*Transition type*/
                                AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].intTranType = intTranType;
                                /*Energy difference between two levels in Rydbergs*/
                                AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].EnergyDiff = fEnergyDiff;
                                /*Scaling parameter C*/
                                AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].ScalingParam = fScalingParam;

                                /*Once the spline points have been filled,fill the second derivatives structure*/
                                /*Creating spline points array*/
                                xs = (double *)MALLOC((unsigned long)intsplinepts*sizeof(double));
                                spl = (double *)MALLOC((unsigned long)intsplinepts*sizeof(double));
                                spl2 = (double *)MALLOC((unsigned long)intsplinepts*sizeof(double));

                                // should be able to just loop to intsplinepts -- ASSERT above protects
                                if(intsplinepts == CHIANTI_SPLINE_MIN)
                                {
                                        for(intxs=0;intxs<CHIANTI_SPLINE_MIN;intxs++)
                                        {
                                                xs[intxs] = 0.25*intxs;
                                                spl[intxs] = AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline[intxs];
                                        }
                                }
                                else if(intsplinepts == CHIANTI_SPLINE_MAX)
                                {
                                        for(intxs=0;intxs<CHIANTI_SPLINE_MAX;intxs++)
                                        {
                                                xs[intxs] = 0.125*intxs;
                                                spl[intxs] = AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline[intxs];
                                        }
                                }
                                else
                                        TotalInsanity();

                                spline(xs, spl,intsplinepts,2e31,2e31,spl2);

                                /*Filling the second derivative structure*/
                                for( long i=0; i<intsplinepts; i++)
                                {
                                        AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].SplineSecDer[i] = spl2[i];
                                }

                                free( xs );
                                free( spl );
                                free( spl2 );
                                splupsstream.ignore(INT_MAX,'\n');
                        }
                        splupsstream.close();
                }
        }

        // free some memory
        free( dBaseStatesEnergy );
        free( intNewIndex );
        if( atmdat.lgChiantiExp)
        {
                free( intExperIndex );
        }

        // close open file handles
        fclose( ioElecCollData );
        if( lgProtonData )
                fclose( ioProtCollData );

        return;
}

---