1 /* 2 * Copyright 2013 Daniel Warner <contact@danrw.com> 3 * 4 * Licensed under the Apache License, Version 2.0 (the "License"); 5 * you may not use this file except in compliance with the License. 6 * You may obtain a copy of the License at 7 * 8 * http://www.apache.org/licenses/LICENSE-2.0 9 * 10 * Unless required by applicable law or agreed to in writing, software 11 * distributed under the License is distributed on an "AS IS" BASIS, 12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 13 * See the License for the specific language governing permissions and 14 * limitations under the License. 15 */ 16 17 #include <keplerian_toolbox/third_party/libsgp4/OrbitalElements.h> 18 #include <keplerian_toolbox/third_party/libsgp4/Tle.h> 19 OrbitalElements(const Tle & tle)20OrbitalElements::OrbitalElements(const Tle &tle) 21 { 22 /* 23 * extract and format tle data 24 */ 25 mean_anomoly_ = tle.MeanAnomaly(false); 26 ascending_node_ = tle.RightAscendingNode(false); 27 argument_perigee_ = tle.ArgumentPerigee(false); 28 eccentricity_ = tle.Eccentricity(); 29 inclination_ = tle.Inclination(false); 30 mean_motion_ = tle.MeanMotion() * kTWOPI / kMINUTES_PER_DAY; 31 bstar_ = tle.BStar(); 32 epoch_ = tle.Epoch(); 33 34 /* 35 * recover original mean motion (xnodp) and semimajor axis (aodp) 36 * from input elements 37 */ 38 const double a1 = pow(kXKE / MeanMotion(), kTWOTHIRD); 39 const double cosio = cos(Inclination()); 40 const double theta2 = cosio * cosio; 41 const double x3thm1 = 3.0 * theta2 - 1.0; 42 const double eosq = Eccentricity() * Eccentricity(); 43 const double betao2 = 1.0 - eosq; 44 const double betao = sqrt(betao2); 45 const double temp = (1.5 * kCK2) * x3thm1 / (betao * betao2); 46 const double del1 = temp / (a1 * a1); 47 const double a0 = a1 * (1.0 - del1 * (1.0 / 3.0 + del1 * (1.0 + del1 * 134.0 / 81.0))); 48 const double del0 = temp / (a0 * a0); 49 50 recovered_mean_motion_ = MeanMotion() / (1.0 + del0); 51 /* 52 * alternative way to calculate 53 * doesnt affect final results 54 * recovered_semi_major_axis_ = pow(XKE / RecoveredMeanMotion(), TWOTHIRD); 55 */ 56 recovered_semi_major_axis_ = a0 / (1.0 - del0); 57 58 /* 59 * find perigee and period 60 */ 61 perigee_ = (RecoveredSemiMajorAxis() * (1.0 - Eccentricity()) - kAE) * kXKMPER; 62 period_ = kTWOPI / RecoveredMeanMotion(); 63 } 64