Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-02-14 13:31:24 

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   // define vectors of IndexedDoubleVectors
0018   IndexedDoubleVector XBVector;
0019   std::vector<IndexedDoubleVector> XBValues;
0020 
0021   // read file and copy to a vector
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;  //,poten;
0028       stringstream linestr;
0029       linestr << line;
0030       linestr >> x1 >> x2 >> x3 >> Bx >> By >> Bz >> perm;
0031       ;  // >> poten;
0032       if (file) {
0033         XBVector.putV6(x1, x2, x3, Bx, By, Bz);
0034         //  if (nLines<5) {
0035         //    double tx1, tx2, tx3, tBx, tBy, tBz;
0036         //    XBVector.getV6(tx1, tx2, tx3, tBx, tBy, tBz);
0037         //    cout << "Read: " <<  setprecision(12) << Bx << " " << tBx << endl;
0038         //  }
0039         XBValues.push_back(XBVector);
0040         ++nLines;
0041       }
0042     }
0043   }
0044 
0045   // compare point by point
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   // check, whether coordinates are Cartesian or cylindrical
0093   std::string::size_type ibeg, iend;
0094   ibeg = name.find('-');   // first occurance of "-"
0095   iend = name.rfind('-');  // last  occurance of "-"
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   // define vectors of IndexedDoubleVectors
0106   IndexedDoubleVector XBVector;
0107   std::vector<IndexedDoubleVector> XBValues;
0108   std::vector<IndexedDoubleVector> IXBValues;
0109   std::vector<IndexedDoubleVector> XBArray;
0110 
0111   // vectors for pre analysis
0112   std::vector<double> valX[3];
0113   std::vector<int> numX[3];
0114 
0115   // copy ASCII file to the vector
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         //  cerr << fixed << x1 << " " <<  x2 << " " <<  x3 << " " <<  Bx << " " <<  By << " " <<  Bz << " " <<  perm << " " <<  poten << " " << endl;
0153 
0154         // Check that the table is homogeneous (ie does not mix perm=1 with perm>>1)
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         // pre analyze file content
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   // PART I: FOR SIMPLE COORDINATES
0196 
0197   // begin simple grid (EasyCoordinate) analysis
0198   for (int i = 0; i < 3; ++i) {
0199     // sort the different values per one coordinate
0200     sort(valX[i].begin(), valX[i].end());
0201     // calculate number of points and constant step size
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   // define indices for easy coordinates
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   // sort according to defined indices
0240   sort(IXBValues.begin(), IXBValues.end());
0241 
0242   // check, if structure is known at this point
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   // get first point = ReferencePoint (and last point for structure checking)
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   // recheck all coordinates
0271   KnownStructure = true;
0272   for (int i = 0; i < 3; ++i) {
0273     ReferencePoint[i] = firstListPoint[i];
0274     // calculate step size for maximal indices
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     // calculate difference between the two reference points minus the total difference (=0)
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       // print result
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   // END OF PART I: FOR SIMPLE COORDINATES
0319   //
0320   // PART II: FOR MORE COMPLICATED COORDINATES
0321   bool goodIndex[3] = {true, true, true};
0322 
0323   if (!KnownStructure) {
0324     // analyze structure of missing coordinates
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     // try to recover info of missing coordiates
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     // try to recover info of last missing coordiates
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       // reset wrong indices and sort
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     // common part of missing coordiate(s) recovery (1 or 2 missing)
0480     if (nEasy > 0) {
0481       // sort without resetting of indices
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       // basic distances (sorry, I have found no smarter solution so far)
0499       // recalculate constant terms for step size from miniGrid[3][2][2][2]
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       // now calculate linear terms for step size from miniGrid[3][2][2][2]
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       // now calculate linear terms offsets from miniGrid[3][2][2][2]
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       // set default values of BasicDistance0 in case of < 2 points
0533       for (int i = 0; i < 3; ++i) {
0534         if (NumberOfPoints[i] < 2)
0535           BasicDistance0[i] = 9e9;
0536       }
0537     }
0538 
0539     // get first point = ReferencePoint (and last point for structure checking)
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     // recheck all coordinates
0554     KnownStructure = true;
0555     for (int i = 0; i < 3; ++i) {
0556       ReferencePoint[i] = firstListPoint[i];
0557       nSteps[i] = NumberOfPoints[i] - 1;
0558       // calculate step size for maximal indices
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       // calculate difference between the two reference points minus the total difference (=0)
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         // print result
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   // END OF PART II: FOR MORE COMPLICATED COORDINATES
0627 
0628   // save coordinates and field values
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   // GridType = 5; // 1/sin(phi) coordinates. Obsolete, not used anymore.
0644   GridType = 6;  // 1/cos(phi) coordinates.
0645 
0646   double phi = (rotateSector - 1) * Geom::pi() / 6.;
0647   double sphi = sin(-phi);
0648   double cphi = cos(-phi);
0649 
0650   // check, whether coordinates are Cartesian or cylindrical
0651   std::string::size_type ibeg, iend;
0652   ibeg = name.find('-');   // first occurance of "-"
0653   iend = name.rfind('-');  // last  occurance of "-"
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   // define vectors of IndexedDoubleVectors
0664   IndexedDoubleVector XBVector;
0665   std::vector<IndexedDoubleVector> XBValues;
0666   std::vector<IndexedDoubleVector> IXBValues;
0667   std::vector<IndexedDoubleVector> XBArray;
0668 
0669   // vectors for pre analysis
0670   std::vector<double> valX[3];
0671   std::vector<int> numX[3];
0672 
0673   // copy ASCII file to the vector
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         // pre analyze file content
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   // begin simple grid (EasyCoordinate) analysis
0741   for (int i = 0; i < 3; ++i) {
0742     // sort the different values per one coordinate
0743     sort(valX[i].begin(), valX[i].end());
0744     // calculate number of points and constant step size
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   // define indices for easy coordinates
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   // sort according to defined indices
0785   sort(IXBValues.begin(), IXBValues.end());
0786 
0787   // check, if structure is known at this point
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   // get first point = ReferencePoint (and last point for structure checking)
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   // recheck all coordinates
0816   KnownStructure = true;
0817   for (int i = 0; i < 3; ++i) {
0818     ReferencePoint[i] = firstListPoint[i];
0819     // calculate step size for maximal indices
0820     double stepSize = BasicDistance0[i];
0821     // calculate difference between the two reference points minus the total difference (=0)
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     // analyze structure of missing coordinates
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     // try to recover info of missing coordiates
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     int nGoodInd = 0;
0964     for (int i = 0; i < 3; ++i) {
0965       if (goodIndex[i])
0966         ++nGoodInd;
0967     }
0968     // common part of missing coordiate(s) recovery (1 or 2 missing)
0969     if (nEasy > 0) {
0970       // sort without resetting of indices
0971       sort(XBArray.begin(), XBArray.end());
0972 
0973       for (int i = 0; i < 3; ++i) {
0974         nSteps[i] = NumberOfPoints[i] - 1;
0975       }
0976       std::vector<double> rMinVec;
0977       std::vector<double> rMaxVec;
0978       std::vector<double> phiVec;
0979       for (int iLine = 0; iLine < nLines; ++iLine) {
0980         XBVector = XBArray.operator[](iLine);
0981         XBVector.getI3(index[0], index[1], index[2]);
0982         // in this case:r,phi,z
0983         XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0984         if (index[2] == 0) {
0985           if (index[0] == 0)
0986             rMinVec.push_back(x1);
0987           if (index[0] == nSteps[0])
0988             rMaxVec.push_back(x1);
0989           if (index[0] == nSteps[0])
0990             phiVec.push_back(x2);
0991         }
0992       }
0993       for (int j = 0; j < NumberOfPoints[1]; ++j) {
0994         double phi = phiVec.operator[](j);
0995         double rMin = rMinVec.operator[](j);
0996         double rMax = rMaxVec.operator[](j);
0997         double rhoMin, rhoMax;
0998         if (GridType == 6) {
0999           rhoMin = rMin * cos(phi);
1000           rhoMax = rMax * cos(phi);
1001         } else if (GridType == 5) {
1002           rhoMin = rMin * sin(phi);
1003           rhoMax = rMax * sin(phi);
1004         } else {
1005           cout << "unknown grid type!" << endl;
1006           abort();
1007         }
1008 
1009         if (j == 0) {
1010           RParAsFunOfPhi[0] = rMax;
1011           RParAsFunOfPhi[1] = rhoMax;
1012           RParAsFunOfPhi[2] = rMin;
1013           RParAsFunOfPhi[3] = rhoMin;
1014         } else {
1015           if (std::abs(rMax - RParAsFunOfPhi[0]) > EPSILON)
1016             RParAsFunOfPhi[0] = 0.;
1017           if (std::abs(rhoMax - RParAsFunOfPhi[1]) > EPSILON)
1018             RParAsFunOfPhi[1] = 0.;
1019           if (std::abs(rMin - RParAsFunOfPhi[2]) > EPSILON)
1020             RParAsFunOfPhi[2] = 0.;
1021           if (std::abs(rhoMin - RParAsFunOfPhi[3]) > EPSILON)
1022             RParAsFunOfPhi[3] = 0.;
1023         }
1024       }
1025 
1026       // set default values of BasicDistance0 in case of < 2 points
1027       for (int i = 0; i < 3; ++i) {
1028         if (NumberOfPoints[i] < 2)
1029           BasicDistance0[i] = 9e9;
1030       }
1031     }
1032 
1033     // get first point = ReferencePoint (and last point for structure checking)
1034     if (!XBArray.empty())
1035       XBVector = XBArray.front();
1036     XBVector.getV6(x1, x2, x3, Bx, By, Bz);
1037     firstListPoint[0] = x1;
1038     firstListPoint[1] = x2;
1039     firstListPoint[2] = x3;
1040     if (!XBArray.empty())
1041       XBVector = XBArray.back();
1042     XBVector.getV6(x1, x2, x3, Bx, By, Bz);
1043     secondRefPoint[0] = x1;
1044     secondRefPoint[1] = x2;
1045     secondRefPoint[2] = x3;
1046 
1047     // recheck all coordinates
1048     KnownStructure = true;
1049     for (int i = 0; i < 3; ++i) {
1050       ReferencePoint[i] = firstListPoint[i];
1051       nSteps[i] = NumberOfPoints[i] - 1;
1052     }
1053 
1054     double sinPhi;  // either sin or cos depending on grid type
1055     if (GridType == 6) {
1056       sinPhi = cos(secondRefPoint[1]);
1057     } else if (GridType == 5) {
1058       sinPhi = sin(secondRefPoint[1]);
1059     } else {
1060       cout << "unknown GridType!" << endl;
1061       abort();
1062     }
1063 
1064     if (std::abs(sinPhi) < EPSILON) {
1065       cout << " WARNING: unexpected sinPhi parameter = 0" << endl;
1066       sinPhi = EPSILON;
1067     }
1068 
1069     double totStepSize =
1070         RParAsFunOfPhi[0] + RParAsFunOfPhi[1] / sinPhi - RParAsFunOfPhi[2] - RParAsFunOfPhi[3] / sinPhi;
1071     double startingPoint = RParAsFunOfPhi[2] + RParAsFunOfPhi[3] / sinPhi;
1072     double totDiff = secondRefPoint[0] - (startingPoint + totStepSize);
1073     if (std::abs(totDiff) < EPSILON * 10.)
1074       systematicGrid[0] = true;
1075     else
1076       systematicGrid[0] = false;
1077     if (!systematicGrid[0])
1078       KnownStructure = false;
1079     totDiff = secondRefPoint[1] - (ReferencePoint[1] + BasicDistance0[1] * nSteps[1]);
1080     if (std::abs(totDiff) < EPSILON * 10.)
1081       systematicGrid[1] = true;
1082     else
1083       systematicGrid[1] = false;
1084     if (!systematicGrid[1])
1085       KnownStructure = false;
1086     totDiff = secondRefPoint[2] - (ReferencePoint[2] + BasicDistance0[2] * nSteps[2]);
1087     if (std::abs(totDiff) < EPSILON * 10.)
1088       systematicGrid[2] = true;
1089     else
1090       systematicGrid[2] = false;
1091     if (!systematicGrid[2])
1092       KnownStructure = false;
1093 
1094     if (PRINT) {
1095       // print result
1096       cout << "  read " << nLines << " lines -> grid structure:" << endl;
1097       cout << "  # of points:   N1 = " << NumberOfPoints[0] << "   N2 = " << NumberOfPoints[1]
1098            << "   N3 = " << NumberOfPoints[2] << endl;
1099       cout << "  ref. point :   X1 = " << ReferencePoint[0] << "   X2 = " << ReferencePoint[1]
1100            << "   X3 = " << ReferencePoint[2] << endl;
1101       cout << "  step parm0 :   A1 = " << BasicDistance0[0] << "   A2 = " << BasicDistance0[1]
1102            << "   A3 = " << BasicDistance0[2] << endl;
1103       cout << "  r param.s. :   R1 = " << RParAsFunOfPhi[0] << " RHO1 = " << RParAsFunOfPhi[1]
1104            << "   R2 = " << RParAsFunOfPhi[2] << " RHO2 = " << RParAsFunOfPhi[3] << endl;
1105       cout << "  easy grid  :   E1 = " << EasyCoordinate[0] << "   E2 = " << EasyCoordinate[1]
1106            << "   E3 = " << EasyCoordinate[2] << endl;
1107       cout << "             :   I1 = " << goodIndex[0] << "   I2 = " << goodIndex[1] << "   I3 = " << goodIndex[2]
1108            << endl;
1109       cout << "  structure  :   S1 = " << systematicGrid[0] << "   S2 = " << systematicGrid[1]
1110            << "   S3 = " << systematicGrid[2];
1111       if (KnownStructure)
1112         cout << "  -->  VALID   (step II)" << endl;
1113       else {
1114         cout << "  -->  INVALID (step II)" << endl;
1115         cout << "  reason for error: ";
1116         if (NumberOfPoints[0] * NumberOfPoints[1] * NumberOfPoints[2] != nLines) {
1117           cout << endl;
1118           cout << "  N1*N2*N3 =/= N.lines  -->  exiting now ..." << endl;
1119           return;
1120         } else {
1121           cout << "  no idea so far  -->  exiting now ..." << endl;
1122           return;
1123         }
1124       }
1125     }
1126   }
1127 
1128   // save coordinates and field values
1129   SixDPoint GridEntry;
1130   for (unsigned int iLine = 0; iLine < XBArray.size(); ++iLine) {
1131     XBVector = XBArray.operator[](iLine);
1132     XBVector.getV6(x1, x2, x3, Bx, By, Bz);
1133     GridEntry.putP6(x1, x2, x3, Bx, By, Bz);
1134     GridData.push_back(GridEntry);
1135   }
1136   if (!XBArray.empty())
1137     XBArray.clear();
1138 
1139   return;
1140 }
1141 
1142 int prepareMagneticFieldGrid::gridType() {
1143   int type = GridType;
1144   if (PRINT) {
1145     if (type == 0)
1146       cout << "  grid type = " << type << "  -->  not determined" << endl;
1147     if (type == 1)
1148       cout << "  grid type = " << type << "  -->  (x,y,z) cube" << endl;
1149     if (type == 2)
1150       cout << "  grid type = " << type << "  -->  (x,y,z) trapezoid" << endl;
1151     if (type == 3)
1152       cout << "  grid type = " << type << "  -->  (r,phi,z) cube" << endl;
1153     if (type == 4)
1154       cout << "  grid type = " << type << "  -->  (r,phi,z) trapezoid" << endl;
1155     if (type == 5)
1156       cout << "  grid type = " << type << "  -->  (r,phi,z) 1/sin(phi)" << endl;
1157     if (type == 6)
1158       cout << "  grid type = " << type << "  -->  (r,phi,z) 1/cos(phi)" << endl;
1159   }
1160   return type;
1161 }
1162 
1163 void prepareMagneticFieldGrid::validateAllPoints() {
1164   if (!KnownStructure)
1165     return;
1166 
1167   double x1, x2, x3, Bx, By, Bz;
1168   int numberOfErrors = 0;
1169 
1170   // loop over three dimensions
1171   int index[3];
1172 
1173   if (GridType < 5) {
1174     for (index[0] = 0; index[0] < NumberOfPoints[0]; ++index[0]) {
1175       for (index[1] = 0; index[1] < NumberOfPoints[1]; ++index[1]) {
1176         for (index[2] = 0; index[2] < NumberOfPoints[2]; ++index[2]) {
1177           // get reference values
1178           putIndicesGetXAndB(index[0], index[1], index[2], x1, x2, x3, Bx, By, Bz);
1179 
1180           // first check: calculate indices from coordinates and compare with original indices
1181           double pnt[3] = {x1, x2, x3};
1182           int tmpIdx[3] = {999, 999, 999};
1183           putCoordGetIndices((x1 + EPSILON), (x2 + EPSILON), (x3 + EPSILON), tmpIdx[0], tmpIdx[1], tmpIdx[2]);
1184           for (int i = 0; i < 3; ++i) {
1185             if (tmpIdx[i] != index[i])
1186               ++numberOfErrors;
1187           }
1188 
1189           // second check: calculate coordinate from indices and compare with original values
1190           double tmpPnt[3] = {0.0, 0.0, 0.0};
1191           putIndCalcXReturnB(index[0], index[1], index[2], tmpPnt[0], tmpPnt[1], tmpPnt[2], Bx, By, Bz);
1192           for (int i = 0; i < 3; ++i) {
1193             if (std::abs(tmpPnt[i] - pnt[i]) > EPSILON)
1194               ++numberOfErrors;
1195           }
1196         }
1197       }
1198     }
1199   }
1200 
1201   if (GridType == 5 || GridType == 6) {
1202     for (index[0] = 0; index[0] < NumberOfPoints[0]; ++index[0]) {
1203       for (index[1] = 0; index[1] < NumberOfPoints[1]; ++index[1]) {
1204         for (index[2] = 0; index[2] < NumberOfPoints[2]; ++index[2]) {
1205           // get reference values
1206           putIndicesGetXAndB(index[0], index[1], index[2], x1, x2, x3, Bx, By, Bz);
1207 
1208           // first check: calculate indices from coordinates and compare with original indices
1209           double pnt[3] = {x1, x2, x3};
1210           int tmpIdx[3] = {999, 999, 999};
1211           putCoordGetIndices((x1 + EPSILON), (x2 + EPSILON), (x3 + EPSILON), tmpIdx[0], tmpIdx[1], tmpIdx[2]);
1212           for (int i = 0; i < 3; ++i) {
1213             if (tmpIdx[i] != index[i]) {
1214               putCoordGetIndices((x1 - EPSILON), (x2 + EPSILON), (x3 + EPSILON), tmpIdx[0], tmpIdx[1], tmpIdx[2]);
1215               if (tmpIdx[i] != index[i]) {
1216                 putCoordGetIndices((x1 + EPSILON), (x2 - EPSILON), (x3 + EPSILON), tmpIdx[0], tmpIdx[1], tmpIdx[2]);
1217                 if (tmpIdx[i] != index[i]) {
1218                   putCoordGetIndices((x1 - EPSILON), (x2 - EPSILON), (x3 + EPSILON), tmpIdx[0], tmpIdx[1], tmpIdx[2]);
1219                   if (tmpIdx[i] != index[i])
1220                     ++numberOfErrors;
1221                 }
1222               }
1223             }
1224           }
1225 
1226           // second check: calculate coordinate from indices and compare with original values
1227           double tmpPnt[3] = {0.0, 0.0, 0.0};
1228           putIndCalcXReturnB(index[0], index[1], index[2], tmpPnt[0], tmpPnt[1], tmpPnt[2], Bx, By, Bz);
1229           for (int i = 0; i < 3; ++i) {
1230             if (std::abs(tmpPnt[i] - pnt[i]) > EPSILON)
1231               ++numberOfErrors;
1232           }
1233         }
1234       }
1235     }
1236   }
1237 
1238   if (numberOfErrors > 0)
1239     KnownStructure = false;
1240   if (PRINT)
1241     cout << "  grid validation of all points done  -->  # errors = " << numberOfErrors << endl;
1242 
1243   return;
1244 }
1245 
1246 void prepareMagneticFieldGrid::saveGridToFile(const std::string &outName) {
1247   // open output file
1248   binary_ofstream outFile(outName);
1249   // write grid type
1250   outFile << GridType;
1251   // write header (depending on grid type)
1252   convertUnits();  // m->cm
1253   if (GridType == 1) {
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   }
1258   if (GridType == 2) {
1259     outFile << NumberOfPoints[0] << NumberOfPoints[1] << NumberOfPoints[2];
1260     outFile << ReferencePoint[0] << ReferencePoint[1] << ReferencePoint[2];
1261     outFile << BasicDistance0[0] << BasicDistance0[1] << BasicDistance0[2];
1262     outFile << BasicDistance1[0][0] << BasicDistance1[1][0] << BasicDistance1[2][0];
1263     outFile << BasicDistance1[0][1] << BasicDistance1[1][1] << BasicDistance1[2][1];
1264     outFile << BasicDistance1[0][2] << BasicDistance1[1][2] << BasicDistance1[2][2];
1265     outFile << BasicDistance2[0][0] << BasicDistance2[1][0] << BasicDistance2[2][0];
1266     outFile << BasicDistance2[0][1] << BasicDistance2[1][1] << BasicDistance2[2][1];
1267     outFile << BasicDistance2[0][2] << BasicDistance2[1][2] << BasicDistance2[2][2];
1268     outFile << EasyCoordinate[0] << EasyCoordinate[1] << EasyCoordinate[2];
1269   }
1270   if (GridType == 3) {
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   }
1275   if (GridType == 4) {
1276     outFile << NumberOfPoints[0] << NumberOfPoints[1] << NumberOfPoints[2];
1277     outFile << ReferencePoint[0] << ReferencePoint[1] << ReferencePoint[2];
1278     outFile << BasicDistance0[0] << BasicDistance0[1] << BasicDistance0[2];
1279     outFile << BasicDistance1[0][0] << BasicDistance1[1][0] << BasicDistance1[2][0];
1280     outFile << BasicDistance1[0][1] << BasicDistance1[1][1] << BasicDistance1[2][1];
1281     outFile << BasicDistance1[0][2] << BasicDistance1[1][2] << BasicDistance1[2][2];
1282     outFile << BasicDistance2[0][0] << BasicDistance2[1][0] << BasicDistance2[2][0];
1283     outFile << BasicDistance2[0][1] << BasicDistance2[1][1] << BasicDistance2[2][1];
1284     outFile << BasicDistance2[0][2] << BasicDistance2[1][2] << BasicDistance2[2][2];
1285     outFile << EasyCoordinate[0] << EasyCoordinate[1] << EasyCoordinate[2];
1286   }
1287   if (GridType == 5 || GridType == 6) {
1288     outFile << NumberOfPoints[0] << NumberOfPoints[1] << NumberOfPoints[2];
1289     outFile << ReferencePoint[0] << ReferencePoint[1] << ReferencePoint[2];
1290     outFile << BasicDistance0[0] << BasicDistance0[1] << BasicDistance0[2];
1291     outFile << RParAsFunOfPhi[0] << RParAsFunOfPhi[1] << RParAsFunOfPhi[2] << RParAsFunOfPhi[3];
1292   }
1293   int nLines = NumberOfPoints[0] * NumberOfPoints[1] * NumberOfPoints[2];
1294   SixDPoint GridPoint;
1295   // write magnetic field values (3*nPoints)
1296   for (int iLine = 0; iLine < nLines; ++iLine) {
1297     GridPoint = GridData.operator[](iLine);
1298     float Bx = float(GridPoint.bx());
1299     float By = float(GridPoint.by());
1300     float Bz = float(GridPoint.bz());
1301     outFile << Bx << By << Bz;
1302     //    if (iLine<5) cout << setprecision(12) << Bx << " " << GridPoint.bx() << endl;
1303   }
1304   // make end and close output file
1305   const std::string lastEntry = "complete";
1306   outFile << lastEntry;
1307   outFile.close();
1308   if (PRINT)
1309     cout << "  output " << outName << endl;
1310   return;
1311 }
1312 
1313 void prepareMagneticFieldGrid::convertUnits() {
1314   double cm = 100.;  // m->cm (just multiply all lengths with 100)
1315   if (XyzCoordinates) {
1316     for (int i = 0; i < 3; ++i) {
1317       ReferencePoint[i] *= cm;
1318     };
1319     for (int i = 0; i < 3; ++i) {
1320       BasicDistance0[i] *= cm;
1321     };
1322     for (int i = 0; i < 3; ++i) {
1323       for (int j = 0; j < 3; ++j) {
1324         BasicDistance1[i][j] *= cm;
1325       };
1326     };
1327     for (int i = 0; i < 3; ++i) {
1328       for (int j = 0; j < 3; ++j) {
1329         BasicDistance2[i][j] *= cm;
1330       };
1331     };
1332     for (int i = 0; i < 4; ++i) {
1333       RParAsFunOfPhi[i] *= cm;
1334     };
1335   }
1336   double du[3] = {100., 1., 100.};  // m->cm ; rad->rad (unchanged)
1337   if (RpzCoordinates) {
1338     for (int i = 0; i < 3; ++i) {
1339       ReferencePoint[i] *= du[i];
1340     };
1341     for (int i = 0; i < 3; ++i) {
1342       BasicDistance0[i] *= du[i];
1343     };
1344     for (int i = 0; i < 3; ++i) {
1345       for (int j = 0; j < 3; ++j) {
1346         BasicDistance1[i][j] *= du[i];
1347       };
1348     };
1349     for (int i = 0; i < 3; ++i) {
1350       for (int j = 0; j < 3; ++j) {
1351         BasicDistance2[i][j] *= du[i];
1352       };
1353     };
1354     for (int i = 0; i < 4; ++i) {
1355       RParAsFunOfPhi[i] *= cm;
1356     };
1357   }
1358   return;
1359 }
1360 
1361 bool prepareMagneticFieldGrid::isReady() { return KnownStructure; }
1362 
1363 void prepareMagneticFieldGrid::interpolateAtPoint(double X1, double X2, double X3, double &Bx, double &By, double &Bz) {
1364   Bx = By = Bz = 0.0;
1365   if (KnownStructure) {
1366     // define interpolation object
1367     VectorFieldInterpolation MagInterpol;
1368     // calculate indices for "CellPoint000"
1369     int index[3];
1370     putCoordGetIndices(X1, X2, X3, index[0], index[1], index[2]);
1371     int index0[3] = {0, 0, 0};
1372     int index1[3] = {0, 0, 0};
1373     for (int i = 0; i < 3; ++i) {
1374       if (NumberOfPoints[i] > 1) {
1375         index0[i] = std::max(0, index[i]);
1376         if (index0[i] > NumberOfPoints[i] - 2)
1377           index0[i] = NumberOfPoints[i] - 2;
1378         ;
1379         index1[i] = std::max(1, index[i] + 1);
1380         ;
1381         if (index1[i] > NumberOfPoints[i] - 1)
1382           index1[i] = NumberOfPoints[i] - 1;
1383       }
1384     }
1385     double tmpX[3];
1386     double tmpB[3];
1387     // define the corners of interpolation volume
1388     putIndicesGetXAndB(index0[0], index0[1], index0[2], tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1389     MagInterpol.defineCellPoint000(tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1390     putIndicesGetXAndB(index1[0], index0[1], index0[2], tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1391     MagInterpol.defineCellPoint100(tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1392     putIndicesGetXAndB(index0[0], index1[1], index0[2], tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1393     MagInterpol.defineCellPoint010(tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1394     putIndicesGetXAndB(index1[0], index1[1], index0[2], tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1395     MagInterpol.defineCellPoint110(tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1396     putIndicesGetXAndB(index0[0], index0[1], index1[2], tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1397     MagInterpol.defineCellPoint001(tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1398     putIndicesGetXAndB(index1[0], index0[1], index1[2], tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1399     MagInterpol.defineCellPoint101(tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1400     putIndicesGetXAndB(index0[0], index1[1], index1[2], tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1401     MagInterpol.defineCellPoint011(tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1402     putIndicesGetXAndB(index1[0], index1[1], index1[2], tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1403     MagInterpol.defineCellPoint111(tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1404     // interpolate
1405     MagInterpol.putSCoordGetVField(X1, X2, X3, Bx, By, Bz);
1406   }
1407 
1408   return;
1409 }
1410 
1411 void prepareMagneticFieldGrid::putCoordGetIndices(
1412     double X1, double X2, double X3, int &Index1, int &Index2, int &Index3) {
1413   double pnt[3] = {X1, X2, X3};
1414   int index[3] = {0, 0, 0};
1415 
1416   if (GridType < 5) {
1417     for (int i = 0; i < 3; ++i) {
1418       if (EasyCoordinate[i]) {
1419         index[i] = int((pnt[i] - ReferencePoint[i]) / BasicDistance0[i]);
1420       }
1421     }
1422     for (int i = 0; i < 3; ++i) {
1423       if (!EasyCoordinate[i]) {
1424         double stepSize = BasicDistance0[i];
1425         double offset = 0.0;
1426         for (int j = 0; j < 3; ++j) {
1427           stepSize += BasicDistance1[i][j] * index[j];
1428           offset += BasicDistance2[i][j] * index[j];
1429         }
1430         index[i] = int((pnt[i] - (ReferencePoint[i] + offset)) / stepSize);
1431       }
1432     }
1433   }
1434   if (GridType == 5 || GridType == 6) {
1435     double sinPhi;  // either sin or cos depending on grid type
1436     if (GridType == 6) {
1437       sinPhi = cos(pnt[1]);
1438     } else if (GridType == 5) {
1439       sinPhi = sin(pnt[1]);
1440     } else {
1441       abort();
1442     }
1443 
1444     if (std::abs(sinPhi) < EPSILON) {
1445       sinPhi = EPSILON;
1446       cout << "ERROR DIVISION BY ZERO" << endl;
1447     }
1448     double stepSize = RParAsFunOfPhi[0] + RParAsFunOfPhi[1] / sinPhi - RParAsFunOfPhi[2] - RParAsFunOfPhi[3] / sinPhi;
1449     stepSize = stepSize / (NumberOfPoints[0] - 1);
1450     double startingPoint = RParAsFunOfPhi[2] + RParAsFunOfPhi[3] / sinPhi;
1451     index[0] = int((pnt[0] - startingPoint) / stepSize);
1452     index[1] = int((pnt[1] - ReferencePoint[1]) / BasicDistance0[1]);
1453     index[2] = int((pnt[2] - ReferencePoint[2]) / BasicDistance0[2]);
1454   }
1455 
1456   Index1 = index[0];
1457   Index2 = index[1];
1458   Index3 = index[2];
1459 
1460   return;
1461 }
1462 
1463 void prepareMagneticFieldGrid::putIndicesGetXAndB(
1464     int Index1, int Index2, int Index3, double &X1, double &X2, double &X3, double &Bx, double &By, double &Bz) {
1465   SixDPoint GridPoint;
1466   GridPoint = GridData.operator[](lineNumber(Index1, Index2, Index3));
1467   X1 = GridPoint.x1();
1468   X2 = GridPoint.x2();
1469   X3 = GridPoint.x3();
1470   Bx = GridPoint.bx();
1471   By = GridPoint.by();
1472   Bz = GridPoint.bz();
1473 
1474   return;
1475 }
1476 
1477 void prepareMagneticFieldGrid::putIndCalcXReturnB(
1478     int Index1, int Index2, int Index3, double &X1, double &X2, double &X3, double &Bx, double &By, double &Bz) {
1479   int index[3] = {Index1, Index2, Index3};
1480   double pnt[3] = {0., 0., 0.};
1481 
1482   if (GridType < 5) {
1483     for (int i = 0; i < 3; ++i) {
1484       if (EasyCoordinate[i]) {
1485         pnt[i] = ReferencePoint[i] + BasicDistance0[i] * index[i];
1486       } else {
1487         double stepSize = BasicDistance0[i];
1488         double offset = 0.0;
1489         for (int j = 0; j < 3; ++j) {
1490           stepSize += BasicDistance1[i][j] * index[j];
1491           offset += BasicDistance2[i][j] * index[j];
1492         }
1493         pnt[i] = ReferencePoint[i] + offset + stepSize * index[i];
1494       }
1495     }
1496   }
1497   if (GridType == 5 || GridType == 6) {
1498     pnt[2] = ReferencePoint[2] + BasicDistance0[2] * index[2];
1499     pnt[1] = ReferencePoint[1] + BasicDistance0[1] * index[1];
1500 
1501     double sinPhi;  // Set either cos or sin depending on grid type
1502     if (GridType == 6) {
1503       sinPhi = cos(pnt[1]);
1504     } else if (GridType == 5) {
1505       sinPhi = sin(pnt[1]);
1506     } else {
1507       cout << "unknown GridType!" << endl;
1508       abort();
1509     }
1510 
1511     if (std::abs(sinPhi) < EPSILON) {
1512       sinPhi = EPSILON;
1513       cout << "ERROR DIVISION BY ZERO" << endl;
1514     }
1515     double stepSize = RParAsFunOfPhi[0] + RParAsFunOfPhi[1] / sinPhi - RParAsFunOfPhi[2] - RParAsFunOfPhi[3] / sinPhi;
1516     stepSize = stepSize / (NumberOfPoints[0] - 1);
1517     double startingPoint = RParAsFunOfPhi[2] + RParAsFunOfPhi[3] / sinPhi;
1518     pnt[0] = startingPoint + stepSize * index[0];
1519   }
1520 
1521   X1 = pnt[0];
1522   X2 = pnt[1];
1523   X3 = pnt[2];
1524 
1525   SixDPoint GridPoint;
1526   GridPoint = GridData.operator[](lineNumber(Index1, Index2, Index3));
1527   Bx = GridPoint.bx();
1528   By = GridPoint.by();
1529   Bz = GridPoint.bz();
1530 
1531   return;
1532 }
1533 
1534 int prepareMagneticFieldGrid::lineNumber(int Index1, int Index2, int Index3) {
1535   return Index1 * NumberOfPoints[1] * NumberOfPoints[2] + Index2 * NumberOfPoints[2] + Index3;
1536 }
1537 
1538 bool prepareMagneticFieldGrid::IndexedDoubleVector::operator<(const IndexedDoubleVector &x) const {
1539   if (I3[0] < x.I3[0])
1540     return true;
1541   else if (I3[0] == x.I3[0]) {
1542     if (I3[1] < x.I3[1])
1543       return true;
1544     else if (I3[1] == x.I3[1]) {
1545       if (I3[2] < x.I3[2])
1546         return true;
1547       else
1548         return false;
1549     } else
1550       return false;
1551   } else
1552     return false;
1553 }
1554 
1555 void prepareMagneticFieldGrid::IndexedDoubleVector::putI3(int index1, int index2, int index3) {
1556   I3[0] = index1;
1557   I3[1] = index2;
1558   I3[2] = index3;
1559   return;
1560 }
1561 
1562 void prepareMagneticFieldGrid::IndexedDoubleVector::getI3(int &index1, int &index2, int &index3) {
1563   index1 = I3[0];
1564   index2 = I3[1];
1565   index3 = I3[2];
1566   return;
1567 }
1568 
1569 void prepareMagneticFieldGrid::IndexedDoubleVector::putV6(
1570     double X1, double X2, double X3, double Bx, double By, double Bz) {
1571   V6[0] = X1;
1572   V6[1] = X2;
1573   V6[2] = X3;
1574   V6[3] = Bx;
1575   V6[4] = By;
1576   V6[5] = Bz;
1577   return;
1578 }
1579 
1580 void prepareMagneticFieldGrid::IndexedDoubleVector::getV6(
1581     double &X1, double &X2, double &X3, double &Bx, double &By, double &Bz) {
1582   X1 = V6[0];
1583   X2 = V6[1];
1584   X3 = V6[2];
1585   Bx = V6[3];
1586   By = V6[4];
1587   Bz = V6[5];
1588   return;
1589 }
1590 
1591 void prepareMagneticFieldGrid::SixDPoint::putP6(double X1, double X2, double X3, double Bx, double By, double Bz) {
1592   P6[0] = X1;
1593   P6[1] = X2;
1594   P6[2] = X3;
1595   P6[3] = Bx;
1596   P6[4] = By;
1597   P6[5] = Bz;
1598   return;
1599 }
1600 
1601 double prepareMagneticFieldGrid::SixDPoint::x1() { return P6[0]; }
1602 
1603 double prepareMagneticFieldGrid::SixDPoint::x2() { return P6[1]; }
1604 
1605 double prepareMagneticFieldGrid::SixDPoint::x3() { return P6[2]; }
1606 
1607 double prepareMagneticFieldGrid::SixDPoint::bx() { return P6[3]; }
1608 
1609 double prepareMagneticFieldGrid::SixDPoint::by() { return P6[4]; }
1610 
1611 double prepareMagneticFieldGrid::SixDPoint::bz() { return P6[5]; }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!