VERB_code_2.3
DiffusionCoefficient.cpp
Go to the documentation of this file.
1 
11 #include "DiffusionCoefficient.h"
12 
13 #include <math.h>
14 #include <vector>
15 #include "rroots.h"
16 #include "../VariousFunctions/bisection.h"
17 
18 #include "../Logging/Output.h"
19 #include "../Exceptions/error.h"
20 
21 #include <ctime>
22 
23 using namespace std;
24 
25 #define double_zero 1.e-21
26 #define min_Dxx 1.e-21
27 
28 #if defined(_WIN32) || defined(_WIN64)
29 #define strncasecmp _strnicmp
30 #define strcasecmp _stricmp
31 #endif
32 
33 
43  bool Load_Failed = false;
44 
45  if (!arr.initialized) {
46  // Allocating memory for the parent matrix class (the function is from Matrix class)
47  arr.AllocateMemory(grid.L.size, grid.pc.size, grid.alpha.size);
48  } else {
49  arr = 0;
50  }
51 
52  arr.name = DxxParameters.DxxName + " " + DxxParameters.waveName + "";
53  this->DxxParameters = DxxParameters;
54 
55  this->type = DxxParameters.DxxType;
56  this->useScale = DxxParameters.useScale;
57 
58  if (DxxParameters.loadOrCalculate != "LOC_LOAD" && DxxParameters.loadOrCalculate != "LOC_LOADORCALCULATE" && DxxParameters.loadOrCalculate != "LOC_CALCULATE") {
59  throw error_msg("DXX_GET", "Undefined load or calculate choice.");
60  }
61 
62  // loading or calculating diffusion coefficients
63  // if DxxParameters.loadOrCalculate == LOC_LOADORCALCULATE - try to load first, and calculate in case of error
64  if (DxxParameters.loadOrCalculate == "LOC_LOAD" || DxxParameters.loadOrCalculate == "LOC_LOADORCALCULATE") {
65  try {
66  // try to load
67  LoadDiffusionCoefficient(grid.L, grid.epc, grid.alpha, DxxParameters.filename.c_str(), DxxParameters.filetype);
68  } catch (error_msg err) {
69  if (DxxParameters.loadOrCalculate == "LOC_LOADORCALCULATE") {
70  Output::echo(2, "%s", err.what().c_str());
71  Load_Failed = true;
72  } else {
73  err.add("UNKNOWN_ERROR", "Unknown error while loading diffusion coefficient file");
74  throw err;
75  }
76  }
77  }
78 
79  // Calculating if need to calculate or failed to load
80  if (DxxParameters.loadOrCalculate == "LOC_CALCULATE" || (DxxParameters.loadOrCalculate == "LOC_LOADORCALCULATE" && Load_Failed)) {
81  Output::echo(3, "Calculating %s...\n", arr.name.c_str());
82  Calculate(grid.L, grid.epc, grid.alpha, DxxParameters);
83  // Transformation to days^-1
84  arr = arr * (60.0 * 60 * 24);
85  // plotting the file, that we could not open to use it next time
86  arr.writeToFile(DxxParameters.filename.c_str(), grid.L.arr, grid.epc.arr, grid.alpha.arr);
87  }
88 
89  // MLT averaging
90  if (DxxParameters.MLT_averaging != 1) {
91  arr = arr * DxxParameters.MLT_averaging;
92  Output::echo(3, "%s is MLT-averaged as %f %%.\n", arr.name.c_str(), DxxParameters.MLT_averaging*100);
93  }
94  // Multiplicator
95  if (DxxParameters.multiplicator != 1) {
96  arr = arr * DxxParameters.multiplicator;
97  Output::echo(3, "%s is multiplied by %f.\n", arr.name.c_str(), DxxParameters.multiplicator);
98  }
99  is_active = false;
100 
101  // Normalizing diffusion coefficients to 'pc':
102  // Version 1 (slowest):
103  //if (this->type == "DCT_Dpp" || this->type == "DCT_DppLpp") arr = arr * grid.pc.arr * grid.pc.arr;
104  //if (this->type == "DCT_Dpa" || this->type == "DCT_DpaLpp") arr = arr * grid.pc.arr;
105  // Version 2 (faster):
106  if (this->type == "DCT_Dpp" || this->type == "DCT_DppLpp") arr.times_equal(grid.pc.arr.times(grid.pc.arr)); // Dxx *= (pc * pc)
107  if (this->type == "DCT_Dpa" || this->type == "DCT_DpaLpp") arr.times_equal(grid.pc.arr); // Dxx *= pc
108 
109  arr.change_ind = clock();
110 }
111 
112 
120  int il, im, ia;
121  bool positivonly=true;
122  double (*int_Dxx_loc)(double lambda, DxxParameters_structure DxxParameters);
123 
124  // Chousing the function by the type
125  if (DxxParameters.DxxType == "DCT_Daa" || DxxParameters.DxxType == "DCT_DaaLpp") {
126  int_Dxx_loc = int_Daa_loc;
127  } else if (DxxParameters.DxxType == "DCT_Dpp" || DxxParameters.DxxType == "DCT_DppLpp") {
128  int_Dxx_loc = int_Dpp_loc;
129  } else if (DxxParameters.DxxType == "DCT_Dpa" || DxxParameters.DxxType == "DCT_DpaLpp") {
130  int_Dxx_loc = int_Dpa_loc;
131  positivonly = false;
132  } else {
133  throw error_msg("DXX_TYPE", "Unknown diffusion coefficient");
134  }
135 
136  Output::echo(0, "Calculating %s for %s... ", DxxParameters.DxxName.c_str(), DxxParameters.waveName.c_str());
137 
138  for (il = 0; il < L.size ; il++) {
139  for (im = 0; im < epc.size; im++) {
140  for (ia = 0; ia < Alpha.size; ia++) {
141  // caltulating the value in one grid point
142  // /2 is because missing factor of 2 in Sammers formula
143  arr[il][im][ia] = Dxx_ba(L.arr[il][im][ia], epc.arr[il][im][ia]/VC::mc2, Alpha.arr[il][im][ia], int_Dxx_loc, DxxParameters) / 2;
144  // remove zeros and/or negative values
145  if (positivonly && (arr[il][im][ia] < min_Dxx)) {
146  arr[il][im][ia] = min_Dxx;
147  } else if (!(arr[il][im][ia] >= 0) && !(arr[il][im][ia] <= 0)) { // if NAN
148  if (positivonly) {
149  arr[il][im][ia] = min_Dxx;
150  } else {
151  arr[il][im][ia] = 0;
152  }
153  }
154  }
155  }
156  }
157 
158  Output::echo(2, "done!\n");
159 }
160 
161 // LOADING PROCEDURES
165 bool DiffusionCoefficient::LoadDiffusionCoefficient(GridElement &L, GridElement &pc, GridElement &alpha, string D_filename, string filetype) {
166  if (!arr.initialized) {
167  throw error_msg("DXX_UNINITIALIZED_MATRIX", "Matrix has not been initialized.");
168  }
169  Output::echo(2, "Loading %s... ", D_filename.c_str());
170  // using different loading functions according to type of file
171  // like, file with/without grid, different order of grid, etc
172  if (filetype == "IFT_GRID") {
173  LoadDiffusionCoefficientFromFileWithGrid(L, pc, alpha, D_filename);
174  } else if (filetype == "IFT_GRID_LPA") {
175  LoadDiffusionCoefficientFromFileWithGrid(L, pc, alpha, D_filename, "IFT_GRID_LPA");
176  } else if (filetype == "IFT_GRID_APL") {
177  LoadDiffusionCoefficientFromFileWithGrid(L, pc, alpha, D_filename, "IFT_GRID_APL");
178  } else if (filetype == "IFT_PLANE") {
179  // if it is "simple" (just matrix_array) - using corresponding function
180  LoadDiffusionCoefficientFromPlaneFile(L, pc, alpha, D_filename);
181  } else {
182  throw error_msg("DXX_LOAD_TYPE_ERR", "Unknown coefficients file format %s", filetype.c_str());
183  }
184  Output::echo(2, "done!\n");
185  return true;
186 }
187 
191  // iterators
192  int il, im, ia;
193 
194  // making stresm
195  ifstream input_D(D_filename.c_str());
196  // checking if file was opened
197  if (input_D != NULL && !input_D.eof()) {
198  for (il = 0; il < L.size; il++)
199  for (im = 0; im < pc.size; im++)
200  for (ia = 0; ia < alpha.size; ia++) {
201  // load file into D
202  input_D >> arr[il][im][ia];
203  if (input_D == NULL) {
204  throw error_msg("LOAD_ERROR", "Wrong coefficient file %s (to few values)\n", D_filename.c_str());
205  }
206  }
207  } else {
208  throw error_msg("DXX_LOAD_GRID_ERR", "error reading input file %s\n", D_filename.c_str());
209  //getchar();
210  return false;
211  }
212  input_D.close();
213  return true;
214 
215 }
216 
217 
220 bool DiffusionCoefficient::LoadDiffusionCoefficientFromFileWithGrid(GridElement &L, GridElement &epc, GridElement &alpha, string D_filename, string gridOrder) {
221  int il, im, ia;
222  double Load_L, Load_epc, Load_alpha;// temporaraly variables for the grid point
223  double err = 5.e-3;
224  string tmp;
225 
226  // making stresm
227  ifstream input_D(D_filename.c_str());
228  // checking if file was opened
229  if (input_D == NULL) {
230  throw error_msg("DXX_FILE_OPEN_ERR", "Diffusion coefficient file %s not found.", D_filename.c_str());
231  }
232 
233  // Skip header!
234  input_D >> tmp;
235  while (strcasecmp(tmp.c_str(), "ZONE") != 0 && !input_D.eof() ) {
236  input_D >> tmp;
237  Output::echo(2, "%s ", tmp.c_str());
238  if (input_D.eof()) {
239  throw error_msg("DXX_LOAD_GRID_ERR", "Keyword 'ZONE' not found in the diffusion coefficient file: %s\n", D_filename.c_str());
240  }
241  }
242  // read to the end of the line with 'zone'
243  input_D.ignore(9999, '\n');
244 
245  if (!input_D.eof()) {
246  // load grid
247  if (gridOrder == "IFT_GRID_LPA") {
248  for (il = 0; il < L.size; il++)
249  for (im = 0; im < epc.size; im++)
250  for (ia = 0; ia < alpha.size; ia++) {
251  input_D >> Load_L >> Load_epc >> Load_alpha;
252  // load D
253  //input_D >> (*this)[il][im][ia];
254  input_D >> arr[il][im][ia];
255  // check if grid is the same
256  if (fabs(Load_L - L.arr[il][im][ia]) > err || fabs(Load_epc - epc.arr[il][im][ia]) > err || fabs(Load_alpha - alpha.arr[il][im][ia]) > err) {
257  throw error_msg("DXX_LOAD_GRID_ERR", "Loading %s: grid mismatch.\nLoaded: L=%e, epc=%e, alpha=%e\nGrid: L=%e, epc=%e, alpha=%e\n", D_filename.c_str(), Load_L, Load_epc, Load_alpha, L.arr[il][im][ia], epc.arr[il][im][ia], alpha.arr[il][im][ia]);
258  }
259  }
260  } else if (gridOrder == "IFT_GRID_APL") {
261  for (ia = 0; ia < alpha.size; ia++)
262  for (im = 0; im < epc.size; im++)
263  for (il = 0; il < L.size; il++) {
264  input_D >> Load_alpha >> Load_epc >> Load_L ;
265  // load D
266  //input_D >> (*this)[il][im][ia];
267  input_D >> arr[il][im][ia];
268  // check if grid is the same
269  if (fabs(Load_L - L.arr[il][im][ia]) > err || fabs(Load_epc - epc.arr[il][im][ia]) > err || fabs(Load_alpha - alpha.arr[il][im][ia]) > err) {
270  throw error_msg("DXX_LOAD_GRID_ERR", "Loading %s: grid mismatch.\nLoaded: L=%e, epc=%e, alpha=%e\nGrid: L=%e, epc=%e, alpha=%e\n", D_filename.c_str(), Load_L, Load_epc, Load_alpha, L.arr[il][im][ia], epc.arr[il][im][ia], alpha.arr[il][im][ia]);
271  }
272  }
273  } else {
274  throw error_msg("UNKNOWN_GRID_TYPE", "Unknown grid order: %s", gridOrder.c_str());
275  }
276  }
277  input_D.close();
278  return true;
279 
280 }
281 
285 void DiffusionCoefficient::MakeDLL(GridElement &L, GridElement &pc, GridElement &alpha, double Kp, string DLLType) {
286  if (DLLType == "DLLT_BE") {
287  MakeDLL_BE(L, pc, alpha, Kp);
288  }
289  else if (DLLType == "DLLT_B") {
290  MakeDLL_B(L, pc, alpha, Kp);
291  }
292  else if (DLLType == "DLLT_FAKE") {
293  MakeDLL_FAKE(L, pc, alpha, Kp);
294  }
295  else if (DLLType == "DLLT_Ozeke") {
296  MakeDLL_Ozeke(L, pc, alpha, Kp);
297  }
298  else if (DLLType == "DLLT_Ozeke_ME") {
299  MakeDLL_Ozeke_ME(L, pc, alpha, Kp);
300  }
301  else {
302  throw error_msg("DLL_ERROR", "Unknown DLL type");
303  }
304  arr.change_ind = clock();
305 
306 }
307 
312  int il, im, ia;
313  for (il = 0; il < L.size; il++) {
314  for (im = 0; im < pc.size; im++) {
315  for (ia = 0; ia < alpha.size; ia++) {
316  arr[il][im][ia] = Df(L.arr[il][im][ia], Kp);
317  //(*this)[il][im][ia] = 0.0;
318  }
319  }
320  }
321 }
322 
327  int il, im, ia;
328  for (il = 0; il < L.size; il++) {
329  for (im = 0; im < pc.size; im++) {
330  for (ia = 0; ia < alpha.size; ia++) {
331  arr[il][im][ia] = Df(L.arr[il][im][ia], Kp) + 3e3*pow(10.0, 0.506*Kp-9.325)*pow(L.arr[il][im][ia],5);;
332  }
333  }
334  }
335 }
336 
341  int il, im, ia;
342 
343  for (il = 0; il < L.size; il++) {
344  for (im = 0; im < pc.size; im++) {
345  for (ia = 0; ia < alpha.size; ia++) {
346 
347  // electrostatic radial diffusion coefficient in 1/day
348  arr[il][im][ia] = VF::Dfe(alpha.arr[il][im][ia],VF::Kfunc(pc.arr[il][im][ia]), L.arr[il][im][ia],Kp);
349  //(*this)[il][im][ia] = 0.2*VF::Dfe(alpha[il][im][ia],VF::Kfunc(pc[il][im][ia]),L[il][im][ia],Kp);
350 
351  // Add magnetostatic radial diffusion coefficient in 1/day
352  arr[il][im][ia] = arr[il][im][ia] + Df(L.arr[il][im][ia], Kp);
353 
354  }
355  }
356  }
357 }
358 
359 
364  int il, im, ia;
365 
366  for (il = 0; il < L.size; il++) {
367  for (im = 0; im < pc.size; im++) {
368  for (ia = 0; ia < alpha.size; ia++) {
369  arr[il][im][ia] = VF::max(Df(L.arr[il][im][ia], Kp), 1.0/100);
370  }
371  }
372  }
373 
374 }
375 
385  int il, im, ia;
386  for (il = 0; il < L.size; il++) {
387  for (im = 0; im < pc.size; im++) {
388  for (ia = 0; ia < alpha.size; ia++) {
389  arr[il][im][ia] = Df_Ozeke(L.arr[il][im][ia], Kp); // this.arr
390  }
391  }
392  }
393 }
394 
404  int il, im, ia;
405  for (il = 0; il < L.size; il++) {
406  for (im = 0; im < pc.size; im++) {
407  for (ia = 0; ia < alpha.size; ia++) {
408  arr[il][im][ia] = Df_Ozeke(L.arr[il][im][ia], Kp) + Df_Ozeke_E(L.arr[il][im][ia], Kp); // this.arr
409  }
410  }
411  }
412 }
413 
414 
418 bool DiffusionCoefficientsGroup::ActivateAndScale(double time, double Kp) {
419  // The last version of activated stuff.
420  // Activating and scaling diffusion coefficients for different wave types
421  double scaleCoeff = 0;
422  unsigned int Dxx_it;
423  // \todo check, if we need to make any changes, otherwise - just exit
424  // check if we really need to do any changes
425  bool need_changes = false;
426  for (Dxx_it = 0; Dxx_it < DxxList.size(); Dxx_it++) {
427  if (time >= DxxList[Dxx_it].DxxParameters.time_start && time < DxxList[Dxx_it].DxxParameters.time_end) {
428  // if inactive - mark as active
429  if (!DxxList[Dxx_it].is_active) {
430  need_changes = true;
431  break;
432  }
433  // if this Dxx needs to be scaled - scale it
434  if (DxxList[Dxx_it].useScale) {
435  need_changes = true;
436  break;
437  }
438  } else {
439  // if active - mark as inactive
440  if (DxxList[Dxx_it].is_active) {
441  need_changes = true;
442  break;
443  }
444  }
445  }
446  if (need_changes == false) return false;
447 
448  CurrentDxx.arr = 0;
449  // Indicator, that we changed the coefficient (and now we probably need to recalculate the model matrix)
450  CurrentDxx.arr.change_ind = clock();
451  // loop thrue list of the diffusion coefficients
452  for (Dxx_it = 0; Dxx_it < DxxList.size(); Dxx_it++) {
453  // if it is time to work - activate
454  if (time >= DxxList[Dxx_it].DxxParameters.time_start && time < DxxList[Dxx_it].DxxParameters.time_end) {
455  // if inactive - mark as active
456  if (!DxxList[Dxx_it].is_active) {
457  DxxList[Dxx_it].is_active = true;
458  Output::echo(2, "Activated %s.\n", DxxList[Dxx_it].arr.name.c_str());
459  }
460  // if this Dxx needs to be scaled - scale it
461  if (DxxList[Dxx_it].useScale) {
462  // getting scale coefficient
463  scaleCoeff = DxxList[Dxx_it].Scale(Kp);
464  Output::echo(2, "Bw scaled %s by %f.\n", DxxList[Dxx_it].arr.name.c_str(), scaleCoeff);
465  // scalling and adding value to the actualized Dxx which is used in calculations
466  CurrentDxx.arr = CurrentDxx.arr + DxxList[Dxx_it].arr * scaleCoeff;
467  } else {
468  CurrentDxx.arr = CurrentDxx.arr + DxxList[Dxx_it].arr;
469  }
470  } else {
471  // if active - mark as inactive
472  if (DxxList[Dxx_it].is_active) {
473  DxxList[Dxx_it].is_active = false;
474  Output::echo(2, "Deactivated %s.", DxxList[Dxx_it].arr.name.c_str());
475  }
476  }
477  }
478  return false;
479 }
480 
484 double DiffusionCoefficient::Scale(double Kp) {
485  double Bw2=0, Bw2_0=0;
486  if (!this->useScale) return 1;
487  if (DxxParameters.waveType == "WT_CHORUS_DAY" || DxxParameters.waveType == "WT_CHORUS_NIGHT") {
488  // scaling for chorus waves
489  if (DxxParameters.DxxKp <= 2.333) Bw2_0 = 2.0*pow(10, 0.73 + 0.91 * DxxParameters.DxxKp);
490  if (DxxParameters.DxxKp > 2.333) Bw2_0 = 2.0*pow(10, 2.5 + 0.18 * DxxParameters.DxxKp);
491  if (Kp <= 2.333) Bw2 = 2.0*pow(10, 0.73 + 0.91 * Kp);// * 1.e-12 * 1.e8;
492  if (Kp > 2.333) Bw2 = 2.0*pow(10, 2.5 + 0.18 * Kp);// * 1.e-12 * 1.e8;
494  /* Ughm? Commented by Dmitriy
495  if(DxxParameters.waveType == "WT_CHORUS_DAY") {
496  return (Bw2/Bw2_0)*2.5;
497  } else {
498  return (Bw2/Bw2_0)/2.5;
499  }*/
501  return (Bw2/Bw2_0);
502  } else if (DxxParameters.waveType == "WT_HISS_PP" || DxxParameters.waveType == "WT_HISS_PLUME") {
503  // scale hiss diffusion coefficients is a square of Bw scale
504  // return pow(Kp/DxxParameters.DxxKp, 2);
505  return pow(Kp/DxxParameters.DxxKp, 2);
506  }
507  else if (DxxParameters.waveType == "WT_HISS_REAL") {
508  // Daa(Kp)=Daa(Kp=4)*10^{g(Kp)-g(Kp=4)}, where g(Kp)=0.2607*Kp - 0.01547*Kp*Kp.
509  //return pow( (Kp*3 + 26 ) / (DxxParameters.DxxKp*3 + 26), 2);
510  return pow(10, 0.2607*Kp - 0.01547*Kp*Kp - 0.2607*DxxParameters.DxxKp + 0.01547*DxxParameters.DxxKp*DxxParameters.DxxKp);
511  } else if (DxxParameters.waveType == "WT_EMIC_PLUME" ) {
512  // scale EMIC diffusion coefficients as a Bw to the forth power
513  return pow(Kp/DxxParameters.DxxKp, 4);
515  } else if (DxxParameters.waveType == "WT_LIGHTNING") {
516  return pow(Kp/DxxParameters.DxxKp, 2);
517  //
519 
520  } else {
521  //throw error_msg("DXX_SCALING", "No scaling for wave type %s", DxxParameters.waveType.c_str());
522  return 1;
523  }
524 }
525 
538 double Alpha_ne(double pangle, double lambda, double L) {
539  return asin( sqrt( f1(lambda)) * sin(pangle) );
540 }
541 
546 double f1(double lambda) {
547  return sqrt(4.0 - 3.0 * pow(cos(lambda), 2))/(pow(cos(lambda), 6));
548 }
549 
554 double B(double lambda, double L) {
555  return f1(lambda) / pow(L, 3) * VC::B_0;
556 }
557 
561 double func_tmp (double x, double Alpha) {
562  return pow(x, 6) + (3.0 * pow(sin(Alpha), 4)) * x - 4.0 * pow(sin(Alpha), 4);
563 }
564 
571 double F_cap(double x, double y, double b, double s, double epsilon, DxxParameters_structure DxxParameters) {
572  double c1 = 2.0*s*(-1+epsilon);
573  double c2 = 1.0-4.0*epsilon+pow(epsilon,2);
574  double c3 = -s*(-1+epsilon)*(b+4.0*epsilon)/2.0;
575  double c4 = epsilon*(b+epsilon) ;
576  double g = pow(x,4) + c1*pow(x,3) + c2*pow(x,2) + c3*x + c4;
577  return y*(pow((x-s),2))*(pow((x+s*epsilon),2))/(x*g);
578 }
579 
590 double F_cap2(double x, double y, double a, double beta, double mu, double s, double epsilon, double Alpha_star, DxxParameters_structure DxxParameters) {
591  double up = 2.0*(x+a)/(beta*mu);
592  double down = 2.0*x-1.0/Alpha_star*(1.0/(x-s)+DxxParameters.eta1*epsilon/(x+s*epsilon)+DxxParameters.eta2*epsilon/(4.0*x+s*epsilon)+DxxParameters.eta3*epsilon/(16.0*x+s*epsilon))+
593  x*1.0/Alpha_star*(1.0/(pow((x-s), 2))+DxxParameters.eta1*epsilon/(pow((x+s*epsilon), 2))+4.0*epsilon*DxxParameters.eta2/(pow((4.0*x+s*epsilon), 2))+16.0*epsilon*DxxParameters.eta3/(pow((16.0*x+s*epsilon), 2)));
594  return up/down;
595 }
596 
602 double quad1(double (*func)(double lambda, DxxParameters_structure DxxParameters), double a, double b, int M, DxxParameters_structure DxxParameters) {
603  // f -integrated function
604  // a and b upper and lower limits
605  // M number of subintegrals
606  double h = (b - a)/M;
607  double s = func(a, DxxParameters);
608  double x;
609  int k;
610  for (k = 1; k <= (M-1); k++) {
611  x = a + h*k;
612  s = s + func(x, DxxParameters);
613  }
614  s = h * (func(a, DxxParameters) + func(b, DxxParameters))/2.0 + h*s;
615  return s;
616 }
617 
623 double Dxx_ba(double L1, double EMeV, double Alpha, double int_Dxx_loc(double lambda, DxxParameters_structure DxxParameters), DxxParameters_structure DxxParameters) {
624  // Searching for zero of function func_tmp=x.^6+(3*sin(DxxParameters.Alpha_eq_tmp)^4.)*x-4*sin(DxxParameters.Alpha_eq_tmp)^4.;
625  // near the point 0.5
627  DxxParameters.L = L1;
628  DxxParameters.EMeV = EMeV;
629  DxxParameters.Alpha = Alpha;
630 
631  double lambda_m;
632  //lambda_m = acos(sqrt(fzero(func_tmp, 0.5)))*0.99;
633  lambda_m = acos(sqrt(bisection (func_tmp, DxxParameters.Alpha, 0, 1)))*DxxParameters.mirror_point_coeff;
634  if (lambda_m > DxxParameters.lam_max*VC::pi/180) lambda_m = DxxParameters.lam_max*VC::pi/180.;
635  // quad1 - calculating integral
636  return (1.0/(1.38-0.32*(sin(DxxParameters.Alpha) + pow(sin(DxxParameters.Alpha),0.5))) * quad1(int_Dxx_loc, DxxParameters.lam_min*VC::pi/180., lambda_m, DxxParameters.nint, DxxParameters));
637 }
638 
646 double int_Daa_loc(double lambda, DxxParameters_structure DxxParameters) {
647  //double Daa = Daa_local(lambda);
648  double Daa = Dxx_local(lambda, Daa_root, DxxParameters);
649  if (DxxParameters.Alpha*180/VC::pi < 0.001) Daa = 0;
650  return Daa * cos(Alpha_ne(DxxParameters.Alpha, lambda, DxxParameters.L)) * pow(cos(lambda), 7) /(pow(cos(DxxParameters.Alpha), 2));
651 }
652 
660 double int_Dpp_loc(double lambda, DxxParameters_structure DxxParameters) {
661  //double Dpp = Dpp_local(lambda);
662  double Dpp = Dxx_local(lambda, Dpp_root, DxxParameters);
663  if (DxxParameters.Alpha*180/VC::pi < 0.001) Dpp = 0;
664  return Dpp * (cos(lambda) * pow((1 + 3*pow(sin(lambda), 2)), 0.5)) / cos(Alpha_ne(DxxParameters.Alpha, lambda, DxxParameters.L));
665 }
666 
674 double int_Dpa_loc(double lambda, DxxParameters_structure DxxParameters) {
675  //double Dpa = Dpa_local(lambda);
676  double Dpa = Dxx_local(lambda, Dpa_root, DxxParameters);
677  if (DxxParameters.Alpha*180/VC::pi < 0.001) Dpa = 0;
678  return Dpa * sin(Alpha_ne(DxxParameters.Alpha, lambda, DxxParameters.L)) * pow(cos(lambda), 7) / (cos(DxxParameters.Alpha) * sin(DxxParameters.Alpha));
679 }
680 
690 double Dxx_local(double lambda, double Dxx_root(double Omega_e, double x, double mu, double su, double y, double beta, double a, double b, double Alpha_star, double s, double epsilon, double d_x, double x_m, double R, DxxParameters_structure DxxParameters), DxxParameters_structure DxxParameters ) {
691 
692  double Omega_e = VC::qe * B(lambda, DxxParameters.L) / VC::m / VC::c_cgs;
693  double Omega_e_eq = VC::qe * B(0, DxxParameters.L) / VC::m / VC::c_cgs; // Omega_e_eq = B_0/(L^3)
694  double N_0 = 124.0 * pow((3.0 / DxxParameters.L), 4);
695  double Omega_p = sqrt(4.0 * VC::pi * N_0 * pow(VC::qe, 2) / VC::m);
696  double Dxx_f = DxxParameters.f; // to hold previous value
697  DxxParameters.f = Omega_p/Omega_e_eq;
698  //double Omega_p;
699  //if (DxxParameters.f > 0) {
700  //double N_0;
701  double LT;
702 
703  // rearranged if statements to be more concise
704  if (DxxParameters.numberDensity == "ND_CONSTANT") {
705  Omega_p = Dxx_f * Omega_e_eq;
706  } else if (DxxParameters.numberDensity == "ND_CHORUS") {
707  // ignore, since parameters already set above
708 // Output::echo(2, "%e \n", DxxParameters.f);
709  } else if (DxxParameters.numberDensity == "ND_CHORUS_NIGHT" || DxxParameters.numberDensity == "ND_CHORUS_DAY") {
710  if (DxxParameters.numberDensity == "ND_CHORUS_NIGHT")
711  LT = 3;
712  else
713  LT = 9;
714 
715  N_0 += 36.0*pow((3.0/DxxParameters.L),3.5) * cos((LT-(7.7*pow((3.0/DxxParameters.L),2.0)+12))*VC::pi/12);
716 // if (DxxParameters.L >= 4.5) { Output::echo(2, "%e \n", DxxParameters.f); }
717  } else if (DxxParameters.numberDensity == "ND_HISS" || DxxParameters.numberDensity == "ND_PL_CA_STARKS") {
718  N_0 = pow(10, (-0.3145*DxxParameters.L + 3.9043));
719  if (DxxParameters.numberDensity == "ND_PL_CA_STARKS" && DxxParameters.L <= 2)
720  {
721  double N1, N2, L1, L2, A, B;
722  N1=10-6.;N2=8.3-6.;L1=1.2;L2=5.9;
723  A=(N2-N1)/(L2-L1); B=N2-A*L2;
724  N_0 = pow(10, A*DxxParameters.L+B);
725  }
726 // Output::echo(2, "%e \n", DxxParameters.f);
727  } else if (DxxParameters.numberDensity == "ND_HISS2") {
728  N_0 = 1e3*pow(4.0/DxxParameters.L, 4);
729 // Output::echo(2, "%e \n", DxxParameters.f);
730  } else {
731  throw error_msg("DXX_DENSITY_ERROR", "Unknown number density type: %s\n", DxxParameters.numberDensity.c_str());
732  }
733 
734 
735  double epsilon = 1.0/1837;
736  double beta = (pow((DxxParameters.EMeV * (DxxParameters.EMeV + 2)), 0.5)) / (DxxParameters.EMeV + 1);
737  double mu = cos(Alpha_ne(DxxParameters.Alpha, lambda, DxxParameters.L));
738 
739  // sin version of mu, returns sqrt( f1(lambda)) * sin(alpha)
740  double su = sin(Alpha_ne(DxxParameters.Alpha, lambda, DxxParameters.L));
741 
742  double lam_sp = -1; // -1 corresponds to electrons
743  double gamma = DxxParameters.EMeV + 1;
744  double a = DxxParameters.s * lam_sp / gamma;
745 
746  double Alpha_star = pow((Omega_e / Omega_p), 2);
747  double b = (1 + epsilon)/Alpha_star;
748 
749  //double d_omega = d_omega_per * Omega_e_eq;
750  //double Omega_m = Omega_m_per * Omega_e_eq;
751  double d_omega;
752  double Omega_m;
753  double Omega_1; // lower cutoff
754  double Omega_2; // upper cutoff
755 
756  // to determine partial or absolute omega type
757  if (DxxParameters.Omega_mType == "OT_PARTIAL") {
758  Omega_m = DxxParameters.Omega_m * Omega_e_eq;
759  d_omega = DxxParameters.d_omega * Omega_e_eq ;
760  Omega_1 = DxxParameters.omega_lc * Omega_e_eq;
761  Omega_2 = DxxParameters.omega_uc * Omega_e_eq;
762  } else if (DxxParameters.Omega_mType == "OT_ABSOLUTE") {
763  Omega_m = DxxParameters.Omega_m * 2 * VC::pi;
764  d_omega = DxxParameters.d_omega * 2 * VC::pi;
765  Omega_1 = DxxParameters.omega_lc * 2 * VC::pi;
766  Omega_2 = DxxParameters.omega_uc * 2 * VC::pi;
767  } else {
768  throw error_msg("DXX_OMEGA_M_TYPE", "Unknown Omega_m type: ", DxxParameters.Omega_mType.c_str());
769  }
770 
771  // if (DxxParameters.Omega_mType == "partial") {
772  // } else if (DxxParameters.Omega_mType == "constant") {
773  // } else {
774  // }
775 
776  //double R = pow((DxxParameters.dB / B(lambda, DxxParameters.L)), 2);
777  double Bw;
778  if (DxxParameters.BwFromLambda) {
779  Bw = pow(10, 0.75 + 0.04*(lambda*180/VC::pi))*1.e-8 * DxxParameters.Bw; // should be gauss...
780  } else {
781  Bw = DxxParameters.Bw;
782  }
783  double R = pow((Bw / B(lambda, DxxParameters.L)), 2);
784 
785  double x_m=0;
786  double d_x=0;
787  double x_1=0;
788  double x_2=0;
789  // determine particle type
790  if (DxxParameters.particle == "PT_ELECTRON") {
791  x_m = Omega_m / Omega_e;
792  d_x = d_omega / Omega_e;
793  x_1 = Omega_1 / Omega_e;
794  x_2 = Omega_2 / Omega_e;
795  } else if (DxxParameters.particle == "PT_PROTON") {
796  double cfm = 16; // cubic feet per minute(?)
797  x_m = Omega_m / Omega_e * epsilon / cfm;
798  d_x = d_omega / Omega_e * epsilon / cfm;
799  x_1 = Omega_1 / Omega_e * epsilon / cfm;
800  x_2 = Omega_2 / Omega_e * epsilon / cfm;
801  }
802 
803  std::vector<double> xx = rrouts(x_1, x_2, DxxParameters.eta1, DxxParameters.eta2, DxxParameters.eta3, epsilon, beta, mu, Alpha_star, a, DxxParameters);
804  //
805  double x, y;
806  double Dxx = 0, Dxx_tmp = 0;
807  for (unsigned int i = 1; i <= xx.size();i++) {
808  x = xx[i-1];
809  y = (x + a) / (beta * mu);
810  Dxx_tmp = Dxx_root(Omega_e, x, mu, su, y, beta, a, b, Alpha_star, DxxParameters.s, epsilon, d_x, x_m, R, DxxParameters);
811  Dxx = Dxx_tmp + Dxx;
812  }
813  if (xx.size() == 0) {
814  Dxx = 0;
815  }
816 
817  return Dxx;
818 }
819 
830 double Daa_root(double Omega_e, double x, double mu, double su, double y, double beta, double a, double b, double Alpha_star, double s, double epsilon, double d_x, double x_m, double R, DxxParameters_structure DxxParameters) {
831  double fcap = F_cap2(x,y,a,beta,mu,s,epsilon,Alpha_star, DxxParameters);
832  return (VC::pi/2) * (1.0/DxxParameters.nu) * Omega_e * (1.0/(pow((DxxParameters.EMeV + 1), 2))) *
833  ((R * pow((1 - x * mu / (y * beta)), 2) * fabs(fcap))/
834  (d_x * fabs(beta * mu - fcap)))*
835  exp( - (pow(((x - x_m) / d_x), 2)));
836 }
837 
848 double Dpa_root(double Omega_e, double x, double mu, double su, double y, double beta, double a, double b, double Alpha_star, double s, double epsilon, double d_x, double x_m, double R, DxxParameters_structure DxxParameters) {
849  double fcap = F_cap2(x,y,a,beta,mu,s,epsilon,Alpha_star,DxxParameters);
850  return -(VC::pi/2)*(1.0/DxxParameters.nu)* Omega_e*su/beta*(1.0/(pow((DxxParameters.EMeV+1), 2)))*
851  ((R*(1-x/y*mu/beta)*(x/y)* fabs(fcap))/
852  (d_x * fabs(beta * mu - fcap)))*
853  exp( -(pow(((x - x_m) / d_x), 2)));
854 }
855 
866 double Dpp_root(double Omega_e, double x, double mu, double su, double y, double beta, double a, double b, double Alpha_star, double s, double epsilon, double d_x, double x_m, double R, DxxParameters_structure DxxParameters) {
867  double fcap = F_cap2(x,y,a,beta,mu,s,epsilon,Alpha_star,DxxParameters);
868  return (VC::pi/2) * (1.0/DxxParameters.nu) * (1./pow(beta, 2)) * (1-pow(mu, 2)) * Omega_e * (1.0/(pow((DxxParameters.EMeV + 1), 2))) *
869  ((R * (pow(x/y, 2)) * fabs(fcap))/
870  (d_x*fabs(beta*mu - fcap)))*
871  exp(-(pow(((x-x_m)/d_x),2)));
872 }
873 
874 
875 //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) {
879 std::vector<double> rrouts(double x_1, double x_2, double eta1, double eta2, double eta3, double epsilon, double beta, double mu, double Alpha_star, double a, DxxParameters_structure DxxParameters) {
880 
881  //double Omega_1 = DxxParameters.omega_lc * Omega_e_eq * epsilon / DxxParameters.cfm;
882  //double Omega_2 = DxxParameters.omega_uc * Omega_e_eq * epsilon / DxxParameters.cfm;
883  //double x_1 = Omega_1 / Omega_e;
884  //double x_2 = Omega_2 / Omega_e;
885 
886  // bunch of coefficients of different orders
887  double c0 = -DxxParameters.s;
888  double c1 = DxxParameters.s * epsilon;
889  double c2 = 1.0/4. * DxxParameters.s * epsilon;
890  double c3 = 1.0/16. * DxxParameters.s * epsilon;
891 
892  double d0 = 1;
893  double d1 = eta1*epsilon;
894  double d2 = eta2/4.*epsilon;
895  double d3 = eta3/16.*epsilon;
896 
897  double g3 = d0 + d1 + d2 + d3;
898  double g2 = d0*(c1+c2+c3) + d1*(c0+c2+c3) + d2*(c0+c1+c3) + d3*(c0+c1+c2);
899  double g1 = d0*(c1*c2 + c1*c3 + c2*c3) + d1*(c0*c2 + c0*c3 + c2*c3) + d2*(c0*c1 + c0*c3 + c1*c3) + d3*(c0*c1 + c0*c2 + c1*c2);
900  double g0 = d0*c1*c2*c3 + d1*c0*c2*c3 + d2*c0*c1*c3 + d3*c0*c1*c2;
901 
902  double h4 = 1.;
903  double h3 = c0 + c1 + c2 + c3;
904  double h2 = c0*c1 + c0*c2 + c0*c3 + c1*c2 + c1*c3 + c2*c3;
905  double h1 = c0*c1*c2 + c0*c1*c3 + c0*c2*c3 + c1*c2*c3;
906  double h0 = c0*c1*c2*c3;
907 
908  // array of calculated coefficients of polynomial
909  double A[7];
910  A[0] = (pow(beta,2)*pow(mu,2)-1.)*h4;
911  A[1] = (pow(beta,2)*pow(mu,2)-1.)*h3-2*a*h4;
912  A[2] = (pow(beta,2)*pow(mu,2)-1.)*h2-2*a*h3-pow(beta,2)*pow(mu,2)*1./Alpha_star*g3-pow(a,2)*h4;
913  A[3] = (pow(beta,2)*pow(mu,2)-1.)*h1-2*a*h2-pow(beta,2)*pow(mu,2)*1./Alpha_star*g2-pow(a,2)*h3;
914  A[4] = (pow(beta,2)*pow(mu,2)-1.)*h0-2*a*h1-pow(beta,2)*pow(mu,2)*1./Alpha_star*g1-pow(a,2)*h2;
915  A[5] = -2*a*h0-pow(beta,2)*pow(mu,2)*1./Alpha_star*g0-pow(a,2)*h1;
916  A[6] = -pow(a,2)*h0;
917 
918  //a6=(beta^2*mu^2-1.)*h4;
919  //a5=(beta^2*mu^2-1.)*h3-2*a*h4;
920  //a4=(beta^2*mu^2-1.)*h2-2*a*h3-beta^2*mu^2*1./Alpha_star*g3-a^2*h4;
921  //a3=(beta^2*mu^2-1.)*h1-2*a*h2-beta^2*mu^2*1./Alpha_star*g2-a^2*h3;
922  //a2=(beta^2*mu^2-1.)*h0-2*a*h1-beta^2*mu^2*1./Alpha_star*g1-a^2*h2;
923  //a1=-2*a*h0-beta^2*mu^2*1./Alpha_star*g0-a^2*h1;
924  //a0=-a^2*h0;
925 
926  std::vector<double> xx;
927  // r is to store quadratic factors, wr is real roots, wi is imaginary roots, quad used for quadratic calculation
928  double r[7], wr[6],wi[6],quad[2];
929  int n, numr; // number of roots
930  n=6;
931  // initialize estimate for 1st root pair
932  quad[0] = 2.71828e-1;
933  quad[1] = 3.14159e-1;
934  // get roots
935  get_quads(A, n, quad, r);
936  numr = roots(r, n, wr, wi);
937 
938  for (int i = 0; i < numr; i++) {
939  if (fabs(wi[i]) < double_zero && wr[i] < x_2 && wr[i] > x_1) {
940  xx.push_back(wr[i]);
941  }
942  }
943 
944  return xx;
945 
946 }
947 
948 
952 double Df(double L, double Kp) {
953  return pow(10, 0.506*Kp-9.325)*pow(L,10);
954 }
963 double Df_Ozeke(double L, double Kp) {
964 
965  return 6.62e-13 * pow(L,8) * pow(10, (pow(L,2) * -0.0327 + L * 0.625 - pow(Kp,2) * 0.0108 + Kp * 0.499) );
966 }
967 
970 double Df_Ozeke_E(double L, double Kp) {
971 
972  return 2.16e-8 * pow(L,6) * pow(10,( 0.217 * L + 0.461 * Kp ));
973 }
974 
975 
976 
977 /* DEPRECATED
979 // http://irbem.svn.sourceforge.net/viewvc/irbem/extras/opendc/
980 double Df_Ozeke(double L, double Kp) {
981 
982  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 );
983 }
984 
985 double Df_Ozeke_E(double L, double Kp) {
986 
987  return 5.75e-9 * pow(L,6) * pow(10,( 0.208 * L + 0.457 * Kp ));
988 }
989 */
990 
1008  Parameters_structure &parameters, Grid &radialDiffusionGrid, Grid &localDiffusionsGrid) {
1009 
1010  // Diffusion coefficient group (Lists of similar diffusion coefficients, like Daa or Dpcpc...)
1011  DLL.arr.AllocateMemory(radialDiffusionGrid.L.size, radialDiffusionGrid.pc.size, radialDiffusionGrid.alpha.size);
1012  Daa.CurrentDxx.arr.AllocateMemory(radialDiffusionGrid.L.size, radialDiffusionGrid.pc.size, radialDiffusionGrid.alpha.size);
1013  DaaLpp.CurrentDxx.arr.AllocateMemory(radialDiffusionGrid.L.size, radialDiffusionGrid.pc.size, radialDiffusionGrid.alpha.size);
1014  Dpcpc.CurrentDxx.arr.AllocateMemory(radialDiffusionGrid.L.size, radialDiffusionGrid.pc.size, radialDiffusionGrid.alpha.size);
1015  DpcpcLpp.CurrentDxx.arr.AllocateMemory(radialDiffusionGrid.L.size, radialDiffusionGrid.pc.size, radialDiffusionGrid.alpha.size);
1016  Dpca.CurrentDxx.arr.AllocateMemory(radialDiffusionGrid.L.size, radialDiffusionGrid.pc.size, radialDiffusionGrid.alpha.size);
1017  DpcaLpp.CurrentDxx.arr.AllocateMemory(radialDiffusionGrid.L.size, radialDiffusionGrid.pc.size, radialDiffusionGrid.alpha.size);
1018 
1019  // Make it zero
1020  DLL.arr = 0;
1021  Daa.CurrentDxx.arr = 0;
1022  DaaLpp.CurrentDxx.arr = 0;
1023  Dpcpc.CurrentDxx.arr = 0;
1024  DpcpcLpp.CurrentDxx.arr = 0;
1025  Dpca.CurrentDxx.arr = 0;
1026  DpcaLpp.CurrentDxx.arr = 0;
1027 
1028  // Temporary diffusion coefficient
1029  DiffusionCoefficient Dxx_tmp;
1030  // Create the diffusion coefficients for each of the Dxx parameters files
1031  unsigned int Dxx_it;
1032  for (Dxx_it = 0; Dxx_it < parameters.DxxParametersList.size(); Dxx_it++) {
1033 
1034  Dxx_tmp.Create_Dxx(parameters.DxxParametersList[Dxx_it], localDiffusionsGrid);
1035 
1036  // Add the created coefficient to an appropriate group
1037  if (parameters.useAlphaDiffusion) {
1038  if (Dxx_tmp.type == "DCT_Daa")
1039  Daa.DxxList.push_back(Dxx_tmp);
1040  else if (Dxx_tmp.type == "DCT_DaaLpp")
1041  DaaLpp.DxxList.push_back(Dxx_tmp);
1042  }
1043 
1044  if (parameters.useEnergyDiffusion) {
1045  if (Dxx_tmp.type == "DCT_Dpp" || Dxx_tmp.type == "DCT_Dpcpc")
1046  Dpcpc.DxxList.push_back(Dxx_tmp);
1047  else if (Dxx_tmp.type == "DCT_DppLpp" || Dxx_tmp.type == "DCT_DpcpcLpp")
1048  DpcpcLpp.DxxList.push_back(Dxx_tmp);
1049  }
1050 
1051  if (parameters.useEnergyAlphaMixedTerms) {
1052  if (Dxx_tmp.type == "DCT_Dpca" || Dxx_tmp.type == "DCT_Dpa")
1053  Dpca.DxxList.push_back(Dxx_tmp);
1054  else if (Dxx_tmp.type == "DCT_DpcaLpp" || Dxx_tmp.type == "DCT_DpaLpp")
1055  DpcaLpp.DxxList.push_back(Dxx_tmp);
1056 
1057  }
1058  }
1059  // Check if Daa*Dpcpc > Dpca^2 for each Dpca
1060  if (parameters.useEnergyAlphaMixedTerms) {
1061  int Daa_it, Dpcpc_it, Dpca_it;
1062  // above Lpp
1063  for (Dpca_it = 0; Dpca_it < Dpca.DxxList.size(); Dpca_it++) {
1064  // For each Dpca...
1065  // Search for Daa with the same wave-type
1066  for (Daa_it = 0; Daa_it < Daa.DxxList.size(); Daa_it++) {
1067  if (Daa.DxxList[Daa_it].DxxParameters.waveType == Dpca.DxxList[Dpca_it].DxxParameters.waveType)
1068  break;
1069  }
1070  // Search for Dpp with the same wave-type
1071  for (Dpcpc_it = 0; Dpcpc_it < Dpcpc.DxxList.size(); Dpcpc_it++) {
1072  if (Dpcpc.DxxList[Dpcpc_it].DxxParameters.waveType == Dpca.DxxList[Dpca_it].DxxParameters.waveType)
1073  break;
1074  }
1075  // Check if Daa*Dpcpc > Dpca^2 for each Dpca
1076  int iL, im, ia;
1077  for (iL = 0; iL < localDiffusionsGrid.L.size; iL++)
1078  for (im = 0; im < localDiffusionsGrid.pc.size; im++)
1079  for (ia = 0; ia < localDiffusionsGrid.alpha.size; ia++)
1080  if (Daa.DxxList[Daa_it].arr[iL][im][ia] * Dpcpc.DxxList[Dpcpc_it].arr[iL][im][ia] * 0.99 <= Dpca.DxxList[Dpca_it].arr[iL][im][ia] * Dpca.DxxList[Dpca_it].arr[iL][im][ia])
1081  Output::echo(0, "WARNING: Diffusion coefficient error. %s: Dpp * Daa * 0.99 <= Dpa^2. \nAt (%d,%d,%d): L=%f, epc=%f, alpha=%f; Dpp=%e, Daa=%e, Dpa=%e\n",
1082  Dpca.DxxList[Dpca_it].DxxParameters.waveType.c_str(),
1083  iL, im, ia,
1084  localDiffusionsGrid.L.arr[iL][im][ia], localDiffusionsGrid.epc.arr[iL][im][ia], localDiffusionsGrid.alpha.arr[iL][im][ia],
1085  Dpcpc.DxxList[Dpcpc_it].arr[iL][im][ia], Daa.DxxList[Daa_it].arr[iL][im][ia], Dpca.DxxList[Dpca_it].arr[iL][im][ia]);
1086  }
1087 
1088  // Below Lpp
1089  for (Dpca_it = 0; Dpca_it < DpcaLpp.DxxList.size(); Dpca_it++) {
1090  // For each Dpca...
1091  // Search for Daa with the same wave-type
1092  for (Daa_it = 0; Daa_it < DaaLpp.DxxList.size(); Daa_it++) {
1093  if (DaaLpp.DxxList[Daa_it].DxxParameters.waveType == DpcaLpp.DxxList[Dpca_it].DxxParameters.waveType)
1094  break;
1095  }
1096  // Search for Dpp with the same wave-type
1097  for (Dpcpc_it = 0; Dpcpc_it < DpcpcLpp.DxxList.size(); Dpcpc_it++) {
1098  if (DpcpcLpp.DxxList[Dpcpc_it].DxxParameters.waveType == DpcaLpp.DxxList[Dpca_it].DxxParameters.waveType)
1099  break;
1100  }
1101  // Check if Daa*Dpcpc > Dpca^2 for each Dpca
1102  int iL, im, ia;
1103  for (iL = 0; iL < localDiffusionsGrid.L.size; iL++)
1104  for (im = 0; im < localDiffusionsGrid.pc.size; im++)
1105  for (ia = 0; ia < localDiffusionsGrid.alpha.size; ia++)
1106  if (DaaLpp.DxxList[Daa_it].arr[iL][im][ia] * DpcpcLpp.DxxList[Dpcpc_it].arr[iL][im][ia] * 0.99 <= DpcaLpp.DxxList[Dpca_it].arr[iL][im][ia] * DpcaLpp.DxxList[Dpca_it].arr[iL][im][ia])
1107  Output::echo(0, "WARNING: Diffusion coefficient error. %s: Dpp * Daa * 0.99 <= Dpa^2. \nAt (%d,%d,%d): L=%f, epc=%f, alpha=%f; Dpp=%e, Daa=%e, Dpa=%e\n",
1108  DpcaLpp.DxxList[Dpca_it].DxxParameters.waveType.c_str(),
1109  iL, im, ia,
1110  localDiffusionsGrid.L.arr[iL][im][ia], localDiffusionsGrid.epc.arr[iL][im][ia], localDiffusionsGrid.alpha.arr[iL][im][ia],
1111  DpcpcLpp.DxxList[Dpcpc_it].arr[iL][im][ia], DaaLpp.DxxList[Daa_it].arr[iL][im][ia], DpcaLpp.DxxList[Dpca_it].arr[iL][im][ia]);
1112  }
1113 
1114  }
1115 
1116  // calculating DLL If we need to do Radial diffusion
1117  if (parameters.useRadialDiffusion) {
1118  DLL.MakeDLL(radialDiffusionGrid.L, radialDiffusionGrid.pc, radialDiffusionGrid.alpha, parameters.Kp[0], parameters.DLLType);
1119 // DLL.writeToFile("./Debug/DLL.plt", radialDiffusionGrid.L, radialDiffusionGrid.epc, radialDiffusionGrid.alpha);
1120  }
1121 
1122 }
1123 
1124 
1125 
1126 void Output1DHeaders(ofstream &output1D,
1129 
1130  // Title of file (first line)
1131  output1D << "Variables = \"time\", \"Kp\", \"Boundary_fluxes\", \"Lpp\", \"tau\" ";
1132 
1133  // Writing headers for all diffusion coefficients scaling coefficients
1134  unsigned int Dxx_it;
1135  for (Dxx_it = 0; Dxx_it < Daa.DxxList.size(); Dxx_it++) {
1136  output1D << "\"" << Daa.DxxList[Dxx_it].arr.name << "\"" << ", ";
1137  }
1138  for (Dxx_it = 0; Dxx_it < Dpcpc.DxxList.size(); Dxx_it++) {
1139  output1D << "\"" << Dpcpc.DxxList[Dxx_it].arr.name << "\"" << ", ";
1140  }
1141  for (Dxx_it = 0; Dxx_it < Dpca.DxxList.size(); Dxx_it++) {
1142  output1D << "\"" << Dpca.DxxList[Dxx_it].arr.name << "\"" << ", ";
1143  }
1144 
1145  for (Dxx_it = 0; Dxx_it < DaaLpp.DxxList.size(); Dxx_it++) {
1146  output1D << "\"" << DaaLpp.DxxList[Dxx_it].arr.name << "\"" << ", ";
1147  }
1148  for (Dxx_it = 0; Dxx_it < DpcpcLpp.DxxList.size(); Dxx_it++) {
1149  output1D << "\"" << DpcpcLpp.DxxList[Dxx_it].arr.name << "\"" << ", ";
1150  }
1151  for (Dxx_it = 0; Dxx_it < DpcaLpp.DxxList.size(); Dxx_it++) {
1152  output1D << "\"" << DpcaLpp.DxxList[Dxx_it].arr.name << "\"" << ", ";
1153  }
1154 
1155  output1D << endl;
1156  output1D << "ZONE T=\"1d-output\" " << endl;
1157 }
1158 
1159 void Output1DValues(ofstream &output1D,
1162  double &time, Parameters_structure &parameters, int iteration) {
1163 
1164  // Title of file (first line)
1165  output1D << time << "\t" << parameters.Kp[iteration] << "\t" << parameters.Bf[iteration] << "\t" << parameters.Lpp[iteration] << "\t" << parameters.tau[iteration] << "\t"; //parameters.dBw2[iteration] << "\t" <<
1166 
1167  unsigned int Dxx_it;
1168  for (Dxx_it = 0; Dxx_it < Daa.DxxList.size(); Dxx_it++) {
1169  output1D << Daa.DxxList[Dxx_it].Scale(parameters.Kp[iteration]) << "\t";
1170  }
1171  for (Dxx_it = 0; Dxx_it < Dpcpc.DxxList.size(); Dxx_it++) {
1172  output1D << Dpcpc.DxxList[Dxx_it].Scale(parameters.Kp[iteration]) << "\t";
1173  }
1174  for (Dxx_it = 0; Dxx_it < Dpca.DxxList.size(); Dxx_it++) {
1175  output1D << Dpca.DxxList[Dxx_it].Scale(parameters.Kp[iteration]) << "\t";
1176  }
1177 
1178  for (Dxx_it = 0; Dxx_it < DaaLpp.DxxList.size(); Dxx_it++) {
1179  output1D << DaaLpp.DxxList[Dxx_it].Scale(parameters.Kp[iteration]) << "\t";
1180  }
1181  for (Dxx_it = 0; Dxx_it < DpcpcLpp.DxxList.size(); Dxx_it++) {
1182  output1D << DpcpcLpp.DxxList[Dxx_it].Scale(parameters.Kp[iteration]) << "\t";
1183  }
1184  for (Dxx_it = 0; Dxx_it < DpcaLpp.DxxList.size(); Dxx_it++) {
1185  output1D << DpcaLpp.DxxList[Dxx_it].Scale(parameters.Kp[iteration]) << "\t";
1186  }
1187 
1188  output1D << endl;
1189 }
Holds list of instances of DiffusionCoefficient class of same type (like Daa, Dpp, etc), but produced by different waves (Daa_chorus, Daa_EMIC, etc).
GridElement alpha
grid elements to be used
Definition: Grid.h:59
Matrix3D< double > arr
array of grid points
Definition: Grid.h:30
Matrix1D< double > Kp
Kp array.
Definition: Parameters.h:106
double eta2
eta values, calculation parameters
Definition: Parameters.h:66
void add(single_error err)
Definition: error.h:96
double int_Dpa_loc(double lambda, DxxParameters_structure DxxParameters)
static const double m
mass of electron in grams
double Alpha
global
Definition: Parameters.h:78
bool ActivateAndScale(double time, double Kp)
double Daa_root(double Omega_e, double x, double mu, double su, double y, double beta, double a, double b, double Alpha_star, double s, double epsilon, double d_x, double x_m, double R, DxxParameters_structure DxxParameters)
Array of values of coordinate axes.
Definition: Grid.h:28
string filetype
Type of the file (w/wo grid etc)
Definition: Parameters.h:35
bool useScale
Flag, use scaling if equal to true.
Definition: Parameters.h:41
string Omega_mType
Omega_m type. Check StrToVal(string input, Omega_mTypes &place) for known values. ...
Definition: Parameters.h:51
double Dpp_root(double Omega_e, double x, double mu, double su, double y, double beta, double a, double b, double Alpha_star, double s, double epsilon, double d_x, double x_m, double R, DxxParameters_structure DxxParameters)
double max(double v1, double v2)
Return maximum.
void Create_Dxx(DxxParameters_structure DxxParameters, Grid &grid)
Create diffusion coefficients.
int roots(double *a, int n, double *wr, double *wi)
Extract individual real or complex roots from list of quadratic factors.
Definition: rroots.cpp:15
double F_cap(double x, double y, double b, double s, double epsilon, DxxParameters_structure DxxParameters)
void Output1DValues(ofstream &output1D, DiffusionCoefficientsGroup &Daa, DiffusionCoefficientsGroup &Dpcpc, DiffusionCoefficientsGroup &Dpca, DiffusionCoefficientsGroup &DaaLpp, DiffusionCoefficientsGroup &DpcpcLpp, DiffusionCoefficientsGroup &DpcaLpp, double &time, Parameters_structure &parameters, int iteration)
Function used in main() to print values of diffusion coefficients in output file. ...
double Df(double L, double Kp)
bool LoadDiffusionCoefficientFromPlaneFile(GridElement &L, GridElement &pc, GridElement &alpha, string D_filename)
Load Dxx from file without grid.
Holds diffusion coefficient matrix and routines to load and calculate it.
bool useAlphaDiffusion
Using diffusions flags.
Definition: Parameters.h:101
double eta1
eta values, calculation parameters
Definition: Parameters.h:66
double eta3
eta values, calculation parameters
Definition: Parameters.h:66
string particle
Type of particles, produced wave? Ions or electrons. Check StrToVal(string input, ParticleTypes &plac...
Definition: Parameters.h:74
string what()
Definition: error.h:117
double EMeV
global
Definition: Parameters.h:78
bool LoadDiffusionCoefficient(GridElement &L, GridElement &pc, GridElement &alpha, string D_filename, string filetype="IFT_GRID")
double lam_min
Minimum latitude.
Definition: Parameters.h:61
double Alpha_ne(double pangle, double lambda, double L)
DiffusionCoefficient CurrentDxx
flag, indicated that the initialization was passed
static const double c_cgs
speed of light, cm/s
void MakeDLL100(GridElement &L, GridElement &pc, GridElement &alpha, double Kp)
Making DLL.
Matrix1D< double > Bf
Boundary flux array.
Definition: Parameters.h:110
double f
f_ratio
Definition: Parameters.h:71
static const double exp
exponential
void MakeDLL_FAKE(GridElement &L, GridElement &pc, GridElement &alpha, double Kp)
Making DLL.
void echo(int msglvl, const char *format,...)
Basic function! Write the message to the file.
Definition: Output.cpp:44
Struct that holds various parameters to be used for Dxx.
Definition: Parameters.h:27
double func_tmp(double x, double Alpha)
void MakeDLL_Ozeke_ME(GridElement &L, GridElement &pc, GridElement &alpha, double Kp)
Making DLL_ME.
double d_omega
d_omega value.
Definition: Parameters.h:53
void AllocateMemory(int size_x, int size_y, int size_z)
Definition: Matrix.cpp:841
string loadOrCalculate
Load diffusion coefficient or calculate flag.
Definition: Parameters.h:44
double quad1(double(*func)(double lambda, DxxParameters_structure DxxParameters), double a, double b, int M, DxxParameters_structure DxxParameters)
double MLT_averaging
MLT averaging part of the orbit.
Definition: Parameters.h:47
static const double mc2
m*c^2 = 0.511 MeV
double Dxx_local(double lambda, double Dxx_root(double Omega_e, double x, double mu, double su, double y, double beta, double a, double b, double Alpha_star, double s, double epsilon, double d_x, double x_m, double R, DxxParameters_structure DxxParameters), DxxParameters_structure DxxParameters)
void Calculate(GridElement &L, GridElement &epc, GridElement &alpha, DxxParameters_structure DxxParameters)
function calculate Dxx
int nint
Number of points in integral.
Definition: Parameters.h:59
string filename
File name for loading or saving diffusion coefficients array.
Definition: Parameters.h:36
std::vector< double > rrouts(double x_1, double x_2, double eta1, double eta2, double eta3, double epsilon, double beta, double mu, double Alpha_star, double a, DxxParameters_structure DxxParameters)
routs finding routine
static const double pi
π
Matrix3D times(const Matrix3D< T > &M) const
Arraywise multiplication (A.*B), stores result in a new matrix.
Definition: Matrix.cpp:1134
string numberDensity
Number density model. Check StrToVal(string input, NumberDensities &place) for known values...
Definition: Parameters.h:46
double mirror_point_coeff
Coefficient for mirror point (like 0.999, cause we can not integrate all the way to mirror point) ...
Definition: Parameters.h:60
double Df_Ozeke(double L, double Kp)
Matrix1D< double > Lpp
Lpp array.
Definition: Parameters.h:114
bool LoadDiffusionCoefficientFromFileWithGrid(GridElement &L, GridElement &pc, GridElement &alpha, string D_filename, string gridOrder="IFT_GRID_LPA")
Load Dxx from the file with Grid.
string DLLType
DLL type (which method to use to calculate). Check StrToVal(string input, DLLTypes &place) for known ...
Definition: Parameters.h:105
double B(double lambda, double L)
double int_Daa_loc(double lambda, DxxParameters_structure DxxParameters)
string DxxType
Dxx type. Check StrToVal(string input, DiffusionCoefficientTypes &place) for known values...
Definition: Parameters.h:29
double bisection(double(*f)(double x, double Alpha), double Alpha, double x0, double x1, double eps)
Looking for roots of function f for given Alpha on the interval x0-x1.
Definition: bisection.cpp:22
double f1(double lambda)
Main parameters structure that holds smaller structures for individual parameters.
Definition: Parameters.h:94
double Bw
Bw value.
Definition: Parameters.h:57
bool useEnergyDiffusion
Using diffusions flags.
Definition: Parameters.h:102
int size
size of grid element
Definition: Grid.h:39
double Df_Ozeke_E(double L, double Kp)
double int_Dpp_loc(double lambda, DxxParameters_structure DxxParameters)
Matrix3D< double > arr
array of diffusion coefficients
double Dpa_root(double Omega_e, double x, double mu, double su, double y, double beta, double a, double b, double Alpha_star, double s, double epsilon, double d_x, double x_m, double R, DxxParameters_structure DxxParameters)
double Dfe(double alpha, double Ke, double L, double Kp)
Electrostatic radial diffusion coefficient.
string waveName
Wave name.
Definition: Parameters.h:33
void CreateAllDiffusionCoefficients(DiffusionCoefficient &DLL, DiffusionCoefficientsGroup &Daa, DiffusionCoefficientsGroup &Dpcpc, DiffusionCoefficientsGroup &Dpca, DiffusionCoefficientsGroup &DaaLpp, DiffusionCoefficientsGroup &DpcpcLpp, DiffusionCoefficientsGroup &DpcaLpp, Parameters_structure &parameters, Grid &radialDiffusionGrid, Grid &localDiffusionsGrid)
void Output1DHeaders(ofstream &output1D, DiffusionCoefficientsGroup &Daa, DiffusionCoefficientsGroup &Dpcpc, DiffusionCoefficientsGroup &Dpca, DiffusionCoefficientsGroup &DaaLpp, DiffusionCoefficientsGroup &DpcpcLpp, DiffusionCoefficientsGroup &DpcaLpp)
Function used in main() to print header in output file. Header needs to understand the output file...
void get_quads(double *a, int n, double *quad, double *x)
Definition: rroots.cpp:201
void MakeDLL_B(GridElement &L, GridElement &pc, GridElement &alpha, double Kp)
Making DLL.
double multiplicator
Dxx is multiplied to this number after all other operations. Useful for fast Dxx modifications.
Definition: Parameters.h:49
void MakeDLL_Ozeke(GridElement &L, GridElement &pc, GridElement &alpha, double Kp)
Making DLL_M.
GridElement L
grid elements to be used
Definition: Grid.h:59
Parameters_structure parameters
Parameters structure, with all parameters from the parameters.ini file. The default parameters define...
Definition: Main.cpp:187
bool BwFromLambda
Bw lambda dependance flag.
Definition: Parameters.h:58
string type
Type of the diffusion coefficient: Daa, Dpp, Dpa etc... Described in types.h file as an enumeration...
GridElement epc
grid elements to be used
Definition: Grid.h:59
Error message - stack of single_errors.
Definition: error.h:54
void MakeDLL_BE(GridElement &L, GridElement &pc, GridElement &alpha, double Kp)
Making DLL.
double Dxx_ba(double L1, double EMeV, double Alpha, double int_Dxx_loc(double lambda, DxxParameters_structure DxxParameters), DxxParameters_structure DxxParameters)
double Omega_m
Omega_m value.
Definition: Parameters.h:52
Matrix1D< double > tau
Lifetime out of the plasmasphere.
Definition: Parameters.h:122
double lam_max
Maximum latitude.
Definition: Parameters.h:62
double omega_lc
omega lower cutoff
Definition: Parameters.h:55
string DxxName
Dxx name.
Definition: Parameters.h:30
Computational grid composed of 3 different GridElement.
Definition: Grid.h:53
double omega_uc
omega upper cutoff
Definition: Parameters.h:54
bool useRadialDiffusion
Using diffusions flags.
Definition: Parameters.h:100
static const double B_0
magnetic field at the equatorial plane on the Earth surface, Gauss
double F_cap2(double x, double y, double a, double beta, double mu, double s, double epsilon, double Alpha_star, DxxParameters_structure DxxParameters)
static const double qe
electron charge, CGS units
void MakeDLL(double Kp)
Making DLL.
double Kfunc(double pc)
Computation of Kinetic energy from given momentum .
double Scale(double Kp)
Function, that scale diffusion coefficients.
GridElement pc
grid elements to be used
Definition: Grid.h:59
vector< DxxParameters_structure > DxxParametersList
List of diffusion coefficients parameters.
Definition: Parameters.h:269
vector< DiffusionCoefficient > DxxList
List of diffusion coefficients in that group. Actually, it's a list of waves used in the diffusion co...