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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.8 - (show annotations) (download)
Tue May 15 14:31:12 2012 UTC (12 years, 6 months ago) by mocchiut
Branch: MAIN
Changes since 1.7: +6 -0 lines
Reprocessing bugs fixed, ToF bugs fixed, new software versions, new quaternions, IGRF bug fixed, code cleanup

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 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
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 /*Sine::Sine()
79 : TObject()
80 {
81 }
82
83 Sine::~Sine()
84 {
85 }*/
86
87 short int Sign_1(double_t a, Int_t b){
88 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 // .* | (X0,Y0,Z0) **--.___|
113 // _** | / *. /
114 // .* | * *
115 // * ..****|***.. / R *
116 // .* | .*.
117 // .* | / *.
118 // * EARTH | / * YISK
119 // * | /_ _ _*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _|
120 // * / * /
121 // * / .*
122 // *. / .*
123 // **/*******
124 // /
125 // /
126 // /
127 // /
128 // /
129 // /
130 // |/
131 // *--
132 // XISK
133 //
134 //****************************************************************************************************/
135 //****************************************************************************************************/
136
137
138 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
140 double_t a = 360/(2*TMath::Pi());
141
142 TMatrixD Xij(3,3);
143 Xij(0,0) = 1; Xij(0,1) = 0; Xij(0,2) = 0;
144 Xij(1,0) = 0; Xij(1,1) = 0; Xij(1,2) = 1;
145 Xij(2,0) = 0; Xij(2,1) = -1; Xij(2,2) = 0;
146
147 TMatrixD Zij(3,3);
148 Zij(0,0) = 0; Zij(0,1) = 0; Zij(0,2) = -1;
149 Zij(1,0) = -1; Zij(1,1) = 0; Zij(1,2) = 0;
150 Zij(2,0) = 0; Zij(2,1) = 1; Zij(2,2) = 0;
151
152 TMatrixD Pij(3,3);
153 Pij(0,0) = pow(q0,2)+pow(q1,2)-pow(q2,2)-pow(q3,2);
154 Pij(0,1) = /*2*(q1*q2+q0*q3);/*/ 2*(q1*q2-q0*q3);
155 Pij(0,2) = /*2*(q1*q3-q0*q2);/*/ 2*(q1*q3+q0*q2);
156 Pij(1,0) = /*2*(q1*q2-q0*q3);/*/ 2*(q1*q2+q0*q3);
157 Pij(1,1) = pow(q0,2)-pow(q1,2)+pow(q2,2)-pow(q3,2);
158 Pij(1,2) = /*2*(q2*q3+q0*q1);/*/ 2*(q2*q3-q0*q1);
159 Pij(2,0) = /*2*(q1*q3+q0*q2);/*/ 2*(q1*q3-q0*q2);
160 Pij(2,1) = /*2*(q2*q3-q0*q1);/*/ 2*(q2*q3+q0*q1);
161 Pij(2,2) = pow(q0,2)-pow(q1,2)-pow(q2,2)+pow(q3,2);
162
163 TMatrixD Aij(3,3);
164
165 Double_t C1 = y0*Vz0 - z0*Vy0;
166 Double_t C2 = z0*Vx0 - x0*Vz0;
167 Double_t C3 = x0*Vy0 - y0*Vx0;
168 Double_t C = sqrt(pow(C1,2) + pow(C2,2) + pow(C3,2));
169 Double_t V0 = sqrt(pow(Vx0,2)+pow(Vy0,2) + pow(Vz0,2));
170 // Double_t R0 = sqrt(pow(x0,2)+pow(y0,2) + pow(z0,2));
171 Aij(0,0) = /*(C2*z0-C3*y0)/(C*R0);/*/Vx0/V0;
172 Aij(0,1) = C1/C;
173 Aij(0,2) = /*x0/R0;/*/(Vy0*C3-Vz0*C2)/(V0*C);
174 Aij(1,0) = /*(C3*x0-C1*z0)/(C*R0);/*/Vy0/V0;
175 Aij(1,1) = C2/C;
176 Aij(1,2) = /*y0/R0;/*/(Vz0*C1-Vx0*C3)/(V0*C);
177 Aij(2,0) = /*(C1*y0-C2*x0)/(C*R0);/*/Vz0/V0;
178 Aij(2,1) = C3/C;
179 Aij(2,2) = /*x0/R0;/*/(Vx0*C2-Vy0*C1)/(V0*C);
180 Aij.Invert();
181
182 TMatrixD Full_(3,3);
183
184 Full_ = Aij*(Pij*Zij);
185
186 //Double_t u13 = Full_(0,2);
187 //Double_t u23 = Full_(1,2);
188 //Double_t u22 = Full_(1,1);
189 //Double_t u33 = Full_(2,2);
190 //Double_t u21 = Full_(1,0);
191
192 Double_t u13 = Full_(0,0);
193 Double_t u23 = -Full_(1,0);
194 Double_t u22 = Full_(1,1);
195 Double_t u33 = Full_(2,0);
196 Double_t u21 = Full_(1,2);
197
198 Tangazh = a*atan(-u13/u33);
199 Kren = a*atan(-u23/sqrt(1 - pow(u23,2)));
200 Ryskanie = a*atan(u21/u22);
201
202 Full_.Delete();
203 Aij.Delete();
204 Pij.Delete();
205 Zij.Delete();
206 Xij.Delete();
207
208 return ;
209 }
210
211
212 void InclinationInfo::Clear(Option_t *t){
213 //Int_t gyh = 0;
214 }
215
216
217 //ClassImp(McmdItem)
218 ClassImp(InclinationInfoI)
219 ClassImp(Quaternions)
220 ClassImp(InclinationInfo)
221 //ClassImp(Sine)

  ViewVC Help
Powered by ViewVC 1.1.23