sight-reduction-geosat.py

"""
 Sight Reduction

 Similar to the last sight reduction one, but modified to use geostationary satellites.

"""

from __future__ import division
import math
from math import sin, cos, tan, asin, acos, atan 

##############################################################################
# Utility Functions
##############################################################################

def deg2rad(x):
    return math.radians(x)

def rad2deg(x):
    return math.degrees(x)

def dms2dec(degrees, minutes, seconds=0):
    return degrees + minutes/60 + seconds/60/60

def dec2dms(_dec):
    """Takes in a decimal time and converts the number to hours, minutes and
    seconds. Returns the result as a tuple in the form (hour, min, sec)
    """
    decimal = _dec % 1
    hour = round(_dec - decimal)

    _min = decimal * 60
    minutes_decimal = _min % 1
    minutes = round(_min - minutes_decimal)

    _sec = minutes_decimal * 60
    seconds_decimal = _min % 1
    seconds = round(_sec - seconds_decimal)

    return (hour, minutes, seconds)

def sind(a):
    """Returns the sine of angle, a, given in degrees."""
    return math.sin(math.radians(a))

def cosd(a):
    """Returns the cosine of angle, a, given in degrees."""
    return math.cos(math.radians(a))

def tand(a):
    """Returns the tangent of angle, a, given in degrees."""
    return math.tan(math.radians(a))

def asind(a):
    """Returns the arcsine of a in degrees."""
    return math.degrees(math.asin(a))

def acosd(a):
    """Returns the arccosine of a in degrees."""
    return math.degrees(math.acos(a))

def atand(a):
    """Returns the arctangent of a in degrees."""
    return math.degrees(math.atan(a))



##############################################################################
# Utility Classes
##############################################################################


class Time(object):

    def __init__(self, hour, minute, second):
        self.h = hour
        self.m = minute
        self.s = second

    def __sub__(self, other):
        """TODO - Make this method more robust. It's fine for the HW problem."""
        dms1 = dms2dec(self.h, self.m, self.s)
        dms2 = dms2dec(other.h, other.m, other.s)
        diff = dms1 - dms2
        h, m, s = dec2dms(diff)
        return diff


class Angle(object):

    def __init__(self, angle):
        self.set(angle)

    def set(self, value):
        tmp = value
        while tmp > 360:
            tmp -= 360
        while tmp < 0:
            tmp += 360
        self.value = tmp

    def get(self):
        return self.value



class Satellite(object):
    """

    The Satellite class is used to define an individual satellite.

    """

    def __init__(self, name='', ra=None):
        self.name = name
        self.ra = ra
        self.dec = 0
        self.Ho = None


    def alt(self, _lat, _lon):

        # Elevation Angle for Star
        theta_e = deg2rad(_lat)       # Earth station latitude     [radians]
        phi_e   = deg2rad(_lon)       # Earth station longitude    [radians]
        phi_s   = deg2rad(self.ra)    # Satellite longitude        [radians]
        phi_se  = phi_s - phi_e       # Difference in longitude    [radians]
        S = sind(self.dec)
        C = cosd(self.dec)*cos(theta_e)
        eta_s   = asin(S*sin(theta_e) + C*cos(phi_se));

        # Elevation Angle for Satellite
        num = sin(eta_s) - R/r;
        den = cos(eta_s);
        eta = atand(num/den);           # Elevation angle, eta          [deg]

        return eta


    def az(self, _lat, _lon):
        # Elevation Angle for Star
        theta_e = deg2rad(_lat)       # Earth station latitude     [radians]
        phi_e   = deg2rad(_lon)       # Earth station longitude    [radians]
        phi_s   = deg2rad(self.ra)    # Satellite longitude        [radians]
        phi_se  = phi_s - phi_e       # Difference in longitude    [radians]
        S = sind(self.dec)
        C = cosd(self.dec)*cos(theta_e)
        eta_s = asin(S*sin(theta_e) + C*cos(phi_se));

        # Elevation Angle for Satellite
        num = sin(eta_s) - R/r;
        den = cos(eta_s);
        eta = atand(num/den);           # Elevation angle, eta          [deg]

        num = sin(phi_se);
        den = cos(theta_e)*tand(self.dec) - sin(theta_e)*cos(phi_se);
        Az = atand(num/den);            # Azimuth                       [deg]

        return 180 + Az


class Position(object):

    def __init__(self, lat, lon):
        self.lat = lat
        self.lon = lon


##############################################################################
# Main - Iterative sight reduction
##############################################################################

if __name__ == '__main__':

    # Radius of Earth
    R = 6378.1e3


    # Positions
    actual = Position(28.6024, -81.200)
    ap = Position(28.6260,-80.6205)

    # Satellites - 27W and 115 W
    sat1 = Satellite('27W', ra=-27)
    sat2 = Satellite('115W', ra=-115)
    r = 42164e3

    MAX_ITERS = 1
    TOLERANCE = 1.0

    Bf = ap.lat 
    Lf = ap.lon

    for i in range(1, MAX_ITERS + 1):

        ## Intercepts -> p = Ho - Hc
        p1 = sat1.alt(actual.lat, actual.lon) - sat1.alt(ap.lat, ap.lon)
        p2 = sat2.alt(actual.lat, actual.lon) - sat2.alt(ap.lat, ap.lon)

        # Azimuth
        Z1 = sat1.az(actual.lat, actual.lon)
        Z2 = sat2.az(actual.lat, actual.lon)

        # A, B, C ...
        A = cosd(Z1)**2 + cosd(Z2)**2 
        B = cosd(Z1)*sind(Z1) + cosd(Z2)*sind(Z2) 
        C = sind(Z1)**2 + sind(Z2)**2 
        D = p1*cosd(Z1) + p2*cosd(Z2) 
        E = p1*sind(Z1) + p2*sind(Z2)
        G = A*C - B**2 

        # New lat/lon estimates
        Li = Lf + (A*E - B*D) / (G*cosd(Bf))     # Longitude
        Bi = Bf + (C*D - B*E)/G                  # Latitude

        # Distance between estimates
        d = 60 * math.sqrt( (Li-Lf)**2 * cosd(Bf)**2 + (Bi - Bf)**2 )

        # If dist is close enough, break
        if d < TOLERANCE:
            break

        # Prepare for next iteration
        Bf = Bi
        Lf = Li


    print 'Number of iterations: %3d' % i
    print 'Tolerance:            %10.2e Nm' % TOLERANCE
    print 'Distance Error:       %10.2e Nm' % d
    print 'Latitude:             %6.4f     deg' % Bf
    print 'Longitude:            %6.4f     deg' % Lf