cloudy  trunk
 All Data Structures 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-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 */
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 */
56 
57  /* set initial ordered field */
58  /* mag field angle_wrt_los set when ordered field specified */
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 */
85 
86  /* energy density - this is positive */
88 
89  /* option for turbulence in equipartition with B field */
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);*/
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 */
109  }
110  else
111  {
112  magnetic.pressure = 0.;
113  magnetic.energydensity = 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 */
225  if( lgTangled )
226  {
227  /* tangled field case */
228  strcpy( optimize.chVarFmt[optimize.nparm], "MAGNETIC FIELD TANGLED =%f GAMMA= %f" );
230  /* this is not varied */
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 */
244  optimize.vincr[optimize.nparm] = 0.5;
245  ++optimize.nparm;
246  }
247  return;
248 }

Generated for cloudy by doxygen 1.8.3.1