/* 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 */
/*fabden called by dlaw command, returns density for any density law */
#include "cddefines.h"
#include "physconst.h"
#include "dense.h"
#include "rfield.h"
#include "fabden.h"

double fabden(double radius, 
  double depth)
{
	double fabden_v;
	double z , z0 , density , rad_midplane;
#	ifdef DEBUG_FUN
	fputs( "<+>fabden()\n", debug_fp );
#	endif

	/* this routine is called when the DLAW command is given,
	 * and must return the hydrogen density at the current radius or depth
	 * RADIUS is the radius, the distance from the center of symmetry.
	 * DEPTH is the depth into the cloud, from the illuminated face
	 * both are in cm
	 *
	 * radius, depth, and the array DensityLaw are double precision, although
	 * FABDEN itself is not
	 *
	 * this is one way to generate a density
	 * fabden = radius * depth
	 *
	 * this is how to use the parameters on the dlaw command
	 * fabden = DensityLaw(1) + radius * DensityLaw(2)
	 *
	 * following must be removed if this sub is to be used */
	/*fprintf( ioQQQ, " the old version of FABDEN must be deleted, and a new one put in place.  Sorry.\n" );*/

	if( depth > radius )
		TotalInsanity();

	rad_midplane = radius * cos( dense.DensityLaw[0] );
	z = radius * sin( dense.DensityLaw[0] );
	z0 = 0.047 * pow( rad_midplane/AU , 1.25 ) * AU;
	/* density in gm / cm3 */
	density = 1.4e-9 * pow( rad_midplane/AU , -11./4. ) * sexp( POW2(z/z0) );
	/* convert density to cm-3 */
	fabden_v = density/(ATOMIC_MASS_UNIT * 1.2);

	fprintf(ioQQQ, "bug fabden zone %.2f radius %e density %e \n",
		fnzone , radius , fabden_v );

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

#if 0

	/* update stromgren radius to most recent value */
	if( rfield.lgUSphON )
		dense.DensityLaw[2] = rfield.rstrom/PARSEC;

	/* this is the density law used in the Wen & O'Dell, 1995, ApJ 438, 784 paper */
	powexp = MIN2((radius/parsec-dense.DensityLaw[2])/dense.DensityLaw[1],dense.DensityLaw[3]);
	fabden_v = pow(10.,dense.DensityLaw[0])*exp(powexp);

	fabden_v = dense.DensityLaw[0];
	/* just here to stop compilers from flagging unused vars */
	temp = radius + depth ; 
	if( fabden_v == 0. )
	{
		puts( "[Stop in fabden]" );
		cdEXIT(EXIT_FAILURE);
	}
	else
	{
		/* when this routine is used the following branck is the correct exit */
		
#		ifdef DEBUG_FUN
		fputs( " <->fabden()\n", debug_fp );
#		endif
		/* following should return correct value of the density at this position, 
		 * temp is there to trick lint */
		return( fabden_v*temp );
	}
#endif
}

