/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * For conditions of distribution and use see copyright notice in license.txt */
/*CO_PopsEmisCool evaluate rotation levels populations, emission, and cooling */
/*cdCO_colden return column density in H2, negative -1 if cannot find state,
 * header is cddrive */

#include "cddefines.h"
#include "taulines.h"
#include "dense.h"
#include "thermal.h"
#include "hmi.h"
#include "radius.h"
#include "atoms.h"
#include "phycon.h"
#include "rt.h"
#include "lines_service.h"
#include "cddrive.h"
#include "mole.h"
/*lint -e778 constant expression evaluatess to 0 in operation '-' */
/*=================================================================*/

/* will be used to save CO column densities */
static double *col12 , *col13;

/* evaluate CO rotation cooling */
void CO_PopsEmisCool(
		EmLine ** Rotate , 
		long int nRotate ,
		float abundan, 
		const char * chLabel ,
		float * Cooling ,
		float * dCoolingDT  )
{

	/* will need to MALLOC space for these but only on first call */
	static double 
		**AulEscp ,
		**col_str ,
		**AulDest, 
		/* AulPump[low][high] is rate (s^-1) from lower to upper level */
		**AulPump,
		**CollRate,
		*pops,
		*create,
		*destroy,
		*depart,
		/* statistical weight */
		*stat ,
		/* excitation energies in kelvin */
		*excit;

	/*static long int **ipdest;*/

	static int lgFirst=TRUE;
	long int i,
	j,
	ilo , 
	ihi;
	static long int nUsed;
	int lgDeBug,lgNegPop;
	int lgZeroPop;
	double rot_cooling , dCoolDT ;
	static long int ndimMalloced = 0;

#	ifdef DEBUG_FUN
	fputs( "<+>CO_PopsEmisCool()\n", debug_fp );
#	endif

	if ( lgFirst )
	{
		/* will never do this again */
		lgFirst = FALSE;
		/* remember how much space we malloced in case ever called with more needed */
		ndimMalloced = nRotate;
		/* allocate the 1D arrays*/
		if ( (excit = (double *)MALLOC( sizeof(double)*(size_t)(nRotate+1) )) == NULL )
			BadMalloc();
		if ( (stat = (double *)MALLOC( sizeof(double)*(size_t)(nRotate+1) )) == NULL )
			BadMalloc();
		if ( (pops = (double *)MALLOC( sizeof(double)*(size_t)(nRotate+1) )) == NULL )
			BadMalloc();
		if ( (create = (double *)MALLOC( sizeof(double)*(size_t)(nRotate+1) )) == NULL )
			BadMalloc();
		if ( (destroy = (double *)MALLOC( sizeof(double)*(size_t)(nRotate+1) )) == NULL )
			BadMalloc();
		if ( (depart = (double *)MALLOC( sizeof(double)*(size_t)(nRotate+1) )) == NULL )
			BadMalloc();
		/* create space for the 2D arrays */
		if ( (AulPump = ((double **)MALLOC((size_t)(nRotate+1)*sizeof(double *)))) == NULL )
			BadMalloc();
		if ( (CollRate = ((double **)MALLOC((size_t)(nRotate+1)*sizeof(double *)))) == NULL )
			BadMalloc();
		if ( (AulDest = ((double **)MALLOC((size_t)(nRotate+1)*sizeof(double *)))) == NULL )
			BadMalloc();
		if ( (AulEscp = ((double **)MALLOC((size_t)(nRotate+1)*sizeof(double *)))) == NULL )
			BadMalloc();
		if ( (col_str = ((double **)MALLOC((size_t)(nRotate+1)*sizeof(double *)))) == NULL )
			BadMalloc();
		for ( i=0; i<(nRotate+1); ++i )
		{
			if ( (AulPump[i] = ((double *)MALLOC((size_t)(nRotate+1)*sizeof(double )))) == NULL )
				BadMalloc();
			if ( (CollRate[i] = ((double *)MALLOC((size_t)(nRotate+1)*sizeof(double )))) == NULL )
				BadMalloc();
			if ( (AulDest[i] = ((double *)MALLOC((size_t)(nRotate+1)*sizeof(double )))) == NULL )
				BadMalloc();
			if ( (AulEscp[i] = ((double *)MALLOC((size_t)(nRotate+1)*sizeof(double )))) == NULL )
				BadMalloc();
			if ( (col_str[i] = ((double *)MALLOC((size_t)(nRotate+1)*sizeof(double )))) == NULL )
				BadMalloc();
		}
	}
	/* this is test for call with too many rotation levels to handle - logic needs
	 * for largest rotor to be called first */
	if ( nRotate > ndimMalloced )
	{
		fprintf(ioQQQ," CO_PopsEmisCool has been called with the number of rotor levels greater than space allocated.\n");
		puts( "[Stop in CO_PopsEmisCool]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* all elements are used, and must be set to zero if zero */
	for ( i=0; i < (nRotate+1); i++ )
	{
		create[i] = 0.;
		destroy[i] = 0.;
		for ( j=0; j < (nRotate+1); j++ )
		{
			col_str[j][i] = 0.;
			AulEscp[j][i] = 0.;
			AulDest[j][i] = 0.;
			AulPump[j][i] = 0.;
		}
	}

	/* the statistical weights of the levels */
	for ( j=0; j < nRotate; j++ )
	{
		/* statistical weights for each level */
		stat[j] = (*Rotate)[j].gLo;
	}
	/* this is the highest level, which is one more than the highest line */
	stat[nRotate] = (*Rotate)[nRotate-1].gHi;

	/* set up the excitation potentials of each level relative to ground -
	 * the struc saves the energy of the line only */
	excit[0] = 0.;
	for ( j=1; j < nRotate; j++ )
	{
		/* excitation energy of each level relative to ground, in K */
		excit[j] = excit[j-1] + (*Rotate)[j-1].EnergyK;
	}

	/* this is the highest level, which is one more than the highest line */
	excit[nRotate] = excit[nRotate-1] + (*Rotate)[nRotate-1].EnergyK;

	nUsed = nRotate;

	/* this determines the largest molecule that can be inverted at this
	 * temperature.  Need Boltzmann factors to be positive for all levels
	 * nUsed is the index of the highest level, so the number of levels (passed
	 * to the solver) is nUsed+1 
	 * excit[nUsed]/phycon.te cannot be much larger than 20, or matrix inversion will fail */
	/* >>chng 03 sep 18, allow to go to 1, a two level atom */
	/*while ( (excit[nUsed] > phycon.te*20.) && (nUsed > 1) )*/
	/* >>chng 03 oct 03, chng 20 to 25 */
	/*while ( (excit[nUsed] > phycon.te*25.) && (nUsed > 1) )*/
	/* >>chng 03 oct 03, keep at least 5 levels */
	/*while ( (excit[nUsed] > phycon.te*25.) && (nUsed > 1) )*/
	while ( (excit[nUsed] > phycon.te*25.) && (nUsed > 5) )
	{
		--nUsed;
	}

	for ( j=0; j < nRotate; j++ )
	{
		/*data[j][j+1] = (*Rotate)[j].Aul*((*Rotate)[j].Pesc + (*Rotate)[j].Pelec_esc);*/
		AulEscp[j+1][j] = (*Rotate)[j].Aul*((*Rotate)[j].Pesc + (*Rotate)[j].Pelec_esc);
		AulDest[j+1][j] = (*Rotate)[j].Aul*(*Rotate)[j].Pdest;
		/* next 2 not not used by levelN since flag passed saying to use CollRate instead */
		(*Rotate)[j].cs = 1.;
		/*data[j+1][j] = (*Rotate)[j].cs*hmi.Hmolec[ipMH2g]/dense.eden ;*/
		col_str[j+1][j] = (*Rotate)[j].cs*hmi.H2_total/dense.eden ;
		/* photon pumping rate */
		AulPump[j][j+1] = (*Rotate)[j].pump;
		/* the continuum indices are on the f, not c, scale, and will be passed to 
		 * atom_levelN, which works on f, not c, scale */
		/*ipdest[j][j+1] = (*Rotate)[j].ipCont;*/
	}

	/* next loop for collisional transitions that have delta J >1 */
	for ( ilo=0; ilo < nRotate-1; ilo++ )
	{
		/* need to have upper limit to to nRotate because 
		 * number ov levels is 1 gt number of lines*/
		for ( ihi=ilo+2; ihi <= nRotate; ihi++ )
		{
			/* this is not used by levelN since flag passed saying to use CollRate instead */
			/*data[ihi][ilo] = 1. *hmi.Hmolec[ipMH2g]/dense.eden ;*/
			/*data[ihi][ilo] = 1. *hmi.H2_total/dense.eden ;*/
			col_str[ihi][ilo] = 1. *hmi.H2_total/dense.eden ;
			/* these are escape and dest rates, which are zero for a rigid rotator */
			/*data[ilo][ihi] = 0 ;*/
			AulEscp[ihi][ilo] = 0 ;
			AulDest[ihi][ilo] = 0 ;
		}
	}

	/* now evaluate the H2 collision rates */
	/* recall one more level than lines */
	for ( ilo=0; ilo < nRotate; ilo++ )
	{
		/* need to have upper limit to to nRotate because 
		 * number of levels is 1 gt number of lines*/
		for ( ihi=ilo+1; ihi <= MIN2(nRotate,ilo+5); ihi++ )
		{
			/* >>refer	CO-H2	collision	de Jong, T., Chu, S-I., & Dalgarno, A. 1975, ApJ, 199, 69 */
			double a[5]={1.66, 2.80, 1.19, 1.00, 1.30};
			double b[5]={1.67, 1.47, 1.85, 1.55, 2.24};
			/* >>refer	CO-He	collision	McKee, C.F., Storey, J.W.V., Watson, D.M., & Green, S., 1982, ApJ, 259, 647 */
			/* this reference gives He-CO collisions,
			 * their table 1 says He collisions are 1.37 slower than H2 collisions*/
			/* >>chng 03 sep 19, from grnd H2 to total H2 */
			double collid = hmi.H2_total + dense.xIonDense[ipHELIUM][0]/1.37;
			/* de Jong et al. explicitly consider temperatures as low as 20K,
			 * don't extrapolate signifcantly lower than this */
			double TeUsed = MAX2(10., phycon.te );

			/* first do deexcitation rate, equation 17 of deJong et al. */
			CollRate[ihi][ilo] = a[ihi-ilo-1]*1.e-10*(*Rotate)[ilo].gLo/(*Rotate)[ilo].gHi*
				(1.+(excit[ihi]-excit[ilo])/TeUsed) *
				sexp( b[ihi-ilo-1]*sqrt((excit[ihi]-excit[ilo])/TeUsed) )*collid;
			/* this is mainly so that rest of code gets valid cs in case crit dense printed */
			if ( ihi == ilo+1 )
			{
				/* save downward collision rate */
				(*Rotate)[ilo].ColUL = (float)CollRate[ihi][ilo];
				/* convert into fake electron collision strength */
				LineConvRate2CS( &(*Rotate)[ilo] , (*Rotate)[ilo].ColUL);
			}

			/* now get excitation rate */
			CollRate[ilo][ihi] = CollRate[ihi][ilo]*
				sexp( (excit[ihi]-excit[ilo])/phycon.te )*
				(*Rotate)[ilo].gHi / (*Rotate)[ilo].gLo;

			/* debug print statement */
			/*fprintf(ioQQQ," %li %li %.2e %.2e \n",ilo, ihi, 
				CollRate[ihi][ilo]/hmi.H2_total  , CollRate[ilo][ihi]/hmi.H2_total );*/
		}
		/* finish off with zeros */
		for ( ihi=ilo+6; ihi <= nRotate; ihi++ )
		{
			CollRate[ihi][ilo] = 0.;
			CollRate[ilo][ihi] = 0.;
		}
	}

	lgDeBug = FALSE;

	atom_levelN(
		/* number of levels is number of lines plus one */
		/* set nUsed so that CO is evaluated even at very low temperatures */
		nUsed+1, /*nRotate+1,*/
		abundan,
		stat,
		excit,
		'K',
		pops,
		depart,
		&AulEscp,
		&col_str,
		&AulDest,
		&AulPump,
		&CollRate,
		create,
		destroy,
		/* say that we have evaluated the collision rates already */
		TRUE,
		/*&ipdest,*/
		&rot_cooling,
		&dCoolDT,
		chLabel,
		/* lgNegPop positive if negative pops occured, negative if too cold */
		&lgNegPop,
		&lgZeroPop,
		lgDeBug );/* option to print stuff - set to true for debug printout */

	/* put cooling into place where we can use it later */
	*Cooling = (float)rot_cooling;

	/* >>chng 03 sep 18, CO rot cooling temp deriv is too small, and temp
	 * changes too big - incr deriv by 5x 
	*dCoolingDT = (float)(dCoolDT);*/
	*dCoolingDT = (float)(rot_cooling/phycon.te);

#	if 0
	if( lgMainCO )
		fprintf(ioQQQ,"COcool\t%i\t%.2f\t%e\t%e\t%e\t%e\t%e\n",
		nUsed,
		fnzone,
		phycon.te,
		*Cooling/abundan,
		*dCoolingDT ,
		*Cooling,
		thermal.htot );/**/
#	endif

	/* zero out higher populations for case where full CO levels not done */
	for ( i=nUsed+1; i<=nRotate; ++i )
	{
		pops[i] = 0.;
		depart[i] = 0.;
	}

	/* can only define first LIMLEVELN elements, the vector's length */
	for ( i=0; i< MIN2(LIMLEVELN,nRotate) ; ++i )
	{
		atoms.PopLevels[i] = pops[i];
		atoms.DepLTELevels[i] = depart[i];
	}

	if ( lgNegPop > 0 )
	{
		fprintf(ioQQQ,"CO_PopsEmisCool called atom_levelN which returned negative populations.\n");
	}

	/* now establish information that is passed out to rest of code's infrastructure */
	for ( j=0; j<nRotate; ++j )
	{
		double EnrLU, EnrUL;
		/* lower upper populations, stim emission correct population */
		(*Rotate)[j].PopLo = pops[j];
		(*Rotate)[j].PopHi = pops[j+1];
		(*Rotate)[j].PopOpc = (pops[j] - pops[j+1]*(*Rotate)[j].gLo/(*Rotate)[j].gHi);

		/* number of photons in the line */
		(*Rotate)[j].phots = (*Rotate)[j].Aul*((*Rotate)[j].Pesc + (*Rotate)[j].Pelec_esc)*pops[j+1];

#		if 0
		/* >>chng 03 oct 04, moved to RT_OTS */
		/* local ots rates, added to line ots array */
		(*Rotate)[j].ots = (*Rotate)[j].Aul*(*Rotate)[j].Pdest*(float)(*Rotate)[j].PopHi;
		RT_OTS_AddLine( (*Rotate)[j].ots , (*Rotate)[j].ipCont);
#		endif

		/* the intensity in the line */
		(*Rotate)[j].xIntensity = (*Rotate)[j].phots*(*Rotate)[j].EnergyErg;

		/* ratio of collisional to total excitation */
		(*Rotate)[j].ColOvTot = (float)(CollRate[j][j+1]/(CollRate[j][j+1]+(*Rotate)[j].pump) );

		/* two cases - collisionally excited (usual case) or 
		 * radiatively excited - in which case line can be a heat source
		 * following are correct heat exchange, if will mix to get correct deriv */
		EnrLU = (*Rotate)[j].PopLo*CollRate[j][j+1]*(*Rotate)[j].EnergyErg;
		EnrUL = (*Rotate)[j].PopHi*CollRate[j+1][j]*(*Rotate)[j].EnergyErg;
		/* energy exchange due to this level
		 * net cooling due to excit minus part of de-excit */
		(*Rotate)[j].cool = EnrLU - EnrUL*(*Rotate)[j].ColOvTot;
		/* net heating is remainder */
		(*Rotate)[j].heat = EnrUL*(1.f - (*Rotate)[j].ColOvTot);
		/* do not add to cooling, since done with evaluated cooling from atom_levelN */
		/*CoolAdd( chLabel, (long)t10->WLAng , t10->cool);*/
	}

	/* generate flag if co cooling important and highest level is fainter
	 * than second highest level */
	if ( rot_cooling / thermal.ctot > 0.1 && 
		 (*Rotate)[nUsed-1].xIntensity > (*Rotate)[nUsed-2].xIntensity &&
		 /*>>chng 03 oct 03, add check whether molecule has been trimed down 
		  * due to small Boltzmann factor */ 
		 /* >>chng 03 oct 20, following had ; after ), so CoolCaped always true */
		 nUsed == nRotate )
	{
		co.lgCOCoolCaped = TRUE;
	}

#	ifdef DEBUG_FUN
	fputs( " <->CO_PopsEmisCool()\n", debug_fp );
#	endif

	return;
}

/*CO_Colden maintain H2 column densities within X */
void CO_Colden( const char *chLabel )
{
	long int iRot;

	if ( strcmp(chLabel,"ZERO") == 0 )
	{
		static int lgFIRST=TRUE;
		if ( lgFIRST )
		{
			lgFIRST = FALSE;
			if ( (col12 = (double *)MALLOC( sizeof(double)*(size_t)(nCORotate+1) )) == NULL )
				BadMalloc();
			if ( (col13 = (double *)MALLOC( sizeof(double)*(size_t)(nCORotate+1) )) == NULL )
				BadMalloc();
		}
		/* the column density (cm-2) of ortho and para H2 */
		/* zero out formation rates and column densites */
		for ( iRot=0; iRot<=nCORotate; ++iRot )
		{
			/* zero it out */
			col12[iRot] = 0.;
			col13[iRot] = 0.;
		}
	}
	else if ( strcmp(chLabel,"ADD ") == 0 )
	{
		/*  add together column densities */
		for ( iRot=0; iRot<nCORotate; ++iRot )
		{
			/* zero it out */
			col12[iRot] += C12O16Rotate[iRot].PopLo*radius.drad_x_fillfac;
			col13[iRot] += C13O16Rotate[iRot].PopLo*radius.drad_x_fillfac;
		}
	}

	/* we will not print column densities so skip that - if not print then we have a problem */
	else if ( strcmp(chLabel,"PRIN") != 0 )
	{
		fprintf( ioQQQ, " CO_Colden does not understand the label %s\n", 
				 chLabel );
		puts( "[Stop in CO_Colden]" );
		cdEXIT(EXIT_FAILURE);
	}
}

/*cdCO_colden return column density in H2, negative -1 if cannot find state,
 * header is cddrive */
double cdCO_colden( long isotope , long iRot )
{

	/* make sure incoming parameters are ok */
	if ( isotope !=12 && isotope != 13 )
	{
		fprintf(ioQQQ," cdCO_colden can only do 12CO and 13CO\n");
		return -1.;
	}
	if ( iRot < 0 || iRot > nCORotate )
	{
		fprintf(ioQQQ," cdCO_colden - rotation quantum number must be 0 or greater, and less than %li\n", 
				nCORotate);
		return -1.;
	}

	/* the incoming parameters are fully qualified - return the column density */
	if ( isotope == 12 )
	{
		return col12[iRot];
	}
	else
	{
		return col13[iRot];
	}
}

/*CO_OTS - add CO lines to ots fields */
void CO_OTS( void )
{

	long int j;

#	ifdef DEBUG_FUN
	fputs( "<+>CO_OTS()\n", debug_fp );
#	endif

	/* add all CO lines */
	for ( j=0; j<nCORotate; ++j )
	{
		C12O16Rotate[j].ots = C12O16Rotate[j].Aul*C12O16Rotate[j].Pdest*
			C12O16Rotate[j].PopHi;
		RT_OTS_AddLine( C12O16Rotate[j].ots , C12O16Rotate[j].ipCont);

		C13O16Rotate[j].ots = C13O16Rotate[j].Aul*C13O16Rotate[j].Pdest*
			C13O16Rotate[j].PopHi;
		RT_OTS_AddLine( C13O16Rotate[j].ots , C13O16Rotate[j].ipCont);
	}

#	ifdef DEBUG_FUN
	fputs( " <->CO_OTS()\n", debug_fp );
#	endif

	return;
}

