13 #include "../VariousFunctions/variousFunctions.h"
14 #include "../Exceptions/error.h"
20 #if defined(_WIN32) || defined(_WIN64)
21 #define strncasecmp _strnicmp
22 #define strcasecmp _stricmp
58 if (GridElement_parameters.size > 1)
60 if (GridElement_parameters.useLogScale)
61 arr[iteratorL][iteratorPc][iteratorAlpha] = pow(10, log10(GridElement_parameters.min) + gridElementDirection*(log10(GridElement_parameters.max) - log10(GridElement_parameters.min))/(GridElement_parameters.size - 1));
63 arr[iteratorL][iteratorPc][iteratorAlpha] = GridElement_parameters.min + gridElementDirection*(GridElement_parameters.max - GridElement_parameters.min)/(GridElement_parameters.size - 1);
65 arr[iteratorL][iteratorPc][iteratorAlpha] = GridElement_parameters.max;
83 string grid_filename,
string gridType,
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);
97 L.GridElement_parameters = parameters_L;
99 L.size = parameters_L.
size;
102 pc.GridElement_parameters = parameters_pc;
104 pc.size = parameters_pc.
size;
107 alpha.GridElement_parameters = parameters_alpha;
108 alpha.arr.
name = parameters_alpha.
name;
109 alpha.size = parameters_alpha.
size;
112 epc.GridElement_parameters = parameters_epc;
114 epc.size = parameters_epc.
size;
130 double mu, Jc, mu_check, Jc_check;
131 if (gridType ==
"GT_FILE") {
135 ifstream input(grid_filename.c_str());
136 if (input != NULL && !input.eof()) {
139 while (strcasecmp(inBuf.c_str(),
"ZONE") != 0 && !input.eof() ) {
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];
149 input.ignore(9999,
'\n');
151 throw error_msg(
"GRID_INITIALIZATION_ERROR",
"Not enough grid points in file %s - check grid resolution.", grid_filename.c_str());
156 input.ignore(1,
'\n');
158 Output::echo(0,
"Warning: probably too many points in the grid file %s - check grid resolution.", grid_filename.c_str());
164 for (il = 0; il < L.size; il++) {
165 for (im = 0; im < pc.size; im++) {
166 for (ia = 0; ia < alpha.size; ia++) {
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]);
176 throw error_msg(
"GRID_INITIALIZATION_ERROR",
"Grid file %s not found.", grid_filename.c_str());
178 }
else if (gridType ==
"GT_L_PC_ALPHA") {
180 for (il = 0; il < L.size ; il++){
181 for (im = 0; im < pc.size; im++){
182 for (ia = 0; ia < alpha.size; ia++){
184 L.SetRegularGridValue(il, im, ia, il);
185 pc.SetRegularGridValue(il, im, ia, im);
186 alpha.SetRegularGridValue(il, im, ia, ia);
194 for (il = 0; il < L.size; il++) {
195 for (im = 0; im < pc.size; im++) {
196 for (ia = 0; ia < alpha.size; ia++) {
200 Jacobian[il][im][ia] = pow(pc.arr[il][im][ia],2) *
VF::G(alpha.arr[il][im][ia]);
209 }
else if (gridType ==
"GT_L_MU_J") {
211 for (il = 0; il < L.size ; il++) {
212 for (im = 0; im < pc.size; im++) {
213 for (ia = 0; ia < alpha.size; ia++) {
215 L.SetRegularGridValue(il, im, ia, il);
221 for (im = 0; im < pc.size; im++){
222 for (ia = 0; ia < alpha.size; ia++) {
224 pc.SetRegularGridValue(il, im, ia, im);
225 alpha.SetRegularGridValue(il, im, ia, ia);
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);
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));
254 for (il = 0; il < L.size; il++) {
255 for (im = 0; im < pc.size; im++) {
256 for (ia = 0; ia < alpha.size; ia++) {
260 Jacobian[il][im][ia] = pow(L.arr[il][im][ia],-2);
270 }
else if (gridType ==
"GT_L_MU_ALPHA") {
272 for (il = 0; il < L.size ; il++) {
273 for (im = 0; im < pc.size; im++) {
274 for (ia = 0; ia < alpha.size; ia++) {
276 L.SetRegularGridValue(il, im, ia, il);
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);
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);
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]);
302 for (il = 0; il < L.size; il++) {
303 for (im = 0; im < pc.size; im++) {
304 for (ia = 0; ia < alpha.size; ia++) {
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]);
314 }
else if (gridType ==
"GT_PERP_TO_L_MU_J") {
316 throw error_msg(
"GRID_INIT_ERR",
"GT_PERP_TO_L_MU_J with empty second grid");
319 alpha = SecondGrid.
alpha;
322 double pc_local_max, pc_local_min;
323 for (il = 0; il < L.size ; il++){
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];
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++){
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;
343 for (il = 0; il < L.size; il++) {
344 for (im = 0; im < pc.size; im++) {
345 for (ia = 0; ia < alpha.size; ia++) {
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]);
355 throw error_msg(
"GRID_UNKNOWN_TYPE",
"Unknown grid type: %s", gridType.c_str());
358 Jacobian.writeToFile(
"./Debug/Jacobian.plt");
361 Grid_initialized =
true;
372 ofstream output(filename.c_str());
376 output <<
"VARIABLES = ";
377 output <<
"\"L, Earth radius\", \"Energy, MeV\", \"Pitch-angle, deg\", \"pc\" " << endl;
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;
401 double find_alpha(
double RHS,
double alpha_min,
double alpha_max,
double ERR,
int max_it,
int it) {
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;
410 if (yy * y0 < 0) alpha_max = alpha_root;
411 else if (yy * y0 > 0) alpha_min = alpha_root;
412 else return alpha_root;
414 }
while (fabs(alpha_max - alpha_min) > ERR);
GridElement alpha
grid elements to be used
Matrix3D< double > arr
array of grid points
Array of values of coordinate axes.
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.
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())
Grid element parameters structure.
int size
size of grid element
int size_y
size x, size y, size z
string name
Grid element name. Optional.
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)
GridElement L
grid elements to be used
double Y(double alpha)
Y(α) function.
double find_alpha(double RHS, double alpha_min, double alpha_max, double ERR, int max_it, int it)
Error message - stack of single_errors.
bool Grid_initialized
check for initialized status
void Output(string filename)
Computational grid composed of 3 different GridElement.
int size_x
size x, size y, size z
void Kfunc_equal(GridElement arg)
double Kfunc(double pc)
Computation of Kinetic energy from given momentum .
GridElement pc
grid elements to be used
int size_z
size x, size y, size z
double G(double alpha)
G-function.