VERB_code_2.2  2
 All Classes Namespaces Files Functions Variables Typedefs Macros Pages
DiffusionCoefficient.cpp
Go to the documentation of this file.
1 /**
2  * Diffusion coefficients calculation, loading, activating, scaling etc code
3  *
4  * \file DiffusionCoefficient.cpp
5  *
6  * \author Developed by Yuri Shprits
7  */
8 
9 #include "DiffusionCoefficient.h"
10 
11 #include <math.h>
12 #include <vector>
13 #include "rroots.h"
14 #include "../VariousFunctions/bisection.h"
15 
16 #include "../Logging/Output.h"
17 #include "../Exceptions/error.h"
18 
19 #include <ctime>
20 
21 using namespace std;
22 
23 #define double_zero 1.e-21
24 #define min_Dxx 1.e-21
25 
26 #if defined(_WIN32) || defined(_WIN64)
27 #define strncasecmp _strnicmp
28 #define strcasecmp _stricmp
29 #endif
30 
31 
32 /** Function Get - retutn diffusion coefficient by loading or calculating it
33  * Two way of call function:
34  * With or without grid parameter
35  * If it called without grid parameter, it takes grid stored in the class
36  * \todo Remove grid from the class DiffusionCoefficients and remove functions, these assume it there
37  * \param DxxParameters - diffusion coefficients param structure
38  */
40  bool Load_Failed = false;
41 
42  if (!arr.initialized) {
43  // Allocating memory for the parent matrix class (the function is from Matrix class)
44  arr.AllocateMemory(grid.L.size, grid.pc.size, grid.alpha.size);
45  } else {
46  arr = 0;
47  }
48 
49  arr.name = DxxParameters.DxxName + " " + DxxParameters.waveName + "";
50  this->DxxParameters = DxxParameters;
51 
52  this->type = DxxParameters.DxxType;
53  this->useScale = DxxParameters.useScale;
54 
55  if (DxxParameters.loadOrCalculate != "LOC_LOAD" && DxxParameters.loadOrCalculate != "LOC_LOADORCALCULATE" && DxxParameters.loadOrCalculate != "LOC_CALCULATE") {
56  throw error_msg("DXX_GET", "Undefined load or calculate choice.");
57  }
58 
59  // loading or calculating diffusion coefficients
60  // if DxxParameters.loadOrCalculate == LOC_LOADORCALCULATE - try to load first, and calculate in case of error
61  if (DxxParameters.loadOrCalculate == "LOC_LOAD" || DxxParameters.loadOrCalculate == "LOC_LOADORCALCULATE") {
62  try {
63  // try to load
64  LoadDiffusionCoefficient(grid.L, grid.epc, grid.alpha, DxxParameters.filename.c_str(), DxxParameters.filetype);
65  } catch (error_msg err) {
66  if (DxxParameters.loadOrCalculate == "LOC_LOADORCALCULATE") {
67  Output::echo(2, "%s", err.what().c_str());
68  Load_Failed = true;
69  } else {
70  err.add("UNKNOWN_ERROR", "Unknown error while loading diffusion coefficient file");
71  throw err;
72  }
73  }
74  }
75 
76  // Calculating if need to calculate or failed to load
77  if (DxxParameters.loadOrCalculate == "LOC_CALCULATE" || (DxxParameters.loadOrCalculate == "LOC_LOADORCALCULATE" && Load_Failed)) {
78  Output::echo(3, "Calculating %s...\n", arr.name.c_str());
79  Calculate(grid.L, grid.epc, grid.alpha, DxxParameters);
80  // Transformation to days^-1
81  arr = arr * (60.0 * 60 * 24);
82  // plotting the file, that we could not open to use it next time
83  arr.writeToFile(DxxParameters.filename.c_str(), grid.L.arr, grid.epc.arr, grid.alpha.arr);
84  }
85 
86  // MLT averaging
87  if (DxxParameters.MLT_averaging != 1) {
88  arr = arr * DxxParameters.MLT_averaging;
89  Output::echo(3, "%s is MLT-averaged as %f %%.\n", arr.name.c_str(), DxxParameters.MLT_averaging*100);
90  }
91  // Multiplicator
92  if (DxxParameters.multiplicator != 1) {
93  arr = arr * DxxParameters.multiplicator;
94  Output::echo(3, "%s is multiplied by %f.\n", arr.name.c_str(), DxxParameters.multiplicator);
95  }
96  is_active = false;
97 
98  // Normalizing diffusion coefficients to 'pc':
99  // Version 1 (slowest):
100  //if (this->type == "DCT_Dpp" || this->type == "DCT_DppLpp") arr = arr * grid.pc.arr * grid.pc.arr;
101  //if (this->type == "DCT_Dpa" || this->type == "DCT_DpaLpp") arr = arr * grid.pc.arr;
102  // Version 2 (faster):
103  if (this->type == "DCT_Dpp" || this->type == "DCT_DppLpp") arr.times_equal(grid.pc.arr.times(grid.pc.arr)); // Dxx *= (pc * pc)
104  if (this->type == "DCT_Dpa" || this->type == "DCT_DpaLpp") arr.times_equal(grid.pc.arr); // Dxx *= pc
105 
106  arr.change_ind = clock();
107 }
108 
109 
110 /** Calculating the diffusion coefficients
111  * \param &L - Grid element L
112  * \param &epc - Grid element epc
113  * \param &Alpha - Grid element Alpha
114  * \param DxxParameters - Dxx parameters structure
115  */
117  int il, im, ia;
118  bool positivonly=true;
119  double (*int_Dxx_loc)(double lambda, DxxParameters_structure DxxParameters);
120 
121  // Chousing the function by the type
122  if (DxxParameters.DxxType == "DCT_Daa" || DxxParameters.DxxType == "DCT_DaaLpp") {
123  int_Dxx_loc = int_Daa_loc;
124  } else if (DxxParameters.DxxType == "DCT_Dpp" || DxxParameters.DxxType == "DCT_DppLpp") {
125  int_Dxx_loc = int_Dpp_loc;
126  } else if (DxxParameters.DxxType == "DCT_Dpa" || DxxParameters.DxxType == "DCT_DpaLpp") {
127  int_Dxx_loc = int_Dpa_loc;
128  positivonly = false;
129  } else {
130  throw error_msg("DXX_TYPE", "Unknown diffusion coefficient");
131  }
132 
133  Output::echo(0, "Calculating %s for %s... ", DxxParameters.DxxName.c_str(), DxxParameters.waveName.c_str());
134 
135  for (il = 0; il < L.size ; il++) {
136  for (im = 0; im < epc.size; im++) {
137  for (ia = 0; ia < Alpha.size; ia++) {
138  // caltulating the value in one grid point
139  // /2 is because missing factor of 2 in Sammers formula
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;
141  // remove zeros and/or negative values
142  if (positivonly && (arr[il][im][ia] < min_Dxx)) {
143  arr[il][im][ia] = min_Dxx;
144  } else if (!(arr[il][im][ia] >= 0) && !(arr[il][im][ia] <= 0)) { // if NAN
145  if (positivonly) {
146  arr[il][im][ia] = min_Dxx;
147  } else {
148  arr[il][im][ia] = 0;
149  }
150  }
151  }
152  }
153  }
154 
155  Output::echo(2, "done!\n");
156 }
157 
158 // LOADING PROCEDURES
159 /** Loading coefficient from file
160  * runs loading function according file format (filetype)
161  */
162 bool DiffusionCoefficient::LoadDiffusionCoefficient(GridElement &L, GridElement &pc, GridElement &alpha, string D_filename, string filetype) {
163  if (!arr.initialized) {
164  throw error_msg("DXX_UNINITIALIZED_MATRIX", "Matrix has not been initialized.");
165  }
166  Output::echo(2, "Loading %s... ", D_filename.c_str());
167  // using different loading functions according to type of file
168  // like, file with/without grid, different order of grid, etc
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") {
176  // if it is "simple" (just matrix_array) - using corresponding function
177  LoadDiffusionCoefficientFromPlaneFile(L, pc, alpha, D_filename);
178  } else {
179  throw error_msg("DXX_LOAD_TYPE_ERR", "Unknown coefficients file format %s", filetype.c_str());
180  }
181  Output::echo(2, "done!\n");
182  return true;
183 }
184 
185 /** Loading coefficient from plane file
186  */
188  // iterators
189  int il, im, ia;
190 
191  // making stresm
192  ifstream input_D(D_filename.c_str());
193  // checking if file was opened
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++) {
198  // load file into D
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());
202  }
203  }
204  } else {
205  throw error_msg("DXX_LOAD_GRID_ERR", "error reading input file %s\n", D_filename.c_str());
206  //getchar();
207  return false;
208  }
209  input_D.close();
210  return true;
211 
212 }
213 
214 
215 /** Loading coef from file with grid
216  */
217 bool DiffusionCoefficient::LoadDiffusionCoefficientFromFileWithGrid(GridElement &L, GridElement &epc, GridElement &alpha, string D_filename, string gridOrder) {
218  int il, im, ia;
219  double Load_L, Load_epc, Load_alpha;// temporaraly variables for the grid point
220  double err = 5.e-3;
221  string tmp;
222 
223  // making stresm
224  ifstream input_D(D_filename.c_str());
225  // checking if file was opened
226  if (input_D == NULL) {
227  throw error_msg("DXX_FILE_OPEN_ERR", "Diffusion coefficient file %s not found.", D_filename.c_str());
228  }
229 
230  // Skip header!
231  input_D >> tmp;
232  while (strcasecmp(tmp.c_str(), "ZONE") != 0 && !input_D.eof() ) {
233  input_D >> tmp;
234  Output::echo(2, "%s ", tmp.c_str());
235  if (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());
237  }
238  }
239  // read to the end of the line with 'zone'
240  input_D.ignore(9999, '\n');
241 
242  if (!input_D.eof()) {
243  // load grid
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;
249  // load D
250  //input_D >> (*this)[il][im][ia];
251  input_D >> arr[il][im][ia];
252  // check if grid is the same
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]);
255  }
256  }
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 ;
262  // load D
263  //input_D >> (*this)[il][im][ia];
264  input_D >> arr[il][im][ia];
265  // check if grid is the same
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]);
268  }
269  }
270  } else {
271  throw error_msg("UNKNOWN_GRID_TYPE", "Unknown grid order: %s", gridOrder.c_str());
272  }
273  }
274  input_D.close();
275  return true;
276 
277 }
278 
279 /** Making DLL
280  */
281 void DiffusionCoefficient::MakeDLL(GridElement &L, GridElement &pc, GridElement &alpha, double Kp, string DLLType) {
282  if (DLLType == "DLLT_BE") {
283  MakeDLL_BE(L, pc, alpha, Kp);
284  }
285  else if (DLLType == "DLLT_B") {
286  MakeDLL_B(L, pc, alpha, Kp);
287  }
288  else if (DLLType == "DLLT_FAKE") {
289  MakeDLL_FAKE(L, pc, alpha, Kp);
290  }
291  else if (DLLType == "DLLT_Ozeke") {
292  MakeDLL_Ozeke(L, pc, alpha, Kp);
293  }
294  else if (DLLType == "DLLT_Ozeke_ME") {
295  MakeDLL_Ozeke_ME(L, pc, alpha, Kp);
296  }
297  else {
298  throw error_msg("DLL_ERROR", "Unknown DLL type");
299  }
300  arr.change_ind = clock();
301 
302 }
303 
304 /** Making DLL
305  */
307  int il, im, ia;
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);
312  //(*this)[il][im][ia] = 0.0;
313  }
314  }
315  }
316 }
317 
318 /** Making DLL
319  */
321  int il, im, ia;
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);;
326  }
327  }
328  }
329 }
330 
331 /** Other version of DLL making procedure
332  */
334  int il, im, ia;
335 
336  for (il = 0; il < L.size; il++) {
337  for (im = 0; im < pc.size; im++) {
338  for (ia = 0; ia < alpha.size; ia++) {
339 
340  // electrostatic radial diffusion coefficient in 1/day
341  arr[il][im][ia] = VF::Dfe(alpha.arr[il][im][ia],VF::Kfunc(pc.arr[il][im][ia]), L.arr[il][im][ia],Kp);
342  //(*this)[il][im][ia] = 0.2*VF::Dfe(alpha[il][im][ia],VF::Kfunc(pc[il][im][ia]),L[il][im][ia],Kp);
343 
344  // Add magnetostatic radial diffusion coefficient in 1/day
345  arr[il][im][ia] = arr[il][im][ia] + Df(L.arr[il][im][ia], Kp);
346 
347  }
348  }
349  }
350 }
351 
352 
353 /** Other version of DLL making procedure
354  */
356  int il, im, ia;
357 
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);
362  }
363  }
364  }
365 
366 }
367 
368 /** Making DLL Ozeke
369  * WARNING! This code uses prepublication values! WARNING!
370  * Previous version of DLL:
371  * Ozeke, L. G., I. R. Mann, K. R. Murphy, I. J. Rae, D. K. Milling, S. R. Elkington, A. A. A. Chan, and H. J. Singer (2012),
372  * ULF Wave Derived Radiation Belt Radial Diffusion Coefficients,
373  * J. Geophys. Res., doi:10.1029/2011JA017463, in press. (accepted 23 February 2012)
374  * Magnetic component
375  */
377  int il, im, ia;
378  for (il = 0; il < L.size; il++) {
379  for (im = 0; im < pc.size; im++) {
380  for (ia = 0; ia < alpha.size; ia++) {
381  arr[il][im][ia] = Df_Ozeke(L.arr[il][im][ia], Kp); // this.arr
382  }
383  }
384  }
385 }
386 
387 /** Making DLL Ozeke
388  * WARNING! This code uses prepublication values! WARNING!
389  * Previous version of DLL:
390  * Ozeke, L. G., I. R. Mann, K. R. Murphy, I. J. Rae, D. K. Milling, S. R. Elkington, A. A. A. Chan, and H. J. Singer (2012),
391  * ULF Wave Derived Radiation Belt Radial Diffusion Coefficients,
392  * J. Geophys. Res., doi:10.1029/2011JA017463, in press. (accepted 23 February 2012)
393  * Magnetic and Electric component
394  */
396  int il, im, 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++) {
400  arr[il][im][ia] = Df_Ozeke(L.arr[il][im][ia], Kp) + Df_Ozeke_E(L.arr[il][im][ia], Kp); // this.arr
401  }
402  }
403  }
404 }
405 
406 
407 /** Activation (enabling/disabling) and scaling of the diffusion coefficients in one group.
408  * It loops thrue the list of the diffusion coefficients, enable which needs to be enabled, disable, which needs to be disabled and scale all enabled
409  */
410 bool DiffusionCoefficientsGroup::ActivateAndScale(double time, double Kp) {
411  // The last version of activated stuff.
412  // Activating and scaling diffusion coefficients for different wave types
413  double scaleCoeff = 0;
414  unsigned int Dxx_it;
415  // \todo check, if we need to make any changes, otherwise - just exit
416  // check if we really need to do any changes
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) {
420  // if inactive - mark as active
421  if (!DxxList[Dxx_it].is_active) {
422  need_changes = true;
423  break;
424  }
425  // if this Dxx needs to be scaled - scale it
426  if (DxxList[Dxx_it].useScale) {
427  need_changes = true;
428  break;
429  }
430  } else {
431  // if active - mark as inactive
432  if (DxxList[Dxx_it].is_active) {
433  need_changes = true;
434  break;
435  }
436  }
437  }
438  if (need_changes == false) return false;
439 
440  CurrentDxx.arr = 0;
441  // Indicator, that we changed the coefficient (and now we probably need to recalculate the model matrix)
442  CurrentDxx.arr.change_ind = clock();
443  // loop thrue list of the diffusion coefficients
444  for (Dxx_it = 0; Dxx_it < DxxList.size(); Dxx_it++) {
445  // if it is time to work - activate
446  if (time >= DxxList[Dxx_it].DxxParameters.time_start && time < DxxList[Dxx_it].DxxParameters.time_end) {
447  // if inactive - mark as active
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());
451  }
452  // if this Dxx needs to be scaled - scale it
453  if (DxxList[Dxx_it].useScale) {
454  // getting scale coefficient
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);
457  // scalling and adding value to the actualized Dxx which is used in calculations
458  CurrentDxx.arr = CurrentDxx.arr + DxxList[Dxx_it].arr * scaleCoeff;
459  } else {
460  CurrentDxx.arr = CurrentDxx.arr + DxxList[Dxx_it].arr;
461  }
462  } else {
463  // if active - mark as inactive
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());
467  }
468  }
469  }
470  return false;
471 }
472 
473 /** Scaling procedure
474  * \todo Upgrade scaling diffusion coefficients procedure with scaling according to external array
475  */
476 double DiffusionCoefficient::Scale(double Kp) {
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") {
480  // scaling for chorus waves
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);// * 1.e-12 * 1.e8;
484  if (Kp > 2.333) Bw2 = 2.0*pow(10, 2.5 + 0.18 * Kp);// * 1.e-12 * 1.e8;
485  //////////////////////////////////////////////////
486  /* Ughm? Commented by Dmitriy
487  if(DxxParameters.waveType == "WT_CHORUS_DAY") {
488  return (Bw2/Bw2_0)*2.5;
489  } else {
490  return (Bw2/Bw2_0)/2.5;
491  }*/
492  ///////////////////////////////////////////////////
493  return (Bw2/Bw2_0);
494  } else if (DxxParameters.waveType == "WT_HISS_PP" || DxxParameters.waveType == "WT_HISS_PLUME") {
495  // scale hiss diffusion coefficients is a square of Bw scale
496  // return pow(Kp/DxxParameters.DxxKp, 2);
497  return pow(Kp/DxxParameters.DxxKp, 2);
498  } else if (DxxParameters.waveType == "WT_EMIC_PLUME" ) {
499  // scale EMIC diffusion coefficients as a Bw to the forth power
500  return pow(Kp/DxxParameters.DxxKp, 4);
501 //////////////////////////kckim////////////////////////////////////
502  } else if (DxxParameters.waveType == "WT_LIGHTNING") {
503  return pow(Kp/DxxParameters.DxxKp, 2);
504  //
505 ///////////////////////////////////////////////////////////////////////////
506 
507  } else {
508  //throw error_msg("DXX_SCALING", "No scaling for wave type %s", DxxParameters.waveType.c_str());
509  return 1;
510  }
511 }
512 
513 
514 double Alpha_ne(double pangle, double lambda, double L) {
515  return asin( sqrt( f1(lambda)) * sin(pangle) );
516 }
517 
518 double f1(double lambda) {
519  return sqrt(4.0 - 3.0 * pow(cos(lambda), 2))/(pow(cos(lambda), 6));
520 }
521 
522 double B(double lambda, double L) {
523  return f1(lambda) / pow(L, 3) * VC::B_0;
524 }
525 
526 double func_tmp (double x, double Alpha) {
527  return pow(x, 6) + (3.0 * pow(sin(Alpha), 4)) * x - 4.0 * pow(sin(Alpha), 4);
528 }
529 
530 double F_cap(double x, double y, double b, double s, double epsilon, DxxParameters_structure DxxParameters) {
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);
537 }
538 
539 double F_cap2(double x, double y, double a, double beta, double mu, double s, double epsilon, double Alpha_star, DxxParameters_structure DxxParameters) {
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)));
543  return up/down;
544 }
545 
546 double quad1(double (*func)(double lambda, DxxParameters_structure DxxParameters), double a, double b, int M, DxxParameters_structure DxxParameters) {
547  // f -integrated function
548  // a and b upper and lower limits
549  // M number of subintegrals
550  double h = (b - a)/M;
551  double s = func(a, DxxParameters);
552  double x;
553  int k;
554  for (k = 1; k <= (M-1); k++) {
555  x = a + h*k;
556  s = s + func(x, DxxParameters);
557  }
558  s = h * (func(a, DxxParameters) + func(b, DxxParameters))/2.0 + h*s;
559  return s;
560 }
561 
562 double Dxx_ba(double L1, double EMeV, double Alpha, double int_Dxx_loc(double lambda, DxxParameters_structure DxxParameters), DxxParameters_structure DxxParameters) {
563  // Searching for zero of function func_tmp=x.^6+(3*sin(DxxParameters.Alpha_eq_tmp)^4.)*x-4*sin(DxxParameters.Alpha_eq_tmp)^4.;
564  // near the point 0.5
565  DxxParameters.L = L1;
566  DxxParameters.EMeV = EMeV;
567  DxxParameters.Alpha = Alpha;
568 
569  double lambda_m;
570  //lambda_m = acos(sqrt(fzero(func_tmp, 0.5)))*0.99;
571  lambda_m = acos(sqrt(bisection (func_tmp, DxxParameters.Alpha, 0, 1)))*DxxParameters.mirror_point_coeff;
572  if (lambda_m > DxxParameters.lam_max*VC::pi/180) lambda_m = DxxParameters.lam_max*VC::pi/180.;
573  // quad1 - calculating integral
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));
575 }
576 
577 
578 double int_Daa_loc(double lambda, DxxParameters_structure DxxParameters) {
579  //double Daa = Daa_local(lambda);
580  double Daa = Dxx_local(lambda, Daa_root, 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));
583 }
584 
585 double int_Dpp_loc(double lambda, DxxParameters_structure DxxParameters) {
586  //double Dpp = Dpp_local(lambda);
587  double Dpp = Dxx_local(lambda, Dpp_root, DxxParameters);
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));
590 }
591 
592 double int_Dpa_loc(double lambda, DxxParameters_structure DxxParameters) {
593  //double Dpa = Dpa_local(lambda);
594  double Dpa = Dxx_local(lambda, Dpa_root, DxxParameters);
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));
597 }
598 
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 ) {
600 
601  double Omega_e = VC::qe * B(lambda, DxxParameters.L) / VC::m / VC::c_cgs;
602  double Omega_e_eq = VC::qe * B(0, DxxParameters.L) / VC::m / VC::c_cgs;
603  //double N_0 = 124.0 * pow((3.0 / L), 4);
604  //double Omega_p = sqrt(4.0 * pi * N_0 * pow(qe, 2) / m);
605  // f = Omega_p/Omega_e_eq;
606  double Omega_p;
607  //if (DxxParameters.f > 0) {
608  double N_0;
609  double LT;
610  if (DxxParameters.numberDensity == "ND_CONSTANT") {
611  Omega_p = DxxParameters.f * Omega_e_eq;
612  } else if (DxxParameters.numberDensity == "ND_CHORUS") {
613  N_0 = 124.0*(pow((3.0/DxxParameters.L),4));
614  Omega_p = sqrt(4.*VC::pi*N_0*pow(VC::qe,2)/VC::m);
615  DxxParameters.f = Omega_p / Omega_e_eq;
616 // Output::echo(2, "%e \n", DxxParameters.f);
617  } else if (DxxParameters.numberDensity == "ND_CHORUS_NIGHT") {
618  LT = 3;
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);
620  Omega_p = sqrt(4.*VC::pi*N_0*pow(VC::qe,2)/VC::m);
621  DxxParameters.f = Omega_p / Omega_e_eq;
622 // if (DxxParameters.L >= 4.5) { Output::echo(2, "%e \n", DxxParameters.f); }
623  } else if (DxxParameters.numberDensity == "ND_CHORUS_DAY") {
624  LT = 9;
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);
626  Omega_p = sqrt(4.*VC::pi*N_0*pow(VC::qe,2)/VC::m);
627  DxxParameters.f = Omega_p / Omega_e_eq;
628 // if (DxxParameters.L >= 4.5) { Output::echo(2, "%e \n", DxxParameters.f); }
629  } else if (DxxParameters.numberDensity == "ND_HISS") {
630  N_0 = pow(10, (-0.3145*DxxParameters.L + 3.9043));
631  Omega_p = sqrt(4.*VC::pi*N_0*pow(VC::qe,2)/VC::m);
632  DxxParameters.f = Omega_p / Omega_e_eq;
633 // Output::echo(2, "%e \n", DxxParameters.f);
634  } else if (DxxParameters.numberDensity == "ND_HISS2") {
635  N_0 = 1e3*pow(4.0/DxxParameters.L, 4);
636  Omega_p = sqrt(4.*VC::pi*N_0*pow(VC::qe,2)/VC::m);
637  DxxParameters.f = Omega_p / Omega_e_eq;
638 // Output::echo(2, "%e \n", DxxParameters.f);
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));
642  Omega_p = sqrt(4.*VC::pi*N_0*pow(VC::qe,2)/VC::m);
643  DxxParameters.f = Omega_p / Omega_e_eq;
644  }else{
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);
649  Omega_p = sqrt(4.*VC::pi*N_0*pow(VC::qe,2)/VC::m);
650  DxxParameters.f = Omega_p / Omega_e_eq;
651  }
652 
653  } else {
654  throw error_msg("DXX_DENSITY_ERROR", "Unknown number density type: %s\n", DxxParameters.numberDensity.c_str());
655  }
656 
657 
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));
662 
663  double lam_sp = -1; // -1 corresponds to electrons
664  double gamma = DxxParameters.EMeV + 1;
665  double a = DxxParameters.s * lam_sp / gamma;
666 
667  double Alpha_star = pow((Omega_e / Omega_p), 2);
668  double b = (1 + epsilon)/Alpha_star;
669 
670  //double d_omega = d_omega_per * Omega_e_eq;
671  //double Omega_m = Omega_m_per * Omega_e_eq;
672  double d_omega;
673  double Omega_m;
674  double Omega_1;
675  double Omega_2;
676 
677  if (DxxParameters.Omega_mType == "OT_PARTIAL") {
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") {
683  Omega_m = DxxParameters.Omega_m * 2 * VC::pi;
684  d_omega = DxxParameters.d_omega * 2 * VC::pi;
685  Omega_1 = DxxParameters.omega_lc * 2 * VC::pi;
686  Omega_2 = DxxParameters.omega_uc * 2 * VC::pi;
687  } else {
688  throw error_msg("DXX_OMEGA_M_TYPE", "Unknown Omega_m type: ", DxxParameters.Omega_mType.c_str());
689  }
690 
691  // if (DxxParameters.Omega_mType == "partial") {
692  // } else if (DxxParameters.Omega_mType == "constant") {
693  // } else {
694  // }
695 
696  //double R = pow((DxxParameters.dB / B(lambda, DxxParameters.L)), 2);
697  double Bw;
698  if (DxxParameters.BwFromLambda) {
699  Bw = pow(10, 0.75 + 0.04*(lambda*180/VC::pi))*1.e-8 * DxxParameters.Bw; // should be gauss...
700  } else {
701  Bw = DxxParameters.Bw;
702  }
703  double R = pow((Bw / B(lambda, DxxParameters.L)), 2);
704 
705  double x_m=0;
706  double d_x=0;
707  double x_1=0;
708  double x_2=0;
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") {
715  double cfm = 16;
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;
720  }
721 
722  std::vector<double> xx = rrouts(x_1, x_2, DxxParameters.eta1, DxxParameters.eta2, DxxParameters.eta3, epsilon, beta, mu, Alpha_star, a, DxxParameters);
723  //
724  double x, y;
725  double Dxx = 0, Dxx_tmp = 0;
726  for (unsigned int i = 1; i <= xx.size();i++) {
727  x = xx[i-1];
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);
730  Dxx = Dxx_tmp + Dxx;
731  }
732  if (xx.size() == 0) {
733  Dxx = 0;
734  }
735 
736  return Dxx;
737 }
738 
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)));
744 }
745 
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)));
751 }
752 
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))) *
755  //((R * (pow(x/y, 2)) * abs(F_cap(x,y,b,s,epsilon)))/
756  //(d_x*abs(beta *mu-F_cap(x,y,b,s,epsilon))))*
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)));
760 }
761 
762 
763 //std::vector<double> rrouts(double Omega_e_eq, double Omega_e, double eta1, double eta2, double eta3, double epsilon, double beta, double mu, double Alpha_star, double a, DxxParameters_structure DxxParameters) {
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) {
765 
766  //double Omega_1 = DxxParameters.omega_lc * Omega_e_eq * epsilon / DxxParameters.cfm;
767  //double Omega_2 = DxxParameters.omega_uc * Omega_e_eq * epsilon / DxxParameters.cfm;
768  //double x_1 = Omega_1 / Omega_e;
769  //double x_2 = Omega_2 / Omega_e;
770 
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;
775 
776  double d0=1;
777  double d1=eta1*epsilon;
778  double d2=eta2/4.*epsilon;
779  double d3=eta3/16.*epsilon;
780 
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;
785 
786  double h4=1.;
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;
791 
792  double A[7];
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;
799  A[6] = -pow(a,2)*h0;
800 
801  //a6=(beta^2*mu^2-1.)*h4;
802  //a5=(beta^2*mu^2-1.)*h3-2*a*h4;
803  //a4=(beta^2*mu^2-1.)*h2-2*a*h3-beta^2*mu^2*1./Alpha_star*g3-a^2*h4;
804  //a3=(beta^2*mu^2-1.)*h1-2*a*h2-beta^2*mu^2*1./Alpha_star*g2-a^2*h3;
805  //a2=(beta^2*mu^2-1.)*h0-2*a*h1-beta^2*mu^2*1./Alpha_star*g1-a^2*h2;
806  //a1=-2*a*h0-beta^2*mu^2*1./Alpha_star*g0-a^2*h1;
807  //a0=-a^2*h0;
808 
809  std::vector<double> xx;
810  double r[7], wr[6],wi[6],quad[2];
811  int n, numr;
812  n=6;
813  // initialize estimate for 1st root pair
814  quad[0] = 2.71828e-1;
815  quad[1] = 3.14159e-1;
816  // get roots
817  get_quads(A, n, quad, r);
818  numr = roots(r, n, wr, wi);
819 
820  for (int i = 0; i < numr; i++) {
821  if (fabs(wi[i]) < double_zero && wr[i] < x_2 && wr[i] > x_1) {
822  xx.push_back(wr[i]);
823  }
824  }
825 
826  return xx;
827 
828 }
829 
830 
831 /// Radial Diffusion coeficeint computed following [Brautigam and Albet , 2000]
832 double Df (double L, double Kp) {
833  return pow(10, 0.506*Kp-9.325)*pow(L,10);
834 }
835 
836 /// Radial Diffusion coeficeint computed Ozeke,
837 // Ozeke et al, 2013
838 // Analytic Expressions for ULL Wave Radiaion Belt Radial Diffusion Coefficient // (before publishing)
839 /// DLL_B [days^-1]
840 double Df_Ozeke(double L, double Kp) {
841 
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) );
843 }
844 
845 /// DLL_E
846 double Df_Ozeke_E(double L, double Kp) {
847 
848  return 2.16e-8 * pow(L,6) * pow(10,( 0.217 * L + 0.461 * Kp ));
849 }
850 
851 
852 
853 /* DEPRECATED
854 /// Radial Diffusion coeficeint computed Ozeke, 2013 library (odc_DLL_Ozeke2013)
855 // http://irbem.svn.sourceforge.net/viewvc/irbem/extras/opendc/
856 double Df_Ozeke(double L, double Kp) {
857 
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 );
859 }
860 
861 double Df_Ozeke_E(double L, double Kp) {
862 
863  return 5.75e-9 * pow(L,6) * pow(10,( 0.208 * L + 0.457 * Kp ));
864 }
865 */
866 
867 
871  Parameters_structure &parameters, Grid &radialDiffusionGrid, Grid &localDiffusionsGrid) {
872 
873  // Diffusion coefficient group (Lists of similar diffusion coefficients, like Daa or Dpcpc...)
874  DLL.arr.AllocateMemory(radialDiffusionGrid.L.size, radialDiffusionGrid.pc.size, radialDiffusionGrid.alpha.size);
875  Daa.CurrentDxx.arr.AllocateMemory(radialDiffusionGrid.L.size, radialDiffusionGrid.pc.size, radialDiffusionGrid.alpha.size);
876  DaaLpp.CurrentDxx.arr.AllocateMemory(radialDiffusionGrid.L.size, radialDiffusionGrid.pc.size, radialDiffusionGrid.alpha.size);
877  Dpcpc.CurrentDxx.arr.AllocateMemory(radialDiffusionGrid.L.size, radialDiffusionGrid.pc.size, radialDiffusionGrid.alpha.size);
878  DpcpcLpp.CurrentDxx.arr.AllocateMemory(radialDiffusionGrid.L.size, radialDiffusionGrid.pc.size, radialDiffusionGrid.alpha.size);
879  Dpca.CurrentDxx.arr.AllocateMemory(radialDiffusionGrid.L.size, radialDiffusionGrid.pc.size, radialDiffusionGrid.alpha.size);
880  DpcaLpp.CurrentDxx.arr.AllocateMemory(radialDiffusionGrid.L.size, radialDiffusionGrid.pc.size, radialDiffusionGrid.alpha.size);
881 
882  // Make it zero
883  DLL.arr = 0;
884  Daa.CurrentDxx.arr = 0;
885  DaaLpp.CurrentDxx.arr = 0;
886  Dpcpc.CurrentDxx.arr = 0;
887  DpcpcLpp.CurrentDxx.arr = 0;
888  Dpca.CurrentDxx.arr = 0;
889  DpcaLpp.CurrentDxx.arr = 0;
890 
891  // Temporary diffusion coefficient
892  DiffusionCoefficient Dxx_tmp;
893  // Create the diffusion coefficients for each of the Dxx parameters files
894  unsigned int Dxx_it;
895  for (Dxx_it = 0; Dxx_it < parameters.DxxParametersList.size(); Dxx_it++) {
896 
897  Dxx_tmp.Create_Dxx(parameters.DxxParametersList[Dxx_it], localDiffusionsGrid);
898 
899  // Add the created coefficient to an appropriate group
900  if (parameters.useAlphaDiffusion) {
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);
905  }
906 
907  if (parameters.useEnergyDiffusion) {
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);
912  }
913 
914  if (parameters.useEnergyAlphaMixedTerms) {
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);
919 
920  }
921  }
922  // Check if Daa*Dpcpc > Dpca^2 for each Dpca
923  if (parameters.useEnergyAlphaMixedTerms) {
924  int Daa_it, Dpcpc_it, Dpca_it;
925  // above Lpp
926  for (Dpca_it = 0; Dpca_it < Dpca.DxxList.size(); Dpca_it++) {
927  // For each Dpca...
928  // Search for Daa with the same wave-type
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)
931  break;
932  }
933  // Search for Dpp with the same wave-type
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)
936  break;
937  }
938  // Check if Daa*Dpcpc > Dpca^2 for each Dpca
939  int iL, im, ia;
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(),
946  iL, im, ia,
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]);
949  }
950 
951  // Below Lpp
952  for (Dpca_it = 0; Dpca_it < DpcaLpp.DxxList.size(); Dpca_it++) {
953  // For each Dpca...
954  // Search for Daa with the same wave-type
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)
957  break;
958  }
959  // Search for Dpp with the same wave-type
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)
962  break;
963  }
964  // Check if Daa*Dpcpc > Dpca^2 for each Dpca
965  int iL, im, ia;
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(),
972  iL, im, ia,
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]);
975  }
976 
977  }
978 
979  // calculating DLL If we need to do Radial diffusion
980  if (parameters.useRadialDiffusion) {
981  DLL.MakeDLL(radialDiffusionGrid.L, radialDiffusionGrid.pc, radialDiffusionGrid.alpha, parameters.Kp[0], parameters.DLLType);
982 // DLL.writeToFile("./Debug/DLL.plt", radialDiffusionGrid.L, radialDiffusionGrid.epc, radialDiffusionGrid.alpha);
983  }
984 
985 }
986 
987 
988 
989 void Output1DHeaders(ofstream &output1D,
992 
993  // Title of file (first line)
994  output1D << "Variables = \"time\", \"Kp\", \"Boundary_fluxes\", \"Lpp\", \"tau\" ";
995 
996  // Writing headers for all diffusion coefficients scaling coefficients
997  unsigned int Dxx_it;
998  for (Dxx_it = 0; Dxx_it < Daa.DxxList.size(); Dxx_it++) {
999  output1D << "\"" << Daa.DxxList[Dxx_it].arr.name << "\"" << ", ";
1000  }
1001  for (Dxx_it = 0; Dxx_it < Dpcpc.DxxList.size(); Dxx_it++) {
1002  output1D << "\"" << Dpcpc.DxxList[Dxx_it].arr.name << "\"" << ", ";
1003  }
1004  for (Dxx_it = 0; Dxx_it < Dpca.DxxList.size(); Dxx_it++) {
1005  output1D << "\"" << Dpca.DxxList[Dxx_it].arr.name << "\"" << ", ";
1006  }
1007 
1008  for (Dxx_it = 0; Dxx_it < DaaLpp.DxxList.size(); Dxx_it++) {
1009  output1D << "\"" << DaaLpp.DxxList[Dxx_it].arr.name << "\"" << ", ";
1010  }
1011  for (Dxx_it = 0; Dxx_it < DpcpcLpp.DxxList.size(); Dxx_it++) {
1012  output1D << "\"" << DpcpcLpp.DxxList[Dxx_it].arr.name << "\"" << ", ";
1013  }
1014  for (Dxx_it = 0; Dxx_it < DpcaLpp.DxxList.size(); Dxx_it++) {
1015  output1D << "\"" << DpcaLpp.DxxList[Dxx_it].arr.name << "\"" << ", ";
1016  }
1017 
1018  output1D << endl;
1019  output1D << "ZONE T=\"1d-output\" " << endl;
1020 }
1021 
1022 void Output1DValues(ofstream &output1D,
1025  double &time, Parameters_structure &parameters, int iteration) {
1026 
1027  // Title of file (first line)
1028  output1D << time << "\t" << parameters.Kp[iteration] << "\t" << parameters.Bf[iteration] << "\t" << parameters.Lpp[iteration] << "\t" << parameters.tau[iteration] << "\t"; //parameters.dBw2[iteration] << "\t" <<
1029 
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";
1033  }
1034  for (Dxx_it = 0; Dxx_it < Dpcpc.DxxList.size(); Dxx_it++) {
1035  output1D << Dpcpc.DxxList[Dxx_it].Scale(parameters.Kp[iteration]) << "\t";
1036  }
1037  for (Dxx_it = 0; Dxx_it < Dpca.DxxList.size(); Dxx_it++) {
1038  output1D << Dpca.DxxList[Dxx_it].Scale(parameters.Kp[iteration]) << "\t";
1039  }
1040 
1041  for (Dxx_it = 0; Dxx_it < DaaLpp.DxxList.size(); Dxx_it++) {
1042  output1D << DaaLpp.DxxList[Dxx_it].Scale(parameters.Kp[iteration]) << "\t";
1043  }
1044  for (Dxx_it = 0; Dxx_it < DpcpcLpp.DxxList.size(); Dxx_it++) {
1045  output1D << DpcpcLpp.DxxList[Dxx_it].Scale(parameters.Kp[iteration]) << "\t";
1046  }
1047  for (Dxx_it = 0; Dxx_it < DpcaLpp.DxxList.size(); Dxx_it++) {
1048  output1D << DpcaLpp.DxxList[Dxx_it].Scale(parameters.Kp[iteration]) << "\t";
1049  }
1050 
1051  output1D << endl;
1052 }