/home66/gary/public_html/cloudy/c08_branch/source/nemala.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 "taulines.h"
#include "trace.h"
#include "string.h"
#include "thirdparty.h"
#include "dense.h"
#include "atmdat.h"
/*File nemala.cpp was developed by Humeshkar B Nemala as a part of his thesis work during the Summer of 2007*/
/* Initially the code has been developed to read in energy levels,radiative and
 * collisional data from the CHIANTI and LEIDEN databases. The idea is to extend it to more databases.
 * In the case of the Leiden database there is a single .dat file which has the energy levels information,
 * radiative and collisional data, with the data corresponding to each collider coming one after the other.
 * In the case of CHIANTI, the energy levels data, radiative data and collision data are present in seperate files.
 * While LEIDEN gives collisional rate coefficients, CHIANTI gives collisional strengths.
 * In the case of CHIANTI only two colliders are used:electrons and protons. They appear as separate files.
 * The electron collision strengths files are always expected to be there. A flag is set and data processed 
 * if the file on proton collision strengths is available.*/

/* There is an initialization file called species.ini which tells Cloudy what kind of data is to be used */
/* Structures are created separately to hold the transition data,radiative and collisional data */
/* The collisional structures are different for different databases depending upon whether */
/* collisional strengths or collisional rate coefficients are used.Finally a superstructure is constructed to hold */
/* the total collisional rate obtained by considering all the colliders */
/* The colliders considered are electron,proton,Atomic Hydrogen,He,He+,He++,Ortho Molecular Hydrogen,Para Molecular Hydrogen and Molecular Hydrogen */
void database_readin(void);
void states_popfill( void);
void states_nelemfill(void);
void database_prep(int);
void emislines_fillredis(void);
STATIC void states_propprint(void);
STATIC int getAtNo(char []);

#define DEBUGSTATE false

emission *AddLine2Stack( realnum Aul, transition *trans )
{

        DEBUG_ENTRY( "AddLine2Stack()" );

        ASSERT(linesAdded2 < MAX_NUM_LINES);
        atmolEmis[linesAdded2].Aul = Aul; 
        atmolEmis[linesAdded2].tran = trans; 
        linesAdded2++;
        ASSERT(linesAdded2 <= MAX_NUM_LINES);
        return &atmolEmis[linesAdded2-1];
}

