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

Diff of /DarthVader/OrbitalInfo/src/OrientationInfo.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1 by mocchiut, Wed Oct 1 15:26:35 2008 UTC revision 1.2 by pam-mep, Tue Nov 15 09:31:28 2011 UTC
# Line 34  TMatrixD OrientationInfo::QuatoECI(Float Line 34  TMatrixD OrientationInfo::QuatoECI(Float
34  }  }
35    
36  TMatrixD OrientationInfo::ECItoGreenwich(TMatrixD Aij, UInt_t t){  TMatrixD OrientationInfo::ECItoGreenwich(TMatrixD Aij, UInt_t t){
     //t=1154304000+86400*365;  
37      TMatrixD Gij(3,3);      TMatrixD Gij(3,3);
38      Double_t omg = (7.292115e-5)*a; // Earth rotation velosity (Around polar axis);      Double_t omg = (7.292115e-5)*a;                     // Earth rotation velosity (Around polar axis);
     //Double_t t1 = 0;  
     //if(t<=1158883200) t1 = 1127347200+229.2732;                    //absTime at 22/09/2005 + diference between Solar midnight and Greenwich sidereal midnight  
     //if(t>1158883200&&t<=1190419200) t1 = 1158883200+172.3415;//absTime at 22/09/2006 + diference between Solar midnight and Greenwich sidereal midnight  
     //if(t>=1190419200&&t<1222041600) t1 = 1190419200+115.39;  //absTime at 22/09/2007 + diference between Solar midnight and Greenwich sidereal midnight  
     //if(t>=1222041600) t1 = 1222041600 + 294.9361;          //absTime at 22/09/2008 + diference between Solar midnight and Greenwich sidereal midnight  
     //UInt_t Nd = (t-t1)/86400;  
     //Int_t DifSuSt = Nd*236.55;  
39      Double_t d = (t-10957*86400-43200);                      //Number of day, passing from 01/01/2000 12:00:00 to t;      Double_t d = (t-10957*86400-43200);                      //Number of day, passing from 01/01/2000 12:00:00 to t;
40      d = d/86400;      d = d/86400;
     //d = t-da*86400+DifSuSt  
     //cout<<"t = "<<t<<"\n";  
     //cout<<"t - 2000y = "<<t-10957*86400-43200<<"\n";  
     //cout<<"d = "<<d<<"\n";  
     //Int_t tl = t%86400;    //!!!!!!!!!!!!!!!!!!!!!!!!  
41      Double_t T = d/36525;                                    //Number of Julian centuries;      Double_t T = d/36525;                                    //Number of Julian centuries;
42            
     //Double_t tl = t-t1-Nd*86400-DifSuSt;  
43      Double_t Se = 6*3600+41*60+236.555367908*d+0.093104*pow(T,2)-(6.2e-6)*pow(T,3);      Double_t Se = 6*3600+41*60+236.555367908*d+0.093104*pow(T,2)-(6.2e-6)*pow(T,3);
44      //cout<<"Se = "<<Se<<"\n";  
     //cout<<t<<endl<<d<<endl<<tl<<endl<<Se+omg*tl*86400/360<<endl;  
45      Int_t tr = (t-10957*86400)%86400;      Int_t tr = (t-10957*86400)%86400;
46      //cout<<"tr = "<<tr<<endl;  
47      Double_t Somg = (Se+49.077+omg*86400*tr/360)*360/86400;      Double_t Somg = (Se+49.077+omg*86400*tr/360)*360/86400;
48      //cout<<"t1 = "<<(t-10957*86400)%86400<<"\n";  
49      //cout<<"tr = "<<tr<<"\n";      //Somg = 25; //for test transition
50      //cout<<"Somg = "<<Se+omg*86400*tr/360<<"\n";  
     //cout<<"Somg = "<<((Somg-360*6)*86400/360/3600-20)*60<<"\n";  
     //cout<<cos(Somg/a)<<endl;  
51      Gij(0,0) = cos(Somg/a);      Gij(0,0) = cos(Somg/a);
52      Gij(0,1) = -sin(Somg/a);      Gij(0,1) = -sin(Somg/a);
53      Gij(0,2) = 0;      Gij(0,2) = 0;
# Line 75  TMatrixD OrientationInfo::ECItoGreenwich Line 58  TMatrixD OrientationInfo::ECItoGreenwich
58      Gij(2,1) = 0;      Gij(2,1) = 0;
59      Gij(2,2) = 1;      Gij(2,2) = 1;
60      Gij.Invert();      Gij.Invert();
     //SetDirAxisGreenwich(Aij);  
     //cout<<(Somg/a)<<endl<<Aij(0,0)<<" "<<Aij(1,0)<<" "<<Aij(2,0)<<endl;  
61      return Gij*Aij;      return Gij*Aij;
62  }  }
63    
64  TMatrixD OrientationInfo::GreenwichtoGEO(Double_t lat, Double_t lon, TMatrixD Aij){  TMatrixD OrientationInfo::GreenwichtoGEO(Double_t lat, Double_t lon, TMatrixD Aij){
65      //Double_t a = 360/(2*TMath::Pi());  
     //Double_t Re = 6000000;  
66      TMatrixD Gij(3,3);      TMatrixD Gij(3,3);
67      TMatrixD Fij(3,3);      TMatrixD Fij(3,3);
       
     TMatrixD Hij(3,3); //TEST  
     TMatrixD Iij(3,3); //TEST  
       
 //    if((lat<0.1)&&(lat>-0.1)){  
         //cout<<"lon = "<<lon<<" lat = "<<lat<<endl;  
         lon=(-lon)/a; lat=(-lat)/a;  
         //cout<<"lon = "<<lon<<" lat = "<<lat<<endl;  
 //      
 //        cout<<"Quaternions Array"<<endl;  
         //cout<<Aij(0,0)<<" "<<Aij(0,1)<<" "<<Aij(0,2)<<endl;  
         //cout<<Aij(1,0)<<" "<<Aij(1,1)<<" "<<Aij(1,2)<<endl;  
         //cout<<Aij(2,0)<<" "<<Aij(2,1)<<" "<<Aij(2,2)<<endl<<endl;  
 //    }  
     //Double_t x0 = (alt+Re)*sin(lat)*sin(lon);  
     //Double_t y0 = (alt+Re)*sin(lat)*cos(lon);  
     //Double_t Sa = lon-a*atan(y0/x0);  
     //if (y0>0&&x0<0) Sa=-Sa+90;  
     //if (y0<0&&x0<0) Sa=Sa-90;  
     //if (y0>0&&x0==0) Sa=90;  
     //if (y0<0&&x0==0) Sa=-90;  
68    
69      Gij(0,0) = cos(lon);      lon=(-lon)/a; lat=(-lat)/a;   // here has the same result as Gij.Invert() in ECItoGreenwich function
70    
71        Gij(0,0) = cos(lon);        // rotation around z-axis:
72      Gij(0,1) = -sin(lon);      Gij(0,1) = -sin(lon);
73      Gij(0,2) = 0;      Gij(0,2) = 0;               //      | cos(lon)      -sin(lon)       0|
74      Gij(1,0) = sin(lon);      Gij(1,0) = sin(lon);        //      | sin(lon)      cos(lon)        0|
75      Gij(1,1) = cos(lon);      Gij(1,1) = cos(lon);        //      |   0                   0       1|
76      Gij(1,2) = 0;      Gij(1,2) = 0;
77      Gij(2,0) = 0;      Gij(2,0) = 0;
78      Gij(2,1) = 0;      Gij(2,1) = 0;
79      Gij(2,2) = 1;      Gij(2,2) = 1;
   
     //cout<<"First rotation"<<endl;  
     //cout<<Gij(0,0)<<" "<<Gij(0,1)<<" "<<Gij(0,2)<<endl;  
     //cout<<Gij(1,0)<<" "<<Gij(1,1)<<" "<<Gij(1,2)<<endl;  
     //cout<<Gij(2,0)<<" "<<Gij(2,1)<<" "<<Gij(2,2)<<endl<<endl;  
       
     //Gij.Invert();  
80            
81      Fij(0,0) = cos(lat);      Fij(0,0) = cos(lat);        // rotation around y-axis at angle -lat, cause rotation around y from x to z axis is negative
82      Fij(0,1) = 0;      Fij(0,1) = 0;               //
83      Fij(0,2) = -sin(lat);      Fij(0,2) = -sin(lat);       //      |cos(-lat)      0       sin(-lat)|              |cos(lat)       0       -sin(lat)|
84      Fij(1,0) = 0;      Fij(1,0) = 0;               //      |       0       1               0 |     ==>     |       0       1               0 |
85      Fij(1,1) = 1;      Fij(1,1) = 1;               //      |-sin(-lat)     0       cos(-lat)|              |sin(lat)       0       cos(lat) |
86      Fij(1,2) = 0;      Fij(1,2) = 0;
87      Fij(2,0) = sin(lat);      Fij(2,0) = sin(lat);
88      Fij(2,1) = 0;      Fij(2,1) = 0;
89      Fij(2,2) = cos(lat);      Fij(2,2) = cos(lat);
90        
     //Fij.Invert();  
       
     //if((lat<0.1)&&(lat>-0.1)){  
     /*    Hij=Gij*Aij;  //TEST  
       
         cout<<"First rotation"<<endl;  
         cout<<Hij(0,0)<<" "<<Hij(0,1)<<" "<<Hij(0,2)<<endl;  
         cout<<Hij(1,0)<<" "<<Hij(1,1)<<" "<<Hij(1,2)<<endl;  
         cout<<Hij(2,0)<<" "<<Hij(2,1)<<" "<<Hij(2,2)<<endl<<endl;  
       
         Iij = Fij*(Gij*Aij); //TEST  
       
         cout<<"Second rotation"<<endl;  
         cout<<Iij(0,0)<<" "<<Iij(0,1)<<" "<<Iij(0,2)<<endl;  
         cout<<Iij(1,0)<<" "<<Iij(1,1)<<" "<<Iij(1,2)<<endl;  
         cout<<Iij(2,0)<<" "<<Iij(2,1)<<" "<<Iij(2,2)<<endl;  
 //      
         Int_t ret;  
         cin>>ret;*/  
 //    }  
91      return Fij*(Gij*Aij);      return Fij*(Gij*Aij);
92  }  }
93    
94    TMatrixD OrientationInfo::GEOtoGeomag(TMatrixD Aij,Double_t Bnorth, Double_t Beast, Double_t Bup){      //Geomagnetic geodetic reference frame
95            Double_t alpha = 0;
96            if(Beast==0. && Bnorth>0) alpha = 0; else
97            if(Beast==0. && Bnorth<0) alpha = 180.; else{
98                    if(Beast > 0) alpha = TMath::ATan(Bnorth/Beast)*TMath::RadToDeg() - 90.;
99                    if(Beast < 0) alpha = TMath::ATan(Bnorth/Beast)*TMath::RadToDeg() + 90.;
100            }
101            alpha = alpha*TMath::DegToRad();
102            Double_t beta = TMath::ATan(Bup/sqrt(pow(Bnorth,2)+pow(Beast,2)));
103            //if(Bup<0.0) beta = TMath::ATan(TMath::Abs(Bup/sqrt(pow(Bnorth,2)+pow(Beast,2))));
104            //if(Bup>0.0) beta = TMath::ATan(TMath::Abs(sqrt(pow(Bnorth,2)+pow(Beast,2))/Bup));
105            //cout<<"GEOtomag:alpha = "<<alpha*TMath::RadToDeg()<<"\tbeta = "<<beta*TMath::RadToDeg()<<endl;
106            TMatrixD Gij(3,3);
107            TMatrixD Fij(3,3);
108            Gij(0,0) = 1;                   //rotation around x-axis at angle alpha
109            Gij(0,1) = 0;
110            Gij(0,2) = 0;                   //      |1              0               0       |
111            Gij(1,0) = 0;                   //      |0      cos(alpha)      -sin(alpha)     |
112            Gij(1,1) = cos(alpha);  //      |0      sin(alpha)      cos(alpha)      |
113            Gij(1,2) = -sin(alpha);
114            Gij(2,0) = 0;
115            Gij(2,1) = sin(alpha);
116            Gij(2,2) = cos(alpha);
117            Gij.Invert();
118            Fij(0,0) = cos(beta);   //rotation around y-axis at angle beta
119            Fij(0,1) = 0;
120            Fij(0,2) = sin(beta);   //      |cos(beta)      0       sin(beta)|
121            Fij(1,0) = 0;                   //      |       0       1               0 |
122            Fij(1,1) = 1;                   //      |-sin(beta)     0       cos(beta)|
123            Fij(1,2) = 0;
124            Fij(2,0) = -sin(beta);
125            Fij(2,1) = 0;
126            Fij(2,2) = cos(beta);
127            Fij.Invert();
128            //Int_t tri;
129            //cin >> tri;
130            return Fij*(Gij*Aij);
131    }
132    
133  TMatrixD OrientationInfo::PamelatoGEO(TMatrixD Aij, Double_t B1, Double_t B2, Double_t B3){  TMatrixD OrientationInfo::PamelatoGEO(TMatrixD Aij, Double_t B1, Double_t B2, Double_t B3){
134      //TMatrixD Gij(3,3);      //TMatrixD Gij(3,3);
135      TMatrixD Hij(3,1);      TMatrixD Hij(3,1);
# Line 164  TMatrixD OrientationInfo::PamelatoGEO(TM Line 137  TMatrixD OrientationInfo::PamelatoGEO(TM
137      Bij(0,0) = B1;      Bij(0,0) = B1;
138      Bij(1,0) = B2;      Bij(1,0) = B2;
139      Bij(2,0) = B3;      Bij(2,0) = B3;
     //Double_t alfa = TMath::ASin(sqrt(1/((Aij(1,2))/Aij(0,2)+1))) * TMath::RadToDeg();  
     //Gij(0,0) = cos(alfa/a);                    
     //Gij(0,1) = -sin(alfa/a);                    
     //Gij(0,2) = 0;                              
     //Gij(1,0) = 0;                              
     //Gij(1,1) = 1;                              
     //Gij(1,2) = 0;                              
     //Gij(2,0) = sin(alfa/a);                    
     //Gij(2,1) = cos(alfa/a);                    
     //Gij(2,2) = 0;                              
       
140      Hij=Aij*Bij;      Hij=Aij*Bij;
141      return Hij;      return Hij;
     //cout<<0.25-Aij(2,2)/(Aij(2,1)*Aij(2,0))<<endl;  
     //cout<<Hij(0,0)<<endl;//" "<<Hij(0,1)<<" "<<Hij(0,2)<<endl;  
     //cout<<Hij(1,0)<<endl;//" "<<Hij(1,1)<<" "<<Hij(1,2)<<endl;  
     //cout<<Hij(2,0)<<endl;//" "<<Hij(2,1)<<" "<<Hij(2,2)<<endl;  
142  }  }
143    
144  TMatrixD OrientationInfo::ColPermutation(TMatrixD Aij){  TMatrixD OrientationInfo::ColPermutation(TMatrixD Aij){
# Line 191  TMatrixD OrientationInfo::ColPermutation Line 149  TMatrixD OrientationInfo::ColPermutation
149      return Aij*Gij;      return Aij*Gij;
150  }  }
151    
152    Float_t OrientationInfo::Larmor(Float_t Ek,Float_t Bm,Int_t iZ,Float_t xA){  //Ek in MeV, Bm in nT, Pitch-angle, rad
153        Float_t mp = 938.272029;    Float_t amu = 931.494043e0;
154        Float_t cc = 299792458.;
155        Float_t ee = 1.60217653e-19;
156        Float_t kg = 1.7826619e-30;
157        Float_t gam = (Ek+mp)/mp;
158        Float_t mm = mp*kg;
159        Float_t omega = iZ*ee*Bm*1e-9/(gam*mm);
160        Float_t larmor = 1e-3*sqrt(1e0-1e0/pow(gam,2))*cc/omega;
161        larmor = 1e-3*Ek*cc/omega; //Ek here is p or for onecharged particle R;
162        return larmor;
163    }
164    
165    TMatrixD OrientationInfo::GetDirectiontoGirocenter(Float_t R, Float_t Px, Float_t Py){
166        TMatrixD GirDir(3,1);
167        if(R>0){
168            GirDir(0,0) = Py;
169            GirDir(1,0) = -Px;
170        }else{
171            GirDir(0,0) = -Py;
172            GirDir(1,0) = Px;
173        }
174        GirDir(2,0) = 0.;
175        return GirDir;
176    }
177    
178  Double_t OrientationInfo::GetPitchAngle(Double_t x1, Double_t y1, Double_t z1, Double_t x2, Double_t y2, Double_t z2){  Double_t OrientationInfo::GetPitchAngle(Double_t x1, Double_t y1, Double_t z1, Double_t x2, Double_t y2, Double_t z2){
179      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();      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();
180  }  }
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.23