#include <BfldHandlerRect2d.h>
Inheritance diagram for BfldHandlerRect2d:
Public Member Functions | |
| BfldHandlerRect2d () | |
| virtual | ~BfldHandlerRect2d () |
| virtual TVector3 | GetBFieldMeshCoord (const TVector3 &xyz) |
| virtual void | SetInterpMethod (const BfldInterpMethod::InterpMethod_t method) |
|
|
Definition at line 32 of file BfldHandlerRect2d.cxx. 00033 : BfldHandler() 00034 { 00035 00036 }
|
|
|
Definition at line 39 of file BfldHandlerRect2d.cxx. 00040 {
00041
00042 }
|
|
|
Implements BfldHandler. Definition at line 81 of file BfldHandlerRect2d.cxx. References BfldInterpMethod::AsString(), BfldMeshRect2d::GeneratorToIndices(), BfldMap::GetBField(), BfldMeshRect2d::GetGeneratorPosition(), BfldMeshRect2d::GetQuadFlag(), BfldMap::GetVariant(), BfldMeshRect2d::InBounds(), BfldMeshRect2d::InQuadrant(), MSG, BfldMeshRect2d::NearestGenerator(), BfldMeshRect2d::Pick3Generators(), and BfldMeshRect2d::Pick4Generators(). 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 }
|
|
|
Reimplemented from BfldHandler. Definition at line 45 of file BfldHandlerRect2d.cxx. References BfldInterpMethod::AsString(), MSG, and BfldHandler::SetInterpMethod(). 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 }
|
1.3.9.1