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