23 double alc(
double L) {
24 return asin(pow(L, -1.5)*pow(4-3./L,-0.25));
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));
40 double G(
double alpha) {
42 return cos(alpha)*sin(alpha)*
T(alpha);
58 double Y(
double alpha) {
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));
72 double B (
double Lparam) {
73 return 0.31/pow(Lparam, 3);
87 double Dfe (
double alpha,
double Ke,
double L,
double Kp) {
89 double fac,dy,ty,factor,gamma,beta,wd;
90 double Re = 6.378*1e8;
94 double scale_factor = 0.65;
101 E_kp = 0.26 * (Kp-1) + 0.1;
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);
116 beta = sqrt(1.0-1.0/(pow(gamma,2.0)));
119 wd = 1.216 * 1e-3 * gamma * pow(beta,2.0) * L * factor;
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);
126 return scale_factor * DLLe * (60.0*60.0*24.0)/pow(Re,2);
138 double n,tau,scale_factor;
142 n=1e3*pow((4.0/L),4.0);
149 tau=( 3e8*pow((energy*1e3),1.5) ) / n;
151 return scale_factor * tau/(60.0*60.0*24.0);
171 tau = pow((2.0/Kp),2.0) * pow(10.0,pow((L-3.5),4.0)+0.4);
193 if (J_L7_function ==
"J_L7") {
195 }
else if (J_L7_function ==
"J_L7_old") {
197 }
else if (J_L7_function ==
"J_L7_corrected") {
199 }
else if (J_L7_function ==
"J_L7_corrected_0") {
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());
228 double J_L7 (
double K,
double x1,
double y1,
double x2,
double y2) {
237 double bcoef = log(y1/y2)/(x2-x1);
238 double acoef = y2*
exp(bcoef*x2);
239 return acoef*
exp(-bcoef*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};
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]);
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};
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]);
274 flux=8222.6*
exp(-7.068*K);
283 double mu2pc (
double L,
double mu,
double alpha) {
284 return sqrt( mu * 2 *
B(L) *
VC::mc2 ) / sin(alpha);
292 double pc2mu (
double L,
double pc,
double alpha) {
293 return pow(pc, 2) * pow(sin(alpha), 2) / (
B(L) * 2 *
VC::mc2);
301 double Jc_calc(
double L,
double pc,
double Alpha) {
302 return 2*pc*L*
Y(Alpha);
313 double f_interp(
double E,
double f1,
double E1,
double f2,
double E2) {
315 if (E < E1 || E > E2) {
318 if (f1 <= 0 || f2 <= 0)
return 0;
319 b = log(f2/f1)/(E1 - E2);
321 f_epc = a *
exp(-b*E);
329 double mu_calc(
double L,
double pc,
double Alpha) {
330 return pow(pc,2)/
VC::mc2 * pow(sin(Alpha),2) / (2 *
B(L));
338 return pow(10, (-0.3145*L + 3.9043));
355 double max(
double v1,
double v2) {
356 return (v1>v2)?v1:v2;
361 std::ostringstream o;
368 if ((iter % d_iter) == 0) {
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
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.
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.