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