/home66/gary/public_html/cloudy/c08_branch/source/atmdat_lamda.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 */
#include "cddefines.h"
#include "lines_service.h"
#include "physconst.h"
#include "taulines.h"
#include "trace.h"
#include "string.h"
#include "thirdparty.h"

#define DEBUGSTATE false

/*Separate Routine to read in the molecules*/
void atmdat_lamda_readin( long intNS )
{
        DEBUG_ENTRY( "atmdat_lamda_readin()" );

        int ngMolLevs = -10000,intCollIndex = -10000,intLColliderIndex = -10000;
        /* type is set to 0 for non chianti and 1 for chianti*/
        long int i,j,nMolLevs,intlnct,intrtct,intgrtct,intCollPart,
                intDCollPart,intCollTran,nCollTrans,intCollTemp,intcollindex,
                ipLo,ipHi;
        /*intrtct refers to radiative transitions count*/
        FILE *atmolLevDATA;
        realnum  fstatwt,fenergyK,fenergyWN,fenergy,feinsteina;

        char chLine[FILENAME_PATH_LENGTH_2] ,
                chEFilename[FILENAME_PATH_LENGTH_2],
                *chColltemp,*chCollName;

        ASSERT( lgSpeciesMolecule[intNS] );

        /*Load the species name in the Species array structure*/
        if(DEBUGSTATE)
                printf("The name of the %li species is %s \n",intNS+1,Species[intNS].chptrSpName);
        strcpy( chEFilename , Species[intNS].chptrSpName );                     
        strcat( chEFilename , ".dat");
        uncaps( chEFilename );

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

        atmolLevDATA = open_data( chEFilename, "r" );

        nMolLevs = 0;
        ngMolLevs = 0;
        intrtct = 0;
        intgrtct = 0;
        intlnct = 0;
        ipHi = 0;
        ipLo = 0;
        j = 0;
        while(intlnct < 3)
        {
                intlnct++;
                if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
                {
                        fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
                        cdEXIT(EXIT_FAILURE);
                }
        }
        /*Extracting out the molecular weight*/
        if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
        {
                fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
                cdEXIT(EXIT_FAILURE);
        }
        Species[intNS].fmolweight = (realnum)atof(chLine);

        /*Discard this line*/
        if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
        {
                fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
                cdEXIT(EXIT_FAILURE);
        }
        /*Reading in the number of energy levels*/
        if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
        {
                fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
                cdEXIT(EXIT_FAILURE);
        }
        ngMolLevs = atoi(chLine);

        if(ngMolLevs <= 0)
        {
                fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
                cdEXIT(EXIT_FAILURE);
        }
        Species[intNS].numLevels_max = ngMolLevs;
        Species[intNS].numLevels_local = ngMolLevs;
        if(DEBUGSTATE)
        {
                printf("The number of energy levels is %li \n",Species[intNS].numLevels_max);
        }
        /*the above refers to the number of energy levels specified in the data file*/
        /*malloc the States array*/
        atmolStates[intNS] = (quantumState *)MALLOC((size_t)(ngMolLevs)*sizeof(quantumState));
        /*malloc the Transition array*/
        atmolTrans[intNS] = (transition **)MALLOC(sizeof(transition*)*(unsigned long)ngMolLevs);
        atmolTrans[intNS][0] = NULL;
        for( ipHi = 1; ipHi < ngMolLevs; ipHi++)
        {
                atmolTrans[intNS][ipHi] = (transition *)MALLOC(sizeof(transition)*(unsigned long)ipHi);
                for(ipLo = 0; ipLo < ipHi; ipLo++)
                {
                        atmolTrans[intNS][ipHi][ipLo].Lo = &atmolStates[intNS][ipLo];
                        atmolTrans[intNS][ipHi][ipLo].Hi = &atmolStates[intNS][ipHi];
                        atmolTrans[intNS][ipHi][ipLo].Emis = &DummyEmis;
                }
        }

        for( intcollindex = 0; intcollindex <NUM_COLLIDERS; intcollindex++ )
        {
                CollRatesArray[intNS][intcollindex] = NULL;
        }


        if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
        {
                fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
                cdEXIT(EXIT_FAILURE);
        }

        for( nMolLevs=0; nMolLevs<ngMolLevs; nMolLevs++)
        {       
                if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
                {
                        fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
                        cdEXIT(EXIT_FAILURE);
                }
                /*information needed for label*/
                strcpy(atmolStates[intNS][nMolLevs].chLabel, "    ");
                strncpy(atmolStates[intNS][nMolLevs].chLabel,Species[intNS].chptrSpName, 4);
                atmolStates[intNS][nMolLevs].chLabel[4] = '\0';

                long i = 1;
                bool lgEOL;
                long index;

                index = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
                fenergy = (realnum)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
                fstatwt = (realnum)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );

                ASSERT( index == nMolLevs + 1 );

                atmolStates[intNS][nMolLevs].energy = fenergy;
                atmolStates[intNS][nMolLevs].g = fstatwt;

                if (nMolLevs > 0)
                {
                        if (atmolStates[intNS][nMolLevs].energy < atmolStates[intNS][nMolLevs-1].energy)
                        {
                                fprintf( ioQQQ, " The energy levels are not in order");
                                cdEXIT(EXIT_FAILURE);
                        }
                }
                if(DEBUGSTATE)
                {
                        printf("The converted energy is %f \n",atmolStates[intNS][nMolLevs].energy);
                        printf("The value of g is %f \n",atmolStates[intNS][nMolLevs].g);
                }
        }

        /* fill in all transition energies, can later overwrite for specific radiative transitions */
        for( i=1; i<ngMolLevs; i++ )
        {
                for( j=0; j<i; j++ )
                {
                        fenergyWN = atmolStates[intNS][i].energy - atmolStates[intNS][j].energy;
                        fenergyK = (realnum)(fenergyWN*T1CM);

                        atmolTrans[intNS][i][j].EnergyK = fenergyK;
                        atmolTrans[intNS][i][j].EnergyWN = fenergyWN;
                        atmolTrans[intNS][i][j].EnergyErg = (realnum)ERG1CM *fenergyWN;

                        /* protect against division by zero. */
                        atmolTrans[intNS][i][j].WLAng = (realnum)(1e+8f/MAX2(1E-5f,fenergyWN)/
                                RefIndex(fenergyWN));
                }
        }

        if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
        {
                fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
                cdEXIT(EXIT_FAILURE);
        }
        if(chLine[0]!='!')
        {
                fprintf( ioQQQ, " The number of energy levels in file %s is not correct.\n",chEFilename);
                cdEXIT(EXIT_FAILURE);
        }
        if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
        {
                fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
                cdEXIT(EXIT_FAILURE);
        }
        intgrtct = atoi(chLine);
        /*The above gives the number of radiative transitions*/
        if(intgrtct <= 0)
        {
                fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
                cdEXIT(EXIT_FAILURE);
        }
        if(DEBUGSTATE)
        {
                printf("The number of radiative transitions is %li \n",intgrtct);
        }
        if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
        {
                fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
                cdEXIT(EXIT_FAILURE);
        }

        for( intrtct=0; intrtct<intgrtct; intrtct++)
        {
                if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
                {
                        fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
                        cdEXIT(EXIT_FAILURE);
                }

                long i = 1;
                bool lgEOL;
                long index;

                index = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
                ipHi = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ) - 1;
                ipLo = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ) - 1;
                feinsteina = (realnum)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
                /* don't need the energy in GHz, so throw it away. */
                FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
                fenergyK = (realnum)(((atmolStates[intNS][ipHi].energy) -(atmolStates[intNS][ipLo].energy))*T1CM);
                ASSERT( index == intrtct + 1 );

                atmolTrans[intNS][ipHi][ipLo].Emis = AddLine2Stack(feinsteina, &atmolTrans[intNS][ipHi][ipLo]);
                atmolTrans[intNS][ipHi][ipLo].EnergyK = fenergyK;
                fenergyWN = (realnum)((fenergyK)/T1CM);
                atmolTrans[intNS][ipHi][ipLo].EnergyWN = fenergyWN;
                atmolTrans[intNS][ipHi][ipLo].EnergyErg = (realnum)ERG1CM *fenergyWN;
                atmolTrans[intNS][ipHi][ipLo].WLAng = (realnum)(1e+8/fenergyWN/RefIndex(fenergyWN));

                ASSERT( !isnan( atmolTrans[intNS][ipHi][ipLo].EnergyK ) );

                if(DEBUGSTATE)
                {
                        printf("The upper level is %ld \n",ipHi+1);
                        printf("The lower level is %ld \n",ipLo+1);
                        printf("The Einstein A  is %E \n",atmolTrans[intNS][ipHi][ipLo].Emis->Aul);
                        printf("The Energy of the transition is %E \n",atmolTrans[intNS][ipHi][ipLo].EnergyK);
                }
        }

        if(read_whole_line( chLine , (int)sizeof(chLine)  , atmolLevDATA ) == NULL )
        {
                fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
                cdEXIT(EXIT_FAILURE);
        }

        if(chLine[0]!='!')
        {
                fprintf( ioQQQ, " The number of radiative transitions in file %s is not correct.\n",chEFilename);
                cdEXIT(EXIT_FAILURE);
        }

        /*Getting the number of collisional partners*/
        if(read_whole_line( chLine , (int)sizeof(chLine)  , atmolLevDATA ) == NULL )
        {
                fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
                cdEXIT(EXIT_FAILURE);
        }
        else
        {
                intCollPart = atoi(chLine);
        }
        /*Checking the number of collisional partners does not exceed 9*/
        if(intCollPart > NUM_COLLIDERS-1)
        {
                fprintf( ioQQQ, " The number of colliders is greater than what is expected in file %s.\n", chEFilename );
                cdEXIT(EXIT_FAILURE);
        }
        /*Creating the duplicate of the number of collisional partners which is reduced*/
        intDCollPart = intCollPart;
        while(intDCollPart > 0)
        {
                if(read_whole_line( chLine , (int)sizeof(chLine)  , atmolLevDATA ) == NULL )
                {
                        fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
                        cdEXIT(EXIT_FAILURE);
                }

                if(read_whole_line( chLine , (int)sizeof(chLine)  , atmolLevDATA ) == NULL )
                {
                        fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
                        cdEXIT(EXIT_FAILURE);
                }
                /*Extract out the name of the collider*/
                /*The following are the rules expected in the datafiles to extract the names of the collider*/
                /*The line which displays the species and the collider starts with a number*/
                /*This refers to the collider in the Leiden database*/
                /*In the Leiden database 1 referes to H2,2 to para-H2,3 to ortho-H2
                4 to electrons,5 to H and 6 to He*/
                chCollName = strtok(chLine," ");
                /*Leiden Collider Index*/
                intLColliderIndex = atoi(chCollName);
                /*In Cloudy,We assign the following indices for the colliders*/
                /*electron=0,proton=1,atomic hydrogen=2,He=3,He+=4,He++=5,oH2=6,pH2=7,H2=8*/

                if(intLColliderIndex == 1)
                {
                        intCollIndex = 8;
                }
                else if(intLColliderIndex == 2)
                {
                        intCollIndex = 7;
                }
                else if(intLColliderIndex == 3)
                {
                        intCollIndex = 6;
                }
                else if(intLColliderIndex == 4)
                {
                        intCollIndex = 0;
                }
                else if(intLColliderIndex == 5)
                {
                        intCollIndex = 2;
                }
                else if(intLColliderIndex == 6)
                {
                        intCollIndex = 3;
                }
                else
                {
                        TotalInsanity();
                }
                /*This is where we allocate memory if the collider exists*/
                /*Needed to take care of the he collisions*/
                CollRatesArray[intNS][intCollIndex] = (double**)MALLOC((unsigned long)(ngMolLevs) * sizeof(double*));
                for( ipHi = 0; ipHi<ngMolLevs; ipHi++ )
                {
                        CollRatesArray[intNS][intCollIndex][ipHi] = (double*)MALLOC((unsigned long)(ngMolLevs) * sizeof(double));
                        for( ipLo = 0; ipLo<ngMolLevs; ipLo++ )
                        {
                                CollRatesArray[intNS][intCollIndex][ipHi][ipLo] = 0.0;
                        }
                }
                if(read_whole_line( chLine , (int)sizeof(chLine)  , atmolLevDATA ) == NULL )
                {
                        fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
                        cdEXIT(EXIT_FAILURE);
                }
                if(read_whole_line( chLine , (int)sizeof(chLine)  , atmolLevDATA ) == NULL )
                {
                        fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
                        cdEXIT(EXIT_FAILURE);
                }
                /*Number of coll trans*/
                intCollTran = atoi(chLine);

                if(read_whole_line( chLine , (int)sizeof(chLine)  , atmolLevDATA ) == NULL )
                {
                        fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
                        cdEXIT(EXIT_FAILURE);
                }
                if(read_whole_line( chLine , (int)sizeof(chLine)  , atmolLevDATA ) == NULL )
                {
                        fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
                        cdEXIT(EXIT_FAILURE);
                }
                /*Number of coll temps*/
                intCollTemp = atoi(chLine);
                /*Storing the number of collisional temperatures*/
                AtmolCollRateCoeff[intNS][intCollIndex].ntemps = intCollTemp;
                /*Mallocing*/
                AtmolCollRateCoeff[intNS][intCollIndex].temps = 
                        (double *)MALLOC((unsigned long)intCollTemp*sizeof(double));
                AtmolCollRateCoeff[intNS][intCollIndex].collrates = 
                        (double***)MALLOC((unsigned long)(ngMolLevs)*sizeof(double**));
                for( i=0; i<ngMolLevs; i++ )
                {
                        AtmolCollRateCoeff[intNS][intCollIndex].collrates[i] = 
                                (double **)MALLOC((unsigned long)(ngMolLevs)*sizeof(double*));
                        for( j=0; j<ngMolLevs; j++ )
                        {
                                AtmolCollRateCoeff[intNS][intCollIndex].collrates[i][j] = 
                                        (double *)MALLOC((unsigned long)(intCollTemp)*sizeof(double));
                        }
                }

                /*Discard the header line*/
                if(read_whole_line( chLine , (int)sizeof(chLine)  , atmolLevDATA ) == NULL )
                {
                        fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
                        cdEXIT(EXIT_FAILURE);
                }
                /*Getting the collisional Temperatures*/
                if(read_whole_line( chLine , (int)sizeof(chLine)  , atmolLevDATA ) == NULL )
                {
                        fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
                        cdEXIT(EXIT_FAILURE);
                }
                /*Filling the collisional temps array*/
                chColltemp =strtok(chLine," ");
                AtmolCollRateCoeff[intNS][intCollIndex].temps[0] =(realnum) atof(chColltemp);
                for( i=1; i<intCollTemp; i++ )
                {
                        chColltemp =strtok(NULL," ");
                        AtmolCollRateCoeff[intNS][intCollIndex].temps[i] = atof(chColltemp);
                }
                /*Discard the header line*/
                if(read_whole_line( chLine , (int)sizeof(chLine)  , atmolLevDATA ) == NULL )
                {
                        fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
                        cdEXIT(EXIT_FAILURE);
                }
                /*Getting all the collisional transitions data*/
                for( nCollTrans=0; nCollTrans<intCollTran; nCollTrans++ )
                {
                        if(read_whole_line( chLine , (int)sizeof(chLine)  , atmolLevDATA ) == NULL )
                        {
                                fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
                                cdEXIT(EXIT_FAILURE);
                        }

                        long i = 1;
                        bool lgEOL;
                        long index;

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

                        /* Indices between the very highest levels seem to be reversed */
                        if( ipHi < ipLo )
                        {
                                ASSERT( ipLo == ngMolLevs - 1);
                                long temp = ipHi;
                                ipHi = ipLo;
                                ipLo = temp;
                        }

                        for( long j=0; j<intCollTemp; j++ )
                        {
                                AtmolCollRateCoeff[intNS][intCollIndex].collrates[ipHi][ipLo][j] = 
                                        (double)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
                        }

                        ASSERT( index == nCollTrans + 1 );

                        if(DEBUGSTATE)
                        {
                                printf("The values of up and lo are %ld & %ld \n",ipHi,ipLo);
                                printf("The collision rates are");
                                for(i=0;i<intCollTemp;i++)
                                {
                                        printf("\n %e",AtmolCollRateCoeff[intNS][intCollIndex].collrates[ipHi][ipLo][i]);
                                }
                                printf("\n");
                        }
                }

                intDCollPart = intDCollPart -1;
        }

        return;
}

---