File indexing completed on 2024-04-06 12:22:33
0001 #include "prepareMagneticFieldGrid.h"
0002
0003 #include "MagneticField/Interpolation/src/VectorFieldInterpolation.h"
0004 #include "MagneticField/Interpolation/src/binary_ofstream.h"
0005 #include "DataFormats/GeometryVector/interface/Pi.h"
0006 #include <iostream>
0007 #include <iomanip>
0008 #include <sstream>
0009 #include <iostream>
0010
0011 using namespace std;
0012
0013 void prepareMagneticFieldGrid::countTrueNumberOfPoints(const std::string &name) const {
0014 int nLines = 0;
0015 int nPoints = 0;
0016
0017
0018 IndexedDoubleVector XBVector;
0019 std::vector<IndexedDoubleVector> XBValues;
0020
0021
0022 double epsilonRadius = 1e-8;
0023 std::ifstream file(name.c_str());
0024 string line;
0025 if (file.good()) {
0026 while (getline(file, line)) {
0027 double x1, x2, x3, Bx, By, Bz, perm;
0028 stringstream linestr;
0029 linestr << line;
0030 linestr >> x1 >> x2 >> x3 >> Bx >> By >> Bz >> perm;
0031 ;
0032 if (file) {
0033 XBVector.putV6(x1, x2, x3, Bx, By, Bz);
0034
0035
0036
0037
0038
0039 XBValues.push_back(XBVector);
0040 ++nLines;
0041 }
0042 }
0043 }
0044
0045
0046 for (int iLine = 0; iLine < nLines; ++iLine) {
0047 double x1, x2, x3, Bx, By, Bz;
0048 XBVector = XBValues.operator[](iLine);
0049 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0050 double pnt[3] = {x1, x2, x3};
0051 double tmpBz = Bz;
0052 bool isSinglePoint = true;
0053 for (int i = 0; i < nLines; ++i) {
0054 if (i < iLine) {
0055 XBVector = XBValues.operator[](i);
0056 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0057 double distPP =
0058 sqrt((pnt[0] - x1) * (pnt[0] - x1) + (pnt[1] - x2) * (pnt[1] - x2) + (pnt[2] - x3) * (pnt[2] - x3));
0059 if (distPP < epsilonRadius) {
0060 isSinglePoint = false;
0061 cout << "same points(" << iLine << ") with dR = " << distPP << endl;
0062 cout << " point 1: " << iLine << " " << pnt[0] << " " << pnt[1] << " " << pnt[2]
0063 << " Bz: " << tmpBz << endl;
0064 cout << " point 2: " << i << " " << x1 << " " << x2 << " " << x3 << " Bz: " << Bz << endl;
0065 }
0066 }
0067 }
0068 if (isSinglePoint)
0069 ++nPoints;
0070 }
0071
0072 if (PRINT)
0073 cout << " file " << name << endl;
0074 if (PRINT)
0075 cout << " # lines = " << nLines << " # points = " << nPoints;
0076 if (nLines == nPoints) {
0077 if (PRINT)
0078 cout << " --> PASSED" << endl;
0079 } else {
0080 if (PRINT)
0081 cout << " --> FAILED" << endl;
0082 }
0083
0084 return;
0085 }
0086
0087 void prepareMagneticFieldGrid::fillFromFile(const std::string &name) {
0088 double phi = (rotateSector - 1) * Geom::pi() / 6.;
0089 double sphi = sin(-phi);
0090 double cphi = cos(-phi);
0091
0092
0093 std::string::size_type ibeg, iend;
0094 ibeg = name.find('-');
0095 iend = name.rfind('-');
0096 if ((name.substr(ibeg + 1, iend - ibeg - 1)) == "xyz") {
0097 XyzCoordinates = true;
0098 } else if ((name.substr(ibeg + 1, iend - ibeg - 1)) == "rpz") {
0099 RpzCoordinates = true;
0100 } else {
0101 cout << " Unrecognized input file name: valid formats are *-xyz-* or *-rpz-*" << endl;
0102 abort();
0103 }
0104
0105
0106 IndexedDoubleVector XBVector;
0107 std::vector<IndexedDoubleVector> XBValues;
0108 std::vector<IndexedDoubleVector> IXBValues;
0109 std::vector<IndexedDoubleVector> XBArray;
0110
0111
0112 std::vector<double> valX[3];
0113 std::vector<int> numX[3];
0114
0115
0116 int nLines = 0;
0117 int index[3] = {0, 0, 0};
0118 int nSteps[3] = {0, 0, 0};
0119 double x1, x2, x3, Bx, By, Bz, perm, poten;
0120 std::ifstream file(name.c_str());
0121 int n_air = 0;
0122 int n_iron = 0;
0123
0124 if (file.good()) {
0125 while (file.good()) {
0126 file >> x1 >> x2 >> x3 >> Bx >> By >> Bz >> perm >> poten;
0127 if (file) {
0128 if (rotateSector > 1) {
0129 if (RpzCoordinates) {
0130 x2 = x2 - phi;
0131 if (fabs(x2) > 0.78) {
0132 cout << "ERROR: Input coordinates do not belong to sector " << rotateSector << " : " << fabs(x2) << endl;
0133 abort();
0134 }
0135 } else {
0136 double x = x1;
0137 double y = x2;
0138 x1 = x * cphi - y * sphi;
0139 x2 = x * sphi + y * cphi;
0140 if (fabs(atan2(x2, x1)) > 0.78) {
0141 cout << "ERROR: Input coordinates do not belong to sector " << rotateSector << " : " << atan2(x2, x1)
0142 << endl;
0143 abort();
0144 }
0145 }
0146 double Bx0 = Bx;
0147 double By0 = By;
0148 Bx = Bx0 * cphi - By0 * sphi;
0149 By = Bx0 * sphi + By0 * cphi;
0150 }
0151
0152
0153
0154
0155 if (perm > 1.)
0156 ++n_iron;
0157 else
0158 ++n_air;
0159
0160 XBVector.putI3(index[0], index[1], index[2]);
0161 XBVector.putV6(x1, x2, x3, Bx, By, Bz);
0162 XBValues.push_back(XBVector);
0163 ++nLines;
0164
0165 double pnt[3] = {x1, x2, x3};
0166 if (nLines == 1) {
0167 for (int i = 0; i < 3; ++i) {
0168 valX[i].push_back(pnt[i]);
0169 numX[i].push_back(1);
0170 }
0171 } else {
0172 for (int i = 0; i < 3; ++i) {
0173 bool knownValue = false;
0174 for (int j = 0; j < (nSteps[i] + 1); ++j) {
0175 if (std::abs(valX[i].operator[](j) - pnt[i]) < EPSILON)
0176 knownValue = true;
0177 if (std::abs(valX[i].operator[](j) - pnt[i]) < EPSILON)
0178 numX[i].operator[](j) += 1;
0179 }
0180 if (!knownValue) {
0181 valX[i].push_back(pnt[i]);
0182 nSteps[i]++;
0183 numX[i].push_back(1);
0184 }
0185 }
0186 }
0187 }
0188 }
0189 if (n_air && n_iron)
0190 cout << "WARNING: table " << name << " contains " << n_air << " points with perm=1 and " << n_iron
0191 << " with perm>1" << endl;
0192 } else
0193 return;
0194
0195
0196
0197
0198 for (int i = 0; i < 3; ++i) {
0199
0200 sort(valX[i].begin(), valX[i].end());
0201
0202 EasyCoordinate[i] = true;
0203 NumberOfPoints[i] = nSteps[i] + 1;
0204 BasicDistance0[i] = 9e9;
0205 double newVal = 0.0;
0206 double oldVal = 0.0;
0207 for (int j = 0; j < NumberOfPoints[i]; ++j) {
0208 newVal = valX[i].operator[](j);
0209 if (j == 1)
0210 BasicDistance0[i] = newVal - oldVal;
0211 else if (j > 1) {
0212 if (std::abs(newVal - oldVal - BasicDistance0[i]) > EPSILON)
0213 EasyCoordinate[i] = false;
0214 }
0215 oldVal = newVal;
0216 }
0217 }
0218
0219
0220 for (int iLine = 0; iLine < nLines; ++iLine) {
0221 XBVector = XBValues.operator[](iLine);
0222 XBVector.getI3(index[0], index[1], index[2]);
0223 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0224 double pnt[3] = {x1, x2, x3};
0225 for (int i = 0; i < 3; ++i) {
0226 for (int j = 0; j < NumberOfPoints[i]; ++j) {
0227 if (EasyCoordinate[i]) {
0228 if (std::abs(valX[i].operator[](j) - pnt[i]) < EPSILON)
0229 index[i] = j;
0230 }
0231 }
0232 }
0233 XBVector.putI3(index[0], index[1], index[2]);
0234 IXBValues.push_back(XBVector);
0235 }
0236 if (!XBValues.empty())
0237 XBValues.clear();
0238
0239
0240 sort(IXBValues.begin(), IXBValues.end());
0241
0242
0243 bool systematicGrid[3] = {false, false, false};
0244 KnownStructure = true;
0245 if (NumberOfPoints[0] * NumberOfPoints[1] * NumberOfPoints[2] != nLines)
0246 KnownStructure = false;
0247 for (int i = 0; i < 3; ++i) {
0248 systematicGrid[i] = EasyCoordinate[i];
0249 if (!systematicGrid[i])
0250 KnownStructure = false;
0251 }
0252
0253 if (KnownStructure == true) {
0254 for (int iLine = 0; iLine < nLines; ++iLine) {
0255 XBVector = IXBValues.operator[](iLine);
0256 XBArray.push_back(XBVector);
0257 }
0258 }
0259
0260
0261 if (!XBArray.empty())
0262 XBVector = XBArray.front();
0263 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0264 double firstListPoint[3] = {x1, x2, x3};
0265 if (!XBArray.empty())
0266 XBVector = XBArray.back();
0267 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0268 double secondRefPoint[3] = {x1, x2, x3};
0269
0270
0271 KnownStructure = true;
0272 for (int i = 0; i < 3; ++i) {
0273 ReferencePoint[i] = firstListPoint[i];
0274
0275 double stepSize = BasicDistance0[i];
0276 for (int j = 0; j < 3; ++j) {
0277 stepSize += BasicDistance1[i][j] * nSteps[j];
0278 }
0279 double offset = 0.0;
0280 for (int j = 0; j < 3; ++j) {
0281 offset += BasicDistance2[i][j] * nSteps[j];
0282 }
0283
0284 double totDiff = secondRefPoint[i] - (ReferencePoint[i] + offset + stepSize * nSteps[i]);
0285 if (std::abs(totDiff) < EPSILON * 10.)
0286 systematicGrid[i] = true;
0287 else
0288 systematicGrid[i] = false;
0289 if (!systematicGrid[i])
0290 KnownStructure = false;
0291 }
0292
0293 if (KnownStructure) {
0294 if (XyzCoordinates)
0295 GridType = 1;
0296 if (RpzCoordinates)
0297 GridType = 3;
0298 if (PRINT) {
0299
0300 cout << " read " << nLines << " lines -> grid structure:" << endl;
0301 cout << " # of points: N1 = " << NumberOfPoints[0] << " N2 = " << NumberOfPoints[1]
0302 << " N3 = " << NumberOfPoints[2] << endl;
0303 cout << " ref. point : X1 = " << ReferencePoint[0] << " X2 = " << ReferencePoint[1]
0304 << " X3 = " << ReferencePoint[2] << endl;
0305 cout << " step parm0 : A1 = " << BasicDistance0[0] << " A2 = " << BasicDistance0[1]
0306 << " A3 = " << BasicDistance0[2] << endl;
0307 cout << " easy grid : E1 = " << EasyCoordinate[0] << " E2 = " << EasyCoordinate[1]
0308 << " E3 = " << EasyCoordinate[2] << endl;
0309 cout << " structure : S1 = " << systematicGrid[0] << " S2 = " << systematicGrid[1]
0310 << " S3 = " << systematicGrid[2];
0311 if (KnownStructure)
0312 cout << " --> VALID (step I)" << endl;
0313 else
0314 cout << " --> INVALID (step I)" << endl;
0315 }
0316 }
0317
0318
0319
0320
0321 bool goodIndex[3] = {true, true, true};
0322
0323 if (!KnownStructure) {
0324
0325 int nEasy = 0;
0326 for (int i = 0; i < 3; ++i) {
0327 if (EasyCoordinate[i])
0328 ++nEasy;
0329 }
0330 std::vector<double> misValX[3];
0331
0332
0333 if (nEasy > 0) {
0334 bool lastLine = false;
0335 int lowerLineLimit = 0;
0336 int upperLineLimit = 0;
0337 int oldIndex[3] = {0, 0, 0};
0338 for (int iLine = 0; iLine < nLines; ++iLine) {
0339 if (iLine == (nLines - 1))
0340 lastLine = true;
0341 XBVector = IXBValues.operator[](iLine);
0342 XBVector.getI3(index[0], index[1], index[2]);
0343 bool newIndices = false;
0344 for (int i = 0; i < 3; ++i) {
0345 if (oldIndex[i] != index[i]) {
0346 oldIndex[i] = index[i];
0347 newIndices = true;
0348 }
0349 }
0350 if (newIndices) {
0351 for (int i = 0; i < 3; ++i) {
0352 if (!EasyCoordinate[i])
0353 sort(misValX[i].begin(), misValX[i].end());
0354 }
0355 lowerLineLimit = upperLineLimit;
0356 upperLineLimit = iLine;
0357 for (int jLine = lowerLineLimit; jLine < upperLineLimit; ++jLine) {
0358 XBVector = IXBValues.operator[](jLine);
0359 XBVector.getI3(index[0], index[1], index[2]);
0360 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0361 double pnt[3] = {x1, x2, x3};
0362 for (int i = 0; i < 3; ++i) {
0363 if (!EasyCoordinate[i]) {
0364 double thisVal = 0.;
0365 double lastVal = 0.;
0366 double deltVal = 0.;
0367 for (int j = 0; j < NumberOfPoints[i]; ++j) {
0368 thisVal = misValX[i].operator[](j);
0369 if (std::abs(thisVal - pnt[i]) < EPSILON)
0370 index[i] = j;
0371 if (j == 1)
0372 deltVal = thisVal - lastVal;
0373 else if (j > 1) {
0374 if (std::abs(thisVal - lastVal - deltVal) > EPSILON)
0375 goodIndex[i] = false;
0376 }
0377 lastVal = thisVal;
0378 }
0379 }
0380 }
0381 XBVector.putI3(index[0], index[1], index[2]);
0382 XBArray.push_back(XBVector);
0383 for (int i = 0; i < 3; ++i) {
0384 if (!misValX[i].empty())
0385 misValX[i].clear();
0386 }
0387 }
0388 }
0389 XBVector = IXBValues.operator[](iLine);
0390 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0391 double pnt[3] = {x1, x2, x3};
0392 for (int i = 0; i < 3; ++i) {
0393 if (!EasyCoordinate[i]) {
0394 if (iLine == 0) {
0395 NumberOfPoints[i] = 1;
0396 misValX[i].push_back(pnt[i]);
0397 } else if (newIndices) {
0398 NumberOfPoints[i] = 1;
0399 misValX[i].push_back(pnt[i]);
0400 } else {
0401 bool knownValue = false;
0402 for (int j = 0; j < NumberOfPoints[i]; ++j) {
0403 if (std::abs(misValX[i].operator[](j) - pnt[i]) < EPSILON)
0404 knownValue = true;
0405 }
0406 if (!knownValue) {
0407 ++NumberOfPoints[i];
0408 misValX[i].push_back(pnt[i]);
0409 }
0410 }
0411 }
0412 }
0413 if (lastLine) {
0414 for (int i = 0; i < 3; ++i) {
0415 if (!EasyCoordinate[i])
0416 sort(misValX[i].begin(), misValX[i].end());
0417 }
0418 lowerLineLimit = upperLineLimit;
0419 upperLineLimit = nLines;
0420 for (int jLine = lowerLineLimit; jLine < upperLineLimit; ++jLine) {
0421 XBVector = IXBValues.operator[](jLine);
0422 XBVector.getI3(index[0], index[1], index[2]);
0423 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0424 double pnt[3] = {x1, x2, x3};
0425 for (int i = 0; i < 3; ++i) {
0426 if (!EasyCoordinate[i]) {
0427 double thisVal = 0.;
0428 double lastVal = 0.;
0429 double deltVal = 0.;
0430 for (int j = 0; j < NumberOfPoints[i]; ++j) {
0431 thisVal = misValX[i].operator[](j);
0432 if (std::abs(thisVal - pnt[i]) < EPSILON)
0433 index[i] = j;
0434 if (j == 1)
0435 deltVal = thisVal - lastVal;
0436 else if (j > 1) {
0437 if (std::abs(thisVal - lastVal - deltVal) > EPSILON)
0438 goodIndex[i] = false;
0439 }
0440 lastVal = thisVal;
0441 }
0442 }
0443 }
0444 XBVector.putI3(index[0], index[1], index[2]);
0445 XBArray.push_back(XBVector);
0446 for (int i = 0; i < 3; ++i) {
0447 if (!misValX[i].empty())
0448 misValX[i].clear();
0449 }
0450 }
0451 }
0452 }
0453 }
0454 int nGoodInd = 0;
0455 for (int i = 0; i < 3; ++i) {
0456 if (goodIndex[i])
0457 ++nGoodInd;
0458 }
0459
0460 if (nGoodInd < 3) {
0461 cout << endl;
0462 cout << "Neasy = 1: beginning of T E S T area!" << endl;
0463 cout << "# good indices = " << nGoodInd << " (" << goodIndex[0] << goodIndex[1] << goodIndex[2] << ") " << endl;
0464
0465 for (int iLine = 0; iLine < nLines; ++iLine) {
0466 XBVector = XBArray.operator[](iLine);
0467 XBVector.getI3(index[0], index[1], index[2]);
0468 for (int i = 0; i < 3; ++i) {
0469 if (!goodIndex[i])
0470 index[i] = 0;
0471 }
0472 XBVector.putI3(index[0], index[1], index[2]);
0473 }
0474 sort(XBArray.begin(), XBArray.end());
0475
0476 cout << "Neasy = 1: end of T E S T area!" << endl;
0477 cout << endl;
0478 }
0479
0480 if (nEasy > 0) {
0481
0482 sort(XBArray.begin(), XBArray.end());
0483
0484 double miniGrid[3][2][2][2] = {{{{0., 0.}, {0., 0.}}, {{0., 0.}, {0., 0.}}},
0485 {{{0., 0.}, {0., 0.}}, {{0., 0.}, {0., 0.}}},
0486 {{{0., 0.}, {0., 0.}}, {{0., 0.}, {0., 0.}}}};
0487 for (int iLine = 0; iLine < nLines; ++iLine) {
0488 XBVector = XBArray.operator[](iLine);
0489 XBVector.getI3(index[0], index[1], index[2]);
0490 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0491 double pnt[3] = {x1, x2, x3};
0492 if (index[0] < 2 && index[1] < 2 && index[2] < 2) {
0493 for (int i = 0; i < 3; ++i) {
0494 miniGrid[i][index[0]][index[1]][index[2]] = pnt[i];
0495 }
0496 }
0497 }
0498
0499
0500 BasicDistance0[0] = miniGrid[0][1][0][0] - miniGrid[0][0][0][0];
0501 BasicDistance0[1] = miniGrid[1][0][1][0] - miniGrid[1][0][0][0];
0502 BasicDistance0[2] = miniGrid[2][0][0][1] - miniGrid[2][0][0][0];
0503
0504 double disd10 = miniGrid[0][1][1][0] - miniGrid[0][0][1][0];
0505 double disd00 = BasicDistance0[0];
0506 double disd01 = miniGrid[0][1][0][1] - miniGrid[0][0][0][1];
0507 BasicDistance1[0][0] = 0.0;
0508 BasicDistance1[0][1] = disd10 - disd00;
0509 BasicDistance1[0][2] = disd01 - disd00;
0510 double dis1d0 = miniGrid[1][1][1][0] - miniGrid[1][1][0][0];
0511 double dis0d0 = BasicDistance0[1];
0512 double dis0d1 = miniGrid[1][0][1][1] - miniGrid[1][0][0][1];
0513 BasicDistance1[1][0] = dis1d0 - dis0d0;
0514 BasicDistance1[1][1] = 0.0;
0515 BasicDistance1[1][2] = dis0d1 - dis0d0;
0516 double dis10d = miniGrid[2][1][0][1] - miniGrid[2][1][0][0];
0517 double dis00d = BasicDistance0[2];
0518 double dis01d = miniGrid[2][0][1][1] - miniGrid[2][0][1][0];
0519 BasicDistance1[2][0] = dis10d - dis00d;
0520 BasicDistance1[2][1] = dis01d - dis00d;
0521 BasicDistance1[2][2] = 0.0;
0522
0523 BasicDistance2[0][0] = 0.0;
0524 BasicDistance2[0][1] = miniGrid[0][0][1][0] - miniGrid[0][0][0][0];
0525 BasicDistance2[0][2] = miniGrid[0][0][0][1] - miniGrid[0][0][0][0];
0526 BasicDistance2[1][0] = miniGrid[1][1][0][0] - miniGrid[1][0][0][0];
0527 BasicDistance2[1][1] = 0.0;
0528 BasicDistance2[1][2] = miniGrid[1][0][0][1] - miniGrid[1][0][0][0];
0529 BasicDistance2[2][0] = miniGrid[2][1][0][0] - miniGrid[2][0][0][0];
0530 BasicDistance2[2][1] = miniGrid[2][0][1][0] - miniGrid[2][0][0][0];
0531 BasicDistance2[2][2] = 0.0;
0532
0533 for (int i = 0; i < 3; ++i) {
0534 if (NumberOfPoints[i] < 2)
0535 BasicDistance0[i] = 9e9;
0536 }
0537 }
0538
0539
0540 if (!XBArray.empty())
0541 XBVector = XBArray.front();
0542 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0543 firstListPoint[0] = x1;
0544 firstListPoint[1] = x2;
0545 firstListPoint[2] = x3;
0546 if (!XBArray.empty())
0547 XBVector = XBArray.back();
0548 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0549 secondRefPoint[0] = x1;
0550 secondRefPoint[1] = x2;
0551 secondRefPoint[2] = x3;
0552
0553
0554 KnownStructure = true;
0555 for (int i = 0; i < 3; ++i) {
0556 ReferencePoint[i] = firstListPoint[i];
0557 nSteps[i] = NumberOfPoints[i] - 1;
0558
0559 double stepSize = BasicDistance0[i];
0560 for (int j = 0; j < 3; ++j) {
0561 stepSize += BasicDistance1[i][j] * nSteps[j];
0562 }
0563 double offset = 0.0;
0564 for (int j = 0; j < 3; ++j) {
0565 offset += BasicDistance2[i][j] * nSteps[j];
0566 }
0567
0568 double totDiff = secondRefPoint[i] - (ReferencePoint[i] + offset + stepSize * nSteps[i]);
0569 if (std::abs(totDiff) < EPSILON * 10.)
0570 systematicGrid[i] = true;
0571 else
0572 systematicGrid[i] = false;
0573 if (!systematicGrid[i])
0574 KnownStructure = false;
0575 }
0576
0577 if (KnownStructure) {
0578 if (XyzCoordinates)
0579 GridType = 2;
0580 if (RpzCoordinates)
0581 GridType = 4;
0582 if (PRINT) {
0583
0584 cout << " read " << nLines << " lines -> grid structure:" << endl;
0585 cout << " # of points: N1 = " << NumberOfPoints[0] << " N2 = " << NumberOfPoints[1]
0586 << " N3 = " << NumberOfPoints[2] << endl;
0587 cout << " ref. point : X1 = " << ReferencePoint[0] << " X2 = " << ReferencePoint[1]
0588 << " X3 = " << ReferencePoint[2] << endl;
0589 cout << " step parm0 : A1 = " << BasicDistance0[0] << " A2 = " << BasicDistance0[1]
0590 << " A3 = " << BasicDistance0[2] << endl;
0591 cout << " step parm1 : B11 = " << BasicDistance1[0][0] << " B21 = " << BasicDistance1[1][0]
0592 << " B31 = " << BasicDistance1[2][0] << endl;
0593 cout << " : B12 = " << BasicDistance1[0][1] << " B22 = " << BasicDistance1[1][1]
0594 << " B32 = " << BasicDistance1[2][1] << endl;
0595 cout << " : B13 = " << BasicDistance1[0][2] << " B23 = " << BasicDistance1[1][2]
0596 << " B33 = " << BasicDistance1[2][2] << endl;
0597 cout << " offset parm: O11 = " << BasicDistance2[0][0] << " O21 = " << BasicDistance2[1][0]
0598 << " O31 = " << BasicDistance2[2][0] << endl;
0599 cout << " : O12 = " << BasicDistance2[0][1] << " O22 = " << BasicDistance2[1][1]
0600 << " O32 = " << BasicDistance2[2][1] << endl;
0601 cout << " : O13 = " << BasicDistance2[0][2] << " O23 = " << BasicDistance2[1][2]
0602 << " O33 = " << BasicDistance2[2][2] << endl;
0603 cout << " easy grid : E1 = " << EasyCoordinate[0] << " E2 = " << EasyCoordinate[1]
0604 << " E3 = " << EasyCoordinate[2] << endl;
0605 cout << " : I1 = " << goodIndex[0] << " I2 = " << goodIndex[1] << " I3 = " << goodIndex[2]
0606 << endl;
0607 cout << " structure : S1 = " << systematicGrid[0] << " S2 = " << systematicGrid[1]
0608 << " S3 = " << systematicGrid[2];
0609 if (KnownStructure)
0610 cout << " --> VALID (step II)" << endl;
0611 else {
0612 cout << " --> INVALID (step II)" << endl;
0613 cout << " reason for error: ";
0614 if (NumberOfPoints[0] * NumberOfPoints[1] * NumberOfPoints[2] != nLines) {
0615 cout << endl;
0616 cout << " N1*N2*N3 =/= N.lines --> exiting now ..." << endl;
0617 return;
0618 } else {
0619 cout << " no idea so far --> exiting now ..." << endl;
0620 return;
0621 }
0622 }
0623 }
0624 }
0625 }
0626
0627
0628
0629 SixDPoint GridEntry;
0630 for (unsigned int iLine = 0; iLine < XBArray.size(); ++iLine) {
0631 XBVector = XBArray.operator[](iLine);
0632 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0633 GridEntry.putP6(x1, x2, x3, Bx, By, Bz);
0634 GridData.push_back(GridEntry);
0635 }
0636 if (!XBArray.empty())
0637 XBArray.clear();
0638
0639 return;
0640 }
0641
0642 void prepareMagneticFieldGrid::fillFromFileSpecial(const std::string &name) {
0643
0644 GridType = 6;
0645
0646 double phi = (rotateSector - 1) * Geom::pi() / 6.;
0647 double sphi = sin(-phi);
0648 double cphi = cos(-phi);
0649
0650
0651 std::string::size_type ibeg, iend;
0652 ibeg = name.find('-');
0653 iend = name.rfind('-');
0654 if ((name.substr(ibeg + 1, iend - ibeg - 1)) == "xyz") {
0655 XyzCoordinates = true;
0656 } else if ((name.substr(ibeg + 1, iend - ibeg - 1)) == "rpz") {
0657 RpzCoordinates = true;
0658 } else {
0659 cout << " Unrecognized input file name: valid formats are *-xyz-* or *-rpz-*" << endl;
0660 abort();
0661 }
0662
0663
0664 IndexedDoubleVector XBVector;
0665 std::vector<IndexedDoubleVector> XBValues;
0666 std::vector<IndexedDoubleVector> IXBValues;
0667 std::vector<IndexedDoubleVector> XBArray;
0668
0669
0670 std::vector<double> valX[3];
0671 std::vector<int> numX[3];
0672
0673
0674 int nLines = 0;
0675 int index[3] = {0, 0, 0};
0676 int nSteps[3] = {0, 0, 0};
0677 double x1, x2, x3, Bx, By, Bz, perm, poten;
0678 std::ifstream file(name.c_str());
0679 if (file.good()) {
0680 while (file.good()) {
0681 file >> x1 >> x2 >> x3 >> Bx >> By >> Bz >> perm >> poten;
0682 if (file) {
0683 if (rotateSector > 1) {
0684 if (RpzCoordinates) {
0685 x2 = x2 - phi;
0686 if (fabs(x2) > 0.78) {
0687 cout << "ERROR: Input coordinates do not belong to sector " << rotateSector << " : " << x2 << endl;
0688 abort();
0689 }
0690 } else {
0691 double x = x1;
0692 double y = x2;
0693 x1 = x * cphi - y * sphi;
0694 x2 = x * sphi + y * cphi;
0695 if (fabs(atan2(x2, x1)) > 0.78) {
0696 cout << "ERROR: Input coordinates do not belong to sector " << rotateSector << " : " << atan2(x2, x1)
0697 << endl;
0698 abort();
0699 }
0700 }
0701 double Bx0 = Bx;
0702 double By0 = By;
0703 Bx = Bx0 * cphi - By0 * sphi;
0704 By = Bx0 * sphi + By0 * cphi;
0705 }
0706
0707 XBVector.putI3(index[0], index[1], index[2]);
0708 XBVector.putV6(x1, x2, x3, Bx, By, Bz);
0709 XBValues.push_back(XBVector);
0710 ++nLines;
0711
0712
0713 double pnt[3] = {x1, x2, x3};
0714 if (nLines == 1) {
0715 for (int i = 0; i < 3; ++i) {
0716 valX[i].push_back(pnt[i]);
0717 numX[i].push_back(1);
0718 }
0719 } else {
0720 for (int i = 0; i < 3; ++i) {
0721 bool knownValue = false;
0722 for (int j = 0; j < (nSteps[i] + 1); ++j) {
0723 if (std::abs(valX[i].operator[](j) - pnt[i]) < EPSILON)
0724 knownValue = true;
0725 if (std::abs(valX[i].operator[](j) - pnt[i]) < EPSILON)
0726 numX[i].operator[](j) += 1;
0727 }
0728 if (!knownValue) {
0729 valX[i].push_back(pnt[i]);
0730 nSteps[i]++;
0731 numX[i].push_back(1);
0732 }
0733 }
0734 }
0735 }
0736 }
0737 } else
0738 return;
0739
0740
0741 for (int i = 0; i < 3; ++i) {
0742
0743 sort(valX[i].begin(), valX[i].end());
0744
0745 EasyCoordinate[i] = true;
0746 NumberOfPoints[i] = nSteps[i] + 1;
0747 BasicDistance0[i] = 9e9;
0748 double newVal = 0.0;
0749 double oldVal = 0.0;
0750 for (int j = 0; j < NumberOfPoints[i]; ++j) {
0751 newVal = valX[i].operator[](j);
0752 if (j == 1)
0753 BasicDistance0[i] = newVal - oldVal;
0754 else if (j > 1) {
0755 if (std::abs(newVal - oldVal - BasicDistance0[i]) > EPSILON)
0756 EasyCoordinate[i] = false;
0757 }
0758 oldVal = newVal;
0759 }
0760 }
0761
0762 std::cout << std::endl;
0763
0764
0765 for (int iLine = 0; iLine < nLines; ++iLine) {
0766 XBVector = XBValues.operator[](iLine);
0767 XBVector.getI3(index[0], index[1], index[2]);
0768 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0769 double pnt[3] = {x1, x2, x3};
0770 for (int i = 0; i < 3; ++i) {
0771 for (int j = 0; j < NumberOfPoints[i]; ++j) {
0772 if (EasyCoordinate[i]) {
0773 if (std::abs(valX[i].operator[](j) - pnt[i]) < EPSILON)
0774 index[i] = j;
0775 }
0776 }
0777 }
0778 XBVector.putI3(index[0], index[1], index[2]);
0779 IXBValues.push_back(XBVector);
0780 }
0781 if (!XBValues.empty())
0782 XBValues.clear();
0783
0784
0785 sort(IXBValues.begin(), IXBValues.end());
0786
0787
0788 bool systematicGrid[3] = {false, false, false};
0789 KnownStructure = true;
0790 if (NumberOfPoints[0] * NumberOfPoints[1] * NumberOfPoints[2] != nLines)
0791 KnownStructure = false;
0792 for (int i = 0; i < 3; ++i) {
0793 systematicGrid[i] = EasyCoordinate[i];
0794 if (!systematicGrid[i])
0795 KnownStructure = false;
0796 }
0797
0798 if (KnownStructure == true) {
0799 for (int iLine = 0; iLine < nLines; ++iLine) {
0800 XBVector = IXBValues.operator[](iLine);
0801 XBArray.push_back(XBVector);
0802 }
0803 }
0804
0805
0806 if (!XBArray.empty())
0807 XBVector = XBArray.front();
0808 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0809 double firstListPoint[3] = {x1, x2, x3};
0810 if (!XBArray.empty())
0811 XBVector = XBArray.back();
0812 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0813 double secondRefPoint[3] = {x1, x2, x3};
0814
0815
0816 KnownStructure = true;
0817 for (int i = 0; i < 3; ++i) {
0818 ReferencePoint[i] = firstListPoint[i];
0819
0820 double stepSize = BasicDistance0[i];
0821
0822 double totDiff = secondRefPoint[i] - (ReferencePoint[i] + stepSize * nSteps[i]);
0823 if (std::abs(totDiff) < EPSILON * 10.)
0824 systematicGrid[i] = true;
0825 else
0826 systematicGrid[i] = false;
0827 if (!systematicGrid[i])
0828 KnownStructure = false;
0829 }
0830 bool goodIndex[3] = {true, true, true};
0831
0832 if (!KnownStructure) {
0833
0834 int nEasy = 0;
0835 for (int i = 0; i < 3; ++i) {
0836 if (EasyCoordinate[i])
0837 ++nEasy;
0838 }
0839 std::vector<double> misValX[3];
0840
0841
0842 if (nEasy > 0) {
0843 bool lastLine = false;
0844 int lowerLineLimit = 0;
0845 int upperLineLimit = 0;
0846 int oldIndex[3] = {0, 0, 0};
0847 for (int iLine = 0; iLine < nLines; ++iLine) {
0848 if (iLine == (nLines - 1))
0849 lastLine = true;
0850 XBVector = IXBValues.operator[](iLine);
0851 XBVector.getI3(index[0], index[1], index[2]);
0852 bool newIndices = false;
0853 for (int i = 0; i < 3; ++i) {
0854 if (oldIndex[i] != index[i]) {
0855 oldIndex[i] = index[i];
0856 newIndices = true;
0857 }
0858 }
0859 if (newIndices) {
0860 for (int i = 0; i < 3; ++i) {
0861 if (!EasyCoordinate[i])
0862 sort(misValX[i].begin(), misValX[i].end());
0863 }
0864 lowerLineLimit = upperLineLimit;
0865 upperLineLimit = iLine;
0866 for (int jLine = lowerLineLimit; jLine < upperLineLimit; ++jLine) {
0867 XBVector = IXBValues.operator[](jLine);
0868 XBVector.getI3(index[0], index[1], index[2]);
0869 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0870 double pnt[3] = {x1, x2, x3};
0871 for (int i = 0; i < 3; ++i) {
0872 if (!EasyCoordinate[i]) {
0873 double thisVal = 0.;
0874 double lastVal = 0.;
0875 double deltVal = 0.;
0876 for (int j = 0; j < NumberOfPoints[i]; ++j) {
0877 thisVal = misValX[i].operator[](j);
0878 if (std::abs(thisVal - pnt[i]) < EPSILON)
0879 index[i] = j;
0880 if (j == 1)
0881 deltVal = thisVal - lastVal;
0882 else if (j > 1) {
0883 if (std::abs(thisVal - lastVal - deltVal) > EPSILON)
0884 goodIndex[i] = false;
0885 }
0886 lastVal = thisVal;
0887 }
0888 }
0889 }
0890 XBVector.putI3(index[0], index[1], index[2]);
0891 XBArray.push_back(XBVector);
0892 for (int i = 0; i < 3; ++i) {
0893 if (!misValX[i].empty())
0894 misValX[i].clear();
0895 }
0896 }
0897 }
0898 XBVector = IXBValues.operator[](iLine);
0899 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0900 double pnt[3] = {x1, x2, x3};
0901 for (int i = 0; i < 3; ++i) {
0902 if (!EasyCoordinate[i]) {
0903 if (iLine == 0) {
0904 NumberOfPoints[i] = 1;
0905 misValX[i].push_back(pnt[i]);
0906 } else if (newIndices) {
0907 NumberOfPoints[i] = 1;
0908 misValX[i].push_back(pnt[i]);
0909 } else {
0910 bool knownValue = false;
0911 for (int j = 0; j < NumberOfPoints[i]; ++j) {
0912 if (std::abs(misValX[i].operator[](j) - pnt[i]) < EPSILON)
0913 knownValue = true;
0914 }
0915 if (!knownValue) {
0916 ++NumberOfPoints[i];
0917 misValX[i].push_back(pnt[i]);
0918 }
0919 }
0920 }
0921 }
0922 if (lastLine) {
0923 for (int i = 0; i < 3; ++i) {
0924 if (!EasyCoordinate[i])
0925 sort(misValX[i].begin(), misValX[i].end());
0926 }
0927 lowerLineLimit = upperLineLimit;
0928 upperLineLimit = nLines;
0929 for (int jLine = lowerLineLimit; jLine < upperLineLimit; ++jLine) {
0930 XBVector = IXBValues.operator[](jLine);
0931 XBVector.getI3(index[0], index[1], index[2]);
0932 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0933 double pnt[3] = {x1, x2, x3};
0934 for (int i = 0; i < 3; ++i) {
0935 if (!EasyCoordinate[i]) {
0936 double thisVal = 0.;
0937 double lastVal = 0.;
0938 double deltVal = 0.;
0939 for (int j = 0; j < NumberOfPoints[i]; ++j) {
0940 thisVal = misValX[i].operator[](j);
0941 if (std::abs(thisVal - pnt[i]) < EPSILON)
0942 index[i] = j;
0943 if (j == 1)
0944 deltVal = thisVal - lastVal;
0945 else if (j > 1) {
0946 if (std::abs(thisVal - lastVal - deltVal) > EPSILON)
0947 goodIndex[i] = false;
0948 }
0949 lastVal = thisVal;
0950 }
0951 }
0952 }
0953 XBVector.putI3(index[0], index[1], index[2]);
0954 XBArray.push_back(XBVector);
0955 for (int i = 0; i < 3; ++i) {
0956 if (!misValX[i].empty())
0957 misValX[i].clear();
0958 }
0959 }
0960 }
0961 }
0962 }
0963
0964 if (nEasy > 0) {
0965
0966 sort(XBArray.begin(), XBArray.end());
0967
0968 for (int i = 0; i < 3; ++i) {
0969 nSteps[i] = NumberOfPoints[i] - 1;
0970 }
0971 std::vector<double> rMinVec;
0972 std::vector<double> rMaxVec;
0973 std::vector<double> phiVec;
0974 for (int iLine = 0; iLine < nLines; ++iLine) {
0975 XBVector = XBArray.operator[](iLine);
0976 XBVector.getI3(index[0], index[1], index[2]);
0977
0978 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0979 if (index[2] == 0) {
0980 if (index[0] == 0)
0981 rMinVec.push_back(x1);
0982 if (index[0] == nSteps[0])
0983 rMaxVec.push_back(x1);
0984 if (index[0] == nSteps[0])
0985 phiVec.push_back(x2);
0986 }
0987 }
0988 for (int j = 0; j < NumberOfPoints[1]; ++j) {
0989 double phi = phiVec.operator[](j);
0990 double rMin = rMinVec.operator[](j);
0991 double rMax = rMaxVec.operator[](j);
0992 double rhoMin, rhoMax;
0993 if (GridType == 6) {
0994 rhoMin = rMin * cos(phi);
0995 rhoMax = rMax * cos(phi);
0996 } else if (GridType == 5) {
0997 rhoMin = rMin * sin(phi);
0998 rhoMax = rMax * sin(phi);
0999 } else {
1000 cout << "unknown grid type!" << endl;
1001 abort();
1002 }
1003
1004 if (j == 0) {
1005 RParAsFunOfPhi[0] = rMax;
1006 RParAsFunOfPhi[1] = rhoMax;
1007 RParAsFunOfPhi[2] = rMin;
1008 RParAsFunOfPhi[3] = rhoMin;
1009 } else {
1010 if (std::abs(rMax - RParAsFunOfPhi[0]) > EPSILON)
1011 RParAsFunOfPhi[0] = 0.;
1012 if (std::abs(rhoMax - RParAsFunOfPhi[1]) > EPSILON)
1013 RParAsFunOfPhi[1] = 0.;
1014 if (std::abs(rMin - RParAsFunOfPhi[2]) > EPSILON)
1015 RParAsFunOfPhi[2] = 0.;
1016 if (std::abs(rhoMin - RParAsFunOfPhi[3]) > EPSILON)
1017 RParAsFunOfPhi[3] = 0.;
1018 }
1019 }
1020
1021
1022 for (int i = 0; i < 3; ++i) {
1023 if (NumberOfPoints[i] < 2)
1024 BasicDistance0[i] = 9e9;
1025 }
1026 }
1027
1028
1029 if (!XBArray.empty())
1030 XBVector = XBArray.front();
1031 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
1032 firstListPoint[0] = x1;
1033 firstListPoint[1] = x2;
1034 firstListPoint[2] = x3;
1035 if (!XBArray.empty())
1036 XBVector = XBArray.back();
1037 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
1038 secondRefPoint[0] = x1;
1039 secondRefPoint[1] = x2;
1040 secondRefPoint[2] = x3;
1041
1042
1043 KnownStructure = true;
1044 for (int i = 0; i < 3; ++i) {
1045 ReferencePoint[i] = firstListPoint[i];
1046 nSteps[i] = NumberOfPoints[i] - 1;
1047 }
1048
1049 double sinPhi;
1050 if (GridType == 6) {
1051 sinPhi = cos(secondRefPoint[1]);
1052 } else if (GridType == 5) {
1053 sinPhi = sin(secondRefPoint[1]);
1054 } else {
1055 cout << "unknown GridType!" << endl;
1056 abort();
1057 }
1058
1059 if (std::abs(sinPhi) < EPSILON) {
1060 cout << " WARNING: unexpected sinPhi parameter = 0" << endl;
1061 sinPhi = EPSILON;
1062 }
1063
1064 double totStepSize =
1065 RParAsFunOfPhi[0] + RParAsFunOfPhi[1] / sinPhi - RParAsFunOfPhi[2] - RParAsFunOfPhi[3] / sinPhi;
1066 double startingPoint = RParAsFunOfPhi[2] + RParAsFunOfPhi[3] / sinPhi;
1067 double totDiff = secondRefPoint[0] - (startingPoint + totStepSize);
1068 if (std::abs(totDiff) < EPSILON * 10.)
1069 systematicGrid[0] = true;
1070 else
1071 systematicGrid[0] = false;
1072 if (!systematicGrid[0])
1073 KnownStructure = false;
1074 totDiff = secondRefPoint[1] - (ReferencePoint[1] + BasicDistance0[1] * nSteps[1]);
1075 if (std::abs(totDiff) < EPSILON * 10.)
1076 systematicGrid[1] = true;
1077 else
1078 systematicGrid[1] = false;
1079 if (!systematicGrid[1])
1080 KnownStructure = false;
1081 totDiff = secondRefPoint[2] - (ReferencePoint[2] + BasicDistance0[2] * nSteps[2]);
1082 if (std::abs(totDiff) < EPSILON * 10.)
1083 systematicGrid[2] = true;
1084 else
1085 systematicGrid[2] = false;
1086 if (!systematicGrid[2])
1087 KnownStructure = false;
1088
1089 if (PRINT) {
1090
1091 cout << " read " << nLines << " lines -> grid structure:" << endl;
1092 cout << " # of points: N1 = " << NumberOfPoints[0] << " N2 = " << NumberOfPoints[1]
1093 << " N3 = " << NumberOfPoints[2] << endl;
1094 cout << " ref. point : X1 = " << ReferencePoint[0] << " X2 = " << ReferencePoint[1]
1095 << " X3 = " << ReferencePoint[2] << endl;
1096 cout << " step parm0 : A1 = " << BasicDistance0[0] << " A2 = " << BasicDistance0[1]
1097 << " A3 = " << BasicDistance0[2] << endl;
1098 cout << " r param.s. : R1 = " << RParAsFunOfPhi[0] << " RHO1 = " << RParAsFunOfPhi[1]
1099 << " R2 = " << RParAsFunOfPhi[2] << " RHO2 = " << RParAsFunOfPhi[3] << endl;
1100 cout << " easy grid : E1 = " << EasyCoordinate[0] << " E2 = " << EasyCoordinate[1]
1101 << " E3 = " << EasyCoordinate[2] << endl;
1102 cout << " : I1 = " << goodIndex[0] << " I2 = " << goodIndex[1] << " I3 = " << goodIndex[2]
1103 << endl;
1104 cout << " structure : S1 = " << systematicGrid[0] << " S2 = " << systematicGrid[1]
1105 << " S3 = " << systematicGrid[2];
1106 if (KnownStructure)
1107 cout << " --> VALID (step II)" << endl;
1108 else {
1109 cout << " --> INVALID (step II)" << endl;
1110 cout << " reason for error: ";
1111 if (NumberOfPoints[0] * NumberOfPoints[1] * NumberOfPoints[2] != nLines) {
1112 cout << endl;
1113 cout << " N1*N2*N3 =/= N.lines --> exiting now ..." << endl;
1114 return;
1115 } else {
1116 cout << " no idea so far --> exiting now ..." << endl;
1117 return;
1118 }
1119 }
1120 }
1121 }
1122
1123
1124 SixDPoint GridEntry;
1125 for (unsigned int iLine = 0; iLine < XBArray.size(); ++iLine) {
1126 XBVector = XBArray.operator[](iLine);
1127 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
1128 GridEntry.putP6(x1, x2, x3, Bx, By, Bz);
1129 GridData.push_back(GridEntry);
1130 }
1131 if (!XBArray.empty())
1132 XBArray.clear();
1133
1134 return;
1135 }
1136
1137 int prepareMagneticFieldGrid::gridType() {
1138 int type = GridType;
1139 if (PRINT) {
1140 if (type == 0)
1141 cout << " grid type = " << type << " --> not determined" << endl;
1142 if (type == 1)
1143 cout << " grid type = " << type << " --> (x,y,z) cube" << endl;
1144 if (type == 2)
1145 cout << " grid type = " << type << " --> (x,y,z) trapezoid" << endl;
1146 if (type == 3)
1147 cout << " grid type = " << type << " --> (r,phi,z) cube" << endl;
1148 if (type == 4)
1149 cout << " grid type = " << type << " --> (r,phi,z) trapezoid" << endl;
1150 if (type == 5)
1151 cout << " grid type = " << type << " --> (r,phi,z) 1/sin(phi)" << endl;
1152 if (type == 6)
1153 cout << " grid type = " << type << " --> (r,phi,z) 1/cos(phi)" << endl;
1154 }
1155 return type;
1156 }
1157
1158 void prepareMagneticFieldGrid::validateAllPoints() {
1159 if (!KnownStructure)
1160 return;
1161
1162 double x1, x2, x3, Bx, By, Bz;
1163 int numberOfErrors = 0;
1164
1165
1166 int index[3];
1167
1168 if (GridType < 5) {
1169 for (index[0] = 0; index[0] < NumberOfPoints[0]; ++index[0]) {
1170 for (index[1] = 0; index[1] < NumberOfPoints[1]; ++index[1]) {
1171 for (index[2] = 0; index[2] < NumberOfPoints[2]; ++index[2]) {
1172
1173 putIndicesGetXAndB(index[0], index[1], index[2], x1, x2, x3, Bx, By, Bz);
1174
1175
1176 double pnt[3] = {x1, x2, x3};
1177 int tmpIdx[3] = {999, 999, 999};
1178 putCoordGetIndices((x1 + EPSILON), (x2 + EPSILON), (x3 + EPSILON), tmpIdx[0], tmpIdx[1], tmpIdx[2]);
1179 for (int i = 0; i < 3; ++i) {
1180 if (tmpIdx[i] != index[i])
1181 ++numberOfErrors;
1182 }
1183
1184
1185 double tmpPnt[3] = {0.0, 0.0, 0.0};
1186 putIndCalcXReturnB(index[0], index[1], index[2], tmpPnt[0], tmpPnt[1], tmpPnt[2], Bx, By, Bz);
1187 for (int i = 0; i < 3; ++i) {
1188 if (std::abs(tmpPnt[i] - pnt[i]) > EPSILON)
1189 ++numberOfErrors;
1190 }
1191 }
1192 }
1193 }
1194 }
1195
1196 if (GridType == 5 || GridType == 6) {
1197 for (index[0] = 0; index[0] < NumberOfPoints[0]; ++index[0]) {
1198 for (index[1] = 0; index[1] < NumberOfPoints[1]; ++index[1]) {
1199 for (index[2] = 0; index[2] < NumberOfPoints[2]; ++index[2]) {
1200
1201 putIndicesGetXAndB(index[0], index[1], index[2], x1, x2, x3, Bx, By, Bz);
1202
1203
1204 double pnt[3] = {x1, x2, x3};
1205 int tmpIdx[3] = {999, 999, 999};
1206 putCoordGetIndices((x1 + EPSILON), (x2 + EPSILON), (x3 + EPSILON), tmpIdx[0], tmpIdx[1], tmpIdx[2]);
1207 for (int i = 0; i < 3; ++i) {
1208 if (tmpIdx[i] != index[i]) {
1209 putCoordGetIndices((x1 - EPSILON), (x2 + EPSILON), (x3 + EPSILON), tmpIdx[0], tmpIdx[1], tmpIdx[2]);
1210 if (tmpIdx[i] != index[i]) {
1211 putCoordGetIndices((x1 + EPSILON), (x2 - EPSILON), (x3 + EPSILON), tmpIdx[0], tmpIdx[1], tmpIdx[2]);
1212 if (tmpIdx[i] != index[i]) {
1213 putCoordGetIndices((x1 - EPSILON), (x2 - EPSILON), (x3 + EPSILON), tmpIdx[0], tmpIdx[1], tmpIdx[2]);
1214 if (tmpIdx[i] != index[i])
1215 ++numberOfErrors;
1216 }
1217 }
1218 }
1219 }
1220
1221
1222 double tmpPnt[3] = {0.0, 0.0, 0.0};
1223 putIndCalcXReturnB(index[0], index[1], index[2], tmpPnt[0], tmpPnt[1], tmpPnt[2], Bx, By, Bz);
1224 for (int i = 0; i < 3; ++i) {
1225 if (std::abs(tmpPnt[i] - pnt[i]) > EPSILON)
1226 ++numberOfErrors;
1227 }
1228 }
1229 }
1230 }
1231 }
1232
1233 if (numberOfErrors > 0)
1234 KnownStructure = false;
1235 if (PRINT)
1236 cout << " grid validation of all points done --> # errors = " << numberOfErrors << endl;
1237
1238 return;
1239 }
1240
1241 void prepareMagneticFieldGrid::saveGridToFile(const std::string &outName) {
1242
1243 binary_ofstream outFile(outName);
1244
1245 outFile << GridType;
1246
1247 convertUnits();
1248 if (GridType == 1) {
1249 outFile << NumberOfPoints[0] << NumberOfPoints[1] << NumberOfPoints[2];
1250 outFile << ReferencePoint[0] << ReferencePoint[1] << ReferencePoint[2];
1251 outFile << BasicDistance0[0] << BasicDistance0[1] << BasicDistance0[2];
1252 }
1253 if (GridType == 2) {
1254 outFile << NumberOfPoints[0] << NumberOfPoints[1] << NumberOfPoints[2];
1255 outFile << ReferencePoint[0] << ReferencePoint[1] << ReferencePoint[2];
1256 outFile << BasicDistance0[0] << BasicDistance0[1] << BasicDistance0[2];
1257 outFile << BasicDistance1[0][0] << BasicDistance1[1][0] << BasicDistance1[2][0];
1258 outFile << BasicDistance1[0][1] << BasicDistance1[1][1] << BasicDistance1[2][1];
1259 outFile << BasicDistance1[0][2] << BasicDistance1[1][2] << BasicDistance1[2][2];
1260 outFile << BasicDistance2[0][0] << BasicDistance2[1][0] << BasicDistance2[2][0];
1261 outFile << BasicDistance2[0][1] << BasicDistance2[1][1] << BasicDistance2[2][1];
1262 outFile << BasicDistance2[0][2] << BasicDistance2[1][2] << BasicDistance2[2][2];
1263 outFile << EasyCoordinate[0] << EasyCoordinate[1] << EasyCoordinate[2];
1264 }
1265 if (GridType == 3) {
1266 outFile << NumberOfPoints[0] << NumberOfPoints[1] << NumberOfPoints[2];
1267 outFile << ReferencePoint[0] << ReferencePoint[1] << ReferencePoint[2];
1268 outFile << BasicDistance0[0] << BasicDistance0[1] << BasicDistance0[2];
1269 }
1270 if (GridType == 4) {
1271 outFile << NumberOfPoints[0] << NumberOfPoints[1] << NumberOfPoints[2];
1272 outFile << ReferencePoint[0] << ReferencePoint[1] << ReferencePoint[2];
1273 outFile << BasicDistance0[0] << BasicDistance0[1] << BasicDistance0[2];
1274 outFile << BasicDistance1[0][0] << BasicDistance1[1][0] << BasicDistance1[2][0];
1275 outFile << BasicDistance1[0][1] << BasicDistance1[1][1] << BasicDistance1[2][1];
1276 outFile << BasicDistance1[0][2] << BasicDistance1[1][2] << BasicDistance1[2][2];
1277 outFile << BasicDistance2[0][0] << BasicDistance2[1][0] << BasicDistance2[2][0];
1278 outFile << BasicDistance2[0][1] << BasicDistance2[1][1] << BasicDistance2[2][1];
1279 outFile << BasicDistance2[0][2] << BasicDistance2[1][2] << BasicDistance2[2][2];
1280 outFile << EasyCoordinate[0] << EasyCoordinate[1] << EasyCoordinate[2];
1281 }
1282 if (GridType == 5 || GridType == 6) {
1283 outFile << NumberOfPoints[0] << NumberOfPoints[1] << NumberOfPoints[2];
1284 outFile << ReferencePoint[0] << ReferencePoint[1] << ReferencePoint[2];
1285 outFile << BasicDistance0[0] << BasicDistance0[1] << BasicDistance0[2];
1286 outFile << RParAsFunOfPhi[0] << RParAsFunOfPhi[1] << RParAsFunOfPhi[2] << RParAsFunOfPhi[3];
1287 }
1288 int nLines = NumberOfPoints[0] * NumberOfPoints[1] * NumberOfPoints[2];
1289 SixDPoint GridPoint;
1290
1291 for (int iLine = 0; iLine < nLines; ++iLine) {
1292 GridPoint = GridData.operator[](iLine);
1293 float Bx = float(GridPoint.bx());
1294 float By = float(GridPoint.by());
1295 float Bz = float(GridPoint.bz());
1296 outFile << Bx << By << Bz;
1297
1298 }
1299
1300 const std::string lastEntry = "complete";
1301 outFile << lastEntry;
1302 outFile.close();
1303 if (PRINT)
1304 cout << " output " << outName << endl;
1305 return;
1306 }
1307
1308 void prepareMagneticFieldGrid::convertUnits() {
1309 double cm = 100.;
1310 if (XyzCoordinates) {
1311 for (int i = 0; i < 3; ++i) {
1312 ReferencePoint[i] *= cm;
1313 };
1314 for (int i = 0; i < 3; ++i) {
1315 BasicDistance0[i] *= cm;
1316 };
1317 for (int i = 0; i < 3; ++i) {
1318 for (int j = 0; j < 3; ++j) {
1319 BasicDistance1[i][j] *= cm;
1320 };
1321 };
1322 for (int i = 0; i < 3; ++i) {
1323 for (int j = 0; j < 3; ++j) {
1324 BasicDistance2[i][j] *= cm;
1325 };
1326 };
1327 for (int i = 0; i < 4; ++i) {
1328 RParAsFunOfPhi[i] *= cm;
1329 };
1330 }
1331 double du[3] = {100., 1., 100.};
1332 if (RpzCoordinates) {
1333 for (int i = 0; i < 3; ++i) {
1334 ReferencePoint[i] *= du[i];
1335 };
1336 for (int i = 0; i < 3; ++i) {
1337 BasicDistance0[i] *= du[i];
1338 };
1339 for (int i = 0; i < 3; ++i) {
1340 for (int j = 0; j < 3; ++j) {
1341 BasicDistance1[i][j] *= du[i];
1342 };
1343 };
1344 for (int i = 0; i < 3; ++i) {
1345 for (int j = 0; j < 3; ++j) {
1346 BasicDistance2[i][j] *= du[i];
1347 };
1348 };
1349 for (int i = 0; i < 4; ++i) {
1350 RParAsFunOfPhi[i] *= cm;
1351 };
1352 }
1353 return;
1354 }
1355
1356 bool prepareMagneticFieldGrid::isReady() { return KnownStructure; }
1357
1358 void prepareMagneticFieldGrid::interpolateAtPoint(double X1, double X2, double X3, double &Bx, double &By, double &Bz) {
1359 Bx = By = Bz = 0.0;
1360 if (KnownStructure) {
1361
1362 VectorFieldInterpolation MagInterpol;
1363
1364 int index[3];
1365 putCoordGetIndices(X1, X2, X3, index[0], index[1], index[2]);
1366 int index0[3] = {0, 0, 0};
1367 int index1[3] = {0, 0, 0};
1368 for (int i = 0; i < 3; ++i) {
1369 if (NumberOfPoints[i] > 1) {
1370 index0[i] = std::max(0, index[i]);
1371 if (index0[i] > NumberOfPoints[i] - 2)
1372 index0[i] = NumberOfPoints[i] - 2;
1373 ;
1374 index1[i] = std::max(1, index[i] + 1);
1375 ;
1376 if (index1[i] > NumberOfPoints[i] - 1)
1377 index1[i] = NumberOfPoints[i] - 1;
1378 }
1379 }
1380 double tmpX[3];
1381 double tmpB[3];
1382
1383 putIndicesGetXAndB(index0[0], index0[1], index0[2], tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1384 MagInterpol.defineCellPoint000(tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1385 putIndicesGetXAndB(index1[0], index0[1], index0[2], tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1386 MagInterpol.defineCellPoint100(tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1387 putIndicesGetXAndB(index0[0], index1[1], index0[2], tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1388 MagInterpol.defineCellPoint010(tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1389 putIndicesGetXAndB(index1[0], index1[1], index0[2], tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1390 MagInterpol.defineCellPoint110(tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1391 putIndicesGetXAndB(index0[0], index0[1], index1[2], tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1392 MagInterpol.defineCellPoint001(tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1393 putIndicesGetXAndB(index1[0], index0[1], index1[2], tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1394 MagInterpol.defineCellPoint101(tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1395 putIndicesGetXAndB(index0[0], index1[1], index1[2], tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1396 MagInterpol.defineCellPoint011(tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1397 putIndicesGetXAndB(index1[0], index1[1], index1[2], tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1398 MagInterpol.defineCellPoint111(tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1399
1400 MagInterpol.putSCoordGetVField(X1, X2, X3, Bx, By, Bz);
1401 }
1402
1403 return;
1404 }
1405
1406 void prepareMagneticFieldGrid::putCoordGetIndices(
1407 double X1, double X2, double X3, int &Index1, int &Index2, int &Index3) {
1408 double pnt[3] = {X1, X2, X3};
1409 int index[3] = {0, 0, 0};
1410
1411 if (GridType < 5) {
1412 for (int i = 0; i < 3; ++i) {
1413 if (EasyCoordinate[i]) {
1414 index[i] = int((pnt[i] - ReferencePoint[i]) / BasicDistance0[i]);
1415 }
1416 }
1417 for (int i = 0; i < 3; ++i) {
1418 if (!EasyCoordinate[i]) {
1419 double stepSize = BasicDistance0[i];
1420 double offset = 0.0;
1421 for (int j = 0; j < 3; ++j) {
1422 stepSize += BasicDistance1[i][j] * index[j];
1423 offset += BasicDistance2[i][j] * index[j];
1424 }
1425 index[i] = int((pnt[i] - (ReferencePoint[i] + offset)) / stepSize);
1426 }
1427 }
1428 }
1429 if (GridType == 5 || GridType == 6) {
1430 double sinPhi;
1431 if (GridType == 6) {
1432 sinPhi = cos(pnt[1]);
1433 } else if (GridType == 5) {
1434 sinPhi = sin(pnt[1]);
1435 } else {
1436 abort();
1437 }
1438
1439 if (std::abs(sinPhi) < EPSILON) {
1440 sinPhi = EPSILON;
1441 cout << "ERROR DIVISION BY ZERO" << endl;
1442 }
1443 double stepSize = RParAsFunOfPhi[0] + RParAsFunOfPhi[1] / sinPhi - RParAsFunOfPhi[2] - RParAsFunOfPhi[3] / sinPhi;
1444 stepSize = stepSize / (NumberOfPoints[0] - 1);
1445 double startingPoint = RParAsFunOfPhi[2] + RParAsFunOfPhi[3] / sinPhi;
1446 index[0] = int((pnt[0] - startingPoint) / stepSize);
1447 index[1] = int((pnt[1] - ReferencePoint[1]) / BasicDistance0[1]);
1448 index[2] = int((pnt[2] - ReferencePoint[2]) / BasicDistance0[2]);
1449 }
1450
1451 Index1 = index[0];
1452 Index2 = index[1];
1453 Index3 = index[2];
1454
1455 return;
1456 }
1457
1458 void prepareMagneticFieldGrid::putIndicesGetXAndB(
1459 int Index1, int Index2, int Index3, double &X1, double &X2, double &X3, double &Bx, double &By, double &Bz) {
1460 SixDPoint GridPoint;
1461 GridPoint = GridData.operator[](lineNumber(Index1, Index2, Index3));
1462 X1 = GridPoint.x1();
1463 X2 = GridPoint.x2();
1464 X3 = GridPoint.x3();
1465 Bx = GridPoint.bx();
1466 By = GridPoint.by();
1467 Bz = GridPoint.bz();
1468
1469 return;
1470 }
1471
1472 void prepareMagneticFieldGrid::putIndCalcXReturnB(
1473 int Index1, int Index2, int Index3, double &X1, double &X2, double &X3, double &Bx, double &By, double &Bz) {
1474 int index[3] = {Index1, Index2, Index3};
1475 double pnt[3] = {0., 0., 0.};
1476
1477 if (GridType < 5) {
1478 for (int i = 0; i < 3; ++i) {
1479 if (EasyCoordinate[i]) {
1480 pnt[i] = ReferencePoint[i] + BasicDistance0[i] * index[i];
1481 } else {
1482 double stepSize = BasicDistance0[i];
1483 double offset = 0.0;
1484 for (int j = 0; j < 3; ++j) {
1485 stepSize += BasicDistance1[i][j] * index[j];
1486 offset += BasicDistance2[i][j] * index[j];
1487 }
1488 pnt[i] = ReferencePoint[i] + offset + stepSize * index[i];
1489 }
1490 }
1491 }
1492 if (GridType == 5 || GridType == 6) {
1493 pnt[2] = ReferencePoint[2] + BasicDistance0[2] * index[2];
1494 pnt[1] = ReferencePoint[1] + BasicDistance0[1] * index[1];
1495
1496 double sinPhi;
1497 if (GridType == 6) {
1498 sinPhi = cos(pnt[1]);
1499 } else if (GridType == 5) {
1500 sinPhi = sin(pnt[1]);
1501 } else {
1502 cout << "unknown GridType!" << endl;
1503 abort();
1504 }
1505
1506 if (std::abs(sinPhi) < EPSILON) {
1507 sinPhi = EPSILON;
1508 cout << "ERROR DIVISION BY ZERO" << endl;
1509 }
1510 double stepSize = RParAsFunOfPhi[0] + RParAsFunOfPhi[1] / sinPhi - RParAsFunOfPhi[2] - RParAsFunOfPhi[3] / sinPhi;
1511 stepSize = stepSize / (NumberOfPoints[0] - 1);
1512 double startingPoint = RParAsFunOfPhi[2] + RParAsFunOfPhi[3] / sinPhi;
1513 pnt[0] = startingPoint + stepSize * index[0];
1514 }
1515
1516 X1 = pnt[0];
1517 X2 = pnt[1];
1518 X3 = pnt[2];
1519
1520 SixDPoint GridPoint;
1521 GridPoint = GridData.operator[](lineNumber(Index1, Index2, Index3));
1522 Bx = GridPoint.bx();
1523 By = GridPoint.by();
1524 Bz = GridPoint.bz();
1525
1526 return;
1527 }
1528
1529 int prepareMagneticFieldGrid::lineNumber(int Index1, int Index2, int Index3) {
1530 return Index1 * NumberOfPoints[1] * NumberOfPoints[2] + Index2 * NumberOfPoints[2] + Index3;
1531 }
1532
1533 bool prepareMagneticFieldGrid::IndexedDoubleVector::operator<(const IndexedDoubleVector &x) const {
1534 if (I3[0] < x.I3[0])
1535 return true;
1536 else if (I3[0] == x.I3[0]) {
1537 if (I3[1] < x.I3[1])
1538 return true;
1539 else if (I3[1] == x.I3[1]) {
1540 if (I3[2] < x.I3[2])
1541 return true;
1542 else
1543 return false;
1544 } else
1545 return false;
1546 } else
1547 return false;
1548 }
1549
1550 void prepareMagneticFieldGrid::IndexedDoubleVector::putI3(int index1, int index2, int index3) {
1551 I3[0] = index1;
1552 I3[1] = index2;
1553 I3[2] = index3;
1554 return;
1555 }
1556
1557 void prepareMagneticFieldGrid::IndexedDoubleVector::getI3(int &index1, int &index2, int &index3) {
1558 index1 = I3[0];
1559 index2 = I3[1];
1560 index3 = I3[2];
1561 return;
1562 }
1563
1564 void prepareMagneticFieldGrid::IndexedDoubleVector::putV6(
1565 double X1, double X2, double X3, double Bx, double By, double Bz) {
1566 V6[0] = X1;
1567 V6[1] = X2;
1568 V6[2] = X3;
1569 V6[3] = Bx;
1570 V6[4] = By;
1571 V6[5] = Bz;
1572 return;
1573 }
1574
1575 void prepareMagneticFieldGrid::IndexedDoubleVector::getV6(
1576 double &X1, double &X2, double &X3, double &Bx, double &By, double &Bz) {
1577 X1 = V6[0];
1578 X2 = V6[1];
1579 X3 = V6[2];
1580 Bx = V6[3];
1581 By = V6[4];
1582 Bz = V6[5];
1583 return;
1584 }
1585
1586 void prepareMagneticFieldGrid::SixDPoint::putP6(double X1, double X2, double X3, double Bx, double By, double Bz) {
1587 P6[0] = X1;
1588 P6[1] = X2;
1589 P6[2] = X3;
1590 P6[3] = Bx;
1591 P6[4] = By;
1592 P6[5] = Bz;
1593 return;
1594 }
1595
1596 double prepareMagneticFieldGrid::SixDPoint::x1() { return P6[0]; }
1597
1598 double prepareMagneticFieldGrid::SixDPoint::x2() { return P6[1]; }
1599
1600 double prepareMagneticFieldGrid::SixDPoint::x3() { return P6[2]; }
1601
1602 double prepareMagneticFieldGrid::SixDPoint::bx() { return P6[3]; }
1603
1604 double prepareMagneticFieldGrid::SixDPoint::by() { return P6[4]; }
1605
1606 double prepareMagneticFieldGrid::SixDPoint::bz() { return P6[5]; }