1 | #include <math.h>
|
---|
2 | #include "circle.h"
|
---|
3 | //++
|
---|
4 | // Class Circle
|
---|
5 | //
|
---|
6 | // include circle.h math.h
|
---|
7 | //--
|
---|
8 | //++
|
---|
9 | //
|
---|
10 | // Links Parents
|
---|
11 | //
|
---|
12 | // Geometry
|
---|
13 | //
|
---|
14 | //--
|
---|
15 | //++
|
---|
16 | // Titre Constructors
|
---|
17 | //--
|
---|
18 | //++
|
---|
19 | Circle::Circle()
|
---|
20 | //
|
---|
21 | //--
|
---|
22 | {
|
---|
23 | UnitVector temp;
|
---|
24 | SetCircle(temp,M_PI/2.);
|
---|
25 | }
|
---|
26 | //++
|
---|
27 | Circle::Circle(double theta, double phi, double aperture)
|
---|
28 | //
|
---|
29 | //--
|
---|
30 | {
|
---|
31 | UnitVector temp(theta,phi);
|
---|
32 | SetCircle(temp,aperture);
|
---|
33 | }
|
---|
34 | //++
|
---|
35 | Circle::Circle(double x, double y, double z, double aperture)
|
---|
36 | //
|
---|
37 | //--
|
---|
38 | {
|
---|
39 | UnitVector temp(x,y,z);
|
---|
40 | SetCircle(temp,aperture);
|
---|
41 | }
|
---|
42 | //++
|
---|
43 | Circle::Circle(const Vector3d& v, double aperture)
|
---|
44 | //
|
---|
45 | //--
|
---|
46 | {
|
---|
47 | UnitVector temp=v;
|
---|
48 | SetCircle(temp,aperture);
|
---|
49 | }
|
---|
50 | //++
|
---|
51 | Circle::Circle(const Circle& c)
|
---|
52 | //
|
---|
53 | // copy constructor
|
---|
54 | //--
|
---|
55 | {
|
---|
56 | UnitVector temp= c.Omega();
|
---|
57 | SetCircle(temp,c._angouv);
|
---|
58 | }
|
---|
59 | //++
|
---|
60 | // Titre Public Methods
|
---|
61 | //--
|
---|
62 | //++
|
---|
63 | void Circle::SetCircle(const UnitVector& temp, double aperture)
|
---|
64 | //
|
---|
65 | //--
|
---|
66 | {
|
---|
67 | _spinunitaxis= temp;
|
---|
68 |
|
---|
69 | _angouv = aperture;
|
---|
70 | _cangouv= cos(_angouv);
|
---|
71 | _sangouv= sin(_angouv);
|
---|
72 |
|
---|
73 | _spinaxis=_spinunitaxis*fabs(_cangouv);
|
---|
74 |
|
---|
75 | _theta =_spinunitaxis.Theta();
|
---|
76 | _ctheta= cos(_theta);
|
---|
77 | _stheta= sin(_theta);
|
---|
78 |
|
---|
79 | _phi =_spinunitaxis.Phi();
|
---|
80 | _cphi= cos(_phi);
|
---|
81 | _sphi= sin(_phi);
|
---|
82 |
|
---|
83 | _x= _spinunitaxis.X();
|
---|
84 | _y= _spinunitaxis.Y();
|
---|
85 | _z= _spinunitaxis.Z();
|
---|
86 | }
|
---|
87 | //++
|
---|
88 | void Circle::SetSpinAxis(double theta, double phi)
|
---|
89 | //
|
---|
90 | //--
|
---|
91 | {
|
---|
92 | UnitVector temp(theta,phi);
|
---|
93 | SetCircle(temp,_angouv);
|
---|
94 | }
|
---|
95 | //++
|
---|
96 | void Circle::SetSpinAxis(const Vector3d& u)
|
---|
97 | //
|
---|
98 | //--
|
---|
99 | {
|
---|
100 | UnitVector temp=u;
|
---|
101 | SetCircle(temp,_angouv);
|
---|
102 | }
|
---|
103 | //++
|
---|
104 | void Circle::SetSpinAxis(double x, double y, double z)
|
---|
105 | //
|
---|
106 | //--
|
---|
107 | {
|
---|
108 | UnitVector temp(x,y,z);
|
---|
109 | SetCircle(temp,_angouv);
|
---|
110 | }
|
---|
111 | //++
|
---|
112 | void Circle::SetApertureAngle(double aperture)
|
---|
113 | //
|
---|
114 | //--
|
---|
115 | {
|
---|
116 | SetCircle(_spinunitaxis,aperture);
|
---|
117 | }
|
---|
118 | //++
|
---|
119 | void Circle::SetApertureAngle(const Circle& c)
|
---|
120 | //
|
---|
121 | //--
|
---|
122 | {
|
---|
123 | SetCircle(_spinunitaxis,c._angouv);
|
---|
124 | }
|
---|
125 |
|
---|
126 | //++
|
---|
127 | UnitVector Circle::ConvToSphere(double psi) const
|
---|
128 | //
|
---|
129 | // Return UnitVector corresponding to a given position donnee on the circle
|
---|
130 | //--
|
---|
131 | {
|
---|
132 | psi=mod(psi,pi2);
|
---|
133 | double xout, yout, zout;
|
---|
134 | double cosa=cos(_angouv);
|
---|
135 | double sina=sin(_angouv);
|
---|
136 | double cost=cos(_theta);
|
---|
137 | double sint=sin(_theta);
|
---|
138 | double cosphi=cos(_phi);
|
---|
139 | double sinphi=sin(_phi);
|
---|
140 | double cosp=cos(psi);
|
---|
141 | double sinp=sin(psi);
|
---|
142 | xout = cosa*sint*cosphi+sina*(sinphi*sinp-cost*cosphi*cosp);
|
---|
143 | yout = cosa*sint*sinphi-sina*(cosphi*sinp+cost*sinphi*cosp);
|
---|
144 | zout = cosa*cost+sina*sint*cosp;
|
---|
145 | return UnitVector(xout,yout,zout);
|
---|
146 | }
|
---|
147 | //++
|
---|
148 | UnitVector Circle::TanOnCircle(double psi) const
|
---|
149 | //
|
---|
150 | // Return UnitVector corresponding to the tangent to the circle
|
---|
151 | // at given position on the circle.
|
---|
152 | //--
|
---|
153 | {
|
---|
154 | psi=mod(psi,pi2);
|
---|
155 | double xout, yout, zout;
|
---|
156 | double cost=cos(_theta);
|
---|
157 | double sint=sin(_theta);
|
---|
158 | double cosphi=cos(_phi);
|
---|
159 | double sinphi=sin(_phi);
|
---|
160 | double cosp=cos(psi);
|
---|
161 | double sinp=sin(psi);
|
---|
162 | xout = cosp*sinphi+sinp*sint*cosphi;
|
---|
163 | yout = -cosp*cosphi+sinp*sint*sinphi;
|
---|
164 | zout = -sinp*cost;
|
---|
165 | return UnitVector(xout,yout,zout);
|
---|
166 | }
|
---|
167 | //++
|
---|
168 | UnitVector Circle::EPhi(double psi) const
|
---|
169 | //
|
---|
170 | // Return the vector tangent to the sphere in the plane (xy)
|
---|
171 | // at a given position on the circle.
|
---|
172 | //--
|
---|
173 | {
|
---|
174 | psi=mod(psi,pi2);
|
---|
175 | return ConvToSphere(psi).EPhi();
|
---|
176 | }
|
---|
177 | //++
|
---|
178 | UnitVector Circle::ETheta(double psi) const
|
---|
179 | //
|
---|
180 | // Return the other tangent vector( orthogonal to EPhi)--
|
---|
181 | // see previous method
|
---|
182 | //--
|
---|
183 | {
|
---|
184 | psi=mod(psi,pi2);
|
---|
185 | return ConvToSphere(psi).ETheta();
|
---|
186 | }
|
---|
187 | //++
|
---|
188 | double Circle::SepAngleTanEPhi02PI(double psi) const
|
---|
189 | //
|
---|
190 | // Return separation angle in [0,2Pi] at a given position on the
|
---|
191 | // circle and EPhi
|
---|
192 | //--
|
---|
193 | {
|
---|
194 | psi=mod(psi,pi2);
|
---|
195 | UnitVector pol=this->TanOnCircle(psi);
|
---|
196 | UnitVector ephi=this->EPhi(psi);
|
---|
197 | double angle=pol.SepAngle(ephi);
|
---|
198 | if( pol.Z() <= 0 ) angle=pi2-angle;
|
---|
199 | return angle;
|
---|
200 | }
|
---|
201 | //++
|
---|
202 | void Circle::Print(ostream& os) const
|
---|
203 | //
|
---|
204 | //--
|
---|
205 | {
|
---|
206 | os << "1 - Circle - Axe de Spin Unitaire : " << _spinunitaxis << endl;
|
---|
207 | os << "1 - Circle - Axe de Spin : " << _spinaxis << endl;
|
---|
208 | os << "2 - Circle - Angle d'ouverture : " << _angouv << endl;
|
---|
209 | os << "3 - Circle - Theta,Phi : " << _theta << "," << _phi << endl;
|
---|
210 | os << "4 - Circle - x,y,z : " << _x << "," << _y << "," << _z << endl;
|
---|
211 | }
|
---|
212 | //++
|
---|
213 | //
|
---|
214 | // inline double Theta() const
|
---|
215 | // inline double Phi() const
|
---|
216 | // inline double ApertureAngle() const
|
---|
217 | // inline Vector3d Omega() const
|
---|
218 | //--
|
---|
219 | //++
|
---|
220 | // Titre Operators
|
---|
221 | //--
|
---|
222 |
|
---|
223 | Circle& Circle::operator=(const Circle& c)
|
---|
224 | {
|
---|
225 | if( this != &c )
|
---|
226 | {
|
---|
227 | UnitVector temp(c.Omega());
|
---|
228 | SetCircle(temp,c.ApertureAngle());
|
---|
229 | }
|
---|
230 | return *this;
|
---|
231 | }
|
---|
232 | //++
|
---|
233 | bool Circle::operator==(const Circle& c) const
|
---|
234 | //
|
---|
235 | //--
|
---|
236 | {
|
---|
237 | bool flag;
|
---|
238 | if( this == &c ) flag=true;
|
---|
239 | else flag=false;
|
---|
240 | return flag;
|
---|
241 | }
|
---|
242 | //++
|
---|
243 | bool Circle::operator!=(const Circle& c) const
|
---|
244 | //
|
---|
245 | //--
|
---|
246 | {
|
---|
247 | return (bool)(1-(this->operator==(c)));
|
---|
248 | }
|
---|
249 | //
|
---|
250 | bool Circle::Intersection(const Circle& c) const
|
---|
251 | {
|
---|
252 | double alphak= _angouv;
|
---|
253 | double alphal= c._angouv;
|
---|
254 | Vector3d ok= _spinaxis;
|
---|
255 | Vector3d ol= c._spinaxis;
|
---|
256 | double gamma= ok.SepAngle(ol);
|
---|
257 |
|
---|
258 | if(fabs(alphak-alphal) < gamma && gamma <= (alphak+alphal) && this != &c) {
|
---|
259 | return true;
|
---|
260 | } else {
|
---|
261 | return false;
|
---|
262 | }
|
---|
263 | }
|
---|
264 | //
|
---|
265 | bool Circle::Intersection(const Circle& c, double* psi) const
|
---|
266 | {
|
---|
267 | double alphak= _angouv;
|
---|
268 | double alphal= c._angouv;
|
---|
269 | Vector3d ok= _spinaxis;
|
---|
270 | Vector3d ol= c._spinaxis;
|
---|
271 | double gamma= ok.SepAngle(ol);
|
---|
272 |
|
---|
273 | if(fabs(alphak-alphal) < gamma && gamma <= (alphak+alphal) && this != &c) {
|
---|
274 |
|
---|
275 | double sgamma= sin(gamma);
|
---|
276 | double cgamma= cos(gamma);
|
---|
277 |
|
---|
278 | double sdphi= _sphi*c._cphi - _cphi*c._sphi;
|
---|
279 | double cdphi= _cphi*c._cphi + _sphi*c._sphi;
|
---|
280 |
|
---|
281 | double ssk= c._stheta*sdphi/sgamma;
|
---|
282 | double csk= (c._ctheta*_stheta-c._stheta*_ctheta*cdphi)/sgamma;
|
---|
283 |
|
---|
284 | double ssl= -_stheta*sdphi/sgamma;
|
---|
285 | double csl= (_ctheta*c._stheta-_stheta*c._ctheta*cdphi)/sgamma;
|
---|
286 |
|
---|
287 | double ak= atan2(ssk,csk);
|
---|
288 | double al= atan2(ssl,csl);
|
---|
289 | double omegak= acos((c._cangouv-_cangouv*cgamma)/sgamma/_sangouv);
|
---|
290 | double omegal= acos((_cangouv-c._cangouv*cgamma)/sgamma/c._sangouv);
|
---|
291 |
|
---|
292 | psi[0]= fmod(ak-omegak+pi2,pi2);
|
---|
293 | psi[1]= fmod(ak+omegak+pi2,pi2);
|
---|
294 | psi[2]= fmod(al-omegal+pi2,pi2);
|
---|
295 | psi[3]= fmod(al+omegal+pi2,pi2);
|
---|
296 |
|
---|
297 | if(psi[0] > psi[1]) {
|
---|
298 | swap(psi[0],psi[1]);
|
---|
299 | swap(psi[2],psi[3]);
|
---|
300 | }
|
---|
301 | return true;
|
---|
302 |
|
---|
303 | } else {
|
---|
304 | psi[0] = -1.;
|
---|
305 | psi[1] = -1.;
|
---|
306 | psi[2] = -1.;
|
---|
307 | psi[3] = -1.;
|
---|
308 | return false;
|
---|
309 | }
|
---|
310 | }
|
---|