16 #include "../VariousFunctions/bisection.h"
18 #include "../Logging/Output.h"
19 #include "../Exceptions/error.h"
25 #define double_zero 1.e-21
26 #define min_Dxx 1.e-21
28 #if defined(_WIN32) || defined(_WIN64)
29 #define strncasecmp _strnicmp
30 #define strcasecmp _stricmp
43 bool Load_Failed =
false;
45 if (!arr.initialized) {
53 this->DxxParameters = DxxParameters;
55 this->type = DxxParameters.
DxxType;
56 this->useScale = DxxParameters.
useScale;
59 throw error_msg(
"DXX_GET",
"Undefined load or calculate choice.");
73 err.
add(
"UNKNOWN_ERROR",
"Unknown error while loading diffusion coefficient file");
81 Output::echo(3,
"Calculating %s...\n", arr.name.c_str());
82 Calculate(grid.
L, grid.
epc, grid.
alpha, DxxParameters);
84 arr = arr * (60.0 * 60 * 24);
106 if (this->type ==
"DCT_Dpp" || this->type ==
"DCT_DppLpp") arr.times_equal(grid.
pc.
arr.
times(grid.
pc.
arr));
107 if (this->type ==
"DCT_Dpa" || this->type ==
"DCT_DpaLpp") arr.times_equal(grid.
pc.
arr);
109 arr.change_ind = clock();
121 bool positivonly=
true;
125 if (DxxParameters.
DxxType ==
"DCT_Daa" || DxxParameters.
DxxType ==
"DCT_DaaLpp") {
127 }
else if (DxxParameters.
DxxType ==
"DCT_Dpp" || DxxParameters.
DxxType ==
"DCT_DppLpp") {
129 }
else if (DxxParameters.
DxxType ==
"DCT_Dpa" || DxxParameters.
DxxType ==
"DCT_DpaLpp") {
133 throw error_msg(
"DXX_TYPE",
"Unknown diffusion coefficient");
138 for (il = 0; il < L.
size ; il++) {
139 for (im = 0; im < epc.
size; im++) {
140 for (ia = 0; ia < Alpha.
size; ia++) {
143 arr[il][im][ia] =
Dxx_ba(L.
arr[il][im][ia], epc.
arr[il][im][ia]/
VC::mc2, Alpha.
arr[il][im][ia], int_Dxx_loc, DxxParameters) / 2;
145 if (positivonly && (arr[il][im][ia] < min_Dxx)) {
146 arr[il][im][ia] = min_Dxx;
147 }
else if (!(arr[il][im][ia] >= 0) && !(arr[il][im][ia] <= 0)) {
149 arr[il][im][ia] = min_Dxx;
166 if (!arr.initialized) {
167 throw error_msg(
"DXX_UNINITIALIZED_MATRIX",
"Matrix has not been initialized.");
172 if (filetype ==
"IFT_GRID") {
173 LoadDiffusionCoefficientFromFileWithGrid(L, pc, alpha, D_filename);
174 }
else if (filetype ==
"IFT_GRID_LPA") {
175 LoadDiffusionCoefficientFromFileWithGrid(L, pc, alpha, D_filename,
"IFT_GRID_LPA");
176 }
else if (filetype ==
"IFT_GRID_APL") {
177 LoadDiffusionCoefficientFromFileWithGrid(L, pc, alpha, D_filename,
"IFT_GRID_APL");
178 }
else if (filetype ==
"IFT_PLANE") {
180 LoadDiffusionCoefficientFromPlaneFile(L, pc, alpha, D_filename);
182 throw error_msg(
"DXX_LOAD_TYPE_ERR",
"Unknown coefficients file format %s", filetype.c_str());
195 ifstream input_D(D_filename.c_str());
197 if (input_D != NULL && !input_D.eof()) {
198 for (il = 0; il < L.
size; il++)
199 for (im = 0; im < pc.
size; im++)
200 for (ia = 0; ia < alpha.
size; ia++) {
202 input_D >> arr[il][im][ia];
203 if (input_D == NULL) {
204 throw error_msg(
"LOAD_ERROR",
"Wrong coefficient file %s (to few values)\n", D_filename.c_str());
208 throw error_msg(
"DXX_LOAD_GRID_ERR",
"error reading input file %s\n", D_filename.c_str());
222 double Load_L, Load_epc, Load_alpha;
227 ifstream input_D(D_filename.c_str());
229 if (input_D == NULL) {
230 throw error_msg(
"DXX_FILE_OPEN_ERR",
"Diffusion coefficient file %s not found.", D_filename.c_str());
235 while (strcasecmp(tmp.c_str(),
"ZONE") != 0 && !input_D.eof() ) {
239 throw error_msg(
"DXX_LOAD_GRID_ERR",
"Keyword 'ZONE' not found in the diffusion coefficient file: %s\n", D_filename.c_str());
243 input_D.ignore(9999,
'\n');
245 if (!input_D.eof()) {
247 if (gridOrder ==
"IFT_GRID_LPA") {
248 for (il = 0; il < L.
size; il++)
249 for (im = 0; im < epc.
size; im++)
250 for (ia = 0; ia < alpha.
size; ia++) {
251 input_D >> Load_L >> Load_epc >> Load_alpha;
254 input_D >> arr[il][im][ia];
256 if (fabs(Load_L - L.
arr[il][im][ia]) > err || fabs(Load_epc - epc.
arr[il][im][ia]) > err || fabs(Load_alpha - alpha.
arr[il][im][ia]) > err) {
257 throw error_msg(
"DXX_LOAD_GRID_ERR",
"Loading %s: grid mismatch.\nLoaded: L=%e, epc=%e, alpha=%e\nGrid: L=%e, epc=%e, alpha=%e\n", D_filename.c_str(), Load_L, Load_epc, Load_alpha, L.
arr[il][im][ia], epc.
arr[il][im][ia], alpha.
arr[il][im][ia]);
260 }
else if (gridOrder ==
"IFT_GRID_APL") {
261 for (ia = 0; ia < alpha.
size; ia++)
262 for (im = 0; im < epc.
size; im++)
263 for (il = 0; il < L.
size; il++) {
264 input_D >> Load_alpha >> Load_epc >> Load_L ;
267 input_D >> arr[il][im][ia];
269 if (fabs(Load_L - L.
arr[il][im][ia]) > err || fabs(Load_epc - epc.
arr[il][im][ia]) > err || fabs(Load_alpha - alpha.
arr[il][im][ia]) > err) {
270 throw error_msg(
"DXX_LOAD_GRID_ERR",
"Loading %s: grid mismatch.\nLoaded: L=%e, epc=%e, alpha=%e\nGrid: L=%e, epc=%e, alpha=%e\n", D_filename.c_str(), Load_L, Load_epc, Load_alpha, L.
arr[il][im][ia], epc.
arr[il][im][ia], alpha.
arr[il][im][ia]);
274 throw error_msg(
"UNKNOWN_GRID_TYPE",
"Unknown grid order: %s", gridOrder.c_str());
286 if (DLLType ==
"DLLT_BE") {
287 MakeDLL_BE(L, pc, alpha, Kp);
289 else if (DLLType ==
"DLLT_B") {
290 MakeDLL_B(L, pc, alpha, Kp);
292 else if (DLLType ==
"DLLT_FAKE") {
293 MakeDLL_FAKE(L, pc, alpha, Kp);
295 else if (DLLType ==
"DLLT_Ozeke") {
296 MakeDLL_Ozeke(L, pc, alpha, Kp);
298 else if (DLLType ==
"DLLT_Ozeke_ME") {
299 MakeDLL_Ozeke_ME(L, pc, alpha, Kp);
302 throw error_msg(
"DLL_ERROR",
"Unknown DLL type");
304 arr.change_ind = clock();
313 for (il = 0; il < L.
size; il++) {
314 for (im = 0; im < pc.
size; im++) {
315 for (ia = 0; ia < alpha.
size; ia++) {
316 arr[il][im][ia] =
Df(L.
arr[il][im][ia], Kp);
328 for (il = 0; il < L.
size; il++) {
329 for (im = 0; im < pc.
size; im++) {
330 for (ia = 0; ia < alpha.
size; ia++) {
331 arr[il][im][ia] =
Df(L.
arr[il][im][ia], Kp) + 3e3*pow(10.0, 0.506*Kp-9.325)*pow(L.
arr[il][im][ia],5);;
343 for (il = 0; il < L.
size; il++) {
344 for (im = 0; im < pc.
size; im++) {
345 for (ia = 0; ia < alpha.
size; ia++) {
352 arr[il][im][ia] = arr[il][im][ia] +
Df(L.
arr[il][im][ia], Kp);
366 for (il = 0; il < L.
size; il++) {
367 for (im = 0; im < pc.
size; im++) {
368 for (ia = 0; ia < alpha.
size; ia++) {
369 arr[il][im][ia] =
VF::max(
Df(L.
arr[il][im][ia], Kp), 1.0/100);
386 for (il = 0; il < L.
size; il++) {
387 for (im = 0; im < pc.
size; im++) {
388 for (ia = 0; ia < alpha.
size; ia++) {
405 for (il = 0; il < L.
size; il++) {
406 for (im = 0; im < pc.
size; im++) {
407 for (ia = 0; ia < alpha.
size; ia++) {
421 double scaleCoeff = 0;
425 bool need_changes =
false;
426 for (Dxx_it = 0; Dxx_it < DxxList.size(); Dxx_it++) {
427 if (time >= DxxList[Dxx_it].DxxParameters.time_start && time < DxxList[Dxx_it].DxxParameters.time_end) {
429 if (!DxxList[Dxx_it].is_active) {
434 if (DxxList[Dxx_it].useScale) {
440 if (DxxList[Dxx_it].is_active) {
446 if (need_changes ==
false)
return false;
450 CurrentDxx.arr.change_ind = clock();
452 for (Dxx_it = 0; Dxx_it < DxxList.size(); Dxx_it++) {
454 if (time >= DxxList[Dxx_it].DxxParameters.time_start && time < DxxList[Dxx_it].DxxParameters.time_end) {
456 if (!DxxList[Dxx_it].is_active) {
457 DxxList[Dxx_it].is_active =
true;
458 Output::echo(2,
"Activated %s.\n", DxxList[Dxx_it].arr.name.c_str());
461 if (DxxList[Dxx_it].useScale) {
463 scaleCoeff = DxxList[Dxx_it].Scale(Kp);
464 Output::echo(2,
"Bw scaled %s by %f.\n", DxxList[Dxx_it].arr.name.c_str(), scaleCoeff);
466 CurrentDxx.arr = CurrentDxx.arr + DxxList[Dxx_it].arr * scaleCoeff;
468 CurrentDxx.arr = CurrentDxx.arr + DxxList[Dxx_it].arr;
472 if (DxxList[Dxx_it].is_active) {
473 DxxList[Dxx_it].is_active =
false;
474 Output::echo(2,
"Deactivated %s.", DxxList[Dxx_it].arr.name.c_str());
485 double Bw2=0, Bw2_0=0;
486 if (!this->useScale)
return 1;
487 if (DxxParameters.waveType ==
"WT_CHORUS_DAY" || DxxParameters.waveType ==
"WT_CHORUS_NIGHT") {
489 if (DxxParameters.DxxKp <= 2.333) Bw2_0 = 2.0*pow(10, 0.73 + 0.91 * DxxParameters.DxxKp);
490 if (DxxParameters.DxxKp > 2.333) Bw2_0 = 2.0*pow(10, 2.5 + 0.18 * DxxParameters.DxxKp);
491 if (Kp <= 2.333) Bw2 = 2.0*pow(10, 0.73 + 0.91 * Kp);
492 if (Kp > 2.333) Bw2 = 2.0*pow(10, 2.5 + 0.18 * Kp);
502 }
else if (DxxParameters.waveType ==
"WT_HISS_PP" || DxxParameters.waveType ==
"WT_HISS_PLUME") {
505 return pow(Kp/DxxParameters.DxxKp, 2);
507 else if (DxxParameters.waveType ==
"WT_HISS_REAL") {
510 return pow(10, 0.2607*Kp - 0.01547*Kp*Kp - 0.2607*DxxParameters.DxxKp + 0.01547*DxxParameters.DxxKp*DxxParameters.DxxKp);
511 }
else if (DxxParameters.waveType ==
"WT_EMIC_PLUME" ) {
513 return pow(Kp/DxxParameters.DxxKp, 4);
515 }
else if (DxxParameters.waveType ==
"WT_LIGHTNING") {
516 return pow(Kp/DxxParameters.DxxKp, 2);
538 double Alpha_ne(
double pangle,
double lambda,
double L) {
539 return asin( sqrt(
f1(lambda)) * sin(pangle) );
546 double f1(
double lambda) {
547 return sqrt(4.0 - 3.0 * pow(cos(lambda), 2))/(pow(cos(lambda), 6));
554 double B(
double lambda,
double L) {
562 return pow(x, 6) + (3.0 * pow(sin(Alpha), 4)) * x - 4.0 * pow(sin(Alpha), 4);
572 double c1 = 2.0*s*(-1+epsilon);
573 double c2 = 1.0-4.0*epsilon+pow(epsilon,2);
574 double c3 = -s*(-1+epsilon)*(b+4.0*epsilon)/2.0;
575 double c4 = epsilon*(b+epsilon) ;
576 double g = pow(x,4) + c1*pow(x,3) + c2*pow(x,2) + c3*x + c4;
577 return y*(pow((x-s),2))*(pow((x+s*epsilon),2))/(x*g);
591 double up = 2.0*(x+a)/(beta*mu);
592 double down = 2.0*x-1.0/Alpha_star*(1.0/(x-s)+DxxParameters.
eta1*epsilon/(x+s*epsilon)+DxxParameters.
eta2*epsilon/(4.0*x+s*epsilon)+DxxParameters.
eta3*epsilon/(16.0*x+s*epsilon))+
593 x*1.0/Alpha_star*(1.0/(pow((x-s), 2))+DxxParameters.
eta1*epsilon/(pow((x+s*epsilon), 2))+4.0*epsilon*DxxParameters.
eta2/(pow((4.0*x+s*epsilon), 2))+16.0*epsilon*DxxParameters.
eta3/(pow((16.0*x+s*epsilon), 2)));
606 double h = (b - a)/M;
607 double s = func(a, DxxParameters);
610 for (k = 1; k <= (M-1); k++) {
612 s = s + func(x, DxxParameters);
614 s = h * (func(a, DxxParameters) + func(b, DxxParameters))/2.0 + h*s;
627 DxxParameters.
L = L1;
628 DxxParameters.
EMeV = EMeV;
629 DxxParameters.
Alpha = Alpha;
636 return (1.0/(1.38-0.32*(sin(DxxParameters.
Alpha) + pow(sin(DxxParameters.
Alpha),0.5))) *
quad1(int_Dxx_loc, DxxParameters.
lam_min*
VC::pi/180., lambda_m, DxxParameters.
nint, DxxParameters));
649 if (DxxParameters.
Alpha*180/
VC::pi < 0.001) Daa = 0;
650 return Daa * cos(
Alpha_ne(DxxParameters.
Alpha, lambda, DxxParameters.
L)) * pow(cos(lambda), 7) /(pow(cos(DxxParameters.
Alpha), 2));
663 if (DxxParameters.
Alpha*180/
VC::pi < 0.001) Dpp = 0;
664 return Dpp * (cos(lambda) * pow((1 + 3*pow(sin(lambda), 2)), 0.5)) / cos(
Alpha_ne(DxxParameters.
Alpha, lambda, DxxParameters.
L));
677 if (DxxParameters.
Alpha*180/
VC::pi < 0.001) Dpa = 0;
678 return Dpa * sin(
Alpha_ne(DxxParameters.
Alpha, lambda, DxxParameters.
L)) * pow(cos(lambda), 7) / (cos(DxxParameters.
Alpha) * sin(DxxParameters.
Alpha));
690 double Dxx_local(
double lambda,
double Dxx_root(
double Omega_e,
double x,
double mu,
double su,
double y,
double beta,
double a,
double b,
double Alpha_star,
double s,
double epsilon,
double d_x,
double x_m,
double R,
DxxParameters_structure DxxParameters),
DxxParameters_structure DxxParameters ) {
694 double N_0 = 124.0 * pow((3.0 / DxxParameters.
L), 4);
696 double Dxx_f = DxxParameters.
f;
697 DxxParameters.
f = Omega_p/Omega_e_eq;
705 Omega_p = Dxx_f * Omega_e_eq;
715 N_0 += 36.0*pow((3.0/DxxParameters.
L),3.5) * cos((LT-(7.7*pow((3.0/DxxParameters.
L),2.0)+12))*
VC::pi/12);
718 N_0 = pow(10, (-0.3145*DxxParameters.
L + 3.9043));
719 if (DxxParameters.
numberDensity ==
"ND_PL_CA_STARKS" && DxxParameters.
L <= 2)
721 double N1, N2, L1, L2, A,
B;
722 N1=10-6.;N2=8.3-6.;L1=1.2;L2=5.9;
723 A=(N2-N1)/(L2-L1); B=N2-A*L2;
724 N_0 = pow(10, A*DxxParameters.
L+B);
728 N_0 = 1e3*pow(4.0/DxxParameters.
L, 4);
731 throw error_msg(
"DXX_DENSITY_ERROR",
"Unknown number density type: %s\n", DxxParameters.
numberDensity.c_str());
735 double epsilon = 1.0/1837;
736 double beta = (pow((DxxParameters.
EMeV * (DxxParameters.
EMeV + 2)), 0.5)) / (DxxParameters.
EMeV + 1);
737 double mu = cos(
Alpha_ne(DxxParameters.
Alpha, lambda, DxxParameters.
L));
740 double su = sin(
Alpha_ne(DxxParameters.
Alpha, lambda, DxxParameters.
L));
743 double gamma = DxxParameters.
EMeV + 1;
744 double a = DxxParameters.
s * lam_sp / gamma;
746 double Alpha_star = pow((Omega_e / Omega_p), 2);
747 double b = (1 + epsilon)/Alpha_star;
758 Omega_m = DxxParameters.
Omega_m * Omega_e_eq;
759 d_omega = DxxParameters.
d_omega * Omega_e_eq ;
760 Omega_1 = DxxParameters.
omega_lc * Omega_e_eq;
761 Omega_2 = DxxParameters.
omega_uc * Omega_e_eq;
762 }
else if (DxxParameters.
Omega_mType ==
"OT_ABSOLUTE") {
768 throw error_msg(
"DXX_OMEGA_M_TYPE",
"Unknown Omega_m type: ", DxxParameters.
Omega_mType.c_str());
779 Bw = pow(10, 0.75 + 0.04*(lambda*180/
VC::pi))*1.e-8 * DxxParameters.
Bw;
781 Bw = DxxParameters.
Bw;
783 double R = pow((Bw /
B(lambda, DxxParameters.
L)), 2);
790 if (DxxParameters.
particle ==
"PT_ELECTRON") {
791 x_m = Omega_m / Omega_e;
792 d_x = d_omega / Omega_e;
793 x_1 = Omega_1 / Omega_e;
794 x_2 = Omega_2 / Omega_e;
795 }
else if (DxxParameters.
particle ==
"PT_PROTON") {
797 x_m = Omega_m / Omega_e * epsilon / cfm;
798 d_x = d_omega / Omega_e * epsilon / cfm;
799 x_1 = Omega_1 / Omega_e * epsilon / cfm;
800 x_2 = Omega_2 / Omega_e * epsilon / cfm;
803 std::vector<double> xx =
rrouts(x_1, x_2, DxxParameters.
eta1, DxxParameters.
eta2, DxxParameters.
eta3, epsilon, beta, mu, Alpha_star, a, DxxParameters);
806 double Dxx = 0, Dxx_tmp = 0;
807 for (
unsigned int i = 1; i <= xx.size();i++) {
809 y = (x + a) / (beta * mu);
810 Dxx_tmp = Dxx_root(Omega_e, x, mu, su, y, beta, a, b, Alpha_star, DxxParameters.
s, epsilon, d_x, x_m, R, DxxParameters);
813 if (xx.size() == 0) {
830 double Daa_root(
double Omega_e,
double x,
double mu,
double su,
double y,
double beta,
double a,
double b,
double Alpha_star,
double s,
double epsilon,
double d_x,
double x_m,
double R,
DxxParameters_structure DxxParameters) {
831 double fcap =
F_cap2(x,y,a,beta,mu,s,epsilon,Alpha_star, DxxParameters);
832 return (
VC::pi/2) * (1.0/DxxParameters.
nu) * Omega_e * (1.0/(pow((DxxParameters.
EMeV + 1), 2))) *
833 ((R * pow((1 - x * mu / (y * beta)), 2) * fabs(fcap))/
834 (d_x * fabs(beta * mu - fcap)))*
835 exp( - (pow(((x - x_m) / d_x), 2)));
848 double Dpa_root(
double Omega_e,
double x,
double mu,
double su,
double y,
double beta,
double a,
double b,
double Alpha_star,
double s,
double epsilon,
double d_x,
double x_m,
double R,
DxxParameters_structure DxxParameters) {
849 double fcap =
F_cap2(x,y,a,beta,mu,s,epsilon,Alpha_star,DxxParameters);
850 return -(
VC::pi/2)*(1.0/DxxParameters.
nu)* Omega_e*su/beta*(1.0/(pow((DxxParameters.
EMeV+1), 2)))*
851 ((R*(1-x/y*mu/beta)*(x/y)* fabs(fcap))/
852 (d_x * fabs(beta * mu - fcap)))*
853 exp( -(pow(((x - x_m) / d_x), 2)));
866 double Dpp_root(
double Omega_e,
double x,
double mu,
double su,
double y,
double beta,
double a,
double b,
double Alpha_star,
double s,
double epsilon,
double d_x,
double x_m,
double R,
DxxParameters_structure DxxParameters) {
867 double fcap =
F_cap2(x,y,a,beta,mu,s,epsilon,Alpha_star,DxxParameters);
868 return (
VC::pi/2) * (1.0/DxxParameters.
nu) * (1./pow(beta, 2)) * (1-pow(mu, 2)) * Omega_e * (1.0/(pow((DxxParameters.
EMeV + 1), 2))) *
869 ((R * (pow(x/y, 2)) * fabs(fcap))/
870 (d_x*fabs(beta*mu - fcap)))*
871 exp(-(pow(((x-x_m)/d_x),2)));
879 std::vector<double>
rrouts(
double x_1,
double x_2,
double eta1,
double eta2,
double eta3,
double epsilon,
double beta,
double mu,
double Alpha_star,
double a,
DxxParameters_structure DxxParameters) {
887 double c0 = -DxxParameters.
s;
888 double c1 = DxxParameters.
s * epsilon;
889 double c2 = 1.0/4. * DxxParameters.
s * epsilon;
890 double c3 = 1.0/16. * DxxParameters.
s * epsilon;
893 double d1 = eta1*epsilon;
894 double d2 = eta2/4.*epsilon;
895 double d3 = eta3/16.*epsilon;
897 double g3 = d0 + d1 + d2 + d3;
898 double g2 = d0*(c1+c2+c3) + d1*(c0+c2+c3) + d2*(c0+c1+c3) + d3*(c0+c1+c2);
899 double g1 = d0*(c1*c2 + c1*c3 + c2*c3) + d1*(c0*c2 + c0*c3 + c2*c3) + d2*(c0*c1 + c0*c3 + c1*c3) + d3*(c0*c1 + c0*c2 + c1*c2);
900 double g0 = d0*c1*c2*c3 + d1*c0*c2*c3 + d2*c0*c1*c3 + d3*c0*c1*c2;
903 double h3 = c0 + c1 + c2 + c3;
904 double h2 = c0*c1 + c0*c2 + c0*c3 + c1*c2 + c1*c3 + c2*c3;
905 double h1 = c0*c1*c2 + c0*c1*c3 + c0*c2*c3 + c1*c2*c3;
906 double h0 = c0*c1*c2*c3;
910 A[0] = (pow(beta,2)*pow(mu,2)-1.)*h4;
911 A[1] = (pow(beta,2)*pow(mu,2)-1.)*h3-2*a*h4;
912 A[2] = (pow(beta,2)*pow(mu,2)-1.)*h2-2*a*h3-pow(beta,2)*pow(mu,2)*1./Alpha_star*g3-pow(a,2)*h4;
913 A[3] = (pow(beta,2)*pow(mu,2)-1.)*h1-2*a*h2-pow(beta,2)*pow(mu,2)*1./Alpha_star*g2-pow(a,2)*h3;
914 A[4] = (pow(beta,2)*pow(mu,2)-1.)*h0-2*a*h1-pow(beta,2)*pow(mu,2)*1./Alpha_star*g1-pow(a,2)*h2;
915 A[5] = -2*a*h0-pow(beta,2)*pow(mu,2)*1./Alpha_star*g0-pow(a,2)*h1;
926 std::vector<double> xx;
928 double r[7], wr[6],wi[6],quad[2];
932 quad[0] = 2.71828e-1;
933 quad[1] = 3.14159e-1;
936 numr =
roots(r, n, wr, wi);
938 for (
int i = 0; i < numr; i++) {
939 if (fabs(wi[i]) < double_zero && wr[i] < x_2 && wr[i] > x_1) {
952 double Df(
double L,
double Kp) {
953 return pow(10, 0.506*Kp-9.325)*pow(L,10);
965 return 6.62e-13 * pow(L,8) * pow(10, (pow(L,2) * -0.0327 + L * 0.625 - pow(Kp,2) * 0.0108 + Kp * 0.499) );
972 return 2.16e-8 * pow(L,6) * pow(10,( 0.217 * L + 0.461 * Kp ));
1031 unsigned int Dxx_it;
1038 if (Dxx_tmp.
type ==
"DCT_Daa")
1039 Daa.
DxxList.push_back(Dxx_tmp);
1040 else if (Dxx_tmp.
type ==
"DCT_DaaLpp")
1041 DaaLpp.
DxxList.push_back(Dxx_tmp);
1045 if (Dxx_tmp.
type ==
"DCT_Dpp" || Dxx_tmp.
type ==
"DCT_Dpcpc")
1046 Dpcpc.
DxxList.push_back(Dxx_tmp);
1047 else if (Dxx_tmp.
type ==
"DCT_DppLpp" || Dxx_tmp.
type ==
"DCT_DpcpcLpp")
1048 DpcpcLpp.
DxxList.push_back(Dxx_tmp);
1052 if (Dxx_tmp.
type ==
"DCT_Dpca" || Dxx_tmp.
type ==
"DCT_Dpa")
1053 Dpca.
DxxList.push_back(Dxx_tmp);
1054 else if (Dxx_tmp.
type ==
"DCT_DpcaLpp" || Dxx_tmp.
type ==
"DCT_DpaLpp")
1055 DpcaLpp.
DxxList.push_back(Dxx_tmp);
1061 int Daa_it, Dpcpc_it, Dpca_it;
1063 for (Dpca_it = 0; Dpca_it < Dpca.
DxxList.size(); Dpca_it++) {
1066 for (Daa_it = 0; Daa_it < Daa.
DxxList.size(); Daa_it++) {
1067 if (Daa.
DxxList[Daa_it].DxxParameters.waveType == Dpca.
DxxList[Dpca_it].DxxParameters.waveType)
1071 for (Dpcpc_it = 0; Dpcpc_it < Dpcpc.
DxxList.size(); Dpcpc_it++) {
1072 if (Dpcpc.
DxxList[Dpcpc_it].DxxParameters.waveType == Dpca.
DxxList[Dpca_it].DxxParameters.waveType)
1077 for (iL = 0; iL < localDiffusionsGrid.
L.
size; iL++)
1078 for (im = 0; im < localDiffusionsGrid.
pc.
size; im++)
1079 for (ia = 0; ia < localDiffusionsGrid.
alpha.
size; ia++)
1080 if (Daa.
DxxList[Daa_it].arr[iL][im][ia] * Dpcpc.
DxxList[Dpcpc_it].arr[iL][im][ia] * 0.99 <= Dpca.
DxxList[Dpca_it].arr[iL][im][ia] * Dpca.
DxxList[Dpca_it].arr[iL][im][ia])
1081 Output::echo(0,
"WARNING: Diffusion coefficient error. %s: Dpp * Daa * 0.99 <= Dpa^2. \nAt (%d,%d,%d): L=%f, epc=%f, alpha=%f; Dpp=%e, Daa=%e, Dpa=%e\n",
1082 Dpca.
DxxList[Dpca_it].DxxParameters.waveType.c_str(),
1084 localDiffusionsGrid.
L.
arr[iL][im][ia], localDiffusionsGrid.
epc.
arr[iL][im][ia], localDiffusionsGrid.
alpha.
arr[iL][im][ia],
1085 Dpcpc.
DxxList[Dpcpc_it].arr[iL][im][ia], Daa.
DxxList[Daa_it].arr[iL][im][ia], Dpca.
DxxList[Dpca_it].arr[iL][im][ia]);
1089 for (Dpca_it = 0; Dpca_it < DpcaLpp.
DxxList.size(); Dpca_it++) {
1092 for (Daa_it = 0; Daa_it < DaaLpp.
DxxList.size(); Daa_it++) {
1093 if (DaaLpp.
DxxList[Daa_it].DxxParameters.waveType == DpcaLpp.
DxxList[Dpca_it].DxxParameters.waveType)
1097 for (Dpcpc_it = 0; Dpcpc_it < DpcpcLpp.
DxxList.size(); Dpcpc_it++) {
1098 if (DpcpcLpp.
DxxList[Dpcpc_it].DxxParameters.waveType == DpcaLpp.
DxxList[Dpca_it].DxxParameters.waveType)
1103 for (iL = 0; iL < localDiffusionsGrid.
L.
size; iL++)
1104 for (im = 0; im < localDiffusionsGrid.
pc.
size; im++)
1105 for (ia = 0; ia < localDiffusionsGrid.
alpha.
size; ia++)
1106 if (DaaLpp.
DxxList[Daa_it].arr[iL][im][ia] * DpcpcLpp.
DxxList[Dpcpc_it].arr[iL][im][ia] * 0.99 <= DpcaLpp.
DxxList[Dpca_it].arr[iL][im][ia] * DpcaLpp.
DxxList[Dpca_it].arr[iL][im][ia])
1107 Output::echo(0,
"WARNING: Diffusion coefficient error. %s: Dpp * Daa * 0.99 <= Dpa^2. \nAt (%d,%d,%d): L=%f, epc=%f, alpha=%f; Dpp=%e, Daa=%e, Dpa=%e\n",
1108 DpcaLpp.
DxxList[Dpca_it].DxxParameters.waveType.c_str(),
1110 localDiffusionsGrid.
L.
arr[iL][im][ia], localDiffusionsGrid.
epc.
arr[iL][im][ia], localDiffusionsGrid.
alpha.
arr[iL][im][ia],
1111 DpcpcLpp.
DxxList[Dpcpc_it].arr[iL][im][ia], DaaLpp.
DxxList[Daa_it].arr[iL][im][ia], DpcaLpp.
DxxList[Dpca_it].arr[iL][im][ia]);
1118 DLL.
MakeDLL(radialDiffusionGrid.
L, radialDiffusionGrid.
pc, radialDiffusionGrid.
alpha, parameters.
Kp[0], parameters.
DLLType);
1131 output1D <<
"Variables = \"time\", \"Kp\", \"Boundary_fluxes\", \"Lpp\", \"tau\" ";
1134 unsigned int Dxx_it;
1135 for (Dxx_it = 0; Dxx_it < Daa.
DxxList.size(); Dxx_it++) {
1136 output1D <<
"\"" << Daa.
DxxList[Dxx_it].arr.name <<
"\"" <<
", ";
1138 for (Dxx_it = 0; Dxx_it < Dpcpc.
DxxList.size(); Dxx_it++) {
1139 output1D <<
"\"" << Dpcpc.
DxxList[Dxx_it].arr.name <<
"\"" <<
", ";
1141 for (Dxx_it = 0; Dxx_it < Dpca.
DxxList.size(); Dxx_it++) {
1142 output1D <<
"\"" << Dpca.
DxxList[Dxx_it].arr.name <<
"\"" <<
", ";
1145 for (Dxx_it = 0; Dxx_it < DaaLpp.
DxxList.size(); Dxx_it++) {
1146 output1D <<
"\"" << DaaLpp.
DxxList[Dxx_it].arr.name <<
"\"" <<
", ";
1148 for (Dxx_it = 0; Dxx_it < DpcpcLpp.
DxxList.size(); Dxx_it++) {
1149 output1D <<
"\"" << DpcpcLpp.
DxxList[Dxx_it].arr.name <<
"\"" <<
", ";
1151 for (Dxx_it = 0; Dxx_it < DpcaLpp.
DxxList.size(); Dxx_it++) {
1152 output1D <<
"\"" << DpcaLpp.
DxxList[Dxx_it].arr.name <<
"\"" <<
", ";
1156 output1D <<
"ZONE T=\"1d-output\" " << endl;
1165 output1D << time <<
"\t" << parameters.
Kp[iteration] <<
"\t" << parameters.
Bf[iteration] <<
"\t" << parameters.
Lpp[iteration] <<
"\t" << parameters.
tau[iteration] <<
"\t";
1167 unsigned int Dxx_it;
1168 for (Dxx_it = 0; Dxx_it < Daa.
DxxList.size(); Dxx_it++) {
1169 output1D << Daa.
DxxList[Dxx_it].Scale(parameters.
Kp[iteration]) <<
"\t";
1171 for (Dxx_it = 0; Dxx_it < Dpcpc.
DxxList.size(); Dxx_it++) {
1172 output1D << Dpcpc.
DxxList[Dxx_it].Scale(parameters.
Kp[iteration]) <<
"\t";
1174 for (Dxx_it = 0; Dxx_it < Dpca.
DxxList.size(); Dxx_it++) {
1175 output1D << Dpca.
DxxList[Dxx_it].Scale(parameters.
Kp[iteration]) <<
"\t";
1178 for (Dxx_it = 0; Dxx_it < DaaLpp.
DxxList.size(); Dxx_it++) {
1179 output1D << DaaLpp.
DxxList[Dxx_it].Scale(parameters.
Kp[iteration]) <<
"\t";
1181 for (Dxx_it = 0; Dxx_it < DpcpcLpp.
DxxList.size(); Dxx_it++) {
1182 output1D << DpcpcLpp.
DxxList[Dxx_it].Scale(parameters.
Kp[iteration]) <<
"\t";
1184 for (Dxx_it = 0; Dxx_it < DpcaLpp.
DxxList.size(); Dxx_it++) {
1185 output1D << DpcaLpp.
DxxList[Dxx_it].Scale(parameters.
Kp[iteration]) <<
"\t";
Holds list of instances of DiffusionCoefficient class of same type (like Daa, Dpp, etc), but produced by different waves (Daa_chorus, Daa_EMIC, etc).
GridElement alpha
grid elements to be used
Matrix3D< double > arr
array of grid points
Matrix1D< double > Kp
Kp array.
double eta2
eta values, calculation parameters
void add(single_error err)
double int_Dpa_loc(double lambda, DxxParameters_structure DxxParameters)
static const double m
mass of electron in grams
bool ActivateAndScale(double time, double Kp)
double Daa_root(double Omega_e, double x, double mu, double su, double y, double beta, double a, double b, double Alpha_star, double s, double epsilon, double d_x, double x_m, double R, DxxParameters_structure DxxParameters)
Array of values of coordinate axes.
string filetype
Type of the file (w/wo grid etc)
bool useScale
Flag, use scaling if equal to true.
string Omega_mType
Omega_m type. Check StrToVal(string input, Omega_mTypes &place) for known values. ...
double Dpp_root(double Omega_e, double x, double mu, double su, double y, double beta, double a, double b, double Alpha_star, double s, double epsilon, double d_x, double x_m, double R, DxxParameters_structure DxxParameters)
double max(double v1, double v2)
Return maximum.
void Create_Dxx(DxxParameters_structure DxxParameters, Grid &grid)
Create diffusion coefficients.
int roots(double *a, int n, double *wr, double *wi)
Extract individual real or complex roots from list of quadratic factors.
double F_cap(double x, double y, double b, double s, double epsilon, DxxParameters_structure DxxParameters)
void Output1DValues(ofstream &output1D, DiffusionCoefficientsGroup &Daa, DiffusionCoefficientsGroup &Dpcpc, DiffusionCoefficientsGroup &Dpca, DiffusionCoefficientsGroup &DaaLpp, DiffusionCoefficientsGroup &DpcpcLpp, DiffusionCoefficientsGroup &DpcaLpp, double &time, Parameters_structure ¶meters, int iteration)
Function used in main() to print values of diffusion coefficients in output file. ...
double Df(double L, double Kp)
bool LoadDiffusionCoefficientFromPlaneFile(GridElement &L, GridElement &pc, GridElement &alpha, string D_filename)
Load Dxx from file without grid.
Holds diffusion coefficient matrix and routines to load and calculate it.
bool useAlphaDiffusion
Using diffusions flags.
double eta1
eta values, calculation parameters
double eta3
eta values, calculation parameters
string particle
Type of particles, produced wave? Ions or electrons. Check StrToVal(string input, ParticleTypes &plac...
bool LoadDiffusionCoefficient(GridElement &L, GridElement &pc, GridElement &alpha, string D_filename, string filetype="IFT_GRID")
double lam_min
Minimum latitude.
double Alpha_ne(double pangle, double lambda, double L)
DiffusionCoefficient CurrentDxx
flag, indicated that the initialization was passed
static const double c_cgs
speed of light, cm/s
void MakeDLL100(GridElement &L, GridElement &pc, GridElement &alpha, double Kp)
Making DLL.
Matrix1D< double > Bf
Boundary flux array.
static const double exp
exponential
void MakeDLL_FAKE(GridElement &L, GridElement &pc, GridElement &alpha, double Kp)
Making DLL.
void echo(int msglvl, const char *format,...)
Basic function! Write the message to the file.
Struct that holds various parameters to be used for Dxx.
double func_tmp(double x, double Alpha)
void MakeDLL_Ozeke_ME(GridElement &L, GridElement &pc, GridElement &alpha, double Kp)
Making DLL_ME.
double d_omega
d_omega value.
void AllocateMemory(int size_x, int size_y, int size_z)
string loadOrCalculate
Load diffusion coefficient or calculate flag.
double quad1(double(*func)(double lambda, DxxParameters_structure DxxParameters), double a, double b, int M, DxxParameters_structure DxxParameters)
double MLT_averaging
MLT averaging part of the orbit.
static const double mc2
m*c^2 = 0.511 MeV
double Dxx_local(double lambda, double Dxx_root(double Omega_e, double x, double mu, double su, double y, double beta, double a, double b, double Alpha_star, double s, double epsilon, double d_x, double x_m, double R, DxxParameters_structure DxxParameters), DxxParameters_structure DxxParameters)
void Calculate(GridElement &L, GridElement &epc, GridElement &alpha, DxxParameters_structure DxxParameters)
function calculate Dxx
int nint
Number of points in integral.
string filename
File name for loading or saving diffusion coefficients array.
std::vector< double > rrouts(double x_1, double x_2, double eta1, double eta2, double eta3, double epsilon, double beta, double mu, double Alpha_star, double a, DxxParameters_structure DxxParameters)
routs finding routine
Matrix3D times(const Matrix3D< T > &M) const
Arraywise multiplication (A.*B), stores result in a new matrix.
string numberDensity
Number density model. Check StrToVal(string input, NumberDensities &place) for known values...
double mirror_point_coeff
Coefficient for mirror point (like 0.999, cause we can not integrate all the way to mirror point) ...
double Df_Ozeke(double L, double Kp)
Matrix1D< double > Lpp
Lpp array.
bool LoadDiffusionCoefficientFromFileWithGrid(GridElement &L, GridElement &pc, GridElement &alpha, string D_filename, string gridOrder="IFT_GRID_LPA")
Load Dxx from the file with Grid.
string DLLType
DLL type (which method to use to calculate). Check StrToVal(string input, DLLTypes &place) for known ...
double B(double lambda, double L)
double int_Daa_loc(double lambda, DxxParameters_structure DxxParameters)
string DxxType
Dxx type. Check StrToVal(string input, DiffusionCoefficientTypes &place) for known values...
double bisection(double(*f)(double x, double Alpha), double Alpha, double x0, double x1, double eps)
Looking for roots of function f for given Alpha on the interval x0-x1.
Main parameters structure that holds smaller structures for individual parameters.
bool useEnergyDiffusion
Using diffusions flags.
int size
size of grid element
bool useEnergyAlphaMixedTerms
double Df_Ozeke_E(double L, double Kp)
double int_Dpp_loc(double lambda, DxxParameters_structure DxxParameters)
Matrix3D< double > arr
array of diffusion coefficients
double Dpa_root(double Omega_e, double x, double mu, double su, double y, double beta, double a, double b, double Alpha_star, double s, double epsilon, double d_x, double x_m, double R, DxxParameters_structure DxxParameters)
double Dfe(double alpha, double Ke, double L, double Kp)
Electrostatic radial diffusion coefficient.
string waveName
Wave name.
void CreateAllDiffusionCoefficients(DiffusionCoefficient &DLL, DiffusionCoefficientsGroup &Daa, DiffusionCoefficientsGroup &Dpcpc, DiffusionCoefficientsGroup &Dpca, DiffusionCoefficientsGroup &DaaLpp, DiffusionCoefficientsGroup &DpcpcLpp, DiffusionCoefficientsGroup &DpcaLpp, Parameters_structure ¶meters, Grid &radialDiffusionGrid, Grid &localDiffusionsGrid)
void Output1DHeaders(ofstream &output1D, DiffusionCoefficientsGroup &Daa, DiffusionCoefficientsGroup &Dpcpc, DiffusionCoefficientsGroup &Dpca, DiffusionCoefficientsGroup &DaaLpp, DiffusionCoefficientsGroup &DpcpcLpp, DiffusionCoefficientsGroup &DpcaLpp)
Function used in main() to print header in output file. Header needs to understand the output file...
void get_quads(double *a, int n, double *quad, double *x)
void MakeDLL_B(GridElement &L, GridElement &pc, GridElement &alpha, double Kp)
Making DLL.
double multiplicator
Dxx is multiplied to this number after all other operations. Useful for fast Dxx modifications.
void MakeDLL_Ozeke(GridElement &L, GridElement &pc, GridElement &alpha, double Kp)
Making DLL_M.
GridElement L
grid elements to be used
Parameters_structure parameters
Parameters structure, with all parameters from the parameters.ini file. The default parameters define...
bool BwFromLambda
Bw lambda dependance flag.
string type
Type of the diffusion coefficient: Daa, Dpp, Dpa etc... Described in types.h file as an enumeration...
GridElement epc
grid elements to be used
Error message - stack of single_errors.
void MakeDLL_BE(GridElement &L, GridElement &pc, GridElement &alpha, double Kp)
Making DLL.
double Dxx_ba(double L1, double EMeV, double Alpha, double int_Dxx_loc(double lambda, DxxParameters_structure DxxParameters), DxxParameters_structure DxxParameters)
double Omega_m
Omega_m value.
Matrix1D< double > tau
Lifetime out of the plasmasphere.
double lam_max
Maximum latitude.
double omega_lc
omega lower cutoff
Computational grid composed of 3 different GridElement.
double omega_uc
omega upper cutoff
bool useRadialDiffusion
Using diffusions flags.
static const double B_0
magnetic field at the equatorial plane on the Earth surface, Gauss
double F_cap2(double x, double y, double a, double beta, double mu, double s, double epsilon, double Alpha_star, DxxParameters_structure DxxParameters)
static const double qe
electron charge, CGS units
void MakeDLL(double Kp)
Making DLL.
double Kfunc(double pc)
Computation of Kinetic energy from given momentum .
double Scale(double Kp)
Function, that scale diffusion coefficients.
GridElement pc
grid elements to be used
vector< DxxParameters_structure > DxxParametersList
List of diffusion coefficients parameters.
vector< DiffusionCoefficient > DxxList
List of diffusion coefficients in that group. Actually, it's a list of waves used in the diffusion co...