13 #include "../VariousFunctions/variousFunctions.h"
14 #include "../Exceptions/error.h"
20 #if defined(_WIN32) || defined(_WIN64)
21 #define strncasecmp _strnicmp
22 #define strcasecmp _stricmp
60 if (GridElement_parameters.size > 1)
62 if (GridElement_parameters.useLogScale)
63 arr[iteratorL][iteratorPc][iteratorAlpha] = pow(10, log10(GridElement_parameters.min) + gridElementDirection*(log10(GridElement_parameters.max) - log10(GridElement_parameters.min))/(GridElement_parameters.size - 1));
65 arr[iteratorL][iteratorPc][iteratorAlpha] = GridElement_parameters.min + gridElementDirection*(GridElement_parameters.max - GridElement_parameters.min)/(GridElement_parameters.size - 1);
67 arr[iteratorL][iteratorPc][iteratorAlpha] = GridElement_parameters.max;
75 string grid_filename,
string gridType,
82 L.arr.AllocateMemory(parameters_L.
size, parameters_pc.
size, parameters_alpha.
size);
83 pc.arr.AllocateMemory(parameters_L.
size, parameters_pc.
size, parameters_alpha.
size);
84 alpha.arr.AllocateMemory(parameters_L.
size, parameters_pc.
size, parameters_alpha.
size);
85 epc.arr.AllocateMemory(parameters_L.
size, parameters_pc.
size, parameters_alpha.
size);
86 Jacobian.AllocateMemory(parameters_L.
size, parameters_pc.
size, parameters_alpha.
size);
89 L.GridElement_parameters = parameters_L;
91 L.size = parameters_L.
size;
94 pc.GridElement_parameters = parameters_pc;
96 pc.size = parameters_pc.
size;
99 alpha.GridElement_parameters = parameters_alpha;
100 alpha.arr.
name = parameters_alpha.
name;
101 alpha.size = parameters_alpha.
size;
104 epc.GridElement_parameters = parameters_epc;
106 epc.size = parameters_epc.
size;
122 double mu, Jc, mu_check, Jc_check;
123 if (gridType ==
"GT_FILE") {
127 ifstream input(grid_filename.c_str());
128 if (input != NULL && !input.eof()) {
131 while (strcasecmp(inBuf.c_str(),
"ZONE") != 0 && !input.eof() ) {
135 input.ignore(9999,
'\n');
136 for (il = 0; il < L.size; il++) {
137 for (im = 0; im < epc.size; im++) {
138 for (ia = 0; ia < alpha.size; ia++) {
139 input >> L.arr[il][im][ia] >> epc.arr[il][im][ia] >> alpha.arr[il][im][ia] >> pc.arr[il][im][ia];
141 input.ignore(9999,
'\n');
143 throw error_msg(
"GRID_INITIALIZATION_ERROR",
"Not enough grid points in file %s - check grid resolution.", grid_filename.c_str());
148 input.ignore(1,
'\n');
150 Output::echo(0,
"Warning: probably too many points in the grid file %s - check grid resolution.", grid_filename.c_str());
156 for (il = 0; il < L.size; il++) {
157 for (im = 0; im < pc.size; im++) {
158 for (ia = 0; ia < alpha.size; ia++) {
160 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]);
168 throw error_msg(
"GRID_INITIALIZATION_ERROR",
"Grid file %s not found.", grid_filename.c_str());
170 }
else if (gridType ==
"GT_L_PC_ALPHA") {
172 for (il = 0; il < L.size ; il++){
173 for (im = 0; im < pc.size; im++){
174 for (ia = 0; ia < alpha.size; ia++){
176 L.SetRegularGridValue(il, im, ia, il);
177 pc.SetRegularGridValue(il, im, ia, im);
178 alpha.SetRegularGridValue(il, im, ia, ia);
186 for (il = 0; il < L.size; il++) {
187 for (im = 0; im < pc.size; im++) {
188 for (ia = 0; ia < alpha.size; ia++) {
192 Jacobian[il][im][ia] = pow(pc.arr[il][im][ia],2) *
VF::G(alpha.arr[il][im][ia]);
201 }
else if (gridType ==
"GT_L_MU_J") {
203 for (il = 0; il < L.size ; il++) {
204 for (im = 0; im < pc.size; im++) {
205 for (ia = 0; ia < alpha.size; ia++) {
207 L.SetRegularGridValue(il, im, ia, il);
213 for (im = 0; im < pc.size; im++){
214 for (ia = 0; ia < alpha.size; ia++) {
216 pc.SetRegularGridValue(il, im, ia, im);
217 alpha.SetRegularGridValue(il, im, ia, ia);
222 for (im = 0; im < pc.size; im++) {
223 for (ia = 0; ia < alpha.size; ia++) {
224 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]);
225 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]);
226 for (il = 0; il < L.size-1 ; il++) {
227 if (alpha.arr[L.size-1][im][ia] <
VC::pi/2)
228 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);
229 else if (alpha.arr[L.size-1][im][ia] >
VC::pi/2)
230 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);
232 alpha.arr[il][im][ia] =
VC::pi/2;
233 pc.arr[il][im][ia] = sqrt(2.0 * mu *
VC::mc2 *
VF::B(L.arr[il][im][ia])) / sin(alpha.arr[il][im][ia]);
234 mu_check =
VF::mu_calc(L.arr[il][im][ia], pc.arr[il][im][ia], alpha.arr[il][im][ia]);
235 Jc_check =
VF::Jc_calc(L.arr[il][im][ia], pc.arr[il][im][ia], alpha.arr[il][im][ia]);
236 if (fabs(mu - mu_check) >
maxERR || fabs(Jc - Jc_check) >
maxERR) {
237 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 for (il = 0; il < L.size; il++) {
247 for (im = 0; im < pc.size; im++) {
248 for (ia = 0; ia < alpha.size; ia++) {
252 Jacobian[il][im][ia] = pow(L.arr[il][im][ia],-2);
262 }
else if (gridType ==
"GT_L_MU_ALPHA") {
264 for (il = 0; il < L.size ; il++) {
265 for (im = 0; im < pc.size; im++) {
266 for (ia = 0; ia < alpha.size; ia++) {
268 L.SetRegularGridValue(il, im, ia, il);
274 for (im = 0; im < pc.size; im++){
275 for (ia = 0; ia < alpha.size; ia++) {
276 pc.SetRegularGridValue(il, im, ia, im);
277 alpha.SetRegularGridValue(il, im, ia, ia);
281 for (im = 0; im < pc.size; im++) {
282 for (ia = 0; ia < alpha.size; ia++) {
283 for (il = 0; il < L.size-1 ; il++) {
284 alpha.SetRegularGridValue(il, im, ia, ia);
286 mu =
VF::pc2mu(L.arr[L.size-1][im][ia], pc.arr[L.size-1][im][ia], alpha.arr[L.size-1][im][ia]);
287 pc.arr[il][im][ia] =
VF::mu2pc(L.arr[il][im][ia], mu, alpha.arr[il][im][ia]);
294 for (il = 0; il < L.size; il++) {
295 for (im = 0; im < pc.size; im++) {
296 for (ia = 0; ia < alpha.size; ia++) {
299 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]);
306 }
else if (gridType ==
"GT_PERP_TO_L_MU_J") {
308 throw error_msg(
"GRID_INIT_ERR",
"GT_PERP_TO_L_MU_J with empty second grid");
311 alpha = SecondGrid.
alpha;
314 double pc_local_max, pc_local_min;
315 for (il = 0; il < L.size ; il++){
317 pc_local_max = SecondGrid.
pc.
arr[il][SecondGrid.
pc.
size-1][alpha.size-1];
318 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];
320 pc_local_min = SecondGrid.
pc.
arr[il][0][alpha.size-1];
321 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];
322 for (im = 0; im < pc.size; im++){
323 for (ia = 0; ia < alpha.size; ia++){
325 if (!pc.GridElement_parameters.useLogScale) pc.arr[il][im][ia] = pc_local_min + im*(pc_local_max - pc_local_min)/(pc.size - 1);
326 else pc.arr[il][im][ia] = pow(10, log10(pc_local_min) + im*(log10(pc_local_max) - log10(pc_local_min))/(pc.size - 1));
327 }
else pc.arr[il][im][ia] = pc.GridElement_parameters.max;
335 for (il = 0; il < L.size; il++) {
336 for (im = 0; im < pc.size; im++) {
337 for (ia = 0; ia < alpha.size; ia++) {
339 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]);
347 throw error_msg(
"GRID_UNKNOWN_TYPE",
"Unknown grid type: %s", gridType.c_str());
350 Jacobian.writeToFile(
"./Debug/Jacobian.plt");
353 Grid_initialized =
true;
364 ofstream output(filename.c_str());
368 output <<
"VARIABLES = ";
369 output <<
"\"L, Earth radius\", \"Energy, MeV\", \"Pitch-angle, deg\", \"pc\" " << endl;
372 output <<
"ZONE T=\"Grid\" " <<
" I=" << alpha.size <<
", J=" << epc.size <<
", K=" << L.size << endl;
373 output.setf(ios_base::scientific, ios_base::floatfield);
374 for (il = 0; il < L.size; il++) {
375 for (im = 0; im < epc.size; im++) {
376 for (ia = 0; ia < alpha.size; ia++) {
377 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;
396 double find_alpha(
double RHS,
double alpha_min,
double alpha_max,
double ERR,
int max_it,
int it) {
402 alpha_root = (alpha_min + alpha_max) / 2;
403 double y0 =
VF::Y(alpha_min)/sin(alpha_min) - RHS, yy =
VF::Y(alpha_root)/sin(alpha_root) - RHS;
405 if (yy * y0 < 0) alpha_max = alpha_root;
406 else if (yy * y0 > 0) alpha_min = alpha_root;
407 else return alpha_root;
409 }
while (fabs(alpha_max - alpha_min) > ERR);