/* 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 */
/*ParseHDEN parse the hden command */
#include "cddefines.h"
#include "input.h"
#include "dense.h"
#include "optimize.h"
#include "parse.h"

void ParseHDEN(char *chCard )
{
	int lgEOL;
	long int i;

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

	/* log of hydrogen density */
	i = 5;
	dense.gas_phase[ipHYDROGEN] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
	if( lgEOL )
	{
		fprintf( ioQQQ, " log(n) MUST be entered with this command. STOP\n" );
		puts( "[Stop in ParseHDEN]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* check for further options */
	if( lgMatch("LINE",chCard) )
	{
		/* linear option */
		if( dense.gas_phase[ipHYDROGEN] <= 0. )
		{
			fprintf( ioQQQ, " A negative density is not physical.\n" );
			puts( "[Stop in ParseHDEN]" );
			cdEXIT(EXIT_FAILURE);
		}
		else
		{
			dense.gas_phase[ipHYDROGEN] = (float)log10(dense.gas_phase[ipHYDROGEN]);
		}
	}
	/* this is the linear initial density */
	dense.den0 = (float)pow(10.f,dense.gas_phase[ipHYDROGEN]);

	/* check if exponent entered */
	dense.DensityPower = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
	/* this branch when exponent was entered - do something with it */
	if( !lgEOL )
	{
		/*  not constant density
		 *  some sort of power law density distribution */
		if( lgMatch("COLU",chCard) )
		{
			/* density will depend on column density to a power
			 * number entered is col den, convert to scale radius
			 * at this point HDEN is LOG10 of hydrogen density */
			dense.rscale = (float)pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));
			if( lgEOL )
			{
				fprintf( ioQQQ, " The column density MUST be set if the col den option is to be used.\n" );
				puts( "[Stop in ParseHDEN]" );
				cdEXIT(EXIT_FAILURE);
			}
			strcpy( dense.chDenseLaw, "POWC" );
		}
		else if( lgMatch("DEPT",chCard) )
		{
			/* depth option, sets scale radius, log cm */
			dense.rscale = (float)pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));
			if( lgEOL )
			{
				fprintf( ioQQQ, " The scale depth MUST be set if the depth option is to be used.\n" );
				puts( "[Stop in ParseHDEN]" );
				cdEXIT(EXIT_FAILURE);
			}
			strcpy( dense.chDenseLaw, "POWD" );
		}
		else
		{
			/* radius option, will be relative to inner radius */
			strcpy( dense.chDenseLaw, "POWR" );
		}
	}

	/* vary option */
	if( optimize.lgVarOn )
	{
		/*  pointer to where to write */
		optimize.nvfpnt[optimize.nparm] = input.nRead;
		optimize.vparm[0][optimize.nparm] = (float)dense.gas_phase[ipHYDROGEN];
		optimize.vincr[optimize.nparm] = 1.;

		/* these are varios options for density laws, 
		 * first is constant density or pressre*/
		if( strcmp(dense.chDenseLaw ,"CDEN") == 0 || 
			 strcmp(dense.chDenseLaw ,"CPRE") == 0 )
		{
			strcpy( optimize.chVarFmt[optimize.nparm], "HDEN=%f" );
			optimize.nvarxt[optimize.nparm] = 1;
		}

		/* power law density distrution */
		else if( strcmp(dense.chDenseLaw,"POWR") == 0 )
		{
			strcpy( optimize.chVarFmt[optimize.nparm], "HDEN=%f, power=%f" );
			optimize.vparm[1][optimize.nparm] = dense.DensityPower;
			optimize.nvarxt[optimize.nparm] = 2;
		}

		/* power law with density scale depending on column density */
		else if( strcmp(dense.chDenseLaw,"POWC") == 0 )
		{
			strcpy( optimize.chVarFmt[optimize.nparm], "HDEN=%f, power=%f, column=%f" );
			optimize.nvarxt[optimize.nparm] = 3;
			optimize.vparm[1][optimize.nparm] = dense.DensityPower;
			optimize.vparm[2][optimize.nparm] = (float)log10(dense.rscale);
		}

		/* power law with density scale depending on depth */
		else if( strcmp(dense.chDenseLaw,"POWD") == 0 )
		{
			strcpy( optimize.chVarFmt[optimize.nparm], "HDEN=%f, power=%f, depth=%f" );
			optimize.nvarxt[optimize.nparm] = 3;
			optimize.vparm[1][optimize.nparm] = dense.DensityPower;
			optimize.vparm[2][optimize.nparm] = (float)log10(dense.rscale);
		}

		/* could not idensity an option */
		else
		{
			fprintf( ioQQQ, " Internal error in HDEN\n" );
			puts( "[Stop in ParseHDEN]" );
			cdEXIT(EXIT_FAILURE);
		}
		++optimize.nparm;
	}

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

