/* 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 */
/*ParseFluc parse the fluctuations command, which affects either density or abundances */
#include "cddefines.h"
#include "dense.h"
#include "parse.h"

void ParseFluc(char *chCard )
{
	int lgEOL;
	long int i;
	double flmax, 
	  flmin, 
	  period, 
	  temp;

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

	/* rapid density fluctuations
	 * first parameter is log of period, 2 is log den max, 3 log Nmin */
	if( lgMatch("ABUN",chCard) )
	{
		/* abundances varied, not density */
		dense.lgDenFlucOn = FALSE;
	}
	else
	{
		/* density is varied */
		dense.lgDenFlucOn = TRUE;
	}

	/* optional keyword COLUMN makes sin over column density rather than radius */
	if( lgMatch("COLU",chCard) )
	{
		/* found key, not fluc over radius, over col den instead */
		dense.lgDenFlucRadius = FALSE;
	}
	else
	{
		/* no key, use default of radius */
		dense.lgDenFlucRadius = TRUE;
	}

	i = 5;
	/* 1st number log of period in centimeters */
	period = pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));
	dense.flong = (float)(6.2831853/period);
	temp = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

	/* 2nd number log of max hydrogen density */
	flmax = pow(10.,temp);

	/* 3rd number log of min hydrogen density */
	flmin = pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));
	if( flmax/flmin > 100. )
	{
		fprintf( ioQQQ, "This range of density probably will not work.\n" );
	}
	if( flmax > 1e15 )
	{
		fprintf( ioQQQ, "These parameters look funny to me.  Please check Hazy.\n" );
	}
	if( lgEOL || (flmin > flmax) )
	{
		fprintf( ioQQQ, "There MUST be three numbers on this line.\n" );
		fprintf( ioQQQ, "These must be the period(cm), max, min densities, all logs, in that order.\n" );
		puts( "[Stop in getfluc]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* this is optional phase shift for the command */
	dense.flcPhase = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

	/* FacAbunSav = (cfirst * COS( depth*flong+flcPhase ) + csecnd) */
	dense.cfirst = (float)((flmax - flmin)/2.);
	dense.csecnd = (float)((flmax + flmin)/2.);
	/* these will be added together with the first mult by sin - which goes to
	 * -1 - must not have a negative density */
	ASSERT( dense.cfirst < dense.csecnd );
	/* >>chng 96 jul 13 moved depset to SetAbundances fac
	 * if( lgDenFlucOn ) then
	 * this is a pressure law
	 * chCPres = 'SINE'
	 * else
	 * this is the metallicity of the gas
	 * do i=3,limelm
	 * depset(i) = flmax
	 * end do
	 * endif
	 *
	 * now get density if this is density option (not abundances) */
	if( dense.lgDenFlucOn )
	{
		strcpy( dense.chDenseLaw, "SINE" );

		/* depth is zero for first zone */
		dense.gas_phase[ipHYDROGEN] = dense.cfirst*(float)cos(dense.flcPhase) + dense.csecnd;
		dense.gas_phase[ipHYDROGEN] = (float)log10(dense.gas_phase[ipHYDROGEN]);
	}

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

