/* 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 */
/*ParseTLaw parse parameters on the tlaw command to set temperature as function of depth,
 * currently only does Bertoldi & Draine simple T law */
#include "cddefines.h"
#include "phycon.h"
#include "thermal.h"
#include "tfidle.h"
#include "parse.h"

void ParseTLaw(char *chCard )
{

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

	/* this says that some type of temperature law has been specified */
	thermal.lgTLaw = TRUE;
	thermal.lgTSetOn = TRUE;

	if( lgMatch("DB96",chCard) )
	{
		/* this is to simulate the temperature law given by equation 41 in 
		 * >>refer	H2	temperature law	Draine, B.T., & Bertoldi, Frank, 1996, ApJ, 468, 269-289 */
		thermal.lgTeBD96 = TRUE;

		/* this is the inital temperature for the BD96 temperature law */
		thermal.T0BD96 = 500.f;
		phycon.te = thermal.T0BD96;
		tfidle(FALSE);

		/* the coefficient on column density for temp dropoff */
		thermal.SigmaBD96 = 6e-22f;
	}
	else if( lgMatch("SN99",chCard) )
	{
		/* this is to simulate the temperature law given by equation 16 in 
		 * >>refer	H2	temperature law	Sternberg, A., & Neufeld, D.A. 1999, ApJ, 516, 371-380 */
		thermal.lgTeSN99 = TRUE;

		/* this is the inital temperature for the BD96 temperature law */
		thermal.T0SN99 = 500.f;
		phycon.te = thermal.T0SN99;
		tfidle(FALSE);
	}
	else
	{
		fprintf(ioQQQ," There must be a keyword on this command.  The one I know about is BD96\n");
		puts( "[Stop in ParseTLaw]" );
		cdEXIT(EXIT_FAILURE);
	}

#if 0
#include "dense.h"
#include "optimize.h"
#include "input.h"
	int lgEOL, 
	  lgEnd;
	long int i, 
	  j;
	char chCAP[5];
	/* all remainder is currently dead code, a copy of DLAW command,
	 * which could be activated if needs arose */
	/* call fcn FABDEN(RADIUS) which uses the ten parameters
	 * N.B.; existing version of FABDEN must be deleted
	 * >>chng 96 nov 29, added table option */
	if( lgMatch("TABL",chCard) )
	{
		/* when called, read in densities from input stream */
		strcpy( dense.chDenseLaw, "DLW2" );
		if( lgMatch("DEPT",chCard) )
		{
			dense.lgDLWDepth = TRUE;
		}
		else
		{
			dense.lgDLWDepth = FALSE;
		}

		input_readarray(chCard,&lgEnd);
		i = 1;
		dense.frad[0] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		dense.fhden[0] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		if( lgEOL )
		{
			fprintf( ioQQQ, " No pairs entered - can\'t interpolate.\n Sorry.\n" );
			puts( "[Stop in getdlaw]" );
			cdEXIT(EXIT_FAILURE);
		}

		dense.nvals = 2;
		lgEnd = FALSE;

		/* read pairs of numbers until we find line starting with END */
		while( !lgEnd && dense.nvals < LIMTABDLAW )
		{
			input_readarray(chCard,&lgEnd);
			if( !lgEnd )
			{
				cap4(chCAP , chCard);
				if( strncmp(chCAP , "END" , 3 ) == 0 )
					lgEnd = TRUE;
			}

			if( !lgEnd )
			{
				i = 1;
				dense.frad[dense.nvals-1] = (float)FFmtRead(chCard
				  ,&i,INPUT_LINE_LENGTH,&lgEOL);
				dense.fhden[dense.nvals-1] = (float)FFmtRead(chCard
				  ,&i,INPUT_LINE_LENGTH,&lgEOL);
				dense.nvals += 1;
			}
		}
		--dense.nvals;

		for( i=1; i < dense.nvals; i++ )
		{
			/* the radius values are assumed to be strictly increasing */
			if( dense.frad[i] <= dense.frad[i-1] )
			{
				fprintf( ioQQQ, " density.in radii must be in increasing order\n" );
				puts( "[Stop in getdlaw]" );
				cdEXIT(EXIT_FAILURE);
			}
		}
	}
	else
	{
		/* this is usual case, call fabden to get density */
		i = 4;
		for( j=0; j < 10; j++ )
		{
			dense.DensityLaw[j] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		}

		/* set flag so we know which law to use later */
		strcpy( dense.chDenseLaw, "DLW1" );

		/* vary option */
		if( optimize.lgVarOn )
		{
			/* NB - there are 5 = LIMEXT numbers on this line - if LIMEXT ever changes,
			 * chnage this too */
			strcpy( optimize.chVarFmt[optimize.nparm], "DLAW %f %f %f %f %f " );

			/* index for where to write */
			optimize.nvfpnt[optimize.nparm] = input.nRead;
			for( j=0; j<LIMEXT; ++j )
			{
				optimize.vparm[j][optimize.nparm] = (float)dense.DensityLaw[j];
			}
			optimize.vincr[optimize.nparm] = 0.5;
			optimize.nvarxt[optimize.nparm] = LIMEXT;
			++optimize.nparm;
		}
	}
#	endif

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