cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
magnetic.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2017 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*ParseMagnet parse magnetic field command */
4 /*Magnetic_init initialize magnetic field parameters */
5 /*Magnetic_reinit - reinitialized magnetic field at start of new iteration */
6 /*Magnetic_evaluate evaluate some parameters to do with magnetic field */
7 #include "cddefines.h"
8 #include "physconst.h"
9 #include "dense.h"
10 #include "doppvel.h"
11 #include "optimize.h"
12 #include "input.h"
13 #include "wind.h"
14 #include "magnetic.h"
15 #include "parser.h"
16 
18 
19 /* the initial magnetic field */
20 static double Btangl_init;
21 
22 /* this is logical var, set in zero, which says whether the magnetic field has been
23  * initialized */
24 static bool lgBinitialized;
25 
26 /* the current magnetic field */
27 static double Btangl_here;
28 
29 /* the initial parallel and tangential fields for ordered case */
30 static double Bpar_init, Btan_init;
31 
32 /* the current parallel and tangential fields for ordered case */
33 static double Bpar_here, Btan_here;
34 
35 /* this is the gamma power law index, default is 4. / 3. */
36 static double gamma_mag;
37 
38 /*Magnetic_evaluate evaluate some parameters to do with magnetic field */
40 {
41 
42  DEBUG_ENTRY( "Magnetic_evaluate()" );
43 
44  /* this flag set true if magnetic field is specified */
45  if( magnetic.lgB )
46  {
47  static double density_initial,
48  /* the square of the Alfven velocity at illuminated face */
49  v_A;
50 
51  /* does magnetic field need to be initialized for this iteration?
52  * flag is reset false at init of code, and at start of every iteration */
53  if( !lgBinitialized )
54  {
55  lgBinitialized = true;
56 
57  /* set initial tangled field */
59 
60  /* set initial ordered field */
61  /* mag field angle_wrt_los set when ordered field specified */
64 
65  /* XMassDensity was set above, safe to use this on first call */
66  density_initial = dense.xMassDensity;
67 
68  /* this is used to update tangential field */
69  v_A = POW2(Bpar_init) / (PI4 * density_initial );
70  }
71 
72  /* now update parameters in tangled field case */
73  /* magnetic pressure is a function of the local density, use gamma law */
74  Btangl_here = Btangl_init * pow(dense.xMassDensity/density_initial, gamma_mag/2.);
75 
76  /* ordered components of field - parallel field is always constant - find tangential component -
77  * but only in wind case */
78  if( !wind.lgStatic() )
79  {
80  /* N B - must preserve sign in this equation - will blow if product of wind speeds is equal to v_A) */
81  /* wind.windv*wind.windv0 == v_A should not occur since mag pressure goes to inf */
82  Btan_here = Btan_init * (POW2(wind.windv0) - v_A)/ (wind.windv*wind.windv0-v_A);
83  }
84 
85  /* magnetic pressure - sum of two fields - can have negative pressure (tension)
86  * is ordered field dominates */
87  magnetic.pressure = POW2(Btangl_here)/PI8 + (POW2(Btan_here) - POW2(Bpar_here)) / PI8;
88 
89  /* energy density - this is positive */
90  magnetic.energydensity = POW2(Btangl_here)/PI8 + (POW2(Btan_here) + POW2(Bpar_here)) / PI8;
91 
92  /* option for turbulence in equipartition with B field */
94  {
95  /* >>chng 05 jan 26, as per Robin Williams email,
96  * evaluate energydensity above, which is +ve, and use that for
97  * velocity here - had used pressure but could not evaluate when negative */
98  /* turbulent velocity is mag pres over density */
99  /* >>chng 06 apr 19, use DoppVel.Heiles_Troland_F */
100  /*DoppVel.TurbVel = (realnum)sqrt(2.*magnetic.energydensity/dense.xMassDensity);*/
103  /* >>chng 06 apr 19, do not double mag pressure, really count turbulence as pressure */
104  /* double magnetic pressure to account for ram pressure due to turbulence,
105  * which is not counted elsewhere
106  magnetic.pressure *= 2.;*/
107  }
108 
109  /* input parser made sure gamma != 1, default magnetic gamma is 4/3 */
110  magnetic.EnthalpyDensity = gamma_mag/(gamma_mag-1.) *
111  POW2(Btangl_here)/PI8 + (POW2(Btan_here) + POW2(Bpar_here)) / PI4;
112  }
113  else
114  {
115  magnetic.pressure = 0.;
116  magnetic.energydensity = 0.;
117  magnetic.EnthalpyDensity = 0.;
118  }
119  return;
120 }
121 
122 /*Magnetic_reinit - reinitialized magnetic field at start of new iteration */
123 void Magnetic_reinit(void)
124 {
125  DEBUG_ENTRY( "Magnetic_reinit()" );
126 
127  /* this says whether B has been initialized in this run */
128  lgBinitialized = false;
129  return;
130 }
131 
132 /* t_magnetic::zero initialize magnetic field parameters */
134 {
135 
136  DEBUG_ENTRY( "t_magnetic::zero()" );
137 
138  gamma_mag = 4./3.;
139  lgB = false;
140  /* this says whether B has been initialized in this run */
141  lgBinitialized = false;
142  /* the initial tangled and ordered fields */
143  Btangl_init = 0.;
144  Btangl_here = DBL_MAX;
145  pressure = DBL_MAX;
146  energydensity = DBL_MAX;
147  Bpar_init = 0.;
148  Btan_init = 0.;
149  Bpar_here = DBL_MAX;
150  Btan_here = DBL_MAX;
151  EnthalpyDensity = DBL_MAX;
152  return;
153 }
154 
155 /*ParseMagnet parse magnetic field command */
157 {
158  bool lgTangled;
159  double angle_wrt_los=-1. , Border_init=-1.;
160 
161  DEBUG_ENTRY( "ParseMagnet()" );
162 
163  /* flag saying B turned on */
164  magnetic.lgB = true;
165 
166  /* check whether ordered is specified - if not this is tangled */
167  if( p.nMatch("ORDE") )
168  {
169  /* ordered field case */
170  lgTangled = false;
171 
172  /* field is ordered, not isotropic - need field direction - spcified
173  * by angle away from beam of light - 0 => parallel to beam */
174 
175  /* this is the log of the ordered field strength */
176  Border_init = p.getNumberCheckAlwaysLog("ordered field");
177 
178  /* this is the angle (default in degrees) with respect to the line of sight */
179  angle_wrt_los = p.getNumberCheck("LOS angle");
180 
181  /* if radians is on the line then the number already is in radians -
182  * else convert to radians */
183  if( !p.nMatch("RADI") )
184  angle_wrt_los *= PI2 / 360.;
185 
186  /* now get initial components that we will work with */
187  Bpar_init = Border_init*cos(angle_wrt_los);
188  Btan_init = Border_init*sin(angle_wrt_los);
189 
190  }
191  else
192  {
193  /* tangled field case */
194  lgTangled = true;
195  /* this is the log of the tangled field strength */
196  Btangl_init = p.getNumberCheckAlwaysLog("tangled field");
197 
198  /* optional gamma for dependence on pressure */
199  gamma_mag = p.getNumberDefault("field gamma law", 4./3.);
200 
201  if( gamma_mag != 0. && gamma_mag <= 1. )
202  {
203  /* impossible value for gamma */
204  fprintf( ioQQQ,
205  " This value of gamma (%.3e) is impossible. Must have gamma = 0 or > 1.\n Sorry.\n",
206  gamma_mag );
208  }
209  }
210 
211  /* vary option */
212  if( optimize.lgVarOn )
213  {
214  /* number of parameters */
216  if( lgTangled )
217  {
218  /* tangled field case */
219  // keyword LOG is not needed, but its presence is checked elsewhere
220  strcpy( optimize.chVarFmt[optimize.nparm], "MAGNETIC FIELD TANGLED= %f LOG GAMMA= %f" );
222  /* this is not varied */
224  }
225  else
226  {
227  /* ordered field case */
228  // keyword LOG is not needed, but its presence is checked elsewhere
229  strcpy( optimize.chVarFmt[optimize.nparm], "MAGNETIC FIELD ORDERED= %f LOG ANGLE RADIANS= %f" );
230  optimize.vparm[0][optimize.nparm] = (realnum)log10( Border_init );
231  /* this is not varied */
232  optimize.vparm[1][optimize.nparm] = (realnum)angle_wrt_los;
233  }
234 
235  /* pointer to where to write */
237  optimize.vincr[optimize.nparm] = 0.5;
238  ++optimize.nparm;
239  }
240  return;
241 }
bool nMatch(const char *chKey) const
Definition: parser.h:140
static double gamma_mag
Definition: magnetic.cpp:36
t_input input
Definition: input.cpp:12
realnum Heiles_Troland_F
Definition: doppvel.h:37
long int nvfpnt[LIMPAR]
Definition: optimize.h:198
realnum windv0
Definition: wind.h:11
static double Btangl_init
Definition: magnetic.cpp:20
t_magnetic magnetic
Definition: magnetic.cpp:17
long int nRead
Definition: input.h:62
realnum TurbVel
Definition: doppvel.h:21
char chVarFmt[LIMPAR][FILENAME_PATH_LENGTH_2]
Definition: optimize.h:267
double pressure
Definition: magnetic.h:41
FILE * ioQQQ
Definition: cddefines.cpp:7
realnum vparm[LIMEXT][LIMPAR]
Definition: optimize.h:192
t_DoppVel DoppVel
Definition: doppvel.cpp:5
Definition: parser.h:42
bool lgTurbEquiMag
Definition: doppvel.h:46
bool lgVarOn
Definition: optimize.h:207
t_dense dense
Definition: global.cpp:15
static double Btan_init
Definition: magnetic.cpp:30
double EnthalpyDensity
Definition: magnetic.h:38
Wind wind
Definition: wind.cpp:5
double energydensity
Definition: magnetic.h:44
static double Bpar_here
Definition: magnetic.cpp:33
bool lgB
Definition: magnetic.h:35
static bool lgBinitialized
Definition: magnetic.cpp:24
#define POW2
Definition: cddefines.h:983
void ParseMagnet(Parser &p)
Definition: magnetic.cpp:156
long int nparm
Definition: optimize.h:204
float realnum
Definition: cddefines.h:124
#define EXIT_FAILURE
Definition: cddefines.h:168
#define cdEXIT(FAIL)
Definition: cddefines.h:484
t_optimize optimize
Definition: optimize.cpp:6
realnum vincr[LIMPAR]
Definition: optimize.h:195
void zero(void)
Definition: magnetic.cpp:133
static double Btangl_here
Definition: magnetic.cpp:27
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
realnum xMassDensity
Definition: dense.h:101
double getNumberCheckAlwaysLog(const char *chDesc)
Definition: parser.cpp:393
static double Btan_here
Definition: magnetic.cpp:33
double getNumberDefault(const char *chDesc, double fdef)
Definition: parser.cpp:367
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1217
double pow(double x, int i)
Definition: cddefines.h:782
bool lgStatic(void) const
Definition: wind.h:24
static double Bpar_init
Definition: magnetic.cpp:30
void Magnetic_reinit(void)
Definition: magnetic.cpp:123
long int nvarxt[LIMPAR]
Definition: optimize.h:198
realnum windv
Definition: wind.h:18
double getNumberCheck(const char *chDesc)
Definition: parser.cpp:358
void Magnetic_evaluate(void)
Definition: magnetic.cpp:39