/* 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 */
/*ParseCosmicRays parse the cosmic rays command */
#include "cddefines.h"
#include "optimize.h"
#include "hextra.h"
#include "ionbal.h"
#include "input.h"
#include "parse.h"

/*ParseCosmicRays parse the cosmic rays command */
void ParseCosmicRays( char *chCard )
{
	long int i;
	int lgEOL;
	float a;
	double var;

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

	/* cosmic ray density, log of rate relative to background, log of H0 rate in neutral gas,
	 * or density of rel. electrons,
	 * quantity is log unless keyword linear appears */
	i = 5;
	/* if no number is present FFmtRead returns zero */
	a = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
	if( lgEOL )
		a = 0.;

	/* if keyword LINEAR not present, then log, and make linear */
	if( !lgMatch("LINE",chCard) )
		a = (float)pow(10.f,a);
	/* a is now linear scale factor, or linear density */

	/* default is cosmic ray ionization rate relative to galactic background, but can
		* also give density, which was the only option originally */
	if( lgMatch("DENS",chCard) )
	{
		if( lgEOL )
		{
			fprintf(ioQQQ,
				" Either the log of the cosmic ray density must appear when density given.\n");
			NoNumb(chCard);
		}
		hextra.cryden = a;

		/*  optional power law density  */
		hextra.crpowr = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

		/*  option to specify a temp for non-rel electrons - but only when a density */
		hextra.crtemp = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		if( lgEOL )
		{
			/* relavistic limit (Balbus and McKee) */
			hextra.crtemp = 2.6e9;
		}
		else
		{
			var = pow(10.f,hextra.crtemp);
			hextra.crtemp = (float)MIN2(var,2.6e9);
		}
	}
	else if( lgMatch( "RATE" , chCard ) )
	{
		/* this sets rate - use stored density and rate for background to set
		 * new density since code works with density */
		ASSERT( a > 0. );
		hextra.cryden = hextra.background_density * a / hextra.background_rate;
		hextra.crtemp = 2.6e9f;
	}
	else
	{
		/* cr relative to galactic background BACK - no check on string since default */
		/* >>chng 04 mar 10, background is now 
			* >>refer	cr	ion	Williams, J.P., Bergin, E.A., Caseli, P., Myers, P.C., & Plume, R. 1998, ApJ, 503, 689 */
		/* galactic background cosmic ray density to produce
			* secondary ionization rate quoted by Tielens and Hollenbach */
		/* hextra.cryden = 2e-9f;*/
		/* >>chng 99 jun 24, slight change to value
			* quoted by 
			* >>refer	cosmic ray	ionization rate	McKee, C.M., 1999, astro-ph 9901370
			* this will produce a total
			* secondary ionization rate of 2.5e-17 s^-1, as tested in 
			* tsuite secondary.in.  If each ionization produces 2.4 eV of heat,
			* the background heating rate should be 9.6e-29 * n*/
		/* >>chng 00 nov 28, changed density to 4.9e-9 to reproduce TH85a
			* when photoionization is turned off. 
			>>refer	cosmic ray	ionization rate	Tielens, A.G.G.M., & Hollenbach, D., 1998, ApJ, 291, 722
			*/
		/* hextra.cryden = 7.07e-9f;*/
		/* this value reproduces the TH cr ionization rate when the factor
			* of 0.46 is included.  This will directly go onto the h ionization rate
			* without the factor of 0.46 there.  this is necessary for the more
			* general case where cr ionization is actually self-consistently determined
			* from rate hot electrons injected into the plasma */
		/*hextra.cryden = 2.25e-9f;*/
		ASSERT( a > 0. );
		hextra.cryden = hextra.background_density * a;
		hextra.crtemp = 2.6e9f;
	}

	/* this is current cosmic ray density divided by background - used in
	 * a few chemical reactions */
	hextra.cryden_ov_background = hextra.cryden / hextra.background_density;
	/* >>chng 05 jan 05, 
	 * set the cr ionization rate to very rough value, before we have enough
	 * information to evaluate it - may be needed in initial guess of H and He ionization*/
	ionbal.CosRayIonRate = hextra.cryden_ov_background * 2.5e-17;

	/* vary option */
	if( optimize.lgVarOn )
	{
		/* will be one parameter */
		optimize.nvarxt[optimize.nparm] = 1;
		strcpy( optimize.chVarFmt[optimize.nparm], "COSMic rays LOG BACKROUND=%f" );
		/* array index for where to write */
		optimize.nvfpnt[optimize.nparm] = input.nRead;
		/* log of cosmic rays rates realtive to background */
		optimize.vparm[0][optimize.nparm] = (float)log10(a);
		/* the increment in the first steps away from the original value */
		optimize.vincr[optimize.nparm] = 0.2f;
		++optimize.nparm;
	}

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

	return;
}

