A script used to solve sight reduction homework problems for navigation class.

sight_reduction.py

"""
 Sight Reduction

 Calculate position on 2000 June 21 2100 UTC 
               
 Star     Observed     Altitude
 -------  -----------  --------
 Regulus  20h 39m 23s  37.4204
 Antares  20h 45m 47s  20.3226
 Kochab   21h 10m 34s  47.2050

 Ship travelling at constant speed of 20 knots, heading 325.
 Position known to nearest whole degree 32N, 15W.

NOTE: Change the class variables in the Star class for the specific problem.

"""

from __future__ import division
import math


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

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


##############################################################################
# Star Class - This is where the bulk of the work happens
##############################################################################

class Star(object):
    """
    The Star class is used to define an individual star's data as well
    as class data common to all stars for the given problem. 

    Class Variables:
        REF_TIME    - The time for which we want to know our position.
        GHA         - A dictionary containing GHA Aries lookup data surrounding
                      the reference time.
    """

    REF_TIME = Time(21, 0, 0)

    # GHA Aries lookup table for date 
    GHA = {
        20: dms2dec(223, 1),
        21: dms2dec(238, 3.4),
        22: dms2dec(253, 5.9)
    }

    def __init__(self, name=''):
        self.name = name
        self.sha = None
        self.dec = None
        self.time = None  # Time observed
        self.Ho = None


    def look(self, _lat, _lon):

        # Time
        t = self.time - Star.REF_TIME
        self.t = t

        # Estimate lat, lon
        _lon = _lon + t*(V/60)*sind(T)/cosd(_lat)
        _lat = _lat + t*(V/60)*cosd(T)
        self.lat = _lat
        self.lon = _lon

        # Interpolation factor - fraction of an hour
        interp_factor = dms2dec(0, self.time.m, self.time.s) / dms2dec(0, 60)
        self.interp_factor = interp_factor

        # Last and next whole hour
        prev_hour = self.time.h
        next_hour = self.time.h + 1 # TODO (jk) - won't work after 11pm, go to be early.
        
        # Interpolate GHA Aries
        hourly_diff = Star.GHA[next_hour] - Star.GHA[prev_hour]
        GHA_Aries = Star.GHA[prev_hour] + interp_factor*(hourly_diff)
        self.GHA_Aries = GHA_Aries

        #GHA_Aries = Star.GHA[_gha_hour]
        GHA = GHA_Aries + self.sha
        LHA = Angle(GHA + _lon).get()

        # Elevation Angle
        S = sind(self.dec);
        C = cosd(self.dec)*cosd(LHA);
        Hc = asind(S*sind(_lat) + C*cosd(_lat));

        # Azimuth
        X = (S*cosd(_lat) - C*sind(_lat)) / cosd(Hc);
        if X > 1:
            X = 1
        elif X < -1:
            X = -1
        A = acosd(X)

        if (LHA > 180):
            Z = A
        else:
            Z = 360 - A

        self.azimuth = Angle(Z).get()
        self.Hc = Hc

    @property
    def p(self):
        return self.Ho - self.Hc

    @property
    def Z(self):
        return self.azimuth


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

if __name__ == '__main__':

    ########## Given ##########
    Lf = -15   # Ship position - estimated longitude
    Bf = 32    # Ship position - estimated latitude
    V = 20     # Ship course - velocity
    T = 325    # Ship course - heading

    # Regulus 
    regulus = Star('Regulus')
    regulus.sha = dms2dec(207, 40.9)
    regulus.dec = dms2dec(11, 52.9)
    regulus.time = Time(20, 39, 23)
    regulus.Ho = 26.8543

    # Antares
    antares = Star('Antares')
    antares.sha = dms2dec(112, 22.6)
    antares.dec = dms2dec(-26, -28.1)
    antares.time = Time(20, 45, 47)
    antares.Ho = 26.0611

    # Kochab
    kochab = Star('Kochab')
    kochab.sha = dms2dec(137, 19.6)
    kochab.dec = dms2dec(74, 5.5)
    kochab.time = Time(21, 10, 34)
    kochab.Ho = 47.5292


    MAX_ITERS = 20
    TOLERANCE = 1.0

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

        # Calculate Az/Alt of each star
        regulus.look(Bf, Lf)
        antares.look(Bf, Lf)
        kochab.look(Bf, Lf)

        ## Intercepts -> p = Ho - Hc
        p1 = regulus.p;
        p2 = antares.p;
        p3 = kochab.p;

        # Azimuth
        Z1 = regulus.Z;
        Z2 = antares.Z;
        Z3 = kochab.Z;

        # A, B, C ...
        A = cosd(Z1)**2 + cosd(Z2)**2 + cosd(Z3)**2
        B = cosd(Z1)*sind(Z1) + cosd(Z2)*sind(Z2) + cosd(Z3)*sind(Z3)
        C = sind(Z1)**2 + sind(Z2)**2 + sind(Z3)**2
        D = p1*cosd(Z1) + p2*cosd(Z2) + p3*cosd(Z3)
        E = p1*sind(Z1) + p2*sind(Z2) + p3*sind(Z3)
        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.2f     deg' % Bf
    print 'Longitude:            %6.2f     deg' % Lf