/[PAMELA software]/DarthVader/OrbitalInfo/src/InclinationInfo.cpp
ViewVC logotype

Annotation of /DarthVader/OrbitalInfo/src/InclinationInfo.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.9 - (hide annotations) (download)
Fri Mar 28 20:47:15 2014 UTC (10 years, 8 months ago) by pam-mep
Branch: MAIN
CVS Tags: v10REDr01, v10RED, HEAD
Changes since 1.8: +132 -25 lines
new OrbitalInfo code

1 pamelaprod 1.1 /***************************************************************************
2     * Copyright (C) 2006 by pamelaprod *
3     * pamelaprod@P1.pamela *
4     * *
5     * This program is free software; you can redistribute it and/or modify *
6     * it under the terms of the GNU General Public License as published by *
7     * the Free Software Foundation; either version 2 of the License, or *
8     * (at your option) any later version. *
9     * *
10     * This program is distributed in the hope that it will be useful, *
11     * but WITHOUT ANY WARRANTY; without even the implied warranty of *
12     * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
13     * GNU General Public License for more details. *
14     * *
15     * You should have received a copy of the GNU General Public License *
16     * along with this program; if not, write to the *
17     * Free Software Foundation, Inc., *
18     * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
19     ***************************************************************************/
20     #include <InclinationInfo.h>
21    
22     using namespace std;
23    
24     // InclinationInfoI()::InclinationInfoI() {
25     // // memset(time,0,6*sizeof(double));
26     // // memset(quad,0,6*4*sizeof(double));
27     // };
28    
29     void InclinationInfoI::fill(TArrayC* data){
30     short extIndex = 0;
31     short innIndex = 0;
32     long tempData = 0;
33     for (int i = 0; i < 6; i++){
34     extIndex = 20*i;
35     time[i] = (((data->At(extIndex) << 24) & 0xFF000000) +
36     ((data->At(extIndex + 1) << 16) & 0x00FF0000) + ((data->At(extIndex + 2) << 8) & 0x0000FF00) +
37     (data->At(extIndex + 3) & 0x000000FF))/128.0;
38     for (int j = 0; j < 4; j++){
39     innIndex = extIndex + 4*j;
40     tempData = ((data->At(innIndex + 4) << 24) & 0xFF000000) + ((data->At(innIndex + 5) << 16) & 0x00FF0000) + ((data->At(innIndex + 6) << 8) & 0x0000FF00) + (data->At(innIndex + 7) & 0x000000FF);
41     if (data->At(innIndex + 4) >> 8) {
42     quat[i][j] = (~tempData * -1.0)/1073741824.0;
43     } else {
44     quat[i][j] = tempData / 1073741824.0;
45     }
46     }
47     }
48     }
49    
50 mocchiut 1.4 void InclinationInfoI::clear() {
51     for(UInt_t i = 0; i < 6; i++){
52     time[i]=0;
53     for(UInt_t j = 0; j < 4; j++) quat[i][j]=0;
54     }
55     return ;
56     }
57 pamelaprod 1.1
58    
59     Quaternions::Quaternions()
60     : InclinationInfoI()
61     {
62     }
63    
64    
65     Quaternions::~Quaternions()
66     {
67     }
68    
69     InclinationInfo::InclinationInfo()
70     : TObject()
71     {
72     }
73    
74     InclinationInfo::~InclinationInfo()
75     {
76     }
77    
78 mocchiut 1.7 /*Sine::Sine()
79     : TObject()
80     {
81     }
82    
83     Sine::~Sine()
84     {
85     }*/
86    
87 mocchiut 1.4 short int Sign_1(double_t a, Int_t b){
88 pamelaprod 1.1 if(a>0){b=1;}
89     if(a<0){b=-1;}
90     else{b=0;}
91     return b;
92     }
93    
94    
95     /******************************************************************************************************************/
96     /******************************************************************************************************************/
97     //********************* ***************************************************************/
98     //********************* COORDINATE SYSTEMS ***************************************************************/
99     //********************* ***************************************************************/
100     //*****************************************************************************************************************/
101     //*****************************************************************************************************************/
102     //
103     // ZISK
104     // +
105     // / \ YOSK ZOSK (Directed by Radius)
106     // | _ _.
107     // | |\ /|
108     // | \ /
109     // | \ /
110     // |.__..__ \ /
111     // Orbit _._.***| **.\/_ XOSK (Directed by velocity)
112 mocchiut 1.4 // .* | (X0,Y0,Z0) **--.___|
113 pamelaprod 1.1 // _** | / *. /
114     // .* | * *
115     // * ..****|***.. / R *
116     // .* | .*.
117     // .* | / *.
118     // * EARTH | / * YISK
119 mocchiut 1.4 // * | /_ _ _*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _|
120 pamelaprod 1.1 // * / * /
121     // * / .*
122     // *. / .*
123     // **/*******
124     // /
125     // /
126     // /
127     // /
128     // /
129     // /
130     // |/
131     // *--
132     // XISK
133     //
134     //****************************************************************************************************/
135     //****************************************************************************************************/
136    
137    
138 mocchiut 1.4 void InclinationInfo::TransAngle(Double_t x0, Double_t y0, Double_t z0, Double_t Vx0, Double_t Vy0, Double_t Vz0, Double_t q0, Double_t q1, Double_t q2, Double_t q3){
139 pamelaprod 1.1
140 pam-mep 1.9 cout.precision(12);
141    
142 mocchiut 1.4 double_t a = 360/(2*TMath::Pi());
143 pamelaprod 1.1
144 mocchiut 1.4 TMatrixD Xij(3,3);
145     Xij(0,0) = 1; Xij(0,1) = 0; Xij(0,2) = 0;
146     Xij(1,0) = 0; Xij(1,1) = 0; Xij(1,2) = 1;
147     Xij(2,0) = 0; Xij(2,1) = -1; Xij(2,2) = 0;
148    
149     TMatrixD Zij(3,3);
150 pam-mep 1.9 Zij(0,0) = 0.0; Zij(0,1) = 0.0; Zij(0,2) = -1.0;
151     Zij(1,0) = -1.0; Zij(1,1) = 0.0; Zij(1,2) = 0.0;
152     Zij(2,0) = 0.0; Zij(2,1) = 1.0; Zij(2,2) = 0.0;
153 pamelaprod 1.1
154     TMatrixD Pij(3,3);
155     Pij(0,0) = pow(q0,2)+pow(q1,2)-pow(q2,2)-pow(q3,2);
156 mocchiut 1.4 Pij(0,1) = /*2*(q1*q2+q0*q3);/*/ 2*(q1*q2-q0*q3);
157     Pij(0,2) = /*2*(q1*q3-q0*q2);/*/ 2*(q1*q3+q0*q2);
158     Pij(1,0) = /*2*(q1*q2-q0*q3);/*/ 2*(q1*q2+q0*q3);
159 pamelaprod 1.1 Pij(1,1) = pow(q0,2)-pow(q1,2)+pow(q2,2)-pow(q3,2);
160 mocchiut 1.4 Pij(1,2) = /*2*(q2*q3+q0*q1);/*/ 2*(q2*q3-q0*q1);
161     Pij(2,0) = /*2*(q1*q3+q0*q2);/*/ 2*(q1*q3-q0*q2);
162     Pij(2,1) = /*2*(q2*q3-q0*q1);/*/ 2*(q2*q3+q0*q1);
163 pamelaprod 1.1 Pij(2,2) = pow(q0,2)-pow(q1,2)-pow(q2,2)+pow(q3,2);
164    
165     TMatrixD Aij(3,3);
166 mocchiut 1.4
167 pamelaprod 1.1 Double_t C1 = y0*Vz0 - z0*Vy0;
168     Double_t C2 = z0*Vx0 - x0*Vz0;
169     Double_t C3 = x0*Vy0 - y0*Vx0;
170     Double_t C = sqrt(pow(C1,2) + pow(C2,2) + pow(C3,2));
171     Double_t V0 = sqrt(pow(Vx0,2)+pow(Vy0,2) + pow(Vz0,2));
172 mocchiut 1.4 // Double_t R0 = sqrt(pow(x0,2)+pow(y0,2) + pow(z0,2));
173 pam-mep 1.9 Aij(0,0) = Vx0/V0;
174 pamelaprod 1.1 Aij(0,1) = C1/C;
175 pam-mep 1.9 Aij(0,2) = (Vy0*C3-Vz0*C2)/(V0*C);
176     Aij(1,0) = Vy0/V0;
177 pamelaprod 1.1 Aij(1,1) = C2/C;
178 pam-mep 1.9 Aij(1,2) = (Vz0*C1-Vx0*C3)/(V0*C);
179     Aij(2,0) = Vz0/V0;
180 pamelaprod 1.1 Aij(2,1) = C3/C;
181 pam-mep 1.9 Aij(2,2) = (Vx0*C2-Vy0*C1)/(V0*C);
182    
183     TMatrixD Bij(3,3);
184     Bij(0,0) = Vx0/V0;
185     Bij(1,0) = C1/C;
186     Bij(2,0) = (Vy0*C3-Vz0*C2)/(V0*C);
187     Bij(0,1) = Vy0/V0;
188     Bij(1,1) = C2/C;
189     Bij(2,1) = (Vz0*C1-Vx0*C3)/(V0*C);
190     Bij(0,2) = Vz0/V0;
191     Bij(1,2) = C3/C;
192     Bij(2,2) = (Vx0*C2-Vy0*C1)/(V0*C);
193     /*
194     cout<<"Coordinates:"<<endl;
195     cout<<x0<<"\t"<<y0<<"\t"<<z0<<"\t"<<Vx0/V0<<"\t"<<Vy0/V0<<"\t"<<Vz0/V0<<endl;
196    
197     cout<<"Pij"<<endl;
198     cout<<Pij(0,0)<<"\t"<<Pij(0,1)<<"\t"<<Pij(0,2)<<endl;
199     cout<<Pij(1,0)<<"\t"<<Pij(1,1)<<"\t"<<Pij(1,2)<<endl;
200     cout<<Pij(2,0)<<"\t"<<Pij(2,1)<<"\t"<<Pij(2,2)<<endl;
201     */
202     //Aij.Invert();
203 pamelaprod 1.1
204 mocchiut 1.4 TMatrixD Full_(3,3);
205 pamelaprod 1.1
206 pam-mep 1.9 Full_ = Bij*(Pij*Zij);
207     /*/temprary
208     TMatrixD Tmp(3,3);
209     Tmp=Pij*Zij;
210    
211     cout<<"Quaternion matrix (Tmp)"<<endl;
212     cout<<Tmp(0,0)<<"\t"<<Tmp(0,1)<<"\t"<<Tmp(0,2)<<endl;
213     cout<<Tmp(1,0)<<"\t"<<Tmp(1,1)<<"\t"<<Tmp(1,2)<<endl;
214     cout<<Tmp(2,0)<<"\t"<<Tmp(2,1)<<"\t"<<Tmp(2,2)<<endl;
215    
216     cout<<"Orientation based on Velocity"<<endl;
217     cout<<Aij(0,0)<<"\t"<<Aij(0,1)<<"\t"<<Aij(0,2)<<endl;
218     cout<<Aij(1,0)<<"\t"<<Aij(1,1)<<"\t"<<Aij(1,2)<<endl;
219     cout<<Aij(2,0)<<"\t"<<Aij(2,1)<<"\t"<<Aij(2,2)<<endl;
220    
221     cout<<"Satellite Axis in Velocity Reference frame (Full_):"<<endl;
222     cout<<Full_(0,0)<<"\t"<<Full_(0,1)<<"\t"<<Full_(0,2)<<endl;
223     cout<<Full_(1,0)<<"\t"<<Full_(1,1)<<"\t"<<Full_(1,2)<<endl;
224     cout<<Full_(2,0)<<"\t"<<Full_(2,1)<<"\t"<<Full_(2,2)<<endl;
225     */
226     Double_t u00 = Full_(0,0);
227     Double_t u10 = Full_(1,0);
228     Double_t u11 = Full_(1,1);
229     Double_t u20 = Full_(2,0);
230     Double_t u12 = Full_(1,2);
231    
232     //Double_t u13 = Full_(0,0);
233     //Double_t u23 = -Full_(1,0);
234 mocchiut 1.4 //Double_t u22 = Full_(1,1);
235 pam-mep 1.9 //Double_t u33 = Full_(2,0);
236     //Double_t u21 = Full_(1,2);
237 mocchiut 1.4
238 pam-mep 1.9 Tangazh = a*atan(-u00/u20);
239     Kren = a*atan(u10/sqrt(1 - pow(u10,2)));
240     Ryskanie = a*atan(u12/u11);
241    
242     //Tangazh = a*atan(-u13/u33);
243     //Kren = a*atan(-u23/sqrt(1 - pow(u23,2)));
244     //Ryskanie = a*atan(u21/u22);
245     // end temprary
246    
247     //10RED CHECK
248     /*
249     u10 = tan(Kren*TMath::DegToRad())/sqrt(pow(tan(Kren*TMath::DegToRad()),2)+1);
250     u11 = -sqrt((1-pow(u10,2))/(1+pow(tan(Ryskanie*TMath::DegToRad()),2)));
251     u12 = u11*tan(Ryskanie*TMath::DegToRad());
252     u20 = -sqrt((1-pow(u10,2))/(1+pow(tan(Tangazh*TMath::DegToRad()),2)));
253     u00 = -u20*tan(Tangazh*TMath::DegToRad());
254    
255     Double_t aa = 1+pow((u20/u00),2);
256     Double_t by = 2*u10*u11*u20/pow(u00,2);
257     Double_t cy = (1+pow(u10/u00,2))*pow(u11,2)-1;
258     Double_t bz = 2*u10*u12*u20/pow(u00,2);
259     Double_t cz = (1+pow(u10/u00,2))*pow(u12,2)-1;
260    
261     Int_t uj = TMath::Sign(1.,Ryskanie)*TMath::Sign(1.,Tangazh);
262     Double_t u21 = (-by+uj*sqrt(pow(by,2)-4*aa*cy))/(2*aa);
263     Double_t u21s = -TMath::Sign(1.,Kren)*TMath::Abs(u21);
264     Double_t u01 = TMath::Sign(1.,Ryskanie)*TMath::Abs((u10*u11+u20*u21)/u00);
265    
266     Int_t fj=1;
267     if(TMath::Sign(1.,Tangazh)>0 && TMath::Sign(1.,Ryskanie)>0) fj=-1;
268    
269     Double_t u22 = (-bz+fj*sqrt(pow(bz,2)-4*aa*cz))/(2*aa);
270     Double_t u22s = -TMath::Sign(1.,Tangazh)*TMath::Abs(u22);
271     Double_t u02 = -TMath::Abs((u10*u12+u20*u22)/u00);
272    
273     TMatrixD Dij(3,3);
274     Dij(0,0) = u00; Dij(0,1) = u01; Dij(0,2) = u02;
275     Dij(1,0) = u10; Dij(1,1) = u11; Dij(1,2) = u12;
276     Dij(2,0) = u20; Dij(2,1) = u21s; Dij(2,2) = u22s;
277    
278     //cout<<"Dij"<<endl;
279     //cout<<Dij(0,0)<<"\t"<<Dij(0,1)<<"\t"<<Dij(0,2)<<endl;
280     //cout<<Dij(1,0)<<"\t"<<Dij(1,1)<<"\t"<<Dij(1,2)<<endl;
281     //cout<<Dij(2,0)<<"\t"<<Dij(2,1)<<"\t"<<Dij(2,2)<<endl;
282    
283     //Aij.Invert();
284     //Zij.Invert();
285     TMatrixD Shij(3,3);
286     TMatrixD Usij(3,3);
287     Usij = (Aij*Dij);
288     Usij.Invert();
289     Shij = Zij*Usij;
290     Shij.Invert();
291    
292     cout<<"Full_ matrix having got from Euler angles"<<endl;
293     cout<<Shij(0,0)<<"\t"<<Shij(0,1)<<"\t"<<Shij(0,2)<<endl;
294     cout<<Shij(1,0)<<"\t"<<Shij(1,1)<<"\t"<<Shij(1,2)<<endl;
295     cout<<Shij(2,0)<<"\t"<<Shij(2,1)<<"\t"<<Shij(2,2)<<endl;
296    
297     cout<<"Bank = "<<Kren<<"\tSPitch = "<<Tangazh<<"\tYaw = "<<Ryskanie<<endl;
298     if(TMath::Abs(Kren)>10.0){
299     // if(Vz0/V0>0.99){
300     Int_t Fer;
301     cin>>Fer;
302     }
303 pamelaprod 1.1
304 mocchiut 1.8 Full_.Delete();
305     Aij.Delete();
306     Pij.Delete();
307     Zij.Delete();
308     Xij.Delete();
309 pam-mep 1.9 Dij.Delete();
310     Shij.Delete();
311     Usij.Delete();
312 mocchiut 1.8
313 pam-mep 1.9 // END 10RED CHECK
314     */
315 pamelaprod 1.1 return ;
316     }
317    
318    
319 mocchiut 1.5 void InclinationInfo::Clear(Option_t *t){
320 mocchiut 1.4 //Int_t gyh = 0;
321     }
322    
323    
324 pamelaprod 1.1 //ClassImp(McmdItem)
325     ClassImp(InclinationInfoI)
326     ClassImp(Quaternions)
327     ClassImp(InclinationInfo)
328 mocchiut 1.7 //ClassImp(Sine)

  ViewVC Help
Powered by ViewVC 1.1.23