VERB_code_2.3
Grid.cpp
Go to the documentation of this file.
1 
8 #include <iostream>
9 #include <iomanip>
10 
11 #include "Grid.h"
12 #include <math.h>
13 #include "../VariousFunctions/variousFunctions.h"
14 #include "../Exceptions/error.h"
15 
16 using namespace std;
17 
18 #define maxERR 1.e-6
19 
20 #if defined(_WIN32) || defined(_WIN64)
21 #define strncasecmp _strnicmp
22 #define strcasecmp _stricmp
23 #endif
24 
31  arr.AllocateMemory(pc.arr.size_x, pc.arr.size_y, pc.arr.size_z);
32  int x, y, z;
33  for (x = 0; x < pc.arr.size_x; x++)
34  for (y = 0; y < pc.arr.size_y; y++)
35  for (z = 0; z < pc.arr.size_z; z++)
36  arr[x][y][z] = VF::Kfunc(pc.arr[x][y][z]);
37 }
38 
57 void GridElement::SetRegularGridValue(int iteratorL, int iteratorPc, int iteratorAlpha, int gridElementDirection) {
58  if (GridElement_parameters.size > 1)
59  // in case of 3d
60  if (GridElement_parameters.useLogScale) // in case of log scale
61  arr[iteratorL][iteratorPc][iteratorAlpha] = pow(10, log10(GridElement_parameters.min) + gridElementDirection*(log10(GridElement_parameters.max) - log10(GridElement_parameters.min))/(GridElement_parameters.size - 1));
62  else // in case of regular scale
63  arr[iteratorL][iteratorPc][iteratorAlpha] = GridElement_parameters.min + gridElementDirection*(GridElement_parameters.max - GridElement_parameters.min)/(GridElement_parameters.size - 1);
64  else // in case of 2d
65  arr[iteratorL][iteratorPc][iteratorAlpha] = GridElement_parameters.max;
66 }
67 
81  Parameters_structure::GridElement parameters_alpha,
82  Parameters_structure::GridElement parameters_epc,
83  string grid_filename, string gridType,
84  Grid SecondGrid) {
85 
86  //if (!L.initialized) {
87  //throw error_msg("GRID_INITIALIZATION_ERROR", "Grid was not initialized");
88 
89  // allocating memory for all member classes
90  L.arr.AllocateMemory(parameters_L.size, parameters_pc.size, parameters_alpha.size);
91  pc.arr.AllocateMemory(parameters_L.size, parameters_pc.size, parameters_alpha.size);
92  alpha.arr.AllocateMemory(parameters_L.size, parameters_pc.size, parameters_alpha.size);
93  epc.arr.AllocateMemory(parameters_L.size, parameters_pc.size, parameters_alpha.size);
94  Jacobian.AllocateMemory(parameters_L.size, parameters_pc.size, parameters_alpha.size);
95 
96  //L.Initialize(parameters_L);
97  L.GridElement_parameters = parameters_L;
98  L.arr.name = parameters_L.name;
99  L.size = parameters_L.size;
100 
101  //pc.Initialize(parameters_pc);
102  pc.GridElement_parameters = parameters_pc;
103  pc.arr.name = parameters_pc.name;
104  pc.size = parameters_pc.size;
105 
106  //alpha.Initialize(parameters_alpha);
107  alpha.GridElement_parameters = parameters_alpha;
108  alpha.arr.name = parameters_alpha.name;
109  alpha.size = parameters_alpha.size;
110 
111  //epc.Initialize(parameters_epc);
112  epc.GridElement_parameters = parameters_epc;
113  epc.arr.name = parameters_epc.name;
114  epc.size = parameters_epc.size;
115 
116  type = gridType; // save the grid type for future use
117 
118  //LSize = L.size;
119  //pcSize = pc.size;
120  //alphaSize = alpha.size;
121  //epcSize = epc.size;
122 
123  //}
124 
125  // iterators, cold be useful later
126  int il, im, ia;
127  // temparary matrix_arrays
128  Matrix3D<double> mu_arr(L.size, pc.size, alpha.size);
129  // temporary variables
130  double mu, Jc, mu_check, Jc_check;
131  if (gridType == "GT_FILE") {
132  int il, im, ia;
133  string inBuf;
134  // making stresm
135  ifstream input(grid_filename.c_str());
136  if (input != NULL && !input.eof()) {
137  // Skipping header.
138  input >> inBuf;
139  while (strcasecmp(inBuf.c_str(), "ZONE") != 0 && !input.eof() ) {
140  input >> inBuf;
141  }
142  // read to the end of the line with 'zone'
143  input.ignore(9999, '\n');
144  for (il = 0; il < L.size; il++) {
145  for (im = 0; im < epc.size; im++) {
146  for (ia = 0; ia < alpha.size; ia++) {
147  input >> L.arr[il][im][ia] >> epc.arr[il][im][ia] >> alpha.arr[il][im][ia] >> pc.arr[il][im][ia];
148  // skip till the end of the line
149  input.ignore(9999, '\n');
150  if (input.eof()) {
151  throw error_msg("GRID_INITIALIZATION_ERROR", "Not enough grid points in file %s - check grid resolution.", grid_filename.c_str());
152  }
153  }
154  }
155  }
156  input.ignore(1, '\n');
157  if (!input.eof()) {
158  Output::echo(0, "Warning: probably too many points in the grid file %s - check grid resolution.", grid_filename.c_str());
159  }
160 
161  // calculating Jacobian for the grid
162  // \todo Calculate Jacobian from the grid data
163  int il, im, ia;
164  for (il = 0; il < L.size; il++) {
165  for (im = 0; im < pc.size; im++) {
166  for (ia = 0; ia < alpha.size; ia++) {
167  // remove all constant - they don't matter
168  Jacobian[il][im][ia] = pow(pc.arr[il][im][ia],2) * pow(L.arr[il][im][ia],2) * VF::G(alpha.arr[il][im][ia]);
169  // Jacobian[il][im][ia] = -4.0 * pow(pc[il][im][ia],2) * pow(L[il][im][ia],2) * VF::G(alpha[il][im][ia])/VF::B(L[il][im][ia]);
170  //psd.G[q0][q1][q2] = -4.0 * pow(psd.pc[q0][q1][q2],2) * pow(psd.L[q0][q1][q2],2) * G(psd.alpha[q0][q1][q2])/B(psd.L[q0][q1][q2]);
171  }
172  }
173  }
174 
175  } else {
176  throw error_msg("GRID_INITIALIZATION_ERROR", "Grid file %s not found.", grid_filename.c_str());
177  }
178  } else if (gridType == "GT_L_PC_ALPHA") {
179  // making grid
180  for (il = 0; il < L.size ; il++){
181  for (im = 0; im < pc.size; im++){
182  for (ia = 0; ia < alpha.size; ia++){
183  // You see, instead of that multiline formula for grid, appears everywherem we have just one.
184  L.SetRegularGridValue(il, im, ia, il);
185  pc.SetRegularGridValue(il, im, ia, im);
186  alpha.SetRegularGridValue(il, im, ia, ia);
187  }
188  }
189  }
190  //epc = pc.Kfunc();
191  epc.Kfunc_equal(pc);
192 
193  // calculating Jacobian for the grid
194  for (il = 0; il < L.size; il++) {
195  for (im = 0; im < pc.size; im++) {
196  for (ia = 0; ia < alpha.size; ia++) {
197  // remove all constant - they don't matter
198  // !!!!!!!!!!!!!Jacobian[il][im][ia] = pow(pc.arr[il][im][ia],2) * pow(L.arr[il][im][ia],2) * VF::G(alpha.arr[il][im][ia]);
199  // Jacobian[il][im][ia] = -4.0 * pow(pc[il][im][ia],2) * pow(L[il][im][ia],2) * VF::G(alpha[il][im][ia])/VF::B(L[il][im][ia]);
200  Jacobian[il][im][ia] = pow(pc.arr[il][im][ia],2) * VF::G(alpha.arr[il][im][ia]);
201 
202  }
203  }
204  }
205  //Jacobian = 1;
206  //Output::echo(0, "Jacobian == 1 !!!!!!!!!!!!!!!!\n");
207 
208 
209  } else if (gridType == "GT_L_MU_J") {
210  // Make L regular grid
211  for (il = 0; il < L.size ; il++) {
212  for (im = 0; im < pc.size; im++) {
213  for (ia = 0; ia < alpha.size; ia++) {
214  // You see, instead of that multiline formula for grid, appears everywherem we have just one.
215  L.SetRegularGridValue(il, im, ia, il);
216  }
217  }
218  }
219  // Making a 2-D grid on the top level at L=Lmax
220  il = L.size-1;
221  for (im = 0; im < pc.size; im++){
222  for (ia = 0; ia < alpha.size; ia++) {
223  // You see, instead of that multiline formula for grid, appears everywherem we have just one.
224  pc.SetRegularGridValue(il, im, ia, im);
225  alpha.SetRegularGridValue(il, im, ia, ia);
226  }
227  }
228  // Making a grid which conserves 1st J1 and second J2 adiabatic invariant
229  if (L.size > 1) {
230  for (im = 0; im < pc.size; im++) {
231  for (ia = 0; ia < alpha.size; ia++) {
232  mu = VF::mu_calc(L.arr[L.size-1][im][ia], pc.arr[L.size-1][im][ia], alpha.arr[L.size-1][im][ia]);
233  Jc = VF::Jc_calc(L.arr[L.size-1][im][ia], pc.arr[L.size-1][im][ia], alpha.arr[L.size-1][im][ia]);
234  for (il = 0; il < L.size-1 ; il++) {
235  if (alpha.arr[L.size-1][im][ia] < VC::pi/2)
236  alpha.arr[il][im][ia] = find_alpha(Jc * sqrt(L.arr[il][im][ia]) / sqrt(8.0 * VF::B(1) * VC::mc2 * mu), alpha.GridElement_parameters.min, VC::pi/2);
237  else if (alpha.arr[L.size-1][im][ia] > VC::pi/2)
238  alpha.arr[il][im][ia] = find_alpha(Jc * sqrt(L.arr[il][im][ia]) / sqrt(8.0 * VF::B(1) * VC::mc2 * mu), VC::pi/2, alpha.GridElement_parameters.max);
239  else
240  alpha.arr[il][im][ia] = VC::pi/2;
241  pc.arr[il][im][ia] = sqrt(2.0 * mu * VC::mc2 * VF::B(L.arr[il][im][ia])) / sin(alpha.arr[il][im][ia]);
242  mu_check = VF::mu_calc(L.arr[il][im][ia], pc.arr[il][im][ia], alpha.arr[il][im][ia]);
243  Jc_check = VF::Jc_calc(L.arr[il][im][ia], pc.arr[il][im][ia], alpha.arr[il][im][ia]);
244  if (fabs(mu - mu_check) > maxERR || fabs(Jc - Jc_check) > maxERR) {
245  throw error_msg("L_MU_J_GRID_ERROR", "Grid computation error (L, pc, a) = (%f, %f, %f) :\n mu=%e, Jc=%e\n mu_chk=%e, Jc_chk=%e\n mu_err=%e, Jc_err=%e.\n", L.arr[il][im][ia], pc.arr[il][im][ia], alpha.arr[il][im][ia], mu, Jc, mu_check, Jc_check, fabs(mu - mu_check), fabs(Jc - Jc_check));
246  }
247  }
248  }
249  }
250  }
251  epc.Kfunc_equal(pc);
252 
253  // calculating Jacobian for the grid
254  for (il = 0; il < L.size; il++) {
255  for (im = 0; im < pc.size; im++) {
256  for (ia = 0; ia < alpha.size; ia++) {
257  // remove all constant - they don't matter
258 
259 // Jacobian[il][im][ia] = pow(pc[il][im][ia],2) * pow(L[il][im][ia],-2) * VF::G(alpha[il][im][ia]);;
260  Jacobian[il][im][ia] = pow(L.arr[il][im][ia],-2); // <--
261 // Jacobian[il][im][ia] = 1;
262 // Jacobian[il][im][ia] = pow(L[il][im][ia],2);
263 
264  // Jacobian[il][im][ia] = -4.0 * pow(pc[il][im][ia],2) * pow(L[il][im][ia],2) * VF::G(alpha[il][im][ia])/VF::B(L[il][im][ia]);
265  //psd.G[q0][q1][q2] = -4.0 * pow(psd.pc[q0][q1][q2],2) * pow(psd.L[q0][q1][q2],2) * G(psd.alpha[q0][q1][q2])/B(psd.L[q0][q1][q2]);
266  }
267  }
268  }
269 
270  } else if (gridType == "GT_L_MU_ALPHA") {
271  // Make L regular grid
272  for (il = 0; il < L.size ; il++) {
273  for (im = 0; im < pc.size; im++) {
274  for (ia = 0; ia < alpha.size; ia++) {
275  // You see, instead of that multiline formula for grid, appears everywherem we have just one.
276  L.SetRegularGridValue(il, im, ia, il);
277  }
278  }
279  }
280  // Making a 2-D grid on the top level at L=Lmax
281  il = L.size-1;
282  for (im = 0; im < pc.size; im++){
283  for (ia = 0; ia < alpha.size; ia++) {
284  pc.SetRegularGridValue(il, im, ia, im);
285  alpha.SetRegularGridValue(il, im, ia, ia);
286 // mu_arr[il][im][ia] = VF::pc2mu(L[il][im][ia], pc[il][im][ia], alpha[il][im][ia]);
287  }
288  }
289  for (im = 0; im < pc.size; im++) {
290  for (ia = 0; ia < alpha.size; ia++) {
291  for (il = 0; il < L.size-1 ; il++) {
292  alpha.SetRegularGridValue(il, im, ia, ia);
293  // mu on the top
294  mu = VF::pc2mu(L.arr[L.size-1][im][ia], pc.arr[L.size-1][im][ia], alpha.arr[L.size-1][im][ia]);
295  pc.arr[il][im][ia] = VF::mu2pc(L.arr[il][im][ia], mu, alpha.arr[il][im][ia]);
296  }
297  }
298  }
299  epc.Kfunc_equal(pc);
300 
301  // calculating Jacobian for the grid
302  for (il = 0; il < L.size; il++) {
303  for (im = 0; im < pc.size; im++) {
304  for (ia = 0; ia < alpha.size; ia++) {
305  // remove all constant - they don't matter
306  // L^-2?? If it's on (L, pc, alpha) grid, I think it should be L^2. It does not matter cause it cancels anyway, but still.
307  Jacobian[il][im][ia] = pow(pc.arr[il][im][ia],2) * pow(L.arr[il][im][ia],-2) * VF::G(alpha.arr[il][im][ia]);
308  // Jacobian[il][im][ia] = -4.0 * pow(pc[il][im][ia],2) * pow(L[il][im][ia],2) * VF::G(alpha[il][im][ia])/VF::B(L[il][im][ia]);
309  //psd.G[q0][q1][q2] = -4.0 * pow(psd.pc[q0][q1][q2],2) * pow(psd.L[q0][q1][q2],2) * G(psd.alpha[q0][q1][q2])/B(psd.L[q0][q1][q2]);
310  }
311  }
312  }
313 
314  } else if (gridType == "GT_PERP_TO_L_MU_J") {
315  if (!SecondGrid.Grid_initialized) {
316  throw error_msg("GRID_INIT_ERR", "GT_PERP_TO_L_MU_J with empty second grid");
317  }
318  L = SecondGrid.L;
319  alpha = SecondGrid.alpha;
320 
321  // make new pc
322  double pc_local_max, pc_local_min;
323  for (il = 0; il < L.size ; il++){
324  // maximum from the upper ENERGY boundary
325  pc_local_max = SecondGrid.pc.arr[il][SecondGrid.pc.size-1][alpha.size-1];
326  for (ia = 0; ia < alpha.size; ia++) if (pc_local_max < SecondGrid.pc.arr[il][SecondGrid.pc.size-1][ia]) pc_local_max = SecondGrid.pc.arr[il][SecondGrid.pc.size-1][ia];
327  // MAXIMUM from the LOWER energy boundary
328  pc_local_min = SecondGrid.pc.arr[il][0][alpha.size-1];
329  for (ia = 0; ia < alpha.size; ia++) if (pc_local_min < SecondGrid.pc.arr[il][0][ia]) pc_local_min = SecondGrid.pc.arr[il][0][ia];
330  for (im = 0; im < pc.size; im++){
331  for (ia = 0; ia < alpha.size; ia++){
332  if (pc.size > 1) {
333  if (!pc.GridElement_parameters.useLogScale) pc.arr[il][im][ia] = pc_local_min + im*(pc_local_max - pc_local_min)/(pc.size - 1);
334  else pc.arr[il][im][ia] = pow(10, log10(pc_local_min) + im*(log10(pc_local_max) - log10(pc_local_min))/(pc.size - 1));
335  } else pc.arr[il][im][ia] = pc.GridElement_parameters.max;
336  }
337  }
338  }
339 
340  epc.Kfunc_equal(pc);
341 
342  // calculating Jacobian for the grid
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  // remove all constant - they don't matter
347  Jacobian[il][im][ia] = pow(pc.arr[il][im][ia],2) * pow(L.arr[il][im][ia],2) * VF::G(alpha.arr[il][im][ia]);
348  // Jacobian[il][im][ia] = -4.0 * pow(pc[il][im][ia],2) * pow(L[il][im][ia],2) * VF::G(alpha[il][im][ia])/VF::B(L[il][im][ia]);
349  //psd.G[q0][q1][q2] = -4.0 * pow(psd.pc[q0][q1][q2],2) * pow(psd.L[q0][q1][q2],2) * G(psd.alpha[q0][q1][q2])/B(psd.L[q0][q1][q2]);
350  }
351  }
352  }
353 
354  } else {
355  throw error_msg("GRID_UNKNOWN_TYPE", "Unknown grid type: %s", gridType.c_str());
356  }
357 
358  Jacobian.writeToFile("./Debug/Jacobian.plt");
359 
360 
361  Grid_initialized = true;
362 
363 }
364 
365 
371 void Grid::Output(string filename) {
372  ofstream output(filename.c_str());
373  int il, im, ia;
374 
375  // output header of the file
376  output << "VARIABLES = ";
377  output << "\"L, Earth radius\", \"Energy, MeV\", \"Pitch-angle, deg\", \"pc\" " << endl;
378 
379  // output values...
380  output << "ZONE T=\"Grid\" " << " I=" << alpha.size << ", J=" << epc.size << ", K=" << L.size << endl;
381  output.setf(ios_base::scientific, ios_base::floatfield);
382  for (il = 0; il < L.size; il++) {
383  for (im = 0; im < epc.size; im++) {
384  for (ia = 0; ia < alpha.size; ia++) {
385  output << L.arr[il][im][ia] << "\t" << epc.arr[il][im][ia] << "\t" << alpha.arr[il][im][ia] << "\t" << VF::pfunc(epc.arr[il][im][ia]) << endl;
386  }
387  }
388  }
389 }
390 
401 double find_alpha(double RHS, double alpha_min, double alpha_max, double ERR, int max_it, int it) {
402 
403  double alpha_root;
404  do {
405  it++;
406 
407  alpha_root = (alpha_min + alpha_max) / 2;
408  double y0 = VF::Y(alpha_min)/sin(alpha_min) - RHS, yy = VF::Y(alpha_root)/sin(alpha_root) - RHS;
409 
410  if (yy * y0 < 0) alpha_max = alpha_root;
411  else if (yy * y0 > 0) alpha_min = alpha_root;
412  else return alpha_root;
413 
414  } while (fabs(alpha_max - alpha_min) > ERR);
415 
416  return alpha_root;
417 
418 }
GridElement alpha
grid elements to be used
Definition: Grid.h:59
Matrix3D< double > arr
array of grid points
Definition: Grid.h:30
Array of values of coordinate axes.
Definition: Grid.h:28
double Jc_calc(double L, double pc, double Alpha)
Caltulating J*c.
double mu2pc(double L, double mu, double alpha)
Convert μ to pc.
void echo(int msglvl, const char *format,...)
Basic function! Write the message to the file.
Definition: Output.cpp:44
double pfunc(double K)
Computation of moumentum from Kinetic energy .
static const double mc2
m*c^2 = 0.511 MeV
void Create_Grid(Parameters_structure::GridElement parameters_L, Parameters_structure::GridElement parameters_pc, Parameters_structure::GridElement parameters_alpha, Parameters_structure::GridElement parameters_epc, string grid_filename, string gridType, Grid SecondGrid=Grid())
Definition: Grid.cpp:79
static const double pi
π
Grid element parameters structure.
Definition: Parameters.h:151
int size
size of grid element
Definition: Grid.h:39
int size_y
size x, size y, size z
Definition: Matrix.h:225
string name
Grid element name. Optional.
Definition: Parameters.h:153
double pc2mu(double L, double pc, double alpha)
Convert pc to μ.
double mu_calc(double L, double pc, double Alpha)
Caltulating μ.
double B(double Lparam)
Dipole magnetic field.
void SetRegularGridValue(int il, int im, int ia, int gridElementDirection)
Definition: Grid.cpp:57
GridElement L
grid elements to be used
Definition: Grid.h:59
double Y(double alpha)
Y(α) function.
double find_alpha(double RHS, double alpha_min, double alpha_max, double ERR, int max_it, int it)
Definition: Grid.cpp:401
Error message - stack of single_errors.
Definition: error.h:54
bool Grid_initialized
check for initialized status
Definition: Grid.h:67
void Output(string filename)
Definition: Grid.cpp:371
Computational grid composed of 3 different GridElement.
Definition: Grid.h:53
int size_x
size x, size y, size z
Definition: Matrix.h:225
void Kfunc_equal(GridElement arg)
Definition: Grid.cpp:30
double Kfunc(double pc)
Computation of Kinetic energy from given momentum .
GridElement pc
grid elements to be used
Definition: Grid.h:59
int size_z
size x, size y, size z
Definition: Matrix.h:225
double G(double alpha)
G-function.