/* 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 */
/*atom_level2 do level population and cooling for two level atom,
 * side effects:
 * set elements of EmLine struc
 * cooling via 	CoolAdd( chLab, (long)t->WLAng , t->cool);
 * cooling derivative */
#include "cddefines.h"
#include "phycon.h"
#include "lines_service.h"
#include "dense.h"
#include "rfield.h"
#include "thermal.h"
#include "cooling.h"
#include "atoms.h"

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

#	ifdef DEBUG_FUN
	fputs( "<+>atom_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 CoolAdd
	 * CDSQTE is EDEN / SQRTE * 8.629E-6
	 */

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

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

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

	/* approximate boltzmann factor to see if results zero */
	boltz = rfield.ContBoltz[ip-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*phycon.cdsqte + t->pump)/(phycon.cdsqte + t->Aul);

	/* following first test needed for 32 bit cpu on search phase
	 * >>chng 96 jul 02, was 1e-30 crashed on dec, change to 1e-25
	 * if( AbunxIon.lt.1e-25 .or. boltz.gt.30. ) then
	 * >>chng 96 jul 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 = AbunxIon;
		t->PopOpc = AbunxIon;
		t->PopHi = 0.;
		t->xIntensity = 0.;
		t->cool = 0.;
		t->phots = 0.;
		t->ots = 0.;
		t->ColOvTot = 0.;
		t->heat = 0.;
		/* level populations */
		atoms.PopLevels[0] = AbunxIon;
		atoms.PopLevels[1] = 0.;
		atoms.DepLTELevels[0] = 1.;
		atoms.DepLTELevels[1] = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->atom_level2()\n", debug_fp );
#		endif
		return;
	}

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

	/* 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. );

	/* this lacks the upper statistical weight */
	col21 = phycon.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.);
	atoms.PopLevels[1] = pfs2;
	t->PopHi = pfs2;

	/* pop of ground */
	pfs1 = pfs2*r;
	atoms.PopLevels[0] = pfs1;

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

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

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

	/* ratio of collisional to total (collisional + pumped) 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, 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 */

	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 -
	 * note that ColOvTot cancels out in the sum heat - cool */
	coolng = Enr12 - Enr21*t->ColOvTot;
	t->cool = coolng;

	/* net heating is remainder */
	t->heat = 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 thermal.heating[0][22]
	 * in CoolSum */
	CoolAdd( chLab, t->WLAng , t->cool);

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

	/* destroyed part of line */
	/* >>chng 00 oct 31, diffuse continuum driving lines cause ionization
	 * runaway with atom_level2 lines close to H lyman edge.  Added factor
	 * of ColOvTot which should have been there all along. */
	/*t->ots = (float)(pfs2) * t->Aul * t->Pdest;*/
	/* >>chng 03 apr 26, add mean of old and new ots rates, NII 2040 drove
	 * oscillation in hizqso */
#	if 0
	otsold = t->ots;
#define	FRACNEW	0.5f
	t->ots = (float)(pfs2) * t->Aul * t->Pdest * t->ColOvTot;
	if( nzone )
		t->ots = t->ots*FRACNEW + otsold*(1.f-FRACNEW);
#undef FRACNEW
#	endif
#	if 0
	/*>>chng 03 oct 04, move to RT_OTS */
	/*  now add this line to ots field, ip is on the f not c scale */
	RT_OTS_AddLine( t->ots , ip);
#	endif
	 
	/* compute ratio Aul/(Aul+Cul) */
	/* level population with correction for stimulated emission */
	t->PopLo = atoms.PopLevels[0];

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

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

