/*popexc do level population for simple two level atom, no radiative transfer */
#include "cddefines.h"
#include "edsqte.h"
#include "teinv.h"
#include "showme.h"
#include "popexc.h"

double popexc(double omega, 
  double g1, 
  double g2, 
  double a21, 
  double bltz, 
  double abund)
{
	double boltz, 
	  popexc_v, 
	  q12, 
	  q21, 
	  r;

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

	/* result is density (cm-3) of excited state times a21
	 * result normalized to n1+n2=abund
	 * cdsqte is eden / sqrte * 8.629e-6
	 * */
	boltz = bltz*teinvCom.teinv;
	if( abund == 0. || boltz > 15. )
	{
		popexc_v = 0.;
		
#	ifdef DEBUG_FUN
		fputs( " <->popexc()\n", debug_fp );
#	endif
		return( popexc_v );
	}

	/*begin sanity check */
	if( omega <= 0. )
	{
		fprintf( ioQQQ, " POPEXC called with non-positive collision strength, =%10.2e\n", 
		  omega );
		fprintf( ioQQQ, " bltz=%10.2e\n", bltz );
		fprintf( ioQQQ, " args were omega , g1 , g2 , a21 , bltz , abund %12.5e%12.5e%12.5e%12.5e%12.5e%12.5e\n", 
		  omega, g1, g2, a21, bltz, abund );
		ShowMe();
		puts( "[Stop in popexc]" );
		exit(1);
	}
	/*end sanity check
	 * */
	q21 = edsqte.cdsqte*omega;
	q12 = q21/g1*exp(-boltz);
	q21 /= g2;
	r = (a21 + q21)/q12;
	popexc_v = abund*a21/(r + 1.);


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