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

Last change on this file since 2963 was 2615, checked in by cmv, 21 years ago

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