/[PAMELA software]/ParDirCalc/TrAng.cpp
ViewVC logotype

Annotation of /ParDirCalc/TrAng.cpp

Parent Directory Parent Directory | Revision Log Revision Log

Revision 1.1 - (hide annotations) (download)
Thu Jan 17 16:15:13 2008 UTC (16 years, 10 months ago) by cafagna
Branch: MAIN
Branch point for: first 
Initial revision
1 cafagna 1.1 #include <TObject.h>
2     #include <TrAng.h>
3     #include "TString.h"
4     #include "TMatrixD.h"
5     #include <iostream>
6     #include <stdio.h>
7    
8     using namespace std;
9    
10     PamelaOrientation::PamelaOrientation() : TObject(){
11     a = 360/(2*TMath::Pi());
12     Re = 6000000;
13     }
14    
15     PamelaOrientation::~PamelaOrientation(){
16     }
17    
18     TMatrixD PamelaOrientation::QuatoECI(Float_t q0, Float_t q1, Float_t q2, Float_t q3){
19     TMatrixD Pij(3,3);
20     Pij(0,0) = pow(q0,2)+pow(q1,2)-pow(q2,2)-pow(q3,2);
21     Pij(0,1) = /*2*(q1*q2+q0*q3);/*/ 2*(q1*q2-q0*q3);
22     Pij(0,2) = /*2*(q1*q3-q0*q2);/*/ 2*(q1*q3+q0*q2);
23     Pij(1,0) = /*2*(q1*q2-q0*q3);/*/ 2*(q1*q2+q0*q3);
24     Pij(1,1) = pow(q0,2)-pow(q1,2)+pow(q2,2)-pow(q3,2);
25     Pij(1,2) = /*2*(q2*q3+q0*q1);/*/ 2*(q2*q3-q0*q1);
26     Pij(2,0) = /*2*(q1*q3+q0*q2);/*/ 2*(q1*q3-q0*q2);
27     Pij(2,1) = /*2*(q2*q3-q0*q1);/*/ 2*(q2*q3+q0*q1);
28     Pij(2,2) = pow(q0,2)-pow(q1,2)-pow(q2,2)+pow(q3,2);
29     return Pij;
30     }
31    
32     TMatrixD PamelaOrientation::ECItoGreenwich(TMatrixD Aij, UInt_t t){
33     //t=1154304000+86400*365;
34     TMatrixD Gij(3,3);
35     Double_t omg = (7.292115e-5)*a; // Earth rotation velosity (Around polar axis);
36     //Double_t t1 = 0;
37     //if(t<=1158883200) t1 = 1127347200+229.2732; //absTime at 22/09/2005 + diference between Solar midnight and Greenwich sidereal midnight
38     //if(t>1158883200&&t<=1190419200) t1 = 1158883200+172.3415;//absTime at 22/09/2006 + diference between Solar midnight and Greenwich sidereal midnight
39     //if(t>=1190419200&&t<1222041600) t1 = 1190419200+115.39; //absTime at 22/09/2007 + diference between Solar midnight and Greenwich sidereal midnight
40     //if(t>=1222041600) t1 = 1222041600 + 294.9361; //absTime at 22/09/2008 + diference between Solar midnight and Greenwich sidereal midnight
41     //UInt_t Nd = (t-t1)/86400;
42     //Int_t DifSuSt = Nd*236.55;
43     Double_t d = (t-10957*86400-43200); //Number of day, passing from 01/01/2000 12:00:00 to t;
44     d = d/86400;
45     //d = t-da*86400+DifSuSt
46     //cout<<"t = "<<t<<"\n";
47     //cout<<"t - 2000y = "<<t-10957*86400-43200<<"\n";
48     //cout<<"d = "<<d<<"\n";
49     //Int_t tl = t%86400; //!!!!!!!!!!!!!!!!!!!!!!!!
50     Double_t T = d/36525; //Number of Julian centuries;
51    
52     //Double_t tl = t-t1-Nd*86400-DifSuSt;
53     Double_t Se = 6*3600+41*60+236.555367908*d+0.093104*pow(T,2)-(6.2e-6)*pow(T,3);
54     //cout<<"Se = "<<Se<<"\n";
55     //cout<<t<<endl<<d<<endl<<tl<<endl<<Se+omg*tl*86400/360<<endl;
56     Int_t tr = (t-10957*86400)%86400;
57     //cout<<"tr = "<<tr<<endl;
58     Double_t Somg = (Se+49.077+omg*86400*tr/360)*360/86400;
59     //cout<<"t1 = "<<(t-10957*86400)%86400<<"\n";
60     //cout<<"tr = "<<tr<<"\n";
61     //cout<<"Somg = "<<Se+omg*86400*tr/360<<"\n";
62     //cout<<"Somg = "<<((Somg-360*6)*86400/360/3600-20)*60<<"\n";
63     //cout<<cos(Somg/a)<<endl;
64     Gij(0,0) = cos(Somg/a);
65     Gij(0,1) = -sin(Somg/a);
66     Gij(0,2) = 0;
67     Gij(1,0) = sin(Somg/a);
68     Gij(1,1) = cos(Somg/a);
69     Gij(1,2) = 0;
70     Gij(2,0) = 0;
71     Gij(2,1) = 0;
72     Gij(2,2) = 1;
73     Gij.Invert();
74     //SetDirAxisGreenwich(Aij);
75     //cout<<(Somg/a)<<endl<<Aij(0,0)<<" "<<Aij(1,0)<<" "<<Aij(2,0)<<endl;
76     return Gij*Aij;
77     }
78    
79     TMatrixD PamelaOrientation::GreenwichtoGEO(Double_t lat, Double_t lon, TMatrixD Aij){
80     //Double_t a = 360/(2*TMath::Pi());
81     //Double_t Re = 6000000;
82     TMatrixD Gij(3,3);
83     TMatrixD Fij(3,3);
84    
85     TMatrixD Hij(3,3); //TEST
86     TMatrixD Iij(3,3); //TEST
87    
88     // if((lat<0.1)&&(lat>-0.1)){
89     //cout<<"lon = "<<lon<<" lat = "<<lat<<endl;
90     lon=(-lon)/a; lat=(-lat)/a;
91     //cout<<"lon = "<<lon<<" lat = "<<lat<<endl;
92     //
93     // cout<<"Quaternions Array"<<endl;
94     //cout<<Aij(0,0)<<" "<<Aij(0,1)<<" "<<Aij(0,2)<<endl;
95     //cout<<Aij(1,0)<<" "<<Aij(1,1)<<" "<<Aij(1,2)<<endl;
96     //cout<<Aij(2,0)<<" "<<Aij(2,1)<<" "<<Aij(2,2)<<endl<<endl;
97     // }
98     //Double_t x0 = (alt+Re)*sin(lat)*sin(lon);
99     //Double_t y0 = (alt+Re)*sin(lat)*cos(lon);
100     //Double_t Sa = lon-a*atan(y0/x0);
101     //if (y0>0&&x0<0) Sa=-Sa+90;
102     //if (y0<0&&x0<0) Sa=Sa-90;
103     //if (y0>0&&x0==0) Sa=90;
104     //if (y0<0&&x0==0) Sa=-90;
105    
106     Gij(0,0) = cos(lon);
107     Gij(0,1) = -sin(lon);
108     Gij(0,2) = 0;
109     Gij(1,0) = sin(lon);
110     Gij(1,1) = cos(lon);
111     Gij(1,2) = 0;
112     Gij(2,0) = 0;
113     Gij(2,1) = 0;
114     Gij(2,2) = 1;
115    
116     //cout<<"First rotation"<<endl;
117     //cout<<Gij(0,0)<<" "<<Gij(0,1)<<" "<<Gij(0,2)<<endl;
118     //cout<<Gij(1,0)<<" "<<Gij(1,1)<<" "<<Gij(1,2)<<endl;
119     //cout<<Gij(2,0)<<" "<<Gij(2,1)<<" "<<Gij(2,2)<<endl<<endl;
120    
121     //Gij.Invert();
122    
123     Fij(0,0) = cos(lat);
124     Fij(0,1) = 0;
125     Fij(0,2) = -sin(lat);
126     Fij(1,0) = 0;
127     Fij(1,1) = 1;
128     Fij(1,2) = 0;
129     Fij(2,0) = sin(lat);
130     Fij(2,1) = 0;
131     Fij(2,2) = cos(lat);
132    
133     //Fij.Invert();
134    
135     //if((lat<0.1)&&(lat>-0.1)){
136     /* Hij=Gij*Aij; //TEST
137    
138     cout<<"First rotation"<<endl;
139     cout<<Hij(0,0)<<" "<<Hij(0,1)<<" "<<Hij(0,2)<<endl;
140     cout<<Hij(1,0)<<" "<<Hij(1,1)<<" "<<Hij(1,2)<<endl;
141     cout<<Hij(2,0)<<" "<<Hij(2,1)<<" "<<Hij(2,2)<<endl<<endl;
142    
143     Iij = Fij*(Gij*Aij); //TEST
144    
145     cout<<"Second rotation"<<endl;
146     cout<<Iij(0,0)<<" "<<Iij(0,1)<<" "<<Iij(0,2)<<endl;
147     cout<<Iij(1,0)<<" "<<Iij(1,1)<<" "<<Iij(1,2)<<endl;
148     cout<<Iij(2,0)<<" "<<Iij(2,1)<<" "<<Iij(2,2)<<endl;
149     //
150     Int_t ret;
151     cin>>ret;*/
152     // }
153     return Fij*(Gij*Aij);
154     }
155    
156     TMatrixD PamelaOrientation::PamelatoGEO(TMatrixD Aij, Double_t B1, Double_t B2, Double_t B3){
157     //TMatrixD Gij(3,3);
158     TMatrixD Hij(3,1);
159     TMatrixD Bij(3,1);
160     Bij(0,0) = B1;
161     Bij(1,0) = B2;
162     Bij(2,0) = B3;
163     //Double_t alfa = TMath::ASin(sqrt(1/((Aij(1,2))/Aij(0,2)+1))) * TMath::RadToDeg();
164     //Gij(0,0) = cos(alfa/a);
165     //Gij(0,1) = -sin(alfa/a);
166     //Gij(0,2) = 0;
167     //Gij(1,0) = 0;
168     //Gij(1,1) = 1;
169     //Gij(1,2) = 0;
170     //Gij(2,0) = sin(alfa/a);
171     //Gij(2,1) = cos(alfa/a);
172     //Gij(2,2) = 0;
173    
174     Hij=Aij*Bij;
175     return Hij;
176     //cout<<0.25-Aij(2,2)/(Aij(2,1)*Aij(2,0))<<endl;
177     //cout<<Hij(0,0)<<endl;//" "<<Hij(0,1)<<" "<<Hij(0,2)<<endl;
178     //cout<<Hij(1,0)<<endl;//" "<<Hij(1,1)<<" "<<Hij(1,2)<<endl;
179     //cout<<Hij(2,0)<<endl;//" "<<Hij(2,1)<<" "<<Hij(2,2)<<endl;
180     }
181    
182     TMatrixD PamelaOrientation::ColPermutation(TMatrixD Aij){
183     TMatrixD Gij(3,3);
184     Gij(0,0) = 1; Gij(0,1) = 0; Gij(0,2) = 0;
185     Gij(1,0) = 0; Gij(1,1) = 0; Gij(1,2) = 1;
186     Gij(2,0) = 0; Gij(2,1) = -1; Gij(2,2) = 0;
187     return Aij*Gij;
188     }
189    
190     Double_t PamelaOrientation::GetPitchAngle(Double_t x1, Double_t y1, Double_t z1, Double_t x2, Double_t y2, Double_t z2){
191     return TMath::ACos((x1*x2 + y1*y2 + z1*z2)/(sqrt(pow(x1,2)+pow(y1,2)+pow(z1,2))*sqrt(pow(x2,2)+pow(y2,2)+pow(z2,2)))) * TMath::RadToDeg();
192     }
  ViewVC Help
Powered by ViewVC 1.1.23