/* 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 "atmdat.h"
#include "lines_service.h"
#include "physconst.h"
#include "rfield.h"
#include "taulines.h"
#include "trace.h"
#include "string.h"
#include "thirdparty.h"
#include "mole.h"
#include "rfield.h"

#define DEBUGSTATE false

/*Separate Routine to read in the molecules*/
void atmdat_LAMDA_readin( long intNS, char *chEFilename )
{
	DEBUG_ENTRY( "atmdat_LAMDA_readin()" );

	int nMolLevs = -10000;
	long int intlnct,intrtct,intgrtct;

	/*intrtct refers to radiative transitions count*/
	FILE *ioLevData;
	realnum  fstatwt,fenergyK,fenergyWN,fenergy,feinsteina;

	char chLine[FILENAME_PATH_LENGTH_2];

	ASSERT( intNS >= 0 );

	const int MAX_NUM_LEVELS = 70;

	dBaseSpecies[intNS].lgMolecular = true;
	dBaseSpecies[intNS].lgLAMDA = true;

	/*Load the species name in the dBaseSpecies array structure*/
	if(DEBUGSTATE)
		printf("The name of the %li species is %s \n",intNS+1,dBaseSpecies[intNS].chLabel);

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

	ioLevData = open_data( chEFilename, "r" );

	nMolLevs = 0;
	intrtct = 0;
	intgrtct = 0;
	intlnct = 0;
	while(intlnct < 3)
	{
		intlnct++;
		if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == 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) , ioLevData ) == NULL )
	{
		fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
		cdEXIT(EXIT_FAILURE);
	}
	dBaseSpecies[intNS].fmolweight = (realnum)atof(chLine);

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

	// sections of the file are separated by line that begin with "!"
	ASSERT( chLine[0] == '!' );

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

	long HighestIndexInFile = nMolLevs;
	/* truncate model to preset max */
	nMolLevs = MIN2( nMolLevs, MAX_NUM_LEVELS );

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

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

	/*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() = "LAMDA";

	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;
		}
	}

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

	// sections of the file are separated by line that begin with "!"
	ASSERT( chLine[0] == '!' );

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

		// skip these levels
		if( ipLev >= nMolLevs )
			continue;

		/*information needed for label*/
		strcpy(dBaseStates[intNS][ipLev].chLabel(), "    ");
		strncpy(dBaseStates[intNS][ipLev].chLabel(),dBaseSpecies[intNS].chLabel, 4);

		// this is a hack to get ^13CO to print as 13CO in line lists	
		if( nMatch( "^13C", dBaseStates[intNS][ipLev].chLabel() ) )
			strcpy(dBaseStates[intNS][ipLev].chLabel(), "13CO");

		dBaseStates[intNS][ipLev].chLabel()[4] = '\0';
		// pad label to exactly four characters.
		if( dBaseStates[intNS][ipLev].chLabel()[2]=='\0' )
		{
			dBaseStates[intNS][ipLev].chLabel()[2]=' ';
			dBaseStates[intNS][ipLev].chLabel()[3]=' ';
		}
		else if( dBaseStates[intNS][ipLev].chLabel()[3]=='\0' )
		{
			dBaseStates[intNS][ipLev].chLabel()[3]=' ';
		}

		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 == ipLev + 1 );
		dBaseStates[intNS][ipLev].energy().set(fenergy,"cm^-1");
		dBaseStates[intNS][ipLev].g() = fstatwt;

		if( ipLev > 0 )
		{
			// Use volatile variables to ensure normalization to
			// standard precision before comparison
			volatile realnum enlev = dBaseStates[intNS][ipLev].energy().Ryd();
			volatile realnum enprev = dBaseStates[intNS][ipLev-1].energy().Ryd();
			if (enlev < enprev)
			{
				fprintf( ioQQQ, " The energy levels are not in order in species %s at index %li.\n",
					dBaseSpecies[intNS].chLabel, ipLev );
				cdEXIT(EXIT_FAILURE);
			}
		}
		if(DEBUGSTATE)
		{
			printf("The converted energy is %f \n",dBaseStates[intNS][ipLev].energy().WN());
			printf("The value of g is %f \n",dBaseStates[intNS][ipLev].g());
		}
	}

	/* 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( 1.01*rfield.emm*RYD_INF, dBaseStates[intNS][ipHi].energy().WN() - dBaseStates[intNS][ipLo].energy().WN() );
		fenergyWN = (realnum)(dBaseStates[intNS][ipHi].energy().WN() - dBaseStates[intNS][ipLo].energy().WN());
		fenergyK = (realnum)(fenergyWN*T1CM);
		
		(*tr).EnergyWN() = fenergyWN;

		/* there are OH hyperfine levels where i+1 and i have exactly
		 * the same energy.  The refractive index routine will FPE with
		 * an energy of zero - so we do this test */
		if( fenergyWN>SMALLFLOAT )
			(*tr).WLAng() = (realnum)(1e+8f/fenergyWN/RefIndex(fenergyWN));
		else
			(*tr).WLAng() = 1e30f;
	}

	if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == 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, expected to find line starting with!.\n",chEFilename);
		fprintf( ioQQQ , "%s\n", chLine );
		cdEXIT(EXIT_FAILURE);
	}
	if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == 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) , ioLevData ) == 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) , ioLevData ) == NULL )
		{
			fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
			cdEXIT(EXIT_FAILURE);
		}

		long i = 1;
		bool lgEOL;
		long index, ipHi, ipLo;

		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;

		ASSERT( ipHi >= 0 );
		ASSERT( ipLo >= 0 );
		ASSERT( index == intrtct + 1 );

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

		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 );
		fenergyWN = (realnum)(dBaseStates[intNS][ipHi].energy().WN() -dBaseStates[intNS][ipLo].energy().WN());
		fenergyWN = MAX2( fenergyWN, 1.01 * RYD_INF * rfield.emm );
		fenergyK = (realnum)(fenergyWN*T1CM);

		TransitionList::iterator tr = dBaseTrans[intNS].begin()+ipdBaseTrans[intNS][ipHi][ipLo];
		ASSERT(!(*tr).hasEmis()); // Check for duplicates
		(*tr).AddLine2Stack();
		(*tr).Emis().Aul() = feinsteina;
		ASSERT( !isnan( (*tr).EnergyK() ) );
		fenergyWN = (realnum)((fenergyK)/T1CM);
		(*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)
		{
			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",(*tr).Emis().Aul());
			printf("The Energy of the transition is %E \n",(*tr).EnergyK());
		}
	}

	if(read_whole_line( chLine , (int)sizeof(chLine)  , ioLevData ) == 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);
	}

	long nCollPartners = -1;
	/*Getting the number of collisional partners*/
	if(read_whole_line( chLine , (int)sizeof(chLine)  , ioLevData ) == NULL )
	{
		fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
		cdEXIT(EXIT_FAILURE);
	}
	else
	{
		nCollPartners = atoi(chLine);
	}
	/* setting the number of colliders to a negative value is a flag to do the molecule in LTE */
	dBaseSpecies[intNS].lgLTE = ( nCollPartners < 0 );
	/*Checking the number of collisional partners does not exceed 9*/
	if( nCollPartners > ipNCOLLIDER )
	{
		fprintf( ioQQQ, " The number of colliders is greater than what is expected in file %s.\n", chEFilename );
		cdEXIT(EXIT_FAILURE);
	}

	FunctPtr func = new FunctLAMDA();
	// loop over partners	
	for( long ipPartner = 0; ipPartner < nCollPartners; ++ ipPartner )
	{
		if(read_whole_line( chLine , (int)sizeof(chLine)  , ioLevData ) == NULL )
		{
			fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
			cdEXIT(EXIT_FAILURE);
		}

		ASSERT( chLine[0] == '!' );

		if(read_whole_line( chLine , (int)sizeof(chLine)  , ioLevData ) == 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*/
		char *chCollName = strtok(chLine," ");
		/*Leiden Collider Index*/
		int intLColliderIndex = atoi(chCollName);
		int intCollIndex = -1;
		/*In Cloudy,We assign the following indices for the colliders*/
		/*electron=0,proton=1,He+=2,He++=3,atomic hydrogen=4,He=5,oH2=6,pH2=7,H2=8*/

		if(intLColliderIndex == 1)
		{
			intCollIndex = ipH2;
		}
		else if(intLColliderIndex == 2)
		{
			intCollIndex = ipH2_PARA;
		}
		else if(intLColliderIndex == 3)
		{
			intCollIndex = ipH2_ORTHO;
		}
		else if(intLColliderIndex == 4)
		{
			intCollIndex = ipELECTRON;
		}
		else if(intLColliderIndex == 5)
		{
			intCollIndex = ipATOM_H;
		}
		else if(intLColliderIndex == 6)
		{
			intCollIndex = ipATOM_HE;
		}
		else
		{
			// this happens for some LAMDA files (as of Jan 20, 2009) because there is no integer
			// in the datafile indicating which collider the subsequent table is for
			TotalInsanity();
		}

		ASSERT( intCollIndex < ipNCOLLIDER );

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

		ASSERT( chLine[0] == '!' );

		if(read_whole_line( chLine , (int)sizeof(chLine)  , ioLevData ) == NULL )
		{
			fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
			cdEXIT(EXIT_FAILURE);
		}
		// Get the number of transitions 
		long nCollTrans = atoi(chLine);

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

		ASSERT( chLine[0] == '!' );

		if(read_whole_line( chLine , (int)sizeof(chLine)  , ioLevData ) == NULL )
		{
			fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
			cdEXIT(EXIT_FAILURE);
		}
		// Get the number of temperatures
		long intCollTemp = atoi(chLine);

		// Read and store the table of rate coefficients
		ReadCollisionRateTable( AtmolCollRateCoeff[intNS][intCollIndex],
			ioLevData, func, nMolLevs, intCollTemp, nCollTrans );
	}
	delete func;

	fclose( ioLevData );

	return;
}