void Nemala_Start( void )
{
        int i,j,los,intNoSp;

        FILE *atmolDATA;
        char *chToken,*chAtNo,*chIonStg;

        char chLine[FILENAME_PATH_LENGTH_2],
                chDLine[FILENAME_PATH_LENGTH_2];

        static int nCalled = 0;

        DEBUG_ENTRY( "Nemala_Start()" );

        linesAdded2 = 0;

        /* only do this once. */
        if( nCalled > 0 )
        {
                return;
        }

        /* this is first call, increment the nCalled counterso never do this again */
        ++nCalled;

        // Reading in the species.ini file: Humeshkar B Nemala
        if( trace.lgTrace )
                fprintf( ioQQQ," Nemala_Start opening species.ini:");

        atmolDATA = open_data( "species.ini", "r" );

        if( read_whole_line( chLine , (int)sizeof(chLine) , atmolDATA ) == NULL )
        {
                fprintf( ioQQQ, " Nemala_Start could not read first line of species.ini.\n");
                cdEXIT(EXIT_FAILURE);
        }

        /* count how many lines are in the file, ignoring all lines
         * starting with '#':This would give the number of molecules */
        nSpecies = 0;
        while( read_whole_line( chLine , (int)sizeof(chLine) , atmolDATA ) != NULL )
        {
                /* we want to count the lines that do not start with %
                 * since these contain data */
                if( (chLine[0] != '#')&&(chLine[0] != '\n')&&(chLine[0] != ' '))
                        ++nSpecies;
        }
        if(DEBUGSTATE)
                printf("The number of species is %li \n",nSpecies);

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

        /*Once we know the number of species we initialize a pointer to a one 
         * dimensional array which contains information about the DB Type */
        lgSpeciesMolecule = (bool *)MALLOC( (unsigned long)(nSpecies)*sizeof(bool));

        /*Initialization of the Species Structure*/
        Species = (species *)MALLOC( (unsigned long)nSpecies*sizeof(species));

        /*Initialization of the collisional rates array structure*/
        /*Assume that there are 9 colliders*/
        CollRatesArray = (double****)MALLOC((unsigned long)nSpecies * sizeof(double***));
        AtmolCollRateCoeff = (CollRateCoeffArray **)MALLOC((unsigned long)nSpecies * sizeof(CollRateCoeffArray*));
        AtmolCollSplines = (CollSplinesArray****)MALLOC((unsigned long)nSpecies *sizeof(CollSplinesArray***));

        /*Mallocing here takes care of the number of colliders*/
        for( i=0; i<nSpecies; i++ )
        {
                AtmolCollRateCoeff[i] = (CollRateCoeffArray *)MALLOC(NUM_COLLIDERS * sizeof(CollRateCoeffArray));
                //AtmolCollSplines[i] = (CollSplinesArray***)MALLOC(NUM_COLLIDERS *sizeof(CollSplinesArray**));
                CollRatesArray[i] = (double***)MALLOC(NUM_COLLIDERS * sizeof(double**));
        }
        /*Initializing all the elements of the matrix to 0*/
        /*For each species and collider*/
        for( i=0; i<nSpecies; i++ )
        {
                for( j=0; j<NUM_COLLIDERS; j++ )
                {
                        AtmolCollRateCoeff[i][j].ntemps = 0;
                        AtmolCollRateCoeff[i][j].temps = NULL;
                        AtmolCollRateCoeff[i][j].collrates = NULL;
                }
        }

        /*****************************************/
        i=0;
        while( read_whole_line( chLine , (int)sizeof(chLine) , atmolDATA ) != NULL )
        {       
                caps( chLine );
                if ((chLine[0]!='#') && (chLine[0]!='\n')&&(chLine[0]!='\t')&&(chLine[0]!='\r'))
                {       
                        ASSERT(i<nSpecies);
                        if(DEBUGSTATE)
                                printf("%s",chLine);
                        strcpy(chDLine, chLine);
                        chToken = strtok(chDLine," ");
                        los = (int)strlen(chToken);
                        if(DEBUGSTATE)
                                printf("The length of the string is %d\n",los);
                        Species[i].chptrSpName = (char *)MALLOC((unsigned long)(los+1) * sizeof(char));
                        strcpy(Species[i].chptrSpName,chToken);
                        Species[i].intAtNo = 0;
                        Species[i].intIonStage = 0;
                        if(DEBUGSTATE)
                                printf("The name of the species is%s \n",Species[i].chptrSpName);
                        if( nMatch("LEID",chLine) )
                        {
                                lgSpeciesMolecule[i] = true;
                        }
                        else if( nMatch("CHIA",chLine))
                        {
                                lgSpeciesMolecule[i] = false;
                                chAtNo = strtok(chToken,"_");
                                Species[i].intAtNo = getAtNo(chAtNo);
                                chIonStg = strtok(NULL,"_");
                                Species[i].intIonStage = atoi(chIonStg)-1;
                        }
                        else
                                TotalInsanity();

                        if(DEBUGSTATE)
                                printf("The species is molecular? %c \n",TorF(lgSpeciesMolecule[i]) );
                        ++i;
                }
        }

        database_readin();

        /*Setting the population of states to 0*/
        states_popfill();
        /*Setting the values of nelem arbitrarily to 1*/
        states_nelemfill();

        /*Setting nelem of the states to an arbitrary value*/
        /*Loop over species*/
        for( intNoSp=0; intNoSp<nSpecies; intNoSp++ )
        {
                database_prep(intNoSp);
        }

        /*Looping over the emission lines and filling in the redistribution function value*/
        emislines_fillredis();

        /*To print the states*/
        if(DEBUGSTATE)
                states_propprint();
        return;
}

/*Filling the redistribution function value*/
void emislines_fillredis(void)
{
        DEBUG_ENTRY("emislines_fillredis()");
        for( long i=0; i<linesAdded2; i++ )
        {
                atmolEmis[i].iRedisFun = ipPRD;
        }
        return;
}

/*This function fills the nelem value arbitrarily*/
void states_nelemfill(void)
{
        DEBUG_ENTRY( "states_nelemfill()" );

        for( long i=0; i<nSpecies; i++ )
        {
                for( long j=0; j<Species[i].numLevels_max; j++ )
                {
                        if( lgSpeciesMolecule[i] )
                        {
                                fixit(); /* this is needed because right now, all lines test on 
                                                  * dense.xIonDense[nelem][IonStg] */
                                atmolStates[i][j].nelem = 1;
                                atmolStates[i][j].IonStg = 1;
                        }
                        else
                        {
                                atmolStates[i][j].nelem = Species[i].intAtNo + 1;
                                atmolStates[i][j].IonStg = Species[i].intIonStage;
                        }
                }
        }
        return;
}

