From 49bb0961268c48d189d6ed45f09787d533613d11 Mon Sep 17 00:00:00 2001 From: Andreas Kristiansen <6977712+AndreasK79@users.noreply.github.com> Date: Sun, 21 Jul 2024 13:32:11 +0200 Subject: [PATCH] [Sat planner] First version of sat planner added --- application/controllers/Satellite.php | 99 ++ application/views/interface_assets/header.php | 4 +- application/views/satellite/pass.php | 649 ++++++++++ application/views/satellite/passtable.php | 47 + assets/js/sections/satpasses.js | 23 + predict/Predict.php | 877 ++++++++++++++ predict/Predict/DeepArg.php | 35 + predict/Predict/DeepStatic.php | 83 ++ predict/Predict/Exception.php | 5 + predict/Predict/Geodetic.php | 13 + predict/Predict/Math.php | 196 +++ predict/Predict/ObsSet.php | 12 + predict/Predict/Pass.php | 27 + predict/Predict/PassDetail.php | 38 + predict/Predict/QTH.php | 20 + predict/Predict/SGPObs.php | 156 +++ predict/Predict/SGPSDP.php | 1061 +++++++++++++++++ predict/Predict/SGSDPStatic.php | 39 + predict/Predict/Sat.php | 324 +++++ predict/Predict/Solar.php | 115 ++ predict/Predict/TLE.php | 232 ++++ predict/Predict/Time.php | 222 ++++ predict/Predict/Vector.php | 13 + predict/README.md | 35 + predict/composer.json | 19 + predict/generatePackage.php | 63 + predict/package.xml | 182 +++ 27 files changed, 4588 insertions(+), 1 deletion(-) create mode 100644 application/views/satellite/pass.php create mode 100644 application/views/satellite/passtable.php create mode 100644 assets/js/sections/satpasses.js create mode 100644 predict/Predict.php create mode 100644 predict/Predict/DeepArg.php create mode 100644 predict/Predict/DeepStatic.php create mode 100644 predict/Predict/Exception.php create mode 100644 predict/Predict/Geodetic.php create mode 100644 predict/Predict/Math.php create mode 100644 predict/Predict/ObsSet.php create mode 100644 predict/Predict/Pass.php create mode 100644 predict/Predict/PassDetail.php create mode 100644 predict/Predict/QTH.php create mode 100644 predict/Predict/SGPObs.php create mode 100644 predict/Predict/SGPSDP.php create mode 100644 predict/Predict/SGSDPStatic.php create mode 100644 predict/Predict/Sat.php create mode 100644 predict/Predict/Solar.php create mode 100644 predict/Predict/TLE.php create mode 100644 predict/Predict/Time.php create mode 100644 predict/Predict/Vector.php create mode 100644 predict/README.md create mode 100644 predict/composer.json create mode 100644 predict/generatePackage.php create mode 100644 predict/package.xml diff --git a/application/controllers/Satellite.php b/application/controllers/Satellite.php index 9679c1ce8..7ebfe720f 100644 --- a/application/controllers/Satellite.php +++ b/application/controllers/Satellite.php @@ -173,4 +173,103 @@ class Satellite extends CI_Controller { echo json_encode($satellite_data, JSON_FORCE_OBJECT); } + public function pass() { + $this->load->model('satellite_model'); + $this->load->model('stations'); + $active_station_id = $this->stations->find_active(); + $pageData['activegrid'] = $this->stations->gridsquare_from_station($active_station_id); + + $pageData['satellites'] = $this->satellite_model->get_all_satellites_with_tle(); + + $footerData = []; + $footerData['scripts'] = [ + 'assets/js/sections/satpasses.js?' . filemtime(realpath(__DIR__ . "/../../assets/js/sections/satpasses.js")), + ]; + + // Render Page + $pageData['page_title'] = "Satellite pass"; + $this->load->view('interface_assets/header', $pageData); + $this->load->view('satellite/pass'); + $this->load->view('interface_assets/footer', $footerData); + } + + public function searchpasses() { + try { + $result = $this->get_tle_for_predict(); + $this->calcpass($result); + } + catch (Exception $e) { + header("Content-type: application/json"); + echo json_encode(['ok' => 'Error', 'message' => $e->getMessage() . $e->getCode()]); + } + } + + public function get_tle_for_predict() { + $sat = $this->security->xss_clean($this->input->post('sat')); + $this->load->model('satellite_model'); + return $this->satellite_model->get_tle($sat); + } + + function calcpass($sat_tle) { + require_once realpath(__DIR__ . "/../../predict/Predict.php"); + require_once realpath(__DIR__ . "/../../predict/Predict/Sat.php"); + require_once realpath(__DIR__ . "/../../predict/Predict/QTH.php"); + require_once realpath(__DIR__ . "/../../predict/Predict/Time.php"); + require_once realpath(__DIR__ . "/../../predict/Predict/TLE.php"); + + // Track execution time of this script + $start = microtime(true); + + // The observer or groundstation is called QTH in ham radio terms + $predict = new Predict(); + $qth = new Predict_QTH(); + $qth->alt = $this->security->xss_clean($this->input->post('altitude')); // Altitude in meters + + $strQRA = $this->security->xss_clean($this->input->post('yourgrid')); + + if ((strlen($strQRA) % 2 == 0) && (strlen($strQRA) <= 10)) { // Check if QRA is EVEN (the % 2 does that) and smaller/equal 8 + $strQRA = strtoupper($strQRA); + if (strlen($strQRA) == 4) $strQRA .= "LL"; // Only 4 Chars? Fill with center "LL" as only A-R allowed + if (strlen($strQRA) == 6) $strQRA .= "55"; // Only 6 Chars? Fill with center "55" + if (strlen($strQRA) == 8) $strQRA .= "LL"; // Only 8 Chars? Fill with center "LL" as only A-R allowed + + if (!preg_match('/^[A-R]{2}[0-9]{2}[A-X]{2}[0-9]{2}[A-X]{2}$/', $strQRA)) { + return false; + } + } + + if(!$this->load->is_loaded('Qra')) { + $this->load->library('Qra'); + } + $homecoordinates = $this->qra->qra2latlong($this->security->xss_clean($this->input->post('yourgrid'))); + + $qth->lat = $homecoordinates[0]; + $qth->lon = $homecoordinates[1]; + + $temp = preg_split('/\n/', $sat_tle->tle); + + $tle = new Predict_TLE($sat_tle->satellite, $temp[0], $temp[1]); // Instantiate it + $sat = new Predict_Sat($tle); // Load up the satellite data + + $now = Predict_Time::get_current_daynum(); // get the current time as Julian Date (daynum) + + // You can modify some preferences in Predict(), the defaults are below + // + $predict->minEle = $this->security->xss_clean($this->input->post('minelevation')); // Minimum elevation for a pass + // $predict->timeRes = 10; // Pass details: time resolution in seconds + // $predict->numEntries = 20; // Pass details: number of entries per pass + // $predict->threshold = -6; // Twilight threshold (sun must be at this lat or lower) + + // Get the passes and filter visible only, takes about 4 seconds for 10 days + $results = $predict->get_passes($sat, $qth, $now, 10); + $filtered = $predict->filterVisiblePasses($results); + + $zone = $this->security->xss_clean($this->input->post('timezone')); + $format = 'm-d-Y H:i:s'; // Time format from PHP's date() function + + $data['filtered'] = $filtered; + $data['zone'] = $zone; + $data['format'] = $format; + $this->load->view('satellite/passtable', $data); + } } diff --git a/application/views/interface_assets/header.php b/application/views/interface_assets/header.php index faafb4446..c2ee0d12f 100644 --- a/application/views/interface_assets/header.php +++ b/application/views/interface_assets/header.php @@ -271,6 +271,8 @@
+ +
@@ -404,7 +406,7 @@ if (!($this->config->item('disable_oqrs') ?? false)) { $oqrs_requests = $this->oqrs_model->oqrs_requests($location_list); ?> -
  • +
  • 0) { echo "" . $oqrs_requests . ""; } ?>
    • diff --git a/application/views/satellite/pass.php b/application/views/satellite/pass.php new file mode 100644 index 000000000..5990dc924 --- /dev/null +++ b/application/views/satellite/pass.php @@ -0,0 +1,649 @@ +
    +

    +
    +
    +
    +
    + + +
    +
    + + +
    +
    + + +
    +
    + + +
    +
    + + +
    +
    + + +
    +
    + + +
    +
    + + +
    +
    + + +
    +
    + + +
    +
    + +
    +
    + +
    +
    +
    + diff --git a/application/views/satellite/passtable.php b/application/views/satellite/passtable.php new file mode 100644 index 000000000..9850e3a7a --- /dev/null +++ b/application/views/satellite/passtable.php @@ -0,0 +1,47 @@ + + + Satellite + AOS Time + Duration + AOS az + AOS el + Max El + LOS Time + LOS Az + LOS El + '; + foreach ($filtered as $pass) { + echo ''; + echo '' . $pass->satname . ''; + echo '' . Predict_Time::daynum2readable($pass->visible_aos, $zone, $format) . ''; + echo '' . returntimediff(Predict_Time::daynum2readable($pass->visible_aos, $zone, $format), Predict_Time::daynum2readable($pass->visible_los, $zone, $format)) . ''; + echo '' . round($pass->visible_aos_az) . ' (' . azDegreesToDirection($pass->visible_aos_az) . ')'; + echo '' . round($pass->visible_aos_el) . ''; + echo '' . round($pass->visible_max_el) . ''; + echo '' . Predict_Time::daynum2readable($pass->visible_los, $zone, $format) . ''; + echo '' . round($pass->visible_los_az) . ' (' . azDegreesToDirection($pass->visible_los_az) . ')'; + echo '' . round($pass->visible_los_el) . ''; + echo ''; + } + echo ''; +} + +function returntimediff($start, $end) { + $datetime1 = DateTime::createFromFormat('m-d-Y H:i:s', $end); + $datetime2 = DateTime::createFromFormat('m-d-Y H:i:s', $start); + $interval = $datetime1->diff($datetime2); + + $minutesDifference = ($interval->days * 24 * 60) + ($interval->h * 60) + $interval->i + ($interval->s / 60); + + return round($minutesDifference) . ' min'; +} + +function azDegreesToDirection($az = 0) { + $i = floor($az / 22.5); + $m = (22.5 * (2 * $i + 1)) / 2; + $i = ($az >= $m) ? $i + 1 : $i; + + return trim(substr('N NNENE ENEE ESESE SSES SSWSW WSWW WNWNW NNWN ', $i * 3, 3)); +} diff --git a/assets/js/sections/satpasses.js b/assets/js/sections/satpasses.js new file mode 100644 index 000000000..7515befd1 --- /dev/null +++ b/assets/js/sections/satpasses.js @@ -0,0 +1,23 @@ +function searchpasses() { + $.ajax({ + url: base_url + 'index.php/satellite/searchpasses', + type: 'post', + data: {'sat': $("#satlist").val(), + 'yourgrid': $("#yourgrid").val(), + 'minelevation': $("#minelevation").val(), + 'minazimuth': $("#minazimuth").val(), + 'maxazimuth': $("#maxazimuth").val(), + 'altitude': $("#altitude").val(), + 'timezone': $("#timezone").val(), + 'date': $("#date").val(), + 'mintime': $("#mintime").val(), + 'maxtime': $("#maxtime").val(), + }, + success: function (html) { + $("#resultpasses").html(html); + }, + error: function(e) { + modalloading=false; + } + }); +} diff --git a/predict/Predict.php b/predict/Predict.php new file mode 100644 index 000000000..591ea49d3 --- /dev/null +++ b/predict/Predict.php @@ -0,0 +1,877 @@ + + John A. Magliacane, KD2BD. + + Comments, questions and bugreports should be submitted via + http://sourceforge.net/projects/gpredict/ + More details can be found at the project home page: + + http://gpredict.oz9aec.net/ + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, visit http://www.fsf.org/ +*/ + +require_once realpath(__DIR__ . "/../predict/Predict/Time.php"); +require_once realpath(__DIR__ . "/../predict/Predict/Math.php"); +require_once realpath(__DIR__ . "/../predict/Predict/Pass.php"); +require_once realpath(__DIR__ . "/../predict/Predict/PassDetail.php"); +require_once realpath(__DIR__ . "/../predict/Predict/Vector.php"); +require_once realpath(__DIR__ . "/../predict/Predict/Geodetic.php"); +require_once realpath(__DIR__ . "/../predict/Predict/ObsSet.php"); +require_once realpath(__DIR__ . "/../predict/Predict/Solar.php"); +require_once realpath(__DIR__ . "/../predict/Predict/SGPSDP.php"); +require_once realpath(__DIR__ . "/../predict/Predict/SGPSDP.php"); + +/** + * The main Predict class. Contains constants for use by other classes, as well as + * the prediction logic. + */ +class Predict +{ + const de2ra = 1.74532925E-2; /* Degrees to Radians */ + const pi = 3.1415926535898; /* Pi */ + const pio2 = 1.5707963267949; /* Pi/2 */ + const x3pio2 = 4.71238898; /* 3*Pi/2 */ + const twopi = 6.2831853071796; /* 2*Pi */ + const e6a = 1.0E-6; + const tothrd = 6.6666667E-1; /* 2/3 */ + const xj2 = 1.0826158E-3; /* J2 Harmonic */ + const xj3 = -2.53881E-6; /* J3 Harmonic */ + const xj4 = -1.65597E-6; /* J4 Harmonic */ + const xke = 7.43669161E-2; + const xkmper = 6.378135E3; /* Earth radius km */ + const xmnpda = 1.44E3; /* Minutes per day */ + const km2mi = 0.621371; /* Kilometers per Mile */ + const ae = 1.0; + const ck2 = 5.413079E-4; + const ck4 = 6.209887E-7; + const __f = 3.352779E-3; + const ge = 3.986008E5; + const __s__ = 1.012229; + const qoms2t = 1.880279E-09; + const secday = 8.6400E4; /* Seconds per day */ + const omega_E = 1.0027379; + const omega_ER = 6.3003879; + const zns = 1.19459E-5; + const c1ss = 2.9864797E-6; + const zes = 1.675E-2; + const znl = 1.5835218E-4; + const c1l = 4.7968065E-7; + const zel = 5.490E-2; + const zcosis = 9.1744867E-1; + const zsinis = 3.9785416E-1; + const zsings = -9.8088458E-1; + const zcosgs = 1.945905E-1; + const zcoshs = 1; + const zsinhs = 0; + const q22 = 1.7891679E-6; + const q31 = 2.1460748E-6; + const q33 = 2.2123015E-7; + const g22 = 5.7686396; + const g32 = 9.5240898E-1; + const g44 = 1.8014998; + const g52 = 1.0508330; + const g54 = 4.4108898; + const root22 = 1.7891679E-6; + const root32 = 3.7393792E-7; + const root44 = 7.3636953E-9; + const root52 = 1.1428639E-7; + const root54 = 2.1765803E-9; + const thdt = 4.3752691E-3; + const rho = 1.5696615E-1; + const mfactor = 7.292115E-5; + const __sr__ = 6.96000E5; /*Solar radius - kilometers (IAU 76)*/ + const AU = 1.49597870E8; /*Astronomical unit - kilometers (IAU 76)*/ + + /* visibility constants */ + const SAT_VIS_NONE = 0; + const SAT_VIS_VISIBLE = 1; + const SAT_VIS_DAYLIGHT = 2; + const SAT_VIS_ECLIPSED = 3; + + /* preferences */ + public $minEle = 10; // Minimum elevation + public $timeRes = 10; // Pass details: time resolution + public $numEntries = 20; // Pass details: number of entries + public $threshold = -6; // Twilight threshold + + /** + * Predict the next pass. + * + * This function simply wraps the get_pass function using the current time + * as parameter. + * + * Note: the data in sat will be corrupt (future) and must be refreshed + * by the caller, if the caller will need it later on (eg. if the caller + * is GtkSatList). + * + * @param Predict_Sat $sat The satellite data. + * @param Predict_QTH $qth The observer data. + * @param int $maxdt The maximum number of days to look ahead. + * + * @return Predict_Pass Pointer instance or NULL if no pass can be + * found. + */ + public function get_next_pass(Predict_Sat $sat, Predict_QTH $qth, $maxdt) + { + /* get the current time and call the get_pass function */ + $now = Predict_Time::get_current_daynum(); + + return $this->get_pass($sat, $qth, $now, $maxdt); + } + + /** Predict first pass after a certain time. + * + * @param Predict_Sat $sat The satellite data. + * @param Predict_QTH $qth The observer's location data. + * @param float $start Starting time. + * @param int $maxdt The maximum number of days to look ahead (0 for no limit). + * + * @return Predict_Pass or NULL if there was an error. + * + * This function will find the first upcoming pass with AOS no earlier than + * t = start and no later than t = (start+maxdt). + * + * note For no time limit use maxdt = 0.0 + * + * note the data in sat will be corrupt (future) and must be refreshed + * by the caller, if the caller will need it later on + */ + public function get_pass(Predict_Sat $sat_in, Predict_QTH $qth, $start, $maxdt) + { + + $aos = 0.0; /* time of AOS */ + $tca = 0.0; /* time of TCA */ + $los = 0.0; /* time of LOS */ + $dt = 0.0; /* time diff */ + $step = 0.0; /* time step */ + $t0 = $start; + $tres = 0.0; /* required time resolution */ + $max_el = 0.0; /* maximum elevation */ + $pass = null; + $detail = null; + $done = false; + $iter = 0; /* number of iterations */ + /* FIXME: watchdog */ + + /*copy sat_in to a working structure*/ + $sat = clone $sat_in; + $sat_working = clone $sat_in; + + /* get time resolution; sat-cfg stores it in seconds */ + $tres = $this->timeRes / 86400.0; + + /* loop until we find a pass with elevation > SAT_CFG_INT_PRED_MIN_EL + or we run out of time + FIXME: we should have a safety break + */ + while (!$done) { + /* Find los of next pass or of current pass */ + $los = $this->find_los($sat, $qth, $t0, $maxdt); // See if a pass is ongoing + $aos = $this->find_aos($sat, $qth, $t0, $maxdt); + /* sat_log_log(SAT_LOG_LEVEL_MSG, "%s:%s:%d: found aos %f and los %f for t0=%f", */ + /* __FILE__, */ + /* __FUNCTION__, */ + /* __LINE__, */ + /* aos, */ + /* los, */ + /* t0); */ + if ($aos > $los) { + // los is from an currently happening pass, find previous aos + $aos = $this->find_prev_aos($sat, $qth, $t0); + } + + /* aos = 0.0 means no aos */ + if ($aos == 0.0) { + $done = true; + } else if (($maxdt > 0.0) && ($aos > ($start + $maxdt)) ) { + /* check whether we are within time limits; + maxdt = 0 mean no time limit. + */ + $done = true; + } else { + //los = find_los (sat, qth, aos + 0.001, maxdt); // +1.5 min later + $dt = $los - $aos; + + /* get time step, which will give us the max number of entries */ + $step = $dt / $this->numEntries; + + /* but if this is smaller than the required resolution + we go with the resolution + */ + if ($step < $tres) { + $step = $tres; + } + + /* create a pass_t entry; FIXME: g_try_new in 2.8 */ + $pass = new Predict_Pass(); + + $pass->aos = $aos; + $pass->los = $los; + $pass->max_el = 0.0; + $pass->aos_az = 0.0; + $pass->los_az = 0.0; + $pass->maxel_az = 0.0; + $pass->vis = '---'; + $pass->satname = $sat->nickname; + $pass->details = array(); + + /* iterate over each time step */ + for ($t = $pass->aos; $t <= $pass->los; $t += $step) { + + /* calculate satellite data */ + $this->predict_calc($sat, $qth, $t); + + /* in the first iter we want to store + pass->aos_az + */ + if ($t == $pass->aos) { + $pass->aos_az = $sat->az; + $pass->orbit = $sat->orbit; + } + + /* append details to sat->details */ + $detail = new Predict_PassDetail(); + $detail->time = $t; + $detail->pos->x = $sat->pos->x; + $detail->pos->y = $sat->pos->y; + $detail->pos->z = $sat->pos->z; + $detail->pos->w = $sat->pos->w; + $detail->vel->x = $sat->vel->x; + $detail->vel->y = $sat->vel->y; + $detail->vel->z = $sat->vel->z; + $detail->vel->w = $sat->vel->w; + $detail->velo = $sat->velo; + $detail->az = $sat->az; + $detail->el = $sat->el; + $detail->range = $sat->range; + $detail->range_rate = $sat->range_rate; + $detail->lat = $sat->ssplat; + $detail->lon = $sat->ssplon; + $detail->alt = $sat->alt; + $detail->ma = $sat->ma; + $detail->phase = $sat->phase; + $detail->footprint = $sat->footprint; + $detail->orbit = $sat->orbit; + $detail->vis = $this->get_sat_vis($sat, $qth, $t); + + /* also store visibility "bit" */ + switch ($detail->vis) { + case self::SAT_VIS_VISIBLE: + $pass->vis[0] = 'V'; + break; + case self::SAT_VIS_DAYLIGHT: + $pass->vis[1] = 'D'; + break; + case self::SAT_VIS_ECLIPSED: + $pass->vis[2] = 'E'; + break; + default: + break; + } + + // Using an array, no need to prepend and reverse the list + // as gpredict does + $pass->details[] = $detail; + + // Look up apparent magnitude if this is a visible pass + if ($detail->vis === self::SAT_VIS_VISIBLE) { + $apmag = $sat->calculateApparentMagnitude($t, $qth); + if ($pass->max_apparent_magnitude === null || $apmag < $pass->max_apparent_magnitude) { + $pass->max_apparent_magnitude = $apmag; + } + } + + /* store elevation if greater than the + previously stored one + */ + if ($sat->el > $max_el) { + $max_el = $sat->el; + $tca = $t; + $pass->maxel_az = $sat->az; + } + + /* g_print ("TIME: %f\tAZ: %f\tEL: %f (MAX: %f)\n", */ + /* t, sat->az, sat->el, max_el); */ + } + + /* calculate satellite data */ + $this->predict_calc($sat, $qth, $pass->los); + /* store los_az, max_el and tca */ + $pass->los_az = $sat->az; + $pass->max_el = $max_el; + $pass->tca = $tca; + + /* check whether this pass is good */ + if ($max_el >= $this->minEle) { + $done = true; + } else { + $done = false; + $t0 = $los + 0.014; // +20 min + $pass = null; + } + + $iter++; + } + } + + return $pass; + } + + /** + * Calculate satellite visibility. + * + * @param Predict_Sat $sat The satellite structure. + * @param Predict_QTH $qth The QTH + * @param float $jul_utc The time at which the visibility should be calculated. + * + * @return int The visiblity constant, 0, 1, 2, or 3 (see above) + */ + public function get_sat_vis(Predict_Sat $sat, Predict_QTH $qth, $jul_utc) + { + /* gboolean sat_sun_status; + gdouble sun_el; + gdouble threshold; + gdouble eclipse_depth; + sat_vis_t vis = SAT_VIS_NONE; */ + + $eclipse_depth = 0.0; + $zero_vector = new Predict_Vector(); + $obs_geodetic = new Predict_Geodetic(); + + /* Solar ECI position vector */ + $solar_vector = new Predict_Vector(); + + /* Solar observed az and el vector */ + $solar_set = new Predict_ObsSet(); + + /* FIXME: could be passed as parameter */ + $obs_geodetic->lon = $qth->lon * self::de2ra; + $obs_geodetic->lat = $qth->lat * self::de2ra; + $obs_geodetic->alt = $qth->alt / 1000.0; + $obs_geodetic->theta = 0; + + Predict_Solar::Calculate_Solar_Position($jul_utc, $solar_vector); + Predict_SGPObs::Calculate_Obs($jul_utc, $solar_vector, $zero_vector, $obs_geodetic, $solar_set); + + if (Predict_Solar::Sat_Eclipsed($sat->pos, $solar_vector, $eclipse_depth)) { + /* satellite is eclipsed */ + $sat_sun_status = false; + } else { + /* satellite in sunlight => may be visible */ + $sat_sun_status = true; + } + + if ($sat_sun_status) { + $sun_el = Predict_Math::Degrees($solar_set->el); + + if ($sun_el <= $this->threshold && $sat->el >= 0.0) { + $vis = self::SAT_VIS_VISIBLE; + } else { + $vis = self::SAT_VIS_DAYLIGHT; + } + } else { + $vis = self::SAT_VIS_ECLIPSED; + } + + return $vis; + } + + /** Find the AOS time of the next pass. + * @author Alexandru Csete, OZ9AEC + * @author John A. Magliacane, KD2BD + * @param Predict_Sat $sat The satellite data. + * @param Predict_QTH $qth The observer's location (QTH) data. + * @param float $start The julian date where calculation should start. + * @param int $maxdt The upper time limit in days (0.0 = no limit) + * @return The julain date of the next AOS or 0.0 if the satellite has no AOS. + * + * This function finds the time of AOS for the first coming pass taking place + * no earlier that start. + * If the satellite is currently within range, the function first calls + * find_los to get the next LOS time. Then the calculations are done using + * the new start time. + * + */ + public function find_aos(Predict_Sat $sat, Predict_QTH $qth, $start, $maxdt) + { + $t = $start; + $aostime = 0.0; + + + /* make sure current sat values are + in sync with the time + */ + $this->predict_calc($sat, $qth, $start); + + /* check whether satellite has aos */ + if (($sat->otype == Predict_SGPSDP::ORBIT_TYPE_GEO) || + ($sat->otype == Predict_SGPSDP::ORBIT_TYPE_DECAYED) || + !$this->has_aos($sat, $qth)) { + + return 0.0; + } + + if ($sat->el > 0.0) { + $t = $this->find_los($sat, $qth, $start, $maxdt) + 0.014; // +20 min + } + + /* invalid time (potentially returned by find_los) */ + if ($t < 0.1) { + return 0.0; + } + + /* update satellite data */ + $this->predict_calc($sat, $qth, $t); + + /* use upper time limit */ + if ($maxdt > 0.0) { + + /* coarse time steps */ + while (($sat->el < -1.0) && ($t <= ($start + $maxdt))) { + $t -= 0.00035 * ($sat->el * (($sat->alt / 8400.0) + 0.46) - 2.0); + $this->predict_calc($sat, $qth, $t); + } + + /* fine steps */ + while (($aostime == 0.0) && ($t <= ($start + $maxdt))) { + + if (abs($sat->el) < 0.005) { + $aostime = $t; + } else { + $t -= $sat->el * sqrt($sat->alt) / 530000.0; + $this->predict_calc($sat, $qth, $t); + } + } + } else { + /* don't use upper time limit */ + + /* coarse time steps */ + while ($sat->el < -1.0) { + + $t -= 0.00035 * ($sat->el * (($sat->alt / 8400.0) + 0.46) - 2.0); + $this->predict_calc($sat, $qth, $t); + } + + /* fine steps */ + while ($aostime == 0.0) { + + if (abs($sat->el) < 0.005) { + $aostime = $t; + } else { + $t -= $sat->el * sqrt($sat->alt) / 530000.0; + $this->predict_calc($sat, $qth, $t); + } + + } + } + + return $aostime; + } + + /** SGP4SDP4 driver for doing AOS/LOS calculations. + * @param Predict_Sat $sat The satellite data. + * @param Predict_QTH $qth The QTH observer location data. + * @param float $t The time for calculation (Julian Date) + * + */ + public function predict_calc(Predict_Sat $sat, Predict_QTH $qth, $t) + { + $obs_set = new Predict_ObsSet(); + $sat_geodetic = new Predict_Geodetic(); + $obs_geodetic = new Predict_Geodetic(); + + $obs_geodetic->lon = $qth->lon * self::de2ra; + $obs_geodetic->lat = $qth->lat * self::de2ra; + $obs_geodetic->alt = $qth->alt / 1000.0; + $obs_geodetic->theta = 0; + + $sat->jul_utc = $t; + $sat->tsince = ($sat->jul_utc - $sat->jul_epoch) * self::xmnpda; + + /* call the norad routines according to the deep-space flag */ + $sgpsdp = Predict_SGPSDP::getInstance($sat); + if ($sat->flags & Predict_SGPSDP::DEEP_SPACE_EPHEM_FLAG) { + $sgpsdp->SDP4($sat, $sat->tsince); + } else { + $sgpsdp->SGP4($sat, $sat->tsince); + } + + Predict_Math::Convert_Sat_State($sat->pos, $sat->vel); + + /* get the velocity of the satellite */ + $sat->vel->w = sqrt($sat->vel->x * $sat->vel->x + $sat->vel->y * $sat->vel->y + $sat->vel->z * $sat->vel->z); + $sat->velo = $sat->vel->w; + Predict_SGPObs::Calculate_Obs($sat->jul_utc, $sat->pos, $sat->vel, $obs_geodetic, $obs_set); + Predict_SGPObs::Calculate_LatLonAlt($sat->jul_utc, $sat->pos, $sat_geodetic); + + while ($sat_geodetic->lon < -self::pi) { + $sat_geodetic->lon += self::twopi; + } + + while ($sat_geodetic->lon > (self::pi)) { + $sat_geodetic->lon -= self::twopi; + } + + $sat->az = Predict_Math::Degrees($obs_set->az); + $sat->el = Predict_Math::Degrees($obs_set->el); + $sat->range = $obs_set->range; + $sat->range_rate = $obs_set->range_rate; + $sat->ssplat = Predict_Math::Degrees($sat_geodetic->lat); + $sat->ssplon = Predict_Math::Degrees($sat_geodetic->lon); + $sat->alt = $sat_geodetic->alt; + $sat->ma = Predict_Math::Degrees($sat->phase); + $sat->ma *= 256.0 / 360.0; + $sat->phase = Predict_Math::Degrees($sat->phase); + + /* same formulas, but the one from predict is nicer */ + //sat->footprint = 2.0 * xkmper * acos (xkmper/sat->pos.w); + $sat->footprint = 12756.33 * acos(self::xkmper / (self::xkmper + $sat->alt)); + $age = $sat->jul_utc - $sat->jul_epoch; + $sat->orbit = floor(($sat->tle->xno * self::xmnpda / self::twopi + + $age * $sat->tle->bstar * self::ae) * $age + + $sat->tle->xmo / self::twopi) + $sat->tle->revnum - 1; + } + + /** Find the LOS time of the next pass. + * @author Alexandru Csete, OZ9AEC + * @author John A. Magliacane, KD2BD + * @param Predict_Sat $sat The satellite data. + * @param Predict_QTH $qth The QTH observer location data. + * @param float $start The time where calculation should start. (Julian Date) + * @param int $maxdt The upper time limit in days (0.0 = no limit) + * @return The time (julian date) of the next LOS or 0.0 if the satellite has no LOS. + * + * This function finds the time of LOS for the first coming pass taking place + * no earlier that start. + * If the satellite is currently out of range, the function first calls + * find_aos to get the next AOS time. Then the calculations are done using + * the new start time. + * The function has a built-in watchdog to ensure that we don't end up in + * lengthy loops. + * + */ + public function find_los(Predict_Sat $sat, Predict_QTH $qth, $start, $maxdt) + { + $t = $start; + $lostime = 0.0; + + + $this->predict_calc($sat, $qth, $start); + + /* check whether satellite has aos */ + if (($sat->otype == Predict_SGPSDP::ORBIT_TYPE_GEO) || + ($sat->otype == Predict_SGPSDP::ORBIT_TYPE_DECAYED) || + !$this->has_aos ($sat, $qth)) { + + return 0.0; + } + + if ($sat->el < 0.0) { + $t = $this->find_aos($sat, $qth, $start, $maxdt) + 0.001; // +1.5 min + } + + /* invalid time (potentially returned by find_aos) */ + if ($t < 0.01) { + return 0.0; + } + + /* update satellite data */ + $this->predict_calc($sat, $qth, $t); + + /* use upper time limit */ + if ($maxdt > 0.0) { + + /* coarse steps */ + while (($sat->el >= 1.0) && ($t <= ($start + $maxdt))) { + $t += cos(($sat->el - 1.0) * self::de2ra) * sqrt($sat->alt) / 25000.0; + $this->predict_calc($sat, $qth, $t); + } + + /* fine steps */ + while (($lostime == 0.0) && ($t <= ($start + $maxdt))) { + + $t += $sat->el * sqrt($sat->alt) / 502500.0; + $this->predict_calc($sat, $qth, $t); + + if (abs($sat->el) < 0.005) { + $lostime = $t; + } + } + } else { + /* don't use upper limit */ + + /* coarse steps */ + while ($sat->el >= 1.0) { + $t += cos(($sat->el - 1.0) * self::de2ra) * sqrt($sat->alt) / 25000.0; + $this->predict_calc($sat, $qth, $t); + } + + /* fine steps */ + while ($lostime == 0.0) { + + $t += $sat->el * sqrt($sat->alt) / 502500.0; + $this->predict_calc($sat, $qth, $t); + + if (abs($sat->el) < 0.005) + $lostime = $t; + } + } + + return $lostime; + } + + /** Find AOS time of current pass. + * @param Predict_Sat $sat The satellite to find AOS for. + * @param Predict_QTH $qth The ground station. + * @param float $start Start time, prefereably now. + * @return The time of the previous AOS or 0.0 if the satellite has no AOS. + * + * This function can be used to find the AOS time in the past of the + * current pass. + */ + public function find_prev_aos(Predict_Sat $sat, Predict_QTH $qth, $start) + { + $aostime = $start; + + /* make sure current sat values are + in sync with the time + */ + $this->predict_calc($sat, $qth, $start); + + /* check whether satellite has aos */ + if (($sat->otype == Predict_SGPSDP::ORBIT_TYPE_GEO) || + ($sat->otype == Predict_SGPSDP::ORBIT_TYPE_DECAYED) || + !$this->has_aos($sat, $qth)) { + + return 0.0; + } + + while ($sat->el >= 0.0) { + $aostime -= 0.0005; // 0.75 min + $this->predict_calc($sat, $qth, $aostime); + } + + return $aostime; + } + + /** Determine whether satellite ever reaches AOS. + * @author John A. Magliacane, KD2BD + * @author Alexandru Csete, OZ9AEC + * @param Predict_Sat $sat The satellite data. + * @param Predict_QTH $qth The observer's location data + * @return bool true if the satellite will reach AOS, false otherwise. + * + */ + public function has_aos(Predict_Sat $sat, Predict_QTH $qth) + { + $retcode = false; + + /* FIXME */ + if ($sat->meanmo == 0.0) { + $retcode = false; + } else { + + /* xincl is already in RAD by select_ephemeris */ + $lin = $sat->tle->xincl; + if ($lin >= self::pio2) { + $lin = self::pi - $lin; + } + + $sma = 331.25 * exp(log(1440.0 / $sat->meanmo) * (2.0 / 3.0)); + $apogee = $sma * (1.0 + $sat->tle->eo) - self::xkmper; + + if ((acos(self::xkmper / ($apogee + self::xkmper)) + ($lin)) > abs($qth->lat * self::de2ra)) { + $retcode = true; + } else { + $retcode = false; + } + } + + return $retcode; + } + + /** Predict passes after a certain time. + * + * + * This function calculates num upcoming passes with AOS no earlier + * than t = start and not later that t = (start+maxdt). The function will + * repeatedly call get_pass until + * the number of predicted passes is equal to num, the time has reached + * limit or the get_pass function returns NULL. + * + * note For no time limit use maxdt = 0.0 + * + * note the data in sat will be corrupt (future) and must be refreshed + * by the caller, if the caller will need it later on (eg. if the caller + * is GtkSatList). + * + * note Prepending to a singly linked list is much faster than appending. + * Therefore, the elements are prepended whereafter the GSList is + * reversed + * + * + * @param Predict_Sat $sat The satellite data + * @param Predict_QTH $qth The observer's location data + * @param float $start The start julian date + * @param int $maxdt The max # of days to look + * @param int $num The max # of passes to get + * @return array of Predict_Pass instances if found, empty array otherwise + */ + public function get_passes(Predict_Sat $sat, Predict_QTH $qth, $start, $maxdt, $num = 0) + { + $passes = array(); + + /* if no number has been specified + set it to something big */ + if ($num == 0) { + $num = 100; + } + + $t = $start; + + for ($i = 0; $i < $num; $i++) { + $pass = $this->get_pass($sat, $qth, $t, $maxdt); + + if ($pass != null) { + $passes[] = $pass; + $t = $pass->los + 0.014; // +20 min + + /* if maxdt > 0.0 check whether we have reached t = start+maxdt + if yes finish predictions + */ + if (($maxdt > 0.0) && ($t >= ($start + $maxdt))) { + $i = $num; + } + } else { + /* we can't get any more passes */ + $i = $num; + } + } + + return $passes; + } + + /** + * Filters out visible passes and adds the visible aos, tca, los, and + * corresponding az and ele for each. + * + * @param array $passes The passes returned from get_passes() + * + * @author Bill Shupp + * @return array + */ + public function filterVisiblePasses(array $passes) + { + $filtered = array(); + + foreach ($passes as $result) { + // Dummy check + if ($result->vis[0] != 'V') { + continue; + } + + $aos = false; + $aos_az = false; + $aos = false; + $tca = false; + $los_az = false; + $aos_el = 0; + $max_el = 0; + + foreach ($result->details as $detail) { + if ($detail->vis != Predict::SAT_VIS_VISIBLE) { + continue; + } + if ($detail->el < $this->minEle) { + continue; + } + + if ($aos == false) { + $aos = $detail->time; + $aos_az = $detail->az; + $aos_el = $detail->el; + $tca = $detail->time; + $los = $detail->time; + $los_az = $detail->az; + $los_el = $detail->el; + $max_el = $detail->el; + $max_el_az = $detail->el; + continue; + } + $los = $detail->time; + $los_az = $detail->az; + $los_el = $detail->el; + + if ($detail->el > $max_el) { + $tca = $detail->time; + $max_el = $detail->el; + $max_el_az = $detail->az; + } + } + + if ($aos === false) { + // Does not reach minimum elevation, skip + continue; + } + + $result->visible_aos = $aos; + $result->visible_aos_az = $aos_az; + $result->visible_aos_el = $aos_el; + $result->visible_tca = $tca; + $result->visible_max_el = $max_el; + $result->visible_max_el_az = $max_el_az; + $result->visible_los = $los; + $result->visible_los_az = $los_az; + $result->visible_los_el = $los_el; + + $filtered[] = $result; + } + + return $filtered; + } + + /** + * Translates aziumuth degrees to compass direction: + * + * N (0°), NNE (22.5°), NE (45°), ENE (67.5°), E (90°), ESE (112.5°), + * SE (135°), SSE (157.5°), S (180°), SSW (202.5°), SW (225°), + * WSW (247.5°), W (270°), WNW (292.5°), NW (315°), NNW (337.5°) + * + * @param int $az The azimuth in degrees, defaults to 0 + * + * @return string + */ + public function azDegreesToDirection($az = 0) + { + $i = floor($az / 22.5); + $m = (22.5 * (2 * $i + 1)) / 2; + $i = ($az >= $m) ? $i + 1 : $i; + + return trim(substr('N NNENE ENEE ESESE SSES SSWSW WSWW WNWNW NNWN ', $i * 3, 3)); + } +} diff --git a/predict/Predict/DeepArg.php b/predict/Predict/DeepArg.php new file mode 100644 index 000000000..f29178fca --- /dev/null +++ b/predict/Predict/DeepArg.php @@ -0,0 +1,35 @@ + 0 ) { + return 1; + } else if ($arg < 0 ) { + return -1; + } else { + return 0; + } + } + + /* Returns the arcsine of the argument */ + public static function ArcSin($arg) + { + if (abs($arg) >= 1 ) { + return (self::Sign($arg) * Predict::pio2); + } else { + return(atan($arg / sqrt(1 - $arg * $arg))); + } + } + + /* Returns arccosine of rgument */ + public static function ArcCos($arg) + { + return Predict::pio2 - self::ArcSin($arg); + } + + /* Adds vectors v1 and v2 together to produce v3 */ + public static function Vec_Add(Predict_Vector $v1, Predict_Vector $v2, Predict_Vector $v3) + { + $v3->x = $v1->x + $v2->x; + $v3->y = $v1->y + $v2->y; + $v3->z = $v1->z + $v2->z; + + $v3->w = sqrt($v3->x * $v3->x + $v3->y * $v3->y + $v3->z * $v3->z); + } + + /* Subtracts vector v2 from v1 to produce v3 */ + public static function Vec_Sub(Predict_Vector $v1, Predict_Vector $v2, Predict_Vector $v3) + { + $v3->x = $v1->x - $v2->x; + $v3->y = $v1->y - $v2->y; + $v3->z = $v1->z - $v2->z; + + $v3->w = sqrt($v3->x * $v3->x + $v3->y * $v3->y + $v3->z * $v3->z); + } + + /* Multiplies the vector v1 by the scalar k to produce the vector v2 */ + public static function Scalar_Multiply($k, Predict_Vector $v1, Predict_Vector $v2) + { + $v2->x = $k * $v1->x; + $v2->y = $k * $v1->y; + $v2->z = $k * $v1->z; + $v2->w = abs($k) * $v1->w; + } + + /* Multiplies the vector v1 by the scalar k */ + public static function Scale_Vector($k, Predict_Vector $v) + { + $v->x *= $k; + $v->y *= $k; + $v->z *= $k; + + $v->w = sqrt($v->x * $v->x + $v->y * $v->y + $v->z * $v->z); + } + + /* Returns the dot product of two vectors */ + public static function Dot(Predict_Vector $v1, Predict_Vector $v2) + { + return ($v1->x * $v2->x + $v1->y * $v2->y + $v1->z * $v2->z); + } + + /* Calculates the angle between vectors v1 and v2 */ + public static function Angle(Predict_Vector $v1, Predict_Vector $v2) + { + $v1->w = sqrt($v1->x * $v1->x + $v1->y * $v1->y + $v1->z * $v1->z); + $v2->w = sqrt($v2->x * $v2->x + $v2->y * $v2->y + $v2->z * $v2->z); + return (self::ArcCos(self::Dot($v1, $v2) / ($v1->w * $v2->w))); + } + + /* Produces cross product of v1 and v2, and returns in v3 */ + public static function Cross(Predict_Vector $v1, Predict_Vector $v2 ,Predict_Vector $v3) + { + $v3->x = $v1->y * $v2->z - $v1->z * $v2->y; + $v3->y = $v1->z * $v2->x - $v1->x * $v2->z; + $v3->z = $v1->x * $v2->y - $v1->y * $v2->x; + + $v3->w = sqrt($v3->x * $v3->x + $v3->y * $v3->y + $v3->z * $v3->z); + } + + /* Normalizes a vector */ + public static function Normalize(Predict_Vector $v ) + { + $v->x /= $v->w; + $v->y /= $v->w; + $v->z /= $v->w; + } + + /* Four-quadrant arctan function */ + public static function AcTan($sinx, $cosx) + { + if ($cosx == 0) { + if ($sinx > 0) { + return Predict::pio2; + } else { + return Predict::x3pio2; + } + } else { + if ($cosx > 0) { + if ($sinx > 0) { + return atan($sinx / $cosx); + } else { + return Predict::twopi + atan($sinx / $cosx); + } + } else { + return Predict::pi + atan($sinx / $cosx); + } + } + } + + /* Returns mod 2pi of argument */ + public static function FMod2p($x) + { + $ret_val = $x; + $i = (int) ($ret_val / Predict::twopi); + $ret_val -= $i * Predict::twopi; + + if ($ret_val < 0) { + $ret_val += Predict::twopi; + } + + return $ret_val; + } + + /* Returns arg1 mod arg2 */ + public static function Modulus($arg1, $arg2) + { + $ret_val = $arg1; + $i = (int) ($ret_val / $arg2); + $ret_val -= $i * $arg2; + + if ($ret_val < 0) { + $ret_val += $arg2; + } + + return $ret_val; + } + + /* Returns fractional part of double argument */ + public static function Frac($arg) + { + return $arg - floor($arg); + } + + /* Converts the satellite's position and velocity */ + /* vectors from normalised values to km and km/sec */ + public static function Convert_Sat_State(Predict_Vector $pos, Predict_Vector $vel) + { + self::Scale_Vector(Predict::xkmper, $pos); + self::Scale_Vector(Predict::xkmper * Predict::xmnpda / Predict::secday, $vel); + } + + /* Returns angle in radians from arg in degrees */ + public static function Radians($arg) + { + return $arg * Predict::de2ra; + } + + /* Returns angle in degrees from arg in rads */ + public static function Degrees($arg) + { + return $arg / Predict::de2ra; + } +} diff --git a/predict/Predict/ObsSet.php b/predict/Predict/ObsSet.php new file mode 100644 index 000000000..9470b93d4 --- /dev/null +++ b/predict/Predict/ObsSet.php @@ -0,0 +1,12 @@ +pos = new Predict_Vector(); + $this->vel = new Predict_Vector(); + } +} diff --git a/predict/Predict/QTH.php b/predict/Predict/QTH.php new file mode 100644 index 000000000..5ccd5d365 --- /dev/null +++ b/predict/Predict/QTH.php @@ -0,0 +1,20 @@ +lat); /* Only run sin($geodetic->lat) once */ + + $geodetic->theta = Predict_Math::FMod2p(Predict_Time::ThetaG_JD($_time) + $geodetic->lon);/*LMST*/ + $c = 1 / sqrt(1 + Predict::__f * (Predict::__f - 2) * $sinGeodeticLat * $sinGeodeticLat); + $sq = (1 - Predict::__f) * (1 - Predict::__f) * $c; + $achcp = (Predict::xkmper * $c + $geodetic->alt) * cos($geodetic->lat); + $obs_pos->x = $achcp * cos($geodetic->theta); /*kilometers*/ + $obs_pos->y = $achcp * sin($geodetic->theta); + $obs_pos->z = (Predict::xkmper * $sq + $geodetic->alt) * $sinGeodeticLat; + $obs_vel->x = -Predict::mfactor * $obs_pos->y; /*kilometers/second*/ + $obs_vel->y = Predict::mfactor * $obs_pos->x; + $obs_vel->z = 0; + $obs_pos->w = sqrt($obs_pos->x * $obs_pos->x + $obs_pos->y * $obs_pos->y + $obs_pos->z * $obs_pos->z); + $obs_vel->w = sqrt($obs_vel->x * $obs_vel->x + $obs_vel->y * $obs_vel->y + $obs_vel->z * $obs_vel->z); + } + + /* Procedure Calculate_LatLonAlt will calculate the geodetic */ + /* position of an object given its ECI position pos and time. */ + /* It is intended to be used to determine the ground track of */ + /* a satellite. The calculations assume the earth to be an */ + /* oblate spheroid as defined in WGS '72. */ + public static function Calculate_LatLonAlt($_time, Predict_Vector $pos, Predict_Geodetic $geodetic) + { + /* Reference: The 1992 Astronomical Almanac, page K12. */ + + /* double r,e2,phi,c; */ + + $geodetic->theta = Predict_Math::AcTan($pos->y, $pos->x); /*radians*/ + $geodetic->lon = Predict_Math::FMod2p($geodetic->theta - Predict_Time::ThetaG_JD($_time)); /*radians*/ + $r = sqrt(($pos->x * $pos->x) + ($pos->y * $pos->y)); + $e2 = Predict::__f * (2 - Predict::__f); + $geodetic->lat = Predict_Math::AcTan($pos->z, $r); /*radians*/ + + do { + $phi = $geodetic->lat; + $sinPhi = sin($phi); + $c = 1 / sqrt(1 - $e2 * ($sinPhi * $sinPhi)); + $geodetic->lat = Predict_Math::AcTan($pos->z + Predict::xkmper * $c * $e2 * $sinPhi, $r); + } while (abs($geodetic->lat - $phi) >= 1E-10); + + $geodetic->alt = $r / cos($geodetic->lat) - Predict::xkmper * $c;/*kilometers*/ + + if ($geodetic->lat > Predict::pio2) { + $geodetic->lat -= Predict::twopi; + } + } + + /* The procedures Calculate_Obs and Calculate_RADec calculate */ + /* the *topocentric* coordinates of the object with ECI position, */ + /* {pos}, and velocity, {vel}, from location {geodetic} at {time}. */ + /* The {obs_set} returned for Calculate_Obs consists of azimuth, */ + /* elevation, range, and range rate (in that order) with units of */ + /* radians, radians, kilometers, and kilometers/second, respectively. */ + /* The WGS '72 geoid is used and the effect of atmospheric refraction */ + /* (under standard temperature and pressure) is incorporated into the */ + /* elevation calculation; the effect of atmospheric refraction on */ + /* range and range rate has not yet been quantified. */ + + /* The {obs_set} for Calculate_RADec consists of right ascension and */ + /* declination (in that order) in radians. Again, calculations are */ + /* based on *topocentric* position using the WGS '72 geoid and */ + /* incorporating atmospheric refraction. */ + public static function Calculate_Obs($_time, Predict_Vector $pos, Predict_Vector $vel, Predict_Geodetic $geodetic, Predict_ObsSet $obs_set) + { + $obs_pos = new Predict_Vector(); + $obs_vel = new Predict_Vector(); + $range = new Predict_Vector(); + $rgvel = new Predict_Vector(); + + self::Calculate_User_PosVel($_time, $geodetic, $obs_pos, $obs_vel); + + $range->x = $pos->x - $obs_pos->x; + $range->y = $pos->y - $obs_pos->y; + $range->z = $pos->z - $obs_pos->z; + + $rgvel->x = $vel->x - $obs_vel->x; + $rgvel->y = $vel->y - $obs_vel->y; + $rgvel->z = $vel->z - $obs_vel->z; + + $range->w = sqrt($range->x * $range->x + $range->y * $range->y + $range->z * $range->z); + + $sin_lat = sin($geodetic->lat); + $cos_lat = cos($geodetic->lat); + $sin_theta = sin($geodetic->theta); + $cos_theta = cos($geodetic->theta); + $top_s = $sin_lat * $cos_theta * $range->x + + $sin_lat * $sin_theta * $range->y + - $cos_lat * $range->z; + $top_e = -$sin_theta * $range->x + + $cos_theta * $range->y; + $top_z = $cos_lat * $cos_theta * $range->x + + $cos_lat * $sin_theta * $range->y + + $sin_lat * $range->z; + $azim = atan(-$top_e / $top_s); /*Azimuth*/ + if ($top_s > 0) { + $azim = $azim + Predict::pi; + } + if ($azim < 0 ) { + $azim = $azim + Predict::twopi; + } + $el = Predict_Math::ArcSin($top_z / $range->w); + $obs_set->az = $azim; /* Azimuth (radians) */ + $obs_set->el = $el; /* Elevation (radians)*/ + $obs_set->range = $range->w; /* Range (kilometers) */ + + /* Range Rate (kilometers/second)*/ + $obs_set->range_rate = Predict_Math::Dot($range, $rgvel) / $range->w; + + /* Corrections for atmospheric refraction */ + /* Reference: Astronomical Algorithms by Jean Meeus, pp. 101-104 */ + /* Correction is meaningless when apparent elevation is below horizon */ + // obs_set->el = obs_set->el + Radians((1.02/tan(Radians(Degrees(el)+ + // 10.3/(Degrees(el)+5.11))))/60); + if ($obs_set->el < 0) { + $obs_set->el = $el; /*Reset to true elevation*/ + } + } +} diff --git a/predict/Predict/SGPSDP.php b/predict/Predict/SGPSDP.php new file mode 100644 index 000000000..41a7de548 --- /dev/null +++ b/predict/Predict/SGPSDP.php @@ -0,0 +1,1061 @@ +flags & self::SGP4_INITIALIZED_FLAG) { + $sat->flags |= self::SGP4_INITIALIZED_FLAG; + + /* Recover original mean motion (xnodp) and */ + /* semimajor axis (aodp) from input elements. */ + $a1 = pow(Predict::xke / $sat->tle->xno, Predict::tothrd); + $sat->sgps->cosio = cos($sat->tle->xincl); + $theta2 = $sat->sgps->cosio * $sat->sgps->cosio; + $sat->sgps->x3thm1 = 3 * $theta2 - 1.0; + $eosq = $sat->tle->eo * $sat->tle->eo; + $betao2 = 1 - $eosq; + $betao = sqrt($betao2); + $del1 = 1.5 * Predict::ck2 * $sat->sgps->x3thm1 / ($a1 * $a1 * $betao * $betao2); + $ao = $a1 * (1 - $del1 * (0.5 * Predict::tothrd + $del1 * (1 + 134.0 / 81.0 * $del1))); + $delo = 1.5 * Predict::ck2 * $sat->sgps->x3thm1 / ($ao * $ao * $betao * $betao2); + $sat->sgps->xnodp = $sat->tle->xno / (1.0 + $delo); + $sat->sgps->aodp = $ao / (1.0 - $delo); + + /* For perigee less than 220 kilometers, the "simple" flag is set */ + /* and the equations are truncated to linear variation in sqrt a */ + /* and quadratic variation in mean anomaly. Also, the c3 term, */ + /* the delta omega term, and the delta m term are dropped. */ + if (($sat->sgps->aodp * (1.0 - $sat->tle->eo) / Predict::ae) < (220.0 / Predict::xkmper + Predict::ae)) { + $sat->flags |= self::SIMPLE_FLAG; + } else { + $sat->flags &= ~self::SIMPLE_FLAG; + } + + /* For perigee below 156 km, the */ + /* values of s and qoms2t are altered. */ + $s4 = Predict::__s__; + $qoms24 = Predict::qoms2t; + $perige = ($sat->sgps->aodp * (1 - $sat->tle->eo) - Predict::ae) * Predict::xkmper; + if ($perige < 156.0) { + if ($perige <= 98.0) { + $s4 = 20.0; + } else { + $s4 = $perige - 78.0; + } + $qoms24 = pow((120.0 - $s4) * Predict::ae / Predict::xkmper, 4); + $s4 = $s4 / Predict::xkmper + Predict::ae; + }; /* FIXME FIXME: End of if(perige <= 98) NO WAY!!!! */ + + $pinvsq = 1.0 / ($sat->sgps->aodp * $sat->sgps->aodp * $betao2 * $betao2); + $tsi = 1.0 / ($sat->sgps->aodp - $s4); + $sat->sgps->eta = $sat->sgps->aodp * $sat->tle->eo * $tsi; + $etasq = $sat->sgps->eta * $sat->sgps->eta; + $eeta = $sat->tle->eo * $sat->sgps->eta; + $psisq = abs(1.0 - $etasq); + $coef = $qoms24 * pow($tsi, 4); + $coef1 = $coef / pow($psisq, 3.5); + $c2 = $coef1 * $sat->sgps->xnodp * ($sat->sgps->aodp * + (1.0 + 1.5 * $etasq + $eeta * (4.0 + $etasq)) + + 0.75 * Predict::ck2 * $tsi / $psisq * $sat->sgps->x3thm1 * + (8.0 + 3.0 * $etasq * (8 + $etasq))); + $sat->sgps->c1 = $c2 * $sat->tle->bstar; + $sat->sgps->sinio = sin($sat->tle->xincl); + $a3ovk2 = -Predict::xj3 / Predict::ck2 * pow(Predict::ae, 3); + $c3 = $coef * $tsi * $a3ovk2 * $sat->sgps->xnodp * Predict::ae * $sat->sgps->sinio / $sat->tle->eo; + $sat->sgps->x1mth2 = 1.0 - $theta2; + $sat->sgps->c4 = 2.0 * $sat->sgps->xnodp * $coef1 * $sat->sgps->aodp * $betao2 * + ($sat->sgps->eta * (2.0 + 0.5 * $etasq) + + $sat->tle->eo * (0.5 + 2.0 * $etasq) - + 2.0 * Predict::ck2 * $tsi / ($sat->sgps->aodp * $psisq) * + (-3.0 * $sat->sgps->x3thm1 * (1.0 - 2.0 * $eeta + $etasq * (1.5 - 0.5 * $eeta)) + + 0.75 * $sat->sgps->x1mth2 * (2.0 * $etasq - $eeta * (1.0 + $etasq)) * + cos(2.0 * $sat->tle->omegao))); + $sat->sgps->c5 = 2.0 * $coef1 * $sat->sgps->aodp * $betao2 * + (1.0 + 2.75 * ($etasq + $eeta) + $eeta * $etasq); + $theta4 = $theta2 * $theta2; + $temp1 = 3.0 * Predict::ck2 * $pinvsq * $sat->sgps->xnodp; + $temp2 = $temp1 * Predict::ck2 * $pinvsq; + $temp3 = 1.25 * Predict::ck4 * $pinvsq * $pinvsq * $sat->sgps->xnodp; + $sat->sgps->xmdot = $sat->sgps->xnodp + 0.5 * $temp1 * $betao * $sat->sgps->x3thm1 + + 0.0625 * $temp2 * $betao * (13.0 - 78.0 * $theta2 + 137.0 * $theta4); + $x1m5th = 1.0 - 5.0 * $theta2; + $sat->sgps->omgdot = -0.5 * $temp1 * $x1m5th + + 0.0625 * $temp2 * (7.0 - 114.0 * $theta2 + 395.0 * $theta4) + + $temp3 * (3.0 - 36.0 * $theta2 + 49.0 * $theta4); + $xhdot1 = -$temp1 * $sat->sgps->cosio; + $sat->sgps->xnodot = $xhdot1 + (0.5 * $temp2 * (4.0 - 19.0 * $theta2) + + 2.0 * $temp3 * (3.0 - 7.0 * $theta2)) * $sat->sgps->cosio; + $sat->sgps->omgcof = $sat->tle->bstar * $c3 * cos($sat->tle->omegao); + $sat->sgps->xmcof = -Predict::tothrd * $coef * $sat->tle->bstar * Predict::ae / $eeta; + $sat->sgps->xnodcf = 3.5 * $betao2 * $xhdot1 * $sat->sgps->c1; + $sat->sgps->t2cof = 1.5 * $sat->sgps->c1; + $sat->sgps->xlcof = 0.125 * $a3ovk2 * $sat->sgps->sinio * + (3.0 + 5.0 * $sat->sgps->cosio) / (1.0 + $sat->sgps->cosio); + $sat->sgps->aycof = 0.25 * $a3ovk2 * $sat->sgps->sinio; + $sat->sgps->delmo = pow(1.0 + $sat->sgps->eta * cos($sat->tle->xmo), 3); + $sat->sgps->sinmo = sin($sat->tle->xmo); + $sat->sgps->x7thm1 = 7.0 * $theta2 - 1.0; + if (~$sat->flags & self::SIMPLE_FLAG) { + $c1sq = $sat->sgps->c1 * $sat->sgps->c1; + $sat->sgps->d2 = 4.0 * $sat->sgps->aodp * $tsi * $c1sq; + $temp = $sat->sgps->d2 * $tsi * $sat->sgps->c1 / 3.0; + $sat->sgps->d3 = (17.0 * $sat->sgps->aodp + $s4) * $temp; + $sat->sgps->d4 = 0.5 * $temp * $sat->sgps->aodp * $tsi * + (221.0 * $sat->sgps->aodp + 31.0 * $s4) * $sat->sgps->c1; + $sat->sgps->t3cof = $sat->sgps->d2 + 2.0 * $c1sq; + $sat->sgps->t4cof = 0.25 * (3.0 * $sat->sgps->d3 + $sat->sgps->c1 * + (12.0 * $sat->sgps->d2 + 10.0 * $c1sq)); + $sat->sgps->t5cof = 0.2 * (3.0 * $sat->sgps->d4 + + 12.0 * $sat->sgps->c1 * $sat->sgps->d3 + + 6.0 * $sat->sgps->d2 * $sat->sgps->d2 + + 15.0 * $c1sq * (2.0 * $sat->sgps->d2 + $c1sq)); + }; /* End of if (isFlagClear(SIMPLE_FLAG)) */ + }; /* End of SGP4() initialization */ + + /* Update for secular gravity and atmospheric drag. */ + $xmdf = $sat->tle->xmo + $sat->sgps->xmdot * $tsince; + $omgadf = $sat->tle->omegao + $sat->sgps->omgdot * $tsince; + $xnoddf = $sat->tle->xnodeo + $sat->sgps->xnodot * $tsince; + $omega = $omgadf; + $xmp = $xmdf; + $tsq = $tsince * $tsince; + $xnode = $xnoddf + $sat->sgps->xnodcf * $tsq; + $tempa = 1.0 - $sat->sgps->c1 * $tsince; + $tempe = $sat->tle->bstar * $sat->sgps->c4 * $tsince; + $templ = $sat->sgps->t2cof * $tsq; + if (~$sat->flags & self::SIMPLE_FLAG) { + $delomg = $sat->sgps->omgcof * $tsince; + $delm = $sat->sgps->xmcof * (pow(1 + $sat->sgps->eta * cos($xmdf), 3) - $sat->sgps->delmo); + $temp = $delomg + $delm; + $xmp = $xmdf + $temp; + $omega = $omgadf - $temp; + $tcube = $tsq * $tsince; + $tfour = $tsince * $tcube; + $tempa = $tempa - $sat->sgps->d2 * $tsq - $sat->sgps->d3 * $tcube - $sat->sgps->d4 * $tfour; + $tempe = $tempe + $sat->tle->bstar * $sat->sgps->c5 * (sin($xmp) - $sat->sgps->sinmo); + $templ = $templ + $sat->sgps->t3cof * $tcube + $tfour * + ($sat->sgps->t4cof + $tsince * $sat->sgps->t5cof); + }; /* End of if (isFlagClear(SIMPLE_FLAG)) */ + + $a = $sat->sgps->aodp * pow($tempa, 2); + $e = $sat->tle->eo - $tempe; + $xl = $xmp + $omega + $xnode + $sat->sgps->xnodp * $templ; + $beta = sqrt(1.0 - ($e * $e)); + $xn = Predict::xke / pow($a, 1.5); + + /* Long period periodics */ + $axn = $e * cos($omega); + $temp = 1.0 / ($a * $beta * $beta); + $xll = $temp * $sat->sgps->xlcof * $axn; + $aynl = $temp * $sat->sgps->aycof; + $xlt = $xl + $xll; + $ayn = $e * sin($omega) + $aynl; + + /* Solve Kepler's' Equation */ + $capu = Predict_Math::FMod2p($xlt - $xnode); + $temp2 = $capu; + + $i = 0; + do { + $sinepw = sin($temp2); + $cosepw = cos($temp2); + $temp3 = $axn * $sinepw; + $temp4 = $ayn * $cosepw; + $temp5 = $axn * $cosepw; + $temp6 = $ayn * $sinepw; + $epw = ($capu - $temp4 + $temp3 - $temp2) / (1.0 - $temp5 - $temp6) + $temp2; + if (abs($epw - $temp2) <= Predict::e6a) { + break; + } + $temp2 = $epw; + } while ($i++ < 10); + + /* Short period preliminary quantities */ + $ecose = $temp5 + $temp6; + $esine = $temp3 - $temp4; + $elsq = $axn * $axn + $ayn * $ayn; + $temp = 1.0 - $elsq; + $pl = $a * $temp; + $r = $a * (1.0 - $ecose); + $temp1 = 1.0 / $r; + $rdot = Predict::xke * sqrt($a) * $esine * $temp1; + $rfdot = Predict::xke * sqrt($pl) * $temp1; + $temp2 = $a * $temp1; + $betal = sqrt($temp); + $temp3 = 1.0 / (1.0 + $betal); + $cosu = $temp2 * ($cosepw - $axn + $ayn * $esine * $temp3); + $sinu = $temp2 * ($sinepw - $ayn - $axn * $esine * $temp3); + $u = Predict_Math::AcTan($sinu, $cosu); + $sin2u = 2.0 * $sinu * $cosu; + $cos2u = 2.0 * $cosu * $cosu - 1.0; + $temp = 1.0 / $pl; + $temp1 = Predict::ck2 * $temp; + $temp2 = $temp1 * $temp; + + /* Update for short periodics */ + $rk = $r * (1.0 - 1.5 * $temp2 * $betal * $sat->sgps->x3thm1) + + 0.5 * $temp1 * $sat->sgps->x1mth2 * $cos2u; + $uk = $u - 0.25 * $temp2 * $sat->sgps->x7thm1 * $sin2u; + $xnodek = $xnode + 1.5 * $temp2 * $sat->sgps->cosio * $sin2u; + $xinck = $sat->tle->xincl + 1.5 * $temp2 * $sat->sgps->cosio * $sat->sgps->sinio * $cos2u; + $rdotk = $rdot - $xn * $temp1 * $sat->sgps->x1mth2 * $sin2u; + $rfdotk = $rfdot + $xn * $temp1 * ($sat->sgps->x1mth2 * $cos2u + 1.5 * $sat->sgps->x3thm1); + + + /* Orientation vectors */ + $sinuk = sin($uk); + $cosuk = cos($uk); + $sinik = sin($xinck); + $cosik = cos($xinck); + $sinnok = sin($xnodek); + $cosnok = cos($xnodek); + $xmx = -$sinnok * $cosik; + $xmy = $cosnok * $cosik; + $ux = $xmx * $sinuk + $cosnok * $cosuk; + $uy = $xmy * $sinuk + $sinnok * $cosuk; + $uz = $sinik * $sinuk; + $vx = $xmx * $cosuk - $cosnok * $sinuk; + $vy = $xmy * $cosuk - $sinnok * $sinuk; + $vz = $sinik * $cosuk; + + /* Position and velocity */ + $sat->pos->x = $rk * $ux; + $sat->pos->y = $rk * $uy; + $sat->pos->z = $rk * $uz; + $sat->vel->x = $rdotk * $ux + $rfdotk * $vx; + $sat->vel->y = $rdotk * $uy + $rfdotk * $vy; + $sat->vel->z = $rdotk * $uz + $rfdotk * $vz; + + $sat->phase = $xlt - $xnode - $omgadf + Predict::twopi; + if ($sat->phase < 0) { + $sat->phase += Predict::twopi; + } + $sat->phase = Predict_Math::FMod2p($sat->phase); + + $sat->tle->omegao1 = $omega; + $sat->tle->xincl1 = $xinck; + $sat->tle->xnodeo1 = $xnodek; + + } /*SGP4*/ + + /* SDP4 */ + /* This function is used to calculate the position and velocity */ + /* of deep-space (period > 225 minutes) satellites. tsince is */ + /* time since epoch in minutes, tle is a pointer to a tle_t */ + /* structure with Keplerian orbital elements and pos and vel */ + /* are vector_t structures returning ECI satellite position and */ + /* velocity. Use Convert_Sat_State() to convert to km and km/s. */ + public function SDP4(Predict_Sat $sat, $tsince) + { + /* Initialization */ + if (~$sat->flags & self::SDP4_INITIALIZED_FLAG) { + + $sat->flags |= self::SDP4_INITIALIZED_FLAG; + + /* Recover original mean motion (xnodp) and */ + /* semimajor axis (aodp) from input elements. */ + $a1 = pow(Predict::xke / $sat->tle->xno, Predict::tothrd); + $sat->deep_arg->cosio = cos($sat->tle->xincl); + $sat->deep_arg->theta2 = $sat->deep_arg->cosio * $sat->deep_arg->cosio; + $sat->sgps->x3thm1 = 3.0 * $sat->deep_arg->theta2 - 1.0; + $sat->deep_arg->eosq = $sat->tle->eo * $sat->tle->eo; + $sat->deep_arg->betao2 = 1.0 - $sat->deep_arg->eosq; + $sat->deep_arg->betao = sqrt($sat->deep_arg->betao2); + $del1 = 1.5 * Predict::ck2 * $sat->sgps->x3thm1 / + ($a1 * $a1 * $sat->deep_arg->betao * $sat->deep_arg->betao2); + $ao = $a1 * (1.0 - $del1 * (0.5 * Predict::tothrd + $del1 * (1.0 + 134.0 / 81.0 * $del1))); + $delo = 1.5 * Predict::ck2 * $sat->sgps->x3thm1 / + ($ao * $ao * $sat->deep_arg->betao * $sat->deep_arg->betao2); + $sat->deep_arg->xnodp = $sat->tle->xno / (1.0 + $delo); + $sat->deep_arg->aodp = $ao / (1.0 - $delo); + + /* For perigee below 156 km, the values */ + /* of s and qoms2t are altered. */ + $s4 = Predict::__s__; + $qoms24 = Predict::qoms2t; + $perige = ($sat->deep_arg->aodp * (1.0 - $sat->tle->eo) - Predict::ae) * Predict::xkmper; + if ($perige < 156.0) { + if ($perige <= 98.0) { + $s4 = 20.0; + } else { + $s4 = $perige - 78.0; + } + $qoms24 = pow((120.0 - $s4) * Predict::ae / Predict::xkmper, 4); + $s4 = $s4 / Predict::xkmper + Predict::ae; + } + $pinvsq = 1.0 / ($sat->deep_arg->aodp * $sat->deep_arg->aodp * + $sat->deep_arg->betao2 * $sat->deep_arg->betao2); + $sat->deep_arg->sing = sin($sat->tle->omegao); + $sat->deep_arg->cosg = cos($sat->tle->omegao); + $tsi = 1.0 / ($sat->deep_arg->aodp - $s4); + $eta = $sat->deep_arg->aodp * $sat->tle->eo * $tsi; + $etasq = $eta * $eta; + $eeta = $sat->tle->eo * $eta; + $psisq = abs(1.0 - $etasq); + $coef = $qoms24 * pow($tsi, 4); + $coef1 = $coef / pow($psisq, 3.5); + $c2 = $coef1 * $sat->deep_arg->xnodp * ($sat->deep_arg->aodp * + (1.0 + 1.5 * $etasq + $eeta * + (4.0 + $etasq)) + 0.75 * Predict::ck2 * $tsi / $psisq * + $sat->sgps->x3thm1 * (8.0 + 3.0 * $etasq * + (8.0 + $etasq))); + $sat->sgps->c1 = $sat->tle->bstar * $c2; + $sat->deep_arg->sinio = sin($sat->tle->xincl); + $a3ovk2 = -Predict::xj3 / Predict::ck2 * pow(Predict::ae, 3); + $sat->sgps->x1mth2 = 1.0 - $sat->deep_arg->theta2; + $sat->sgps->c4 = 2.0 * $sat->deep_arg->xnodp * $coef1 * + $sat->deep_arg->aodp * $sat->deep_arg->betao2 * + ($eta * (2.0 + 0.5 * $etasq) + $sat->tle->eo * + (0.5 + 2.0 * $etasq) - 2.0 * Predict::ck2 * $tsi / + ($sat->deep_arg->aodp * $psisq) * (-3.0 * $sat->sgps->x3thm1 * + (1.0 - 2.0 * $eeta + $etasq * + (1.5 - 0.5 * $eeta)) + + 0.75 * $sat->sgps->x1mth2 * + (2.0 * $etasq - $eeta * (1.0 + $etasq)) * + cos(2.0 * $sat->tle->omegao))); + $theta4 = $sat->deep_arg->theta2 * $sat->deep_arg->theta2; + $temp1 = 3.0 * Predict::ck2 * $pinvsq * $sat->deep_arg->xnodp; + $temp2 = $temp1 * Predict::ck2 * $pinvsq; + $temp3 = 1.25 * Predict::ck4 * $pinvsq * $pinvsq * $sat->deep_arg->xnodp; + $sat->deep_arg->xmdot = $sat->deep_arg->xnodp + 0.5 * $temp1 * $sat->deep_arg->betao * + $sat->sgps->x3thm1 + 0.0625 * $temp2 * $sat->deep_arg->betao * + (13.0 - 78.0 * $sat->deep_arg->theta2 + 137.0 * $theta4); + $x1m5th = 1.0 - 5.0 * $sat->deep_arg->theta2; + $sat->deep_arg->omgdot = -0.5 * $temp1 * $x1m5th + 0.0625 * $temp2 * + (7.0 - 114.0 * $sat->deep_arg->theta2 + 395.0 * $theta4) + + $temp3 * (3.0 - 36.0 * $sat->deep_arg->theta2 + 49.0 * $theta4); + $xhdot1 = -$temp1 * $sat->deep_arg->cosio; + $sat->deep_arg->xnodot = $xhdot1 + (0.5 * $temp2 * (4.0 - 19.0 * $sat->deep_arg->theta2) + + 2.0 * $temp3 * (3.0 - 7.0 * $sat->deep_arg->theta2)) * + $sat->deep_arg->cosio; + $sat->sgps->xnodcf = 3.5 * $sat->deep_arg->betao2 * $xhdot1 * $sat->sgps->c1; + $sat->sgps->t2cof = 1.5 * $sat->sgps->c1; + $sat->sgps->xlcof = 0.125 * $a3ovk2 * $sat->deep_arg->sinio * + (3.0 + 5.0 * $sat->deep_arg->cosio) / (1.0 + $sat->deep_arg->cosio); + $sat->sgps->aycof = 0.25 * $a3ovk2 * $sat->deep_arg->sinio; + $sat->sgps->x7thm1 = 7.0 * $sat->deep_arg->theta2 - 1.0; + + /* initialize Deep() */ + $this->Deep(self::dpinit, $sat); + }; /*End of SDP4() initialization */ + + /* Update for secular gravity and atmospheric drag */ + $xmdf = $sat->tle->xmo + $sat->deep_arg->xmdot * $tsince; + $sat->deep_arg->omgadf = $sat->tle->omegao + $sat->deep_arg->omgdot * $tsince; + $xnoddf = $sat->tle->xnodeo + $sat->deep_arg->xnodot * $tsince; + $tsq = $tsince * $tsince; + $sat->deep_arg->xnode = $xnoddf + $sat->sgps->xnodcf * $tsq; + $tempa = 1.0 - $sat->sgps->c1 * $tsince; + $tempe = $sat->tle->bstar * $sat->sgps->c4 * $tsince; + $templ = $sat->sgps->t2cof * $tsq; + $sat->deep_arg->xn = $sat->deep_arg->xnodp; + + /* Update for deep-space secular effects */ + $sat->deep_arg->xll = $xmdf; + $sat->deep_arg->t = $tsince; + + $this->Deep(self::dpsec, $sat); + + $xmdf = $sat->deep_arg->xll; + $a = pow(Predict::xke / $sat->deep_arg->xn, Predict::tothrd) * $tempa * $tempa; + $sat->deep_arg->em = $sat->deep_arg->em - $tempe; + $xmam = $xmdf + $sat->deep_arg->xnodp * $templ; + + /* Update for deep-space periodic effects */ + $sat->deep_arg->xll = $xmam; + + $this->Deep(self::dpper, $sat); + + $xmam = $sat->deep_arg->xll; + $xl = $xmam + $sat->deep_arg->omgadf + $sat->deep_arg->xnode; + $beta = sqrt(1.0 - $sat->deep_arg->em * $sat->deep_arg->em); + $sat->deep_arg->xn = Predict::xke / pow($a, 1.5); + + /* Long period periodics */ + $axn = $sat->deep_arg->em * cos($sat->deep_arg->omgadf); + $temp = 1.0 / ($a * $beta * $beta); + $xll = $temp * $sat->sgps->xlcof * $axn; + $aynl = $temp * $sat->sgps->aycof; + $xlt = $xl + $xll; + $ayn = $sat->deep_arg->em * sin($sat->deep_arg->omgadf) + $aynl; + + /* Solve Kepler's Equation */ + $capu = Predict_Math::FMod2p ($xlt - $sat->deep_arg->xnode); + $temp2 = $capu; + + $i = 0; + do { + $sinepw = sin($temp2); + $cosepw = cos($temp2); + $temp3 = $axn * $sinepw; + $temp4 = $ayn * $cosepw; + $temp5 = $axn * $cosepw; + $temp6 = $ayn * $sinepw; + $epw = ($capu - $temp4 + $temp3 - $temp2) / (1.0 - $temp5 - $temp6) + $temp2; + if (abs($epw - $temp2) <= Predict::e6a) { + break; + } + $temp2 = $epw; + } while ($i++ < 10); + + /* Short period preliminary quantities */ + $ecose = $temp5 + $temp6; + $esine = $temp3 - $temp4; + $elsq = $axn * $axn + $ayn * $ayn; + $temp = 1.0 - $elsq; + $pl = $a * $temp; + $r = $a * (1.0 - $ecose); + $temp1 = 1.0 / $r; + $rdot = Predict::xke * sqrt($a) * $esine * $temp1; + $rfdot = Predict::xke * sqrt($pl) * $temp1; + $temp2 = $a * $temp1; + $betal = sqrt($temp); + $temp3 = 1.0 / (1.0 + $betal); + $cosu = $temp2 * ($cosepw - $axn + $ayn * $esine * $temp3); + $sinu = $temp2 * ($sinepw - $ayn - $axn * $esine * $temp3); + $u = Predict_Math::AcTan($sinu, $cosu); + $sin2u = 2.0 * $sinu * $cosu; + $cos2u = 2.0 * $cosu * $cosu - 1.0; + $temp = 1.0 / $pl; + $temp1 = Predict::ck2 * $temp; + $temp2 = $temp1 * $temp; + + /* Update for short periodics */ + $rk = $r * (1.0 - 1.5 * $temp2 * $betal * $sat->sgps->x3thm1) + + 0.5 * $temp1 * $sat->sgps->x1mth2 * $cos2u; + $uk = $u - 0.25 * $temp2 * $sat->sgps->x7thm1 * $sin2u; + $xnodek = $sat->deep_arg->xnode + 1.5 * $temp2 * $sat->deep_arg->cosio * $sin2u; + $xinck = $sat->deep_arg->xinc + 1.5 * $temp2 * + $sat->deep_arg->cosio * $sat->deep_arg->sinio * $cos2u; + $rdotk = $rdot - $sat->deep_arg->xn * $temp1 * $sat->sgps->x1mth2 * $sin2u; + $rfdotk = $rfdot + $sat->deep_arg->xn * $temp1 * + ($sat->sgps->x1mth2 * $cos2u + 1.5 * $sat->sgps->x3thm1); + + /* Orientation vectors */ + $sinuk = sin($uk); + $cosuk = cos($uk); + $sinik = sin($xinck); + $cosik = cos($xinck); + $sinnok = sin($xnodek); + $cosnok = cos($xnodek); + $xmx = -$sinnok * $cosik; + $xmy = $cosnok * $cosik; + $ux = $xmx * $sinuk + $cosnok * $cosuk; + $uy = $xmy * $sinuk + $sinnok * $cosuk; + $uz = $sinik * $sinuk; + $vx = $xmx * $cosuk - $cosnok * $sinuk; + $vy = $xmy * $cosuk - $sinnok * $sinuk; + $vz = $sinik * $cosuk; + + /* Position and velocity */ + $sat->pos->x = $rk * $ux; + $sat->pos->y = $rk * $uy; + $sat->pos->z = $rk * $uz; + $sat->vel->x = $rdotk * $ux + $rfdotk * $vx; + $sat->vel->y = $rdotk * $uy + $rfdotk * $vy; + $sat->vel->z = $rdotk * $uz + $rfdotk * $vz; + + /* Phase in rads */ + $sat->phase = $xlt - $sat->deep_arg->xnode - $sat->deep_arg->omgadf + Predict::twopi; + if ($sat->phase < 0.0) { + $sat->phase += Predict::twopi; + } + $sat->phase = Predict_Math::FMod2p ($sat->phase); + + $sat->tle->omegao1 = $sat->deep_arg->omgadf; + $sat->tle->xincl1 = $sat->deep_arg->xinc; + $sat->tle->xnodeo1 = $sat->deep_arg->xnode; + } /* SDP4 */ + + + /* DEEP */ + /* This function is used by SDP4 to add lunar and solar */ + /* perturbation effects to deep-space orbit objects. */ + public function Deep($ientry, Predict_Sat $sat) + { + switch ($ientry) { + case self::dpinit : /* Entrance for deep space initialization */ + $sat->dps->thgr = Predict_Time::ThetaG($sat->tle->epoch, $sat->deep_arg); + $eq = $sat->tle->eo; + $sat->dps->xnq = $sat->deep_arg->xnodp; + $aqnv = 1.0 / $sat->deep_arg->aodp; + $sat->dps->xqncl = $sat->tle->xincl; + $xmao = $sat->tle->xmo; + $xpidot = $sat->deep_arg->omgdot + $sat->deep_arg->xnodot; + $sinq = sin($sat->tle->xnodeo); + $cosq = cos($sat->tle->xnodeo); + $sat->dps->omegaq = $sat->tle->omegao; + $sat->dps->preep = 0; + + /* Initialize lunar solar terms */ + $day = $sat->deep_arg->ds50 + 18261.5; /* Days since 1900 Jan 0.5 */ + if ($day != $sat->dps->preep) { + $sat->dps->preep = $day; + $xnodce = 4.5236020 - 9.2422029E-4 * $day; + $stem = sin($xnodce); + $ctem = cos($xnodce); + $sat->dps->zcosil = 0.91375164 - 0.03568096 * $ctem; + $sat->dps->zsinil = sqrt(1.0 - $sat->dps->zcosil * $sat->dps->zcosil); + $sat->dps->zsinhl = 0.089683511 * $stem / $sat->dps->zsinil; + $sat->dps->zcoshl = sqrt(1.0 - $sat->dps->zsinhl * $sat->dps->zsinhl); + $c = 4.7199672 + 0.22997150 * $day; + $gam = 5.8351514 + 0.0019443680 * $day; + $sat->dps->zmol = Predict_Math::FMod2p($c - $gam); + $zx = 0.39785416 * $stem / $sat->dps->zsinil; + $zy = $sat->dps->zcoshl * $ctem + 0.91744867 * $sat->dps->zsinhl * $stem; + $zx = Predict_Math::AcTan($zx, $zy); + $zx = $gam + $zx - $xnodce; + $sat->dps->zcosgl = cos($zx); + $sat->dps->zsingl = sin($zx); + $sat->dps->zmos = 6.2565837 + 0.017201977 * $day; + $sat->dps->zmos = Predict_Math::FMod2p($sat->dps->zmos); + } /* End if(day != preep) */ + + /* Do solar terms */ + $sat->dps->savtsn = 1E20; + $zcosg = Predict::zcosgs; + $zsing = Predict::zsings; + $zcosi = Predict::zcosis; + $zsini = Predict::zsinis; + $zcosh = $cosq; + $zsinh = $sinq; + $cc = Predict::c1ss; + $zn = Predict::zns; + $ze = Predict::zes; + $zmo = $sat->dps->zmos; + $xnoi = 1.0 / $sat->dps->xnq; + + /* Loop breaks when Solar terms are done a second */ + /* time, after Lunar terms are initialized */ + for(;;) { + /* Solar terms done again after Lunar terms are done */ + $a1 = $zcosg * $zcosh + $zsing * $zcosi * $zsinh; + $a3 = -$zsing * $zcosh + $zcosg * $zcosi * $zsinh; + $a7 = -$zcosg * $zsinh + $zsing * $zcosi * $zcosh; + $a8 = $zsing * $zsini; + $a9 = $zsing * $zsinh + $zcosg * $zcosi * $zcosh; + $a10 = $zcosg * $zsini; + $a2 = $sat->deep_arg->cosio * $a7 + $sat->deep_arg->sinio * $a8; + $a4 = $sat->deep_arg->cosio * $a9 + $sat->deep_arg->sinio * $a10; + $a5 = -$sat->deep_arg->sinio * $a7 + $sat->deep_arg->cosio * $a8; + $a6 = -$sat->deep_arg->sinio * $a9 + $sat->deep_arg->cosio * $a10; + $x1 = $a1 * $sat->deep_arg->cosg + $a2 * $sat->deep_arg->sing; + $x2 = $a3 * $sat->deep_arg->cosg + $a4 * $sat->deep_arg->sing; + $x3 = -$a1 * $sat->deep_arg->sing + $a2 * $sat->deep_arg->cosg; + $x4 = -$a3 * $sat->deep_arg->sing + $a4 * $sat->deep_arg->cosg; + $x5 = $a5 * $sat->deep_arg->sing; + $x6 = $a6 * $sat->deep_arg->sing; + $x7 = $a5 * $sat->deep_arg->cosg; + $x8 = $a6 * $sat->deep_arg->cosg; + $z31 = 12 * $x1 * $x1 - 3 * $x3 * $x3; + $z32 = 24 * $x1 * $x2 - 6 * $x3 * $x4; + $z33 = 12 * $x2 * $x2 - 3 * $x4 * $x4; + $z1 = 3 * ($a1 * $a1 + $a2 * $a2) + $z31 * $sat->deep_arg->eosq; + $z2 = 6 * ($a1 * $a3 + $a2 * $a4) + $z32 * $sat->deep_arg->eosq; + $z3 = 3 * ($a3 * $a3 + $a4 * $a4) + $z33 * $sat->deep_arg->eosq; + $z11 = -6 * $a1 * $a5 + $sat->deep_arg->eosq * (-24 * $x1 * $x7 - 6 * $x3 * $x5); + $z12 = -6 * ($a1 * $a6 + $a3 * $a5) + $sat->deep_arg->eosq * + (-24 * ($x2 * $x7 + $x1 * $x8) - 6 * ($x3 * $x6 + $x4 * $x5)); + $z13 = -6 * $a3 * $a6 + $sat->deep_arg->eosq * (-24 * $x2 * $x8 - 6 * $x4 * $x6); + $z21 = 6 * $a2 * $a5 + $sat->deep_arg->eosq * (24 * $x1 * $x5 - 6 * $x3 * $x7); + $z22 = 6 * ($a4 * $a5 + $a2 * $a6) + $sat->deep_arg->eosq * + (24 * ($x2 * $x5 + $x1 * $x6) - 6 * ($x4 * $x7 + $x3 * $x8)); + $z23 = 6 * $a4 * $a6 + $sat->deep_arg->eosq * (24 * $x2 * $x6 - 6 * $x4 * $x8); + $z1 = $z1 + $z1 + $sat->deep_arg->betao2 * $z31; + $z2 = $z2 + $z2 + $sat->deep_arg->betao2 * $z32; + $z3 = $z3 + $z3 + $sat->deep_arg->betao2 * $z33; + $s3 = $cc * $xnoi; + $s2 = -0.5 * $s3 / $sat->deep_arg->betao; + $s4 = $s3 * $sat->deep_arg->betao; + $s1 = -15 * $eq * $s4; + $s5 = $x1 * $x3 + $x2 * $x4; + $s6 = $x2 * $x3 + $x1 * $x4; + $s7 = $x2 * $x4 - $x1 * $x3; + $se = $s1 * $zn * $s5; + $si = $s2 * $zn * ($z11 + $z13); + $sl = -$zn * $s3 * ($z1 + $z3 - 14 - 6 * $sat->deep_arg->eosq); + $sgh = $s4 * $zn * ($z31 + $z33 - 6); + $sh = -$zn * $s2 * ($z21 + $z23); + if ($sat->dps->xqncl < 5.2359877E-2) { + $sh = 0; + } + $sat->dps->ee2 = 2 * $s1 * $s6; + $sat->dps->e3 = 2 * $s1 * $s7; + $sat->dps->xi2 = 2 * $s2 * $z12; + $sat->dps->xi3 = 2 * $s2 * ($z13 - $z11); + $sat->dps->xl2 = -2 * $s3 * $z2; + $sat->dps->xl3 = -2 * $s3 * ($z3 - $z1); + $sat->dps->xl4 = -2 * $s3 * (-21 - 9 * $sat->deep_arg->eosq) * $ze; + $sat->dps->xgh2 = 2 * $s4 * $z32; + $sat->dps->xgh3 = 2 * $s4 * ($z33 - $z31); + $sat->dps->xgh4 = -18 * $s4 * $ze; + $sat->dps->xh2 = -2 * $s2 * $z22; + $sat->dps->xh3 = -2 * $s2 * ($z23 - $z21); + + if ($sat->flags & self::LUNAR_TERMS_DONE_FLAG) { + break; + } + + /* Do lunar terms */ + $sat->dps->sse = $se; + $sat->dps->ssi = $si; + $sat->dps->ssl = $sl; + $sat->dps->ssh = $sh / $sat->deep_arg->sinio; + $sat->dps->ssg = $sgh - $sat->deep_arg->cosio * $sat->dps->ssh; + $sat->dps->se2 = $sat->dps->ee2; + $sat->dps->si2 = $sat->dps->xi2; + $sat->dps->sl2 = $sat->dps->xl2; + $sat->dps->sgh2 = $sat->dps->xgh2; + $sat->dps->sh2 = $sat->dps->xh2; + $sat->dps->se3 = $sat->dps->e3; + $sat->dps->si3 = $sat->dps->xi3; + $sat->dps->sl3 = $sat->dps->xl3; + $sat->dps->sgh3 = $sat->dps->xgh3; + $sat->dps->sh3 = $sat->dps->xh3; + $sat->dps->sl4 = $sat->dps->xl4; + $sat->dps->sgh4 = $sat->dps->xgh4; + $zcosg = $sat->dps->zcosgl; + $zsing = $sat->dps->zsingl; + $zcosi = $sat->dps->zcosil; + $zsini = $sat->dps->zsinil; + $zcosh = $sat->dps->zcoshl * $cosq + $sat->dps->zsinhl * $sinq; + $zsinh = $sinq * $sat->dps->zcoshl - $cosq * $sat->dps->zsinhl; + $zn = Predict::znl; + $cc = Predict::c1l; + $ze = Predict::zel; + $zmo = $sat->dps->zmol; + $sat->flags |= self::LUNAR_TERMS_DONE_FLAG; + } /* End of for(;;) */ + + $sat->dps->sse = $sat->dps->sse + $se; + $sat->dps->ssi = $sat->dps->ssi + $si; + $sat->dps->ssl = $sat->dps->ssl + $sl; + $sat->dps->ssg = $sat->dps->ssg + $sgh - $sat->deep_arg->cosio / $sat->deep_arg->sinio * $sh; + $sat->dps->ssh = $sat->dps->ssh + $sh / $sat->deep_arg->sinio; + + /* Geopotential resonance initialization for 12 hour orbits */ + $sat->flags &= ~self::RESONANCE_FLAG; + $sat->flags &= ~self::SYNCHRONOUS_FLAG; + + if (!(($sat->dps->xnq < 0.0052359877) && ($sat->dps->xnq > 0.0034906585))) { + if( ($sat->dps->xnq < 0.00826) || ($sat->dps->xnq > 0.00924) ) { + return; + } + if ($eq < 0.5) { + return; + } + $sat->flags |= self::RESONANCE_FLAG; + $eoc = $eq * $sat->deep_arg->eosq; + $g201 = -0.306 - ($eq - 0.64) * 0.440; + if ($eq <= 0.65) { + $g211 = 3.616 - 13.247 * $eq + 16.290 * $sat->deep_arg->eosq; + $g310 = -19.302 + 117.390 * $eq - 228.419 * + $sat->deep_arg->eosq + 156.591 * $eoc; + $g322 = -18.9068 + 109.7927 * $eq - 214.6334 * + $sat->deep_arg->eosq + 146.5816 * $eoc; + $g410 = -41.122 + 242.694 * $eq - 471.094 * + $sat->deep_arg->eosq + 313.953 * $eoc; + $g422 = -146.407 + 841.880 * $eq - 1629.014 * + $sat->deep_arg->eosq + 1083.435 * $eoc; + $g520 = -532.114 + 3017.977 * $eq - 5740 * + $sat->deep_arg->eosq + 3708.276 * $eoc; + } else { + $g211 = -72.099 + 331.819 * $eq - 508.738 * + $sat->deep_arg->eosq + 266.724 * $eoc; + $g310 = -346.844 + 1582.851 * $eq - 2415.925 * + $sat->deep_arg->eosq + 1246.113 * $eoc; + $g322 = -342.585 + 1554.908 * $eq - 2366.899 * + $sat->deep_arg->eosq + 1215.972 * $eoc; + $g410 = -1052.797 + 4758.686 * $eq - 7193.992 * + $sat->deep_arg->eosq + 3651.957 * $eoc; + $g422 = -3581.69 + 16178.11 * $eq - 24462.77 * + $sat->deep_arg->eosq+ 12422.52 * $eoc; + if ($eq <= 0.715) { + $g520 = 1464.74 - 4664.75 * $eq + 3763.64 * $sat->deep_arg->eosq; + } else { + $g520 = -5149.66 + 29936.92 * $eq - 54087.36 * + $sat->deep_arg->eosq + 31324.56 * $eoc; + } + } /* End if (eq <= 0.65) */ + + if ($eq < 0.7) { + $g533 = -919.2277 + 4988.61 * $eq - 9064.77 * + $sat->deep_arg->eosq + 5542.21 * $eoc; + $g521 = -822.71072 + 4568.6173 * $eq - 8491.4146 * + $sat->deep_arg->eosq + 5337.524 * $eoc; + $g532 = -853.666 + 4690.25 * $eq - 8624.77 * + $sat->deep_arg->eosq + 5341.4 * $eoc; + } + else { + $g533 = -37995.78 + 161616.52 * $eq - 229838.2* + $sat->deep_arg->eosq + 109377.94 * $eoc; + $g521 = -51752.104 + 218913.95 * $eq - 309468.16* + $sat->deep_arg->eosq + 146349.42 * $eoc; + $g532 = -40023.88 + 170470.89 * $eq - 242699.48* + $sat->deep_arg->eosq + 115605.82 * $eoc; + } /* End if (eq <= 0.7) */ + + $sini2 = $sat->deep_arg->sinio * $sat->deep_arg->sinio; + $f220 = 0.75 * (1 + 2 * $sat->deep_arg->cosio + $sat->deep_arg->theta2); + $f221 = 1.5 * $sini2; + $f321 = 1.875 * $sat->deep_arg->sinio * (1 - 2 * + $sat->deep_arg->cosio - 3 * $sat->deep_arg->theta2); + $f322 = -1.875 * $sat->deep_arg->sinio * (1 + 2* + $sat->deep_arg->cosio - 3 * $sat->deep_arg->theta2); + $f441 = 35 * $sini2 * $f220; + $f442 = 39.3750 * $sini2 * $sini2; + $f522 = 9.84375 * $sat->deep_arg->sinio * ($sini2 * (1 - 2 * $sat->deep_arg->cosio - 5 * + $sat->deep_arg->theta2) + 0.33333333 * (-2 + 4 * $sat->deep_arg->cosio + + 6 * $sat->deep_arg->theta2)); + $f523 = $sat->deep_arg->sinio * (4.92187512 * $sini2 * (-2 - 4 * + $sat->deep_arg->cosio + 10 * $sat->deep_arg->theta2) + 6.56250012 + * (1 + 2 * $sat->deep_arg->cosio - 3 * $sat->deep_arg->theta2)); + $f542 = 29.53125 * $sat->deep_arg->sinio * (2 - 8 * + $sat->deep_arg->cosio + $sat->deep_arg->theta2 * + (-12 + 8 * $sat->deep_arg->cosio + 10 * $sat->deep_arg->theta2)); + $f543 = 29.53125 * $sat->deep_arg->sinio * (-2 - 8 * $sat->deep_arg->cosio + + $sat->deep_arg->theta2 * (12 + 8 * $sat->deep_arg->cosio - 10 * + $sat->deep_arg->theta2)); + $xno2 = $sat->dps->xnq * $sat->dps->xnq; + $ainv2 = $aqnv * $aqnv; + $temp1 = 3 * $xno2 * $ainv2; + $temp = $temp1 * Predict::root22; + $sat->dps->d2201 = $temp * $f220 * $g201; + $sat->dps->d2211 = $temp * $f221 * $g211; + $temp1 = $temp1 * $aqnv; + $temp = $temp1 * Predict::root32; + $sat->dps->d3210 = $temp * $f321 * $g310; + $sat->dps->d3222 = $temp * $f322 * $g322; + $temp1 = $temp1 * $aqnv; + $temp = 2 * $temp1 * Predict::root44; + $sat->dps->d4410 = $temp * $f441 * $g410; + $sat->dps->d4422 = $temp * $f442 * $g422; + $temp1 = $temp1 * $aqnv; + $temp = $temp1 * Predict::root52; + $sat->dps->d5220 = $temp * $f522 * $g520; + $sat->dps->d5232 = $temp * $f523 * $g532; + $temp = 2 * $temp1 * Predict::root54; + $sat->dps->d5421 = $temp * $f542 * $g521; + $sat->dps->d5433 = $temp * $f543 * $g533; + $sat->dps->xlamo = $xmao + $sat->tle->xnodeo + $sat->tle->xnodeo - $sat->dps->thgr - $sat->dps->thgr; + $bfact = $sat->deep_arg->xmdot + $sat->deep_arg->xnodot + + $sat->deep_arg->xnodot - Predict::thdt - Predict::thdt; + $bfact = $bfact + $sat->dps->ssl + $sat->dps->ssh + $sat->dps->ssh; + } else { + $sat->flags |= self::RESONANCE_FLAG; + $sat->flags |= self::SYNCHRONOUS_FLAG; + /* Synchronous resonance terms initialization */ + $g200 = 1 + $sat->deep_arg->eosq * (-2.5 + 0.8125 * $sat->deep_arg->eosq); + $g310 = 1 + 2 * $sat->deep_arg->eosq; + $g300 = 1 + $sat->deep_arg->eosq * (-6 + 6.60937 * $sat->deep_arg->eosq); + $f220 = 0.75 * (1 + $sat->deep_arg->cosio) * (1 + $sat->deep_arg->cosio); + $f311 = 0.9375 * $sat->deep_arg->sinio * $sat->deep_arg->sinio * + (1 + 3 * $sat->deep_arg->cosio) - 0.75 * (1 + $sat->deep_arg->cosio); + $f330 = 1 + $sat->deep_arg->cosio; + $f330 = 1.875 * $f330 * $f330 * $f330; + $sat->dps->del1 = 3 * $sat->dps->xnq * $sat->dps->xnq * $aqnv * $aqnv; + $sat->dps->del2 = 2 * $sat->dps->del1 * $f220 * $g200 * Predict::q22; + $sat->dps->del3 = 3 * $sat->dps->del1 * $f330 * $g300 * Predict::q33 * $aqnv; + $sat->dps->del1 = $sat->dps->del1 * $f311 * $g310 * Predict::q31 * $aqnv; + $sat->dps->fasx2 = 0.13130908; + $sat->dps->fasx4 = 2.8843198; + $sat->dps->fasx6 = 0.37448087; + $sat->dps->xlamo = $xmao + $sat->tle->xnodeo + $sat->tle->omegao - $sat->dps->thgr; + $bfact = $sat->deep_arg->xmdot + $xpidot - Predict::thdt; + $bfact = $bfact + $sat->dps->ssl + $sat->dps->ssg + $sat->dps->ssh; + } /* End if( !(xnq < 0.0052359877) && (xnq > 0.0034906585) ) */ + + $sat->dps->xfact = $bfact - $sat->dps->xnq; + + /* Initialize integrator */ + $sat->dps->xli = $sat->dps->xlamo; + $sat->dps->xni = $sat->dps->xnq; + $sat->dps->atime = 0; + $sat->dps->stepp = 720; + $sat->dps->stepn = -720; + $sat->dps->step2 = 259200; + /* End case self::dpinit: */ + return; + + case self::dpsec: /* Entrance for deep space secular effects */ + $sat->deep_arg->xll = $sat->deep_arg->xll + $sat->dps->ssl * $sat->deep_arg->t; + $sat->deep_arg->omgadf = $sat->deep_arg->omgadf + $sat->dps->ssg * $sat->deep_arg->t; + $sat->deep_arg->xnode = $sat->deep_arg->xnode + $sat->dps->ssh * $sat->deep_arg->t; + $sat->deep_arg->em = $sat->tle->eo + $sat->dps->sse * $sat->deep_arg->t; + $sat->deep_arg->xinc = $sat->tle->xincl + $sat->dps->ssi * $sat->deep_arg->t; + if ($sat->deep_arg->xinc < 0) { + $sat->deep_arg->xinc = -$sat->deep_arg->xinc; + $sat->deep_arg->xnode = $sat->deep_arg->xnode + Predict::pi; + $sat->deep_arg->omgadf = $sat->deep_arg->omgadf - Predict::pi; + } + if(~$sat->flags & self::RESONANCE_FLAG ) { + return; + } + + do { + if ( ($sat->dps->atime == 0) || + (($sat->deep_arg->t >= 0) && ($sat->dps->atime < 0)) || + (($sat->deep_arg->t < 0) && ($sat->dps->atime >= 0)) ) { + /* Epoch restart */ + if ($sat->deep_arg->t >= 0) { + $delt = $sat->dps->stepp; + } else { + $delt = $sat->dps->stepn; + } + + $sat->dps->atime = 0; + $sat->dps->xni = $sat->dps->xnq; + $sat->dps->xli = $sat->dps->xlamo; + } else { + if (abs($sat->deep_arg->t) >= abs($sat->dps->atime)) { + if ($sat->deep_arg->t > 0) { + $delt = $sat->dps->stepp; + } else { + $delt = $sat->dps->stepn; + } + } + } + + do { + if (abs($sat->deep_arg->t - $sat->dps->atime) >= $sat->dps->stepp) { + $sat->flags |= self::DO_LOOP_FLAG; + $sat->flags &= ~self::EPOCH_RESTART_FLAG; + } + else { + $ft = $sat->deep_arg->t - $sat->dps->atime; + $sat->flags &= ~self::DO_LOOP_FLAG; + } + + if (abs($sat->deep_arg->t) < abs($sat->dps->atime)) { + if ($sat->deep_arg->t >= 0) { + $delt = $sat->dps->stepn; + } else { + $delt = $sat->dps->stepp; + } + $sat->flags |= (self::DO_LOOP_FLAG | self::EPOCH_RESTART_FLAG); + } + + /* Dot terms calculated */ + if ($sat->flags & self::SYNCHRONOUS_FLAG) { + $xndot = $sat->dps->del1 * sin($sat->dps->xli - $sat->dps->fasx2) + $sat->dps->del2 * sin(2 * ($sat->dps->xli - $sat->dps->fasx4)) + + $sat->dps->del3 * sin(3 * ($sat->dps->xli - $sat->dps->fasx6)); + $xnddt = $sat->dps->del1 * cos($sat->dps->xli - $sat->dps->fasx2) + 2 * $sat->dps->del2 * cos(2 * ($sat->dps->xli - $sat->dps->fasx4)) + + 3 * $sat->dps->del3 * cos(3 * ($sat->dps->xli - $sat->dps->fasx6)); + } else { + $xomi = $sat->dps->omegaq + $sat->deep_arg->omgdot * $sat->dps->atime; + $x2omi = $xomi + $xomi; + $x2li = $sat->dps->xli + $sat->dps->xli; + $xndot = $sat->dps->d2201 * sin($x2omi + $sat->dps->xli - Predict::g22) + + $sat->dps->d2211 * sin($sat->dps->xli - Predict::g22) + + $sat->dps->d3210 * sin($xomi + $sat->dps->xli - Predict::g32) + + $sat->dps->d3222 * sin(-$xomi + $sat->dps->xli - Predict::g32) + + $sat->dps->d4410 * sin($x2omi + $x2li- Predict::g44) + + $sat->dps->d4422 * sin($x2li- Predict::g44) + + $sat->dps->d5220 * sin($xomi + $sat->dps->xli- Predict::g52) + + $sat->dps->d5232 * sin(-$xomi + $sat->dps->xli- Predict::g52) + + $sat->dps->d5421 * sin($xomi + $x2li - Predict::g54) + + $sat->dps->d5433 * sin(-$xomi + $x2li - Predict::g54); + $xnddt = $sat->dps->d2201 * cos($x2omi + $sat->dps->xli- Predict::g22) + + $sat->dps->d2211 * cos($sat->dps->xli - Predict::g22) + + $sat->dps->d3210 * cos($xomi + $sat->dps->xli - Predict::g32) + + $sat->dps->d3222 * cos(-$xomi + $sat->dps->xli - Predict::g32) + + $sat->dps->d5220 * cos($xomi + $sat->dps->xli - Predict::g52) + + $sat->dps->d5232 * cos(-$xomi + $sat->dps->xli - Predict::g52) + + 2 * ($sat->dps->d4410 * cos($x2omi + $x2li - Predict::g44) + + $sat->dps->d4422 * cos($x2li - Predict::g44) + + $sat->dps->d5421 * cos($xomi + $x2li - Predict::g54) + + $sat->dps->d5433 * cos(-$xomi + $x2li - Predict::g54)); + } /* End of if (isFlagSet(SYNCHRONOUS_FLAG)) */ + + $xldot = $sat->dps->xni + $sat->dps->xfact; + $xnddt = $xnddt * $xldot; + + if ($sat->flags & self::DO_LOOP_FLAG) { + $sat->dps->xli = $sat->dps->xli + $xldot * $delt + $xndot * $sat->dps->step2; + $sat->dps->xni = $sat->dps->xni + $xndot * $delt + $xnddt * $sat->dps->step2; + $sat->dps->atime = $sat->dps->atime + $delt; + } + } while (($sat->flags & self::DO_LOOP_FLAG) && + (~$sat->flags & self::EPOCH_RESTART_FLAG)); + } + while (($sat->flags & self::DO_LOOP_FLAG) && ($sat->flags & self::EPOCH_RESTART_FLAG)); + + $sat->deep_arg->xn = $sat->dps->xni + $xndot * $ft + $xnddt * $ft * $ft * 0.5; + $xl = $sat->dps->xli + $xldot * $ft + $xndot * $ft * $ft * 0.5; + $temp = -$sat->deep_arg->xnode + $sat->dps->thgr + $sat->deep_arg->t * Predict::thdt; + + if (~$sat->flags & self::SYNCHRONOUS_FLAG) { + $sat->deep_arg->xll = $xl + $temp + $temp; + } else { + $sat->deep_arg->xll = $xl - $sat->deep_arg->omgadf + $temp; + } + + return; + /* End case dpsec: */ + + case self::dpper: /* Entrance for lunar-solar periodics */ + $sinis = sin($sat->deep_arg->xinc); + $cosis = cos($sat->deep_arg->xinc); + if (abs($sat->dps->savtsn - $sat->deep_arg->t) >= 30) { + $sat->dps->savtsn = $sat->deep_arg->t; + $zm = $sat->dps->zmos + Predict::zns * $sat->deep_arg->t; + $zf = $zm + 2 * Predict::zes * sin($zm); + $sinzf = sin($zf); + $f2 = 0.5 * $sinzf * $sinzf - 0.25; + $f3 = -0.5 * $sinzf * cos($zf); + $ses = $sat->dps->se2 * $f2 + $sat->dps->se3 * $f3; + $sis = $sat->dps->si2 * $f2 + $sat->dps->si3 * $f3; + $sls = $sat->dps->sl2 * $f2 + $sat->dps->sl3 * $f3 + $sat->dps->sl4 * $sinzf; + $sat->dps->sghs = $sat->dps->sgh2 * $f2 + $sat->dps->sgh3 * $f3 + $sat->dps->sgh4 * $sinzf; + $sat->dps->shs = $sat->dps->sh2 * $f2 + $sat->dps->sh3 * $f3; + $zm = $sat->dps->zmol + Predict::znl * $sat->deep_arg->t; + $zf = $zm + 2 * Predict::zel * sin($zm); + $sinzf = sin($zf); + $f2 = 0.5 * $sinzf * $sinzf - 0.25; + $f3 = -0.5 * $sinzf * cos($zf); + $sel = $sat->dps->ee2 * $f2 + $sat->dps->e3 * $f3; + $sil = $sat->dps->xi2 * $f2 + $sat->dps->xi3 * $f3; + $sll = $sat->dps->xl2 * $f2 + $sat->dps->xl3 * $f3 + $sat->dps->xl4 * $sinzf; + $sat->dps->sghl = $sat->dps->xgh2 * $f2 + $sat->dps->xgh3 * $f3 + $sat->dps->xgh4 * $sinzf; + $sat->dps->sh1 = $sat->dps->xh2 * $f2 + $sat->dps->xh3 * $f3; + $sat->dps->pe = $ses + $sel; + $sat->dps->pinc = $sis + $sil; + $sat->dps->pl = $sls + $sll; + } + + $pgh = $sat->dps->sghs + $sat->dps->sghl; + $ph = $sat->dps->shs + $sat->dps->sh1; + $sat->deep_arg->xinc = $sat->deep_arg->xinc + $sat->dps->pinc; + $sat->deep_arg->em = $sat->deep_arg->em + $sat->dps->pe; + + if ($sat->dps->xqncl >= 0.2) { + /* Apply periodics directly */ + $ph = $ph / $sat->deep_arg->sinio; + $pgh = $pgh - $sat->deep_arg->cosio * $ph; + $sat->deep_arg->omgadf = $sat->deep_arg->omgadf + $pgh; + $sat->deep_arg->xnode = $sat->deep_arg->xnode + $ph; + $sat->deep_arg->xll = $sat->deep_arg->xll + $sat->dps->pl; + } else { + /* Apply periodics with Lyddane modification */ + $sinok = sin($sat->deep_arg->xnode); + $cosok = cos($sat->deep_arg->xnode); + $alfdp = $sinis * $sinok; + $betdp = $sinis * $cosok; + $dalf = $ph * $cosok + $sat->dps->pinc * $cosis * $sinok; + $dbet = -$ph * $sinok + $sat->dps->pinc * $cosis * $cosok; + $alfdp = $alfdp + $dalf; + $betdp = $betdp + $dbet; + $sat->deep_arg->xnode = Predict_Math::FMod2p($sat->deep_arg->xnode); + $xls = $sat->deep_arg->xll + $sat->deep_arg->omgadf + $cosis * $sat->deep_arg->xnode; + $dls = $sat->dps->pl + $pgh - $sat->dps->pinc * $sat->deep_arg->xnode * $sinis; + $xls = $xls + $dls; + $xnoh = $sat->deep_arg->xnode; + $sat->deep_arg->xnode = Predict_Math::AcTan($alfdp, $betdp); + + /* This is a patch to Lyddane modification */ + /* suggested by Rob Matson. */ + if(abs($xnoh - $sat->deep_arg->xnode) > Predict::pi) { + if ($sat->deep_arg->xnode < $xnoh) { + $sat->deep_arg->xnode += Predict::twopi; + } else { + $sat->deep_arg->xnode -= Predict::twopi; + } + } + + $sat->deep_arg->xll = $sat->deep_arg->xll + $sat->dps->pl; + $sat->deep_arg->omgadf = $xls - $sat->deep_arg->xll - cos($sat->deep_arg->xinc) * + $sat->deep_arg->xnode; + } /* End case dpper: */ + return; + + } /* End switch(ientry) */ + + } /* End of Deep() */ + + /** + * Singleton + * + * @param Predict_Sat $sat The current satellite data instance + * + * @return Predict_SGPSDP + */ + public static function getInstance(Predict_Sat $sat) + { + static $instances = array(); + $catnr = $sat->tle->catnr; + if (!isset($instances[$catnr])) { + $instances[$catnr] = new self(); + } + return $instances[$catnr]; + } +} +?> diff --git a/predict/Predict/SGSDPStatic.php b/predict/Predict/SGSDPStatic.php new file mode 100644 index 000000000..12a43a4cd --- /dev/null +++ b/predict/Predict/SGSDPStatic.php @@ -0,0 +1,39 @@ +header); + $this->name = $headerParts[0]; + $this->nickname = $this->name; + $this->tle = $tle; + $this->pos = new Predict_Vector(); + $this->vel = new Predict_Vector(); + $this->sgps = new Predict_SGSDPStatic(); + $this->deep_arg = new Predict_DeepArg(); + $this->dps = new Predict_DeepStatic(); + + $this->select_ephemeris(); + $this->sat_data_init_sat($this); + } + + /* Selects the apropriate ephemeris type to be used */ + /* for predictions according to the data in the TLE */ + /* It also processes values in the tle set so that */ + /* they are apropriate for the sgp4/sdp4 routines */ + public function select_ephemeris() + { + /* Preprocess tle set */ + $this->tle->xnodeo *= Predict::de2ra; + $this->tle->omegao *= Predict::de2ra; + $this->tle->xmo *= Predict::de2ra; + $this->tle->xincl *= Predict::de2ra; + $temp = Predict::twopi / Predict::xmnpda / Predict::xmnpda; + + /* store mean motion before conversion */ + $this->meanmo = $this->tle->xno; + $this->tle->xno = $this->tle->xno * $temp * Predict::xmnpda; + $this->tle->xndt2o *= $temp; + $this->tle->xndd6o = $this->tle->xndd6o * $temp / Predict::xmnpda; + $this->tle->bstar /= Predict::ae; + + /* Period > 225 minutes is deep space */ + $dd1 = Predict::xke / $this->tle->xno; + $dd2 = Predict::tothrd; + $a1 = pow($dd1, $dd2); + $r1 = cos($this->tle->xincl); + $dd1 = 1.0 - $this->tle->eo * $this->tle->eo; + $temp = Predict::ck2 * 1.5 * ($r1 * $r1 * 3.0 - 1.0) / pow($dd1, 1.5); + $del1 = $temp / ($a1 * $a1); + $ao = $a1 * (1.0 - $del1 * (Predict::tothrd * 0.5 + $del1 * + ($del1 * 1.654320987654321 + 1.0))); + $delo = $temp / ($ao * $ao); + $xnodp = $this->tle->xno / ($delo + 1.0); + + /* Select a deep-space/near-earth ephemeris */ + if (Predict::twopi / $xnodp / Predict::xmnpda >= .15625) { + $this->flags |= Predict_SGPSDP::DEEP_SPACE_EPHEM_FLAG; + } else { + $this->flags &= ~Predict_SGPSDP::DEEP_SPACE_EPHEM_FLAG; + } + } + + /** Initialise satellite data. + * @param sat The satellite to initialise. + * @param qth Optional QTH info, use (0,0) if NULL. + * + * This function calculates the satellite data at t = 0, ie. epoch time + * The function is called automatically by gtk_sat_data_read_sat. + */ + public function sat_data_init_sat(Predict_Sat $sat, Predict_QTH $qth = null) + { + $obs_geodetic = new Predict_Geodetic(); + $obs_set = new Predict_ObsSet(); + $sat_geodetic = new Predict_Geodetic(); + /* double jul_utc, age; */ + + $jul_utc = Predict_Time::Julian_Date_of_Epoch($sat->tle->epoch); // => tsince = 0.0 + $sat->jul_epoch = $jul_utc; + + /* initialise observer location */ + if ($qth != null) { + $obs_geodetic->lon = $qth->lon * Predict::de2ra; + $obs_geodetic->lat = $qth->lat * Predict::de2ra; + $obs_geodetic->alt = $qth->alt / 1000.0; + $obs_geodetic->theta = 0; + } + else { + $obs_geodetic->lon = 0.0; + $obs_geodetic->lat = 0.0; + $obs_geodetic->alt = 0.0; + $obs_geodetic->theta = 0; + } + + /* execute computations */ + $sdpsgp = Predict_SGPSDP::getInstance($sat); + if ($sat->flags & Predict_SGPSDP::DEEP_SPACE_EPHEM_FLAG) { + $sdpsgp->SDP4($sat, 0.0); + } else { + $sdpsgp->SGP4($sat, 0.0); + } + + /* scale position and velocity to km and km/sec */ + Predict_Math::Convert_Sat_State($sat->pos, $sat->vel); + + /* get the velocity of the satellite */ + $sat->vel->w = sqrt($sat->vel->x * $sat->vel->x + $sat->vel->y * $sat->vel->y + $sat->vel->z * $sat->vel->z); + $sat->velo = $sat->vel->w; + Predict_SGPObs::Calculate_Obs($jul_utc, $sat->pos, $sat->vel, $obs_geodetic, $obs_set); + Predict_SGPObs::Calculate_LatLonAlt($jul_utc, $sat->pos, $sat_geodetic); + + while ($sat_geodetic->lon < -Predict::pi) { + $sat_geodetic->lon += Predict::twopi; + } + + while ($sat_geodetic->lon > Predict::pi) { + $sat_geodetic->lon -= Predict::twopi; + } + + $sat->az = Predict_Math::Degrees($obs_set->az); + $sat->el = Predict_Math::Degrees($obs_set->el); + $sat->range = $obs_set->range; + $sat->range_rate = $obs_set->range_rate; + $sat->ssplat = Predict_Math::Degrees($sat_geodetic->lat); + $sat->ssplon = Predict_Math::Degrees($sat_geodetic->lon); + $sat->alt = $sat_geodetic->alt; + $sat->ma = Predict_Math::Degrees($sat->phase); + $sat->ma *= 256.0 / 360.0; + $sat->footprint = 2.0 * Predict::xkmper * acos (Predict::xkmper/$sat->pos->w); + $age = 0.0; + $sat->orbit = floor(($sat->tle->xno * Predict::xmnpda / Predict::twopi + + $age * $sat->tle->bstar * Predict::ae) * $age + + $sat->tle->xmo / Predict::twopi) + $sat->tle->revnum - 1; + + /* orbit type */ + $sat->otype = $sat->get_orbit_type($sat); + } + + public function get_orbit_type(Predict_Sat $sat) + { + $orbit = Predict_SGPSDP::ORBIT_TYPE_UNKNOWN; + + if ($this->geostationary($sat)) { + $orbit = Predict_SGPSDP::ORBIT_TYPE_GEO; + } else if ($this->decayed($sat)) { + $orbit = Predict_SGPSDP::ORBIT_TYPE_DECAYED; + } else { + $orbit = Predict_SGPSDP::ORBIT_TYPE_UNKNOWN; + } + + return $orbit; + } + + + /** Determinte whether satellite is in geostationary orbit. + * @author John A. Magliacane, KD2BD + * @param sat Pointer to satellite data. + * @return TRUE if the satellite appears to be in geostationary orbit, + * FALSE otherwise. + * + * A satellite is in geostationary orbit if + * + * fabs (sat.meanmotion - 1.0027) < 0.0002 + * + * Note: Appearantly, the mean motion can deviate much more from 1.0027 than 0.0002 + */ + public function geostationary(Predict_Sat $sat) + { + if (abs($sat->meanmo - 1.0027) < 0.0002) { + return true; + } else { + return false; + } + } + + + /** Determine whether satellite has decayed. + * @author John A. Magliacane, KD2BD + * @author Alexandru Csete, OZ9AEC + * @param sat Pointer to satellite data. + * @return TRUE if the satellite appears to have decayed, FALSE otherwise. + * @bug Modified version of the predict code but it is not tested. + * + * A satellite is decayed if + * + * satepoch + ((16.666666 - sat.meanmo) / (10.0*fabs(sat.drag))) < "now" + * + */ + public function decayed(Predict_Sat $sat) + { + /* tle.xndt2o/(twopi/xmnpda/xmnpda) is the value before converted the + value matches up with the value in predict 2.2.3 */ + /*** FIXME decayed is treated as a static quantity. + It is time dependent. Also sat->jul_utc is often zero + when this function is called + ***/ + if ($sat->jul_epoch + ((16.666666 - $sat->meanmo) / + (10.0 * abs($sat->tle->xndt2o / (Predict::twopi / Predict::xmnpda / Predict::xmnpda)))) < $sat->jul_utc) { + return true; + } else { + return false; + } + } + + /** + * Experimental attempt at calculating apparent magnitude. Known intrinsic + * magnitudes are listed inside the function for now. + * + * @param float $time The daynum the satellite is calculated for + * @param Predict_QTH $qth The observer location + * + * @return null on failure, float otherwise + */ + public function calculateApparentMagnitude($time, Predict_QTH $qth) + { + // Recorded intrinsic magnitudes and their respective + // illumination and distance from heavens-above.com + static $intrinsicMagnitudes = array( + '25544' => array( + 'mag' => -1.3, + 'illum' => .5, + 'distance' => 1000, + ) + ); + + // Return null if we don't have a record of the intrinsic mag + if (!isset($intrinsicMagnitudes[$this->tle->catnr])) { + return null; + } + $imag = $intrinsicMagnitudes[$this->tle->catnr]; + + // Convert the observer's geodetic info to radians and km so + // we can compare vectors + $observerGeo = new Predict_Geodetic(); + $observerGeo->lat = Predict_Math::Radians($qth->lat); + $observerGeo->lon = Predict_Math::Radians($qth->lon); + $observerGeo->alt = $qth->alt * 1000; + + // Now determine the sun and observer positions + $observerPos = new Predict_Vector(); + $observerVel = new Predict_Vector(); + $solarVector = new Predict_Vector(); + Predict_Solar::Calculate_Solar_Position($time, $solarVector); + Predict_SGPObs::Calculate_User_PosVel($time, $observerGeo, $observerPos, $observerVel); + + // Determine the solar phase and and thus the percent illumination + $observerSatPos = new Predict_Vector(); + Predict_Math::Vec_Sub($this->pos, $observerPos, $observerSatPos); + $phaseAngle = Predict_Math::Degrees(Predict_Math::Angle($solarVector, $observerSatPos)); + $illum = $phaseAngle / 180; + + $illuminationChange = $illum / $imag['illum']; + $inverseSquareOfDistanceChange = pow(($imag['distance'] / $this->range), 2); + $changeInMagnitude = log( + $illuminationChange * $inverseSquareOfDistanceChange, + self::POGSONS_RATIO + ); + + return $imag['mag'] - $changeInMagnitude; + } +} diff --git a/predict/Predict/Solar.php b/predict/Predict/Solar.php new file mode 100644 index 000000000..8ec351b0f --- /dev/null +++ b/predict/Predict/Solar.php @@ -0,0 +1,115 @@ +x = $R * cos($Lsa); + $solar_vector->y = $R * sin($Lsa) * cos($eps); + $solar_vector->z = $R * sin($Lsa) * sin($eps); + $solar_vector->w = $R; + } + + /* Calculates stellite's eclipse status and depth */ + public static function Sat_Eclipsed(Predict_Vector $pos, Predict_Vector $sol, &$depth) + { + $Rho = new Predict_Vector(); + $earth = new Predict_Vector(); + + /* Determine partial eclipse */ + $sd_earth = Predict_Math::ArcSin(Predict::xkmper / $pos->w); + Predict_Math::Vec_Sub($sol, $pos, $Rho); + $sd_sun = Predict_Math::ArcSin(Predict::__sr__ / $Rho->w); + Predict_Math::Scalar_Multiply(-1, $pos, $earth); + $delta = Predict_Math::Angle($sol, $earth); + $depth = $sd_earth - $sd_sun - $delta; + + if ($sd_earth < $sd_sun) { + return 0; + } else if ($depth >= 0) { + return 1; + } else { + return 0; + } + } + + /** + * Finds the current location of the sun based on the observer location + * + * @param Predict_QTH $qth The observer location + * @param int $daynum The daynum or null to use the current daynum + * + * @return Predict_ObsSet + */ + public static function FindSun(Predict_QTH $qth, $daynum = null) + { + if ($daynum === null) { + $daynum = Predict_Time::get_current_daynum(); + } + + $obs_geodetic = new Predict_Geodetic(); + $obs_geodetic->lon = $qth->lon * Predict::de2ra; + $obs_geodetic->lat = $qth->lat * Predict::de2ra; + $obs_geodetic->alt = $qth->alt / 1000.0; + $obs_geodetic->theta = 0; + + $solar_vector = new Predict_Vector(); + $zero_vector = new Predict_Vector(); + $solar_set = new Predict_ObsSet(); + + self::Calculate_Solar_Position($daynum, $solar_vector); + Predict_SGPObs::Calculate_Obs( + $daynum, + $solar_vector, + $zero_vector, + $obs_geodetic, + $solar_set + ); + + $solar_set->az = Predict_Math::Degrees($solar_set->az); + $solar_set->el = Predict_Math::Degrees($solar_set->el); + + return $solar_set; + } +} diff --git a/predict/Predict/TLE.php b/predict/Predict/TLE.php new file mode 100644 index 000000000..e78e4345c --- /dev/null +++ b/predict/Predict/TLE.php @@ -0,0 +1,232 @@ +Good_Elements($line1, $line2)) { + throw new Predict_Exception('Invalid TLE contents'); + } + + $this->header = $header; + $this->line1 = $line1; + $this->line2 = $line2; + + /** Decode Card 1 **/ + /* Satellite's catalogue number */ + $this->catnr = (int) substr($line1, 2, 5); + + /* International Designator for satellite */ + $this->idesg = substr($line1, 9, 8); + + /* Epoch time; this is the complete, unconverted epoch. */ + /* Replace spaces with 0 before casting, as leading spaces are allowed */ + $this->epoch = (float) str_replace(' ', '0', substr($line1, 18, 14)); + + /* Now, convert the epoch time into year, day + and fraction of day, according to: + + YYDDD.FFFFFFFF + */ + + // Adjust for 2 digit year through 2056 + $this->epoch_year = (int) substr($line1, 18, 2); + if ($this->epoch_year > 56) { + $this->epoch_year = $this->epoch_year + 1900; + } else { + $this->epoch_year = $this->epoch_year + 2000; + } + + /* Epoch day */ + $this->epoch_day = (int) substr($line1, 20, 3); + + /* Epoch fraction of day */ + $this->epoch_fod = (float) substr($line1, 23, 9); + + + /* Satellite's First Time Derivative */ + $this->xndt2o = (float) substr($line1, 33, 10); + + /* Satellite's Second Time Derivative */ + $this->xndd6o = (float) (substr($line1, 44, 1) . '.' . substr($line1, 45, 5) . 'E' . substr($line1, 50, 2)); + + /* Satellite's bstar drag term + FIXME: How about buff[0] ???? + */ + $this->bstar = (float) (substr($line1, 53, 1) . '.' . substr($line1, 54, 5) . 'E' . substr($line1, 59, 2)); + + /* Element Number */ + $this->elset = (int) substr($line1, 64, 4); + + /** Decode Card 2 **/ + /* Satellite's Orbital Inclination (degrees) */ + $this->xincl = (float) substr($line2, 8, 8); + + /* Satellite's RAAN (degrees) */ + $this->xnodeo = (float) substr($line2, 17, 8); + + /* Satellite's Orbital Eccentricity */ + $this->eo = (float) ('.' . substr($line2, 26, 7)); + + /* Satellite's Argument of Perigee (degrees) */ + $this->omegao = (float) substr($line2, 34, 8); + + /* Satellite's Mean Anomaly of Orbit (degrees) */ + $this->xmo = (float) substr($line2, 43, 8); + + /* Satellite's Mean Motion (rev/day) */ + $this->xno = (float) substr($line2, 52, 11); + + /* Satellite's Revolution number at epoch */ + $this->revnum = (float) substr($line2, 63, 5); + } + + /* Calculates the checksum mod 10 of a line from a TLE set and */ + /* returns true if it compares with checksum in column 68, else false.*/ + /* tle_set is a character string holding the two lines read */ + /* from a text file containing NASA format Keplerian elements. */ + /* NOTE!!! The stuff about two lines is not quite true. + The function assumes that tle_set[0] is the begining + of the line and that there are 68 elements - see the consumer + */ + public function Checksum_Good($tle_set) + { + if (strlen($tle_set) < 69) { + return false; + } + + $checksum = 0; + + for ($i = 0; $i < 68; $i++) { + if (($tle_set[$i] >= '0') && ($tle_set[$i] <= '9')) { + $value = $tle_set[$i] - '0'; + } else if ($tle_set[$i] == '-' ) { + $value = 1; + } else { + $value = 0; + } + + $checksum += $value; + } + + $checksum %= 10; + $check_digit = $tle_set[68] - '0'; + + return $checksum == $check_digit; + } + + /* Carries out various checks on a TLE set to verify its validity */ + /* $line1 is the first line of the TLE, $line2 is the second line */ + /* from a text file containing NASA format Keplerian elements. */ + public function Good_Elements($line1, $line2) + { + /* Verify checksum of both lines of a TLE set */ + if (!$this->Checksum_Good($line1) || !$this->Checksum_Good($line2)) { + return false; + } + + /* Check the line number of each line */ + if (($line1[0] != '1') || ($line2[0] != '2')) { + return false; + } + + /* Verify that Satellite Number is same in both lines */ + if (strncmp($line1[2], $line2[2], 5) != 0) { + return false; + } + + /* Check that various elements are in the right place */ + if (($line1[23] != '.') || + ($line1[34] != '.') || + ($line2[11] != '.') || + ($line2[20] != '.') || + ($line2[37] != '.') || + ($line2[46] != '.') || + ($line2[54] != '.') || + (strncmp(substr($line1, 61), ' 0 ', 3) != 0)) { + + return false; + } + + return true; + } + + /** + * A function to allow checksum creation of a line. This is driven by + * the fact that some TLEs from SpaceTrack are missing checksum numbers. + * You can use this to create a checksum for a line, but you should + * probably have confidence that the TLE data itself is good. YMMV. + * + * @throws Predict_Exception if the line is not exactly 68 chars + * @return string + */ + static public function createChecksum($line) + { + if (strlen($line) != 68) { + throw Predict_Exception('Invalid line, needs to e 68 chars'); + } + + $checksum = 0; + + for ($i = 0; $i < 68; $i++) { + if (($line[$i] >= '0') && ($line[$i] <= '9')) { + $value = (int) $line[$i]; + } else if ($line[$i] == '-' ) { + $value = 1; + } else { + $value = 0; + } + + $checksum += $value; + } + + $checksum %= 10; + + return $checksum; + } +} diff --git a/predict/Predict/Time.php b/predict/Predict/Time.php new file mode 100644 index 000000000..bfb515850 --- /dev/null +++ b/predict/Predict/Time.php @@ -0,0 +1,222 @@ +ds50 = $jd - 2433281.5 + $UT; + + return Predict_Math::FMod2p(6.3003880987 * $deep_arg->ds50 + 1.72944494); + } + + /* See the ThetaG doc block above */ + public static function ThetaG_JD($jd) + { + /* Reference: The 1992 Astronomical Almanac, page B6. */ + $UT = Predict_Math::Frac($jd + 0.5); + $jd = $jd - $UT; + $TU = ($jd - 2451545.0) / 36525; + $GMST = 24110.54841 + $TU * (8640184.812866 + $TU * (0.093104 - $TU * 6.2E-6)); + $GMST = Predict_Math::Modulus($GMST + Predict::secday * Predict::omega_E * $UT, Predict::secday); + + return Predict::twopi * $GMST / Predict::secday; + } + + /** + * Read the system clock and return the current Julian day. From phpPredict + * + * @return float + */ + public static function get_current_daynum() { + // Gets the current decimal day number from microtime + + list($usec, $sec) = explode(' ', microtime()); + return self::unix2daynum($sec, $usec); + } + + /** + * Converts a standard unix timestamp and optional + * milliseconds to a daynum + * + * @param int $sec Seconds from the unix epoch + * @param int $usec Optional milliseconds + * + * @return float + */ + public static function unix2daynum($sec, $usec = 0) + { + $time = ((($sec + $usec) / 86400.0) - 3651.0); + return $time + 2444238.5; + } + + /* The function Delta_ET has been added to allow calculations on */ + /* the position of the sun. It provides the difference between UT */ + /* (approximately the same as UTC) and ET (now referred to as TDT).*/ + /* This function is based on a least squares fit of data from 1950 */ + /* to 1991 and will need to be updated periodically. */ + public static function Delta_ET($year) + { + /* Values determined using data from 1950-1991 in the 1990 + Astronomical Almanac. See DELTA_ET.WQ1 for details. */ + + $delta_et = 26.465 + 0.747622 * ($year - 1950) + + 1.886913 * sin(Predict::twopi * ($year - 1975) / 33); + + return $delta_et; + } + + /** + * Converts a daynum to a unix timestamp. From phpPredict. + * + * @param float $dn Julian Daynum + * + * @return float + */ + public static function daynum2unix($dn) { + // Converts a daynum to a UNIX timestamp + + return (86400.0 * ($dn - 2444238.5 + 3651.0)); + } + + /** + * Converts a daynum to a readable time format. + * + * @param float $dn The julian date + * @param string $zone The zone string, defaults to America/Los_Angeles + * @param string $format The date() function's format string. Defaults to m-d-Y H:i:s + * + * @return string + */ + public static function daynum2readable($dn, $zone = 'America/Los_Angeles', $format = 'm-d-Y H:i:s') + { + $unix = self::daynum2unix($dn); + $date = new DateTime("@" . round($unix)); + $dateTimezone = new DateTimezone($zone); + $date->setTimezone($dateTimezone); + return $date->format($format); + } + + /** + * Returns the unix timestamp of a TLE's epoch + * + * @param Predict_TLE $tle The TLE object + * + * @return int + */ + public static function getEpochTimeStamp(Predict_TLE $tle) + { + $year = $tle->epoch_year; + $day = $tle->epoch_day; + $sec = round(86400 * $tle->epoch_fod); + + $zone = new DateTimeZone('GMT'); + $date = new DateTime(); + $date->setTimezone($zone); + $date->setDate($year, 1, 1); + $date->setTime(0, 0, 0); + + return $date->format('U') + (86400 * $day) + $sec - 86400; + } +} diff --git a/predict/Predict/Vector.php b/predict/Predict/Vector.php new file mode 100644 index 000000000..84749058b --- /dev/null +++ b/predict/Predict/Vector.php @@ -0,0 +1,13 @@ +setOptions(array( + 'baseinstalldir' => '/', + 'simpleoutput' => true, + 'packagedirectory' => './', + 'filelistgenerator' => 'file', + 'ignore' => array('generatePackage.php', 'xhprof_lib/*'), + 'dir_roles' => array( + 'tests' => 'test', + 'examples' => 'doc' + ), + 'exceptions' => array('README.md' => 'doc'), +)); + +$packagexml->setPackage('Predict'); +$packagexml->setSummary('A partial port of the Gpredict program for satellite tracking'); +$packagexml->setDescription( + 'Predict is a partial PHP port of the Gpredict (http://gpredict.oz9aec.net/) program that ' + . 'allows real-time tracking and orbit prediction of satellites from two line element sets. ' + . 'It supports the SGP4 and SDP4 models for prediction.' +); + +$packagexml->setChannel('shupp.github.com/pirum'); +$packagexml->setAPIVersion('0.2.2'); +$packagexml->setReleaseVersion('0.2.2'); + +$packagexml->setReleaseStability('alpha'); + +$packagexml->setAPIStability('alpha'); + +$packagexml->setNotes(' +* Addec Predict_TLE::createChecksum() +* Updates to examples +'); +$packagexml->setPackageType('php'); +$packagexml->addRelease(); + +$packagexml->detectDependencies(); + +$packagexml->addMaintainer('lead', + 'shupp', + 'Bill Shupp', + 'shupp@php.net'); +$packagexml->setLicense('GPL v2.1', + 'http://www.opensource.org/licenses/gpl-license.php'); + +$packagexml->setPhpDep('5.2.0'); +$packagexml->setPearinstallerDep('1.4.0b1'); +$packagexml->addExtensionDep('required', 'date'); + +$packagexml->generateContents(); +$packagexml->writePackageFile(); + +?> diff --git a/predict/package.xml b/predict/package.xml new file mode 100644 index 000000000..783eb5477 --- /dev/null +++ b/predict/package.xml @@ -0,0 +1,182 @@ + + + Predict + shupp.github.com/pirum + A partial port of the Gpredict program for satellite tracking + Predict is a partial PHP port of the Gpredict (http://gpredict.oz9aec.net/) program that allows real-time tracking and orbit prediction of satellites from two line element sets. It supports the SGP4 and SDP4 models for prediction. + + Bill Shupp + shupp + shupp@php.net + yes + + 2014-03-06 + + + 0.2.2 + 0.2.2 + + + alpha + alpha + + GPL v2.1 + +* Addec Predict_TLE::createChecksum() +* Updates to examples + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 5.2.0 + + + 1.4.0b1 + + + date + + + + + + + + 0.1.0 + 0.1.0 + + + alpha + alpha + + 2011-08-22 + GPL v2.1 + +* Initial release + + + + + 0.1.1 + 0.1.1 + + + alpha + alpha + + 2011-09-05 + GPL v2.1 + +* Fixed minimum elevation bug in visible pass detection +* Updated precision of Pogson's Ratio, refactored calculate magnitude to be more readble, as well as added comments + + + + + 0.1.2 + 0.1.2 + + + alpha + alpha + + 2011-09-25 + GPL v2.1 + +* Added Predict_Time::getEpochTimeStamp() + + + + + 0.2.0 + 0.2.0 + + + alpha + alpha + + 2012-05-21 + GPL v2.1 + +* Fixed bug in Predict_Time::unix2daynum() +* Updated example iss.tle file + + + + + 0.2.1 + 0.2.1 + + + alpha + alpha + + 2012-07-04 + GPL v2.1 + +* Corrected role of README.md + + + + + 0.2.2 + 0.2.2 + + + alpha + alpha + + 2014-03-06 + GPL v2.1 + +* Addec Predict_TLE::createChecksum() +* Updates to examples + + + +