Back to home page

Project CMSSW displayed by LXR
 
 

    

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   // 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     // common part of missing coordiate(s) recovery (1 or 2 missing)
0964     if (nEasy > 0) {
0965       // sort without resetting of indices
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         // in this case:r,phi,z
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       // set default values of BasicDistance0 in case of < 2 points
1022       for (int i = 0; i < 3; ++i) {
1023         if (NumberOfPoints[i] < 2)
1024           BasicDistance0[i] = 9e9;
1025       }
1026     }
1027 
1028     // get first point = ReferencePoint (and last point for structure checking)
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     // recheck all coordinates
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;  // either sin or cos depending on grid type
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       // print result
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   // save coordinates and field values
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   // loop over three dimensions
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           // get reference values
1173           putIndicesGetXAndB(index[0], index[1], index[2], x1, x2, x3, Bx, By, Bz);
1174 
1175           // first check: calculate indices from coordinates and compare with original indices
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           // second check: calculate coordinate from indices and compare with original values
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           // get reference values
1201           putIndicesGetXAndB(index[0], index[1], index[2], x1, x2, x3, Bx, By, Bz);
1202 
1203           // first check: calculate indices from coordinates and compare with original indices
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           // second check: calculate coordinate from indices and compare with original values
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   // open output file
1243   binary_ofstream outFile(outName);
1244   // write grid type
1245   outFile << GridType;
1246   // write header (depending on grid type)
1247   convertUnits();  // m->cm
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   // write magnetic field values (3*nPoints)
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     //    if (iLine<5) cout << setprecision(12) << Bx << " " << GridPoint.bx() << endl;
1298   }
1299   // make end and close output file
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.;  // m->cm (just multiply all lengths with 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.};  // m->cm ; rad->rad (unchanged)
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     // define interpolation object
1362     VectorFieldInterpolation MagInterpol;
1363     // calculate indices for "CellPoint000"
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     // define the corners of interpolation volume
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     // interpolate
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;  // either sin or cos depending on grid type
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;  // Set either cos or sin depending on grid type
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]; }
[ 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!