<?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