source: Sophya/trunk/SophyaLib/SkyMap/vector3d.cc@ 1827

Last change on this file since 1827 was 1177, checked in by ansari, 25 years ago

bug constructeur Vector3d(theta, phi) corrige.

dominique
File size: 7.6 KB
Line 
1#include "machdefs.h"
2#include <math.h>
3#include "pexceptions.h"
4#include "vector3d.h"
5#include "utilgeom.h"
6//++
7// Class Vector3d
8//
9// include vector3d.h utilgeom.h longlat.h math.h
10//
11//
12// 3-D geometry.
13// All computations are made with angles in radians and with spherical
14// coordinates theta, phi.
15//
16// Concerning Euler's angles, the reference is :
17//
18// "Classical Mechanics" 2nd edition, H. Goldstein, Addison Wesley
19//--
20//++
21// Titre Constructors
22//--
23//++
24Vector3d::Vector3d()
25//
26//--
27{
28 Setxyz(1.,0.,0.);
29}
30//++
31Vector3d::Vector3d(double x, double y, double z)
32//
33//--
34{
35 _x=x;
36 _y=y;
37 _z=z;
38 xyz2ThetaPhi();
39}
40//++
41Vector3d::Vector3d(double theta, double phi)
42//
43//--
44{
45 // _theta=mod(theta,M_PI); // dans [0;pi]
46 // Version precedente fausse: _theta=M_PI est valide. Or mod(M_PI,M_PI)=0!
47 // De plus theta>pi ou <0 n'a pas de sens. Dominique Yvon
48 if( (theta<0.) || (theta>M_PI) )
49 { string exmsg = "Wrong initialisation of theta in Vector3d::Vector3d(double theta, double phi)";
50 throw( ParmError(exmsg) );
51 }
52 _theta=theta; // dans [0;pi]
53 _phi=mod(phi,pi2); // dans [0;2pi]
54 ThetaPhi2xyz();
55}
56//++
57Vector3d::Vector3d(const LongLat& ll)
58//
59//--
60{
61 _theta=ll.Theta(); // dans [0;pi]
62 _phi=ll.Phi(); // dans [0;2pi]
63 ThetaPhi2xyz();
64}
65//++
66Vector3d::Vector3d(const Vector3d& v)
67//
68//--
69{
70 _x=v._x;
71 _y=v._y;
72 _z=v._z;
73 _theta=v._theta;
74 _phi=v._phi;
75}
76//++
77// Titre Public methods
78//--
79//++
80void Vector3d::SetThetaPhi(double theta, double phi)
81//
82//--
83{
84 _theta=mod(theta,M_PI);
85 _phi=mod(phi,pi2);
86 ThetaPhi2xyz();
87}
88//++
89void Vector3d::Setxyz(double x, double y, double z)
90//
91//--
92{
93 _x=x;
94 _y=y;
95 _z=z;
96 xyz2ThetaPhi();
97}
98//++
99void Vector3d::ThetaPhi2xyz()
100//
101//--
102{
103 _x=sin(_theta)*cos(_phi);
104 _y=sin(_theta)*sin(_phi);
105 _z=cos(_theta);
106}
107//++
108void Vector3d::xyz2ThetaPhi()
109//
110//--
111{
112 double norm=this->Norm();
113 if( norm != 0. )
114 {
115 _theta=acos(_z/norm); // dans [0,Pi]
116 if( mod(_theta,M_PI) == 0. ) _phi=0.; // on est sur +-Oz, le vecteur z est en phi=0
117 // else _phi=acos(_x/sin(_theta)/norm)+M_PI*(_y<0);
118 else _phi=scangle(_y/sin(_theta)/norm,_x/sin(_theta)/norm);
119 }
120 else // vecteur nul
121 {
122 _theta=0.;
123 _phi=0.;
124 }
125}
126//++
127Vector3d& Vector3d::Normalize()
128//
129//--
130{
131 double norm=this->Norm();
132 if( norm != 0. ) (*this)/=norm;
133 else cerr << "Division par zero" << endl;
134 return *this;
135}
136//++
137double Vector3d::Norm() const
138//
139//--
140{
141 return sqrt(_x*_x+_y*_y+_z*_z);
142}
143//++
144double Vector3d::Psc(const Vector3d& v) const
145//
146// dot product
147//--
148{
149 return _x*v._x+_y*v._y+_z*v._z;
150}
151//++
152double Vector3d::SepAngle(const Vector3d& v) const
153//
154// angular gap between 2 vectors in [0,Pi]
155//--
156{
157 double n1=this->Norm();
158 double n2=v.Norm();
159 double ret;
160 if( n1!=0. && n2!=0. ) ret=acos((this->Psc(v))/n1/n2);
161 else
162 {
163 cerr << "Division par zero" << endl;
164 ret=0.;
165 }
166 return ret;
167}
168//++
169Vector3d Vector3d::Vect(const Vector3d& v) const
170//
171// vector product
172//--
173{
174 double xo=_y*v._z-_z*v._y;
175 double yo=_z*v._x-_x*v._z;
176 double zo=_x*v._y-_y*v._x;
177 return Vector3d(xo,yo,zo);
178}
179//++
180Vector3d Vector3d::VperpPhi() const
181//
182// perpendicular vector, with equal phi
183//--
184{
185 double theta;
186 if( _theta != pi_over_2 ) theta=_theta+(0.5-(_theta>pi_over_2))*M_PI; // on tourne theta de +-pi/2
187 else theta=0.;
188 return Vector3d(theta,_phi);
189}
190//++
191Vector3d Vector3d::VperpTheta() const
192//
193// perpendicular vector with equal theta
194//--
195{ cerr<< " Erreur in Vector3d::VperpTheta()"<<endl;
196 throw PError("Vector3d::VperpTheta() - Logic Error DY/Reza 20/02/2000");
197 // Bug ??? (D. Yvon, Fevrier 2000)
198 // double phi=mod(_phi+pi_over_2,pi2); // on tourne phi de pi/2
199 // return Vector3d(_theta,phi);
200}
201
202Vector3d Vector3d::EPhi() const
203{
204 Vector3d temp;
205 if ( fabs(_z) == 1. ) // si on est en +- Oz
206 {
207 temp=Vector3d(1.,0.,0.);
208 }
209 else
210 {
211 Vector3d k(0,0,-1);
212 temp=this->Vect(k);
213 temp.Normalize();
214 }
215 return temp;
216}
217//++
218Vector3d Vector3d::ETheta() const
219//
220//--
221{
222 Vector3d temp=this->Vect(EPhi());
223 temp.Normalize();
224 return temp;
225}
226
227//++
228Vector3d Vector3d::Euler(double phi, double theta, double psi) const
229//
230// Euler's rotations
231//--
232{
233 double cpsi=cos(psi);
234 double ctheta=cos(theta);
235 double cphi=cos(phi);
236 double spsi=sin(psi);
237 double stheta=sin(theta);
238 double sphi=sin(phi);
239 double xnew=(cpsi*cphi-ctheta*sphi*spsi)*_x
240 +(cpsi*sphi+ctheta*cphi*spsi)*_y
241 +spsi*stheta*_z;
242 double ynew=(-spsi*cphi-ctheta*sphi*cpsi)*_x
243 +(-spsi*sphi+ctheta*cphi*cpsi)*_y
244 +cpsi*stheta*_z;
245 double znew=stheta*sphi*_x-stheta*cphi*_y+ctheta*_z;
246 return Vector3d(xnew,ynew,znew);
247}
248//++
249Vector3d Vector3d::InvEuler(double phi, double theta, double psi) const
250//
251// inverse rotation
252//--
253{
254 double cpsi=cos(psi);
255 double ctheta=cos(theta);
256 double cphi=cos(phi);
257 double spsi=sin(psi);
258 double stheta=sin(theta);
259 double sphi=sin(phi);
260 double xnew=(cpsi*cphi-ctheta*sphi*spsi)*_x
261 -(spsi*cphi+ctheta*sphi*cpsi)*_y
262 +sphi*stheta*_z;
263 double ynew=(cpsi*sphi+ctheta*cphi*spsi)*_x
264 +(-spsi*sphi+ctheta*cphi*cpsi)*_y
265 -cphi*stheta*_z;
266 double znew=stheta*spsi*_x+stheta*cpsi*_y+ctheta*_z;
267 return Vector3d(xnew,ynew,znew);
268}
269//++
270Vector3d Vector3d::Rotate(const Vector3d& omega, double phi) const
271//
272// rotation of angle phi around an axis omega (Maxwell's rule)
273//--
274{
275 Vector3d rotationaxis=omega;
276 rotationaxis.Normalize();
277 double n=this->Psc(rotationaxis);
278 Vector3d myself=*this;
279 Vector3d rotate=n*rotationaxis+(myself-n*rotationaxis)*cos(phi)-(myself^rotationaxis)*sin(phi);
280 return rotate;
281}
282//++
283void Vector3d::Print(ostream& os) const
284//
285//--
286{
287 os << "Vector3: (X,Y,Z)= (" << _x << ", " << _y << ", " << _z
288 << ") --- Theta,Phi= " << _theta << ", " << _phi << "\n"
289 << "Norme = " << this->Norm() << endl;
290}
291//++
292// Titre Operators
293//--
294//++
295Vector3d& Vector3d::operator += (const Vector3d& v)
296//
297//--
298{
299 *this=*this+v;
300 return *this;
301}
302//++
303Vector3d& Vector3d::operator -= (const Vector3d& v)
304//
305//--
306{
307 *this=*this-v;
308 return *this;
309}
310//++
311Vector3d& Vector3d::operator += (double d)
312//
313//--
314{
315 Setxyz(_x+d,_y+d,_z+d);
316 return *this;
317}
318//++
319Vector3d& Vector3d::operator /= (double d)
320//
321//--
322{
323 if( d != 0. ) Setxyz(_x/d,_y/d,_z/d);
324 else cerr << "Division par zero." << endl;
325 return *this;
326}
327//++
328Vector3d& Vector3d::operator *= (double d)
329//
330//--
331{
332 Setxyz(_x*d,_y*d,_z*d);
333 return *this;
334}
335//++
336Vector3d Vector3d::operator ^ (const Vector3d& v) const
337//
338// vector product
339//--
340{
341 return this->Vect(v);
342}
343//++
344Vector3d Vector3d::operator + (double d) const
345//
346//--
347{
348 return Vector3d(_x+d,_y+d,_z+d);
349}
350//++
351Vector3d Vector3d::operator + (const Vector3d& v) const
352//
353//--
354{
355 return Vector3d(_x+v._x,_y+v._y,_z+v._z);
356}
357//++
358Vector3d Vector3d::operator - (double d) const
359//
360//--
361{
362 return *this+(-d);
363}
364//++
365Vector3d Vector3d::operator - (const Vector3d& v) const
366//
367//--
368{
369 return *this+(v*(-1.));
370}
371//++
372Vector3d Vector3d::operator * (double d) const
373//
374//--
375{
376 return Vector3d(d*_x,d*_y,d*_z);
377}
378//++
379double Vector3d::operator * (const Vector3d& v) const
380//
381// dot product
382//--
383{
384 return this->Psc(v);
385}
386//++
387Vector3d Vector3d::operator / (double d) const
388//
389//--
390{
391 Vector3d ret=*this;
392 if( d != 0. ) ret/=d;
393 else cerr << "Division par zero." << endl;
394 return ret;
395}
396//++
397Vector3d& Vector3d::operator = (const Vector3d& v)
398//
399//--
400{
401 if( this != &v )
402 {
403 _x=v._x;
404 _y=v._y;
405 _z=v._z;
406 _theta=v._theta;
407 _phi=v._phi;
408 }
409 return *this;
410}
411//++
412bool Vector3d::operator == (const Vector3d& v)
413//
414//--
415{
416 return (this==&v);
417}
418
Note: See TracBrowser for help on using the repository browser.