/*This function prints the various properties of states*/
STATIC void states_propprint(void)
{
        DEBUG_ENTRY( "states_propprint()" );

        for( long i=0; i<nSpecies; i++ )
        {
                printf("The species is %s \n",Species[i].chptrSpName);
                printf("The data output is in the following format \n");
                printf("Label Energy St.wt Pop Lifetime\n");

                for( long j=0; j<Species[i].numLevels_max; j++ )
                {
                        printf("This is the %ld state \n",j);
                        printf("%s %f %f %f %e \n",atmolStates[i][j].chLabel,
                                atmolStates[i][j].energy,
                                atmolStates[i][j].g,
                                atmolStates[i][j].Pop,
                                atmolStates[i][j].lifetime);
                }
        }
        return;
}


void database_prep(int intSpIndex)
{
        realnum fsumAs;
        int ipHi,ipLo;

        DEBUG_ENTRY( "database_prep()" );

        /*Get the lifetimes*/
        atmolStates[intSpIndex][0].lifetime= BIGFLOAT;
        for( ipHi=1; ipHi < Species[intSpIndex].numLevels_max; ipHi++ )
        {
                fsumAs = SMALLFLOAT;
                for( ipLo=0; ipLo<ipHi; ipLo++ )
                {
                        if(atmolTrans[intSpIndex][ipHi][ipLo].Emis!=NULL)
                                fsumAs = fsumAs + atmolTrans[intSpIndex][ipHi][ipLo].Emis->Aul;
                }
                atmolStates[intSpIndex][ipHi].lifetime = 1./fsumAs;
        }
        return;
}

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

        long intNS;

        atmolStates = (quantumState **)MALLOC((unsigned long)nSpecies * sizeof(quantumState*));
        atmolTrans = (transition ***)MALLOC((unsigned long)nSpecies * sizeof(transition**));

        /* nSpecies should be global*/
        /*The data array containing the species names should also be global*/
        for( intNS = 0; intNS < nSpecies; intNS++ )
        {               
                if ( lgSpeciesMolecule[intNS] )
                {       
                        atmdat_lamda_readin( intNS );
                }
                else
                {       
                        atmdat_Chianti_readin( intNS );
                }
        }
        return;
}

STATIC int getAtNo(char *p)
{
        DEBUG_ENTRY("getAtNo()");
        if(strcmp(p,"H")==0)
        {
                return(ipHYDROGEN);
        }
        else if(strcmp(p,"HE")== 0)
        {
                return(ipHELIUM);
        }
        else if(strcmp(p,"C")== 0)
        {
                return(ipCARBON);
        }
        else if(strcmp(p,"N")== 0)
        {
                return(ipNITROGEN);
        }
        else if(strcmp(p,"O")== 0)
        {
                return(ipOXYGEN);
        }
        else if(strcmp(p,"NE")== 0)
        {
                return(ipNEON);
        }
        else if(strcmp(p,"NA")== 0)
        {
                return(ipSODIUM);
        }
        else if(strcmp(p,"MG")== 0)
        {
                return(ipMAGNESIUM);
        }
        else if(strcmp(p,"AL")== 0)
        {
                return(ipALUMINIUM);
        }
        else if(strcmp(p,"SI")== 0)
        {
                return(ipSILICON);
        }
        else if(strcmp(p,"P")== 0)
        {
                return(ipPHOSPHORUS);
        }
        else if(strcmp(p,"S")== 0)
        {
                return(ipSULPHUR);
        }
        else if(strcmp(p,"CL")== 0)
        {
                return(ipCHLORINE);
        }
        else if(strcmp(p,"AR")== 0)
        {
                return(ipARGON );
        }
        else if(strcmp(p,"K")== 0)
        {
                return(ipPOTASSIUM);
        }
        else if(strcmp(p,"CA")== 0)
        {
                return(ipCALCIUM);
        }
        else if(strcmp(p,"SC")== 0)
        {
                return(ipSCANDIUM);
        }
        else if(strcmp(p,"TI")== 0)
        {
                return(ipTITANIUM);
        }
        else if(strcmp(p,"V")== 0)
        {
                return(ipVANADIUM);
        }
        else if(strcmp(p,"CR")== 0)
        {
                return(ipCHROMIUM );
        }
        else if(strcmp(p,"MN")== 0)
        {
                return(ipMANGANESE);
        }
        else if(strcmp(p,"FE")== 0)
        {
                return(ipIRON);
        }
        else if(strcmp(p,"CO")== 0)
        {
                return(ipCOBALT);
        }
        else if(strcmp(p,"NI")== 0)
        {
                return(ipNICKEL);
        }
        else if(strcmp(p,"CU")== 0)
        {
                return(ipCOPPER);
        }
        else if(strcmp(p,"ZN")== 0)
        {
                return(ipZINC);
        }
        else
        {
                TotalInsanity();
        }
}

---