/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  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 "optimize.h"
#include "input.h"
#include "wind.h"
#include "magnetic.h"
#include "parser.h"

t_magnetic magnetic;

/* 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 bool lgBinitialized;

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

/* the initial 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)
{

	DEBUG_ENTRY( "Magnetic_evaluate()" );

	/* 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 initial 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.lgStatic() )
		{
			/* 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 */
			/* >>chng 06 apr 19, use DoppVel.Heiles_Troland_F */
			/*DoppVel.TurbVel = (realnum)sqrt(2.*magnetic.energydensity/dense.xMassDensity);*/
			DoppVel.TurbVel = (realnum)sqrt(6.*magnetic.energydensity/dense.xMassDensity/
				DoppVel.Heiles_Troland_F);
			/* >>chng 06 apr 19, do not double mag pressure, really count turbulence as pressure */
			/* double magnetic pressure to account for ram pressure due to turbulence,
			 * which 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.;
	}
	return;
}

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

	/* this says whether B has been initialized in this run */
	lgBinitialized = false;
	return;
}

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

	DEBUG_ENTRY( "Magnetic_init()" );

	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;
	return;
}

/*ParseMagnet parse magnetic field command  */
void ParseMagnet(Parser &p )
{
	bool lgTangled;
	double angle_wrt_los=-1. , Border_init=-1.;

	DEBUG_ENTRY( "ParseMagnet()" );

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

	/* check whether ordered is specified - if not this is tangled */
	if( p.nMatch("ORDE") )
	{
		/* ordered field case */
		lgTangled = false;

		/* field is ordered, not isotropic - need field direction - spcified
		 * by angle away from beam of light - 0 => parallel to beam */

		/* this is the log of the ordered field strength */
		Border_init = p.getNumberCheckAlwaysLog("ordered field");

		/* this is the angle (default in degrees) with respect to the line of sight */
		angle_wrt_los = p.getNumberCheck("LOS angle");

		/* if radians is on the line then the number already is in radians - 
		 * else convert to radians */
		if( !p.nMatch("RADI") )
			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
	{
		/* tangled field case */
		lgTangled = true;
		/* this is the log of the tangled field strength */
		Btangl_init = p.getNumberCheckAlwaysLog("tangled field");

		/* optional gamma for dependence on pressure */
		gamma_mag = p.getNumberDefault("field gamma law", 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 );
			cdEXIT(EXIT_FAILURE);
		}
	}

	/* vary option */
	if( optimize.lgVarOn )
	{
		/* number of parameters */
		optimize.nvarxt[optimize.nparm] = 2;
		if( lgTangled )
		{
			/* tangled field case */
			// keyword LOG is not needed, but its presence is checked elsewhere
			strcpy( optimize.chVarFmt[optimize.nparm], "MAGNETIC FIELD TANGLED= %f LOG GAMMA= %f" );
			optimize.vparm[0][optimize.nparm] = (realnum)log10( Btangl_init );
			/* this is not varied */
			optimize.vparm[1][optimize.nparm] = (realnum)gamma_mag;
		}
		else
		{
			/* ordered field case */
			// keyword LOG is not needed, but its presence is checked elsewhere
			strcpy( optimize.chVarFmt[optimize.nparm], "MAGNETIC FIELD ORDERED= %f LOG ANGLE RADIANS= %f" );
			optimize.vparm[0][optimize.nparm] = (realnum)log10( Border_init );
			/* this is not varied */
			optimize.vparm[1][optimize.nparm] = (realnum)angle_wrt_los;
		}

		/* pointer to where to write */
		optimize.nvfpnt[optimize.nparm] = input.nRead;
		optimize.vincr[optimize.nparm] = 0.5;
		++optimize.nparm;
	}
	return;
}