OSHUN  beta
Arbitrary Order Spherical-Harmonic 1D-3P Vlasov-Fokker-Planck-Maxwell code
formulary.cpp
Go to the documentation of this file.
1 
15 //-------------------------------------------------------------------
16 
17 // Standard libraries
18 #include <iostream>
19 #include <vector>
20 #include <valarray>
21 #include <complex>
22 #include <fstream>
23 #include <iomanip>
24 #include <cstdlib>
25 #include <sstream>
26 #include <string>
27 #include <cstring>
28 #include <math.h>
29 #include <map>
30 
31 // My libraries
32 #include "lib-array.h"
33 
34 // Declerations
35 #include "input.h"
36 #include "formulary.h"
37 
38 
39 //-------------------------------------------------------------------
40 //-------------------------------------------------------------------
41 // List all the units as in "label, conversion_factor"
42 //-------------------------------------------------------------------
43 Formulary::Formulary() : n(Input::List().density_np),
44  wp(5.64*1.0e+4*sqrt(n)),
45  skindepth(cL/wp),
46  B0(-1.0*wp*me/qe),
47  T0(pow(Input::List().pth_ref,2.0)*keVnorm*1e3),
48  Zeta(Input::List().hydrocharge)
49  // ref_nuei(sqrt(2.0/pi)*LOGei(1.0,Input::List().pth_ref,Zeta)/exp(LOGei(1.0,Input::List().pth_ref,Zeta)) * wp)
50  {
51  // Physically: label * d = const in any system
52 
53  // velocity - c
54  D["Velocity"] = units("c",1.0);
55  D["Velocity_cgs"] = units("cm/sec", c );
56  D["Velocity_si"] = units("m/sec", c * 0.01 );
57 
58  D["v"] = units("c",1.0);
59  D["v_cgs"] = units("cm/sec", c );
60  D["v_si"] = units("m/sec", c * 0.01 );
61 
62  D["vx"] = units("c",1.0);
63  D["vx_cgs"] = units("cm/sec", c );
64  D["vx_si"] = units("m/sec", c * 0.01 );
65 
66  D["vy"] = units("c",1.0);
67  D["vy_cgs"] = units("cm/sec", c );
68  D["vy_si"] = units("m/sec", c * 0.01 );
69 
70  D["vz"] = units("c",1.0);
71  D["vz_cgs"] = units("cm/sec", c );
72  D["vz_si"] = units("m/sec", c * 0.01 );
73 
74  // time - 1/wp
75  D["Time"] = units("1/wp",1.0);
76  D["Time_cgs"] = units("sec", 1.0/wp);
77  D["Time_si"] = units("sec", 1.0/wp);
78  D["Time_fs"] = units("fsec", 1.0/wp*1.0e+15);
79  D["Time_ps"] = units("psec", 1.0/wp*1.0e+12);
80 
81  D["t"] = units("1/wp",1.0);
82  D["t_cgs"] = units("sec", 1.0/wp);
83  D["t_si"] = units("sec", 1.0/wp);
84 
85  // space - c/wp
86  D["Space"] = units("c/wp",1.0);
87  D["Space_cgs"] = units("cm", c/wp);
88  D["Space_si"] = units("m", 0.01 * c/wp); // 1 cm = 0.01 m
89 
90  D["x"] = units("c/wp",1.0);
91  D["x_cgs"] = units("cm", c/wp);
92  D["x_si"] = units("m", 0.01 * c/wp); // 1 cm = 0.01 m
93 
94  D["y"] = units("c/wp",1.0);
95  D["y_cgs"] = units("cm", c/wp);
96  D["y_si"] = units("m", 0.01 * c/wp); // 1 cm = 0.01 m
97 
98  D["z"] = units("c/wp",1.0);
99  D["z_cgs"] = units("cm", c/wp);
100  D["z_si"] = units("m", 0.01 * c/wp); // 1 cm = 0.01 m
101 
102  // density - n
103  D["Density"] = units("n",1.0);
104  D["Density_cgs"] = units("cm^-3", n);
105  D["Density_si"] = units("m^-3", 1.0e+6 * n); // cm^-3 = 10^6 m^-3
106 
107  D["n"] = units("n",1.0);
108  D["n_cgs"] = units("cm^-3", n);
109  D["n_si"] = units("m^-3", 1.0e+6 * n); // cm^-3 = 10^6 m^-3
110 
111  // charge - e
112  D["Charge"] = units("e",1.0);
113  D["Charge_cgs"] = units("Fr", e);
114  D["Charge_si"] = units("C", (10.0/c) * e); // 1 Fr = (10/c) * C
115 
116  D["q"] = units("e",1.0);
117  D["q_cgs"] = units("Fr", e);
118  D["q_si"] = units("C", (10.0/c) * e); // 1 Fr = (10/c) * C
119 
120  // momentum - mc
121  D["Momentum"] = units("mc",1.0);
122  D["Momentum_cgs"] = units("gr*cm/sec", m*c );
123  D["Momentum_si"] = units("kg*m/sec", m*c * 1.0e-5 );
124 
125  D["p"] = units("mc",1.0);
126  D["p_cgs"] = units("gr*cm/sec", m*c );
127  D["p_si"] = units("kg*m/sec", m*c * 1.0e-5 );
128 
129  D["px"] = units("mc",1.0);
130  D["px_cgs"] = units("gr*cm/sec", m*c );
131  D["px_si"] = units("kg*m/sec", m*c * 1.0e-5 );
132 
133  D["py"] = units("mc",1.0);
134  D["py_cgs"] = units("gr*cm/sec", m*c );
135  D["py_si"] = units("kg*m/sec", m*c * 1.0e-5 );
136 
137  D["pz"] = units("mc",1.0);
138  D["pz_cgs"] = units("gr*cm/sec", m*c );
139  D["pz_si"] = units("kg*m/sec", m*c * 1.0e-5 );
140 //-----
141 
142  // energy - mc^2
143  D["Energy"] = units("mc^2",1.0);
144  D["Energy_cgs"] = units("erg", m*pow(c,2) );
145  D["Energy_si"] = units("J", 1.0e-7 * m*pow(c,2) ); // erg = 10^-7 J
146  D["Energy_eV"] = units("eV", 1.0e-7 * m*pow(c,2)/((10.0/c) * e) );
147 
148  D["T"] = units("mc^2",1.0);
149  D["T_cgs"] = units("erg", m*pow(c,2) );
150  D["T_si"] = units("J", 1.0e-7 * m*pow(c,2) ); // erg = 10^-7 J
151  D["T_eV"] = units("eV", 1.0e-7 * m*pow(c,2)/((10.0/c) * e) );
152 
153  // current - ewp
154  D["Current"] = units("ewp",1.0);
155  D["Current_cgs"] = units("Fr/s", e*wp );
156  D["Current_si"] = units("A", 10.0/c * e*wp ); // Fr/s = 10/c A
157 
158  D["Jx"] = units("ewp",1.0);
159  D["Jx_cgs"] = units("Fr/s", e*wp );
160  D["Jx_si"] = units("A", 10.0/c * e*wp ); // Fr/s = 10/c A
161 
162  D["Jy"] = units("ewp",1.0);
163  D["Jy_cgs"] = units("Fr/s", e*wp );
164  D["Jy_si"] = units("A", 10.0/c * e*wp ); // Fr/s = 10/c A
165 
166  D["Jz"] = units("ewp",1.0);
167  D["Jz_cgs"] = units("Fr/s", e*wp );
168  D["Jz_si"] = units("A", 10.0/c * e*wp ); // Fr/s = 10/c A
169 
170  // pressure - n*mc^2
171  D["Pressure"] = units("n*mc^2",1.0);
172  D["Pressure_cgs"] = units("Ba", n*m*c*c );
173  D["Pressure_si"] = units("Pa", 0.1 * n*m*c*c ); // Ba = 0.1 Pa
174  D["Pressure_Mbar"] = units("Mbar", 1.0e-12 * n*m*c*c ); // Ba = 10^-12Mbar
175 
176  D["P"] = units("n*mc^2",1.0);
177  D["P_cgs"] = units("Ba", n*m*c*c );
178  D["P_si"] = units("Pa", 0.1 * n*m*c*c ); // Ba = 0.1 Pa
179  D["P_Mbar"] = units("Mbar", 1.0e-12 * n*m*c*c ); // Ba = 10^-12Mbar
180 //-----
181 
182  // Electric field - mcwp/e
183  D["Efield"] = units("mcwp/e",1.0);
184  D["Efield_cgs"] = units("statV/cm", m*c*wp/e);
185  D["Efield_si"] = units("V/m", c*1.0e-6 * m*c*wp/e); // statV/m = 10^-6*c * V/m
186 
187  D["Ex"] = units("mcwp/e",1.0);
188  D["Ex_cgs"] = units("statV/cm", m*c*wp/e);
189  D["Ex_si"] = units("V/m", c*1.0e-6 * m*c*wp/e); // statV/m = 10^-6*c * V/m
190 
191  D["Ey"] = units("mcwp/e",1.0);
192  D["Ey_cgs"] = units("statV/cm", m*c*wp/e);
193  D["Ey_si"] = units("V/m", c*1.0e-6 * m*c*wp/e); // statV/m = 10^-6*c * V/m
194 
195  D["Ez"] = units("mcwp/e",1.0);
196  D["Ez_cgs"] = units("statV/cm", m*c*wp/e);
197  D["Ez_si"] = units("V/m", c*1.0e-6 * m*c*wp/e); // statV/m = 10^-6*c * V/m
198 
199  // Magnetic field - mcwp/e
200  D["Bfield"] = units("mcwp/e",1.0);
201  D["Bfield_cgs"] = units("Gauss", m*c*wp/e);
202  D["Bfield_si"] = units("Tesla", 1.0e-4 * m*c*wp/e); // Gauss = 10^-4 Tesla
203 
204  D["Bx"] = units("mcwp/e",1.0);
205  D["Bx_cgs"] = units("Gauss", m*c*wp/e);
206  D["Bx_si"] = units("Tesla", 1.0e-4 * m*c*wp/e); // Gauss = 10^-4 Tesla
207 
208  D["By"] = units("mcwp/e",1.0);
209  D["By_cgs"] = units("Gauss", m*c*wp/e);
210  D["By_si"] = units("Tesla", 1.0e-4 * m*c*wp/e); // Gauss = 10^-4 Tesla
211 
212  D["Bz"] = units("mcwp/e",1.0);
213  D["Bz_cgs"] = units("Gauss", m*c*wp/e);
214  D["Bz_si"] = units("Tesla", 1.0e-4 * m*c*wp/e); // Gauss = 10^-4 Tesla
215 
216  // force - mcwp
217  D["Force"] = units("mcwp",1.0);
218  D["Force_cgs"] = units("dyne", m*c*wp);
219  D["Force_si"] = units("N", 1.0e-5 * m*c*wp ); // dyne = 10^-5 N
220 
221 }
222 
223 //-------------------------------------------------------------------
224 // Calculate the Coulomb logarithm
225 //-------------------------------------------------------------------
226 double Formulary::LOGee(double ne, double Te) {
227 // ne = density/np, Te = energy/mc^2
228 // Note: we assume nonrelativistic distribution functions
229  if (ne < nmin) return 2.0;
230 
231  // Te /= (3.0*ne);
232  Te *= Units("Energy","eV").d; // Temperature in eV
233  ne *= Units("Density","cgs").d;
234 
235  Te = log(Te);
236  ne = log(ne);
237 
238  return max(2.0,23.5 - 0.5*ne + 1.25*Te - sqrt(0.00001+0.0625*(Te-2.0)*(Te-2.0)));
239 }
240 
241 double Formulary::LOGee(double ne, string un, double Te, string uT) {
242 // where "un" are the units of ne, and uT of Te
243  Te /= Units("Energy",uT).d;
244  ne /= Units("Density",un).d;
245 
246  return LOGee(ne,Te);
247 }
248 
249 double Formulary::LOGei(double ne, double Te, double Z) {
250 // ne = density/np, Te = energy/mc^2
251 // Note: we assume nonrelativistic distribution functions
252  if (ne < nmin) return 2.0;
253 
254  Te *= Units("Energy","eV").d; // Temperature in eV
255  ne *= Units("Density","cgs").d;
256 
257  if ( Te > 10.0*Z*Z ) {
258  return max(2.0, 24.0 - 0.5*log(ne) + log(Te));
259  }
260  return max(2.0, 23.0 - 0.5*log(ne) + 1.5*log(Te) - log(Z));
261 }
262 
263 double Formulary::LOGei(double ne, string un, double Te, string uT, double Z) {
264 // where "un" are the units of ne, and uT of Te
265  Te /= Units("Energy",uT).d; // Temperature in eV
266  ne /= Units("Density",un).d;
267 
268  return LOGei(ne,Te,Z);
269 }
270 
271 double Formulary::LOGii(double m1, double Z1, double n1, double T1,
272  double m2, double Z2, double n2, double T2) {
273 // ne = density/np, Te = energy/mc^2
274 // Note: we assume nonrelativistic distribution functions
275  double mp=1;
276  if (n1 < nmin || n2 < nmin) return 2.0;
277 
278  T1 *= Units("Energy","eV").d; // Temperature in eV
279  T2 *= Units("Energy","eV").d; // Temperature in eV
280  n1 *= Units("Density","cgs").d;
281  n2 *= Units("Density","cgs").d;
282 
283  return max(2.0, 23.0 - log(Z1*Z2*(m1/mp+m2/mp)/(m1/mp*T2+m2/mp*T1)*sqrt(n1*Z1*Z1/T1+n2*Z2*Z2/T2)));
284 }
285 
286 //-------------------------------------------------------------------
287 // Simple formulary expresions
288 //-------------------------------------------------------------------
289 
290 // Thermal velocity: NRL Formulary 2009, p.29
291 double Formulary::vth(double Te){
292  Te *= Units("Energy","eV").d;
293  return sqrt(Te)*4.19e+7/c;
294 }
295 double Formulary::vth(double Te, string uT){
296  Te /= Units("Energy", uT).d;
297  return vth(Te);
298 }
299 // Electron-ion self-collision time: NRL Formulary 2009, p.37
300 double Formulary::Tau_i(double ne, double Te, double Zeta){
301  return 2.09e+7*(pow(sqrt(Te/Zeta*Units("Energy","eV").d),3)/
302  (ne*Units("Density","cgs").d*LOGei(ne,Te,Zeta))) /sqrt(me_over_mp);
303 }
304 
305 // Electron self-collision time: NRL Formulary 2009, p.37
306 double Formulary::Tau_e(double ne, double Te){
307  return 3.44e+5*(pow(sqrt(Te*Units("Energy","eV").d),3)/
308  (ne*Units("Density","cgs").d*LOGee(ne,Te))) ; //
310 }
311 double Formulary::Tau_e(double ne, string un, double Te, string uT){
312  Te /= Units("Energy",uT).d;
313  ne /= Units("Density",un).d;
314  return Tau_e(ne,Te);
315 }
316 
317 // Electron mean free path
318 double Formulary::MFP(double ne, double Te){
319  return vth(Te)*Tau_e(ne,Te);
320 }
321 double Formulary::MFP(double ne, string un, double Te, string uT){
322  Te /= Units("Energy",uT).d;
323  ne /= Units("Density",un).d;
324  return MFP(ne,Te);
325 }
326 //-------------------------------------------------------------------
327 
328 Formulary& formulary() {
329  static Formulary f;
330  return f;
331 }
332 //-------------------------------------------------------------------
333 //*******************************************************************
formulary.h
Plasma Formulary and Units - Declarations.
Formulary::Tau_e
double Tau_e(double ne, double Te)
Definition: formulary.cpp:306
Formulary::D
map< string, units > D
Definition: formulary.h:124
Formulary::m
static constexpr double m
Definition: formulary.h:121
input.h
Input reader - Declarations.
lib-array.h
Underlying data structures.
Formulary::wp
double wp
Definition: formulary.h:100
Formulary::Zeta
double Zeta
Definition: formulary.h:104
units::d
double d
Definition: formulary.h:52
Formulary::LOGee
double LOGee(double ne, double Te)
Definition: formulary.cpp:226
Formulary::Tau_i
double Tau_i(double ne, double Te, double Zeta)
Definition: formulary.cpp:300
Formulary::e
static constexpr double e
Definition: formulary.h:120
Formulary::LOGei
double LOGei(double ne, double Te, double Z)
Definition: formulary.cpp:249
Formulary::c
static constexpr double c
Definition: formulary.h:119
Formulary::MFP
double MFP(double ne, double Te)
Definition: formulary.cpp:318
Formulary::nmin
static constexpr double nmin
Definition: formulary.h:126
Formulary::n
double n
Definition: formulary.h:99
formulary
Formulary & formulary()
Definition: formulary.cpp:328
Formulary::LOGii
double LOGii(double m1, double Z1, double n1, double T1, double m2, double Z2, double n2, double T2)
Definition: formulary.cpp:271
Formulary::Units
units Units(string key)
Definition: formulary.h:74
Input::List
Input_List & List()
Definition: input.cpp:1585
Formulary::vth
double vth(double Te)
Definition: formulary.cpp:291
Formulary::me_over_mp
static constexpr double me_over_mp
Definition: formulary.h:116
Input
Definition: input.h:18