--- quicklook/SatelliteInclination/src/InclinationInfo.cpp 2007/02/08 00:49:35 1.1.1.1 +++ quicklook/SatelliteInclination/src/InclinationInfo.cpp 2007/04/17 07:46:48 1.2 @@ -1,544 +1,162 @@ -/*************************************************************************** - * 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 -#include -#include "TMatrixD.h" - -using namespace std; - -Quaternions::Quaternions() - : InclinationInfoItem() -{ -} - - -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; -} - -void InclinationInfo::QuaternionstoAngle(Quaternions Qua){ - - double_t a11 = pow(Qua.quat[0][0],2.)+pow(Qua.quat[0][1],2.)-pow(Qua.quat[0][2],2.)-pow(Qua.quat[0][3],2.); - double_t a12 = 2*(Qua.quat[0][1]*Qua.quat[0][2]+Qua.quat[0][0]*Qua.quat[0][3]); - double_t a13 = 2*(Qua.quat[0][1]*Qua.quat[0][3]-Qua.quat[0][0]*Qua.quat[0][2]); - double_t a21 = 2*(Qua.quat[0][1]*Qua.quat[0][2]-Qua.quat[0][0]*Qua.quat[0][3]); - double_t a22 = pow(Qua.quat[0][0],2.)+pow(Qua.quat[0][2],2.)-pow(Qua.quat[0][1],2.)-pow(Qua.quat[0][3],2.); - double_t a23 = 2*(Qua.quat[0][2]*Qua.quat[0][3]+Qua.quat[0][0]*Qua.quat[0][1]); - double_t a31 = 2*(Qua.quat[0][1]*Qua.quat[0][3]+Qua.quat[0][0]*Qua.quat[0][2]); - double_t a32 = 2*(Qua.quat[0][2]*Qua.quat[0][3]-Qua.quat[0][0]*Qua.quat[0][1]); - double_t a33 = pow(Qua.quat[0][0],2.)+pow(Qua.quat[0][3],2.)-pow(Qua.quat[0][1],2.)-pow(Qua.quat[0][2],2.); - double_t a = 360/(2*TMath::Pi()); - double_t eksi = 0.0000001; - double_t eteta = 0.0000001; - double_t ksteta = a22*a22/(a12*a12+a22*a22); - double_t ksksi = a33*a33/(a33*a33+a31*a31); - - Int_t kj1; - if (a33<0){kj1=1; - } else {kj1=0;}; - Int_t kj2; - if (ksksi>eksi){kj2=1; - } else {kj2=0;}; - Int_t kj3; - if (ksksi<=eksi){kj3=1; - } else {kj3=0;}; - Int_t kj4; - if (a22<0){kj4=1; - } else {kj4=0;}; - Int_t kj5; - if (ksteta>eteta){kj5=1; - } else {kj5=0;}; - Int_t kj6; - if (ksteta<=eteta){kj6=1; - } else {kj6=0;}; - if (abs((int)a32)>1){exit(1);}; - Int_t fr; - - Double_t gamar = -atan(a32/sqrt(1-pow(a32,2.))); - Double_t ksir = (-atan(a31/a33)-TMath::Pi()*kj1*Sign_1(a31, fr))*kj2-0.5*TMath::Pi()*kj3*Sign_1(a31, fr); - Double_t tetar = -(-atan(a12/a22)-TMath::Pi()*kj4*Sign_1(a12, fr))*kj5+0.5*TMath::Pi()*kj6*Sign_1(a12, fr); -// if (gamar<0){A11=gamar*a+360;}else{A11=gamar*a;}; -// if (ksir<0){A11=ksir*a+360;}else{A11=ksir*a;}; -// if (tetar<0){A13=tetar*a+360;}else{A13=tetar*a;}; - -// gamar = acos(pow(Qua.quat[0][0],2.)+pow(Qua.quat[0][3],2.)-pow(Qua.quat[0][1],2.)-pow(Qua.quat[0][2],2.)); -// tetar = atan((Qua.quat[0][2]*Qua.quat[0][3]+Qua.quat[0][0]*Qua.quat[0][2])/(Qua.quat[0][2]*Qua.quat[0][3]+Qua.quat[0][1]*Qua.quat[0][0])); -// ksir = atan((Qua.quat[0][1]*Qua.quat[0][3]-Qua.quat[0][0]*Qua.quat[0][2])/(Qua.quat[0][2]*Qua.quat[0][3]-Qua.quat[0][1]*Qua.quat[0][0])); - - - A13=tetar*a; - A12=ksir*a; - A11=gamar*a; - -return ; -} - -/******************************************************************************************************************/ -/******************************************************************************************************************/ -//********************* ***************************************************************/ -//********************* COORDINATE SYSTEMS ***************************************************************/ -//********************* ***************************************************************/ -//*****************************************************************************************************************/ -//*****************************************************************************************************************/ -// -// ZISK -// + -// / \ YOSK ZOSK (Directed by Radius) -// | _ _. -// | |\ /| -// | \ / -// | \ / -// |.__..__ \ / -// Orbit _._.***| **.\/_ XOSK (Directed by velocity) -// .* | (X0,Y0,Z0) **--.___\ -// _** | / *. / -// .* | * * -// * ..****|***.. / R * -// .* | .*. -// .* | / *. -// * EARTH | / * YISK -// * | /_ _ _*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _\ -// * / * / -// * / .* -// *. / .* -// **/******* -// / -// / -// / -// / -// / -// / -// |/ -// *-- -// XISK -// -//****************************************************************************************************/ -//****************************************************************************************************/ - -//void OrbitalInfo::coefplane(Double_t x1,Double_t y1,Double_t z1,Double_t x2, Double_t y2, Double_t z2){ -// k1 = ((z1/y1)*((x2*y1 - x1*y2)/(z2*y1 - z1*y2)) - x1/y1); -// k2 = (x1*y2 - x2*y1)/(z2*y1 - z1*y2); -// } - -//Double_t OrbitalInfo::AngBetAxes(Double_t x1,Double_t y1,Double_t z1,Double_t x2, Double_t y2, Double_t z2){ -// return 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))); -// } - -//Double_t OrbitalInfo::ValueT(Double_t x1, Double_t y1, Double_t z1, Double_t n1, Double_t n2){ -// return -(x1+n1*y1+n2*z1)/(1+pow(n1,2)+pow(n2,2)); -// } - -//Double_t OrbitalInfo::AngBetPlan(Double_t x1, Double_t y1, Double_t z1, Double_t n1, Double_t n2){ -// return asin((x1+n1*y1+n2*z1)/(sqrt(pow(x1,2)+pow(y1,2)+pow(z1,2))*sqrt(1+pow(n1,2)+pow(n2,2)))); -// } - -void InclinationInfo::TransAngle(Double_t x0, Double_t y0, Double_t z0, Double_t Vx0, Double_t Vy0, Double_t Vz0, Double_t gamar, Double_t ksir, Double_t tetar, Double_t q0, Double_t q1, Double_t q2, Double_t q3){ - - double_t a = 360/(2*TMath::Pi()); - -// Points on three axes of Resurs' coordinate system (RCS) - Int_t XRCS[3]; Int_t YRCS[3]; Int_t ZRCS[3]; - -// Angles between our Axes(RCS) and planes of Orbital Coordinate System (OCS); -// Double_t AboAa0ZX[3]; -// Double_t AboAa0XY[3]; -// Double_t AboAa0YZ[3]; - -// Angles between our Axes(RCS) and Axes of OCS -// Double_t AboA0X[3]; -// Double_t AboA0Y[3]; -// Double_t AboA0Z[3]; - -//Angles between Proection of our axes on every plane of OCS and axes of it plane. -// Double_t AbPoAaAoP0ZX[3]; -// Double_t AbPoAaAoP0XY[3]; -// Double_t AbPoAaAoP0YZ[3]; - - XRCS[0] = 1; YRCS[0] = 0; ZRCS[0] = 0; // Points on X-axis RCS. - XRCS[1] = 0; YRCS[1] = 1; ZRCS[1] = 0; // Points on Y-axis RCS. - XRCS[2] = 0; YRCS[2] = 0; ZRCS[2] = 1; // Points on Z-axis - -// Transition matrix RCS -> Inertial Coordinate System (ICS) - TMatrixD Bij(3,3); - Bij(0,0) = cos(tetar)*cos(ksir)-sin(tetar)*sin(gamar)*sin(ksir); - Bij(0,1) = -sin(tetar)*cos(gamar); - Bij(0,2) = cos(tetar)*sin(ksir)+sin(tetar)*sin(gamar)*cos(ksir); - Bij(1,0) = sin(tetar)*cos(ksir)+cos(tetar)*sin(gamar)*sin(ksir); - Bij(1,1) = cos(tetar)*cos(gamar); - Bij(1,2) = sin(tetar)*sin(ksir)-cos(tetar)*sin(gamar)*cos(ksir); - Bij(2,0) = -sin(ksir)*cos(gamar); - Bij(2,1) = sin(gamar); - Bij(2,2) = cos(ksir)*cos(gamar); - -//********************************************************************************************/ -//*********************** OTHER METHOD OF GETING COORDINATE IN OCS****************************/ -//********************************************************************************************/ - - 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); - Pij(0,2) = 2*(q1*q3-q0*q2); - Pij(1,0) = 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); - Pij(2,0) = 2*(q1*q3+q0*q2); - Pij(2,1) = 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 AbPoAaAoP0ZX_2m[3]; Double_t AboAa0ZX_2m[3]; -// Double_t AbPoAaAoP0XY_2m[3]; Double_t AboAa0XY_2m[3]; -// Double_t AbPoAaAoP0YZ_2m[3]; Double_t AboAa0YZ_2m[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)); - //cout<<"C1= "<<(Vy0*C3-Vz0*C2)/(V0*C)<<", C2= "<<(Vz0*C1-Vx0*C3)/(V0*C)<<", C3="<<(Vx0*C2-Vy0*C1)/(V0*C)<<"\n"; - 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(); -// Double_t Xnew = Aij(0,0)*(X-x0)+Aij(0,1)*(Y-y0)+Aij(0,2)*(Z-z0); -// Double_t Ynew = Aij(1,0)*(X-x0)+Aij(1,1)*(Y-y0)+Aij(1,2)*(Z-z0); -// Double_t Znew = Aij(2,0)*(X-x0)+Aij(2,1)*(Y-y0)+Aij(2,2)*(Z-z0); - -/*********************************************************************************************/ - - Double_t Azim = atan(R0*C3/(y0*C1-x0*C2)); - Double_t Sa = sin(Azim); Double_t Ca = cos(Azim); - Double_t R1 = sqrt(pow(x0,2)+pow(y0,2)); - Double_t Sb = z0/R0; Double_t Cb = R1/R0; - Double_t Sl = y0/R1; Double_t Cl = x0/R1; - - TMatrixD Tij(3,3); - Tij(0,0) = -Cl*Sb*Ca-Sa*Sl; - Tij(0,1) = Sa*Cl-Ca*Sl*Sb; - Tij(0,2) = Ca*Cb; - Tij(1,0) = Ca*Sl-Sa*Sb*Cl; - Tij(1,1) = -Sa*Sl*Sb-Ca*Cl; - Tij(1,2) = Sa*Cb; - Tij(2,0) = Cb*Cl; - Tij(2,1) = Cb*Sl; - Tij(2,2) = Sb; - //cout<<"Tij\n"; - //cout< x+k1*y+k2*z=0; -// Double_t k1y; -// Double_t k2y; - //cout<>uchu; - -// Next we must find equation of Y-axis of OCS. For it we must find equation of line passing through -// point (0,0,0) perpendicularly by 0ZX plane of OCS -// generaly equation is: -// x = x0 + At; -// y = y0 + Bt; -// z = z0 + Ct; -// But we have point (x0,y0,z0) is (0,0,0) and other plane equation. For us it's: -// x = t; -// y = (B/A)*t => y = (B/A)*x; or x/x1=y/y1=z/z1 where -// z = (C/A)*t z = (C/A)*x; y1=B/A,z1=C/A and x1 we must find - -// if ((YonX<0 && ZonZ>0)||(YonX>0 && ZonZ<0)) XonY = 1; -// if ((YonX>0 && ZonZ>0)||(YonX>0 && ZonZ<0)) XonY = 1; -// Double_t XonY = 1; Double_t YonY = -k1; Double_t ZonY = k2; - -// coefficients for equations of 0XY plane of OCS. -// coefplane(XonX,YonX,ZonX,XonY,YonY,ZonY); -// Double_t k1XY = k1; Double_t k2XY = k2; - //cout<<"P3= "<<(XonY+k1XY*YonY+k2XY*ZonY)<<"\n"; - //cout<<"P3= "<<(XonX+k1XY*YonX+k2XY*ZonX)<<"\n"; - //cout<<"k1XY= "<>yuip; - - for (Int_t i = 0; i<3; i++) { -// Values of points on axes of RCS in ICS - Double_t XICS = Pij(0,0)*XRCS[i] + Pij(0,1)*YRCS[i] + Pij(0,2)*ZRCS[i];// + x0; - Double_t YICS = Pij(1,0)*XRCS[i] + Pij(1,1)*YRCS[i] + Pij(1,2)*ZRCS[i];// + y0; - Double_t ZICS = Pij(2,0)*XRCS[i] + Pij(2,1)*YRCS[i] + Pij(2,2)*ZRCS[i];// + z0; - //cout<<"XICS= "<>oiu; - -// Angles between our Axis and Z,Y,X-axes of OCS -// AboA0Z[i] = AngBetAxes(XICS,YICS,ZICS,XonZ,YonZ,ZonZ); -// AboA0Y[i] = AngBetAxes(XICS,YICS,ZICS,XonY,YonY,ZonY); -// AboA0X[i] = AngBetAxes(XICS,YICS,ZICS,XonX,YonX,ZonX); - -//Find coordinate of our point in OCS -// Double_t XOCS; -// Double_t YOCS; -// Double_t ZOCS; -// Double_t T = ValueT(XICS,YICS,ZICS,k1y,k2y); -// Double_t XonXZ = XICS + T; -// Double_t YonXZ = YICS + k1y*T; -// Double_t ZonXZ = ZICS + k2y*T; -// Double_t R = T*sqrt(1+pow(k1y,2)+pow(k2y,2)); -// YOCS = R; - //cout<<"CHECK= "<>alsdj; - -//Find Angles between RCS-axes and OCS-planes; -// AboAa0ZX[i] = AngBetPlan(XICS,YICS,ZICS,k1y,k2y); -// AboAa0XY[i] = AngBetPlan(XICS,YICS,ZICS,k1XY,k2XY); -// AboAa0YZ[i] = AngBetPlan(XICS,YICS,ZICS,k1YZ,k2YZ); - //AbPoAaAoP0ZX[i] = atan(ZOCS/XOCS); AbPoAaAoP0ZX_2m[i] = atan(ZOCS_2m/XOCS_2m); - //AbPoAaAoP0XY[i] = atan(XOCS/YOCS); AbPoAaAoP0XY_2m[i] = atan(YOCS_2m/XOCS_2m); - //AbPoAaAoP0YZ[i] = atan(ZOCS/YOCS); AbPoAaAoP0YZ_2m[i] = atan(ZOCS_2m/YOCS_2m); - -// if (XOCS_2m>0) AbPoAaAoP0ZX_2m[i] = atan(ZOCS_2m/XOCS_2m); -// if ((XOCS_2m<0)&&(ZOCS_2m>=0)) AbPoAaAoP0ZX_2m[i] = 180/a + atan(ZOCS_2m/XOCS_2m); -// if ((XOCS_2m<0)&&(ZOCS_2m<0)) AbPoAaAoP0ZX_2m[i] = atan(ZOCS_2m/XOCS_2m) - 180/a; -// if ((XOCS_2m=0)&&(ZOCS_2m>0)) AbPoAaAoP0ZX_2m[i] = 90/a; -// if ((XOCS_2m=0)&&(ZOCS_2m<0)) AbPoAaAoP0ZX_2m[i] = -90/a; -// if ((XOCS_2m=0)&&(ZOCS_2m=0)) AbPoAaAoP0ZX_2m[i] = 0; - -// if (XOCS_2m>0) AbPoAaAoP0XY_2m[i] = atan(YOCS_2m/XOCS_2m); -// if ((XOCS_2m<0)&&(YOCS_2m>=0)) AbPoAaAoP0XY_2m[i] = 180/a + atan(YOCS_2m/XOCS_2m); -// if ((XOCS_2m<0)&&(YOCS_2m<0)) AbPoAaAoP0XY_2m[i] = atan(YOCS_2m/XOCS_2m) - 180/a; -// if ((XOCS_2m=0)&&(YOCS_2m>0)) AbPoAaAoP0XY_2m[i] = 90/a; -// if ((XOCS_2m=0)&&(YOCS_2m<0)) AbPoAaAoP0XY_2m[i] = -90/a; -// if ((XOCS_2m=0)&&(YOCS_2m=0)) AbPoAaAoP0XY_2m[i] = 0; - -// if (YOCS_2m>0) AbPoAaAoP0YZ_2m[i] = atan(ZOCS_2m/YOCS_2m); -// if ((YOCS_2m<0)&&(ZOCS_2m>=0)) AbPoAaAoP0YZ_2m[i] = 180/a + atan(ZOCS_2m/YOCS_2m); -// if ((YOCS_2m<0)&&(ZOCS_2m<0)) AbPoAaAoP0YZ_2m[i] = atan(ZOCS_2m/YOCS_2m) - 180/a; -// if ((YOCS_2m=0)&&(ZOCS_2m>0)) AbPoAaAoP0YZ_2m[i] = 90/a; -// if ((YOCS_2m=0)&&(ZOCS_2m<0)) AbPoAaAoP0YZ_2m[i] = -90/a; -// if ((YOCS_2m=0)&&(ZOCS_2m=0)) AbPoAaAoP0YZ_2m[i] = 0; - - //if (i==0) cout<<"AbPoAaAoP0ZX_2m[i]"<0) Ryskanie = a*atan(u21/u22); - if ((u22<0)&&(u21>=0)) Ryskanie = 180 + a*atan(u21/u22); - if ((u22<0)&&(u21<0)) Ryskanie = a*atan(u21/u22) - 180; - if ((u22=0)&&(u21>0)) Ryskanie = 90; - if ((u22=0)&&(u21<0)) Ryskanie = -90; - if ((u22=0)&&(u21=0)) Ryskanie = 0; - -// AXrXo = AboA0X[0]*a; AXrZXo = AboAa0ZX[0]*a; ApXrZXoZ = AbPoAaAoP0ZX[0]*a; -// AXrYo = AboA0Y[0]*a; AXrXYo = AboAa0XY[0]*a; ApXrXYoZ = AbPoAaAoP0XY[0]*a; -// AXrZo = AboA0Z[0]*a; AXrYZo = AboAa0YZ[0]*a; ApXrYZoZ = AbPoAaAoP0YZ[0]*a; -// AYrXo = AboA0X[1]*a; AYrZXo = AboAa0ZX[1]*a; ApYrZXoZ = AbPoAaAoP0ZX[1]*a; -// AYrYo = AboA0Y[1]*a; AYrXYo = AboAa0XY[1]*a; ApYrXYoZ = AbPoAaAoP0XY[1]*a; -// AYrZo = AboA0Z[1]*a; AYrYZo = AboAa0YZ[1]*a; ApYrYZoZ = AbPoAaAoP0YZ[1]*a; -// AZrXo = AboA0X[2]*a; AZrZXo = AboAa0ZX[2]*a; ApZrZXoZ = AbPoAaAoP0ZX[2]*a; -// AZrYo = AboA0Y[2]*a; AZrXYo = AboAa0XY[2]*a; ApZrXYoZ = AbPoAaAoP0XY[2]*a; -// AZrZo = AboA0Z[2]*a; AZrYZo = AboAa0YZ[2]*a; ApZrYZoZ = AbPoAaAoP0YZ[2]*a; - -// AXrZXo_2m = AboAa0ZX_2m[0]*a; ApXrZXoZ_2m = AbPoAaAoP0ZX_2m[0]*a; -// AXrXYo_2m = AboAa0XY_2m[0]*a; ApXrXYoZ_2m = AbPoAaAoP0XY_2m[0]*a; -// AXrYZo_2m = AboAa0YZ_2m[0]*a; ApXrYZoZ_2m = AbPoAaAoP0YZ_2m[0]*a; -// AYrZXo_2m = AboAa0ZX_2m[1]*a; ApYrZXoZ_2m = AbPoAaAoP0ZX_2m[1]*a; -// AYrXYo_2m = AboAa0XY_2m[1]*a; ApYrXYoZ_2m = AbPoAaAoP0XY_2m[1]*a; -// AYrYZo_2m = AboAa0YZ_2m[1]*a; ApYrYZoZ_2m = AbPoAaAoP0YZ_2m[1]*a; -// AZrZXo_2m = AboAa0ZX_2m[2]*a; ApZrZXoZ_2m = AbPoAaAoP0ZX_2m[2]*a; -// AZrXYo_2m = AboAa0XY_2m[2]*a; ApZrXYoZ_2m = AbPoAaAoP0XY_2m[2]*a; -// AZrYZo_2m = AboAa0YZ_2m[2]*a; ApZrYZoZ_2m = AbPoAaAoP0YZ_2m[2]*a; - -/* - //Int_t Y=2;Int_t X=2;Int_t Z=4;Int_t Vx=5;Int_t Vz=9;Int_t Vy=5; - - Double_t X[2]; Double_t Y[2]; Double_t Z[2]; Double_t Vx[2]; Double_t Vy[2]; Double_t Vz[2]; - - TMatrixD Aij(3,3); - TMatrixD Bij(3,3); - - Bij(0,0) = cos(tetar)*cos(ksir)-sin(tetar)*sin(gamar)*sin(ksir); - Bij(0,1) = -sin(tetar)*cos(gamar); - Bij(0,2) = cos(tetar)*sin(ksir)+sin(tetar)*sin(gamar)*cos(ksir); - Bij(1,0) = sin(tetar)*cos(ksir)+cos(tetar)*sin(gamar)*sin(ksir); - Bij(1,1) = cos(tetar)*cos(gamar); - Bij(1,2) = sin(tetar)*sin(ksir)-cos(tetar)*sin(gamar)*cos(ksir); - Bij(2,0) = -sin(ksir)*cos(gamar); - Bij(2,1) = sin(gamar); - Bij(2,2) = cos(ksir)*cos(gamar); - - Double_t C1 = Y[0]*Vz[0] - Z[0]*Vy[0]; - Double_t C2 = Z[0]*Vx[0] - X[0]*Vz[0]; - Double_t C3 = X[0]*Vy[0] - Y[0]*Vx[0]; - 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)); - Aij(0,0) = Vx0/V0; - Aij(0,1) = C1/C; - Aij(0,2) = (Vy0*C3-Vz0*C2)/(V0*C); - Aij(1,0) = Vy0/V0; - Aij(1,1) = C2/C; - Aij(1,2) = (Vz0*C1-Vx0*C3)/(V0*C); - Aij(2,0) = Vz0/V0; - Aij(2,1) = C3/C; - Aij(2,2) = (Vx0*C2-Vy0*C1)/(V0*C); - Aij.Invert(); - Double_t Xnew = Aij(0,0)*(X-x0)+Aij(0,1)*(Y-y0)+Aij(0,2)*(Z-z0); - Double_t Ynew = Aij(1,0)*(X-x0)+Aij(1,1)*(Y-y0)+Aij(1,2)*(Z-z0); - Double_t Znew = Aij(2,0)*(X-x0)+Aij(2,1)*(Y-y0)+Aij(2,2)*(Z-z0); - */ - //A21 = NewTetar; - //A22 = NewGamar; - //A23 = NewKsir; - -return ; -} - +/*************************************************************************** + * 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 +#include +#include "TMatrixD.h" + +using namespace std; + +Quaternions::Quaternions() + : InclinationInfoItem() +{ +} + + +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 ; +} +