14 #include "../VariousFunctions/bisection.h"
16 #include "../Logging/Output.h"
17 #include "../Exceptions/error.h"
23 #define double_zero 1.e-21
24 #define min_Dxx 1.e-21
26 #if defined(_WIN32) || defined(_WIN64)
27 #define strncasecmp _strnicmp
28 #define strcasecmp _stricmp
40 bool Load_Failed =
false;
42 if (!arr.initialized) {
50 this->DxxParameters = DxxParameters;
52 this->type = DxxParameters.
DxxType;
53 this->useScale = DxxParameters.
useScale;
56 throw error_msg(
"DXX_GET",
"Undefined load or calculate choice.");
70 err.
add(
"UNKNOWN_ERROR",
"Unknown error while loading diffusion coefficient file");
78 Output::echo(3,
"Calculating %s...\n", arr.name.c_str());
79 Calculate(grid.
L, grid.
epc, grid.
alpha, DxxParameters);
81 arr = arr * (60.0 * 60 * 24);
103 if (this->type ==
"DCT_Dpp" || this->type ==
"DCT_DppLpp") arr.times_equal(grid.
pc.
arr.
times(grid.
pc.
arr));
104 if (this->type ==
"DCT_Dpa" || this->type ==
"DCT_DpaLpp") arr.times_equal(grid.
pc.
arr);
106 arr.change_ind = clock();
118 bool positivonly=
true;
122 if (DxxParameters.
DxxType ==
"DCT_Daa" || DxxParameters.
DxxType ==
"DCT_DaaLpp") {
124 }
else if (DxxParameters.
DxxType ==
"DCT_Dpp" || DxxParameters.
DxxType ==
"DCT_DppLpp") {
126 }
else if (DxxParameters.
DxxType ==
"DCT_Dpa" || DxxParameters.
DxxType ==
"DCT_DpaLpp") {
130 throw error_msg(
"DXX_TYPE",
"Unknown diffusion coefficient");
135 for (il = 0; il < L.
size ; il++) {
136 for (im = 0; im < epc.
size; im++) {
137 for (ia = 0; ia < Alpha.
size; ia++) {
140 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;
142 if (positivonly && (arr[il][im][ia] <
min_Dxx)) {
144 }
else if (!(arr[il][im][ia] >= 0) && !(arr[il][im][ia] <= 0)) {
163 if (!arr.initialized) {
164 throw error_msg(
"DXX_UNINITIALIZED_MATRIX",
"Matrix has not been initialized.");
169 if (filetype ==
"IFT_GRID") {
170 LoadDiffusionCoefficientFromFileWithGrid(L, pc, alpha, D_filename);
171 }
else if (filetype ==
"IFT_GRID_LPA") {
172 LoadDiffusionCoefficientFromFileWithGrid(L, pc, alpha, D_filename,
"IFT_GRID_LPA");
173 }
else if (filetype ==
"IFT_GRID_APL") {
174 LoadDiffusionCoefficientFromFileWithGrid(L, pc, alpha, D_filename,
"IFT_GRID_APL");
175 }
else if (filetype ==
"IFT_PLANE") {
177 LoadDiffusionCoefficientFromPlaneFile(L, pc, alpha, D_filename);
179 throw error_msg(
"DXX_LOAD_TYPE_ERR",
"Unknown coefficients file format %s", filetype.c_str());
192 ifstream input_D(D_filename.c_str());
194 if (input_D != NULL && !input_D.eof()) {
195 for (il = 0; il < L.
size; il++)
196 for (im = 0; im < pc.
size; im++)
197 for (ia = 0; ia < alpha.
size; ia++) {
199 input_D >> arr[il][im][ia];
200 if (input_D == NULL) {
201 throw error_msg(
"LOAD_ERROR",
"Wrong coefficient file %s (to few values)\n", D_filename.c_str());
205 throw error_msg(
"DXX_LOAD_GRID_ERR",
"error reading input file %s\n", D_filename.c_str());
219 double Load_L, Load_epc, Load_alpha;
224 ifstream input_D(D_filename.c_str());
226 if (input_D == NULL) {
227 throw error_msg(
"DXX_FILE_OPEN_ERR",
"Diffusion coefficient file %s not found.", D_filename.c_str());
232 while (strcasecmp(tmp.c_str(),
"ZONE") != 0 && !input_D.eof() ) {
236 throw error_msg(
"DXX_LOAD_GRID_ERR",
"Keyword 'ZONE' not found in the diffusion coefficient file: %s\n", D_filename.c_str());
240 input_D.ignore(9999,
'\n');
242 if (!input_D.eof()) {
244 if (gridOrder ==
"IFT_GRID_LPA") {
245 for (il = 0; il < L.
size; il++)
246 for (im = 0; im < epc.
size; im++)
247 for (ia = 0; ia < alpha.
size; ia++) {
248 input_D >> Load_L >> Load_epc >> Load_alpha;
251 input_D >> arr[il][im][ia];
253 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) {
254 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]);
257 }
else if (gridOrder ==
"IFT_GRID_APL") {
258 for (ia = 0; ia < alpha.
size; ia++)
259 for (im = 0; im < epc.
size; im++)
260 for (il = 0; il < L.
size; il++) {
261 input_D >> Load_alpha >> Load_epc >> Load_L ;
264 input_D >> arr[il][im][ia];
266 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) {
267 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]);
271 throw error_msg(
"UNKNOWN_GRID_TYPE",
"Unknown grid order: %s", gridOrder.c_str());
282 if (DLLType ==
"DLLT_BE") {
283 MakeDLL_BE(L, pc, alpha, Kp);
285 else if (DLLType ==
"DLLT_B") {
286 MakeDLL_B(L, pc, alpha, Kp);
288 else if (DLLType ==
"DLLT_FAKE") {
289 MakeDLL_FAKE(L, pc, alpha, Kp);
291 else if (DLLType ==
"DLLT_Ozeke") {
292 MakeDLL_Ozeke(L, pc, alpha, Kp);
294 else if (DLLType ==
"DLLT_Ozeke_ME") {
295 MakeDLL_Ozeke_ME(L, pc, alpha, Kp);
298 throw error_msg(
"DLL_ERROR",
"Unknown DLL type");
300 arr.change_ind = clock();
308 for (il = 0; il < L.
size; il++) {
309 for (im = 0; im < pc.
size; im++) {
310 for (ia = 0; ia < alpha.
size; ia++) {
311 arr[il][im][ia] =
Df(L.
arr[il][im][ia], Kp);
322 for (il = 0; il < L.
size; il++) {
323 for (im = 0; im < pc.
size; im++) {
324 for (ia = 0; ia < alpha.
size; ia++) {
325 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);;
336 for (il = 0; il < L.
size; il++) {
337 for (im = 0; im < pc.
size; im++) {
338 for (ia = 0; ia < alpha.
size; ia++) {
345 arr[il][im][ia] = arr[il][im][ia] +
Df(L.
arr[il][im][ia], Kp);
358 for (il = 0; il < L.
size; il++) {
359 for (im = 0; im < pc.
size; im++) {
360 for (ia = 0; ia < alpha.
size; ia++) {
361 arr[il][im][ia] =
VF::max(
Df(L.
arr[il][im][ia], Kp), 1.0/100);
378 for (il = 0; il < L.
size; il++) {
379 for (im = 0; im < pc.
size; im++) {
380 for (ia = 0; ia < alpha.
size; ia++) {
397 for (il = 0; il < L.
size; il++) {
398 for (im = 0; im < pc.
size; im++) {
399 for (ia = 0; ia < alpha.
size; ia++) {
413 double scaleCoeff = 0;
417 bool need_changes =
false;
418 for (Dxx_it = 0; Dxx_it < DxxList.size(); Dxx_it++) {
419 if (time >= DxxList[Dxx_it].DxxParameters.time_start && time < DxxList[Dxx_it].DxxParameters.time_end) {
421 if (!DxxList[Dxx_it].is_active) {
426 if (DxxList[Dxx_it].useScale) {
432 if (DxxList[Dxx_it].is_active) {
438 if (need_changes ==
false)
return false;
442 CurrentDxx.arr.change_ind = clock();
444 for (Dxx_it = 0; Dxx_it < DxxList.size(); Dxx_it++) {
446 if (time >= DxxList[Dxx_it].DxxParameters.time_start && time < DxxList[Dxx_it].DxxParameters.time_end) {
448 if (!DxxList[Dxx_it].is_active) {
449 DxxList[Dxx_it].is_active =
true;
450 Output::echo(2,
"Activated %s.\n", DxxList[Dxx_it].arr.name.c_str());
453 if (DxxList[Dxx_it].useScale) {
455 scaleCoeff = DxxList[Dxx_it].Scale(Kp);
456 Output::echo(2,
"Bw scaled %s by %f.\n", DxxList[Dxx_it].arr.name.c_str(), scaleCoeff);
458 CurrentDxx.arr = CurrentDxx.arr + DxxList[Dxx_it].arr * scaleCoeff;
460 CurrentDxx.arr = CurrentDxx.arr + DxxList[Dxx_it].arr;
464 if (DxxList[Dxx_it].is_active) {
465 DxxList[Dxx_it].is_active =
false;
466 Output::echo(2,
"Deactivated %s.", DxxList[Dxx_it].arr.name.c_str());
477 double Bw2=0, Bw2_0=0;
478 if (!this->useScale)
return 1;
479 if (DxxParameters.waveType ==
"WT_CHORUS_DAY" || DxxParameters.waveType ==
"WT_CHORUS_NIGHT") {
481 if (DxxParameters.DxxKp <= 2.333) Bw2_0 = 2.0*pow(10, 0.73 + 0.91 * DxxParameters.DxxKp);
482 if (DxxParameters.DxxKp > 2.333) Bw2_0 = 2.0*pow(10, 2.5 + 0.18 * DxxParameters.DxxKp);
483 if (Kp <= 2.333) Bw2 = 2.0*pow(10, 0.73 + 0.91 * Kp);
484 if (Kp > 2.333) Bw2 = 2.0*pow(10, 2.5 + 0.18 * Kp);
494 }
else if (DxxParameters.waveType ==
"WT_HISS_PP" || DxxParameters.waveType ==
"WT_HISS_PLUME") {
497 return pow(Kp/DxxParameters.DxxKp, 2);
498 }
else if (DxxParameters.waveType ==
"WT_EMIC_PLUME" ) {
500 return pow(Kp/DxxParameters.DxxKp, 4);
502 }
else if (DxxParameters.waveType ==
"WT_LIGHTNING") {
503 return pow(Kp/DxxParameters.DxxKp, 2);
514 double Alpha_ne(
double pangle,
double lambda,
double L) {
515 return asin( sqrt(
f1(lambda)) * sin(pangle) );
518 double f1(
double lambda) {
519 return sqrt(4.0 - 3.0 * pow(cos(lambda), 2))/(pow(cos(lambda), 6));
522 double B(
double lambda,
double L) {
527 return pow(x, 6) + (3.0 * pow(sin(Alpha), 4)) * x - 4.0 * pow(sin(Alpha), 4);
531 double c1 = 2.0*s*(-1+epsilon);
532 double c2 = 1.0-4.0*epsilon+pow(epsilon,2);
533 double c3 = -s*(-1+epsilon)*(b+4.0*epsilon)/2.0;
534 double c4 = epsilon*(b+epsilon) ;
535 double g = pow(x,4) + c1*pow(x,3) + c2*pow(x,2) + c3*x + c4;
536 return y*(pow((x-s),2))*(pow((x+s*epsilon),2))/(x*g);
540 double up = 2.0*(x+a)/(beta*mu);
541 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))+
542 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)));
550 double h = (b - a)/M;
551 double s = func(a, DxxParameters);
554 for (k = 1; k <= (M-1); k++) {
556 s = s + func(x, DxxParameters);
558 s = h * (func(a, DxxParameters) + func(b, DxxParameters))/2.0 + h*s;
565 DxxParameters.
L = L1;
566 DxxParameters.
EMeV = EMeV;
567 DxxParameters.
Alpha = Alpha;
574 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));
581 if (DxxParameters.
Alpha*180/
VC::pi < 0.001) Daa = 0;
582 return Daa * cos(
Alpha_ne(DxxParameters.
Alpha, lambda, DxxParameters.
L)) * pow(cos(lambda), 7) /(pow(cos(DxxParameters.
Alpha), 2));
588 if (DxxParameters.
Alpha*180/
VC::pi < 0.001) Dpp = 0;
589 return Dpp * (cos(lambda) * pow((1 + 3*pow(sin(lambda), 2)), 0.5)) / cos(
Alpha_ne(DxxParameters.
Alpha, lambda, DxxParameters.
L));
595 if (DxxParameters.
Alpha*180/
VC::pi < 0.001) Dpa = 0;
596 return Dpa * sin(
Alpha_ne(DxxParameters.
Alpha, lambda, DxxParameters.
L)) * pow(cos(lambda), 7) / (cos(DxxParameters.
Alpha) * sin(DxxParameters.
Alpha));
599 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 ) {
611 Omega_p = DxxParameters.
f * Omega_e_eq;
613 N_0 = 124.0*(pow((3.0/DxxParameters.
L),4));
615 DxxParameters.
f = Omega_p / Omega_e_eq;
617 }
else if (DxxParameters.
numberDensity ==
"ND_CHORUS_NIGHT") {
619 N_0 = 124.0*(pow((3.0/DxxParameters.
L),4)) + 36.0*pow((3.0/DxxParameters.
L),3.5) * cos((LT-(7.7*pow((3.0/DxxParameters.
L),2.0)+12))*
VC::pi/12);
621 DxxParameters.
f = Omega_p / Omega_e_eq;
625 N_0 = 124.0*(pow((3.0/DxxParameters.
L),4)) + 36.0*pow((3.0/DxxParameters.
L),3.5) * cos((LT-(7.7*pow((3.0/DxxParameters.
L),2.0)+12))*
VC::pi/12);
627 DxxParameters.
f = Omega_p / Omega_e_eq;
630 N_0 = pow(10, (-0.3145*DxxParameters.
L + 3.9043));
632 DxxParameters.
f = Omega_p / Omega_e_eq;
635 N_0 = 1e3*pow(4.0/DxxParameters.
L, 4);
637 DxxParameters.
f = Omega_p / Omega_e_eq;
639 }
else if (DxxParameters.
numberDensity ==
"ND_PL_CA_STARKS") {
640 if (DxxParameters.
L>2) {
641 N_0 = pow(10, (-0.3145*DxxParameters.
L + 3.9043));
643 DxxParameters.
f = Omega_p / Omega_e_eq;
645 double N1, N2, L1, L2, A,
B;
646 N1=10-6.;N2=8.3-6.;L1=1.2;L2=5.9;
647 A=(N2-N1)/(L2-L1);B=N2-A*L2;
648 N_0 = pow(10, A*DxxParameters.
L+B);
650 DxxParameters.
f = Omega_p / Omega_e_eq;
654 throw error_msg(
"DXX_DENSITY_ERROR",
"Unknown number density type: %s\n", DxxParameters.
numberDensity.c_str());
658 double epsilon = 1.0/1837;
659 double beta = (pow((DxxParameters.
EMeV * (DxxParameters.
EMeV + 2)), 0.5)) / (DxxParameters.
EMeV + 1);
660 double mu = cos(
Alpha_ne(DxxParameters.
Alpha, lambda, DxxParameters.
L));
661 double su = sin(
Alpha_ne(DxxParameters.
Alpha, lambda, DxxParameters.
L));
664 double gamma = DxxParameters.
EMeV + 1;
665 double a = DxxParameters.
s * lam_sp / gamma;
667 double Alpha_star = pow((Omega_e / Omega_p), 2);
668 double b = (1 + epsilon)/Alpha_star;
678 Omega_m = DxxParameters.
Omega_m * Omega_e_eq;
679 d_omega = DxxParameters.
d_omega * Omega_e_eq ;
680 Omega_1 = DxxParameters.
omega_lc * Omega_e_eq;
681 Omega_2 = DxxParameters.
omega_uc * Omega_e_eq;
682 }
else if (DxxParameters.
Omega_mType ==
"OT_ABSOLUTE") {
688 throw error_msg(
"DXX_OMEGA_M_TYPE",
"Unknown Omega_m type: ", DxxParameters.
Omega_mType.c_str());
699 Bw = pow(10, 0.75 + 0.04*(lambda*180/
VC::pi))*1.e-8 * DxxParameters.
Bw;
701 Bw = DxxParameters.
Bw;
703 double R = pow((Bw /
B(lambda, DxxParameters.
L)), 2);
709 if (DxxParameters.
particle ==
"PT_ELECTRON") {
710 x_m = Omega_m / Omega_e;
711 d_x = d_omega / Omega_e;
712 x_1 = Omega_1 / Omega_e;
713 x_2 = Omega_2 / Omega_e;
714 }
else if (DxxParameters.
particle ==
"PT_PROTON") {
716 x_m = Omega_m / Omega_e * epsilon / cfm;
717 d_x = d_omega / Omega_e * epsilon / cfm;
718 x_1 = Omega_1 / Omega_e * epsilon / cfm;
719 x_2 = Omega_2 / Omega_e * epsilon / cfm;
722 std::vector<double> xx =
rrouts(x_1, x_2, DxxParameters.
eta1, DxxParameters.
eta2, DxxParameters.
eta3, epsilon, beta, mu, Alpha_star, a, DxxParameters);
725 double Dxx = 0, Dxx_tmp = 0;
726 for (
unsigned int i = 1; i <= xx.size();i++) {
728 y = (x + a) / (beta * mu);
729 Dxx_tmp = Dxx_root(Omega_e, x, mu, su, y, beta, a, b, Alpha_star, DxxParameters.
s, epsilon, d_x, x_m, R, DxxParameters);
732 if (xx.size() == 0) {
739 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) {
740 return (
VC::pi/2) * (1.0/DxxParameters.
nu) * Omega_e * (1.0/(pow((DxxParameters.
EMeV + 1), 2))) *
741 ((R * pow((1 - x * mu / (y * beta)), 2) * fabs(
F_cap2(x,y,a,beta,mu,s,epsilon,Alpha_star, DxxParameters)))/
742 (d_x * fabs(beta * mu -
F_cap2(x,y,a,beta,mu,s,epsilon,Alpha_star, DxxParameters))))*
743 exp( - (pow(((x - x_m) / d_x), 2)));
746 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) {
747 return -(
VC::pi/2)*(1.0/DxxParameters.
nu)* Omega_e*su/beta*(1.0/(pow((DxxParameters.
EMeV+1), 2)))*
748 ((R*(1-x/y*mu/beta)*(x/y)* fabs(
F_cap2(x,y,a,beta,mu,s,epsilon,Alpha_star,DxxParameters)))/
749 (d_x * fabs(beta * mu -
F_cap2(x,y,a,beta,mu,s,epsilon,Alpha_star, DxxParameters))))*
750 exp( -(pow(((x - x_m) / d_x), 2)));
753 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) {
754 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))) *
757 ((R * (pow(x/y, 2)) * fabs(
F_cap2(x,y,a,beta,mu,s,epsilon,Alpha_star,DxxParameters)))/
758 (d_x*fabs(beta *mu-
F_cap2(x,y,a,beta,mu,s,epsilon,Alpha_star, DxxParameters))))*
759 exp(-(pow(((x-x_m)/d_x),2)));
764 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) {
771 double c0=-DxxParameters.
s;
772 double c1=DxxParameters.
s*epsilon;
773 double c2=1.0/4.*DxxParameters.
s*epsilon;
774 double c3=1.0/16.*DxxParameters.
s*epsilon;
777 double d1=eta1*epsilon;
778 double d2=eta2/4.*epsilon;
779 double d3=eta3/16.*epsilon;
781 double g3=d0+d1+d2+d3;
782 double g2=d0*(c1+c2+c3)+d1*(c0+c2+c3)+d2*(c0+c1+c3)+d3*(c0+c1+c2);
783 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);
784 double g0=d0*c1*c2*c3+d1*c0*c2*c3+d2*c0*c1*c3+d3*c0*c1*c2;
787 double h3=c0+c1+c2+c3;
788 double h2=c0*c1+c0*c2+c0*c3+c1*c2+c1*c3+c2*c3;
789 double h1=c0*c1*c2+c0*c1*c3+c0*c2*c3+c1*c2*c3;
790 double h0=c0*c1*c2*c3;
793 A[0] = (pow(beta,2)*pow(mu,2)-1.)*h4;
794 A[1] = (pow(beta,2)*pow(mu,2)-1.)*h3-2*a*h4;
795 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;
796 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;
797 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;
798 A[5] = -2*a*h0-pow(beta,2)*pow(mu,2)*1./Alpha_star*g0-pow(a,2)*h1;
809 std::vector<double> xx;
810 double r[7], wr[6],wi[6],quad[2];
814 quad[0] = 2.71828e-1;
815 quad[1] = 3.14159e-1;
818 numr =
roots(r, n, wr, wi);
820 for (
int i = 0; i < numr; i++) {
821 if (fabs(wi[i]) <
double_zero && wr[i] < x_2 && wr[i] > x_1) {
832 double Df (
double L,
double Kp) {
833 return pow(10, 0.506*Kp-9.325)*pow(L,10);
842 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) );
848 return 2.16e-8 * pow(L,6) * pow(10,( 0.217 * L + 0.461 * Kp ));
855 // http://irbem.svn.sourceforge.net/viewvc/irbem/extras/opendc/
856 double Df_Ozeke(double L, double Kp) {
858 return 1.94e-6 * pow(L,8) * pow(10, pow(L,2) * 0.112 + L * Kp * 0.0824 - L * 1.32 + pow(Kp,2) * 0.0291 - Kp * 0.270 );
861 double Df_Ozeke_E(double L, double Kp) {
863 return 5.75e-9 * pow(L,6) * pow(10,( 0.208 * L + 0.457 * Kp ));
901 if (Dxx_tmp.
type ==
"DCT_Daa")
902 Daa.
DxxList.push_back(Dxx_tmp);
903 else if (Dxx_tmp.
type ==
"DCT_DaaLpp")
904 DaaLpp.
DxxList.push_back(Dxx_tmp);
908 if (Dxx_tmp.
type ==
"DCT_Dpp" || Dxx_tmp.
type ==
"DCT_Dpcpc")
909 Dpcpc.
DxxList.push_back(Dxx_tmp);
910 else if (Dxx_tmp.
type ==
"DCT_DppLpp" || Dxx_tmp.
type ==
"DCT_DpcpcLpp")
911 DpcpcLpp.
DxxList.push_back(Dxx_tmp);
915 if (Dxx_tmp.
type ==
"DCT_Dpca" || Dxx_tmp.
type ==
"DCT_Dpa")
916 Dpca.
DxxList.push_back(Dxx_tmp);
917 else if (Dxx_tmp.
type ==
"DCT_DpcaLpp" || Dxx_tmp.
type ==
"DCT_DpaLpp")
918 DpcaLpp.
DxxList.push_back(Dxx_tmp);
924 int Daa_it, Dpcpc_it, Dpca_it;
926 for (Dpca_it = 0; Dpca_it < Dpca.
DxxList.size(); Dpca_it++) {
929 for (Daa_it = 0; Daa_it < Daa.
DxxList.size(); Daa_it++) {
930 if (Daa.
DxxList[Daa_it].DxxParameters.waveType == Dpca.
DxxList[Dpca_it].DxxParameters.waveType)
934 for (Dpcpc_it = 0; Dpcpc_it < Dpcpc.
DxxList.size(); Dpcpc_it++) {
935 if (Dpcpc.
DxxList[Dpcpc_it].DxxParameters.waveType == Dpca.
DxxList[Dpca_it].DxxParameters.waveType)
940 for (iL = 0; iL < localDiffusionsGrid.
L.
size; iL++)
941 for (im = 0; im < localDiffusionsGrid.
pc.
size; im++)
942 for (ia = 0; ia < localDiffusionsGrid.
alpha.
size; ia++)
943 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])
944 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",
945 Dpca.
DxxList[Dpca_it].DxxParameters.waveType.c_str(),
947 localDiffusionsGrid.
L.
arr[iL][im][ia], localDiffusionsGrid.
epc.
arr[iL][im][ia], localDiffusionsGrid.
alpha.
arr[iL][im][ia],
948 Dpcpc.
DxxList[Dpcpc_it].arr[iL][im][ia], Daa.
DxxList[Daa_it].arr[iL][im][ia], Dpca.
DxxList[Dpca_it].arr[iL][im][ia]);
952 for (Dpca_it = 0; Dpca_it < DpcaLpp.
DxxList.size(); Dpca_it++) {
955 for (Daa_it = 0; Daa_it < DaaLpp.
DxxList.size(); Daa_it++) {
956 if (DaaLpp.
DxxList[Daa_it].DxxParameters.waveType == DpcaLpp.
DxxList[Dpca_it].DxxParameters.waveType)
960 for (Dpcpc_it = 0; Dpcpc_it < DpcpcLpp.
DxxList.size(); Dpcpc_it++) {
961 if (DpcpcLpp.
DxxList[Dpcpc_it].DxxParameters.waveType == DpcaLpp.
DxxList[Dpca_it].DxxParameters.waveType)
966 for (iL = 0; iL < localDiffusionsGrid.
L.
size; iL++)
967 for (im = 0; im < localDiffusionsGrid.
pc.
size; im++)
968 for (ia = 0; ia < localDiffusionsGrid.
alpha.
size; ia++)
969 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])
970 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",
971 DpcaLpp.
DxxList[Dpca_it].DxxParameters.waveType.c_str(),
973 localDiffusionsGrid.
L.
arr[iL][im][ia], localDiffusionsGrid.
epc.
arr[iL][im][ia], localDiffusionsGrid.
alpha.
arr[iL][im][ia],
974 DpcpcLpp.
DxxList[Dpcpc_it].arr[iL][im][ia], DaaLpp.
DxxList[Daa_it].arr[iL][im][ia], DpcaLpp.
DxxList[Dpca_it].arr[iL][im][ia]);
981 DLL.
MakeDLL(radialDiffusionGrid.
L, radialDiffusionGrid.
pc, radialDiffusionGrid.
alpha, parameters.
Kp[0], parameters.
DLLType);
994 output1D <<
"Variables = \"time\", \"Kp\", \"Boundary_fluxes\", \"Lpp\", \"tau\" ";
998 for (Dxx_it = 0; Dxx_it < Daa.
DxxList.size(); Dxx_it++) {
999 output1D <<
"\"" << Daa.
DxxList[Dxx_it].arr.name <<
"\"" <<
", ";
1001 for (Dxx_it = 0; Dxx_it < Dpcpc.
DxxList.size(); Dxx_it++) {
1002 output1D <<
"\"" << Dpcpc.
DxxList[Dxx_it].arr.name <<
"\"" <<
", ";
1004 for (Dxx_it = 0; Dxx_it < Dpca.
DxxList.size(); Dxx_it++) {
1005 output1D <<
"\"" << Dpca.
DxxList[Dxx_it].arr.name <<
"\"" <<
", ";
1008 for (Dxx_it = 0; Dxx_it < DaaLpp.
DxxList.size(); Dxx_it++) {
1009 output1D <<
"\"" << DaaLpp.
DxxList[Dxx_it].arr.name <<
"\"" <<
", ";
1011 for (Dxx_it = 0; Dxx_it < DpcpcLpp.
DxxList.size(); Dxx_it++) {
1012 output1D <<
"\"" << DpcpcLpp.
DxxList[Dxx_it].arr.name <<
"\"" <<
", ";
1014 for (Dxx_it = 0; Dxx_it < DpcaLpp.
DxxList.size(); Dxx_it++) {
1015 output1D <<
"\"" << DpcaLpp.
DxxList[Dxx_it].arr.name <<
"\"" <<
", ";
1019 output1D <<
"ZONE T=\"1d-output\" " << endl;
1028 output1D << time <<
"\t" << parameters.
Kp[iteration] <<
"\t" << parameters.
Bf[iteration] <<
"\t" << parameters.
Lpp[iteration] <<
"\t" << parameters.
tau[iteration] <<
"\t";
1030 unsigned int Dxx_it;
1031 for (Dxx_it = 0; Dxx_it < Daa.
DxxList.size(); Dxx_it++) {
1032 output1D << Daa.
DxxList[Dxx_it].Scale(parameters.
Kp[iteration]) <<
"\t";
1034 for (Dxx_it = 0; Dxx_it < Dpcpc.
DxxList.size(); Dxx_it++) {
1035 output1D << Dpcpc.
DxxList[Dxx_it].Scale(parameters.
Kp[iteration]) <<
"\t";
1037 for (Dxx_it = 0; Dxx_it < Dpca.
DxxList.size(); Dxx_it++) {
1038 output1D << Dpca.
DxxList[Dxx_it].Scale(parameters.
Kp[iteration]) <<
"\t";
1041 for (Dxx_it = 0; Dxx_it < DaaLpp.
DxxList.size(); Dxx_it++) {
1042 output1D << DaaLpp.
DxxList[Dxx_it].Scale(parameters.
Kp[iteration]) <<
"\t";
1044 for (Dxx_it = 0; Dxx_it < DpcpcLpp.
DxxList.size(); Dxx_it++) {
1045 output1D << DpcpcLpp.
DxxList[Dxx_it].Scale(parameters.
Kp[iteration]) <<
"\t";
1047 for (Dxx_it = 0; Dxx_it < DpcaLpp.
DxxList.size(); Dxx_it++) {
1048 output1D << DpcaLpp.
DxxList[Dxx_it].Scale(parameters.
Kp[iteration]) <<
"\t";