BfldHandlerRect2d Class Reference

#include <BfldHandlerRect2d.h>

Inheritance diagram for BfldHandlerRect2d:
BfldHandler

List of all members.

Public Member Functions

 BfldHandlerRect2d ()
virtual ~BfldHandlerRect2d ()
virtual TVector3 GetBFieldMeshCoord (const TVector3 &xyz)
virtual void SetInterpMethod (const BfldInterpMethod::InterpMethod_t method)

Detailed Description

Definition at line 17 of file BfldHandlerRect2d.h.


Constructor & Destructor Documentation

BfldHandlerRect2d::BfldHandlerRect2d (  ) 

Definition at line 32 of file BfldHandlerRect2d.cxx.

00033    : BfldHandler()
00034 {
00035 
00036 }

BfldHandlerRect2d::~BfldHandlerRect2d (  )  [virtual]

Definition at line 39 of file BfldHandlerRect2d.cxx.

00040 {
00041 
00042 }


Member Function Documentation

TVector3 BfldHandlerRect2d::GetBFieldMeshCoord ( const TVector3 &  xyz  )  [virtual]

Implements BfldHandler.

Definition at line 81 of file BfldHandlerRect2d.cxx.

References MuELoss::a, bfld::AsString(), BfldHandler::fInterpMethod, BfldHandler::fMap, BfldHandler::fMesh, BfldMeshRect2d::GeneratorToIndices(), BfldMap::GetBField(), BfldMeshRect2d::GetGeneratorPosition(), BfldMeshRect2d::GetQuadFlag(), BfldMap::GetVariant(), BfldMeshRect2d::InBounds(), BfldMeshRect2d::InQuadrant(), BfldInterpMethod::kBilinear, BfldInterpMethod::kClosest, BfldInterpMethod::kPlanar, BfldInterpMethod::kPlanarVec, Msg::kVerbose, Msg::kWarning, MSG, BfldMeshRect2d::NearestGenerator(), BfldMeshRect2d::Pick3Generators(), BfldMeshRect2d::Pick4Generators(), and zeroBfld().

