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) {