/* 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 */
/*ParseMagnet parse magnetic field command  */
/*Magnetic_init initialize magnetic field parameters */
/*Magnetic_reinit - reinitialized magnetic field at start of new iteration */
/*Magnetic_evaluate evaluate some parameters to do with magnetic field */
#include "cddefines.h"
#include "physconst.h"
#include "dense.h"
#include "doppvel.h"
#include "wind.h"
#include "magnetic.h"

/* the initial magnetic field */
static double Btangl_init;

/* this is logical var, set in zero, which says whether the magnetic field has been
	* initialized */
static int lgBinitialized;

/* the current magnetic field */
static double Btangl_here;

/* the inital parallel and tangential fields for ordered case */
static double Bpar_init, Btan_init;

/* the current parallel and tangential fields for ordered case */
static double Bpar_here, Btan_here;

/* this is the gamma power law index, default is 4. / 3. */
static double gamma_mag;

/*Magnetic_evaluate evaluate some parameters to do with magnetic field */
void Magnetic_evaluate(void)
{

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

	/* this flag set TRUE if magnetic field is specified */
	if( magnetic.lgB )
	{
		static double density_initial, 
			/* the square of the alfven velocity at illuminated face */
			v_A;

		/* does magnetic field need to be initialized for this iteration? 
		 * flag is reset false at init of code, and at start of every iteration */
		if( !lgBinitialized )
		{
			lgBinitialized = TRUE;

			/* set initial tangled field */
			Btangl_here = Btangl_init;

			/* set inital ordered field */
			/* mag field angle_wrt_los set when ordered field specified */
			Bpar_here = Bpar_init;
			Btan_here = Btan_init;

			/* XMassDensity was set above, safe to use this on first call */
			density_initial = dense.xMassDensity;

			/* this is used to update tangential field */
			v_A = POW2(Bpar_init) / (PI4 * density_initial );
		}

		/* now update parameters in tangled field case */
		/* magnetic pressure is a function of the local density, use gamma law */
		Btangl_here = Btangl_init * pow(dense.xMassDensity/density_initial, gamma_mag/2.);

		/* ordered components of field - parallel field is always constant - find tangential component -
		 * but only in wind case */
		if( wind.windv0 != 0. )
		{
			/* N B - must preserve sign in this equation - will blow if product of wind speeds is equal to v_A) */
			/* wind.windv*wind.windv0 == v_A should not occur since mag pressure goes to inf */
			Btan_here = Btan_init * (POW2(wind.windv0) - v_A)/ (wind.windv*wind.windv0-v_A);
		}

		/* magnetic pressure - sum of two fields - can have negative pressure (tension)
		 * is ordered field dominates */
		magnetic.pressure = POW2(Btangl_here)/PI8 + (POW2(Btan_here) - POW2(Bpar_here)) / PI8;

		/* energy density - this is positive */
		magnetic.energydensity = POW2(Btangl_here)/PI8 + (POW2(Btan_here) + POW2(Bpar_here)) / PI8;

		/* option for turbulence in equipartition with B field */
		if( DoppVel.lgTurbEquiMag )
		{
			/* >>chng 05 jan 26, as per Robin Williams email,
			 * evaluate energydensity above, which is +ve, and use that for
			 * velocity here - had used pressure but could not evaluate when negative */
			/* turbulent velocity is mag pres over density */
			DoppVel.TurbVel = (float)sqrt(2.*magnetic.energydensity/dense.xMassDensity);
			/* double magnetic pressure to account for ram pressure due to turbulence,
			 * wihch is not counted elsewhere */
			magnetic.pressure *= 2.;
		}

		/* input parser made sure gamma != 1, default magnetic gamma is 4/3 */
		magnetic.EnthalpyDensity = gamma_mag/(gamma_mag-1.) *
			POW2(Btangl_here)/PI8 + (POW2(Btan_here) + POW2(Bpar_here)) / PI4;
	}
	else
	{
		magnetic.pressure = 0.;
		magnetic.energydensity = 0.;
		magnetic.EnthalpyDensity = 0.;
	}

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

	return;
}

/*Magnetic_reinit - reinitialized magnetic field at start of new iteration */
void Magnetic_reinit(void)
{

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

	/* this says whether B has been initialized in this run */
	lgBinitialized = FALSE;

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

	return;
}

/* Magnetic_init initialize magnetic field parameters */
void Magnetic_init(void)
{

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

	gamma_mag = 4./3.;
	magnetic.lgB = FALSE;
	/* this says whether B has been initialized in this run */
	lgBinitialized = FALSE;
	/* the initial tangled and ordered fields */
	Btangl_init = 0.;
	Btangl_here = DBL_MAX;
	magnetic.pressure = DBL_MAX;
	magnetic.energydensity = DBL_MAX;
	Bpar_init = 0.;
	Btan_init = 0.;
	Bpar_here = DBL_MAX;
	Btan_here = DBL_MAX;
	magnetic.EnthalpyDensity = DBL_MAX;

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

	return;
}

/*ParseMagnet parse magnetic field command  */
void ParseMagnet(char *chCard )
{
	int lgEOL;
	long int i;

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

	/* flag saying B turned on */
	magnetic.lgB = TRUE;

	/* check whether ordered is specified - if not this is tangled */
	if( lgMatch("ORDE" , chCard ) )
	{
		double angle_wrt_los , Border_init;
		/* field is ordered, not isotropic - need field direction - spcified
		 * by angle away from beam of light - 0 => parallel to beam */

		i = 5;
		/* this is the log of the ordered field strength */
		Border_init = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		Border_init = pow(10.,Border_init );

		/* this is the angle (default in degrees) with respect to the line of sight */
		angle_wrt_los = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		if( lgEOL )
			NoNumb(chCard);

		/* if radians is on the line then the number already is in radians - 
		 * else convert to radians */
		if( !lgMatch("RADI" , chCard ) )
			angle_wrt_los *= PI2 / 360.;

		/* now get initial components that we will work with */
		Bpar_init = Border_init*cos(angle_wrt_los);
		Btan_init = Border_init*sin(angle_wrt_los);

	}
	else
	{
		i = 5;
		/* this is the log of the tangled field strength */
		Btangl_init = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		if( lgEOL )
			NoNumb(chCard);
		Btangl_init = pow(10.,Btangl_init );

		/* optional gamma for dependence on pressure */
		gamma_mag = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		if( lgEOL )
			gamma_mag = 4./3.;

		if( gamma_mag !=0. && gamma_mag <=1. )
		{
			/* impossible value for gamma */
			fprintf( ioQQQ, 
				" This value of gamma (%.3e) is impossible.  Must have gamma = 0 or > 1.\n Sorry.\n",
				gamma_mag );
			puts( "[Stop in ParseMagnet]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

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

	return;
}

