cloudy trunk
Loading...
Searching...
No Matches
magnetic.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 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 */
20static double Btangl_init;
21
22/* this is logical var, set in zero, which says whether the magnetic field has been
23 * initialized */
24static bool lgBinitialized;
25
26/* the current magnetic field */
27static double Btangl_here;
28
29/* the initial parallel and tangential fields for ordered case */
30static double Bpar_init, Btan_init;
31
32/* the current parallel and tangential fields for ordered case */
33static double Bpar_here, Btan_here;
34
35/* this is the gamma power law index, default is 4. / 3. */
36static 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 */
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 */
93 if( DoppVel.lgTurbEquiMag )
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);*/
101 DoppVel.TurbVel = (realnum)sqrt(6.*magnetic.energydensity/dense.xMassDensity/
102 DoppVel.Heiles_Troland_F);
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.) *
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 */
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/* Magnetic_init initialize magnetic field parameters */
134{
135
136 DEBUG_ENTRY( "Magnetic_init()" );
137
138 gamma_mag = 4./3.;
139 magnetic.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 magnetic.pressure = DBL_MAX;
146 magnetic.energydensity = DBL_MAX;
147 Bpar_init = 0.;
148 Btan_init = 0.;
149 Bpar_here = DBL_MAX;
150 Btan_here = DBL_MAX;
151 magnetic.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 */
215 optimize.nvarxt[optimize.nparm] = 2;
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" );
221 optimize.vparm[0][optimize.nparm] = (realnum)log10( Btangl_init );
222 /* this is not varied */
223 optimize.vparm[1][optimize.nparm] = (realnum)gamma_mag;
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 */
236 optimize.nvfpnt[optimize.nparm] = input.nRead;
237 optimize.vincr[optimize.nparm] = 0.5;
238 ++optimize.nparm;
239 }
240 return;
241}
FILE * ioQQQ
Definition cddefines.cpp:7
#define POW2
Definition cddefines.h:929
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
bool nMatch(const char *chKey) const
Definition parser.h:135
double getNumberCheckAlwaysLog(const char *chDesc)
Definition parser.cpp:308
double getNumberCheck(const char *chDesc)
Definition parser.cpp:273
double getNumberDefault(const char *chDesc, double fdef)
Definition parser.cpp:282
t_dense dense
Definition dense.cpp:24
t_DoppVel DoppVel
Definition doppvel.cpp:5
t_input input
Definition input.cpp:12
static double Btan_here
Definition magnetic.cpp:33
t_magnetic magnetic
Definition magnetic.cpp:17
void Magnetic_init(void)
Definition magnetic.cpp:133
static double Btangl_here
Definition magnetic.cpp:27
void Magnetic_evaluate(void)
Definition magnetic.cpp:39
static double gamma_mag
Definition magnetic.cpp:36
static double Btangl_init
Definition magnetic.cpp:20
static bool lgBinitialized
Definition magnetic.cpp:24
static double Bpar_here
Definition magnetic.cpp:33
static double Btan_init
Definition magnetic.cpp:30
void Magnetic_reinit(void)
Definition magnetic.cpp:123
static double Bpar_init
Definition magnetic.cpp:30
void ParseMagnet(Parser &p)
Definition magnetic.cpp:156
t_optimize optimize
Definition optimize.cpp:5
UNUSED const double PI4
Definition physconst.h:35
UNUSED const double PI2
Definition physconst.h:32
UNUSED const double PI8
Definition physconst.h:38
Wind wind
Definition wind.cpp:5