/*************************************************************************** * Copyright (C) 2006 by pamelaprod * * pamelaprod@P1.pamela * * * * This program is free software; you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation; either version 2 of the License, or * * (at your option) any later version. * * * * This program is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * GNU General Public License for more details. * * * * You should have received a copy of the GNU General Public License * * along with this program; if not, write to the * * Free Software Foundation, Inc., * * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * ***************************************************************************/ #include using namespace std; // InclinationInfoI()::InclinationInfoI() { // // memset(time,0,6*sizeof(double)); // // memset(quad,0,6*4*sizeof(double)); // }; void InclinationInfoI::fill(TArrayC* data){ short extIndex = 0; short innIndex = 0; long tempData = 0; for (int i = 0; i < 6; i++){ extIndex = 20*i; time[i] = (((data->At(extIndex) << 24) & 0xFF000000) + ((data->At(extIndex + 1) << 16) & 0x00FF0000) + ((data->At(extIndex + 2) << 8) & 0x0000FF00) + (data->At(extIndex + 3) & 0x000000FF))/128.0; for (int j = 0; j < 4; j++){ innIndex = extIndex + 4*j; tempData = ((data->At(innIndex + 4) << 24) & 0xFF000000) + ((data->At(innIndex + 5) << 16) & 0x00FF0000) + ((data->At(innIndex + 6) << 8) & 0x0000FF00) + (data->At(innIndex + 7) & 0x000000FF); if (data->At(innIndex + 4) >> 8) { quat[i][j] = (~tempData * -1.0)/1073741824.0; } else { quat[i][j] = tempData / 1073741824.0; } } } } void InclinationInfoI::clear() { for(UInt_t i = 0; i < 6; i++){ time[i]=0; for(UInt_t j = 0; j < 4; j++) quat[i][j]=0; } return ; } Quaternions::Quaternions() : InclinationInfoI() { } Quaternions::~Quaternions() { } InclinationInfo::InclinationInfo() : TObject() { } InclinationInfo::~InclinationInfo() { } short int Sign_1(double_t a, Int_t b){ if(a>0){b=1;} if(a<0){b=-1;} else{b=0;} return b; } /******************************************************************************************************************/ /******************************************************************************************************************/ //********************* ***************************************************************/ //********************* COORDINATE SYSTEMS ***************************************************************/ //********************* ***************************************************************/ //*****************************************************************************************************************/ //*****************************************************************************************************************/ // // ZISK // + // / \ YOSK ZOSK (Directed by Radius) // | _ _. // | |\ /| // | \ / // | \ / // |.__..__ \ / // Orbit _._.***| **.\/_ XOSK (Directed by velocity) // .* | (X0,Y0,Z0) **--.___| // _** | / *. / // .* | * * // * ..****|***.. / R * // .* | .*. // .* | / *. // * EARTH | / * YISK // * | /_ _ _*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _| // * / * / // * / .* // *. / .* // **/******* // / // / // / // / // / // / // |/ // *-- // XISK // //****************************************************************************************************/ //****************************************************************************************************/ 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){ double_t a = 360/(2*TMath::Pi()); TMatrixD Xij(3,3); Xij(0,0) = 1; Xij(0,1) = 0; Xij(0,2) = 0; Xij(1,0) = 0; Xij(1,1) = 0; Xij(1,2) = 1; Xij(2,0) = 0; Xij(2,1) = -1; Xij(2,2) = 0; TMatrixD Zij(3,3); Zij(0,0) = 0; Zij(0,1) = 0; Zij(0,2) = -1; Zij(1,0) = -1; Zij(1,1) = 0; Zij(1,2) = 0; Zij(2,0) = 0; Zij(2,1) = 1; Zij(2,2) = 0; TMatrixD Pij(3,3); Pij(0,0) = pow(q0,2)+pow(q1,2)-pow(q2,2)-pow(q3,2); Pij(0,1) = /*2*(q1*q2+q0*q3);/*/ 2*(q1*q2-q0*q3); Pij(0,2) = /*2*(q1*q3-q0*q2);/*/ 2*(q1*q3+q0*q2); Pij(1,0) = /*2*(q1*q2-q0*q3);/*/ 2*(q1*q2+q0*q3); Pij(1,1) = pow(q0,2)-pow(q1,2)+pow(q2,2)-pow(q3,2); Pij(1,2) = /*2*(q2*q3+q0*q1);/*/ 2*(q2*q3-q0*q1); Pij(2,0) = /*2*(q1*q3+q0*q2);/*/ 2*(q1*q3-q0*q2); Pij(2,1) = /*2*(q2*q3-q0*q1);/*/ 2*(q2*q3+q0*q1); Pij(2,2) = pow(q0,2)-pow(q1,2)-pow(q2,2)+pow(q3,2); TMatrixD Aij(3,3); Double_t C1 = y0*Vz0 - z0*Vy0; Double_t C2 = z0*Vx0 - x0*Vz0; Double_t C3 = x0*Vy0 - y0*Vx0; Double_t C = sqrt(pow(C1,2) + pow(C2,2) + pow(C3,2)); Double_t V0 = sqrt(pow(Vx0,2)+pow(Vy0,2) + pow(Vz0,2)); // Double_t R0 = sqrt(pow(x0,2)+pow(y0,2) + pow(z0,2)); Aij(0,0) = /*(C2*z0-C3*y0)/(C*R0);/*/Vx0/V0; Aij(0,1) = C1/C; Aij(0,2) = /*x0/R0;/*/(Vy0*C3-Vz0*C2)/(V0*C); Aij(1,0) = /*(C3*x0-C1*z0)/(C*R0);/*/Vy0/V0; Aij(1,1) = C2/C; Aij(1,2) = /*y0/R0;/*/(Vz0*C1-Vx0*C3)/(V0*C); Aij(2,0) = /*(C1*y0-C2*x0)/(C*R0);/*/Vz0/V0; Aij(2,1) = C3/C; Aij(2,2) = /*x0/R0;/*/(Vx0*C2-Vy0*C1)/(V0*C); Aij.Invert(); TMatrixD Full_(3,3); Full_ = Aij*(Pij*Zij); //Double_t u13 = Full_(0,2); //Double_t u23 = Full_(1,2); //Double_t u22 = Full_(1,1); //Double_t u33 = Full_(2,2); //Double_t u21 = Full_(1,0); Double_t u13 = Full_(0,0); Double_t u23 = -Full_(1,0); Double_t u22 = Full_(1,1); Double_t u33 = Full_(2,0); Double_t u21 = Full_(1,2); Tangazh = a*atan(-u13/u33); Kren = a*atan(-u23/sqrt(1 - pow(u23,2))); Ryskanie = a*atan(u21/u22); return ; } void InclinationInfo::Clear(Option_t *t){ //Int_t gyh = 0; } //ClassImp(McmdItem) ClassImp(InclinationInfoI) ClassImp(Quaternions) ClassImp(InclinationInfo)