/*pop3 solve three level atom without radiative transfer */
#include "cddefines.h"
#include "phycon.h"
#include "pop3.h"
#include "edsqte.h"

double pop3(double g1, 
  double g2, 
  double g3, 
  double o12, 
  double o13, 
  double o23, 
  double a21, 
  double a31, 
  double a32, 
  double ex12, 
  double ex23, 
  float *pop2, 
  double abund, 
  double gam2)
{
	double alf, 
	  b12, 
	  b13, 
	  b23, 
	  bet, 
	  c12, 
	  c13, 
	  c21, 
	  c23, 
	  c31, 
	  c32, 
	  ex, 
	  fac, 
	  pop3_v;

#	ifdef DEBUG_FUN
	fputs( "<+>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 = ex12/phycon.te;
	if( (abund <= 0.) || (ex > 20.) )
	{
		pop3_v = 0.;
		*pop2 = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->pop3()\n", debug_fp );
#		endif
		return( pop3_v );
	}

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

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

	c12 = edsqte.cdsqte*o12/g1*b12;
	c13 = edsqte.cdsqte*o13/g1*b13;
	c23 = edsqte.cdsqte*o23/g2*b23;
	c32 = edsqte.cdsqte*o23/g3;
	c31 = edsqte.cdsqte*o13/g3;
	c21 = edsqte.cdsqte*o12/g2;

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

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

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