/*ParseHDEN parse the hden command */
#include <string.h>
#include "cddefines.h"
#include "input.h"
#include "pressure.h"
#include "called.h"
#include "phycon.h"
#include "denpow.h"
#include "varypar.h"
#include "ffmtread.h"
#include "lgmatch.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;
	phycon.hden = (float)FFmtRead(chCard,&i,76,&lgEOL);
	if( lgEOL )
	{
		if( called.lgTalk )
		{
			fprintf( ioQQQ, " log(n) MUST be entered with this command. STOP\n" );
		}
		puts( "[Stop in ParseHDEN]" );
		exit(1);
	}

	/* check for further options */
	if( lgMatch("LINE",chCard) )
	{
		/* linear option */
		if( phycon.hden <= 0. )
		{
			fprintf( ioQQQ, " A negative density is not physical.\n" );
			puts( "[Stop in ParseHDEN]" );
			exit(1);
		}
		else
		{
			phycon.hden = (float)log10(phycon.hden);
		}
	}

	/* check if exponent entered */
	denpow.DensityPower = (float)FFmtRead(chCard,&i,76,&lgEOL);
	denpow.den0 = (float)pow(10.,phycon.hden);
	/* 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 */
			denpow.rscale = (float)pow(10.,FFmtRead(chCard,&i,76,&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]" );
				exit(1);
			}
			strcpy( pressure.chCPres, "POWC" );
		}
		else if( lgMatch("DEPT",chCard) )
		{
			/* depth option, sets scale radius, log cm */
			denpow.rscale = (float)pow(10.,FFmtRead(chCard,&i,76,&lgEOL));
			if( lgEOL )
			{
				fprintf( ioQQQ, " The scale depth MUST be set if the depth option is to be used.\n" );
				puts( "[Stop in ParseHDEN]" );
				exit(1);
			}
			strcpy( pressure.chCPres, "POWD" );
		}
		else
		{
			/* radius option, will be relative to inner radius */
			strcpy( pressure.chCPres, "POWR" );
		}
	}

	/* vary option */
	if( VaryPar.lgVarOn )
	{
		/*  pointer to where to write */
		VaryPar.nvfpnt[VaryPar.nparm] = input.nRead;
		VaryPar.vparm[0][VaryPar.nparm] = (float)phycon.hden;
		VaryPar.vincr[VaryPar.nparm] = 1.;

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

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

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

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

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


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