/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * For conditions of distribution and use see copyright notice in license.txt */
/*RT_line_driving derive radiative acceleration due to line absorption of incident continuum */
#include "cddefines.h"
#include "physconst.h"
#include "iso.h"
#include "dense.h"
#include "taulines.h"
#include "h2.h"
#include "atomfeii.h"
#include "rt.h"

double RT_line_driving(void)
{
	long int i, 
	  ipHi, 
	  nelem, 
	  ipLo,
	  ipISO;

	double AllHeavy, 
	  AllRest, 
	  OneLine, 
	  fe2drive, 
	  forlin_v, 
	  h2drive,
	  accel_iso[NISO];

	/* following used for debugging */
	/* double 
	  RestMax, 
	  HeavMax, 
	  hydromax;
	  long int 
	  ipRestMax, 
	  ihmax; */

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

	/* this function finds the total rate the gas absorbs energy
	 * this result is divided by the calling routine to find the
	 * momentum absorbed by the gas, and eventually the radiative acceleration
	 *
	 * the energy absorbed by the line is
	 * Abundance * energy * A *(g_up/g_lo) * occnum * escape prob
	 * where occnum is the photon occupation number, and the g's are
	 * the ratios of statistical weights */

	/* total energy absorbed in this zone, per cubic cm
	 * do hydrogen first */

	for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
	{
		accel_iso[ipISO] = 0;
		for( nelem=ipISO; nelem < LIMELM; nelem++ )
		{
			if( (dense.IonHigh[nelem] >= nelem + 1-ipISO)  )
			{
				/* do not put in highest level since its not real */
				for( ipLo=0; ipLo < (iso.numLevels_max[ipISO][nelem] - 2); ipLo++ )
				{
					/* must not include the two-photon continuum in the folipLoing*/
					for( ipHi=ipLo+1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
					{
						/* do not include bogus lines */
						if( EmisLines[ipISO][nelem][ipHi][ipLo].ipCont > 0 )
						{
							OneLine = EmisLines[ipISO][nelem][ipHi][ipLo].pump*
							EmisLines[ipISO][nelem][ipHi][ipLo].EnergyErg*
							EmisLines[ipISO][nelem][ipHi][ipLo].PopOpc*
							dense.xIonDense[nelem][nelem+1-ipISO];

							accel_iso[ipISO] += OneLine;
						}
					}
				}
			}
		}
		accel_iso[ipISO] *= EN1RYD;
	}

	/* all heavy element lines in calculation of cooling 
	 * these are the level 1 lines */
	AllHeavy = 0.;
	for( i=1; i <= nLevel1; i++ )
	{
		OneLine = 
		  TauLines[i].pump*
		  TauLines[i].EnergyErg*
		  TauLines[i].PopOpc;
		AllHeavy += OneLine;
	}

	/* all heavy element lines treated with g-bar 
	 * these are the level 2 lines, f should be ok */
	AllRest = 0.;
	for( i=0; i < nWindLine; i++ )
	{
		OneLine = 
			TauLine2[i].pump*
			TauLine2[i].EnergyErg*
			TauLine2[i].PopOpc;
		AllRest += OneLine;
	}
	for( i=0; i < nUTA; i++ )
	{
		if( UTALines[i].Aul > 0. )
		{
			OneLine = 
				UTALines[i].pump*
				UTALines[i].EnergyErg*
				UTALines[i].PopOpc;
			AllRest += OneLine;
		}
	}
	for( i=0; i < nHFLines; i++ )
	{
		OneLine = 
			HFLines[i].pump*
			HFLines[i].EnergyErg*
			HFLines[i].PopOpc;
		AllRest += OneLine;
	}
	for( i=0; i < nCORotate; i++ )
	{
		OneLine = 
			C12O16Rotate[i].pump*
			C12O16Rotate[i].EnergyErg*
			C12O16Rotate[i].PopOpc;
		AllRest += OneLine;
		OneLine = 
			C13O16Rotate[i].pump*
			C13O16Rotate[i].EnergyErg*
			C13O16Rotate[i].PopOpc;
		AllRest += OneLine;
	}

	/* the H2 molecule */
	h2drive = H2_Accel();

	/* The large model FeII atom */
	fe2drive = 0.;
	FeIIAccel(&fe2drive);

	/*lint -e771 accel_iso not init */
	forlin_v = AllHeavy + accel_iso[ipH_LIKE] + accel_iso[ipHE_LIKE] + 
		fe2drive + h2drive + AllRest;
	/*lint +e771 accel_iso not init */

	/*fprintf(ioQQQ," wind te %e %e %e %e %e\n", 	
		AllHeavy , HydroAccel , fe2drive , he1l , AllRest );*/

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

