/*level2 do level population and cooling for two level atom */
#include <string.h>
#include "cddefines.h"
#include "taulines.h"
#include "ionfracs.h"
#include "phycon.h"
#include "edsqte.h"
#include "tcool.h"
#include "poplevls.h"
#include "rfield.h"
#include "chionlbl.h"
#include "sexp.h"
#include "addotslin.h"
#include "coladd.h"
#include "level2.h"

void level2(EmLine * t)
{
	char chLab[5];
	long int ion, 
	  ipCont, 
	  nelem;
	double AbunxIon, 
	  Enr12, 
	  Enr21, 
	  a21, 
	  aem, 
	  boltz, 
	  col12, 
	  col21, 
	  coolng, 
	  g1, 
	  g2, 
	  omega, 
	  pfs1,
	  pfs2,
	  plower, 
	  r, 
	  rate12, 
	  ri21;

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

	/* result is density (cm-3) of excited state times A21
	 * result normalized to N1+N2=ABUND
	 * routine increments dCooldT, call coladd
	 * CDSQTE is EDEN / SQRTE * 8.629E-6
	 */

	ion = t->IonStg ;
	nelem = t->nelem ;

	/* xIonFracs(nelem,i) is density of ith ionization stage (cm^-3) */
	AbunxIon = xIonFracs[ion][nelem-1];

	/* continuum pointer */
	ipCont = t->ipCont;

	/* approximate boltzmann factor to see if results zero */
	boltz = rfield.ContBoltz[ipCont-1];

	/* collision strength for this transition, omega is zero for hydrogenic
	 * species which are done in special hydro routines */
	omega = t->cs;

	/* ROUGH check whether upper level has significant population,*/
	r = (boltz*edsqte.cdsqte + t->pump)/(edsqte.cdsqte + t->Aul);

	/* following first test needed for 32 bit cpu on search phase
	 * >>chng 96 july 2, was 1e-30 crashed on dec, change to 1e-25
	 * if( AbunxIon.lt.1e-25 .or. boltz.gt.30. ) then
	 * >>chng 96 july 11, to below, since can be strong pumping when
	 * boltzman factor essentially zero */
	/* omega in following is zero for hydrogenic species, since done
	 * in hydro routines, so this should cause us to quit on this test */
	/* >>chng 99 nov 29, from 1e-25 to 1e-30 to keep same result for
	 * very low density models, where AbunxIon is very small but still significant*/
	/*if( omega*AbunxIon < 1e-25 || r < 1e-25 )*/
	if( omega*AbunxIon < 1e-30 || r < 1e-25 )
	{
		/* put in pop since possible just too cool */
		t->PopLo = (float)AbunxIon;
		t->PopOpc = (float)AbunxIon;
		t->PopHi = 0.;
		t->xIntensity = 0.;
		t->cool = 0.;
		t->phots = 0.;
		t->ots = 0.;
		/* >>>chng 99 nov 29, from 0 to 1 since very low density models have
		 * no coll excit but Putline came up with significant pump */
		t->AovTot = 0.;
		t->ColOvTot = 0.;
		t->heat = 0.;
		/* level populations */
		PopLevls.PopLevels[0] = (float)AbunxIon;
		PopLevls.PopLevels[1] = 0.;
		PopLevls.DepLTELevels[0] = 1.;
		PopLevls.DepLTELevels[1] = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->level2()\n", debug_fp );
#		endif
		return;
	}

	/* net rate down A21*(escape + destruction) */
	a21 = t->Aul*(t->Pesc + t->Pdest);

	/* put null terminated line label into chLab using line array*/
	chIonLbl(chLab,t);

	/* statistical weights of upper and lower levels */
	g1 = t->gLo;
	g2 = t->gHi;

	/* now get real Boltzmann factor */
	boltz = t->EnergyK/phycon.te;

	assert( boltz > 0. );
	boltz = sexp(boltz);

	assert( g1 > 0. && g2 > 0. );

	col21 = edsqte.cdsqte*omega;

	/* rate 1 to 2 is both collisions and pumping */
	col12 = col21/g1*boltz;
	col21 /= g2;

	/* the total excitation rate from lower to upper, collisional and pumping */
	rate12 = col12 + t->pump;

	/* induced emissions down */
	ri21 = t->pump*g1/g2;

	/* this is the ratio of lower to upper level populations */
	r = (a21 + col21 + ri21)/rate12;

	/* upper level pop */
	pfs2 = AbunxIon/(r + 1.);
	PopLevls.PopLevels[1] = (float)pfs2;
	t->PopHi = (float)pfs2;

	/* pop of ground */
	pfs1 = pfs2*r;
	PopLevls.PopLevels[0] = (float)pfs1;

	/* departure coef of excited state rel to ground */
	PopLevls.DepLTELevels[0] = 1.;
	if( boltz > 1e-20 && PopLevls.PopLevels[1] > 1e-20 )
	{
		/* this line could have zero boltz factor but radiatively excited
		 * dec alpha does not obey () in fast mode?? */
		PopLevls.DepLTELevels[1] = (float)((PopLevls.PopLevels[1]/PopLevls.PopLevels[0])/
		  (boltz*g2/g1));
	}
	else
	{
		PopLevls.DepLTELevels[1] = 0.;
	}

	/* number of escaping line photons, used elsewhere for outward beam */
	t->phots = t->Aul*t->Pesc*(float)pfs2;

	/* intensity of line */
	t->xIntensity = t->phots*t->EnergyErg;
	plower = AbunxIon - pfs2;

	/* ratio of collisional to total excitation */
	t->ColOvTot = (float)(col12/rate12);

	/* 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 */
	Enr12 = plower*col12*t->EnergyErg;
	Enr21 = pfs2*col21*t->EnergyErg;

	/* energy exchange due to this level
	 * net cooling due to excit minus part of de-excit */
	/*t->cool = Enr12 - Enr21*t->ColOvTot;*/
	coolng = Enr12 - Enr21*t->ColOvTot;
	/*coolng = t->cool;*/
	t->cool = (float)coolng;

	/* net heating is remainder */
	t->heat = (float)(Enr21*(1. - t->ColOvTot));

	/* expression pre jul 3 95, changed for case where line heating dominates
	 * coolng = (plower*col12 - pfs2*col21)*t->t(ipLnEnrErg)
	 * t->t(ipLnCool) = coolng */

	/* add to cooling - heating part was taken out above,
	 * and is not added in here - it will be added to heating(1,23)
	 * in CoolSum */
	coladd( chLab, (long)t->WLAng , coolng);

	/* derivative of cooling function */
	tcool.dCooldT += coolng * (t->EnergyK * tcool.tsq1 - tcool.halfte );

	/* destroyed part of line */
	t->ots = (float)(pfs2) * t->Aul * t->Pdest;

	/*  now add this line to ots field */
	AddOTSLin( t->ots , ipCont);
	 
	/* compute ratio Aul/(Aul+Cul) */
	aem = t->Aul*t->Pesc;
	t->AovTot = (float)(aem/(aem + col21 + ri21));

	/* level population with correction for stimulated emission */
	t->PopLo = PopLevls.PopLevels[0];

	t->PopOpc = 
		PopLevls.PopLevels[0] - PopLevls.PopLevels[1]*(float)(g1/g2);

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