001 /*
002 * $Id: Kepler.java,v 1.3 2012/02/19 17:35:39 davep Exp $
003 *
004 * This file is part of McIDAS-V
005 *
006 * Copyright 2007-2012
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 */
052 package edu.wisc.ssec.mcidasv.data.adde.sgp4;
053
054 //import jsattrak.utilities.StateVector;
055
056 public 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