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

Annotation of /ParDirCalc/CheckGeom/TrAng.cpp

Parent Directory Parent Directory | Revision Log Revision Log

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