/* 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_pop3 solve 3-level atom without radiative transfer, returns pops of level 2 and 3 */
#include "cddefines.h"
#include "phycon.h"
#include "atoms.h"

/*atom_pop3 return value is population for 3-level atom, cm^-3 */
double atom_pop3(
	/* statictical weights of levels 1, 2, and 3 */
	double g1, double g2, double g3, 

	/* collision strengths between three levels */
	double o12, double o13, double o23, 

	/* transition probabilities between three levels */
	double a21, double a31, double a32, 

	/* excitation energy in Kelvin */
	double Tex12, double Tex23, 

	/* returned population of level 2, cm^-3 */
	float *pop2, 

	/* incoming total abundance of ion */
	double abund, 

	/* possible photodestruction of level 2, normally 0 */
	double gam2,
	
	/* excitation rates (s-1) due to "other" processes,
	 * these are not included in the energy exchange */
	double r12, 
	double r13 )
{
	double alf, 
	  b12, 
	  b13, 
	  b23, 
	  bet, 
	  c21, 
	  c23, 
	  c31, 
	  c32, 
	  ex, 
	  fac, 
	  pop3_v;

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

	/* computes level populations for 3 level atom, all col rad coupling
	 * results are populations of levels 2 and 3, (cm^-3) no A included */
	ex = Tex12/phycon.te;
	if( (abund <= 0.) || (ex > 20. && r12<SMALLFLOAT ) )
	{
		pop3_v = 0.;
		*pop2 = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->atom_pop3()\n", debug_fp );
#		endif
		return( pop3_v );
	}

	b12 = exp(-ex);
	b23 = exp(-Tex23/phycon.te);

	b13 = b12*b23;
	if( b13 == 0. && r12<SMALLFLOAT )
	{
		pop3_v = 0.;
		*pop2 = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->atom_pop3()\n", debug_fp );
#		endif
		return( pop3_v );
	}

	/* these rates have units s-1 */
	atoms.c12 = phycon.cdsqte*o12/g1*b12 + r12;
	atoms.c13 = phycon.cdsqte*o13/g1*b13 + r13;
	c23 = phycon.cdsqte*o23/g2*b23;
	c32 = phycon.cdsqte*o23/g3;
	c31 = phycon.cdsqte*o13/g3;
	c21 = phycon.cdsqte*o12/g2;

	alf = a21 + c21 + c23 + gam2;
	bet = a31 + a32 + c31 + c32;
	*pop2 = (float)((atoms.c13/bet + atoms.c12/(c32 + a32))/(alf/(c32 + a32) - c23/bet));
	pop3_v = (atoms.c13 + *pop2*c23)/bet;

	/* renorm to 1+pop2+atom_pop3=1 */
	fac = abund/(1. + *pop2 + pop3_v);
	*pop2 *= (float)fac;
	pop3_v *= fac;

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

	return( pop3_v );
}

