/* 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 "taulines.h"
#include "trace.h"
#include "string.h"
#include "input.h"
#include "thirdparty.h"
#include "dense.h"
#include "atmdat.h"
#include "mole.h"
#include "elementnames.h"
#include "version.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 */
STATIC void states_popfill(void);
STATIC void states_nelemfill(void);
STATIC void database_prep(int);
STATIC void set_fractionation( species *sp );
STATIC void states_propprint(void);
/*SpeciesJunk set all elements of species struc to dangerous values */
STATIC void SpeciesJunk( species *sp );

#define DEBUGSTATE false
void database_readin( void )
{
	int i,intNoSp;

	FILE *ioMASTERLIST, *ioVERSION;

	char *chToken;

	char chLine[FILENAME_PATH_LENGTH_2],
		chDLine[FILENAME_PATH_LENGTH_2],
		chPath[FILENAME_PATH_LENGTH_2] = "";

	const int MAX_NUM_SPECIES = 1000;

	char chLabels[MAX_NUM_SPECIES][CHARS_SPECIES];
	char chLabelsOrig[MAX_NUM_SPECIES][CHARS_SPECIES];
	char chPaths[MAX_NUM_SPECIES][FILENAME_PATH_LENGTH_2];

	static int nCalled = 0;
	long nSpeciesLAMDA, nSpeciesSTOUT, nSpeciesCHIANTI;

	DEBUG_ENTRY( "database_readin()" );

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

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

	// read masterlists, count number of species
	nSpecies = 0;

	////////////////////////////////////
	//  
	// Read LAMDA masterlist 
	//
	//////////////////////////////////

	/* count how many lines are in the file, ignoring all lines
	 * starting with '#':This would give the number of molecules */
	nSpeciesLAMDA = 0;

	if( atmdat.lgLamdaOn )
	{
		long numModelsNotUsed = 0;
		strcpy( chPath, "lamda" );
		strcat( chPath, input.chDelimiter );
		strcat( chPath, "masterlist" );
		ioMASTERLIST = open_data( chPath, "r" );

		if( read_whole_line( chLine , (int)sizeof(chLine) , ioMASTERLIST ) == NULL )
		{
			fprintf( ioQQQ, " database_readin could not read first line of LAMDA masterlist.\n");
			cdEXIT(EXIT_FAILURE);
		}

		do
		{	
			if ((chLine[0]!='#') && (chLine[0]!='\n')&&(chLine[0]!='\t')&&(chLine[0]!='\r'))
			{	
				strcpy(chDLine, chLine);
				chToken = strtok(chDLine," \t\n");
				if( findspecies( chToken ) != null_mole  || 
					( chToken[1]=='-' && findspecies( chToken+2 ) != null_mole ) )
				{
					ASSERT( nSpecies + 1 <= MAX_NUM_SPECIES );
					ASSERT( nSpeciesLAMDA + 1 <= MAX_NUM_SPECIES );
					ASSERT( strlen(chToken) < CHARS_SPECIES );
					strcpy( chLabels[nSpecies], chToken );
					chLabels[nSpecies][CHARS_SPECIES-1] = '\0';

					// path is, for example, LAMDA/no.dat
					strcpy( chPaths[nSpecies], "lamda" );
					strcat( chPaths[nSpecies], input.chDelimiter );
					chToken = strtok( NULL," \t\n" );
					strcat( chPaths[nSpecies], chToken );
					++nSpecies;
					++nSpeciesLAMDA;
				}
				else
					++numModelsNotUsed;
			}
		}
		while( read_whole_line( chLine , (int)sizeof(chLine) , ioMASTERLIST ) != NULL );

		/* \todo 1 - save this and stuff as note since not really a "PROBLEM" but worth reporting */
		//if( !t_version::Inst().lgRelease && numModelsNotUsed > 0 )
		//	fprintf( ioQQQ, "\n PROBLEM - %li LAMDA models could not be found in chemistry network.\n\n\n", numModelsNotUsed );

		fclose(ioMASTERLIST);
	}

	//////////////////////////////////
	//  
	// Read CDMS/JPL masterlist
	//
	// These data files are in LAMDA format
	//
	//////////////////////////////////

	if( atmdat.lgCalpgmOn )
	{
		long numModelsNotUsed = 0;
		strcpy( chPath, "cdms+jpl" );
		strcat( chPath, input.chDelimiter );
		strcat( chPath, "masterlist" );
		ioMASTERLIST = open_data( chPath, "r" );

		if( read_whole_line( chLine , (int)sizeof(chLine) , ioMASTERLIST ) == NULL )
		{
			fprintf( ioQQQ, " database_readin could not read first line of CDMS/JPL masterlist.\n");
			cdEXIT(EXIT_FAILURE);
		}

		do
		{	
			if ((chLine[0]!='#') && (chLine[0]!='\n')&&(chLine[0]!='\t')&&(chLine[0]!='\r'))
			{	
				strcpy(chDLine, chLine);
				chToken = strtok(chDLine," \t\n");
				// hacks for alternative dialects...
				if( strcmp( chToken, "SH" ) == 0 )
					strcpy( chToken, "HS" );
				if( strcmp( chToken, "SH+" ) == 0 )
					strcpy( chToken, "HS+" );
				if( strcmp( chToken, "CCH" ) == 0 )
					strcpy( chToken, "C2H" );
				if( findspecies( chToken ) != null_mole  || 
					( chToken[1]=='-' && findspecies( chToken+2 ) != null_mole ) )
				{
					ASSERT( nSpecies + 1 <= MAX_NUM_SPECIES );
					ASSERT( nSpeciesLAMDA + 1 <= MAX_NUM_SPECIES );
					strcpy( chLabels[nSpecies], chToken );
					chLabels[nSpecies][CHARS_SPECIES-1] = '\0';

					strcpy( chPaths[nSpecies], "cdms+jpl" );
					strcat( chPaths[nSpecies], input.chDelimiter );
					chToken = strtok( NULL," \t\n" );
					strcat( chPaths[nSpecies], chToken );
					++nSpecies;
					++nSpeciesLAMDA;
				}
				else
					++numModelsNotUsed;
			}
		}
		while( read_whole_line( chLine , (int)sizeof(chLine) , ioMASTERLIST ) != NULL );

		// if( !t_version::Inst().lgRelease && numModelsNotUsed > 0 )
		//	fprintf( ioQQQ, "\nPROBLEM - %li CDMS/JPL models could not be found in chemistry network.\n\n",
		//		 numModelsNotUsed );

		fclose(ioMASTERLIST);
	}

	////////////////////////////////////
	//
	// Read STOUT masterlist and VERSION
	//
	///////////////////////////////////
	nSpeciesSTOUT = 0;
	if( atmdat.lgStoutOn )
	{
		// default location of Stout masterlist file
		strcpy( chPath, "stout" );
		strcat( chPath, input.chDelimiter );
		strcat( chPath, "masterlist" );
		strcat( chPath, input.chDelimiter );

		strcat( chPath, atmdat.chStoutFile );

		// first try local directory, then data/SED
		if( (ioMASTERLIST = fopen( atmdat.chStoutFile , "r" ) ) == NULL )
		{
			ioMASTERLIST = open_data( chPath, "r" );
		}

		// magic number
		if( read_whole_line( chLine , (int)sizeof(chLine) , ioMASTERLIST ) == NULL )
		{
			fprintf( ioQQQ, " database_readin could not read first line of stout.ini.\n");
			cdEXIT(EXIT_FAILURE);
		}

		bool lgEOLST;
		long int ipST = 1;
		long int nYrRdST = (long)FFmtRead(chLine,&ipST,sizeof(chLine),&lgEOLST);
		long int nMonRdST = (long)FFmtRead(chLine,&ipST,sizeof(chLine),&lgEOLST);
		long int nDayRdST = (long)FFmtRead(chLine,&ipST,sizeof(chLine),&lgEOLST);

		static long int nYrST =11 , nMonST = 10, nDayST = 25;
		if( ( nYrRdST != nYrST ) || ( nMonRdST != nMonST ) || ( nDayRdST != nDayST ) )
		{
			fprintf( ioQQQ,
				" I expected to find the number %2.2li %2.2li %2.2li and got %2.2li %2.2li %2.2li instead.\n" ,
				nYrST , nMonST , nDayST , nYrRdST , nMonRdST , nDayRdST );
			fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
			cdEXIT(EXIT_FAILURE);
		}
		if( read_whole_line( chLine , (int)sizeof(chLine) , ioMASTERLIST ) == NULL )
		{
			fprintf( ioQQQ, " database_readin could not read first line of CHIANTI masterlist.\n");
			cdEXIT(EXIT_FAILURE);
		}

		do
		{
			strcpy(chDLine, chLine);
			chToken = strtok(chDLine," \n");
			if ((chLine[0]!='#') && (chLine[0]!='\n')&&(chLine[0]!='\t')&&(chLine[0]!='\r'))
			{
				ASSERT( nSpecies + 1 <= MAX_NUM_SPECIES );
				ASSERT( nSpeciesSTOUT + 1 <= MAX_NUM_SPECIES );
				strcpy( chLabels[nSpecies], chToken );
				strcpy( chLabelsOrig[nSpecies], chLabels[nSpecies] );

				char *chElement, chTokenTemp[7];
				strcpy( chTokenTemp, chToken );
				chElement = strtok(chTokenTemp," \n");
				chElement = strtok(chTokenTemp,"_");
				uncaps( chElement );

//					printf("STOUT:%s\n",chLabels[nSpecies]);

				// path is, for example, CHIANTI/ar/ar_10/ar_10
				// we will append extensions later
				strcpy( chPaths[nSpecies], "stout" );
				strcat( chPaths[nSpecies], input.chDelimiter );
				strcat( chPaths[nSpecies], chElement );
				strcat( chPaths[nSpecies], input.chDelimiter );
				strcat( chPaths[nSpecies], chLabels[nSpecies] );
				strcat( chPaths[nSpecies], input.chDelimiter );
				strcat( chPaths[nSpecies], chLabels[nSpecies] );

				ASSERT( isalpha(chToken[0]) );
				long cursor=0;
				chLabels[nSpecies][0] = chToken[0];
				if( isalpha(chToken[1]) )
				{
					chLabels[nSpecies][1] = chToken[1];
					cursor = 2;
				}
				else
				{
					chLabels[nSpecies][1] = ' ';
					cursor = 1;
				}

				ASSERT( chToken[cursor++]=='_' );
				ASSERT( isdigit(chToken[cursor]) );

				if( isdigit(chToken[cursor+1]) )
				{
					chLabels[nSpecies][2] = chToken[cursor++];
					chLabels[nSpecies][3] =	chToken[cursor++];
				}
				else
				{
					chLabels[nSpecies][2] = ' ';
					chLabels[nSpecies][3] =	chToken[cursor++];
				}
				chLabels[nSpecies][4] = '\0';
				ASSERT( chToken[cursor]=='\0' || chToken[cursor]=='d' );

				// now capitalize the first letter
				chLabels[nSpecies][0] = toupper( chLabels[nSpecies][0] );
				++nSpecies;
				++nSpeciesSTOUT;
			}
		}
		while( read_whole_line( chLine , (int)sizeof(chLine) , ioMASTERLIST ) != NULL );
		fclose(ioMASTERLIST);
	}


	////////////////////////////////////
	//  
	// Read CHIANTI masterlist and VERSION
	//
	///////////////////////////////////

	nSpeciesCHIANTI = 0;

	if( atmdat.lgChiantiOn )
	{
		char chPathSave[FILENAME_PATH_LENGTH_2];
		strcpy( chPath, "chianti" );
		strcat( chPath, input.chDelimiter );
		//Preserve the path /chianti/ with chPathSave
		//Start reading in the chianti version number
		strcpy( chPathSave , chPath );
		strcat(chPath,"VERSION");
		ioVERSION = open_data(chPath,"r");
		if( read_whole_line( chLine , (int)sizeof(chLine) , ioVERSION ) == NULL )
		{
			fprintf( ioQQQ, " database_readin could not read first line of the Chianti VERSION.\n");
			cdEXIT(EXIT_FAILURE);
		}
		fclose(ioVERSION);
		// chianti version - string since can contain letters
		strncpy(atmdat.chVersion,chLine,atmdat.iVersionLength);
		// remove newline we may captured
		long len = strlen(atmdat.chVersion);
		if( atmdat.chVersion[len-1] == '\n' )
			atmdat.chVersion[len-1] = '\0';
		// may have been null earlier in string, but make sure terminated at limit
		atmdat.chVersion[atmdat.iVersionLength-1] = '\0';
		//Restore the previous chPath
		strcpy(chPath,chPathSave);
		// Read in the masterlist
		strcat( chPath, "masterlist" );
		strcat( chPath, input.chDelimiter );
		// save copy
		strcpy( chPathSave , chPath );

		// our subset of Chianti
		strcat( chPath, atmdat.chCloudyChiantiFile );

		// first try local directory, then data/chianti
		if( (ioMASTERLIST = fopen( atmdat.chCloudyChiantiFile , "r" ) ) == NULL )
		{
			ioMASTERLIST = open_data( chPath, "r" );
		}

		// magic number
		if( read_whole_line( chLine , (int)sizeof(chLine) , ioMASTERLIST ) == NULL )
		{
			fprintf( ioQQQ, " database_readin could not read first line of CloudyChianti.ini.\n");
			cdEXIT(EXIT_FAILURE);
		}

		bool lgEOL;
		long int ip = 1;
		long int nYrRd = (long)FFmtRead(chLine,&ip,sizeof(chLine),&lgEOL);
		long int nMonRd = (long)FFmtRead(chLine,&ip,sizeof(chLine),&lgEOL);
		long int nDayRd = (long)FFmtRead(chLine,&ip,sizeof(chLine),&lgEOL);

		static long int nYr=11 , nMon = 10, nDay = 3;
		if( ( nYrRd != nYr ) || ( nMonRd != nMon ) || ( nDayRd != nDay ) )
		{
			fprintf( ioQQQ,
				" database_readin: the version of CloudyChianti.ini is not the current version.\n" );
			fprintf( ioQQQ,
				" database_readin obtain the current version from the Cloudy web site.\n" );
			fprintf( ioQQQ,
				" I expected to find the number %2.2li %2.2li %2.2li and got %2.2li %2.2li %2.2li instead.\n" ,
				nYr , nMon , nDay , nYrRd , nMonRd , nDayRd );
			fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
			cdEXIT(EXIT_FAILURE);
		}

		if( read_whole_line( chLine , (int)sizeof(chLine) , ioMASTERLIST ) == NULL )
		{
			fprintf( ioQQQ, " database_readin could not read first line of CHIANTI masterlist.\n");
			cdEXIT(EXIT_FAILURE);
		}

		do
		{

			if ((chLine[0]!='#') && (chLine[0]!='\n')&&(chLine[0]!='\t')&&(chLine[0]!='\r'))
			{
				strcpy(chDLine, chLine);
				chToken = strtok(chDLine," \n");

				fixit();  //insert logic here to exclude some ions (for example, iso sequences)
				// exclude for now the satellite lines (denoted by a "d" after the label
				//if( chToken[3]=='d' || chToken[4]=='d' || chToken[5]=='d' )
				if( chToken[3]!='d' && chToken[4]!='d' && chToken[5]!='d' )
				{
					ASSERT( nSpecies + 1 <= MAX_NUM_SPECIES );
					ASSERT( nSpeciesCHIANTI + 1 <= MAX_NUM_SPECIES );
					strcpy( chLabels[nSpecies], chToken );
					strcpy( chLabelsOrig[nSpecies], chLabels[nSpecies]);

					bool skipSpecies = false;

					//Check for duplicate species with Stout
					for( int j = nSpeciesLAMDA; j < nSpecies; j++)
					{
						if( strcmp( chLabelsOrig[j], chLabelsOrig[nSpecies] ) == 0)
						{
							printf("Skipping the Chianti version of %s, using Stout version\n",chLabels[nSpecies]);
							skipSpecies = true;
							break;
						}
					}
					if( skipSpecies)
						continue;

					char *chElement, chTokenTemp[7];
					strcpy( chTokenTemp, chToken );
					chElement = strtok(chTokenTemp," \n");
					chElement = strtok(chTokenTemp,"_");
					uncaps( chElement );

					// path is, for example, CHIANTI/ar/ar_10/ar_10
					// we will append extensions later
					strcpy( chPaths[nSpecies], "chianti" );
					strcat( chPaths[nSpecies], input.chDelimiter );
					strcat( chPaths[nSpecies], chElement );
					strcat( chPaths[nSpecies], input.chDelimiter );
					strcat( chPaths[nSpecies], chLabels[nSpecies] );
					strcat( chPaths[nSpecies], input.chDelimiter );
					strcat( chPaths[nSpecies], chLabels[nSpecies] );

					ASSERT( isalpha(chToken[0]) );
					long cursor=0;
					chLabels[nSpecies][0] = chToken[0];
					if( isalpha(chToken[1]) )
					{
						chLabels[nSpecies][1] = chToken[1];
						cursor = 2;
					}
					else
					{
						chLabels[nSpecies][1] = ' ';
						cursor = 1;
					}

					ASSERT( chToken[cursor++]=='_' );
					ASSERT( isdigit(chToken[cursor]) );

					if( isdigit(chToken[cursor+1]) )
					{
						chLabels[nSpecies][2] = chToken[cursor++];
						chLabels[nSpecies][3] =	chToken[cursor++];
					}
					else
					{
						chLabels[nSpecies][2] = ' ';
						chLabels[nSpecies][3] =	chToken[cursor++];
					}
					chLabels[nSpecies][4] = '\0';
					ASSERT( chToken[cursor]=='\0' || chToken[cursor]=='d' );

					// now capitalize the first letter
					chLabels[nSpecies][0] = toupper( chLabels[nSpecies][0] );
					++nSpecies;
					++nSpeciesCHIANTI;
				}
			}
		}
		while( read_whole_line( chLine , (int)sizeof(chLine) , ioMASTERLIST ) != NULL );

		fclose(ioMASTERLIST);
	}

	/* no species found, nothing to do */
	if( nSpecies==0 )
		return;

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

	/*Initialization of the collisional rates array structure*/
	AtmolCollRateCoeff.reserve( nSpecies );
	AtmolCollSplines = (CollSplinesArray****)MALLOC((unsigned long)nSpecies *sizeof(CollSplinesArray***));
	StoutCollData = (StoutColls****)MALLOC((unsigned long)nSpecies *sizeof(StoutColls***));

	/*Mallocing here takes care of the number of colliders*/
	for( i=0; i<nSpecies; i++ )
	{
		AtmolCollRateCoeff.reserve( i, ipNCOLLIDER );
	}
	AtmolCollRateCoeff.alloc();

	// malloc state and transition arrays
	dBaseStates.resize(nSpecies);
	ipdBaseTrans.resize(nSpecies);

	for( i = 0; i < nSpecies; i++ )
	{
		dBaseTrans.push_back(TransitionList("dBaseTrans",&dBaseStates[i]));
		SpeciesJunk( &dBaseSpecies[i] );
		// label should be a minimum of 4 characters long
		size_t los = max(4,strlen(chLabels[i]));
		ASSERT( los >= 4 && los <= 7 );
		dBaseSpecies[i].chLabel = new char[los+1];
		strcpy(dBaseSpecies[i].chLabel,chLabels[i]);
		dBaseSpecies[i].lgActive = true;

		/* set type and isotopologue fractions */
		set_fractionation( &dBaseSpecies[i] );

		// set_fractionation trims off "p-","o-", etc.  Now have set label.  Check size.
		los = (int)strlen( dBaseSpecies[i].chLabel );
		ASSERT( los < CHARS_SPECIES );

		// pad label to at least four characters.
		if( dBaseSpecies[i].chLabel[2]=='\0' )
		{
			dBaseSpecies[i].chLabel[2]=' ';
			dBaseSpecies[i].chLabel[3]=' ';
			dBaseSpecies[i].chLabel[4]='\0';
		}
		else if( dBaseSpecies[i].chLabel[3]=='\0' )
		{
			dBaseSpecies[i].chLabel[3]=' ';
			dBaseSpecies[i].chLabel[4]='\0';
		}

		if( i<nSpeciesLAMDA )
		{
			// Read in LAMDA data files
			atmdat_LAMDA_readin( i, chPaths[i] );
		}
		else if( i < nSpeciesLAMDA + nSpeciesSTOUT )
		{
			atmdat_STOUT_readin( i, chPaths[i] );
		}
		else if( i < nSpeciesLAMDA + nSpeciesSTOUT + nSpeciesCHIANTI )
		{
			// Read in CHIANTI data files
			atmdat_CHIANTI_readin( i, chPaths[i] );
		}
		else
			TotalInsanity();
	}

	states_popfill();
	states_nelemfill();

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

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

STATIC void set_fractionation( species *sp )
{
	DEBUG_ENTRY("set_fractionation()");

	char chToken[3];

	sp->fracIsotopologue = 1.f;
	//types include "p-", "o-", "e-", and "a-"
	strncpy( chToken, sp->chLabel, 2 );
	chToken[2] = '\0';
	if( strcmp( "p-", chToken )==0 )
		sp->fracType = 0.25f;
	else if( strcmp( "o-", chToken )==0 )
		sp->fracType = 0.75f;
	else if( strcmp( "e-", chToken )==0 )
		sp->fracType = 0.5f;
	else if( strcmp( "a-", chToken )==0 ) 
		sp->fracType = 0.5f;
	else
		sp->fracType = 1.0f;

	fixit(); // what fraction should e-type and a-type Methanol have?  Assume 50/50 for now.

	// Now scrape the type specifier off the label.
	if( sp->chLabel[1]=='-')
		memmove(sp->chLabel,sp->chLabel+2,strlen(sp->chLabel+2)+1);

	return;
}

/*This function zeros the population of all states */
STATIC void states_popfill( void)
{
	DEBUG_ENTRY( "states_popfill()" );

	for( long i=0; i<nSpecies; i++)
	{
		for( long j=0; j<dBaseSpecies[i].numLevels_max; j++)
		{
			dBaseStates[i][j].Pop() = 0.;
		}
	}
	return;
}

/*This function fills the nelem and IonStg fields */
STATIC void states_nelemfill(void)
{
	DEBUG_ENTRY( "states_nelemfill()" );

	for( long i=0; i<nSpecies; i++ )
	{
		long nelem = 0, IonStg;
		char chLabelChemical[CHARS_SPECIES];

		if( dBaseSpecies[i].lgMolecular )
		{
			fixit(); 
			/* these should never be used if lgMolecular
			 *set to dangerous values instead of unity. */
			nelem = -1;
			IonStg = -1;
			strcpy( chLabelChemical, dBaseSpecies[i].chLabel );
		}
		else
		{
			char chToken[3];
			strncpy( chToken, dBaseSpecies[i].chLabel, 2 );
			chToken[2] = '\0';
			strcpy( chLabelChemical, chToken );
			if( chLabelChemical[1]==' ' )
				chLabelChemical[1] = '\0';
			for( long ipElement=0; ipElement<LIMELM; ipElement++ )
			{
				if( strcmp( elementnames.chElementSym[ipElement], chToken )==0 )
				{
					nelem = ipElement + 1;
					break;
				}
			}
			ASSERT( nelem > 0 && nelem <= LIMELM );
			strncpy( chToken, dBaseSpecies[i].chLabel + 2, 2 );
			IonStg = atoi(chToken);
			char chStage[5] = {'\0'};
			if( IonStg==2 )
				sprintf( chStage, "+" );
			else if( IonStg>1 )
				sprintf( chStage, "+%li", IonStg-1 );
			strcat( chLabelChemical, chStage );
			ASSERT( IonStg >= 1 && IonStg <= nelem+1 );
			//Prevent importing of iso-sequences from Chianti
			if( nelem - IonStg < NISO )
			{
				fprintf(ioQQQ, " PROBLEM: Cannot use Chianti model for %s%li\n",elementnames.chElementSym[nelem-1],IonStg);
				fprintf(ioQQQ, "  Iso-sequences are handled by our own model.\n");
				cdEXIT(EXIT_FAILURE);
			}
			dBaseSpecies[i].fmolweight = dense.AtomicWeight[nelem-1];
			// do not evaluate our cooling if we are using Chianti for this species

			if( dBaseTrans[i].chLabel() == "Chianti" )
			{
				dense.lgIonChiantiOn[nelem-1][IonStg-1] = true;
			}
			else if( dBaseTrans[i].chLabel() == "Stout" )
			{
				dense.lgIonStoutOn[nelem-1][IonStg-1] = true;
			}
			else
			{
				TotalInsanity();
			}

			if( atmdat.lgChiantiHybrid || atmdat.lgStoutHybrid )
			{
				// used in cool_dima to indicate whether to include line
				// with shorter wl than these databases
				dense.maxWN[nelem-1][IonStg-1] = dBaseSpecies[i].maxWN;
			}
			else
			{
				dense.maxWN[nelem-1][IonStg-1] = 0.;
			}
		}

		molecule *sp = findspecies(chLabelChemical);
		if( sp == null_mole )
		{
			dBaseSpecies[i].index = INT_MAX;
			if( nelem-1 >= ipHYDROGEN && dense.lgElmtOn[nelem-1] )
				fprintf(ioQQQ," PROBLEM: could not find species %li - %s\n",i,
					chLabelChemical );
		}
		else
		{
			dBaseSpecies[i].index = sp->index;
			mole.species[ sp->index ].levels = &dBaseStates[i];
			mole.species[ sp->index ].lines = &dBaseTrans[i];
		}

		for( long j=0; j<dBaseSpecies[i].numLevels_max; j++ )
		{
			dBaseStates[i][j].nelem() = nelem;
			dBaseStates[i][j].IonStg() = IonStg;
		}
	}
	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",dBaseSpecies[i].chLabel);
		printf("The data output is in the following format \n");
		printf("Label Energy St.wt Pop Lifetime\n");

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


STATIC void database_prep(int intSpIndex)
{
	vector<realnum> fsumAs(dBaseSpecies[intSpIndex].numLevels_max,SMALLFLOAT);

	DEBUG_ENTRY( "database_prep()" );

	/*Get the lifetimes*/
	for( EmissionList::iterator em = dBaseTrans[intSpIndex].Emis().begin();
		  em != dBaseTrans[intSpIndex].Emis().end(); ++em) 
	{
		fsumAs[(*em).Tran().ipHi()] += (*em).Aul();
		(*em).iRedisFun() = ipPRD;
	}
	
	dBaseStates[intSpIndex][0].lifetime()= BIGFLOAT;
	for( int ipHi=1; ipHi < dBaseSpecies[intSpIndex].numLevels_max; ipHi++ )
	{
		dBaseStates[intSpIndex][ipHi].lifetime() = 1./fsumAs[ipHi];
	}
	return;
}

/*SpeciesJunk set all elements of species struc to dangerous values */
STATIC void SpeciesJunk( species *sp )
{
	sp->chLabel = NULL;
	set_NaN(sp->fmolweight);
	set_NaN(sp->fracIsotopologue );
	set_NaN(sp->fracType);
	sp->lgMolecular = false;
	sp->numLevels_local = -INT_MAX;
	sp->numLevels_max = -INT_MAX;

	return;
}