<?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, "/sat/");
}
//------------------------------------------------------------------------------------------//
// Print the results inside an HTML table //
//------------------------------------------------------------------------------------------//
print_head();
print_table();
// -------------------------------- END of Table -----------------------------------------//
// ------------------------------- Print "Functions" --------------------------------------//
function print_head() // HTTP header is included here in order to set the cookie BEFORE.
{
print "<!DOCTYPE HTML PUBLIC '-//W3C//DTD HTML 4.01 Transitional//EN''http://www.w3.org/TR/html4/loose.dtd'>";
print "<HTML lang='en'><HEAD><META http-equiv='content-type' content='text/html; charset=iso-8859-1'>
<TITLE>Sat Calc PHP</TITLE>";
if ($GLOBALS["bw"] == "TRUE") print "<LINK href='bw.css' rel='stylesheet' type='text/css'></HEAD><BODY>"; else print "<LINK href='sat_1.css' rel='stylesheet' type='text/css'></HEAD><BODY>";
print "<!-- generated at http://www.golombek.com/sat -->";
print "<H1>Clarke Belt Snapshot</H1><H3>Generated on <strong>";
print date("r"); print "</strong></h3>";
if ($GLOBALS["do_look_angle"] == "TRUE")
{print "<H3>Look angles for {$GLOBALS["observer_site"]} ( ";
if ($GLOBALS["observer_lat_sign"] == North) {printf ("%03.2f", $GLOBALS["observer_lat"]); print "° North, ";}
if ($GLOBALS["observer_lat_sign"] == South) {printf ("%03.2f", -$GLOBALS["observer_lat"]); print "° South, ";}
if ($GLOBALS["observer_long_sign"] == West) {printf ("%03.2f", -$GLOBALS["observer_long"]); print "° West, ";}
if ($GLOBALS["observer_long_sign"] == East) {printf ("%03.2f", $GLOBALS["observer_long"]); print "° East, ";}
print "Alt. {$GLOBALS["observer_alt"]} km. )</H3>";
}
}
function print_table_head ()
{
Print "<br><table width='auto%' border='0' align='center' cellpadding='2' cellspacing='2'>
<TR> <TD class='sath'>Satellite</TD><TD class='lonh'>Longitude</TD><TD class='lonh'>Inclination</TD>";
if ($GLOBALS["do_look_angle"])
{ print "<TD class='azh'> Azimuth </TD>
<TD class='azh'> Elevation </TD>";
if ($GLOBALS["calc_range"]) print "<TD class='ranh'> Range</TD>";
}
print "</TR>";
}
function print_end()
{
if ($GLOBALS["do_look_angle"]) print "<H3>Only satellites with elevation > 2° are shown</H3>";
print "<H3>{$GLOBALS["no_of_sat"]}";
print " Satellites in file <i>{$GLOBALS["tle_url"]} </i> - ";
print "Average TLE Age: <b>{$GLOBALS["average_age"]}</b> days</H3><br>";
print "<TABLE width='100%' border='0' cellspacing='0' cellpadding='5'>
<TR><TD align='left'><H2>file generated at
<A href='http://www.golombek.com/sat'>http://www.golombek.com/sat</A></H2></TD>
<TD align='right'><a href='http://validator.w3.org/check/referer'>
<IMG src='images/html401.png' width='80' height='15' border='0' alt='valid HTML 4.1'></A>
<A href='http://jigsaw.w3.org/css-validator/check/referer'>
<IMG src='css.png' width='80' height='15' border='0' alt='valid CSS1'></A></TD>
</TR></TABLE></BODY></HTML>";
// ------------------------------- Google AdSense --------------------------------------//
print "<div align='center'>";
print "<script type='text/javascript'><!--
google_ad_client = 'pub-8601536057106122';
google_ad_width = 728;
google_ad_height = 90;
google_ad_format = '728x90_as';
google_ad_type = 'text_image';
google_ad_channel ='';
google_color_border = '6666FF';
google_color_bg = 'CCCCCC';
google_color_link = '0000DD';
google_color_url = '008000';
google_color_text = '000000';
//--></script>
<script type='text/javascript'
src='http://pagead2.googlesyndication.com/pagead/show_ads.js'>
</script>";
print "</div>";
print "<script src='http://www.google-analytics.com/urchin.js' type='text/javascript'></script>
<script type='text/javascript'>_uacct = 'UA-513682-1';
urchinTracker();
</script>";
}
////////////////////////////////////////////////////////////////////////////////////////////////
function print_table()
{
print_table_head ();
$lcl_res_arr = $GLOBALS["lista"];
$do_look_angle = $_GET ['calc_ang'];
$calc_range = $_GET ['calc_range'];
$range_units = $_GET ['range_units'];
$symbol_rate = $_GET ['sym_rate'];
if ($_GET ['round_trip'] == "TRUE") $round_trip = 2; else $round_trip = 1;
for ($i=0; $i<($GLOBALS["no_of_sat"] - 1); $i++)
{
$sub_sat_long = $lcl_res_arr[$i]['sub_sat_long'];
$obs_elevation = $lcl_res_arr[$i]['Elevation'];
$obs_azimuth = $lcl_res_arr[$i]['Azimuth'];
$range_magnitude = $lcl_res_arr[$i]['Range_km'];
if (($obs_elevation > 1.5) || (!$do_look_angle))
{
print "<tr><TD class='sat'>";
print "{$lcl_res_arr [$i]['Sat']}</TD>";
print "<TD class='lon'>";
if ($sub_sat_long >= 0) {printf ("%03.2f", $sub_sat_long); print "° East";}
if ($sub_sat_long < 0) {printf ("%03.2f", -$sub_sat_long); print "° West";}
print "</td>";
if ( $lcl_res_arr [$i]['Inclination'] > 2 ) print "<TD class='incl'>";
if ( $lcl_res_arr [$i]['Inclination'] <= 2 ) print "<TD class='lon'>";
printf ("%1.3f", $lcl_res_arr [$i]['Inclination']);
print "°</td>";
if ($do_look_angle) // Prints the Look angles
{
print "<TD class='az'>";
printf ("%03.2f", $obs_azimuth);
print "°</td>";
print "<TD class='el'>";
printf ("%02.2f", $obs_elevation);
print "°</td>";
if ($calc_range) // Print the 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) ;
print "<TD class='ran'>";
printf ("%.2f", $range_print);
if ( $range_units == 'Kilometers' ) print" km</TD>";
if ($range_units == 'milliseconds' ) print" ms</TD>";;
if ($range_units == 'symbols' ) print" sym</TD>";
}
}
print "</tr>";
}
}
print "</table>";
print_end();
}
function nofile($file)
{
print "<TITLE>File Error</TITLE>
<LINK href='sat_1.css' rel='stylesheet' type='text/css'></HEAD><BODY>";
print "<H1>Sorry, an error has been encountered</H1><br>";
print "The TLE file at <a href=$file>$file<a> could not be opened, it appears that it doesn't exist.Check that the URL is correct.<BR><BR>";
print "<DIV align='center'><IMG SRC='error404.png'></DIV>";
// ------------------------------- Google Analytics --------------------------------------//
print "<script src='http://www.google-analytics.com/urchin.js' type='text/javascript'></script>
<script type='text/javascript'>_uacct = 'UA-513682-1';
urchinTracker();
</script></BODY></HTML>";}
function array_sort($array, $key)
{
for ($i = 0; $i < sizeof($array); $i++)
{ $sort_values[$i] = $array[$i][$key]; }
asort ($sort_values);
reset ($sort_values);
while (list ($arr_key, $arr_val) = each ($sort_values))
{ $sorted_arr[] = $array[$arr_key]; }
return $sorted_arr;
}
function rarray_sort($array, $key)
{
for ($i = 0; $i < sizeof($array); $i++)
{ $sort_values[$i] = $array[$i][$key]; }
arsort ($sort_values);
reset ($sort_values);
while (list ($arr_key, $arr_val) = each ($sort_values))
{ $sorted_arr[] = $array[$arr_key]; }
return $sorted_arr;
}
?>