VERB_code_2.3
variousFunctions.cpp
Go to the documentation of this file.
1 
10 #include "variousFunctions.h"
11 #include <math.h>
12 #include <valarray>
13 #include <iostream>
14 #include <sstream>
15 
16 // namespace described in header file
17 namespace VF {
18 
23  double alc(double L) {
24  return asin(pow(L, -1.5)*pow(4-3./L,-0.25));
25  }
26 
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 
40  double G(double alpha) {
41  // cos(x)*sin(x)*T(y);
42  return cos(alpha)*sin(alpha)*T(alpha);
43  }
44 
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 
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 
67 
72  double B (double Lparam) {
73  return 0.31/pow(Lparam, 3);
74  }
75 
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 
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 
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 
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 
212  double pfunc(double K) {//, double mc2) {
213  return sqrt( pow( K / VC::mc2 + 1, 2) - 1) * VC::mc2;
214  }
215 
220  double Kfunc(double pc) {// , double mc2) {
221  return ( sqrt( 1.0 + pow( pc / VC::mc2, 2) ) - 1) * VC::mc2;
222  }
223 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
355  double max(double v1, double v2) {
356  return (v1>v2)?v1:v2;
357  }
358 
360  string dtostr(double n) {
361  std::ostringstream o;
362  o << n;
363  return o.str();
364  }
365 
367  bool check_time(int iter, int d_iter) {
368  if ((iter % d_iter) == 0) {
369  return true;
370  }
371  return false;
372  }
373 }
double max(double v1, double v2)
Return maximum.
double Precipitation(double Kp, double L, double energy)
Empirical electron lifetime.
double f_interp(double E, double f1, double E1, double f2, double E2)
Interpolation f on specific energy.
double Jc_calc(double L, double pc, double Alpha)
Caltulating J*c.
static const double c_si
speed of light, km/s
double mu2pc(double L, double mu, double alpha)
Convert μ to pc.
double Outer_Boundary_Flux_at_L7(double E7, string J_L7_function)
Call various flux model of the outer boundary.
static const double exp
exponential
string dtostr(double n)
Convert double to string.
double J_L7_corrected(double K)
Outer boundary spectrum (accurate)
double pfunc(double K)
Computation of moumentum from Kinetic energy .
static const double mc2
m*c^2 = 0.511 MeV
Various functions.
double f1(double lambda)
double Coulomb(double L, double energy)
Coulomb scattering electron lifetime.
double J_L7(double K, double x1, double y1, double x2, double y2)
Interpolate spectrum.
double bounce_time_new(double L, double pc, double alpha)
bounce time
double Dfe(double alpha, double Ke, double L, double Kp)
Electrostatic radial diffusion coefficient.
Various functions definition, see namespace VF for details.
double pc2mu(double L, double pc, double alpha)
Convert pc to μ.
double mu_calc(double L, double pc, double Alpha)
Caltulating μ.
double B(double Lparam)
Dipole magnetic field.
double T(double alpha)
Function for bounce time.
double alc(double L)
Loss cone calculations.
double J_L7_corrected_0(double K)
Outer boundary spectrum (accurate with 0 end the end)
double Y(double alpha)
Y(α) function.
Error message - stack of single_errors.
Definition: error.h:54
bool check_time(int iter, int d_iter)
Small routine for checking if it is right iteration for the next output.
double density(double L)
Chorus density model.
static const double zero_f
Minimum value for PSD.
static const double B_0
magnetic field at the equatorial plane on the Earth surface, Gauss
double J_L7_old(double K)
Outer boundary spectrum (old)
static const double RE
Radius fo the Earch in km.
double Kfunc(double pc)
Computation of Kinetic energy from given momentum .
double G(double alpha)
G-function.