<?php
// -------------------- Constants ----------------------------------------------------- //
define ("earth_grav_const" , 3.986e5); // Earth's Graviational Constant (Km^3/s^2)
define ("eg_4pi" , 10096.66709265246); // Earth's Graviational Constant / 4*pi^2
define ("sideral_day_sec" , 86164.0984); // One Earth Rotation in seconds
define ("equatorial_radius", 6378.137); // WGS-84
define ("flattening" , 298.257223563); // 1/flattening per WGS-84
define ("ghaa_deg" , 99.4033); // The angle between the Greenwich meridian
define ("date_of_GHAA" , "1/1/1990 00:00"); // and 'Aries' at this date in Degrees (Vernal Equinox)
define ("Tropical_year" , 365.242197); // Tropical Year in days
define ("c" , 2.99792458e5); // spped of light in Km/s
// http://www.colorado.edu/geography/gcraft/notes/datum/elist.html
// The amount the earth rotates during the 3 minutes and 5.59 seconds
// difference between the 24 hrs (Solar) day and the sideral day
$extra_earth_rot_per_day = (2 * M_PI) / Tropical_year;
// The total earth rotation in one Solar Day = 1 sideral day + the above figure
$earth_rot_rad_sec = ($extra_earth_rot_per_day + (2 * M_PI))/86400; // Earth Rotation Rate (radians/sec)
// Initializing the Current Date Array
$curr_date_array = getdate();
$curr_year = $curr_date_array ['year'];
$curr_day = $curr_date_array ['yday'] + 1; // php counts days from 0, TLE epoch counts from 1
$curr_sec =($curr_date_array ['hours'] * 3600) +
($curr_date_array ['minutes'] * 60) +
$curr_date_array ['seconds']; // - 3600; // Mistyrious One Hour shift (20.10.03 - 26.10.03) Possible DST bug
$curr_time_unix = $curr_date_array [0];
$curr_frac_day = $curr_sec / 86400; // This is how the sideral time appears in the TLEs, fraction of solar day
// Bringing the RA in sync with the current time.
// Seconds since reference GHAA was specified
$delta_ghaa_sec = $curr_time_unix - (strtotime (date_of_GHAA)); // time is in sedconds since 1.1.70 (UNIX time)
// The current Angle to Aries from the Greenwich meridian. This is the X axis of the Geocentric Inetial coordinates
$current_ghaa_rad = deg2rad(ghaa_deg) + ($delta_ghaa_sec * $earth_rot_rad_sec);
$cos_ghaa = cos(-$current_ghaa_rad);
$sin_ghaa = sin(-$current_ghaa_rad);
// Get the data from the form for Look Angles and range
$do_look_angle = $_GET ['calc_ang'];
$calc_range = $_GET ['calc_range'];
$range_units = $_GET ['range_units'];
$symbol_rate = $_GET ['sym_rate'];
$sort_order = $_GET ['sort_order'];
$use_cookie = $_GET ['use_cookie'];
$bw = $_GET ['print'];
if ($_GET ['round_trip'] == "TRUE") $round_trip = 2; else $round_trip = 1;
$observer_site = $_GET ['loc_name'];
$observer_lat = $_GET ['loc_lat'];
$observer_long = $_GET ['loc_lon'];
$observer_alt = $_GET ['loc_alt'];
$observer_alt = $observer_alt/1000;
$observer_lat_sign = $_GET ['lat_menu'];
$observer_long_sign = $_GET ['lon_menu'];
if ($observer_lat_sign == South) $observer_lat = -$observer_lat;
if ($observer_long_sign == West ) $observer_long = -$observer_long;
// Observer's UNIT vector In Geocentric Equatorial Coordinates
// These are NOT Inertial coordinates, they are referenced to Greenwich
$obs_lat_rad = deg2rad($observer_lat);
$obs_long_rad = deg2rad($observer_long);
$observer_up_x = cos($obs_lat_rad) * cos($obs_long_rad);
$observer_up_y = cos($obs_lat_rad) * sin($obs_long_rad);
$observer_up_z = sin($obs_lat_rad);
$observer_east_x = -sin($obs_long_rad);
$observer_east_y = cos($obs_long_rad);
$observer_north_x= -sin($obs_lat_rad) * cos($obs_long_rad);
$observer_north_y= -sin($obs_lat_rad) * sin($obs_long_rad);
$observer_north_z= cos($obs_lat_rad);
//Observer's XYZ coordinates on the surface of the earth.
$polar_radius = equatorial_radius * (1 - (1/flattening));
$rx = ( pow(equatorial_radius,2) / sqrt( (pow(equatorial_radius,2) * pow(cos($obs_lat_rad),2) + (pow($polar_radius,2) * pow(sin($obs_lat_rad ),2))))) + $observer_alt;
$rz = ( pow($polar_radius,2) / sqrt( (pow(equatorial_radius,2) * pow(cos($obs_lat_rad),2) + (pow($polar_radius,2) * pow(sin($obs_lat_rad ),2))))) + $observer_alt;
$observer_x = $observer_up_x * $rx;
$observer_y = $observer_up_y * $rx;
$observer_z = $observer_up_z * $rz;
if (!do_look_angle) $obs_elevation = 2; // To force printing if Look angles aren't calculated
// Open the TLE file specified in the form from the page that invoked
// the script, pasre it and put it into an intermediate Array
$tle_url= $_GET ['tle_url'];
$tle_data = fopen ($tle_url,"r");
if (!$tle_data) nofile($tle_url); // File doesn't exist
if (!$tle_data) return;
$no_of_sat=0;
while ((!feof ($tle_data)) and ($tle_data))
{
$buffer1 = fgets($tle_data, 128);
$buffer2 = fgets($tle_data, 128);
$buffer3 = fgets($tle_data, 128);
$raw_sat_ele[$no_of_sat] = array ('Name' => $buffer1,
'line1' => $buffer2,
'line2' => $buffer3);
$no_of_sat++;
}
// -- Parsing the raw TLE data into elements -- //
for ($i=0; $i<$no_of_sat; $i++)
{
$sat_name = substr ($raw_sat_ele [$i]['Name'], 0, 26 );
// $sat_cata = substr ($raw_sat_ele [$i]['line1'], 2, 5 );
$sat_launch_year = substr ($raw_sat_ele [$i]['line1'], 9, 2 );
// $sat_launch_no = substr ($raw_sat_ele [$i]['line1'], 11, 3 );
// $sat_payload_no = substr ($raw_sat_ele [$i]['line1'], 14, 1 );
$epoch_year = substr ($raw_sat_ele [$i]['line1'], 18, 2 );
$epoch_sider_day = substr ($raw_sat_ele [$i]['line1'], 20, 3 );
$epoch_sider_time = substr ($raw_sat_ele [$i]['line1'], 23, 9 );
$sat_inclination = substr ($raw_sat_ele [$i]['line2'], 8, 8 );
$sat_ra_asc_node = substr ($raw_sat_ele [$i]['line2'], 17, 8 );
$sat_eccentricity = substr ($raw_sat_ele [$i]['line2'], 26, 7 );
$sat_arg_perigee = substr ($raw_sat_ele [$i]['line2'], 34, 8 );
$sat_mean_anomaly = substr ($raw_sat_ele [$i]['line2'], 43, 8 );
$sat_mean_motion = substr ($raw_sat_ele [$i]['line2'], 52, 11 );
$sat_rev_at_epoch = substr ($raw_sat_ele [$i]['line2'], 63, 5 );
// And passing it to an "elements" Array
$sat_elements [$i] = array ( 'Name' => $sat_name,
'Eccentricity' => $sat_eccentricity,
'Inclination' => $sat_inclination,
'Mean_anomaly' => $sat_mean_anomaly,
'Rev_at_epoch' => $sat_rev_at_epoch,
'Mean_motion' => $sat_mean_motion,
'Arg_perigee' => $sat_arg_perigee,
'RA_asc_node' => $sat_ra_asc_node,
'epoch_year' => $epoch_year,
'epoch_day' => $epoch_sider_day,
'epoch_time' => $epoch_sider_time );
};
// --------------------- Main Loop, passes through every TLE -------------------------- //
for ($i=0; $i<$no_of_sat-1; $i++)
{
// Calculating elapsed time since epoch
$ep_yr = $sat_elements [$i]['epoch_year'];
if ($ep_yr>60) $ep_yr = $ep_yr + 1900;
if ($ep_yr<60) $ep_yr = $ep_yr + 2000;
$ep_dy = $sat_elements [$i]['epoch_day'];
$ep_tm = $sat_elements [$i]['epoch_time'];
$ep_tm_sec = $ep_tm * 86400;
$ep_unix_y = getdate( mktime (0, 0, 0, 1, 1, $ep_yr) ); // Used to calcluate time since epoch taking
$ep_unix_s = $ep_unix_y[0] + (($ep_dy + $ep_tm) * 86400); // into account leap years
$delta_sec = $curr_time_unix - $ep_unix_s;
$delta_frac_yr = $delta_sec / (86400 * Tropical_year);
// Calculating Axes of the elliptic orbit
$mean_motion = $sat_elements [$i]['Mean_motion']; // Mean Motion in the TLE is given in Revolution / Day
$seconds_pre_rev = $mean_motion * sideral_day_sec; // Number of seconds per revolution
$eccentricity = $sat_elements [$i]['Eccentricity'];
$eccentricity = $eccentricity/1e7; // Eccentricity in the TLEs assumes decimal point
$semi_major_axis = pow((eg_4pi*pow($seconds_pre_rev, 2)),(1/3));
$semi_minor_axis = $semi_major_axis * sqrt(1-pow($eccentricity,2));
// Find the 'Eccentric Anomaly' by iterating (Newton's Method)
// Mean Anomaly = Ecc Anomaly - (Eccentricity * Sin (Ecc Anomaly))
$mean_mot_rad_sec = ($mean_motion * 2 * M_PI) / 86400; // Mean Motion in RAD / sec
$mean_anomaly_deg = $sat_elements [$i]['Mean_anomaly']; // Mean Anomaly in the TLE is given in Degrees
$mean_anomaly_rad = deg2rad($mean_anomaly_deg);
$rads_since_epoch = ($mean_mot_rad_sec * $delta_sec) + $mean_anomaly_rad; // (Mean Motion * Elapsed time) + Mean Anomaly
$frac_rev = $rads_since_epoch - (2 * M_PI * ( floor( $rads_since_epoch / (2 * M_PI) ) ));
// Do the iterations with this initisl value
$ecc_an = $frac_rev;
do
{
$cos_ea = cos ($ecc_an);
$sin_ea = sin ($ecc_an);
$denom = 1 - ($cos_ea * $eccentricity);
$iter = ( $ecc_an - ( $eccentricity * $sin_ea) - $rads_since_epoch ) / $denom;
$ecc_an = $ecc_an - $iter;
} while (abs($iter) > 0.0001);
$sat_range = $semi_major_axis * $denom; // Satellite range from the CENTER of the earth
// Calculating Satellite position vector on the Orbital Plane
$sat_orb_plane_X = $semi_major_axis * ($cos_ea - $eccentricity);
$sat_orb_plane_Y = $semi_minor_axis * $sin_ea;
// Partial Rotation Matrix to transform from the Orbital Plane to Inertial (Celestial) Coordinates
$incl = $sat_elements [$i]['Inclination'];
$argpg = $sat_elements [$i]['Arg_perigee'];
$raan = $sat_elements [$i]['RA_asc_node'];
$cos_arg_per = cos (deg2rad($argpg));
$sin_arg_per = sin (deg2rad($argpg));
$cos_raan = cos (deg2rad($raan));
$sin_raan = sin (deg2rad($raan));
$cos_incl = cos (deg2rad($incl));
$sin_incl = sin (deg2rad($incl));
$cel_x_x = ( $cos_arg_per * $cos_raan) - ($sin_arg_per * $sin_raan * $cos_incl); $cel_x_y = (-$sin_arg_per * $cos_raan) - ($cos_arg_per * $sin_raan * $cos_incl);
$cel_y_x = ( $cos_arg_per * $sin_raan) + ($sin_arg_per * $cos_raan * $cos_incl); $cel_y_y = (-$sin_arg_per * $sin_raan) + ($cos_arg_per * $cos_raan * $cos_incl);
$cel_z_x = ( $sin_arg_per * $sin_incl); $cel_z_y = ( $cos_arg_per * $sin_incl);
// Calculatins Satellite position vector in Celestial Coordinates
$sat_celestial_X = ($sat_orb_plane_X * $cel_x_x) + ($sat_orb_plane_Y * $cel_x_y);
$sat_celestial_Y = ($sat_orb_plane_X * $cel_y_x) + ($sat_orb_plane_Y * $cel_y_y);
$sat_celestial_Z = ($sat_orb_plane_X * $cel_z_x) + ($sat_orb_plane_Y * $cel_z_y);
// Satellite Coordinates in Geocentric Equatorial Coordinates (from RA to LONG, etc.)
$sat_geoc_x = ( $sat_celestial_X * $cos_ghaa ) - ( $sat_celestial_Y * $sin_ghaa );
$sat_geoc_y = ( $sat_celestial_X * $sin_ghaa ) + ( $sat_celestial_Y * $cos_ghaa );
$sat_geoc_z = $sat_celestial_Z ;
// Calculate Subsatellite Point
$sub_sat_long = rad2deg(atan2($sat_geoc_y, $sat_geoc_x)); // East / West
$sub_sat_lat = rad2deg(asin ($sat_geoc_z/$sat_range )); // North / Souht
if ($do_look_angle) // Calculating Range Vector
{
$range_x = $sat_geoc_x - $observer_x;
$range_y = $sat_geoc_y - $observer_y;
$range_z = $sat_geoc_z - $observer_z;
$range_magnitude = sqrt(pow($range_x,2)+pow($range_y,2)+pow($range_z,2));
$range_norm_x = $range_x / $range_magnitude; //Normalized Range Vector
$range_norm_y = $range_y / $range_magnitude;
$range_norm_z = $range_z / $range_magnitude;
$range_up = ($range_norm_x * $observer_up_x) + ($range_norm_y * $observer_up_y) + ($range_norm_z * $observer_up_z);
$range_east = ($range_norm_x * $observer_east_x) + ($range_norm_y * $observer_east_y);
$range_north= ($range_norm_x * $observer_north_x)+ ($range_norm_y * $observer_north_y) + ($range_norm_z * $observer_north_z);
// Calculating the Look Angles
$obs_elevation = rad2deg(asin($range_up));
$obs_azimuth = rad2deg(atan2($range_east,$range_north));
if ( $obs_azimuth < 0 ) $obs_azimuth = 360 + $obs_azimuth;
// Calculating the range according to the units
if ($calc_range)
{
if ($range_units == 'Kilometers' ) $range_print = round (($range_magnitude * $round_trip),2);
if ($range_units == 'milliseconds' ) $range_print = round ((($range_magnitude * $round_trip * 1000) / c) , 2) ;
if ($range_units == 'symbols' ) $range_print = round ((($range_magnitude * $round_trip * 1000 * $symbol_rate) / c) , 2) ;
}
}
//------------------------ Fill the results array for use later-----------------------------//
$results [$i] = array ( 'Sat' => $sat_elements [$i]['Name'] ,
'Epoch_age_s' => $delta_sec,
'sub_sat_long' => $sub_sat_long,
'Inclination' => $sat_elements [$i]['Inclination'],
'Azimuth' => $obs_azimuth,
'Elevation' => $obs_elevation,
'Range_km' => $range_magnitude
);
$total_epoch_age = $total_epoch_age + $results [$i] ['Epoch_age_s'];
}
// Finalizing
fclose ($tle_data);
$average_age = round((($total_epoch_age / $no_of_sat) / 86400),2);
$lista = $results;
if ($sort_order == "Alphabetic") $lista = array_sort ($results, 'Sat');
if ($sort_order == "Longitude East to West") $lista = rarray_sort ($results, 'sub_sat_long');
if ($sort_order == "Longitude West to East") $lista = array_sort ($results, 'sub_sat_long');
// -- Save The Location in a cookie
if ($use_cookie)
{
setcookie ("loc_cook", $observer_site, time()+60*60*24*90, "/sat/");
setcookie ("lat_cook", $observer_lat, time()+60*60*24*90, "/sat/");
setcookie ("lon_cook", $observer_long, time()+60*60*24*90, "/sat/");
setcookie ("alt_cook", $observer_alt, time()+60*60*24*90