VERB_code_2.2  2
 All Classes Namespaces Files Functions Variables Typedefs Macros Pages
variousFunctions.cpp
Go to the documentation of this file.
1 /**
2  * \file variousFunctions.cpp
3  * \brief Variouse functions description, see namespace VF for details
4  *
5  * \ref deprecated
6  *
7  * \author Developed by Yuri Shprits, Dmitry Subbotin, Kung-Chang Kim
8  */
9 
10 #include "variousFunctions.h"
11 #include <math.h>
12 #include <valarray>
13 #include <iostream>
14 #include <sstream>
15 
16 // namespace descriped in header file
17 namespace VF {
18 
19  /**
20  * Loss cone in rad.
21  * \f[ sin(\alpha)^2 = \frac{L^{-3}}{\sqrt{4 - \frac{3}{L}}} \f]
22  */
23  double alc(double L) {
24  return asin(pow(L, -1.5)*pow(4-3./L,-0.25));
25  }
26 
27  /** T(&alpha;) - function, related to bounce time <br>
28  * \f[ y = sin(\alpha) \f]
29  * \f[ T = 1.38 + 0.055 \times y^{1/3} - 0.32 \times y^{1/2} - 0.037 \times y^{2/3} - 0.394 \times y + 0.056 \times y^{4/3} \f]
30  * [Orlova and Shprits, On the bounce-averaging of scattering rates and the calculation of bounce period ,2011]
31  */
32  double T(double alpha) {
33  double y = sin(alpha);
34  return (1.38 + 0.055 * pow(y, 1.0/3) - 0.32 * pow(y, 1.0/2) - 0.037 * pow(y, 2.0/3) - 0.394 * y + 0.056 * pow(y, 4.0/3));
35  }
36 
37  /** This function used in Fokker-Plank equation. \f$ G = \cos(\alpha) \times \sin(\alpha) \f$ from Jacobian from %mu;,J to %alpha;, p and bounce avarage.
38  * G-function \f$ G = \cos(\alpha) \times \sin(\alpha) \times T(\alpha)\f$
39  */
40  double G(double alpha) {
41  // cos(x)*sin(x)*T(y);
42  return cos(alpha)*sin(alpha)*T(alpha);
43  }
44 
45  /**
46  * Bounce time, in days.
47  * \f[ \tau = \frac{4.0 \times L \times R_E \times mc^2}{pc \times c \times T(\alpha) \times ( 60 \times 60 \times 24 )}\f]
48  */
49  double bounce_time_new(double L, double pc, double alpha) {
50  return 4.0 * L * VC::RE * VC::mc2 / pc / VC::c_si * T(alpha) / 60 / 60 / 24;
51  }
52 
53  /**
54  * Y(&alpha;) function. More accurate approximation.
55  * \f[ y = sin(\alpha); T_0, T_1 = constant; \f]
56  * \f[ Y = 2(1 - y) \times T_0 + (T_0 - T_1) \times (y \log{y} + 2y - 2 \sqrt{y}) \f]
57  */
58  double Y(double alpha) {
59  double T0 = 1.3802;
60  double T1 = 0.7405;
61  double y = sin(alpha);
62  return 2.0*(1.0 - y)*T0 + (T0 - T1)*(y*log(y) + 2.0*y - 2.0*sqrt(y));
63  }
64 
65 
66  ////////////////////////////
67 
68  /**
69  * Computation of dipole magnetic field in Gauss.
70  * \f[ B = \frac{0.31}{L^3} \f]
71  */
72  double B (double Lparam) {
73  return 0.31/pow(Lparam, 3);
74  }
75 
76  /**
77  * ---------------------------------------------------------------------------------- <br>
78  * Electrostatic radial diffusion coefficient [see Brautigam and Albert(2000),JGR] <br>
79  * ---------------------------------------------------------------------------------- <br>
80  * This model seems to overestimate radial diffusion rate, so we need to scale down <br>
81  * by adjusting a parameter "scale_factor". <br>
82  * ---------------------------------------------------------------------------------- <br>
83  * (1) Inputs : alpha in radian, Ke in MeV, L and Kp in dimensionless unit <br>
84  * (2) Outputs : DLLe in 1/day or Re^2/day <br>
85  * ---------------------------------------------------------------------------------- <br>
86  */
87  double Dfe (double alpha, double Ke, double L, double Kp) {
88  double E_kp,DLLe;
89  double fac,dy,ty,factor,gamma,beta,wd;
90  double Re = 6.378*1e8; // Earth radius(cm)
91 
92  // Sometime E-field effect is too strong, so we need to scale down
93  // by adjusting this parameter until more accurate DLLe model is found
94  double scale_factor = 0.65;
95 
96  // exponenetial decay time in sec
97  double T0 = 2700.0;
98 
99  // convective electric field in mV/m, Kp=1 to 6;
100  if(Kp > 6) Kp=6;
101  E_kp = 0.26 * (Kp-1) + 0.1;
102 
103  // convert energy(MeV) to mu(MeV/G)
104  //pc = Ke*(Ke+2.0*VC::mc2); // pc in MeV
105  //mu = (pc*sin(alpha)*sin(alpha)) / (VF::B(L)*VC::mc2);
106 
107  // angular drift velocity wd (rad/s) averaged over a bounce in a dipole
108  // [see Schulz and Lanzerotti, 1974 and Schulz, 1991]
109 
110  fac = pow(sin(alpha),(3.0/4.0));
111  dy = (5.520692-2.357194*sin(alpha)+1.279385*fac) / 12.0;
112  ty = 1.380173-(0.639693*fac);
113  factor = dy/ty; // pitch angle dependency
114 
115  gamma = (Ke + VC::mc2)/VC::mc2;
116  beta = sqrt(1.0-1.0/(pow(gamma,2.0)));
117  //gamma = sqrt(1.0+2.0*VF::B(L)*mu/VC::mc2);
118  //wd = (1.475*1e-3*mu) / gamma / pow(L,2.0) * factor; // wd(rad/s),mu(MeV/G)
119  wd = 1.216 * 1e-3 * gamma * pow(beta,2.0) * L * factor;
120 
121  //wd = (1.475*1e-3*Ke) / VF::B(L) / pow(L,2.0) * factor; // wd(rad/s),Ke(MeV),B(G)
122 
123  // electrostatic radial diffusion coefficient in cm^2/sec
124  DLLe = (T0*pow(E_kp,2.0)*1e6) / (4.0*pow(VC::B_0,2)) / (1.0+pow((wd*T0/2.0),2.0)) * pow(L,6); // E_kp(mV/m),wd(1/s),T0(s),B_0(G)
125 
126  return scale_factor * DLLe * (60.0*60.0*24.0)/pow(Re,2); // DLLe in Re^2/day or 1/day
127  }
128 
129  /**
130  * -------------------------------------------------------------------- <br>
131  * Coulomb scattering electron lifetime [Lyons and Thorne, 1973, JGR] <br>
132  * -------------------------------------------------------------------- <br>
133  * (1) Inputs : L in dimensionless, energy in MeV <br>
134  * (2) Outputs: tau in days <br>
135  * -------------------------------------------------------------------- <br>
136  */
137  double Coulomb(double L,double energy) {
138  double n,tau,scale_factor;
139 
140  scale_factor = 1.0; // 1.8
141  // n is the equatorial electron density in #/cm3
142  n=1e3*pow((4.0/L),4.0); // original model
143 
144  // Carpenter and Anderson model, 1992, JGR
145  // If we apply the following model, tau model also should be changed.
146  //n=pow(10, (-0.3145*L + 3.9043));
147 
148  // electron lifetime due to coulomb scattering in sec
149  tau=( 3e8*pow((energy*1e3),1.5) ) / n;
150 
151  return scale_factor * tau/(60.0*60.0*24.0); // in days
152  }
153 
154  /**
155  * --------------------------------------------------------------------- <br>
156  * Empirical electron lifetime related to the resonant scattering by <br>
157  * plasmaspheric hiss, lightining-generated whistlers <br>
158  * and VLF transmitter signals, which is valid only in the slot region <br>
159  * between 2.5 and Lpp (plasmapause location) [see Su et al, 2010, JGR] <br>
160  * --------------------------------------------------------------------- <br>
161  * (1) Inputs : energy in MeV, Kp and L in dimensionless <br>
162  * (2) Outputs : tau in days <br>
163  * --------------------------------------------------------------------- <br>
164  * This model is valid up to Lpp, so class "PSD::SourcesAndLosses" checks if input L is larger than Lpp or not.
165  * For L<2.5 \f$ \tau = 10^99 \f$, for L>=2.5
166  * \f[ \tau = \frac{2}{Kp}^2 \times 10^{(L-3.5)^4}+0.4 \f]
167  */
168  double Precipitation(double Kp,double L,double energy) {
169  double n,tau;
170  if(L>=2.5) { // this model is valid up to Lpp
171  tau = pow((2.0/Kp),2.0) * pow(10.0,pow((L-3.5),4.0)+0.4); // in days
172  } else {
173  tau =1e99;
174  }
175 
176  return tau; // in days
177  }
178 
179  /**
180  * -------------------------------------------------------------------- <br>
181  * Subroutine to call various flux model of the outer boundary at L=7 <br>
182  * -------------------------------------------------------------------- <br>
183  * (1) Inputs : EL - energy(MeV) at L=7 <br>
184  * J_L7_function - name of boundary flux function <br>
185  * (2) Outputs : JL - flux at L=7 <br>
186  * -------------------------------------------------------------------- <br>
187  */
188  double Outer_Boundary_Flux_at_L7(double E7,string J_L7_function) {
189 
190  double J7;
191 
192  // boundary flux model at L=7
193  if (J_L7_function == "J_L7") {
194  J7 = VF::J_L7( E7,0.2,2.e3,1,7 );
195  } else if (J_L7_function == "J_L7_old") {
196  J7 = VF::J_L7_old( E7 );
197  } else if (J_L7_function == "J_L7_corrected") {
198  J7 = VF::J_L7_corrected( E7 );
199  } else if (J_L7_function == "J_L7_corrected_0") {
200  J7 = VF::J_L7_corrected_0( E7 );
201  } else {
202  throw error_msg("INITIAL_PSD", "Unknown boundary function J_L7: %s. Should be 'J_L7' or 'J_L7_corrected'.", J_L7_function.c_str());
203  }
204 
205  return J7;
206  }
207 
208  /**
209  * Computation of moumentum from Kinetic energy.
210  * \f[ pc = (\sqrt{ { \frac{K}{mc^2} + 1}^2 } - 1) \times mc^2 \f]
211  */
212  double pfunc(double K) {//, double mc2) {
213  return sqrt( pow( K / VC::mc2 + 1, 2) - 1) * VC::mc2;
214  }
215 
216  /**
217  * Computation of Kinetic energy from given momentum.
218  * \f[ K = (\sqrt{ 1 + {\frac{pc}{mc^2}}^2 } - 1) \times mc^2 \f]
219  */
220  double Kfunc(double pc) {// , double mc2) {
221  return ( sqrt( 1.0 + pow( pc / VC::mc2, 2) ) - 1) * VC::mc2;
222  }
223 
224  /**
225  * Interpolate outer boundary.
226  * This function used by VF::J_L7_corrected, VF::J_L7_corrected_0 ... The output is interpolation of the initial spectrum for specific point.
227  */
228  double J_L7 (double K, double x1, double y1, double x2, double y2) {
229  //double x1=0.2;
230  //double y1=2.e3;
231  //double x2=1.;
232  //double y2=7.;
233  // x1=0.2
234  // y1=2.e3
235  // x2=1.
236  // y2=7.
237  double bcoef = log(y1/y2)/(x2-x1);
238  double acoef = y2*exp(bcoef*x2);
239  return acoef*exp(-bcoef*K);
240  }
241 
242  /**
243  * Specifying outer boundary, more accurate function.
244  */
245  double J_L7_corrected (double K) {
246  double x[10] = {1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 2e-1, 1e0, 3e0, 1e1, 2e1};
247  double y[10] = {3e9, 2e7, 9e6, 3e5, 1e4, 2e3, 7e0, 3e-3, 1e-3, 1e-3};
248 
249  unsigned i=0;
250 
251  while (x[++i] < K && i < sizeof(x)/8-1);
252  double result = J_L7 (K, x[i-1], y[i-1], x[i], y[i]);
253  return result;
254  }
255 
256  /**
257  * Specifying outer boundary, more accurate function, zero at infinity
258  */
259  double J_L7_corrected_0 (double K) {
260  double x[10] = {1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 2e-1, 1e0, 3e0, 1e1, 1e2};
261  double y[10] = {3e9, 2e7, 9e6, 3e5, 1e4, 2e3, 7e0, 3e-3, 1e-3, VC::zero_f};
262 
263  unsigned i=0;
264 
265  while (x[++i] < K && i < sizeof(x)/8-1);
266  double result = J_L7 (K, x[i-1], y[i-1], x[i], y[i]);
267  return result;
268  }
269 
270  /** old version \f$ J = 8222.6 \times \exp(-7.068 \times K); \f$ */
271  double J_L7_old (double K) {
272  double flux;
273 
274  flux=8222.6*exp(-7.068*K); // K(kinetic energy) in MeV
275 
276  return flux;
277  }
278 
279  /**
280  * Calculating pc from L, &mu;, alpha.
281  * \f[ pc = \frac{\sqrt{ \mu \times 2 \times B(L) \times mc^2 }}{\sin(\alpha)} \f]
282  */
283  double mu2pc (double L, double mu, double alpha) {
284  return sqrt( mu * 2 * B(L) * VC::mc2 ) / sin(alpha);
285  //return sqrt(mu*B(L)); // alpha == 90'
286  }
287 
288  /**
289  * Calculating &mu; from L, pc, alpha.
290  * \f[ \mu = \frac{pc^2 \times \sin(\alpha)^2}{B(L) \times 2 \times mc^2 } \f]
291  */
292  double pc2mu (double L, double pc, double alpha) {
293  return pow(pc, 2) * pow(sin(alpha), 2) / (B(L) * 2 * VC::mc2);
294  }
295 
296  /**
297  * Caltulating J*c from L, pc, &alpha;.
298  * \warning Determined up to a constant.
299  * \f[ J = 2 \times pc \times L \times Y(\alpha) \f]
300  */
301  double Jc_calc(double L, double pc, double Alpha) {
302  return 2*pc*L*Y(Alpha);
303  }
304 
305  // double Jc2Alpha (double L, double pc, double Jc) {
306  // double a = 1;
307  // return sqrt(acos(Jc/2/0.7405/a/L/pc));
308  // }
309 
310  /**
311  * Getting f on specific energy. Linear interpolation of logarithm.
312  */
313  double f_interp(double E, double f1, double E1, double f2, double E2) {
314  double a, b, f_epc;
315  if (E < E1 || E > E2) {
316  return -1;
317  }
318  if (f1 <= 0 || f2 <= 0) return 0;
319  b = log(f2/f1)/(E1 - E2);
320  a = f1/exp(-b*E1);
321  f_epc = a * exp(-b*E);
322  return f_epc;
323  }
324 
325  /**
326  * Caltulating &mu; from L, pc, &alpha;.
327  * \f[ \mu = \frac{pc^2}{mc^2} \times \frac{ \sin(\alpha)^2 } { 2 \times B(L) }; \f]
328  */
329  double mu_calc(double L, double pc, double Alpha) {
330  return pow(pc,2)/VC::mc2 * pow(sin(Alpha),2) / (2 * B(L));
331  }
332 
333  /**
334  * Chorus density model? (deprecated?)
335  * \f$ D = 10^{-0.3145L + 3.9043} \f$
336  */
337  double density(double L) {
338  return pow(10, (-0.3145*L + 3.9043));
339  }
340 
341  /*
342  * Chorus density model?
343  */
344  /*Matrix3D<double> density(Matrix3D<double> &L) {
345  Matrix3D<double> tmp;
346  tmp.AllocateMemory(L.size_x, L.size_y, L.size_z);
347  for (int x = 0; x < tmp.size_x; x++)
348  for (int y = 0; y < tmp.size_y; y++)
349  for (int z = 0; z < tmp.size_z; z++)
350  tmp[x][y][z] = pow(10, (-0.3145)*L[x][y][z] + 3.9043);
351  return tmp;
352  }*/
353 
354  /// Return maximum.
355  double max(double v1, double v2) {
356  return (v1>v2)?v1:v2;
357  }
358 
359  /// Convert double to string
360  string dtostr(double n) {
361  std::ostringstream o;
362  o << n;
363  return o.str();
364  }
365 
366  /// Small routine for checking if it is right iteration for the next output
367  bool check_time(int iter, int d_iter) {
368  if ((iter % d_iter) == 0) {
369  return true;
370  }
371  return false;
372  }
373 }