Transformation between Coordinate Systems

coord_transfer.py

# Created                     July 4, 2016
# Last Major Modification     July 15, 2016
# Author                      Sibo W

import datetime
import pytz
import math

def get_lst(long, time):
    '''Takes in: long   observer's longitude in degrees, east as positive
                 time   the time to be calculated in datetime type
                        use 'NOW' if current time is wanted
    Returns Local (Mean) Sidereal Time.
    The observer's east longitude has been neglected.'''
    # Define subclsss UTC of tzinfo
    UTC = pytz.timezone('UTC')

    # Initialize current time if necessary
    if time == 'NOW':
        time = datetime.datetime.utcnow().replace(tzinfo = UTC)

    # Calculate number of days since J2000
    J2000 = datetime.datetime(2000, 1, 1, 11, 58, 55, 816000, UTC)
    # Note that J2000.0 = January 1, 2000, 12:00:00.000 Terrestrial Time
    #                   = January 1, 2000, 11:58:55.816 UTC
    days = (time - J2000).total_seconds() / 86400
    gst = days * 360.98562628 + 280.46061837
    # Note that the earth rotates 366.25 * 360 / 365.25 = 360.98562628
    # degrees each solar day
    # Add the sidereal angle at J2000.0 at Prime Meridian 280.46061837
    lst = (gst + long) % 360
    return lst


def equatorial2horizontal(ra, dec, lst, lat):
    '''Given   Right Ascension       ra    in degrees
               Declinition           dec   in degrees
               Local Sidereal Time   lst   in degrees
               Observer's latitude   lat0  in degrees
       Returns tuple (altitude, azimuth), both in degrees
    '''
    ra = math.radians(ra)
    dec = math.radians(dec)
    lst = math.radians(lst)
    lat = math.radians(lat)
    hour_angle = (lst - ra)
    # calculates hour angle (lst - ra)
    temp = math.sin(dec) * math.sin(lat) + \
           math.cos(dec) * math.cos(lat) * math.cos(hour_angle)
    # where temp = sin(altitude)
    altitude = math.asin(temp)
    # gets altitude and converts it into degrees
    temp = -math.sin(hour_angle) * math.cos(dec) / math.cos(altitude)
    # where temp = sin(azimuth)
    altitude = math.degrees(altitude)
    azimuth = math.degrees(math.asin(temp))
    # gets altitude and converts it into degrees
    azimuth %= 360
    return altitude, azimuth


def horizontal2equatorial(alt, azi, lst, lat):
    '''Given   Azimuth               azi   in degrees
               Altitude              alt   in degrees
               Local Sidereal Time   lst   in hours
               Observer's latitude   lat0  in degrees
       Returns tuple (Right Ascension, declination), both in degrees
    '''
    azi = math.radians(azi)
    alt = math.radians(alt)
    lst = math.radians(lst)
    lat = math.radians(lat)
    temp = math.sin(alt) * math.sin(lat) + \
           math.cos(alt) * math.cos(lat) * math.cos(azi)
    # where temp = sin(declination)
    declination = math.asin(temp)
    # gets declination and converts it into degrees
    temp = -math.sin(azi) * math.cos(alt) / math.cos(declination)
    # where temp = sin(hour angle)
    declination = math.degrees(declination)
    hour_angle = math.degrees(math.asin(temp))
    # gets hour angle and converts it into degrees
    ra = math.degrees(lst) - hour_angle
    # gets Right Ascension from hour angle and lst
    return (ra, declination)



# TEST CASES

# UTC = pytz.timezone('UTC')
# test = datetime.datetime(2016, 7, 15, 7, 0, 0, 0, UTC)
# lst = get_lst(172, test)
# print(lst)
#
# print(equatorial2horizontal(120, 20, lst, 33))
#
# print(horizontal2equatorial(10.17, 287.33, lst, 33))