00082 {
00083    // Search for the right location and then interpolate
00084 
00085    Double_t x = xyz.X();
00086    Double_t y = xyz.Y();
00087 
00088    // negative map number implies left-right flip with no
00089    // change in the flux circulation around the coil
00090    if (fMap->GetVariant() < 0) x = -x;
00091 
00092    // cast up so we can use the specific methods for this mesh/map type
00093    BfldMeshRect2d *fMeshRect2d = dynamic_cast<BfldMeshRect2d*>(fMesh);
00094    assert(fMeshRect2d);
00095    BfldMapRect2d *fMapRect2d   = dynamic_cast<BfldMapRect2d*>(fMap);
00096    assert(fMapRect2d);
00097 
00098    // special to these maps ..  out of range ==> zero field
00099    // do NOT extrapolate (since some far fields touch the border)
00100    if ( ! fMeshRect2d->InBounds(x,y) ) return zeroBfld;
00101 
00102    // special case if mesh is only for a quadrant
00103    if ( fMeshRect2d->InQuadrant(x,y) != 0) {
00104      x = TMath::Abs(x);
00105      y = TMath::Abs(y);
00106    }
00107 
00108    // this will form the basis for the return TVector3
00109    Double_t bvec[3] = { 0, 0, 0 };
00110 
00111    Int_t generator;
00112 
00113    Double_t a[3],b[3],c[3],d[3];
00114 
00115    static MsgFormat i8("i8");
00116    static MsgFormat f8p3("f8.3");
00117 
00118    switch (fInterpMethod) {
00119    case BfldInterpMethod::kClosest:
00120      {
00121        //
00122        // Closest/Nearest Neighbor "Interpolation"
00123        //
00124        //MAXMSG("Bfld",Msg::kInfo,10) << "BfldHandlerRect2d kClosest" << endl;
00125        generator = fMeshRect2d->NearestGenerator(x,y);
00126        
00127        // For now no interpolation
00128        if (generator<0) {
00129          // no field outside map region
00130        } else {
00131          // nearest neightbor
00132          TVector3 bnear = fMap->GetBField(generator);
00133          bvec[0] = bnear.X();
00134          bvec[1] = bnear.Y();
00135          bvec[2] = bnear.Z();
00136        }
00137      }
00138      break;
00139 
00140    case BfldInterpMethod::kPlanarVec:
00141      {
00142        //
00143        // Planar Interpolation
00144        // based on linear interpolation of three points
00145        //
00146        //MAXMSG("Bfld",Msg::kInfo,10) << "BfldHandlerRect2d kPlanarVec" << endl;
00147        Int_t gen_nearest, gen_cw, gen_ccw;
00148        fMeshRect2d->Pick3Generators(x,y,gen_nearest,gen_cw,gen_ccw);
00149        
00150        TVector3 xyz1 = fMeshRect2d->GetGeneratorPosition(gen_nearest);
00151        TVector3 xyz2 = fMeshRect2d->GetGeneratorPosition(gen_cw);
00152        TVector3 xyz3 = fMeshRect2d->GetGeneratorPosition(gen_ccw);
00153        
00154        TVector3 b1  = fMap->GetBField(gen_nearest);
00155        TVector3 b2  = fMap->GetBField(gen_cw);
00156        TVector3 b3  = fMap->GetBField(gen_ccw);
00157        
00158        for (Int_t i = 0; i < 3; i++) {
00159          // Fake vectors (x,y,B_i)
00160          TVector3 p1 = xyz1;
00161          TVector3 p2 = xyz2;
00162          TVector3 p3 = xyz3;
00163          // third component is the Bfield component we're interpolating
00164          p1[2] = b1[i];
00165          p2[2] = b2[i];
00166          p3[2] = b3[i];
00167          
00168          // differences used to compute normal vector
00169          TVector3 p2p1 = p2 - p1;
00170          TVector3 p3p1 = p3 - p1;
00171          
00172          TVector3 normal = p2p1.Cross(p3p1);
00173          a[i] = normal[0];
00174          b[i] = normal[1];
00175          c[i] = normal[2];
00176          d[i] = -normal.Dot(p1);
00177 
00178          if ( c[i] == 0 ) c[i] = FLT_EPSILON;
00179          bvec[i] = -(a[i]*x + b[i]*y + d[i]) / c[i];  
00180          
00181 #ifdef VERBOSE
00182          MSG("Bfld",Msg::kVerbose) 
00183            << " Normal[" << i << "]  (" << f8p3(normal[0]) << "," 
00184            << f8p3(normal[1]) << "," << f8p3(normal[2]) << ")" << endl;
00185          MSG("Bfld",Msg::kVerbose) 
00186            << " p2p1[" << i << "]  (" << f8p3(p2p1[0]) << "," 
00187            << f8p3(p2p1[1]) << "," << f8p3(p2p1[2]) << ")" << endl;
00188          MSG("Bfld",Msg::kVerbose) 
00189            << " p3p1[" << i << "]  (" << f8p3(p3p1[0]) << "," 
00190            << f8p3(p3p1[1]) << "," << f8p3(p3p1[2]) << ")" << endl;
00191 #endif         
00192        }
00193 
00194 #ifdef VERBOSE
00195        if (TMath::Abs(bvec[0]) < 100. ) {
00196          // okay value
00197        } else {
00198          // bad value
00199          MSG("Bfld",Msg::kWarning) 
00200            << "BfldHandlerRect2d::GetBField planar interp (x,y)=" << endl
00201            << "(" << f8p3(x) << "," << f8p3(y) << ") yields "
00202            << "(" << f8p3(bvec[0]) << "," << f8p3(bvec[1]) << "," << f8p3(bvec[2]) << ")" << endl
00203            << "  nearest " << i8(gen_nearest) << " (" << f8p3(xyz1[0]) << "," 
00204            << f8p3(xyz1[1]) << ") where B="
00205            << "(" << f8p3(b1[0]) << "," << f8p3(b1[1]) << "," << f8p3(b1[2]) << ")" << endl
00206            << "  cw      " << i8(gen_cw)      << " (" << f8p3(xyz2[0]) << "," 
00207            << f8p3(xyz2[1]) << ") where B="
00208            << "(" << f8p3(b2[0]) << "," << f8p3(b2[1]) << "," << f8p3(b2[2]) << ")" << endl
00209            << "  ccw     " << i8(gen_ccw)     << " (" << f8p3(xyz3[0]) << "," 
00210            << f8p3(xyz3[1]) << ") where B="
00211            << "(" << f8p3(b3[0]) << "," << f8p3(b3[1]) << "," << f8p3(b3[2]) << ")" << endl
00212            << " calculate a " << a[0] << "," << a[1] << "," << a[2] << endl
00213            << " calculate b " << b[0] << "," << b[1] << "," << b[2] << endl
00214            << " calculate c " << c[0] << "," << c[1] << "," << c[2] << endl
00215            << " calculate d " << d[0] << "," << d[1] << "," << d[2] << endl;
00216        }
00217 #endif
00218 
00219      }
00220      break;
00221 
00222    case BfldInterpMethod::kPlanar:
00223      {
00224        //
00225        // Planar Interpolation
00226        // based on linear interpolation of three points
00227        // using "components" rather than TVector3 to do calculation
00228        //
00229        //MAXMSG("Bfld",Msg::kInfo,10) << "BfldHandlerRect2d kPlanar" << endl;
00230        Int_t gen_nearest, gen_cw, gen_ccw;
00231        fMeshRect2d->Pick3Generators(x,y,gen_nearest,gen_cw,gen_ccw);
00232        
00233        TVector3 xyz1 = fMeshRect2d->GetGeneratorPosition(gen_nearest);
00234        TVector3 xyz2 = fMeshRect2d->GetGeneratorPosition(gen_cw);
00235        TVector3 xyz3 = fMeshRect2d->GetGeneratorPosition(gen_ccw);
00236        
00237        TVector3 b1  = fMap->GetBField(gen_nearest);
00238        TVector3 b2  = fMap->GetBField(gen_cw);
00239        TVector3 b3  = fMap->GetBField(gen_ccw);
00240        
00241        for (Int_t i = 0; i < 3; i++) {
00242          a[i] =   (xyz2.Y()-xyz1.Y()) * (b3[i]-b1[i])
00243            - (xyz3.Y()-xyz1.Y()) * (b2[i]-b1[i]);
00244          
00245          b[i] =   (xyz3.X()-xyz1.X()) * (b2[i]-b1[i])
00246            - (xyz2.X()-xyz1.X()) * (b3[i]-b1[i]);
00247          
00248          c[i] =   (xyz2.X()-xyz1.X()) * (xyz3.Y()-xyz1.Y())
00249            - (xyz3.X()-xyz1.X()) * (xyz2.Y()-xyz1.Y());
00250          
00251          Double_t mag = TMath::Sqrt(a[i]*a[i] + b[i]*b[i] + c[i]*c[i]);
00252          a[i] = a[i]/mag;
00253          b[i] = b[i]/mag;
00254          c[i] = c[i]/mag;
00255          
00256          d[i] = - (xyz1.X()*a[i] + xyz1.Y()*b[i] + b1[i]*c[i]);
00257          
00258          bvec[i] = -(a[i]*x + b[i]*y + d[i]) / c[i];
00259          
00260        }
00261 
00262 #ifdef VERBOSE
00263        static MsgFormat i8("i8");
00264        static MsgFormat f8p3("f8.3");
00265        MSG("Bfld",Msg::kVerbose) 
00266          << "BfldHandlerRect2d::GetBField planar interp " << endl
00267          << "(" << f8p3(x) << "," << f8p3(y) << ") yields " 
00268          << "(" << f8p3(bvec[0]) << "," << f8p3(bvec[1]) << "," << f8p3(bvec[2]) << ")" << endl
00269          << "  nearest " << i8(gen_nearest) << " (" << f8p3(xyz1[0]) << "," 
00270          << f8p3(xyz1[1]) << ") where B="
00271          << "(" << f8p3(b1[0]) << "," << f8p3(b1[1]) << "," << f8p3(b1[2]) << ")" << endl
00272          << "  cw      " << i8(gen_cw)      << " (" << f8p3(xyz2[0]) << "," 
00273          << f8p3(xyz2[1]) << ") where B="
00274          << "(" << f8p3(b2[0]) << "," << f8p3(b2[1]) << "," << f8p3(b2[2]) << ")" << endl
00275          << "  ccw     " << i8(gen_ccw)     << " (" << f8p3(xyz3[0]) << "," 
00276          << f8p3(xyz3[1]) << ") where B="
00277          << "(" << f8p3(b3[0]) << "," << f8p3(b3[1]) << "," << f8p3(b3[2]) << ")" << endl
00278          << " calculate a " << a[0] << "," << a[1] << "," << a[2] << endl
00279          << " calculate b " << b[0] << "," << b[1] << "," << b[2] << endl
00280          << " calculate c " << c[0] << "," << c[1] << "," << c[2] << endl
00281          << " calculate d " << d[0] << "," << d[1] << "," << d[2] << endl;
00282 #endif
00283      }
00284      break;
00285 
00286    case BfldInterpMethod::kBilinear:
00287      {
00288        //
00289        // Bi-linear Interpolation
00290        //
00291        //MAXMSG("Bfld",Msg::kInfo,10) << "BfldHandlerRect2d kBilinear" << endl;
00292        Int_t gen_ll, gen_lr, gen_ur, gen_ul;
00293        fMeshRect2d->Pick4Generators(x,y,gen_ll,gen_lr,gen_ur,gen_ul);
00294       
00295 #ifdef DEBUG_RWH 
00296        // RWH
00297        int ix, iy;
00298        fMeshRect2d->GeneratorToIndices(gen_ll,ix,iy);
00299        cout << " x " << setw(10) << setprecision(6) << x << " ix " << ix 
00300             << " y " << setw(10) << setprecision(6) << y << " iy " << iy 
00301             << resetiosflags(ios::floatfield)  // reset to "%g" format
00302             << endl;
00303 #endif
00304 
00305        TVector3 xyz_ll = fMeshRect2d->GetGeneratorPosition(gen_ll);
00306        //     TVector3 xyz_lr = fMeshRect2d->GetGeneratorPosition(gen_lr);
00307        TVector3 xyz_ur = fMeshRect2d->GetGeneratorPosition(gen_ur);
00308        //     TVector3 xyz_ul = fMeshRect2d->GetGeneratorPosition(gen_ul);
00309        
00310        TVector3 b_ll  = fMap->GetBField(gen_ll);
00311        TVector3 b_lr  = fMap->GetBField(gen_lr);
00312        TVector3 b_ur  = fMap->GetBField(gen_ur);
00313        TVector3 b_ul  = fMap->GetBField(gen_ul);
00314        
00315        Double_t dx = xyz_ur.X() - xyz_ll.X();
00316        Double_t dy = xyz_ur.Y() - xyz_ll.Y();
00317 
00318        if ( 0.0 == dx || 0.0 == dy ) {
00319          static Int_t nwarn = 10;
00320          if ( nwarn > 0) {
00321            nwarn--;
00322            MSG("Bfld",Msg::kWarning)
00323              << "BfldHandlerRect2d::GetBField bilinear interp "
00324              << " dx = " << dx << " dy = " << dy
00325              << " generator indx ll " << gen_ll << " ur " << gen_ur << endl
00326              << " for x = " << x << " y = " << y
00327              << ((nwarn==0) ? "\n   ....last message of this type." : " ")
00328              << endl;
00329          }
00330          if ( 0.0 == dx ) dx = 1.0;
00331          if ( 0.0 == dy ) dy = 1.0;
00332        }
00333                                            
00334        Double_t dx1 = (x - xyz_ll.X())/dx;
00335        Double_t dx2 = (xyz_ur.X() - x)/dx;
00336        
00337        Double_t dy1 = (y - xyz_ll.Y())/dy;
00338        Double_t dy2 = (xyz_ur.Y() - y)/dy;
00339        
00340 #ifdef DEBUG_RWH
00341        // RWH
00342        cout << " dx " << dx1 << " " << dx2
00343             << " dy " << dy1 << " " << dy2 << endl;
00344        cout << " LL bx " << b_ll.x() << " by " << b_ll.y() << endl;
00345        cout << " LR bx " << b_lr.x() << " by " << b_lr.y() << endl;
00346        cout << " UL bx " << b_ul.x() << " by " << b_ul.y() << endl;
00347        cout << " UR bx " << b_ur.x() << " by " << b_ur.y() << endl;
00348 #endif
00349 
00350        Double_t a = dx2*dy2;
00351        Double_t b = dx1*dy2;
00352        Double_t c = dx2*dy1;
00353        Double_t d = dx1*dy1;
00354        bvec[0] = a*b_ll.x()+b*b_lr.x()+c*b_ul.x()+d*b_ur.x();
00355        bvec[1] = a*b_ll.y()+b*b_lr.y()+c*b_ul.y()+d*b_ur.y();
00356        bvec[2] = a*b_ll.z()+b*b_lr.z()+c*b_ul.z()+d*b_ur.z();
00357 
00358 #ifdef DEBUG_RWH
00359        // RWH
00360        cout << " b result bx " << bvec[0] << " by " << bvec[1] << endl;
00361 #endif
00362      }
00363      break;
00364 
00365    default:
00366       MSG("Bfld",Msg::kWarning)
00367          << "BfldHandlerRect2d::GetBField " << endl
00368          << "       unimplemented interpolation method: "
00369          << BfldInterpMethod::AsString(fInterpMethod) << endl;
00370    }
00371 
00372    // handle quadrant symmetry
00373    //     1 | 0
00374    //   ----+----   + centered at (x0,y0)
00375    //     2 | 3
00376    //
00377    // flip signs and *then* swap (order matters)
00378    Int_t iflag = fMeshRect2d->GetQuadFlag(x,y);
00379    if (iflag & 1) bvec[0] = -bvec[0];
00380    if (iflag & 2) bvec[1] = -bvec[1];
00381    if (iflag & 4) {
00382      Double_t temp = bvec[0];
00383      bvec[0] = bvec[1];
00384      bvec[1] = temp;
00385    }
00386 
00387    // negative map number implies left-right flip with no
00388    // change in the flux circulation around the coil
00389    if (fMap->GetVariant() < 0) {
00390       //bvec[0] = +bvec[0];
00391       bvec[1] = -bvec[1];
00392       bvec[2] = -bvec[2];
00393    }
00394 
00395    return TVector3(bvec[0],bvec[1],bvec[2]);
00396 }

