/*HydroEinstA calculates Einstein A's from  osillator strengths*/
#include "cddefines.h"
#include "hydrooscilstr.h"
#include "hydroeinsta.h"

double HydroEinstA(long int lower, 
	  long int iupper)
{
	long int lo;
	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
	 * */
	if( lower*iupper == 2 )
	{
		/* 2s to 2p radiative transition */
		EinstA_v = 1e-12;
	}
	else if( lower == 0 && iupper == 1 )
	{
		/* 2s two photon */
		EinstA_v = 8.23;
	}
	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 )
		{
			/* 2s level */
			lo = 2;
			BranchingRatio = 0.25;
		}
		else if( lower == 2 )
		{
			/* 2p level */
			lo = 2;
			BranchingRatio = 0.75;
		}
		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);
		EinstA_v = HydroOscilStr(xl,xu)/(1.499e-8)/(POW2(xmicron))*xl*xl/
		  xu/xu*BranchingRatio;
	}

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