/*forlin derive radiative acceleration due to line absorption of incident continuum */
#include <string.h>
#include "cddefines.h"
#include "physconst.h"
#include "nhe1lvl.h"
#include "hydrogenic.h"
#include "rfield.h"
#include "taulines.h"
#include "ionfracs.h"
#include "he1tau.h"
#include "phe1lv.h"
#include "he.h"
#include "dumpline.h"
#include "iphe1l.h"
#include "ionrange.h"
#include "pop371.h"
#include "forlin.h"

double forlin(void)
{
	long int i, 
	  ipHi, 
	  ipZ, 
	  ipLo;

	double AllHeavy, 
	  AllRest, 
	  OneLine, 
	  fe2drive, 
	  forlin_v, 
	  he1l, 
	  HydroAccel;

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

#	ifdef DEBUG_FUN
	fputs( "<+>forlin()\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 */

	/* 666 error!
	 * put helium triplets in */

	/* total energy absorbed in this zone, per cubic cm
	 * do hydrogen first */
	HydroAccel = 0;

	/* do not put in highest level since its not real */
	for( ipZ=0; ipZ < LIMELM; ipZ++ )
	{
		if( (IonRange.IonHigh[ipZ] == ipZ + 2)  )
		{
			for( ipLo=IP1S; ipLo <= (nhlevel - 2); ipLo++ )
			{
				/* must not include the two-photon continuum in the folipLoing*/
				for( ipHi=MAX2(IP2P,ipLo+1); ipHi <= (nhlevel - 1); ipHi++ )
				{
					OneLine = HydroLines[ipZ][ipHi][ipLo].pump*
					  HydroLines[ipZ][ipHi][ipLo].EnergyErg*
					  HydroLines[ipZ][ipHi][ipLo].PopOpc*
					  /* >>>chng 99 sep 19, this did not have the following
					   * factor and caused build-up of radiation pressure in
					   * constant pressure models that went totally neutral */
					  xIonFracs[ipZ+2][ipZ];
					HydroAccel += OneLine;
				}
			}
		}
	}

	HydroAccel *= EN1RYD;

	/* total energy absorbed in this zone, per cubic cm
	 * helium */
	he1l = 0;
	for( ipLo=0; ipLo < 7; ipLo++ )
	{
		for( ipHi=ipLo + 1; ipHi < 8; ipHi++ )
		{
			he1l += phe1lv.he1n[ipLo]*he1tauCOM.he1con[ipLo][ipHi]*
			  rfield.anu[iphe1lCom.iphe1l[ipLo][ipHi]-1];
		}
	}
	he1l *= EN1RYD*he.heii;

	/* all heavy element lines in calculation of cooling 
	 * these are the level 1 lines */
	AllHeavy = 0.;
	/*HeavMax = 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;
	}

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

	forlin_v = AllHeavy + HydroAccel + fe2drive + he1l + AllRest;

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

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