/* This file is part of Cloudy and is copyright (C) 1978-2004 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,
  /* string used to identify calling program in case of error */
  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
#	if 0
	double 
	  data[N_SEQ_BORON][N_SEQ_BORON], 
	  dest[N_SEQ_BORON][N_SEQ_BORON], 
	  /* pump[low][high] is rate (s^-1) from lower to upper level */
	  pump[N_SEQ_BORON][N_SEQ_BORON],
	  pops[N_SEQ_BORON],
	  /* statistical weight */
	  stat[N_SEQ_BORON] ,
	  /* excitation energies in kelvin */
	  excit[N_SEQ_BORON];
	long int ipdest[N_SEQ_BORON][N_SEQ_BORON];
#	endif
	static double 
	  **data, 
	  **dest, 
	  /* pump[low][high] is rate (s^-1) from lower to upper level */
	  **pump,
	  **CollRate,
	  *pops,
	  *create,
	  *destroy,
	  *depart,
	  /* statistical weight */
	  *stat ,
	  /* excitation energies in kelvin */
	  *excit;
	/*static long int **ipdest;*/

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

	static int lgFirst=TRUE;
	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( (pump = ((double **)MALLOC((N_SEQ_BORON)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (CollRate = ((double **)MALLOC((N_SEQ_BORON)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (dest = ((double **)MALLOC((N_SEQ_BORON)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (data = ((double **)MALLOC((N_SEQ_BORON)*sizeof(double *)))) == NULL )
			BadMalloc();
		for( i=0; i<(N_SEQ_BORON); ++i )
		{
			if( (pump[i] = ((double *)MALLOC((N_SEQ_BORON)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (CollRate[i] = ((double *)MALLOC((N_SEQ_BORON)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (dest[i] = ((double *)MALLOC((N_SEQ_BORON)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (data[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.;

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

		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;
			dest[j][i] = 0.;
			pump[j][i] = 0.;
		}
		data[i][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 kelvins */
	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);
	dest[0][1] = t10->Aul*t10->Pdest;
	data[1][0] = t10->cs;
	pump[0][1] = t10->pump;
	/* these indices are on the f, not c, scale, and will be passed to 
	 * atom_levelN, which works on f, not c, scale 
	ipdest[0][1] = t10->ipCont;*/

	data[0][2] = t20->Aul*(t20->Pesc + t20->Pelec_esc);
	dest[0][2] = t20->Aul*t20->Pdest;
	data[2][0] = t20->cs;
	pump[0][2] = t20->pump;
	/*ipdest[0][2] = t20->ipCont;*/

	data[0][3] = t30->Aul*(t30->Pesc + t30->Pelec_esc);
	dest[0][3] = t30->Aul*t30->Pdest;
	data[3][0] = t30->cs;
	pump[0][3] = t30->pump;
	/*ipdest[0][3] = t30->ipCont;*/

	data[0][4] = 1e-8;/* made up trans prob */
	dest[0][4] = 0.;
	data[4][0] = cs40;
	pump[0][4] = 0.;
	/*ipdest[0][4] = 0;*/

	data[1][2] = t21->Aul*(t21->Pesc + t21->Pelec_esc);
	dest[1][2] = t21->Aul*t21->Pdest;
	data[2][1] = t21->cs;
	pump[1][2] = t21->pump;
	/*ipdest[1][2] = t21->ipCont;*/

	data[1][3] = t31->Aul*(t31->Pesc + t31->Pelec_esc);
	dest[1][3] = t31->Aul*t31->Pdest;
	data[3][1] = t31->cs;
	pump[1][3] = t31->pump;
	/*ipdest[1][3] = t31->ipCont;*/

	data[1][4] = t41->Aul*(t41->Pesc + t41->Pelec_esc);
	dest[1][4] = t41->Aul*t41->Pdest;
	data[4][1] = t41->cs;
	pump[1][4] = t41->pump;
	/*ipdest[1][4] = t41->ipCont;*/

	data[2][3] = 1e-8;/* made up trans prob */
	dest[2][3] = 0.;
	data[3][2] = cs32;
	pump[2][3] = 0.;
	/*ipdest[2][3] = 0;*/

	data[2][4] = 1e-8;/* made up trans prob */
	dest[2][4] = 0.;
	data[4][2] = cs42;
	pump[2][4] = 0.;
	/*ipdest[2][4] = 0;*/

	data[3][4] = 1e-8;/* made up trans prob */
	dest[3][4] = 0.;
	data[4][3] = cs43;
	pump[3][4] = 0.;
	/*ipdest[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,
		pops,
		depart,
		&data,
		&dest,
		&pump,
		&CollRate,
		create,
		destroy,
		FALSE,/* say atom_levelN should evaluate coll rates from cs */
		/*&ipdest,*/
		&b_cooling,
		&dCoolDT,
		chLabel,
		&lgNegPop,
		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];

#		if 0
		/* >>chng 03 oct 03, moved to RT_OTS */
		t10->ots = t10->Aul*t10->Pdest*(float)pops[1];
		t20->ots = t20->Aul*t20->Pdest*(float)pops[2];
		t30->ots = t30->Aul*t30->Pdest*(float)pops[3];
		t21->ots = t21->Aul*t21->Pdest*(float)pops[2];
		t31->ots = t31->Aul*t31->Pdest*(float)pops[3];
		t41->ots = t41->Aul*t41->Pdest*(float)pops[4];

		/*  now add this line to ots field */
		RT_OTS_AddLine( t10->ots , t10->ipCont);
		RT_OTS_AddLine( t20->ots , t20->ipCont);
		RT_OTS_AddLine( t30->ots , t30->ipCont);
		RT_OTS_AddLine( t21->ots , t21->ipCont);
		RT_OTS_AddLine( t31->ots , t31->ipCont);
		RT_OTS_AddLine( t41->ots , t41->ipCont);
#		endif

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

