CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

eulerTest.cc
Go to the documentation of this file.
1 // eulerTest.cc
2 
3 // Test extreme cases for Euler angles --
4 // We perturb the matrix slightly before taking Euler angles.
5 // A test will have failed if the Rotation resulting from taking euler angles
6 // and forming a rotation from them, is ever significantly different than the
7 // origninal rotation.
8 
9 
10 #include "CLHEP/Vector/Rotation.h"
11 #include "CLHEP/Vector/EulerAngles.h"
12 #include "CLHEP/Units/PhysicalConstants.h"
13 
14 #include <iostream>
15 #include <math.h>
16 
17 using std::cout;
18 using namespace CLHEP;
19 
20 const int Nperturbs = 24;
21 
22 void perturb ( int i, const HepRotation & r, HepRotation & rp, double & del );
23 bool compareR ( const HepRotation & r1, const HepRotation & r2, double tol );
24 
25 bool test (double phi, double theta, double psi) {
26  HepRotation r (phi, theta, psi);
27  HepRotation rp;
28  HepRotation rpe;
30  bool retval = true;
31  double del;
32  cout << "\n\n -------------------------------- \n\n";
33  for ( int i = 0; i < Nperturbs; ++i ) {
34  perturb ( i, r, rp, del );
35  e = rp.eulerAngles();
36  cout << "(" << phi <<", " << theta << ", " << psi << ") --> "<< e << "\n";
37  cout << " e.phi() = " << e.phi()
38  << " rp.phi() = " << rp.phi()
39  << " difference is " << (e.phi() - rp.phi())/del << " * del\n";
40  if ( std::fabs (e.phi() - rp.phi()) > 200*del ) {
41  cout << "phi discrepancy in " << rp << "\n";
42  cout << " e.phi() = " << e.phi()
43  << " rp.phi() = " << rp.phi()
44  << " difference is " << e.phi() - rp.phi() << "\n";
45  if (del < 1.0e-4) cout << "??????????\n";
46  retval = false;
47  }
48  if ( std::fabs (e.theta() - rp.theta()) > 200*del ) {
49  cout << "theta discrepancy in " << rp << "\n";
50  cout << " e.theta() = " << e.theta()
51  << " rp.theta() = " << rp.theta()
52  << " difference is " << e.theta() - rp.theta() << "\n";
53  if (del < 1.0e-4) cout << "??????????\n";
54  retval = false;
55  }
56  cout << " e.psi() = " << e.psi()
57  << " rp.psi() = " << rp.psi()
58  << " difference is " << (e.psi() - rp.psi())/del << " * del\n";
59  if ( std::fabs (e.psi() - rp.psi()) > 200*del ) {
60  cout << "psi discrepancy in " << rp << "\n";
61  cout << " e.psi() = " << e.psi()
62  << " rp.psi() = " << rp.psi()
63  << " difference is " << e.psi() - rp.psi() << "\n";
64  if (del < 1.0e-4) cout << "??????????\n";
65  retval = false;
66  }
67  rpe.set(e);
68  retval &= compareR (rpe, rp, del);
69  }
70  return retval;
71 }
72 
73 int main () {
74 
75  bool res = true;
76  double PI = CLHEP::pi;
77 
78  // Some cases not in the potentially unstable region:
79  res &= test ( .05, PI/5, .1 );
80  res &= test ( .4, PI/7, -.35 );
81  res &= test ( PI/3, PI/6, -PI/3 );
82  res &= test ( PI/5, PI/2, -PI/2.5 );
83  res &= test ( -PI/5, PI/2, PI/3 );
84  res &= test ( -4*PI/5, PI/2, PI/2.5 );
85  res &= test ( 4*PI/5, PI/2, 2*PI/3 );
86  res &= test ( -3*PI/4, PI/2, -PI/3 );
87  res &= test ( 5*PI/6, PI/2, -4*PI/5 );
88  res &= test ( 5*PI/6, PI/2, -PI/3 );
89 
90  // Specialized cases
91  res &= test ( .05, 0, .1 );
92  res &= test ( .2, 0, 1.1 );
93  res &= test ( -.4, 0, .4 );
94  res &= test ( -2.4, 0, 2.0 );
95  res &= test ( -2.4, 0, -2.0 );
96  res &= test ( -2.2, 0, 1.8 );
97  res &= test ( -2.2, 0, -1.8 );
98  res &= test ( .05, PI, .1 );
99  res &= test ( .2, PI, 1.1 );
100  res &= test ( -.4, PI, .4 );
101  res &= test ( -2.4, PI, 2.0 );
102  res &= test ( -2.4, PI, -2.0 );
103  res &= test ( -2.2, PI, 1.8 );
104  res &= test ( -2.2, PI, -1.8 );
105 
106  // Cases near zero
107  res &= test ( .1, .0000000004, .5 );
108  res &= test ( -1.2, .0000000004, .5 );
109  res &= test ( .7, .0000000004, -.6 );
110  res &= test ( 1.5, .0000000004, -1.1 );
111  res &= test ( 1.4, .0000000004, -1.5 );
112  res &= test ( -.1, .0000000000028, .5 );
113  res &= test ( -1.2, .0000000000028, -.5 );
114  res &= test ( .7, .0000000000028, -.6 );
115  res &= test ( -1.5, .0000000000028, -1.1 );
116  res &= test ( 1.4, .0000000000028, 1.5 );
117 
118  // Cases near PI
119  double nearPI = PI - .00000002;
120  res &= test ( .1, nearPI, .5 );
121  res &= test ( -1.2, nearPI, .5 );
122  res &= test ( .7, nearPI, -.6 );
123  res &= test ( 1.5, nearPI, -1.1 );
124  res &= test ( 1.4, nearPI, -1.5 );
125  res &= test ( 2.4, nearPI, -1.6 );
126  res &= test ( 2.3, nearPI, 1.9 );
127  res &= test ( -2.8, nearPI, .6 );
128  res &= test ( -.4, nearPI, -3.1 );
129  nearPI = PI - .000000000009;
130  res &= test ( .1, nearPI, -.5 );
131  res &= test ( 1.2, nearPI, .5 );
132  res &= test ( .7, nearPI, -.6 );
133  res &= test ( 1.5, nearPI, 1.1 );
134  res &= test ( -1.4, nearPI, -1.5 );
135  res &= test ( 2.1, nearPI, -1.2 );
136  res &= test ( 2.9, nearPI, .9 );
137  res &= test ( -2.8, nearPI, 1.6 );
138  res &= test ( -.4, nearPI, -3.0 );
139 
140  if (!res) return -1;
141  return 0;
142 }
143 
144 bool compareR ( const HepRotation & r1, const HepRotation & r2, double tol ) {
145  HepRep3x3 m1 = r1.rep3x3();
146  HepRep3x3 m2 = r2.rep3x3();
147  double flaw = 0;
148  flaw = max (flaw, (m1.xx_ - m2.xx_));
149  flaw = max (flaw, (m1.xy_ - m2.xy_));
150  flaw = max (flaw, (m1.xz_ - m2.xz_));
151  flaw = max (flaw, (m1.yx_ - m2.yx_));
152  flaw = max (flaw, (m1.yy_ - m2.yy_));
153  flaw = max (flaw, (m1.yz_ - m2.yz_));
154  flaw = max (flaw, (m1.zx_ - m2.zx_));
155  flaw = max (flaw, (m1.zy_ - m2.zy_));
156  flaw = max (flaw, (m1.zz_ - m2.zz_));
157  if (flaw > 20*std::sqrt(tol)) {
158  cout << "???????? comparison flaw at level of " << flaw << "\n"
159  << r1 << r2;
160  }
161  cout << "flaw size is " << flaw << " (" << flaw/std::sqrt(tol) << ")\n\n";
162  return (flaw <= tol);
163 }
164 
165 void perturb ( int i, const HepRotation & r, HepRotation & rp, double & del ) {
166 
167  HepRep3x3 p0 ( 1, 3, -2,
168  -1, -2, 4,
169  2, -1, -1 );
170 
171  HepRep3x3 p1 ( 1, -1, -2,
172  1, 3, -1,
173  2, -1, -3 );
174 
175  HepRep3x3 p2 ( 5, 1, -5,
176  -3, -2, 3,
177  -1, 4, -1 );
178 
179  HepRep3x3 p3 ( -2, -2, 1,
180  -1, -2, -4,
181  4, 2, -2 );
182 
183  HepRep3x3 p[4];
184  p[0] = p0;
185  p[1] = p1;
186  p[2] = p2;
187  p[3] = p3;
188 
189  int cycle = i/4;
190  double q;
191  switch (cycle){
192  case 0: q = 1.0e-14; break;
193  case 1: q = 1.0e-12; break;
194  case 2: q = 1.0e-10; break;
195  case 3: q = 1.0e-8; break;
196  case 4: q = 1.0e-6; break;
197  case 5: q = 1.0e-4; break;
198  }
199  HepRep3x3 d = p[i%4];
200  HepRep3x3 m = r.rep3x3();
201  if ((m.zz_ + q*d.zz_) < -1) {
202  d.zz_ = -d.zz_;
203  }
204  cout << "i = " << i << " q is " << q << "\n";
205  rp.set (HepRep3x3 ( m.xx_ + q*d.xx_ , m.xy_ + q*d.xy_ , m.xz_ + q*d.xz_ ,
206  m.yx_ + q*d.yx_ , m.yy_ + q*d.yy_ , m.yz_ + q*d.yz_ ,
207  m.zx_ + q*d.zx_ , m.zy_ + q*d.zy_ , m.zz_ + q*d.zz_ ) );
208  del = q;
209 }
incomplete * p0
double psi() const
Definition: RotationE.cc:113
const int Nperturbs
Definition: eulerTest.cc:20
bool test(double phi, double theta, double psi)
Definition: eulerTest.cc:25
double phi() const
HepRotation & set(const Hep3Vector &axis, double delta)
Definition: RotationA.cc:28
void perturb(int i, const HepRotation &r, HepRotation &rp, double &del)
Definition: eulerTest.cc:165
double theta() const
Definition: RotationE.cc:107
HepRep3x3 rep3x3() const
double phi() const
Definition: RotationE.cc:73
HepEulerAngles eulerAngles() const
Definition: RotationE.cc:206
bool compareR(const HepRotation &r1, const HepRotation &r2, double tol)
Definition: eulerTest.cc:144
int main()
Definition: eulerTest.cc:73