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 /** ----------------------------------------------------------------
029 * Contains functions to read TLE data files and initalize the SGP4 propogator
030 * as well as other routines from sgp4ext.cpp
031 *
032 * <p>
033 * sgp4ext.cpp header information:
034 * <p>
035 *
036 * sgp4ext.cpp
037 *
038 * this file contains extra routines needed for the main test program for sgp4.
039 * these routines are derived from the astro libraries.
040 *
041 * companion code for
042 * fundamentals of astrodynamics and applications
043 * 2007
044 * by david vallado
045 *
046 * (w) 719-573-2600, email dvallado@agi.com
047 *
048 * current :
049 * 7 may 08 david vallado
050 * fix sgn
051 * changes :
052 * 2 apr 07 david vallado
053 * fix jday floor and str lengths
054 * updates for constants
055 * 14 aug 06 david vallado
056 * original baseline
057 * ---------------------------------------------------------------- */
058
059 package edu.wisc.ssec.mcidasv.data.adde.sgp4;
060
061 import java.text.DecimalFormat;
062 import java.text.DecimalFormatSymbols;
063 import java.text.ParsePosition;
064
065 /**
066 * 19 June 2009
067 * @author Shawn E. Gano, shawn@gano.name
068 */
069 public class SGP4utils
070 {
071 public static double pi = SGP4unit.pi;
072 public static char OPSMODE_AFSPC = 'a';
073 public static char OPSMODE_IMPROVED = 'i';
074
075 /**
076 * Reads the data from the TLE and initializes the SGP4 propogator variables and stores them in the SGP4unit.Gravconsttype object
077 * DOES NOT PERFORM ANY INTERNAL CHECK BEYOND BASICS OF THE TLE DATA use other methods to do that if desired.
078 *
079 * @param satName
080 * @param line1 TLE line 1
081 * @param line2 TLE line 2
082 * @param opsmode
083 * @param whichconst which constants to use in propogation
084 * @param satrec object to store the SGP4 data
085 * @return if the sgp4 propogator was initialized properly
086 */
087 public static boolean readTLEandIniSGP4(String satName, String line1, String line2, char opsmode, SGP4unit.Gravconsttype whichconst, SGP4SatData satrec)
088 {
089 final double deg2rad = pi / 180.0; // 0.0174532925199433
090 final double xpdotp = 1440.0 / (2.0 * pi); // 229.1831180523293
091
092 double sec, tumin;//, mu, radiusearthkm, xke, j2, j3, j4, j3oj2;
093 // double startsec, stopsec, startdayofyr, stopdayofyr, jdstart, jdstop;
094 // int startyear, stopyear, startmon, stopmon, startday, stopday,
095 // starthr, stophr, startmin, stopmin;
096 // int cardnumb, j; // numb,
097 //long revnum = 0, elnum = 0;
098 //char classification, intldesg[11], tmpstr[80];
099 int year = 0;
100 int mon, day, hr, minute;//, nexp, ibexp;
101
102 double[] temp = SGP4unit.getgravconst(whichconst);
103 tumin = temp[0];
104 // mu = temp[1];
105 // radiusearthkm = temp[2];
106 // xke = temp[3];
107 // j2 = temp[4];
108 // j3 = temp[5];
109 // j4 = temp[6];
110 // j3oj2 = temp[7];
111
112 satrec.error = 0;
113
114 satrec.name = satName;
115
116 // SEG -- save gravity - moved to SGP4unit.sgp4init fo consistancy
117 //satrec.gravconsttype = whichconst;
118
119 // get variables from the two lines
120 satrec.line1 = line1;
121 try
122 {
123 readLine1(line1, satrec);
124 }
125 catch(Exception e)
126 {
127 System.out.println("Error Reading TLE line 1: " + e.toString());
128 satrec.tleDataOk = false;
129 satrec.error = 7;
130 return false;
131 }
132
133 satrec.line2 = line2;
134 try
135 {
136 readLine2(line2, satrec);
137 }
138 catch(Exception e)
139 {
140 System.out.println("Error Reading TLE line 2: " + e.toString());
141 satrec.tleDataOk = false;
142 satrec.error = 7;
143 return false;
144 }
145
146
147 // ---- find no, ndot, nddot ----
148 satrec.no = satrec.no / xpdotp; //* rad/min
149 satrec.nddot = satrec.nddot * Math.pow(10.0, satrec.nexp);
150 satrec.bstar = satrec.bstar * Math.pow(10.0, satrec.ibexp);
151
152 // ---- convert to sgp4 units ----
153 satrec.a = Math.pow(satrec.no * tumin, (-2.0 / 3.0));
154 satrec.ndot = satrec.ndot / (xpdotp * 1440.0); //* ? * minperday
155 satrec.nddot = satrec.nddot / (xpdotp * 1440.0 * 1440);
156
157 // ---- find standard orbital elements ----
158 satrec.inclo = satrec.inclo * deg2rad;
159 satrec.nodeo = satrec.nodeo * deg2rad;
160 satrec.argpo = satrec.argpo * deg2rad;
161 satrec.mo = satrec.mo * deg2rad;
162
163 satrec.alta = satrec.a * (1.0 + satrec.ecco) - 1.0;
164 satrec.altp = satrec.a * (1.0 - satrec.ecco) - 1.0;
165
166 // ----------------------------------------------------------------
167 // find sgp4epoch time of element set
168 // remember that sgp4 uses units of days from 0 jan 1950 (sgp4epoch)
169 // and minutes from the epoch (time)
170 // ----------------------------------------------------------------
171
172 // ---------------- temp fix for years from 1957-2056 -------------------
173 // --------- correct fix will occur when year is 4-digit in tle ---------
174 if(satrec.epochyr < 57)
175 {
176 year = satrec.epochyr + 2000;
177 }
178 else
179 {
180 year = satrec.epochyr + 1900;
181 }
182
183 // computes the m/d/hr/min/sec from year and epoch days
184 MDHMS mdhms = days2mdhms(year, satrec.epochdays);
185 mon = mdhms.mon;
186 day = mdhms.day;
187 hr = mdhms.hr;
188 minute = mdhms.minute;
189 sec = mdhms.sec;
190 // computes the jd from m/d/...
191 satrec.jdsatepoch = jday(year, mon, day, hr, minute, sec);
192
193 // ---------------- initialize the orbit at sgp4epoch -------------------
194 boolean result = SGP4unit.sgp4init(whichconst, opsmode, satrec.satnum,
195 satrec.jdsatepoch - 2433281.5, satrec.bstar,
196 satrec.ecco, satrec.argpo, satrec.inclo, satrec.mo, satrec.no,
197 satrec.nodeo, satrec);
198
199 return result;
200
201 } // readTLEandIniSGP4
202
203 private static boolean readLine1(String line1, SGP4SatData satrec) throws Exception
204 {
205 String tleLine1 = line1; // first line
206 if(!tleLine1.startsWith("1 "))
207 {
208 throw new Exception("TLE line 1 not valid first line");
209 }
210
211 // satnum
212 satrec.satnum = (int)readFloatFromString(tleLine1.substring(2, 7));
213 // classification
214 satrec.classification = tleLine1.substring(7, 8); // 1 char
215 // intln designator
216 satrec.intldesg = tleLine1.substring(9, 17); // should be within 8
217
218 // epochyr
219 satrec.epochyr = (int)readFloatFromString(tleLine1.substring(18, 20));
220
221 // epoch days
222 satrec.epochdays = readFloatFromString(tleLine1.substring(20, 32));
223
224 // ndot
225 satrec.ndot = readFloatFromString(tleLine1.substring(33, 43));
226
227 // nddot
228 //nexp
229 if((tleLine1.substring(44, 52)).equals(" "))
230 {
231 satrec.nddot = 0;
232 satrec.nexp = 0;
233 }
234 else
235 {
236 satrec.nddot = readFloatFromString(tleLine1.substring(44, 50)) / 1.0E5;
237 //nexp
238 satrec.nexp = (int)readFloatFromString(tleLine1.substring(50, 52));
239 }
240 //bstar
241 satrec.bstar = readFloatFromString(tleLine1.substring(53, 59)) / 1.0E5;
242 //ibex
243 satrec.ibexp = (int)readFloatFromString(tleLine1.substring(59, 61));
244
245 // these last things are not essential so just try to read them, and give a warning - but no error
246 try
247 {
248 // num b.
249 satrec.numb = (int)readFloatFromString(tleLine1.substring(62, 63));
250
251 // elnum
252 satrec.elnum = (long)readFloatFromString(tleLine1.substring(64, 68));
253 }
254 catch(Exception e)
255 {
256 System.out.println("Warning: Error Reading numb or elnum from TLE line 1 sat#:" + satrec.satnum);
257 }
258
259 // checksum
260 //int checksum1 = (int) readFloatFromString(tleLine1.substring(68));
261
262 // if no errors yet everything went ok
263 return true;
264 } // readLine1
265
266 private static boolean readLine2(String line2, SGP4SatData satrec) throws Exception
267 {
268 /* Read the second line of elements. */
269
270 //theLine = aFile.readLine();
271 String tleLine2 = line2; // second line
272 if(!tleLine2.startsWith("2 "))
273 {
274 throw new Exception("TLE line 2 not valid second line");
275 }
276
277 // satnum
278 int satnum = (int)readFloatFromString(tleLine2.substring(2, 7));
279 if(satnum != satrec.satnum)
280 {
281 System.out.println("Warning TLE line 2 Sat Num doesn't match line1 for sat: " + satrec.name);
282 }
283
284 // inclination
285 satrec.inclo = readFloatFromString(tleLine2.substring(8, 17));
286
287 // nodeo
288 satrec.nodeo = readFloatFromString(tleLine2.substring(17, 26));
289
290 //satrec.ecco
291 satrec.ecco = readFloatFromString(tleLine2.substring(26, 34)) / 1.0E7;
292
293 // satrec.argpo
294 satrec.argpo = readFloatFromString(tleLine2.substring(34, 43));
295
296 // satrec.mo
297 satrec.mo = readFloatFromString(tleLine2.substring(43, 52));
298
299 // no
300 satrec.no = readFloatFromString(tleLine2.substring(52, 63));
301
302 // try to read other data
303 try
304 {
305 // revnum
306 satrec.revnum = (long)readFloatFromString(tleLine2.substring(63, 68));
307 }
308 catch(Exception e)
309 {
310 System.out.println("Warning: Error Reading revnum from TLE line 2 sat#:" + satrec.satnum + "\n" + e.toString());
311 satrec.revnum = -1;
312 }
313
314 // // checksum
315 // int checksum2 = (int) readFloatFromString(tleLine2.substring(68));
316
317 return true;
318 } // readLine1
319
320 /**
321 * Read float data from a string
322 * @param inStr
323 * @return
324 * @throws Exception
325 */
326 protected static double readFloatFromString(String inStr) throws Exception
327 {
328 // make sure decimal sparator is '.' so it works in other countries
329 // because of this can't use Double.parse
330 DecimalFormat dformat = new DecimalFormat("#");
331 DecimalFormatSymbols dfs = new DecimalFormatSymbols();
332 dfs.setDecimalSeparator('.');
333 dformat.setDecimalFormatSymbols(dfs);
334
335 // trim white space and if there is a + at the start
336 String trimStr = inStr.trim();
337 if(trimStr.startsWith("+"))
338 {
339 trimStr = trimStr.substring(1);
340 }
341
342 // parse until we hit the end or invalid char
343 ParsePosition pp = new ParsePosition(0);
344 Number num = dformat.parse(trimStr, pp);
345 if(null == num)
346 {
347 throw new Exception("Invalid Float In TLE");
348 }
349
350 return num.doubleValue();
351 } // readFloatFromString
352
353 /** -----------------------------------------------------------------------------
354 *
355 * procedure jday
356 *
357 * this procedure finds the julian date given the year, month, day, and time.
358 * the julian date is defined by each elapsed day since noon, jan 1, 4713 bc.
359 *
360 * algorithm : calculate the answer in one step for efficiency
361 *
362 * author : david vallado 719-573-2600 1 mar 2001
363 *
364 * inputs description range / units
365 * year - year 1900 .. 2100
366 * mon - month 1 .. 12
367 * day - day 1 .. 28,29,30,31
368 * hr - universal time hour 0 .. 23
369 * min - universal time min 0 .. 59
370 * sec - universal time sec 0.0 .. 59.999
371 *
372 * outputs :
373 * jd - julian date days from 4713 bc
374 *
375 * locals :
376 * none.
377 *
378 * coupling :
379 * none.
380 *
381 * references :
382 * vallado 2007, 189, alg 14, ex 3-14
383 *
384 * ---------------------------------------------------------------------------
385 * @param year
386 * @param mon
387 * @param day
388 * @param hr
389 * @param minute
390 * @param sec
391 * @return
392 */
393 public static double jday(
394 int year, int mon, int day, int hr, int minute, double sec//,
395 //double& jd
396 )
397 {
398 double jd;
399 jd = 367.0 * year -
400 Math.floor((7 * (year + Math.floor((mon + 9) / 12.0))) * 0.25) +
401 Math.floor(275 * mon / 9.0) +
402 day + 1721013.5 +
403 ((sec / 60.0 + minute) / 60.0 + hr) / 24.0; // ut in days
404 // - 0.5*sgn(100.0*year + mon - 190002.5) + 0.5;
405
406 return jd;
407 } // end jday
408
409 /* -----------------------------------------------------------------------------
410 *
411 * procedure days2mdhms
412 *
413 * this procedure converts the day of the year, days, to the equivalent month
414 * day, hour, minute and second.
415 *
416 *
417 *
418 * algorithm : set up array for the number of days per month
419 * find leap year - use 1900 because 2000 is a leap year
420 * loop through a temp value while the value is < the days
421 * perform int conversions to the correct day and month
422 * convert remainder into h m s using type conversions
423 *
424 * author : david vallado 719-573-2600 1 mar 2001
425 *
426 * inputs description range / units
427 * year - year 1900 .. 2100
428 * days - julian day of the year 0.0 .. 366.0
429 *
430 * outputs :
431 * mon - month 1 .. 12
432 * day - day 1 .. 28,29,30,31
433 * hr - hour 0 .. 23
434 * min - minute 0 .. 59
435 * sec - second 0.0 .. 59.999
436 *
437 * locals :
438 * dayofyr - day of year
439 * temp - temporary extended values
440 * inttemp - temporary int value
441 * i - index
442 * lmonth[12] - int array containing the number of days per month
443 *
444 * coupling :
445 * none.
446 * --------------------------------------------------------------------------- */
447 // returns MDHMS object with the mdhms variables
448 public static MDHMS days2mdhms(
449 int year, double days//,
450 //int& mon, int& day, int& hr, int& minute, double& sec
451 )
452 {
453 // return variables
454 //int mon, day, hr, minute, sec
455 MDHMS mdhms = new MDHMS();
456
457 int i, inttemp, dayofyr;
458 double temp;
459 int lmonth[] =
460 {
461 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31
462 };
463
464 dayofyr = (int)Math.floor(days);
465 /* ----------------- find month and day of month ---------------- */
466 if((year % 4) == 0) // doesn't work for dates starting 2100 and beyond
467 {
468 lmonth[1] = 29;
469 }
470
471 i = 1;
472 inttemp = 0;
473 while((dayofyr > inttemp + lmonth[i - 1]) && (i < 12))
474 {
475 inttemp = inttemp + lmonth[i - 1];
476 i++;
477 }
478 mdhms.mon = i;
479 mdhms.day = dayofyr - inttemp;
480
481 /* ----------------- find hours minutes and seconds ------------- */
482 temp = (days - dayofyr) * 24.0;
483 mdhms.hr = (int)Math.floor(temp);
484 temp = (temp - mdhms.hr) * 60.0;
485 mdhms.minute = (int)Math.floor(temp);
486 mdhms.sec = (temp - mdhms.minute) * 60.0;
487
488 return mdhms;
489 } // end days2mdhms
490
491 // Month Day Hours Min Sec
492 public static class MDHMS
493 {
494 public int mon = 0;
495 ;
496 public int day = 0;
497 ;
498 public int hr = 0;
499 ;
500 public int minute = 0;
501 ;
502 public double sec = 0;
503 }
504
505 /** -----------------------------------------------------------------------------
506 *
507 * procedure invjday
508 *
509 * this procedure finds the year, month, day, hour, minute and second
510 * given the julian date. tu can be ut1, tdt, tdb, etc.
511 *
512 * algorithm : set up starting values
513 * find leap year - use 1900 because 2000 is a leap year
514 * find the elapsed days through the year in a loop
515 * call routine to find each individual value
516 *
517 * author : david vallado 719-573-2600 1 mar 2001
518 *
519 * inputs description range / units
520 * jd - julian date days from 4713 bc
521 *
522 * outputs :
523 * year - year 1900 .. 2100
524 * mon - month 1 .. 12
525 * day - day 1 .. 28,29,30,31
526 * hr - hour 0 .. 23
527 * min - minute 0 .. 59
528 * sec - second 0.0 .. 59.999
529 *
530 * locals :
531 * days - day of year plus fractional
532 * portion of a day days
533 * tu - julian centuries from 0 h
534 * jan 0, 1900
535 * temp - temporary double values
536 * leapyrs - number of leap years from 1900
537 *
538 * coupling :
539 * days2mdhms - finds month, day, hour, minute and second given days and year
540 *
541 * references :
542 * vallado 2007, 208, alg 22, ex 3-13
543 * ---------------------------------------------------------------------------
544 *
545 * @param jd
546 * @return [year,mon,day,hr,minute,sec]
547 */
548 public static double[] invjday(double jd)
549 {
550 // return vars
551 double year, mon, day, hr, minute, sec;
552
553 int leapyrs;
554 double days, tu, temp;
555
556 /* --------------- find year and days of the year --------------- */
557 temp = jd - 2415019.5;
558 tu = temp / 365.25;
559 year = 1900 + (int)Math.floor(tu);
560 leapyrs = (int)Math.floor((year - 1901) * 0.25);
561
562 // optional nudge by 8.64x10-7 sec to get even outputs
563 days = temp - ((year - 1900) * 365.0 + leapyrs) + 0.00000000001;
564
565 /* ------------ check for case of beginning of a year ----------- */
566 if(days < 1.0)
567 {
568 year = year - 1;
569 leapyrs = (int)Math.floor((year - 1901) * 0.25);
570 days = temp - ((year - 1900) * 365.0 + leapyrs);
571 }
572
573 /* ----------------- find remaing data ------------------------- */
574 //days2mdhms(year, days, mon, day, hr, minute, sec);
575 MDHMS mdhms = days2mdhms((int)year, days);
576 mon = mdhms.mon;
577 day = mdhms.day;
578 hr = mdhms.hr;
579 minute = mdhms.minute;
580 sec = mdhms.sec;
581
582 sec = sec - 0.00000086400; // ?
583
584 return new double[]
585 {
586 year, mon, day, hr, minute, sec
587 };
588 } // end invjday
589
590 /** -----------------------------------------------------------------------------
591 *
592 * function rv2coe
593 *
594 * this function finds the classical orbital elements given the geocentric
595 * equatorial position and velocity vectors.
596 *
597 * author : david vallado 719-573-2600 21 jun 2002
598 *
599 * revisions
600 * vallado - fix special cases 5 sep 2002
601 * vallado - delete extra check in inclination code 16 oct 2002
602 * vallado - add constant file use 29 jun 2003
603 * vallado - add mu 2 apr 2007
604 *
605 * inputs description range / units
606 * r - ijk position vector km
607 * v - ijk velocity vector km / s
608 * mu - gravitational parameter km3 / s2
609 *
610 * outputs :
611 * p - semilatus rectum km
612 * a - semimajor axis km
613 * ecc - eccentricity
614 * incl - inclination 0.0 to pi rad
615 * omega - longitude of ascending node 0.0 to 2pi rad
616 * argp - argument of perigee 0.0 to 2pi rad
617 * nu - true anomaly 0.0 to 2pi rad
618 * m - mean anomaly 0.0 to 2pi rad
619 * arglat - argument of latitude (ci) 0.0 to 2pi rad
620 * truelon - true longitude (ce) 0.0 to 2pi rad
621 * lonper - longitude of periapsis (ee) 0.0 to 2pi rad
622 *
623 * locals :
624 * hbar - angular momentum h vector km2 / s
625 * ebar - eccentricity e vector
626 * nbar - line of nodes n vector
627 * c1 - v**2 - u/r
628 * rdotv - r dot v
629 * hk - hk unit vector
630 * sme - specfic mechanical energy km2 / s2
631 * i - index
632 * e - eccentric, parabolic,
633 * hyperbolic anomaly rad
634 * temp - temporary variable
635 * typeorbit - type of orbit ee, ei, ce, ci
636 *
637 * coupling :
638 * mag - magnitude of a vector
639 * cross - cross product of two vectors
640 * angle - find the angle between two vectors
641 * newtonnu - find the mean anomaly
642 *
643 * references :
644 * vallado 2007, 126, alg 9, ex 2-5
645 * ---------------------------------------------------------------------------
646 *
647 * @param r
648 * @param v
649 * @param mu
650 * @return [p, a, ecc, incl, omega, argp, nu, m, arglat, truelon, lonper]
651 */
652 public static double[] rv2coe(
653 double[] r, double[] v, double mu//,
654 // double& p, double& a, double& ecc, double& incl, double& omega, double& argp,
655 // double& nu, double& m, double& arglat, double& truelon, double& lonper
656 )
657 {
658
659 // return variables
660 double p, a, ecc, incl, omega, argp, nu, m, arglat, truelon, lonper;
661
662 // internal
663
664 double undefined, small, magr, magv, magn, sme,
665 rdotv, infinite, temp, c1, hk, twopi, magh, halfpi, e;
666 double[] hbar = new double[3];
667 double[] nbar = new double[3];
668 double[] ebar = new double[3];
669
670 int i;
671 String typeorbit;
672
673 twopi = 2.0 * Math.PI;
674 halfpi = 0.5 * Math.PI;
675 small = 0.00000001;
676 undefined = 999999.1;
677 infinite = 999999.9;
678
679 // SEG m needs to be ini
680 m = undefined;
681
682 // ------------------------- implementation -----------------
683 magr = mag(r);
684 magv = mag(v);
685
686 // ------------------ find h n and e vectors ----------------
687 cross(r, v, hbar);
688 magh = mag(hbar);
689 if(magh > small)
690 {
691 nbar[0] = -hbar[1];
692 nbar[1] = hbar[0];
693 nbar[2] = 0.0;
694 magn = mag(nbar);
695 c1 = magv * magv - mu / magr;
696 rdotv = dot(r, v);
697 for(i = 0; i <= 2; i++)
698 {
699 ebar[i] = (c1 * r[i] - rdotv * v[i]) / mu;
700 }
701 ecc = mag(ebar);
702
703 // ------------ find a e and semi-latus rectum ----------
704 sme = (magv * magv * 0.5) - (mu / magr);
705 if(Math.abs(sme) > small)
706 {
707 a = -mu / (2.0 * sme);
708 }
709 else
710 {
711 a = infinite;
712 }
713 p = magh * magh / mu;
714
715 // ----------------- find inclination -------------------
716 hk = hbar[2] / magh;
717 incl = Math.acos(hk);
718
719 // -------- determine type of orbit for later use --------
720 // ------ elliptical, parabolic, hyperbolic inclined -------
721 typeorbit = "ei";
722 if(ecc < small)
723 {
724 // ---------------- circular equatorial ---------------
725 if((incl < small) | (Math.abs(incl - Math.PI) < small))
726 {
727 typeorbit = "ce";
728 }
729 else
730 // -------------- circular inclined ---------------
731 {
732 typeorbit = "ci";
733 }
734 }
735 else
736 {
737 // - elliptical, parabolic, hyperbolic equatorial --
738 if((incl < small) | (Math.abs(incl - Math.PI) < small))
739 {
740 typeorbit = "ee";
741 }
742 }
743
744 // ---------- find longitude of ascending node ------------
745 if(magn > small)
746 {
747 temp = nbar[0] / magn;
748 if(Math.abs(temp) > 1.0)
749 {
750 temp = Math.signum(temp);
751 }
752 omega = Math.acos(temp);
753 if(nbar[1] < 0.0)
754 {
755 omega = twopi - omega;
756 }
757 }
758 else
759 {
760 omega = undefined;
761 }
762
763 // ---------------- find argument of perigee ---------------
764 if(typeorbit.equalsIgnoreCase("ei") == true) // 0 = true for cpp strcmp
765 {
766 argp = angle(nbar, ebar);
767 if(ebar[2] < 0.0)
768 {
769 argp = twopi - argp;
770 }
771 }
772 else
773 {
774 argp = undefined;
775 }
776
777 // ------------ find true anomaly at epoch -------------
778 if(typeorbit.startsWith("e"))//typeorbit[0] == 'e' )
779 {
780 nu = angle(ebar, r);
781 if(rdotv < 0.0)
782 {
783 nu = twopi - nu;
784 }
785 }
786 else
787 {
788 nu = undefined;
789 }
790
791 // ---- find argument of latitude - circular inclined -----
792 if(typeorbit.equalsIgnoreCase("ci") == true)
793 {
794 arglat = angle(nbar, r);
795 if(r[2] < 0.0)
796 {
797 arglat = twopi - arglat;
798 }
799 m = arglat;
800 }
801 else
802 {
803 arglat = undefined;
804 }
805
806 // -- find longitude of perigee - elliptical equatorial ----
807 if((ecc > small) && (typeorbit.equalsIgnoreCase("ee") == true))
808 {
809 temp = ebar[0] / ecc;
810 if(Math.abs(temp) > 1.0)
811 {
812 temp = Math.signum(temp);
813 }
814 lonper = Math.acos(temp);
815 if(ebar[1] < 0.0)
816 {
817 lonper = twopi - lonper;
818 }
819 if(incl > halfpi)
820 {
821 lonper = twopi - lonper;
822 }
823 }
824 else
825 {
826 lonper = undefined;
827 }
828
829 // -------- find true longitude - circular equatorial ------
830 if((magr > small) && (typeorbit.equalsIgnoreCase("ce") == true))
831 {
832 temp = r[0] / magr;
833 if(Math.abs(temp) > 1.0)
834 {
835 temp = Math.signum(temp);
836 }
837 truelon = Math.acos(temp);
838 if(r[1] < 0.0)
839 {
840 truelon = twopi - truelon;
841 }
842 if(incl > halfpi)
843 {
844 truelon = twopi - truelon;
845 }
846 m = truelon;
847 }
848 else
849 {
850 truelon = undefined;
851 }
852
853 // ------------ find mean anomaly for all orbits -----------
854 if(typeorbit.startsWith("e"))//typeorbit[0] == 'e' )
855 {
856 double[] tt = newtonnu(ecc, nu);
857 e = tt[0];
858 m = tt[1];
859 }
860 }
861 else
862 {
863 p = undefined;
864 a = undefined;
865 ecc = undefined;
866 incl = undefined;
867 omega = undefined;
868 argp = undefined;
869 nu = undefined;
870 m = undefined;
871 arglat = undefined;
872 truelon = undefined;
873 lonper = undefined;
874 }
875
876 return new double[]
877 {
878 p, a, ecc, incl, omega, argp, nu, m, arglat, truelon, lonper
879 };
880 } // end rv2coe
881
882 /** -----------------------------------------------------------------------------
883 *
884 * function newtonnu
885 *
886 * this function solves keplers equation when the true anomaly is known.
887 * the mean and eccentric, parabolic, or hyperbolic anomaly is also found.
888 * the parabolic limit at 168 is arbitrary. the hyperbolic anomaly is also
889 * limited. the hyperbolic sine is used because it's not double valued.
890 *
891 * author : david vallado 719-573-2600 27 may 2002
892 *
893 * revisions
894 * vallado - fix small 24 sep 2002
895 *
896 * inputs description range / units
897 * ecc - eccentricity 0.0 to
898 * nu - true anomaly -2pi to 2pi rad
899 *
900 * outputs :
901 * e0 - eccentric anomaly 0.0 to 2pi rad 153.02
902 * m - mean anomaly 0.0 to 2pi rad 151.7425
903 *
904 * locals :
905 * e1 - eccentric anomaly, next value rad
906 * sine - sine of e
907 * cose - cosine of e
908 * ktr - index
909 *
910 * coupling :
911 * asinh - arc hyperbolic sine
912 *
913 * references :
914 * vallado 2007, 85, alg 5
915 * ---------------------------------------------------------------------------
916 *
917 * @param ecc
918 * @param nu
919 * @return [e0, m]
920 */
921 public static double[] newtonnu(double ecc, double nu)
922 {
923 // return vars
924 double e0, m;
925
926 // internal
927
928 double small, sine, cose;
929
930 // --------------------- implementation ---------------------
931 e0 = 999999.9;
932 m = 999999.9;
933 small = 0.00000001;
934
935 // --------------------------- circular ------------------------
936 if(Math.abs(ecc) < small)
937 {
938 m = nu;
939 e0 = nu;
940 }
941 else // ---------------------- elliptical -----------------------
942 if(ecc < 1.0 - small)
943 {
944 sine = (Math.sqrt(1.0 - ecc * ecc) * Math.sin(nu)) / (1.0 + ecc * Math.cos(nu));
945 cose = (ecc + Math.cos(nu)) / (1.0 + ecc * Math.cos(nu));
946 e0 = Math.atan2(sine, cose);
947 m = e0 - ecc * Math.sin(e0);
948 }
949 else // -------------------- hyperbolic --------------------
950 if(ecc > 1.0 + small)
951 {
952 if((ecc > 1.0) && (Math.abs(nu) + 0.00001 < Math.PI - Math.acos(1.0 / ecc)))
953 {
954 sine = (Math.sqrt(ecc * ecc - 1.0) * Math.sin(nu)) / (1.0 + ecc * Math.cos(nu));
955 e0 = asinh(sine);
956 m = ecc * Math.sinh(e0) - e0;
957 }
958 }
959 else // ----------------- parabolic ---------------------
960 if(Math.abs(nu) < 168.0 * Math.PI / 180.0)
961 {
962 e0 = Math.tan(nu * 0.5);
963 m = e0 + (e0 * e0 * e0) / 3.0;
964 }
965
966 if(ecc < 1.0)
967 {
968 m = (m % (2.0 * Math.PI));
969 if(m < 0.0)
970 {
971 m = m + 2.0 * Math.PI;
972 }
973 e0 = e0 % (2.0 * Math.PI);
974 }
975
976 return new double[]
977 {
978 e0, m
979 };
980 } // end newtonnu
981
982
983 /* -----------------------------------------------------------------------------
984 *
985 * function asinh
986 *
987 * this function evaluates the inverse hyperbolic sine function.
988 *
989 * author : david vallado 719-573-2600 1 mar 2001
990 *
991 * inputs description range / units
992 * xval - angle value any real
993 *
994 * outputs :
995 * arcsinh - result any real
996 *
997 * locals :
998 * none.
999 *
1000 * coupling :
1001 * none.
1002 *
1003 * --------------------------------------------------------------------------- */
1004 public static double asinh(double xval)
1005 {
1006 return Math.log(xval + Math.sqrt(xval * xval + 1.0));
1007 } // end asinh
1008
1009 /* -----------------------------------------------------------------------------
1010 *
1011 * function mag
1012 *
1013 * this procedure finds the magnitude of a vector. the tolerance is set to
1014 * 0.000001, thus the 1.0e-12 for the squared test of underflows.
1015 *
1016 * author : david vallado 719-573-2600 1 mar 2001
1017 *
1018 * inputs description range / units
1019 * vec - vector
1020 *
1021 * outputs :
1022 * vec - answer stored in fourth component
1023 *
1024 * locals :
1025 * none.
1026 *
1027 * coupling :
1028 * none.
1029 * --------------------------------------------------------------------------- */
1030 public static double mag(double[] x)
1031 {
1032 return Math.sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
1033 } // end mag
1034
1035 /* -----------------------------------------------------------------------------
1036 *
1037 * procedure cross
1038 *
1039 * this procedure crosses two vectors.
1040 *
1041 * author : david vallado 719-573-2600 1 mar 2001
1042 *
1043 * inputs description range / units
1044 * vec1 - vector number 1
1045 * vec2 - vector number 2
1046 *
1047 * outputs :
1048 * outvec - vector result of a x b
1049 *
1050 * locals :
1051 * none.
1052 *
1053 * coupling :
1054 * mag magnitude of a vector
1055 ---------------------------------------------------------------------------- */
1056 public static void cross(double[] vec1, double[] vec2, double[] outvec)
1057 {
1058 outvec[0] = vec1[1] * vec2[2] - vec1[2] * vec2[1];
1059 outvec[1] = vec1[2] * vec2[0] - vec1[0] * vec2[2];
1060 outvec[2] = vec1[0] * vec2[1] - vec1[1] * vec2[0];
1061 } // end cross
1062
1063 /* -----------------------------------------------------------------------------
1064 *
1065 * function dot
1066 *
1067 * this function finds the dot product of two vectors.
1068 *
1069 * author : david vallado 719-573-2600 1 mar 2001
1070 *
1071 * inputs description range / units
1072 * vec1 - vector number 1
1073 * vec2 - vector number 2
1074 *
1075 * outputs :
1076 * dot - result
1077 *
1078 * locals :
1079 * none.
1080 *
1081 * coupling :
1082 * none.
1083 *
1084 * --------------------------------------------------------------------------- */
1085 public static double dot(double[] x, double[] y)
1086 {
1087 return (x[0] * y[0] + x[1] * y[1] + x[2] * y[2]);
1088 } // end dot
1089
1090 /* -----------------------------------------------------------------------------
1091 *
1092 * procedure angle
1093 *
1094 * this procedure calculates the angle between two vectors. the output is
1095 * set to 999999.1 to indicate an undefined value. be sure to check for
1096 * this at the output phase.
1097 *
1098 * author : david vallado 719-573-2600 1 mar 2001
1099 *
1100 * inputs description range / units
1101 * vec1 - vector number 1
1102 * vec2 - vector number 2
1103 *
1104 * outputs :
1105 * theta - angle between the two vectors -pi to pi
1106 *
1107 * locals :
1108 * temp - temporary real variable
1109 *
1110 * coupling :
1111 * dot dot product of two vectors
1112 * --------------------------------------------------------------------------- */
1113 public static double angle(double[] vec1, double[] vec2)
1114 {
1115 double small, undefined, magv1, magv2, temp;
1116 small = 0.00000001;
1117 undefined = 999999.1;
1118
1119 magv1 = mag(vec1);
1120 magv2 = mag(vec2);
1121
1122 if(magv1 * magv2 > small * small)
1123 {
1124 temp = dot(vec1, vec2) / (magv1 * magv2);
1125 if(Math.abs(temp) > 1.0)
1126 {
1127 temp = Math.signum(temp) * 1.0;
1128 }
1129 return Math.acos(temp);
1130 }
1131 else
1132 {
1133 return undefined;
1134 }
1135 } // end angle
1136 }