25 #include "../VariousFunctions/variousFunctions.h"
28 #include "../Exceptions/error.h"
48 if (!arr.initialized) {
63 if (grid.
L.
size != 1) {
64 throw error_msg(
"INITIAL_PSD_2D_METHOD_FOR_3D_GRID",
"Initial PSD error: imposibble to use ORBIT_FLUX_2D if L.size > 1");
70 if (grid.
L.
size != 1) {
71 throw error_msg(
"INITIAL_PSD_2D_METHOD_FOR_3D_GRID",
"Initial PSD error: imposibble to use ORBIT_FLUX_2D if L.size > 1");
73 cout<<
"kim................................"<<endl;
74 Load_initial_f_maxwell(grid.
L, grid.
pc, grid.
alpha);
80 if (grid.
L.
size == 1) {
81 throw error_msg(
"INITIAL_PSD_3D_METHOD_FOR_2D_GRID",
"Initial PSD error: imposibble to use STEADY_STATE if L.size == 1");
90 if (grid.
L.
size == 1) {
91 throw error_msg(
"INITIAL_PSD_3D_METHOD_FOR_2D_GRID",
"Initial PSD error: imposibble to use STEADY_STATE if L.size == 1");
96 else if (parameters.
initial_PSD_Type ==
"IPSDT_STEADY_STATE_FROM_BOUNDARY") {
98 if (grid.
L.
size == 1) {
99 throw error_msg(
"INITIAL_PSD_3D_METHOD_FOR_2D_GRID",
"Initial PSD error: imposibble to use STEADY_STATE if L.size == 1");
102 throw error_msg(
"INITIAL_PSD_ERROR",
"Initial PSD error: upper boundary is not initialized, can not be used for steady state calculation.");
106 Load_initial_f_from_outer_L(grid.
L, grid.
pc, grid.
alpha,
117 for (il = 0; il < grid.
L.
size; il++)
118 for (im = 0; im < grid.
pc.
size; im++)
131 throw error_msg(
"INITIAL_PSD_TYPE_UNKNOWN",
"Unknown initial PSD type");
138 for (il = 0; il < grid.
L.
size; il++)
139 for (im = 0; im < grid.
pc.
size; im++)
140 for (ia = 0; ia < grid.
alpha.
size; ia++) {
147 arr[il][im][ia] = arr[il][im][ia] * pow(sinh(grid.
alpha.
arr[il][im][ia]/
VF::alc(grid.
L.
arr[il][im][ia])),2)/pow(sinh(1.0),2);
157 output_without_grid_file <<
"VARIABLES = \"Phase_Space_Density\" " << endl;
158 output_without_grid_file.setf(ios_base::scientific, ios_base::floatfield);
222 if (parameters.initial_PSD_Type == "IPSDT_STEADY_STATE_TWO_ZONE") {
229 Lsize = (Lmax-Lmin)/dL+1;
230 //Lpp = (5.6 - 0.46*Kp);
232 for (il=0; il<Lsize; il++) L_1d[il] = 1+il*dL; // make 1d grid
234 for (im = 0; im < pc.size; im++){
235 for (il = 0; il < L.size; il++) {
237 Ke1[il] = VF::Kfunc(pc[il][im][alpha.size-1]); // Energy in MeV
238 alpha1 = alpha[il][im][alpha.size-1]; // ~90 deg alpha
239 // cout<<Ke1[il]<<" "<<il<<Lsize<<endl;
240 // make a array of tau
241 //if(L_1d[il] >= Lpp) { // calculate tau
243 //} else { // calculate tauLpp
245 tau_0[il] = VF::Coulomb(L_1d[il],Ke1[il]);
246 } else if (tauLpp == 2e99) {
247 tau_0[il] = VF::Precipitation(Kp,L_1d[il],Ke1[il]);
248 } else if (tauLpp == 3e99) {
249 tau_c = VF::Coulomb(L_1d[il],Ke1[il]);
250 tau_p = VF::Precipitation(Kp,L_1d[il],Ke1[il]);
251 tau_0[il] = 1.0 / (1.0/tau_p+1.0/tau_c);
255 //cout<<tau_0[il]<<" "<<il<<endl;
258 steady_state_two_zone(fL, tau_0, Kp, alpha1, Ke1, L.size, L_1d, fb_out, fb_in);
260 L_upper = L[L.size-1][0][0]; // L-upper boundary we want to change, which is given in parameters.ini file
262 if(L_upper != Lmax) { // if L is not equal to 7
263 unsigned i=0; // interpolate PSD using Lagragian
264 while (L_1d[++i] <= L_upper);
265 scale_fm[il][im] = fL[i-1] * ( (L_upper-L_1d[i])/(L_1d[i-1]-L_1d[i]) )
266 + fL[i] * ( (L_upper-L_1d[i-1])/(L_1d[i]-L_1d[i-1]) );
267 } else { // if L is equal to 7
268 scale_fm[il][im] = fb_out; // scale_f = 1;
271 // scaling PSD ///////////////////////
272 for (il = 0; il < L.size; il++) {
273 for (im = 0; im < pc.size; im++) {
274 for (ia = 0; ia < alpha.size; ia++) {
275 arr[il][im][ia] = arr[il][im][ia] * scale_fm[il][im];
303 double dt,
double Lpp,
304 double tau,
double tauLpp,
double Kp)
309 double tau_p,tau_c,tauLpp0;
312 for (il = 0; il < L.
size; il++) {
313 for (im = 0; im < pc.
size; im++) {
314 for (ia = 0; ia < alpha.
size; ia++) {
332 if (L.
arr[il][im][ia] >= Lpp) {
334 arr[il][im][ia] = arr[il][im][ia] - arr[il][im][ia] * (dt / tau);
340 }
else if (parameters.
usetauLpp ==
"precipitaion") {
342 }
else if (parameters.
usetauLpp ==
"combined") {
345 tauLpp0 = 1.0/(1.0/tau_p+1.0/tau_c);
349 arr[il][im][ia] = arr[il][im][ia] - arr[il][im][ia] * (dt / tauLpp0);
386 string pc_lowerBoundaryCondition_calculationType,
387 string pc_upperBoundaryCondition_calculationType,
388 string alpha_lowerBoundaryCondition_calculationType,
389 string alpha_upperBoundaryCondition_calculationType)
401 double Dpca_, Dpca_a_L, Dpca_a_R, Dpca_m_L, Dpca_m_R;
402 double J_, J_a_L, J_a_R, J_m_L, J_m_R;
403 double Laml, Lmal, Lamr, Lmar;
404 double PSD_m_r, PSD_m_l, PSD_a_r, PSD_a_l, Lam, Lma;
407 if (PSD_parameters.approximationMethod ==
"AM_Split_LR") {
410 for (il = 0; il < L.
size; il++) {
411 for (im = 0; im < pc.
size; im++) {
412 for (ia = 1; ia < alpha.
size-1; ia++) {
413 PSD_Der_a_L[il][im][ia] =
VF::G(alpha.
arr[il][im][ia]) * Dpca.
arr[il][im][ia]*(arr[il][im][ia] - arr[il][im][ia-1])/(alpha.
arr[il][im][ia] - alpha.
arr[il][im][ia-1]);
414 PSD_Der_a_R[il][im][ia] =
VF::G(alpha.
arr[il][im][ia]) * Dpca.
arr[il][im][ia]*(arr[il][im][ia+1] - arr[il][im][ia])/(alpha.
arr[il][im][ia+1] - alpha.
arr[il][im][ia]);
417 for (im = 1; im < pc.
size-1; im++) {
418 for (ia = 0; ia < alpha.
size; ia++) {
419 PSD_Der_m_L[il][im][ia] = pc.
arr[il][im][ia] * pc.
arr[il][im][ia] * Dpca.
arr[il][im][ia]*(arr[il][im][ia] - arr[il][im-1][ia])/(pc.
arr[il][im][ia] - pc.
arr[il][im-1][ia]);
420 PSD_Der_m_R[il][im][ia] = pc.
arr[il][im][ia] * pc.
arr[il][im][ia] * Dpca.
arr[il][im][ia]*(arr[il][im+1][ia] - arr[il][im][ia])/(pc.
arr[il][im+1][ia] - pc.
arr[il][im][ia]);
426 for (il = 0; il < L.
size; il++) {
427 for (im = 1; im < pc.
size-1; im++) {
428 for (ia = 1; ia < alpha.
size-1; ia++) {
430 Dpca_ = Dpca.
arr[il][im][ia];
431 Dpca_a_L = Dpca.
arr[il][im][ia-1];
432 Dpca_a_R = Dpca.
arr[il][im][ia+1];
433 Dpca_m_L = Dpca.
arr[il][im-1][ia];
434 Dpca_m_R = Dpca.
arr[il][im+1][ia];
435 J_ = Jacobian[il][im][ia];
436 J_a_L = Jacobian[il][im][ia-1];
437 J_a_R = Jacobian[il][im][ia+1];
438 J_m_L = Jacobian[il][im-1][ia];
439 J_m_R = Jacobian[il][im+1][ia];
443 Lmar = 1.0/(pc.
arr[il][im][ia] * pc.
arr[il][im][ia]) * (PSD_Der_a_R[il][im][ia] - PSD_Der_a_R[il][im-1][ia])/(pc.
arr[il][im][ia] - pc.
arr[il][im-1][ia]);
444 Lmal = 1.0/(pc.
arr[il][im][ia] * pc.
arr[il][im][ia]) * (PSD_Der_a_L[il][im+1][ia] - PSD_Der_a_L[il][im][ia])/(pc.
arr[il][im+1][ia] - pc.
arr[il][im][ia]);
447 Lamr = 1.0/(
VF::G(alpha.
arr[il][im][ia])) * (PSD_Der_m_R[il][im][ia] - PSD_Der_m_R[il][im][ia-1])/(alpha.
arr[il][im][ia] - alpha.
arr[il][im][ia-1]);
448 Laml = 1.0/(
VF::G(alpha.
arr[il][im][ia])) * (PSD_Der_m_L[il][im][ia+1] - PSD_Der_m_L[il][im][ia])/(alpha.
arr[il][im][ia+1] - alpha.
arr[il][im][ia]);
450 arr[il][im][ia] = arr[il][im][ia] + dt * 0.5 * (Lmar + Lamr + Lmal + Laml);
464 }
else if (PSD_parameters.approximationMethod ==
"AM_Split_C") {
465 for (il = 0; il < L.
size; il++) {
466 for (im = 1; im < pc.
size-1; im++) {
467 for (ia = 1; ia < alpha.
size-1; ia++) {
469 Dpca_ = Dpca.
arr[il][im][ia];
470 Dpca_a_L = Dpca.
arr[il][im][ia-1];
471 Dpca_a_R = Dpca.
arr[il][im][ia+1];
472 Dpca_m_L = Dpca.
arr[il][im-1][ia];
473 Dpca_m_R = Dpca.
arr[il][im+1][ia];
474 J_ = Jacobian[il][im][ia];
475 J_a_L = Jacobian[il][im][ia-1];
476 J_a_R = Jacobian[il][im][ia+1];
477 J_m_L = Jacobian[il][im-1][ia];
478 J_m_R = Jacobian[il][im+1][ia];
481 PSD_m_r = Dpca.
arr[il][im+1][ia]*Jacobian[il][im+1][ia]*(arr[il][im+1][ia+1] - arr[il][im+1][ia-1])/(alpha.
arr[il][im+1][ia+1] - alpha.
arr[il][im+1][ia-1]);
482 PSD_m_l = Dpca.
arr[il][im-1][ia]*Jacobian[il][im-1][ia]*(arr[il][im-1][ia+1] - arr[il][im-1][ia-1])/(alpha.
arr[il][im-1][ia+1] - alpha.
arr[il][im-1][ia-1]);
483 Lma = 1.0/Jacobian[il][im][ia] * (PSD_m_r - PSD_m_l) / (pc.
arr[il][im+1][ia] - pc.
arr[il][im-1][ia]);
486 PSD_a_r = Dpca.
arr[il][im][ia+1]*Jacobian[il][im][ia+1]*(arr[il][im+1][ia+1] - arr[il][im-1][ia+1])/(pc.
arr[il][im+1][ia+1] - pc.
arr[il][im-1][ia+1]);
487 PSD_a_l = Dpca.
arr[il][im][ia-1]*Jacobian[il][im][ia-1]*(arr[il][im+1][ia-1] - arr[il][im-1][ia-1])/(pc.
arr[il][im+1][ia-1] - pc.
arr[il][im-1][ia-1]);
488 Lam = 1.0/Jacobian[il][im][ia] * (PSD_a_r - PSD_a_l) / (alpha.
arr[il][im][ia+1] - alpha.
arr[il][im][ia-1]);
491 arr[il][im][ia] = arr[il][im][ia] + dt * (Lma + Lam);
497 throw error_msg(
"APPROXIMATION",
"Unknown approximation method for mixed terms");
527 string alpha_lowerBoundaryCondition_calculationType,
528 string alpha_upperBoundaryCondition_calculationType)
531 char out_name_matr_a[256], out_name_matr_b[256], out_name_matr_c[256];
541 matr_A.Initialize(alpha.
size, 1, 1, 1);
542 matr_B.Initialize(alpha.
size, 1, 1, 0),
543 matr_C.Initialize(alpha.
size, 1, 1, 0);
560 static Matrix3D<double> Daa_1d(1, 1, alpha.
size), DaaLpp_1d(1, 1, alpha.
size), Jacobian_1d(1, 1, alpha.
size), alpha_1d(1, 1, alpha.
size), L_1d(1, 1, alpha.
size), pc_1d(1, 1, alpha.
size);
562 static Matrix2D<double> alpha_lowerBoundaryCondition_1d(1, 1), alpha_upperBoundaryCondition_1d(1, 1);
566 DiagMatrix::iterator it;
567 for (il = 0; il < L.
size; il++) {
568 for (im = 0; im < pc.
size; im++) {
571 for (ia = 0; ia < alpha.
size; ia++) {
572 slicePSD1D[ia] = arr[il][im][ia];
574 Daa_1d[0][0][ia] = Daa.
arr[il][im][ia];
575 DaaLpp_1d[0][0][ia] = DaaLpp.
arr[il][im][ia];
577 L_1d[0][0][ia] = L.
arr[il][im][ia];
578 pc_1d[0][0][ia] = pc.
arr[il][im][ia];
579 alpha_1d[0][0][ia] = alpha.
arr[il][im][ia];
581 Jacobian_1d[0][0][ia] = Jacobian[il][im][ia];
584 alpha_upperBoundaryCondition_1d[0][0] = alpha_upperBoundaryCondition[il][im];
585 alpha_lowerBoundaryCondition_1d[0][0] = alpha_lowerBoundaryCondition[il][im];
589 matr_A, matr_B, matr_C,
590 L_1d, pc_1d, alpha_1d,
592 Zero_Matrix_2d, Zero_Matrix_2d,
593 Zero_Matrix_2d, Zero_Matrix_2d,
594 alpha_lowerBoundaryCondition_1d, alpha_upperBoundaryCondition_1d,
595 "BCT_CONSTANT_VALUE",
"BCT_CONSTANT_VALUE",
596 "BCT_CONSTANT_VALUE",
"BCT_CONSTANT_VALUE",
597 alpha_lowerBoundaryCondition_calculationType, alpha_upperBoundaryCondition_calculationType,
599 Zero_Matrix, Zero_Matrix,
601 Zero_Matrix, Zero_Matrix,
602 Jacobian_1d, dt, Lpp);
614 for (in = 0; in < alpha.
size; in++) {
617 RHS[in] = matr_C[0][in];
618 for (it = matr_B.begin(); it != matr_B.end(); it++) {
620 if (in + it->first >= 0 && in + it->first < matr_B.total_size) {
621 RHS[in] += it->second[in] * slicePSD1D[in + it->first];
641 for (ia = 0; ia < alpha.
size; ia++) {
642 arr[il][im][ia] = slicePSD1D[ia];
673 string pc_lowerBoundaryCondition_calculationType,
674 string pc_upperBoundaryCondition_calculationType)
676 char out_name_matr_a[256], out_name_matr_b[256], out_name_matr_c[256];
687 matr_A.Initialize(1, pc.
size, 1, 1);
688 matr_B.Initialize(1, pc.
size, 1, 0),
689 matr_C.Initialize(1, pc.
size, 1, 0);
705 static Matrix3D<double> Dpcpc_1d(1, pc.
size, 1), DpcpcLpp_1d(1, pc.
size, 1), Jacobian_1d(1, pc.
size, 1), L_1d(1, pc.
size, 1), pc_1d(1, pc.
size, 1), alpha_1d(1, pc.
size, 1);
707 static Matrix2D<double> pc_lowerBoundaryCondition_1d(1, 1), pc_upperBoundaryCondition_1d(1, 1);
711 DiagMatrix::iterator it;
712 for (il = 0; il < L.
size; il++) {
713 for (ia = 0; ia < alpha.
size; ia++) {
716 for (im = 0; im < pc.
size; im++) {
718 slicePSD1D[im] = arr[il][im][ia];
720 Dpcpc_1d[0][im][0] = Dpcpc.
arr[il][im][ia];
721 DpcpcLpp_1d[0][im][0] = DpcpcLpp.
arr[il][im][ia];
723 L_1d[0][im][0] = L.
arr[il][im][ia];
724 pc_1d[0][im][0] = pc.
arr[il][im][ia];
725 alpha_1d[0][im][0] = alpha.
arr[il][im][ia];
727 Jacobian_1d[0][im][0] = Jacobian[il][im][ia];
730 pc_upperBoundaryCondition_1d[0][0] = pc_upperBoundaryCondition[il][ia];
731 pc_lowerBoundaryCondition_1d[0][0] = pc_lowerBoundaryCondition[il][ia];
735 matr_A, matr_B, matr_C,
736 L_1d, pc_1d, alpha_1d,
738 Zero_Matrix_2d, Zero_Matrix_2d,
739 pc_lowerBoundaryCondition_1d, pc_upperBoundaryCondition_1d,
740 Zero_Matrix_2d, Zero_Matrix_2d,
741 "BCT_CONSTANT_VALUE",
"BCT_CONSTANT_VALUE",
742 pc_lowerBoundaryCondition_calculationType, pc_upperBoundaryCondition_calculationType,
743 "BCT_CONSTANT_VALUE",
"BCT_CONSTANT_VALUE",
745 Dpcpc_1d, DpcpcLpp_1d,
746 Zero_Matrix, Zero_Matrix,
747 Zero_Matrix, Zero_Matrix,
748 Jacobian_1d, dt, Lpp);
760 DiagMatrix::iterator it;
761 for (in = 0; in < pc.
size; in++) {
764 RHS[in] = matr_C[0][in];
765 for (it = matr_B.begin(); it != matr_B.end(); it++) {
767 if (in + it->first >= 0 && in + it->first < matr_B.total_size) {
768 RHS[in] += it->second[in] * slicePSD1D[in + it->first];
788 for (im = 0; im < pc.
size; im++) {
789 arr[il][im][ia] = slicePSD1D[im];
819 string L_lowerBoundaryCondition_calculationType,
820 string L_upperBoundaryCondition_calculationType,
821 double tau,
double tauLpp)
823 char out_name_matr_a[256], out_name_matr_b[256], out_name_matr_c[256];
834 matr_A.Initialize(1, 1, L.
size, 1);
835 matr_B.Initialize(1, 1, L.
size, 0),
836 matr_C.Initialize(1, 1, L.
size, 0);
854 static Matrix2D<double> L_lowerBoundaryCondition_1d(1, 1), L_upperBoundaryCondition_1d(1, 1);
858 DiagMatrix::iterator it;
859 for (ia = 0; ia < alpha.
size; ia++) {
860 for (im = 0; im < pc.
size; im++) {
863 for (il = 0; il < L.
size; il++) {
865 slicePSD1D[il] = arr[il][im][ia];
867 DLL_1d[il][0][0] = DLL.
arr[il][im][ia];
869 L_1d[il][0][0] = L.
arr[il][im][ia];
873 Jacobian_1d[il][0][0] = Jacobian[il][im][ia];
876 L_upperBoundaryCondition_1d[0][0] = L_upperBoundaryCondition[im][ia];
877 L_lowerBoundaryCondition_1d[0][0] = L_lowerBoundaryCondition[im][ia];
881 matr_A, matr_B, matr_C,
882 L_1d, Zero_Matrix, Zero_Matrix,
884 L_lowerBoundaryCondition_1d, L_upperBoundaryCondition_1d,
885 Zero_Matrix_2d, Zero_Matrix_2d,
886 Zero_Matrix_2d, Zero_Matrix_2d,
887 L_lowerBoundaryCondition_calculationType, L_upperBoundaryCondition_calculationType,
888 "BCT_CONSTANT_VALUE",
"BCT_CONSTANT_VALUE",
889 "BCT_CONSTANT_VALUE",
"BCT_CONSTANT_VALUE",
891 Zero_Matrix, Zero_Matrix,
892 Zero_Matrix, Zero_Matrix,
893 Zero_Matrix, Zero_Matrix,
894 Jacobian_1d, dt, Lpp,
909 DiagMatrix::iterator it;
910 for (in = 0; in < L.
size; in++) {
913 RHS[in] = matr_C[0][in];
914 for (it = matr_B.begin(); it != matr_B.end(); it++) {
916 if (in + it->first >= 0 && in + it->first < matr_B.total_size) {
917 RHS[in] += it->second[in] * slicePSD1D[in + it->first];
937 for (il = 0; il < L.
size; il++) {
938 arr[il][im][ia] = slicePSD1D[il];
979 string pc_lowerBoundaryCondition_calculationType,
980 string pc_upperBoundaryCondition_calculationType,
981 string alpha_lowerBoundaryCondition_calculationType,
982 string alpha_upperBoundaryCondition_calculationType)
993 matr_A.Initialize(alpha.
size, pc.
size, 1, 1);
994 matr_B.Initialize(alpha.
size, pc.
size, 1, 0),
995 matr_C.Initialize(alpha.
size, pc.
size, 1, 0);
1007 static Matrix1D<double> slicePSD1D(matr_A.total_size), RHS(matr_A.total_size);
1010 static Matrix3D<double> Dpcpc_slice(1, pc.
size, alpha.
size), DpcpcLpp_slice(1, pc.
size, alpha.
size), Daa_slice(1, pc.
size, alpha.
size), DaaLpp_slice(1, pc.
size, alpha.
size), Dpca_slice(1, pc.
size, alpha.
size), DpcaLpp_slice(1, pc.
size, alpha.
size);
1011 static Matrix3D<double> L_slice(1, pc.
size, alpha.
size), pc_slice(1, pc.
size, alpha.
size), alpha_slice(1, pc.
size, alpha.
size), Jacobian_slice(1, pc.
size, alpha.
size);
1013 static Matrix2D<double> pc_lowerBoundaryCondition_slice(1, alpha.
size), pc_upperBoundaryCondition_slice(1, alpha.
size);
1014 static Matrix2D<double> alpha_lowerBoundaryCondition_slice(1, pc.
size), alpha_upperBoundaryCondition_slice(1, pc.
size);
1018 DiagMatrix::iterator it;
1019 for (il = 0; il < L.
size; il++) {
1021 for (im = 0; im < pc.
size; im++) {
1022 for (ia = 0; ia < alpha.
size; ia++) {
1024 in = matr_A.index1d(ia, im);
1025 slicePSD1D[in] = arr[il][im][ia];
1027 Dpcpc_slice[0][im][ia] = Dpcpc.
arr[il][im][ia];
1028 DpcpcLpp_slice[0][im][ia] = DpcpcLpp.
arr[il][im][ia];
1029 Daa_slice[0][im][ia] = Daa.
arr[il][im][ia];
1030 DaaLpp_slice[0][im][ia] = DaaLpp.
arr[il][im][ia];
1031 Dpca_slice[0][im][ia] = Dpca.
arr[il][im][ia];
1032 DpcaLpp_slice[0][im][ia] = DpcaLpp.
arr[il][im][ia];
1034 L_slice[0][im][ia] = L.
arr[il][im][ia];
1035 pc_slice[0][im][ia] = pc.
arr[il][im][ia];
1036 alpha_slice[0][im][ia] = alpha.
arr[il][im][ia];
1038 Jacobian_slice[0][im][ia] = Jacobian[il][im][ia];
1042 for (ia = 0; ia < alpha.
size; ia++) {
1043 pc_upperBoundaryCondition_slice[0][ia] = pc_upperBoundaryCondition[il][ia];
1044 pc_lowerBoundaryCondition_slice[0][ia] = pc_lowerBoundaryCondition[il][ia];
1046 for (im = 0; im < pc.
size; im++) {
1047 alpha_upperBoundaryCondition_slice[0][im] = alpha_upperBoundaryCondition[il][im];
1048 alpha_lowerBoundaryCondition_slice[0][im] = alpha_lowerBoundaryCondition[il][im];
1053 matr_A, matr_B, matr_C,
1054 L_slice, pc_slice, alpha_slice,
1056 Zero_Matrix_2d, Zero_Matrix_2d,
1057 pc_lowerBoundaryCondition_slice, pc_upperBoundaryCondition_slice,
1058 alpha_lowerBoundaryCondition_slice, alpha_upperBoundaryCondition_slice,
1059 "BCT_CONSTANT_VALUE",
"BCT_CONSTANT_VALUE",
1060 pc_lowerBoundaryCondition_calculationType, pc_upperBoundaryCondition_calculationType,
1061 alpha_lowerBoundaryCondition_calculationType, alpha_upperBoundaryCondition_calculationType,
1063 Dpcpc_slice, DpcpcLpp_slice,
1064 Daa_slice, DaaLpp_slice,
1065 Dpca_slice, DpcaLpp_slice,
1066 Jacobian_slice, dt, Lpp);
1075 DiagMatrix::iterator it;
1076 for (im = 0; im < pc.
size; im++) {
1077 for (ia = 0; ia < alpha.
size; ia++) {
1078 in = matr_B.index1d(ia, im);
1081 RHS[in] = matr_C[0][in];
1082 for (it = matr_B.begin(); it != matr_B.end(); it++) {
1084 if (in + it->first >= 0 && in + it->first < matr_B.total_size) {
1085 RHS[in] += it->second[in] * slicePSD1D[in + it->first];
1100 int modelMatrixLineNumber;
1102 double sum, error, error_max;
1103 static bool recalculate_A =
true;
1104 if (PSD_parameters.solutionMethod ==
"SM_Gauss") {
1106 if (recalculate_A) {
1118 for (it = matr_A.begin(); it != matr_A.end(); it++) {
1119 for (modelMatrixLineNumber = 0; modelMatrixLineNumber < n; modelMatrixLineNumber++) {
1122 if (modelMatrixLineNumber + it->first >= 0 && modelMatrixLineNumber + it->first < n) {
1126 A[modelMatrixLineNumber][modelMatrixLineNumber + it->first] = it->second[modelMatrixLineNumber];
1136 recalculate_A =
true;
1141 slicePSD1D[n-1] = B[n-1] / A[n-1][n-1];
1142 for (k = n - 2; k >= 0; k--) {
1144 for (l = k + 1; l < n; l++)
1145 sum += A[k][l] * slicePSD1D[l];
1146 slicePSD1D[k] = B[k] - sum;
1154 for (in = 0; in < n; in++) {
1156 for (it = matr_A.begin(); it != matr_A.end(); it++) {
1158 if (in + it->first >= 0 && in + it->first < n) {
1159 error -= it->second[in] * slicePSD1D[in + it->first];
1162 error_max =
VF::max( error_max, error );
1167 }
else if (PSD_parameters.solutionMethod ==
"SM_Relaxation") {
1170 Output::echo(5,
"Inverting A matrix (relaxation)... ");
1175 }
else if (PSD_parameters.solutionMethod ==
"SM_Lapack") {
1180 Lapack(matr_A, RHS, slicePSD1D);
1183 }
else if (PSD_parameters.solutionMethod ==
"SM_GMRES") {
1189 gmres_wrapout(matr_A, RHS, slicePSD1D, PSD_parameters.GMRES_parameters);
1192 }
else if (PSD_parameters.solutionMethod ==
"SM_GMRES2") {
1198 gmres2_wrapout(matr_A, RHS, slicePSD1D, PSD_parameters.GMRES_parameters);
1202 }
else if (PSD_parameters.solutionMethod ==
"SM_Tridiag") {
1204 throw error_msg(
"SOLUTIONERR",
"Tridiagonal method can't be used");
1206 throw error_msg(
"SOLUTIONERR",
"Unknown solution method");
1210 for (im = 0; im < pc.
size; im++) {
1211 for (ia = 0; ia < alpha.
size; ia++) {
1212 in = matr_A.index1d(ia, im);
1213 arr[il][im][ia] = slicePSD1D[in];
1230 string pc_lowerBoundaryCondition_calculationType,
1231 string pc_upperBoundaryCondition_calculationType,
1232 string alpha_lowerBoundaryCondition_calculationType,
1233 string alpha_upperBoundaryCondition_calculationType)
1244 matr_A.Initialize(alpha.
size, pc.
size, 1, 1);
1245 matr_B.Initialize(alpha.
size, pc.
size, 1, 0),
1246 matr_C.Initialize(alpha.
size, pc.
size, 1, 0);
1261 static Matrix1D<double> slicePSD1D(matr_A.total_size), RHS(matr_A.total_size);
1264 static Matrix3D<double> Dpcpc_slice(1, pc.
size, alpha.
size), DpcpcLpp_slice(1, pc.
size, alpha.
size), Daa_slice(1, pc.
size, alpha.
size), DaaLpp_slice(1, pc.
size, alpha.
size), Dpca_slice(1, pc.
size, alpha.
size), DpcaLpp_slice(1, pc.
size, alpha.
size);
1265 static Matrix3D<double> L_slice(1, pc.
size, alpha.
size), pc_slice(1, pc.
size, alpha.
size), alpha_slice(1, pc.
size, alpha.
size), Jacobian_slice(1, pc.
size, alpha.
size);
1267 static Matrix2D<double> pc_lowerBoundaryCondition_slice(1, alpha.
size), pc_upperBoundaryCondition_slice(1, alpha.
size);
1268 static Matrix2D<double> alpha_lowerBoundaryCondition_slice(1, pc.
size), alpha_upperBoundaryCondition_slice(1, pc.
size);
1272 DiagMatrix::iterator it;
1273 for (il = 0; il < L.
size; il++) {
1275 for (im = 0; im < pc.
size; im++) {
1276 for (ia = 0; ia < alpha.
size; ia++) {
1278 in = matr_A.index1d(ia, im);
1279 slicePSD1D[in] = arr[il][im][ia];
1281 Dpcpc_slice[0][im][ia] = Dpcpc.
arr[il][im][ia];
1282 DpcpcLpp_slice[0][im][ia] = DpcpcLpp.
arr[il][im][ia];
1283 Daa_slice[0][im][ia] = Daa.
arr[il][im][ia];
1284 DaaLpp_slice[0][im][ia] = DaaLpp.
arr[il][im][ia];
1285 Dpca_slice[0][im][ia] = Dpca.
arr[il][im][ia];
1286 DpcaLpp_slice[0][im][ia] = DpcaLpp.
arr[il][im][ia];
1288 L_slice[0][im][ia] = L.
arr[il][im][ia];
1289 pc_slice[0][im][ia] = pc.
arr[il][im][ia];
1290 alpha_slice[0][im][ia] = alpha.
arr[il][im][ia];
1292 Jacobian_slice[0][im][ia] = Jacobian[il][im][ia];
1296 for (ia = 0; ia < alpha.
size; ia++) {
1297 pc_upperBoundaryCondition_slice[0][ia] = pc_upperBoundaryCondition[il][ia];
1298 pc_lowerBoundaryCondition_slice[0][ia] = pc_lowerBoundaryCondition[il][ia];
1300 for (im = 0; im < pc.
size; im++) {
1301 alpha_upperBoundaryCondition_slice[0][im] = alpha_upperBoundaryCondition[il][im];
1302 alpha_lowerBoundaryCondition_slice[0][im] = alpha_lowerBoundaryCondition[il][im];
1307 matr_A, matr_B, matr_C,
1308 L_slice, pc_slice, alpha_slice,
1310 Zero_Matrix_2d, Zero_Matrix_2d,
1311 pc_lowerBoundaryCondition_slice, pc_upperBoundaryCondition_slice,
1312 alpha_lowerBoundaryCondition_slice, alpha_upperBoundaryCondition_slice,
1313 "BCT_CONSTANT_VALUE",
"BCT_CONSTANT_VALUE",
1314 pc_lowerBoundaryCondition_calculationType, pc_upperBoundaryCondition_calculationType,
1315 alpha_lowerBoundaryCondition_calculationType, alpha_upperBoundaryCondition_calculationType,
1317 Dpcpc_slice, DpcpcLpp_slice,
1318 Daa_slice, DaaLpp_slice,
1319 Dpca_slice, DpcaLpp_slice,
1320 Jacobian_slice, dt, Lpp);
1329 DiagMatrix::iterator it;
1330 for (im = 0; im < pc.
size; im++) {
1331 for (ia = 0; ia < alpha.
size; ia++) {
1332 in = matr_B.index1d(ia, im);
1335 RHS[in] = matr_C[0][in];
1336 for (it = matr_B.begin(); it != matr_B.end(); it++) {
1338 if (in + it->first >= 0 && in + it->first < matr_B.total_size) {
1339 RHS[in] += it->second[in] * slicePSD1D[in + it->first];
1354 int modelMatrixLineNumber;
1356 double sum, error, error_max;
1357 static bool recalculate_A =
true;
1358 if (PSD_parameters.solutionMethod ==
"SM_Gauss") {
1360 if (recalculate_A) {
1372 for (it = matr_A.begin(); it != matr_A.end(); it++) {
1373 for (modelMatrixLineNumber = 0; modelMatrixLineNumber < n; modelMatrixLineNumber++) {
1376 if (modelMatrixLineNumber + it->first >= 0 && modelMatrixLineNumber + it->first < n) {
1380 A[modelMatrixLineNumber][modelMatrixLineNumber + it->first] = it->second[modelMatrixLineNumber];
1390 recalculate_A =
true;
1395 slicePSD1D[n-1] = B[n-1] / A[n-1][n-1];
1396 for (k = n - 2; k >= 0; k--) {
1398 for (l = k + 1; l < n; l++)
1399 sum += A[k][l] * slicePSD1D[l];
1400 slicePSD1D[k] = B[k] - sum;
1408 for (in = 0; in < n; in++) {
1410 for (it = matr_A.begin(); it != matr_A.end(); it++) {
1412 if (in + it->first >= 0 && in + it->first < n) {
1413 error -= it->second[in] * slicePSD1D[in + it->first];
1416 error_max =
VF::max( error_max, error );
1421 }
else if (PSD_parameters.solutionMethod ==
"SM_Relaxation") {
1424 Output::echo(5,
"Inverting A matrix (relaxation)... ");
1429 }
else if (PSD_parameters.solutionMethod ==
"SM_Lapack") {
1434 Lapack(matr_A, RHS, slicePSD1D);
1437 }
else if (PSD_parameters.solutionMethod ==
"SM_GMRES") {
1443 gmres_wrapout(matr_A, RHS, slicePSD1D, PSD_parameters.GMRES_parameters);
1446 }
else if (PSD_parameters.solutionMethod ==
"SM_GMRES2") {
1452 gmres2_wrapout(matr_A, RHS, slicePSD1D, PSD_parameters.GMRES_parameters);
1456 }
else if (PSD_parameters.solutionMethod ==
"SM_Tridiag") {
1458 throw error_msg(
"SOLUTIONERR",
"Tridiagonal method can't be used");
1460 throw error_msg(
"SOLUTIONERR",
"Unknown solution method");
1464 for (im = 0; im < pc.
size; im++) {
1465 for (ia = 0; ia < alpha.
size; ia++) {
1466 in = matr_A.index1d(ia, im);
1467 arr[il][im][ia] = slicePSD1D[in];
1487 if (!(arr[il][im][ia] >= 0 || arr[il][im][ia] <= 0) || arr[il][im][ia] > maxNum || arr[il][im][ia] < -maxNum){
1488 throw error_msg(
"INTERPOLATION_ERROR",
"Interpolation error at [%d][%d][%d]", il, im, ia);
1502 if (!(arr[im] >= 0 || arr[im] <= 0) || arr[im] > maxNum || arr[im] < -maxNum){
1503 throw error_msg(
"INTERPOLATION_ERROR",
"Interpolation error at [%d][%d][%d]", il, im, ia);
1512 double lowerBoundary1d=0, upperBoundary1d=0;
1516 if (interpolationParameters_structure.
type ==
"IT_NONE") {
1518 }
else if (oldGrid.
type == newGrid.
type) {
1521 for (il = 0; il < oldGrid.
L.
size; il++) {
1522 for (ia = 0; ia < oldGrid.
alpha.
size; ia++) {
1523 for (im = 0; im < oldGrid.
pc.
size; im++) {
1525 arr[il][im][ia] = otherPSD.
arr[il][im][ia];
1541 for (il = 0; il < oldGrid.
L.
size; il++) {
1542 for (ia = 0; ia < oldGrid.
alpha.
size; ia++) {
1543 for (im = 0; im < oldGrid.
pc.
size; im++) {
1545 if (newGrid.
L.
arr[il][im][ia] != oldGrid.
L.
arr[il][im][ia] || newGrid.
alpha.
arr[il][im][ia] != oldGrid.
alpha.
arr[il][im][ia]) {
1547 throw error_msg(
"INTERPOLATION_ERROR",
"Interpolation error: grids must have the same L and Alpha");
1555 if (interpolationParameters_structure.
useLog ==
"Log_10") {
1558 matrix_array_1d_old[im] = log10(otherPSD.
arr[il][im][ia]);
1566 lowerBoundary1d = log10(
VF::max(newGrid_pc_lowerBoundaryCondition[il][ia],
VC::zero_f));
1567 upperBoundary1d = log10(
VF::max(newGrid_pc_upperBoundaryCondition[il][ia],
VC::zero_f));
1569 else if (interpolationParameters_structure.
useLog ==
"Log_2") {
1572 matrix_array_1d_old[im] = log(otherPSD.
arr[il][im][ia])/log(2.0);
1575 matrix_array_1d_old[im] = log(
VC::zero_f)/log(2.0);
1582 lowerBoundary1d = log(
VF::max(newGrid_pc_lowerBoundaryCondition[il][ia],
VC::zero_f)) / log(2.0);
1583 upperBoundary1d = log(
VF::max(newGrid_pc_upperBoundaryCondition[il][ia],
VC::zero_f)) / log(2.0);
1585 else if (interpolationParameters_structure.
useLog ==
"Log_E") {
1588 matrix_array_1d_old[im] = log(otherPSD.
arr[il][im][ia]);
1598 lowerBoundary1d = log(
VF::max(newGrid_pc_lowerBoundaryCondition[il][ia],
VC::zero_f));
1599 upperBoundary1d = log(
VF::max(newGrid_pc_upperBoundaryCondition[il][ia],
VC::zero_f));
1601 else if (interpolationParameters_structure.
useLog ==
"NoLog") {
1602 matrix_array_1d_old[im] = otherPSD.
arr[il][im][ia];
1603 lowerBoundary1d =
VF::max(newGrid_pc_lowerBoundaryCondition[il][ia],
VC::zero_f);
1604 upperBoundary1d =
VF::max(newGrid_pc_upperBoundaryCondition[il][ia],
VC::zero_f);
1607 throw error_msg(
"INTERPOLATION_ERROR",
"Usage of log(PSD) is not define properly: %s", interpolationParameters_structure.
useLog.c_str());
1616 grid_1d_old[im] = oldGrid.
pc.
arr[il][im][ia];
1617 grid_1d_new[im] = newGrid.
pc.
arr[il][im][ia];
1621 if (interpolationParameters_structure.
type ==
"IT_LINEAR") {
1623 matrix_array_1d_new.
Interpolate(matrix_array_1d_old, grid_1d_old, grid_1d_new);
1624 }
else if (interpolationParameters_structure.
type ==
"IT_POLYNOMIAL") {
1626 matrix_array_1d_new.
Polilinear(matrix_array_1d_old, grid_1d_old, grid_1d_new, lowerBoundary1d, upperBoundary1d);
1627 }
else if (interpolationParameters_structure.
type ==
"IT_SPLINE") {
1629 matrix_array_1d_new.
Spline(matrix_array_1d_old, grid_1d_old, grid_1d_new, lowerBoundary1d, upperBoundary1d, interpolationParameters_structure.
linearSplineCoef, interpolationParameters_structure.
maxSecondDerivative);
1630 }
else if (interpolationParameters_structure.
type ==
"IT_COPY") {
1632 matrix_array_1d_new = matrix_array_1d_old;
1634 throw error_msg(
"INTERPOLATION_TYPE_ERROR",
"Unknown interpolation type");
1637 for (im = 0; im < newGrid.
pc.
size; im++) {
1644 if (interpolationParameters_structure.
useLog ==
"Log_10") arr[il][im][ia] = pow(10, matrix_array_1d_new[im]);
1645 else if (interpolationParameters_structure.
useLog ==
"Log_2") arr[il][im][ia] = pow(2, matrix_array_1d_new[im]);
1646 else if (interpolationParameters_structure.
useLog ==
"Log_E") arr[il][im][ia] = pow(
VC::exp, matrix_array_1d_new[im]);
1647 else arr[il][im][ia] = matrix_array_1d_new[im];
1671 output_without_grid_file <<
"ZONE T=\"" <<
internal << time <<
"\" I=" << arr.size_z <<
", J=" << arr.size_y <<
", K=" << arr.size_x << endl;
1672 for (il = 0; il < arr.size_x; il++) {
1673 for (im = 0; im < arr.size_y; im++) {
1674 for (ia = 0; ia < arr.size_z; ia++) {
1675 output_without_grid_file << arr[il][im][ia] << endl;
1694 double Load_L=-1, Load_epc=-1, Load_alpha=-1;
1698 ifstream input_f(filename);
1699 if (input_f != NULL && !input_f.eof()) {
1701 getline(input_f, inBuf);
1702 getline(input_f, inBuf);
1703 for (il = 0; il < L.
size; il++) {
1704 for (im = 0; im < epc.
size; im++) {
1705 for (ia = 0; ia < alpha.
size; ia++) {
1706 if (withGrid) input_f >> Load_L >> Load_epc >> Load_alpha;
1707 input_f >> arr[il][im][ia];
1709 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) {
1710 throw error_msg(
"INITIAL_PSD_LOAD_GRID_ERR",
"Loading %s: grid mismatch.\nLoaded: L=%e, epc=%e, alpha=%e\nGrid: L=%e, epc=%e, alpha=%e\n", filename, Load_L, Load_epc, Load_alpha, L.
arr[il][im][ia], epc.
arr[il][im][ia], alpha.
arr[il][im][ia]);
1716 throw error_msg(
"INITIAL_PSD_ERROR",
"Error reading PSD input file %s.\n", filename);
1727 double p1, Lmax, Lmin, scale_factor, Lgrid;
1728 int Lsize_7, addGrid;
1736 int ia_90=alpha.
size-1;
1737 while (alpha.
arr[L.
size-1][0][ia_90] >
VC::pi/2 && ia_90 > 0) ia_90--;
1741 if (L.
arr[L.
size-1][0][ia_90] < 7) {
1745 Lmax = L.
arr[L.
size-1][0][ia_90];
1746 Lmin = L.
arr[0][0][ia_90];
1747 Lgrid = (Lmax-Lmin)/(L.
size-1);
1748 Lsize_7 = (7.0-Lmin)/Lgrid + 1;
1749 addGrid = Lsize_7 - L.
size;
1755 for (il = 0; il < L.
size; il++) {
1756 L_1d[il] = L.
arr[il][0][ia_90];
1758 for (il = L.
size; il < Lsize_7; il++) {
1759 L_1d[il] = Lmin+Lgrid*il;
1766 for (il = 0; il < L.
size; il++) L_1d[il] = L.
arr[il][0][ia_90];
1770 steady_state(fL, tau, Kp, Lsize_7, L_1d, fb_out, fb_in);
1775 if (L_1d[L_1d.
size_x-1] > 7) {
1776 int il_7 = L.
size-1;
1777 while (L_1d[il_7] > 7 && il_7 > 0) il_7--;
1778 Output::echo(5,
"Steady state normalized by fL(L[%d] = %f) = %f.\n", il_7, L_1d[il_7], fL(il_7));
1785 for (
int i = 0; i < fL.
size_x; i++)
if (fL[i] < min_f) fL[i] = min_f;
1798 for (im = 0; im < pc.
size; im++) {
1799 for (il = 0; il < L.
size; il++) {
1800 for (ia = 0; ia < alpha.
size; ia++) {
1821 arr[il][im][ia] = arr[il][im][ia] * pow(sinh(alpha.
arr[il][im][ia]/
VF::alc(L.
arr[il][im][ia])),2)/pow(sinh(1.0),2);
1837 double Lmax, Lmin, Lgrid, Lpp;
1838 int Lsize_7, addGrid;
1839 double tau_p, tau_c, alpha1, p1;
1851 int ia_90 = alpha.
size - 1;
1852 while (alpha.
arr[L.
size - 1][0][ia_90] >
VC::pi / 2 && ia_90 > 0) ia_90--;
1856 if (L.
arr[L.
size - 1][0][ia_90] < 7.0) {
1859 Lmax = L.
arr[L.
size - 1][0][ia_90];
1860 Lmin = L.
arr[0][0][ia_90];
1861 Lgrid = (Lmax - Lmin) / (L.
size - 1);
1862 Lsize_7 = (7.0 - Lmin) / Lgrid + 1;
1863 addGrid = Lsize_7 - L.
size;
1870 for (il = 0; il < L.
size; il++) {
1871 for (im = 0; im < pc.
size; im++) {
1872 L_1d_i[il] = L.
arr[il][0][ia_90];
1873 pc_i[il][im] = pc.
arr[il][im][ia_90];
1878 for (il = L.
size; il < Lsize_7; il++) {
1879 for (im = 0; im < pc.
size; im++) {
1880 L_1d_i[il] = Lmin + Lgrid * il;
1894 for (il = 0; il < L.
size; il++) {
1895 for (im = 0; im < pc.
size; im++) {
1896 L_1d_i[il] = L.
arr[il][0][ia_90];
1897 pc_i[il][im] = pc.
arr[il][im][ia_90];
1904 Lpp = (5.6 - 0.46 * Kp);
1910 for (im = 0; im < pc.
size; im++) {
1911 for (il = 0; il < Lsize_7; il++) {
1914 alpha1 = alpha.
arr[L.
size - 1][0][alpha.
size - 1];
1916 if(PSD_parameters.usetauLpp_SteadyState ==
"coulomb") {
1918 }
else if(PSD_parameters.usetauLpp_SteadyState ==
"precipitation") {
1920 }
else if (PSD_parameters.usetauLpp_SteadyState ==
"combined") {
1923 tau_0[il] = 1.0 / (1.0/tau_p+1.0/tau_c);
1938 for (il = 0; il < Lsize_7; il++) {
1939 if (fL_i[il] < min_f) {
1940 fLm[il][im] = min_f;
1942 fLm[il][im] = fL_i[il];
1948 for (im = 0; im < pc.
size; im++) {
1949 scale_factor[im] = fb_out / fLm[L.
size - 1][im];
1952 for (im = 0; im < pc.
size; im++) {
1953 for (il = 0; il < L.
size; il++) {
1954 for (ia = 0; ia < alpha.
size; ia++) {
1963 int ia_90 = alpha.
size - 1;
1964 while (alpha.
arr[il][im][ia_90] >
VC::pi / 2 && ia_90 > 0)
1968 L.
arr[il][im][ia_90]));
1969 arr[il][im][ia] = sin(alpha.
arr[il][im][ia]) * scale_factor[im]
1971 VF::Kfunc(p1), J_L7_function) * pow(1.0 / p1, 2);
1976 L.
arr[il][im][ia])) {
1977 arr[il][im][ia] = arr[il][im][ia] * pow(
1979 L.
arr[il][im][ia])), 2) / pow(sinh(1.0), 2);
1986 arr.writeToFile(
"./Debug/f_0.plt", L, pc, alpha);
2020 for (i = 0; i <= nx-1; i++) {
2022 Dxx[i] = pow(10, 0.506*Kp-9.325)*pow(L[i],10) +
VF::Dfe(alpha,Ke[i],L[i],Kp);
2027 for (i = 1; i < nx-1; i++) {
2028 alfa[i] = 1.0/((L[i+1] - L[i-1])/2)*pow(L[i],2);
2030 C_1[i] = ((Dxx[i] + Dxx[i-1])/2)/(L[i] - L[i-1])/(pow(((L[i] + L[i-1])/2), 2));
2031 C_2[i] = ((Dxx[i+1] + Dxx[i])/2)/(L[i+1] - L[i])/(pow(((L[i+1] + L[i])/2), 2));
2036 for (i = 1; i < nx-1; i++) {
2037 A[i] = alfa[i] * C_1[i];
2038 B[i] = -alfa[i] * (C_1[i] + C_2[i]) - C_3[i];
2039 C[i] = alfa[i] * C_2[i];
2042 for (i = 0; i < nx; i++) D[i] = 0;;
2064 tridag(&A[0], &B[0], &C[0], &D[0], &f[0], nx);
2081 for (il = 0; il < L.
size; il++) L_1d[il] = L.
arr[il][0][alpha.
size-1];
2083 for (
int i = 0; i < fL.
size_x; i++)
if (fL[i] < min_f) fL[i] = min_f;
2090 for (im = 0; im < pc.
size; im++) {
2091 for (il = 0; il < L.
size; il++) {
2092 for (ia = 0; ia < alpha.
size; ia++) {
2094 arr[il][im][ia] = sin(alpha.
arr[il][im][ia]) * fL[il] * L_UpperBoundaryCondition[im][ia] * pow(1.0/p1,2);
2100 arr.writeToFile(
"./Debug/f_0.plt", L, pc, alpha);
2131 ifstream input_f(filename);
2132 if (input_f != NULL && !input_f.eof()) {
2133 for (k = 0; k <= lmax; k++) {
2134 input_f >> orbit[k] >> dayno[k] >> eventtime[k];
2135 for (j = 0; j <= jmax; j++) {
2136 input_f >> i >> energy[j] >> jperp[k][j];
2137 energy[j] = energy[j]/1000;
2141 throw error_msg(
"INITIAL_PSD_ERROR",
"Error reading initial conditions input file %s.\n", filename);
2149 for (i = 0; i <= jmax; i++) jcopy[i] = jperp[2][i];
2163 for (i = 0; i < pc.
size; i++)
2167 for (i = 0; i < jmax+2; i++) {
2168 jcopy[i] = log10(jcopy[i]);
2170 j_0m.Spline(jcopy, energy, epc1D, jcopy[0], log10(
VC::zero_f), 0, 1);
2171 for (i = 0; i < pc.
size; i++) {
2172 j_0m[i] = pow(10, j_0m[i]);
2176 for (i = 0; i < pc.
size; i++) {
2178 f_0m[i] = j_0m[i]/pow(pc.
arr[0][i][0], 2);
2182 j_0m.writeToFile(
"./Debug/j_0m.plt", epc1D);
2186 for (il = 0; il < L.
size; il++) {
2187 for (im = 0; im < pc.
size; im++) {
2188 for (ia = 0; ia < alpha.
size; ia++) {
2189 arr[il][im][ia] = f_0m[im] * sin( alpha.
arr[il][im][ia] );
2191 arr[il][im][ia] = arr[il][im][ia] * pow(sinh(alpha.
arr[il][im][ia]/
VF::alc(L.
arr[il][im][ia])),2)/pow(sinh(1.0),2);
2211 double Ke, alf, flux;
2213 for (il = 0; il < L.
size; il++) {
2214 for (im = 0; im < pc.
size; im++) {
2215 for (ia = 0; ia < alpha.
size; ia++) {
2218 alf = alpha.
arr[il][im][ia];
2219 flux =
exp(-(Ke-0.2)/0.1);
2221 arr[il][im][ia] = flux / pow(pc.
arr[0][im][0], 2);
2257 for (i = 0; i <= nx-1; i++) {
2259 Dxx[i] = pow(10, 0.506*Kp-9.325)*pow(L[i],10);
2262 for (i = 1; i < nx-1; i++) {
2263 alfa[i] = 1.0/((L[i+1] - L[i-1])/2)*pow(L[i],2);
2265 C_1[i] = ((Dxx[i] + Dxx[i-1])/2)/(L[i] - L[i-1])/(pow(((L[i] + L[i-1])/2), 2));
2266 C_2[i] = ((Dxx[i+1] + Dxx[i])/2)/(L[i+1] - L[i])/(pow(((L[i+1] + L[i])/2), 2));
2271 for (i = 1; i < nx-1; i++) {
2272 A[i] = alfa[i] * C_1[i];
2273 B[i] = -alfa[i] * (C_1[i] + C_2[i]) - C_3[i];
2274 C[i] = alfa[i] * C_2[i];
2277 for (i = 0; i < nx; i++) D[i] = 0;;
2299 tridag(&A[0], &B[0], &C[0], &D[0], &f[0], nx);