/* 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 */
/*AtomSeqBoron compute cooling from 5-level boron sequence model atom */
#include "cddefines.h"
#include "cooling.h"
#include "thermal.h"
#include "dense.h"
#include "atoms.h"

void AtomSeqBoron(
	/* indices for all lines are on the C scale since they will be stuffed into
	 * C arrays.  so, t10 referes to the 2-1 transition */
	EmLine * t10, 
	EmLine * t20, 
	EmLine * t30, 
	EmLine * t21, 
	EmLine * t31, 
	EmLine * t41, 
	double cs40,
	double cs32,
	double cs42,
	double cs43,
	/* pump rate s-1 due to UV permitted lines */
	double pump_rate ,
	/* string used to identify calling program in case of error */
	const char *chLabel
	)
{

	/* this routine has three possible returns:
	 * abundance is zero
	 * too cool for full 5-level atom, but still do ground term 
	 * full solution */

	/* boron sequence is now a five level atom */
#	define	N_SEQ_BORON	5
	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;

	double b_cooling,
	  dCoolDT;
	double EnrLU, EnrUL;
	float abundan;

	static int lgFirst=TRUE,
		lgZeroPop;
	int i , j;
	int 
		/* flag to signal negative level populations */
		lgNegPop,
		/* flag to turn on debug print in atom_levelN */
		lgDeBug;

	if( lgFirst )
	{
		/* will never do this again */
		lgFirst = FALSE;
		/* allocate the 1D arrays*/
		if( (excit = (double *)MALLOC( sizeof(double)*(N_SEQ_BORON) )) == NULL )
			BadMalloc();
		if( (stat = (double *)MALLOC( sizeof(double)*(N_SEQ_BORON) )) == NULL )
			BadMalloc();
		if( (pops = (double *)MALLOC( sizeof(double)*(N_SEQ_BORON) )) == NULL )
			BadMalloc();
		if( (create = (double *)MALLOC( sizeof(double)*(N_SEQ_BORON) )) == NULL )
			BadMalloc();
		if( (destroy = (double *)MALLOC( sizeof(double)*(N_SEQ_BORON) )) == NULL )
			BadMalloc();
		if( (depart = (double *)MALLOC( sizeof(double)*(N_SEQ_BORON) )) == NULL )
			BadMalloc();
		/* create space for the 2D arrays */
		if( (AulPump = ((double **)MALLOC((N_SEQ_BORON)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (CollRate = ((double **)MALLOC((N_SEQ_BORON)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (AulDest = ((double **)MALLOC((N_SEQ_BORON)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (AulEscp = ((double **)MALLOC((N_SEQ_BORON)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (col_str = ((double **)MALLOC((N_SEQ_BORON)*sizeof(double *)))) == NULL )
			BadMalloc();
		for( i=0; i<(N_SEQ_BORON); ++i )
		{
			if( (AulPump[i] = ((double *)MALLOC((N_SEQ_BORON)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (CollRate[i] = ((double *)MALLOC((N_SEQ_BORON)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (AulDest[i] = ((double *)MALLOC((N_SEQ_BORON)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (AulEscp[i] = ((double *)MALLOC((N_SEQ_BORON)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (col_str[i] = ((double *)MALLOC((N_SEQ_BORON)*sizeof(double )))) == NULL )
				BadMalloc();
		}
	}

	/* total abundance of this species */
	abundan = dense.xIonDense[ t10->nelem -1][t10->IonStg-1];
	if( abundan <= 0. )
	{
		/* this branch, no abundance of ion */
		t10->PopLo = 0.;
		t20->PopLo = 0.;
		t30->PopLo = 0.;
		t21->PopLo = 0.;
		t31->PopLo = 0.;
		t41->PopLo = 0.;

		t10->PopOpc = 0.;
		t20->PopOpc = 0.;
		t30->PopOpc = 0.;
		t21->PopOpc = 0.;
		t31->PopOpc = 0.;
		t41->PopOpc = 0.;

		t10->PopHi = 0.;
		t20->PopHi = 0.;
		t30->PopHi = 0.;
		t21->PopHi = 0.;
		t31->PopHi = 0.;
		t41->PopHi = 0.;

		t10->xIntensity = 0.;
		t20->xIntensity = 0.;
		t30->xIntensity = 0.;
		t21->xIntensity = 0.;
		t31->xIntensity = 0.;
		t41->xIntensity = 0.;

		t10->cool = 0.;
		t20->cool = 0.;
		t30->cool = 0.;
		t21->cool = 0.;
		t31->cool = 0.;
		t41->cool = 0.;

		t10->phots = 0.;
		t20->phots = 0.;
		t30->phots = 0.;
		t21->phots = 0.;
		t31->phots = 0.;
		t41->phots = 0.;

		t10->ColOvTot = 0.;
		t20->ColOvTot = 0.;
		t30->ColOvTot = 0.;
		t21->ColOvTot = 0.;
		t31->ColOvTot = 0.;
		t41->ColOvTot = 0.;

		t10->heat = 0.;
		t20->heat = 0.;
		t30->heat = 0.;
		t21->heat = 0.;
		t31->heat = 0.;
		t41->heat = 0.;

		CoolAdd( chLabel, t10->WLAng , 0.);
		CoolAdd( chLabel, t20->WLAng , 0.);
		CoolAdd( chLabel, t30->WLAng , 0.);
		CoolAdd( chLabel, t21->WLAng , 0.);
		CoolAdd( chLabel, t31->WLAng , 0.);
		CoolAdd( chLabel, t41->WLAng , 0.);

		/* level populations */
		/* LIMLEVELN is the dimension of the atoms vectors */
		ASSERT( N_SEQ_BORON <= LIMLEVELN);
		for( i=0; i < N_SEQ_BORON; i++ )
		{
			atoms.PopLevels[i] = 0.;
			atoms.DepLTELevels[i] = 1.;
		}
		
#		ifdef DEBUG_FUN
		fputs( " <->AtomSeqBoron()\n", debug_fp );
#		endif
		return;
	}

	ASSERT( t10->cs > 0.);
	ASSERT( t20->cs > 0.);
	ASSERT( t30->cs > 0.);
	ASSERT( t21->cs > 0.);
	ASSERT( t31->cs > 0.);
	ASSERT( t41->cs > 0.);
	ASSERT( cs40>0.);
	ASSERT( cs32>0.);
	ASSERT( cs42>0.);
	ASSERT( cs43>0.);

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

	/* statistical weights */
	stat[0] = t10->gLo;
	stat[1] = t10->gHi;
	stat[2] = t20->gHi;
	stat[3] = t30->gHi;
	stat[4] = t41->gHi;
	ASSERT( stat[0]>0. && stat[1]>0. &&stat[2]>0. &&stat[3]>0. &&stat[4]>0.);
	ASSERT( fabs(t10->gLo/2.-1.) < FLT_EPSILON);
	ASSERT( fabs(t10->gHi/4.-1.) < FLT_EPSILON);
	ASSERT( fabs(t20->gLo/2.-1.) < FLT_EPSILON);
	ASSERT( fabs(t20->gHi/2.-1.) < FLT_EPSILON);
	ASSERT( fabs(t30->gLo/2.-1.) < FLT_EPSILON);
	ASSERT( fabs(t30->gHi/4.-1.) < FLT_EPSILON);
	ASSERT( fabs(t21->gLo/4.-1.) < FLT_EPSILON);
	ASSERT( fabs(t21->gHi/2.-1.) < FLT_EPSILON);
	ASSERT( fabs(t31->gLo/4.-1.) < FLT_EPSILON);
	ASSERT( fabs(t31->gHi/4.-1.) < FLT_EPSILON);
	ASSERT( fabs(t41->gLo/4.-1.) < FLT_EPSILON);
	ASSERT( fabs(t41->gHi/6.-1.) < FLT_EPSILON);

	/* excitation energy of each level relative to ground, in Kelvin */
	excit[0] = 0.;
	excit[1] = t10->EnergyK;
	excit[2] = t20->EnergyK;
	excit[3] = t30->EnergyK;
	excit[4] = t41->EnergyK + t10->EnergyK;
	ASSERT( excit[1]>0. &&excit[2]>0. &&excit[3]>0. &&excit[4]>0.);

	/* fill in einstein As, collision strengths, pumping rates */
	/*data[0][1] = t10->Aul*(t10->Pesc + t10->Pelec_esc);*/
	AulEscp[1][0] = t10->Aul*(t10->Pesc + t10->Pelec_esc);
	AulDest[1][0] = t10->Aul*t10->Pdest;
	/*data[1][0] = t10->cs;*/
	col_str[1][0] = t10->cs;
	AulPump[0][1] = t10->pump;

	/* add FUV pump transitions to this pump rate */
	AulPump[0][1] += pump_rate;

	/*data[0][2] = t20->Aul*(t20->Pesc + t20->Pelec_esc);*/
	AulEscp[2][0] = t20->Aul*(t20->Pesc + t20->Pelec_esc);
	AulDest[2][0] = t20->Aul*t20->Pdest;
	/*data[2][0] = t20->cs;*/
	col_str[2][0] = t20->cs;
	AulPump[0][2] = t20->pump;

	/*data[0][3] = t30->Aul*(t30->Pesc + t30->Pelec_esc);*/
	AulEscp[3][0] = t30->Aul*(t30->Pesc + t30->Pelec_esc);
	AulDest[3][0] = t30->Aul*t30->Pdest;
	/*data[3][0] = t30->cs;*/
	col_str[3][0] = t30->cs;
	AulPump[0][3] = t30->pump;

	/*data[0][4] = 1e-8;*//* made up trans prob */
	AulEscp[4][0] = 1e-8;/* made up trans prob */
	AulDest[4][0] = 0.;
	/*data[4][0] = cs40;*/
	col_str[4][0] = cs40;
	AulPump[0][4] = 0.;

	/*data[1][2] = t21->Aul*(t21->Pesc + t21->Pelec_esc);*/
	AulEscp[2][1] = t21->Aul*(t21->Pesc + t21->Pelec_esc);
	AulDest[2][1] = t21->Aul*t21->Pdest;
	/*data[2][1] = t21->cs;*/
	col_str[2][1] = t21->cs;
	AulPump[1][2] = t21->pump;

	/*data[1][3] = t31->Aul*(t31->Pesc + t31->Pelec_esc);*/
	AulEscp[3][1] = t31->Aul*(t31->Pesc + t31->Pelec_esc);
	AulDest[3][1] = t31->Aul*t31->Pdest;
	/*data[3][1] = t31->cs;*/
	col_str[3][1] = t31->cs;
	AulPump[1][3] = t31->pump;

	/*data[1][4] = t41->Aul*(t41->Pesc + t41->Pelec_esc);*/
	AulEscp[4][1] = t41->Aul*(t41->Pesc + t41->Pelec_esc);
	AulDest[4][1] = t41->Aul*t41->Pdest;
	/*data[4][1] = t41->cs;*/
	col_str[4][1] = t41->cs;
	AulPump[1][4] = t41->pump;

	/*data[2][3] = 1e-8;*//* made up trans prob */
	AulEscp[3][2] = 1e-8;/* made up trans prob */
	AulDest[3][2] = 0.;
	/*data[3][2] = cs32;*/
	col_str[3][2] = cs32;
	AulPump[2][3] = 0.;

	/*data[2][4] = 1e-8;*//* made up trans prob */
	AulEscp[4][2] = 1e-8;/* made up trans prob */
	AulDest[4][2] = 0.;
	/*data[4][2] = cs42;*/
	col_str[4][2] = cs42;
	AulPump[2][4] = 0.;

	/*data[3][4] = 1e-8;*//* made up trans prob */
	AulEscp[4][3] = 1e-8;/* made up trans prob */
	AulDest[4][3] = 0.;
	/*data[4][3] = cs43;*/
	col_str[4][3] = cs43;
	AulPump[3][4] = 0.;

	lgDeBug = FALSE;
#	if 0
	if( strcmp( chLabel , "N  3" ) == 0)
		lgDeBug = TRUE;
#	endif

	/* lgNegPop positive if negative pops occured, negative if too cold */
	atom_levelN(N_SEQ_BORON,
		abundan,
		stat,
		excit,
		'K',
		pops,
		depart,
		&AulEscp,
		&col_str,
		&AulDest,
		&AulPump,
		&CollRate,
		create,
		destroy,
		FALSE,/* say atom_levelN should evaluate coll rates from cs */
		&b_cooling,
		&dCoolDT,
		chLabel,
		&lgNegPop,
		&lgZeroPop,
		lgDeBug );/* option to print suff - set to true for debug printout */

	/* this returned lgNegPop < 0 in case were gas too cool to excite highest level.
	 * in this case we still want to evaluate the cooling due to the split ground term
	 * in this sequence */
	if( lgNegPop < 0 )
	{
		/* it was too cool for atom_levelN to do whole atom, but let's still do the
		 * split ground term */
		atom_level2(t10);

		/* now put in pop for absorption into higher levels, for line optical depths,
		 * PopLevels was evaluated by atom_level2 */
		t20->PopLo = atoms.PopLevels[0];
		t30->PopLo = atoms.PopLevels[0];
		t21->PopLo = atoms.PopLevels[1];
		t31->PopLo = atoms.PopLevels[1];
		t41->PopLo = atoms.PopLevels[1];

		t20->PopOpc = atoms.PopLevels[0];
		t30->PopOpc = atoms.PopLevels[0];
		t21->PopOpc = atoms.PopLevels[1];
		t31->PopOpc = atoms.PopLevels[1];
		t41->PopOpc = atoms.PopLevels[1];

		t20->PopHi = 0.;
		t30->PopHi = 0.;
		t21->PopHi = 0.;
		t31->PopHi = 0.;
		t41->PopHi = 0.;

		t20->xIntensity = 0.;
		t30->xIntensity = 0.;
		t21->xIntensity = 0.;
		t31->xIntensity = 0.;
		t41->xIntensity = 0.;

		t20->cool = 0.;
		t30->cool = 0.;
		t21->cool = 0.;
		t31->cool = 0.;
		t41->cool = 0.;

		t20->phots = 0.;
		t30->phots = 0.;
		t21->phots = 0.;
		t31->phots = 0.;
		t41->phots = 0.;

#		if 0
		/* >>chng 03 oct 04, move to RT_OTS */
		t20->ots = 0.;
		t30->ots = 0.;
		t21->ots = 0.;
		t31->ots = 0.;
		t41->ots = 0.;
#		endif

		t20->ColOvTot = 0.;
		t30->ColOvTot = 0.;
		t21->ColOvTot = 0.;
		t31->ColOvTot = 0.;
		t41->ColOvTot = 0.;

		t20->heat = 0.;
		t30->heat = 0.;
		t21->heat = 0.;
		t31->heat = 0.;
		t41->heat = 0.;

		CoolAdd( chLabel, t20->WLAng , 0.);
		CoolAdd( chLabel, t30->WLAng , 0.);
		CoolAdd( chLabel, t21->WLAng , 0.);
		CoolAdd( chLabel, t31->WLAng , 0.);
		CoolAdd( chLabel, t41->WLAng , 0.);

		/* level populations */
		/* LIMLEVELN is the dimension of the atoms vectors */
		ASSERT( N_SEQ_BORON <= LIMLEVELN);
		for( i=2; i < N_SEQ_BORON; i++ )
		{
			atoms.PopLevels[i] = 0.;
			atoms.DepLTELevels[i] = 1.;
		}
		
	}
	else
	{
		/* atom_levelN did not evaluate PopLevels, so save pops here */
		/* LIMLEVELN is the dimension of the atoms vectors */
		ASSERT( N_SEQ_BORON <= LIMLEVELN);
		for( i=0; i< N_SEQ_BORON ; ++i )
		{
			atoms.PopLevels[i] = pops[i];
			atoms.DepLTELevels[i] = depart[i];
		}
		/* this branch, we have a full valid solution */
		t10->PopLo = pops[0];
		t20->PopLo = pops[0];
		t30->PopLo = pops[0];
		t21->PopLo = pops[1];
		t31->PopLo = pops[1];
		t41->PopLo = pops[1];

		t10->PopOpc = (pops[0] - pops[1]*t10->gLo/t10->gHi);
		t20->PopOpc = (pops[0] - pops[2]*t20->gLo/t20->gHi);
		t30->PopOpc = (pops[0] - pops[3]*t30->gLo/t30->gHi);
		t21->PopOpc = (pops[1] - pops[2]*t21->gLo/t21->gHi);
		t31->PopOpc = (pops[1] - pops[3]*t31->gLo/t31->gHi);
		t41->PopOpc = (pops[1] - pops[4]*t41->gLo/t41->gHi);

		t10->PopHi = pops[1];
		t20->PopHi = pops[2];
		t30->PopHi = pops[3];
		t21->PopHi = pops[2];
		t31->PopHi = pops[3];
		t41->PopHi = pops[4];

		t10->phots = t10->Aul*(t10->Pesc + t10->Pelec_esc)*pops[1];
		t20->phots = t20->Aul*(t20->Pesc + t20->Pelec_esc)*pops[2];
		t30->phots = t30->Aul*(t30->Pesc + t30->Pelec_esc)*pops[3];
		t21->phots = t21->Aul*(t21->Pesc + t21->Pelec_esc)*pops[2];
		t31->phots = t31->Aul*(t31->Pesc + t31->Pelec_esc)*pops[3];
		t41->phots = t41->Aul*(t41->Pesc + t41->Pelec_esc)*pops[4];

		t10->xIntensity = t10->phots*t10->EnergyErg;
		t20->xIntensity = t20->phots*t20->EnergyErg;
		t30->xIntensity = t30->phots*t30->EnergyErg;
		t21->xIntensity = t21->phots*t21->EnergyErg;
		t31->xIntensity = t31->phots*t31->EnergyErg;
		t41->xIntensity = t41->phots*t41->EnergyErg;

		/* ratio of collisional to total excitation */
		t10->ColOvTot = (float)(CollRate[0][1]/(CollRate[0][1]+t10->pump) );
		t20->ColOvTot = (float)(CollRate[0][2]/(CollRate[0][2]+t20->pump) );
		t30->ColOvTot = (float)(CollRate[0][3]/(CollRate[0][3]+t30->pump) );
		t21->ColOvTot = (float)(CollRate[1][2]/(CollRate[1][2]+t21->pump) );
		t31->ColOvTot = (float)(CollRate[1][3]/(CollRate[1][3]+t31->pump) );
		t41->ColOvTot = (float)(CollRate[1][4]/(CollRate[1][4]+t41->pump) );

		/* derivative of cooling function */
		thermal.dCooldT += dCoolDT;

		/* two cases - collisionally excited (usual case) or 
		 * radiatively excited - in which case line can be a heat source
		 * following are correct heat exchange, they will mix to get correct deriv 
		 * the sum of heat-cool will add up to EnrUL - EnrLU - this is a trick to
		 * keep stable solution by effectively dividing up heating and cooling,
		 * so that negative cooling does not occur */

		EnrLU = t10->PopLo*CollRate[0][1]*t10->EnergyErg;
		EnrUL = t10->PopHi*CollRate[1][0]*t10->EnergyErg;
		/* energy exchange due to this level
		 * net cooling due to excit minus part of de-excit */
		t10->cool = EnrLU - EnrUL*t10->ColOvTot;
		/* net heating is remainder */
		t10->heat = EnrUL*(1. - t10->ColOvTot);
		/* add to cooling */
		CoolAdd( chLabel, t10->WLAng , t10->cool);

		EnrLU = t20->PopLo*CollRate[0][2]*t20->EnergyErg;
		EnrUL = t20->PopHi*CollRate[2][0]*t20->EnergyErg;
		t20->cool = EnrLU - EnrUL*t20->ColOvTot;
		t20->heat = EnrUL*(1. - t20->ColOvTot);
		/* add to cooling */
		CoolAdd( chLabel, t20->WLAng , t20->cool);

		EnrLU = t30->PopLo*CollRate[0][3]*t30->EnergyErg;
		EnrUL = t30->PopHi*CollRate[3][0]*t30->EnergyErg;
		t30->cool = EnrLU - EnrUL*t30->ColOvTot;
		t30->heat = EnrUL*(1. - t30->ColOvTot);
		/* add to cooling */
		CoolAdd( chLabel, t30->WLAng , t30->cool);

		EnrLU = t21->PopLo*CollRate[1][2]*t21->EnergyErg;
		EnrUL = t21->PopHi*CollRate[2][1]*t21->EnergyErg;
		t21->cool = EnrLU - EnrUL*t21->ColOvTot;
		t21->heat = EnrUL*(1. - t21->ColOvTot);
		/* add to cooling */
		CoolAdd( chLabel, t21->WLAng , t21->cool);

		EnrLU = t31->PopLo*CollRate[1][3]*t31->EnergyErg;
		EnrUL = t31->PopHi*CollRate[3][1]*t31->EnergyErg;
		t31->cool = EnrLU - EnrUL*t31->ColOvTot;
		t31->heat = EnrUL*(1. - t31->ColOvTot);
		/* add to cooling */
		CoolAdd( chLabel, t31->WLAng , t31->cool);

		EnrLU = t41->PopLo*CollRate[1][4]*t41->EnergyErg;
		EnrUL = t41->PopHi*CollRate[4][1]*t41->EnergyErg;
		t41->cool = EnrLU - EnrUL*t41->ColOvTot;
		t41->heat = EnrUL*(1. - t41->ColOvTot);
		/* add to cooling */
		CoolAdd( chLabel, t41->WLAng , t41->cool);

	}

	return;
}