void BfldHandlerRect2d::SetInterpMethod ( const BfldInterpMethod::InterpMethod_t  method  )  [virtual]

Reimplemented from BfldHandler.

Definition at line 45 of file BfldHandlerRect2d.cxx.

References bfld::AsString(), BfldInterpMethod::kBilinear, BfldInterpMethod::kDefault, BfldInterpMethod::kNatural, Msg::kWarning, and MSG.

Referenced by BfldCanvasRect2d::BfldCanvasRect2d().

00046 {
00047    // set the interpolation method for this handler
00048    // the Rect2d specialization translates "kDefault" to "kBilinear"
00049    // 
00050    // and should also disallow illegal choices
00051 
00052    static Int_t nwarn = 10;
00053 
00054    BfldInterpMethod::InterpMethod_t chosen = method;
00055 
00056    switch (chosen) {
00057    case BfldInterpMethod::kNatural:
00058       // can't do these methods
00059       if (nwarn > 0) {
00060          nwarn--;
00061          MSG("Bfld",Msg::kWarning)
00062             << "BfldHandlerRect2d::SetInterpMethod '"
00063             << BfldInterpMethod::AsString(chosen) 
00064             << "' not supported"<< endl;
00065       }
00066       // fall through to our choice of default
00067    case BfldInterpMethod::kDefault:
00068       // pick a default method
00069       chosen = BfldInterpMethod::kBilinear;
00070       break;
00071    default:
00072       // no need to adjust anything for most methods
00073       break;
00074    }
00075    
00076    BfldHandler::SetInterpMethod(chosen); // do the normal thing
00077    
00078 }


The documentation for this class was generated from the following files:

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1