00001
00008 #include <iostream>
00009 #include <iomanip>
00010
00011 #include "Grid.h"
00012 #include <math.h>
00013 #include "../VariousFunctions/variousFunctions.h"
00014 #include "../Exceptions/error.h"
00015
00016 using namespace std;
00017
00018 #define maxERR 1.e-6
00019
00020 #if defined(_WIN32) || defined(_WIN64)
00021 #define strncasecmp _strnicmp
00022 #define strcasecmp _stricmp
00023 #endif
00024
00035 GridElement::GridElement(ParamStructure::GridElement parameters, int L_size, int pc_size, int alpha_size) {
00036
00037 AllocateMemory(L_size, pc_size, alpha_size);
00038
00039 Initialize(parameters);
00040 }
00041
00048 void GridElement::Initialize(ParamStructure::GridElement parameters, int L_size, int pc_size, int alpha_size) {
00049
00050 AllocateMemory(L_size, pc_size, alpha_size);
00051
00052 Initialize(parameters);
00053 }
00054
00062 void GridElement::AllocateMemory(int L_size, int pc_size, int alpha_size) {
00063
00064 Matrix3D<double>::AllocateMemory(L_size, pc_size, alpha_size);
00065 }
00066
00072 void GridElement::Initialize(ParamStructure::GridElement parameters) {
00073 GridElement_parameters = parameters;
00074 this->name = parameters.name;
00075 this->size = parameters.size;
00076 }
00077
00083 GridElement& GridElement::operator= (const Matrix3D<double> &M) {
00084
00085 Matrix3D<double>::operator=(M);
00086 return *this;
00087 }
00088
00094 GridElement& GridElement::operator= (double val) {
00095
00096 Matrix3D<double>::operator=(val);
00097 return *this;
00098 }
00099
00103 GridElement GridElement::Kfunc() {
00104 int x, y, z;
00105 GridElement tempGridElement = *this;
00106 for (x = 0; x < this->size_x; x++)
00107 for (y = 0; y < this->size_y; y++)
00108 for (z = 0; z < this->size_z; z++)
00109 tempGridElement[x][y][z] = VF::Kfunc(((*this))[x][y][z]);
00110 return tempGridElement;
00111 }
00112
00120 void GridElement::Kfunc(GridElement arg) {
00121 (*this) = arg;
00122 int x, y, z;
00123 for (x = 0; x < this->size_x; x++)
00124 for (y = 0; y < this->size_y; y++)
00125 for (z = 0; z < this->size_z; z++)
00126 (*this)[x][y][z] = VF::Kfunc(arg[x][y][z]);
00127 }
00128
00147 void GridElement::SetRegularGridValue(int iteratorL, int iteratorPc, int iteratorAlpha, int gridElementDirection) {
00148 if (GridElement_parameters.size > 1)
00149
00150 if (GridElement_parameters.useLogScale)
00151 (*this)[iteratorL][iteratorPc][iteratorAlpha] = pow(10, log10(GridElement_parameters.min) + gridElementDirection*(log10(GridElement_parameters.max) - log10(GridElement_parameters.min))/(GridElement_parameters.size - 1));
00152 else
00153 (*this)[iteratorL][iteratorPc][iteratorAlpha] = GridElement_parameters.min + gridElementDirection*(GridElement_parameters.max - GridElement_parameters.min)/(GridElement_parameters.size - 1);
00154 else
00155 (*this)[iteratorL][iteratorPc][iteratorAlpha] = GridElement_parameters.max;
00156 }
00157
00171
00172 Grid::Grid(ParamStructure::GridElement parameters_L,
00173 ParamStructure::GridElement parameters_pc,
00174 ParamStructure::GridElement parameters_alpha,
00175 ParamStructure::GridElement parameters_epc,
00176 string grid_filename,
00177 string grid_type, Grid SecondGrid) {
00178
00179 Initialize(parameters_L,
00180 parameters_pc,
00181 parameters_alpha,
00182 parameters_epc,
00183 grid_type);
00184
00185 if (grid_type == "GT_FILE") {
00186 int il, im, ia;
00187 string inBuf;
00188
00189 ifstream input(grid_filename.c_str());
00190 if (input != NULL && !input.eof()) {
00191
00192 input >> inBuf;
00193 while (strcasecmp(inBuf.c_str(), "ZONE") != 0 && !input.eof() ) {
00194 input >> inBuf;
00195 }
00196
00197 input.ignore(9999, '\n');
00198 for (il = 0; il < L.size; il++) {
00199 for (im = 0; im < epc.size; im++) {
00200 for (ia = 0; ia < alpha.size; ia++) {
00201 input >> L[il][im][ia] >> epc[il][im][ia] >> alpha[il][im][ia] >> pc[il][im][ia];
00202
00203 input.ignore(9999, '\n');
00204 if (input.eof()) {
00205 throw error_msg("GRID_INITIALIZATION_ERROR", "Not enough grid points in file %s - check grid resolution.", grid_filename.c_str());
00206 }
00207 }
00208 }
00209 }
00210 input.ignore(1, '\n');
00211 if (!input.eof()) {
00212 Output::echo(0, "Warning: probably too many points in the grid file %s - check grid resolution.", grid_filename.c_str());
00213 }
00214
00215
00216
00217 Jacobian.AllocateMemory(L.size, pc.size, alpha.size);
00218 int il, im, ia;
00219 for (il = 0; il < L.size; il++) {
00220 for (im = 0; im < pc.size; im++) {
00221 for (ia = 0; ia < alpha.size; ia++) {
00222
00223 Jacobian[il][im][ia] = pow(pc[il][im][ia],2) * pow(L[il][im][ia],2) * VF::G(alpha[il][im][ia]);
00224
00225
00226 }
00227 }
00228 }
00229
00230 } else {
00231 throw error_msg("GRID_INITIALIZATION_ERROR", "Grid file %s was not found", grid_filename.c_str());
00232 }
00233 } else {
00234 MakeGrid(grid_type, SecondGrid);
00235 }
00236
00237 }
00238
00250 void Grid::Initialize(ParamStructure::GridElement parameters_L,
00251 ParamStructure::GridElement parameters_pc,
00252 ParamStructure::GridElement parameters_alpha,
00253 ParamStructure::GridElement parameters_epc,
00254 string grid_type) {
00255
00256 L.AllocateMemory(parameters_L.size, parameters_pc.size, parameters_alpha.size);
00257 pc.AllocateMemory(parameters_L.size, parameters_pc.size, parameters_alpha.size);
00258 alpha.AllocateMemory(parameters_L.size, parameters_pc.size, parameters_alpha.size);
00259 epc.AllocateMemory(parameters_L.size, parameters_pc.size, parameters_alpha.size);
00260 L.Initialize(parameters_L);
00261 pc.Initialize(parameters_pc);
00262 alpha.Initialize(parameters_alpha);
00263 epc.Initialize(parameters_epc);
00264 type = grid_type;
00265 LSize = L.size;
00266 pcSize = pc.size;
00267 alphaSize = alpha.size;
00268 epcSize = epc.size;
00269 initialized = true;
00270 }
00271
00280
00281 void Grid::MakeGrid(string gridType, Grid SecondGrid) {
00282 if (!initialized) {
00283 throw error_msg("GRID_INITIALIZATION_ERROR", "Grid was not initialized");
00284 }
00285
00286
00287 int il, im, ia;
00288
00289 Matrix3D<double> mu_arr(L.size, pc.size, alpha.size);
00290
00291 double mu, Jc, mu_check, Jc_check;
00292 if (gridType == "GT_L_PC_ALPHA") {
00293
00294 for (il = 0; il < L.size ; il++){
00295 for (im = 0; im < pc.size; im++){
00296 for (ia = 0; ia < alpha.size; ia++){
00297
00298 L.SetRegularGridValue(il, im, ia, il);
00299 pc.SetRegularGridValue(il, im, ia, im);
00300 alpha.SetRegularGridValue(il, im, ia, ia);
00301 }
00302 }
00303 }
00304
00305 epc.Kfunc(pc);
00306
00307
00308 Jacobian.AllocateMemory(L.size, pc.size, alpha.size);
00309 for (il = 0; il < L.size; il++) {
00310 for (im = 0; im < pc.size; im++) {
00311 for (ia = 0; ia < alpha.size; ia++) {
00312
00313 Jacobian[il][im][ia] = pow(pc[il][im][ia],2) * pow(L[il][im][ia],2) * VF::G(alpha[il][im][ia]);
00314
00315
00316 }
00317 }
00318 }
00319
00320 } else if (gridType == "GT_L_MU_J") {
00321
00322 for (il = 0; il < L.size ; il++) {
00323 for (im = 0; im < pc.size; im++) {
00324 for (ia = 0; ia < alpha.size; ia++) {
00325
00326 L.SetRegularGridValue(il, im, ia, il);
00327 }
00328 }
00329 }
00330
00331 il = L.size-1;
00332 for (im = 0; im < pc.size; im++){
00333 for (ia = 0; ia < alpha.size; ia++) {
00334
00335 pc.SetRegularGridValue(il, im, ia, im);
00336 alpha.SetRegularGridValue(il, im, ia, ia);
00337 }
00338 }
00339
00340 if (L.size > 1) {
00341 for (im = 0; im < pc.size; im++) {
00342 for (ia = 0; ia < alpha.size; ia++) {
00343 mu = VF::mu_calc(L[L.size-1][im][ia], pc[L.size-1][im][ia], alpha[L.size-1][im][ia]);
00344 Jc = VF::Jc_calc(L[L.size-1][im][ia], pc[L.size-1][im][ia], alpha[L.size-1][im][ia]);
00345 for (il = 0; il < L.size-1 ; il++) {
00346 if (alpha[L.size-1][im][ia] < VC::pi/2)
00347 alpha[il][im][ia] = find_alpha(Jc * sqrt(L[il][im][ia]) / sqrt(8.0 * VF::B(1) * VC::mc2 * mu), alpha.GridElement_parameters.min, VC::pi/2);
00348 else if (alpha[L.size-1][im][ia] > VC::pi/2)
00349 alpha[il][im][ia] = find_alpha(Jc * sqrt(L[il][im][ia]) / sqrt(8.0 * VF::B(1) * VC::mc2 * mu), VC::pi/2, alpha.GridElement_parameters.max);
00350 else
00351 alpha[il][im][ia] = VC::pi/2;
00352 pc[il][im][ia] = sqrt(2.0 * mu * VC::mc2 * VF::B(L[il][im][ia])) / sin(alpha[il][im][ia]);
00353 mu_check = VF::mu_calc(L[il][im][ia], pc[il][im][ia], alpha[il][im][ia]);
00354 Jc_check = VF::Jc_calc(L[il][im][ia], pc[il][im][ia], alpha[il][im][ia]);
00355 if (fabs(mu - mu_check) > maxERR || fabs(Jc - Jc_check) > maxERR) {
00356 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[il][im][ia], pc[il][im][ia], alpha[il][im][ia], mu, Jc, mu_check, Jc_check, fabs(mu - mu_check), fabs(Jc - Jc_check));
00357 }
00358 }
00359 }
00360 }
00361 }
00362 epc.Kfunc(pc);
00363
00364
00365 Jacobian.AllocateMemory(L.size, pc.size, alpha.size);
00366 for (il = 0; il < L.size; il++) {
00367 for (im = 0; im < pc.size; im++) {
00368 for (ia = 0; ia < alpha.size; ia++) {
00369
00370
00371
00372 Jacobian[il][im][ia] = pow(L[il][im][ia],-2);
00373
00374
00375
00376
00377
00378 }
00379 }
00380 }
00381
00382 } else if (gridType == "GT_L_MU_ALPHA") {
00383
00384 for (il = 0; il < L.size ; il++) {
00385 for (im = 0; im < pc.size; im++) {
00386 for (ia = 0; ia < alpha.size; ia++) {
00387
00388 L.SetRegularGridValue(il, im, ia, il);
00389 }
00390 }
00391 }
00392
00393 il = L.size-1;
00394 for (im = 0; im < pc.size; im++){
00395 for (ia = 0; ia < alpha.size; ia++) {
00396 pc.SetRegularGridValue(il, im, ia, im);
00397 alpha.SetRegularGridValue(il, im, ia, ia);
00398
00399 }
00400 }
00401 for (im = 0; im < pc.size; im++) {
00402 for (ia = 0; ia < alpha.size; ia++) {
00403 for (il = 0; il < L.size-1 ; il++) {
00404 alpha.SetRegularGridValue(il, im, ia, ia);
00405
00406 mu = VF::pc2mu(L[L.size-1][im][ia], pc[L.size-1][im][ia], alpha[L.size-1][im][ia]);
00407 pc[il][im][ia] = VF::mu2pc(L[il][im][ia], mu, alpha[il][im][ia]);
00408 }
00409 }
00410 }
00411 epc.Kfunc(pc);
00412
00413
00414 Jacobian.AllocateMemory(L.size, pc.size, alpha.size);
00415 for (il = 0; il < L.size; il++) {
00416 for (im = 0; im < pc.size; im++) {
00417 for (ia = 0; ia < alpha.size; ia++) {
00418
00419
00420 Jacobian[il][im][ia] = pow(pc[il][im][ia],2) * pow(L[il][im][ia],-2) * VF::G(alpha[il][im][ia]);
00421
00422
00423 }
00424 }
00425 }
00426
00427 } else if (gridType == "GT_PERP_TO_L_MU_J") {
00428 if (SecondGrid.initialized == false) {
00429 throw error_msg("GRID_INIT_ERR", "GT_PERP_TO_L_MU_J with empty second grid");
00430 }
00431 L = SecondGrid.L;
00432 alpha = SecondGrid.alpha;
00433
00434
00435 double pc_local_max, pc_local_min;
00436 for (il = 0; il < L.size ; il++){
00437
00438 pc_local_max = SecondGrid.pc[il][SecondGrid.pc.size-1][alpha.size-1];
00439 for (ia = 0; ia < alpha.size; ia++) if (pc_local_max < SecondGrid.pc[il][SecondGrid.pc.size-1][ia]) pc_local_max = SecondGrid.pc[il][SecondGrid.pc.size-1][ia];
00440
00441 pc_local_min = SecondGrid.pc[il][0][alpha.size-1];
00442 for (ia = 0; ia < alpha.size; ia++) if (pc_local_min < SecondGrid.pc[il][0][ia]) pc_local_min = SecondGrid.pc[il][0][ia];
00443 for (im = 0; im < pc.size; im++){
00444 for (ia = 0; ia < alpha.size; ia++){
00445 if (pc.size > 1) {
00446 if (!pc.GridElement_parameters.useLogScale) pc[il][im][ia] = pc_local_min + im*(pc_local_max - pc_local_min)/(pc.size - 1);
00447 else pc[il][im][ia] = pow(10, log10(pc_local_min) + im*(log10(pc_local_max) - log10(pc_local_min))/(pc.size - 1));
00448 } else pc[il][im][ia] = pc.GridElement_parameters.max;
00449 }
00450 }
00451 }
00452
00453 epc.Kfunc(pc);
00454
00455
00456 Jacobian.AllocateMemory(L.size, pc.size, alpha.size);
00457 for (il = 0; il < L.size; il++) {
00458 for (im = 0; im < pc.size; im++) {
00459 for (ia = 0; ia < alpha.size; ia++) {
00460
00461 Jacobian[il][im][ia] = pow(pc[il][im][ia],2) * pow(L[il][im][ia],2) * VF::G(alpha[il][im][ia]);
00462
00463
00464 }
00465 }
00466 }
00467
00468 } else {
00469 throw error_msg("GRID_UNKNOWN_TYPE", "Unknown grid type");
00470 }
00471
00472 }
00473
00479 void Grid::Output(string filename) {
00480 ofstream output(filename.c_str());
00481 int il, im, ia;
00482
00483
00484 output << "VARIABLES = ";
00485 output << "\"L, Earth radius\", \"Energy, MeV\", \"Pitch-angle, deg\", \"pc\" " << endl;
00486
00487
00488 output << "ZONE T=\"Grid\" " << " I=" << alpha.size << ", J=" << epc.size << ", K=" << L.size << endl;
00489 output.setf(ios_base::scientific, ios_base::floatfield);
00490 for (il = 0; il < L.size; il++) {
00491 for (im = 0; im < epc.size; im++) {
00492 for (ia = 0; ia < alpha.size; ia++) {
00493 output << L[il][im][ia] << "\t" << epc[il][im][ia] << "\t" << alpha[il][im][ia] << "\t" << VF::pfunc(epc[il][im][ia]) << endl;
00494 }
00495 }
00496 }
00497 }
00498
00509 double find_alpha(double RHS, double alpha_min, double alpha_max, double ERR, int max_it, int it) {
00510
00511 double alpha_root;
00512 do {
00513 it++;
00514
00515 alpha_root = (alpha_min + alpha_max) / 2;
00516 double y0 = VF::Y(alpha_min)/sin(alpha_min) - RHS, yy = VF::Y(alpha_root)/sin(alpha_root) - RHS;
00517
00518 if (yy * y0 < 0) alpha_max = alpha_root;
00519 else if (yy * y0 > 0) alpha_min = alpha_root;
00520 else return alpha_root;
00521
00522 } while (fabs(alpha_max - alpha_min) > ERR);
00523
00524 return alpha_root;
00525
00526 }
00527