/* 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 */
/*HydroEinstA calculates Einstein A's from  osillator strengths*/
#include "cddefines.h"
#include "hydrooscilstr.h"
#include "hydroeinsta.h"
#include "physconst.h"

double HydroEinstA(long int n1, 
	  long int n2)
{
	long int lo;
	long int lower, iupper;
	double EinstA_v, 
	  BranchingRatio, 
	  ryd, 
	  xl, 
	  xmicron, 
	  xu;

#	ifdef DEBUG_FUN
	fputs( "<+>HydroEinstA()\n", debug_fp );
#	endif
	/* (lower,upper) of Johnson 1972.  */

	/* lower = 0 is ground, n=1
	 * lower = 1 is 2s, n=2
	 * lower = 2 is 2p, n=2
	 * */

	/* sort out the order of upper and lower, so can be called either way */
	lower = MIN2( n1 , n2 );
	iupper = MAX2( n1, n2 );
	if( lower < 0 || lower == iupper )
	{
		fprintf(ioQQQ," HydroEinstA called with impossible ns, =%li %li\n", lower, iupper);
		cdEXIT(EXIT_FAILURE);
	}

	if( lower*iupper == 2 )
	{
		/* 2s to 2p radiative transition */
		EinstA_v = 1e-12;
	}
	else if( lower == 0 && iupper == 1 )
	{
		/* 2s two photon */
		/* >>chng 01 jan 26, from 8.23 to 8.226 */
		EinstA_v =  8.226;
	}
	else if( lower == 0 && iupper == 2 )
	{
		/* 2p 1s is special since do not assume 2s 2p mixed */
		EinstA_v = 6.265e8;
	}
	else
	{
		/* this is general case for E1 transitions*/
		if( lower == 0 )
		{
			/* 1s level */
			lo = 1;
			BranchingRatio = 1.;
		}
		else if( lower == 1 )
		{
			/* decays from above to 2s level */
			lo = 2;
			/* >>chng 01 sep 20, from 0.25 to 0.33, to get 2nu correct */
			/*BranchingRatio = 0.25;*/
			BranchingRatio = 0.34;
		}
		else if( lower == 2 )
		{
			/* decays from above to 2p level */
			lo = 2;
			/* >>chng 01 sep 20, from 0.75 to 0.67, to get 2nu correct */
			/*BranchingRatio = 0.75;*/
			BranchingRatio = 1.-0.34;
		}
		else
		{
			/* one of the higher levels that are only n states */
			lo = lower;
			BranchingRatio = 1.;
		}

		xl = (double)(MIN2(iupper,lo));
		xu = (double)(MAX2(iupper,lo));
		ryd = 1./POW2(xl) - 1./POW2(xu);
		/* xmicron = 1./(ryd*10.9737312); */
		xmicron = 1.E4/(ryd*RYD_INF);
		EinstA_v = HydroOscilStr(xl,xu)*TRANS_PROB_CONST*1e8f/(POW2(xmicron))*xl*xl/
		  xu/xu*BranchingRatio;
		
	}

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

