/home66/gary/public_html/cloudy/c08_branch/source/magnetic.cpp

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2008 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"

/* 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.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 */
                        /* >>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(char *chCard )
{
        bool lgEOL;
        long int i;
        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( nMatch("ORDE" , chCard ) )
        {
                /* 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 */

                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( !nMatch("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
        {
                /* tangled field case */
                lgTangled = true;
                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 );
                        cdEXIT(EXIT_FAILURE);
                }
        }

        /* vary option */
        if( optimize.lgVarOn )
        {
                /* number of parameters */
                optimize.nvarxt[optimize.nparm] = 2;
                if( lgTangled )
                {
                        /* tangled field case */
                        strcpy( optimize.chVarFmt[optimize.nparm], "MAGNETIC FIELD TANGLED =%f 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 */
                        strcpy( optimize.chVarFmt[optimize.nparm], "MAGNETIC FIELD ORDERED =%f 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;
}

---