/* 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_radiative_acceleration radiative acceleration due to incident continuum, calls rt_line_driving for lines*/
#include "cddefines.h"
#include "physconst.h"
#include "rfield.h"
#include "opacity.h"
#include "wind.h"
#include "dense.h"
#include "rt.h"

void RT_radiative_acceleration(void)
{
	long int i;
	double rforce;

#	ifdef DEBUG_FUN
	fputs( "<+>RT_radiative_acceleration()\n", debug_fp );
#	endif
	/*radiative acceleration due to incident continuum, calls rt_line_driving for lines
	 *used for wind solution derives AccelLine and AccelCont */

	rforce = 0.;
	for( i=0; i < rfield.nflux; i++ )
	{
		rforce += (rfield.flux[i] + rfield.outlin[i] + rfield.outlin_noplot[i]+ rfield.ConInterOut[i])*
		  rfield.anu[i]*(opac.opacity_abs[i] + opac.opacity_sct[i]);
	}

	/* radiative acceleration; xMassDensity is gm per cc, eval when PTOT called */
	wind.AccelLine = (float)(RT_line_driving()/SPEEDLIGHT/dense.xMassDensity);
	wind.AccelCont = (float)(rforce*EN1RYD/SPEEDLIGHT/dense.xMassDensity);
	/* this is numerically unstable */
	wind.AccelPres = 0.;
	/* total acceleration */
	wind.AccelTot = wind.AccelCont + wind.AccelLine + wind.AccelPres;
	/*fprintf(ioQQQ," debuggrad acc\t%.2f %.2e %.2e %.2e\n",
		fnzone,wind.AccelLine, wind.AccelCont,opac.opacity_abs[997] );*/
	/*if( wind.AccelCont > 1. )
	{
		for( i=0; i < rfield.nflux; i++ )
		{
			rforce = (rfield.flux[i] + rfield.outlin[i] + rfield.outlin_noplot[i]+ rfield.ConInterOut[i])*
			rfield.anu[i]*(opac.opacity_abs[i] + opac.opacity_sct[i]);
			fprintf(ioQQQ,"debugg rforce\t%.3e\t%.2e\t%.2e\n",
				rfield.anu[i],
				(rfield.flux[i] + rfield.outlin[i] + rfield.outlin_noplot[i]+ rfield.ConInterOut[i]),
				(opac.opacity_abs[i] + opac.opacity_sct[i]) );
		}
	}*/

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