001/* 002 * $Id: Kepler.java,v 1.2 2011/03/24 16:06:33 davep Exp $ 003 * 004 * This file is part of McIDAS-V 005 * 006 * Copyright 2007-2011 007 * Space Science and Engineering Center (SSEC) 008 * University of Wisconsin - Madison 009 * 1225 W. Dayton Street, Madison, WI 53706, USA 010 * https://www.ssec.wisc.edu/mcidas 011 * 012 * All Rights Reserved 013 * 014 * McIDAS-V is built on Unidata's IDV and SSEC's VisAD libraries, and 015 * some McIDAS-V source code is based on IDV and VisAD source code. 016 * 017 * McIDAS-V is free software; you can redistribute it and/or modify 018 * it under the terms of the GNU Lesser Public License as published by 019 * the Free Software Foundation; either version 3 of the License, or 020 * (at your option) any later version. 021 * 022 * McIDAS-V is distributed in the hope that it will be useful, 023 * but WITHOUT ANY WARRANTY; without even the implied warranty of 024 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 025 * GNU Lesser Public License for more details. 026 * 027 * You should have received a copy of the GNU Lesser Public License 028 * along with this program. If not, see http://www.gnu.org/licenses. 029 */ 030// Shawn E. Gano 031// Kepler routines 032/** 033 * ===================================================================== 034 * Copyright (C) 2009 Shawn E. Gano 035 * 036 * This file is part of JSatTrak. 037 * 038 * JSatTrak is free software: you can redistribute it and/or modify 039 * it under the terms of the GNU Lesser General Public License as published by 040 * the Free Software Foundation, either version 3 of the License, or 041 * (at your option) any later version. 042 * 043 * JSatTrak is distributed in the hope that it will be useful, 044 * but WITHOUT ANY WARRANTY; without even the implied warranty of 045 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 046 * GNU Lesser General Public License for more details. 047 * 048 * You should have received a copy of the GNU Lesser General Public License 049 * along with JSatTrak. If not, see <http://www.gnu.org/licenses/>. 050 * ===================================================================== 051 */ 052package edu.wisc.ssec.mcidasv.data.adde.sgp4; 053 054//import jsattrak.utilities.StateVector; 055 056public class Kepler 057{ 058 059 /** 060 * Overload for state 061 * @param GM Gravitational coefficient (gravitational constant * mass of central body) 062 * @param kepElements Keplarian elements (see other state function) 063 * @param dt Time since epoch (seconds?) 064 * @return State vector (x,y,z,vx,vy,vz) 065 */ 066/* 067 public static double[] state(double GM, double[] kepElements, double dt) 068 { 069 return state(GM, kepElements[0], kepElements[1],kepElements[2], kepElements[3], kepElements[4], kepElements[5], dt); 070 } 071*/ 072 /** 073 * Computes the satellite state vector from osculating Keplerian elements for elliptic orbits 074 * <p>Notes: 075 * <br> The semimajor axis a=Kep(0), dt and GM must be given in consistent units, 076 * e.g. [m], [s] and [m^3/s^2]. The resulting units of length and velocity 077 * are implied by the units of GM, e.g. [m] and [m/s]. 078 * 079 * @param GM Gravitational coefficient (gravitational constant * mass of central body) 080 * @param a Semimajor axis 081 * @param e Eccentricity 082 * @param i Inclination [rad] 083 * @param Omega Longitude of the ascending node [rad] or is this RAAN? (I think it is RAAN) 084 * @param omega Argument of pericenter [rad] 085 * @param M0 Mean anomaly at epoch [rad] 086 * @param dt Time since epoch (seconds?) 087 * @return State vector (x,y,z,vx,vy,vz) 088 */ 089/* 090 public static double[] state(double GM, double a, double e, double i, double Omega, double omega, double M0, double dt) 091 { 092 double M; // mean anomaly 093 double n; 094 double E, cosE, sinE; // eccentric anomaly 095 double fac, R, V; 096 097 // vectors for position (r) and velocity (v) 098 double[] r = new double[3]; 099 double[] v = new double[3]; 100 double[] state = new double[6]; // full state 101 102 // transformation matrices 103 double[][] Rx = new double[3][3]; 104 // double[][] Ry = new double[3][3]; 105 double[][] Rz = new double[3][3]; 106 107 // rotation matrix 108 double[][] PQW = new double[3][3]; 109 110 // Mean anomaly 111 if (dt == 0.0) 112 { 113 M = M0; 114 } else 115 { 116 n = Math.sqrt(GM / (a * a * a)); 117 M = M0 + n * dt; 118 } 119 120 // Eccentric anomaly 121 E = EccAnom(M, e); 122 cosE = Math.cos(E); 123 sinE = Math.sin(E); 124 125 // Perifocal coordinates 126 fac = Math.sqrt((1.0 - e) * (1.0 + e)); 127 128 R = a * (1.0 - e * cosE); // Distance 129 V = Math.sqrt(GM * a) / R; // Velocity 130 131 // r 132 r[0] = a * (cosE - e); 133 r[1] = a * fac * sinE; 134 r[2] = 0.0; 135 136 // v 137 v[0] = -V * sinE; 138 v[1] = +V * fac * cosE; 139 v[2] = 0.0; 140 141 // Transformation to reference system (Gaussian vectors) 142 Rx = MathUtils.R_x(-i); 143 Rz = MathUtils.R_z(-Omega); 144 145 // PQW = R_z(-Omega) * R_x(-i) * R_z(-omega); 146 PQW = MathUtils.mult(Rz, Rx); 147 Rz = MathUtils.R_z(-omega); 148 PQW = MathUtils.mult(PQW, Rz); 149 150 r = MathUtils.mult(PQW, r); 151 v = MathUtils.mult(PQW, v); 152 153 // State vector 154 state[0] = r[0]; 155 state[1] = r[1]; 156 state[2] = r[2]; 157 state[3] = v[0]; 158 state[4] = v[1]; 159 state[5] = v[2]; 160 161 // return Stack(r,v); 162 return state; 163 164 } // state 165*/ 166 167 /** 168 * Computes the eccentric anomaly for elliptic orbits 169 * 170 * @param M Mean anomaly in [rad] 171 * @param e Eccentricity of the orbit [0,1] 172 * @return Eccentric anomaly in [rad] 173 */ 174/* 175 public static double EccAnom(double M, double e) 176 { 177 178 // Constants 179 final int maxit = 15; 180 final double eps = 100.0 * 2.22E-16; // JAVA FIND MACHINE PRECISION // 100.0*eps_mach 181 182 // Variables 183 184 int i = 0; 185 double E, f; 186 187 // Starting value 188 M = MathUtils.Modulo(M, 2.0 * Math.PI); 189 if (e < 0.8) 190 E = M; 191 else 192 E = Math.PI; 193 194 // Iteration 195 do 196 { 197 f = E - e * Math.sin(E) - M; 198 E = E - f / (1.0 - e * Math.cos(E)); 199 ++i; 200 if (i == maxit) 201 { 202 System.out.println(" convergence problems in EccAnom\n"); 203 break; 204 } 205 } while (Math.abs(f) > eps); 206 207 return E; 208 209 } // EccAnom 210*/ 211 212 /** 213 * Computes the osculating Keplerian elements from the satellite state vector for elliptic orbits for Earth 214 * @param StateVector j2K cartesian state vector object (t,x,y,z,dx,dy,dz) 215 * @return Keplerian elements (a,e,i,Omega,omega,M) see state function for element definitions 216 */ 217/* 218 public static double[] SingularOsculatingElementsEarth( StateVector state ) 219 { 220 double[] r = new double[] {state.state[1],state.state[2],state.state[3]}; 221 double[] v = new double[] {state.state[4],state.state[5],state.state[6]}; 222 223 return SingularOsculatingElements( AstroConst.GM_Earth, r, v ); 224 } 225*/ 226 /** 227 * Computes the osculating Keplerian elements from the satellite state vector for elliptic orbits 228 * @param GM Gravitational coefficient (gravitational constant * mass of central body) 229 * @param state j2K cartesian state vector (x,y,z,dx,dy,dz) 230 * @return Keplerian elements (a,e,i,Omega,omega,M) see state function for element definitions 231 */ 232/* 233 public static double[] SingularOsculatingElements( double GM, double[] state ) 234 { 235 double[] r = new double[] {state[0],state[1],state[2]}; 236 double[] v = new double[] {state[3],state[4],state[5]}; 237 238 return SingularOsculatingElements( GM, r, v ); 239 } 240*/ 241 /** 242 * Computes the osculating Keplerian elements from the satellite state vector for elliptic orbits 243 * 244 * @param GM Gravitational coefficient (gravitational constant * mass of central body) 245 * @param r State vector (x,y,z) 246 * @param v State vector (vx,vy,vz) 247 * @return Keplerian elements (a,e,i,Omega,omega,M) see state function for element definitions 248 */ 249/* 250 public static double[] SingularOsculatingElements( double GM, double[] r, double[] v ) 251 { 252 // Variables 253 254 double[] h = new double[3]; 255 double H, u, R; 256 double eCosE, eSinE, e2, E, nu; 257 double a,e,i,Omega,omega,M; 258 259 h = MathUtils.cross(r,v); // Areal velocity 260 H = MathUtils.norm(h); 261 262 Omega = Math.atan2( h[0], -h[1] ); // Long. ascend. node 263 Omega = Omega % (Math.PI*2.0); 264 i = Math.atan2( Math.sqrt(h[0]*h[0]+h[1]*h[1]), h[2] ); // Inclination 265 u = Math.atan2( r[2]*H, -r[0]*h[1]+r[1]*h[0] ); // Arg. of latitude 266 267 R = MathUtils.norm(r); // Distance 268 269 //System.out.println("R="+R); 270 //System.out.println("v dot v="+(MathUtils.dot(v,v))); 271 //System.out.println("GM="+GM); 272 273 a = 1.0 / (2.0/R-MathUtils.dot(v,v)/GM); // Semi-major axis 274 275 eCosE = 1.0-R/a; // e*cos(E) 276 eSinE = MathUtils.dot(r,v)/Math.sqrt(GM*a); // e*sin(E) 277 278 e2 = eCosE*eCosE +eSinE*eSinE; 279 e = Math.sqrt(e2); // Eccentricity 280 E = Math.atan2(eSinE,eCosE); // Eccentric anomaly 281 282 M = (E-eSinE)%(2.0*Math.PI); // Mean anomaly 283 284 nu = Math.atan2(Math.sqrt(1.0-e2)*eSinE, eCosE-e2); // True anomaly 285 286 omega = (u-nu)%(2.0*Math.PI); // Arg. of perihelion 287 288 // Keplerian elements vector 289 290 //System.out.println("Real a = " + a); 291 292 return (new double[] {a,e,i,Omega,omega,M}); 293 } // SingularElements 294*/ 295 296 /** 297 * calculate oculating orbital period of a eliptical orbit from position and velocity 298 * 299 * @param GM Gravitational coefficient (gravitational constant * mass of central body) 300 * @param r State vector (x,y,z) 301 * @param v State vector (vx,vy,vz) 302 * @return oculating orbital period 303 */ 304 public static double CalculatePeriod( double GM, double[] r, double[] v ) 305 { 306 double R = MathUtils.norm(r); // current radius 307 308 double a = 1.0 / (2.0/R-MathUtils.dot(v,v)/GM); // Semi-major axis 309 310 return ( 2.0*Math.PI*Math.sqrt(a*a*a/GM) ); 311 } // CalculatePeriod 312 313} // Kepler