Notes on Celestial Navigation

CN_LIB - I/O and Maths functions

I/O functions

Date and Time functions

Astronomical Algorithms functions

Various functions to calculate solar coordinates using Meeus' lower accuracy algorithms (chapter 25)

Matrix operations

Functions for implementing solution of linear systems (eg Watkins & Janiczek)

(For explanation of notation, conventions etc, see python-programs).

Copyright (C) 2024 Ian Staniforth

This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/.

# I/O functions

from math import *

def stod(s):
    d = (input(s)+",0").split(",")
    return(float(d[0])+copysign(float(d[1]),float(d[0]))/60)

def dtos(d):
    sign = "-" if d < 0 else ""
    deg = int(abs(d))
    min = round((abs(d)%1)*60,1)
    if min >= 60: min-=60; deg+=1;
    return("%s%d\xb0 %#04.1f'" % (sign,deg,min))

def ttos(t):
    sign = "-" if t < 0 else ""
    h = int(abs(t))
    t = (abs(t)%1)*60
    m = int(t)
    s = round((t%1)*60)
    if s >= 60: s-=60; m+=1;
    if m >= 60: m-=60; h+=1;
    return("%s%#02dh %#02dm %#02ds" % (sign,h,m,s))

def ttosm(t):
    sign = "-" if t < 0 else ""
    m = int(abs(t))
    s = round((abs(t)%1)*60)
    if s >= 60: s-=60; m+=1;
    return("%s%#02dm %#02ds" % (sign,m,s))

def stot(s):
    t = (input(s)+":0:0").split(":")
    return(float(t[0])+copysign(float(t[1]),float(t[0]))/60+
           copysign(float(t[2]),float(t[0]))/3600)

# Trigonometrical functions (in degrees)

def Sin(x): return sin(radians(x))
def Cos(x): return cos(radians(x))
def Tan(x): return tan(radians(x))
def Cot(x): return 1/tan(radians(x))
def aSin(x): return degrees(asin(x))
def aCos(x): return degrees(acos(x))
def aTan(x): return degrees(atan(x))
def aTan2(x,y): return degrees(atan2(x,y))

# Date and Time functions

def julianday(y,m,d):
    if (m <= 2):
        y = y-1
        m = m+12
    a = int(y/100)
    b = 2-a+int(a/4)
    jd = int(365.25*(y+4716))+int(30.6001*(m+1))+d+b-1524.5
    return(jd)

def dayofyear(year,month,day):
    dsum = [0,31,59,90,120,151,181,212,243,273,304,334]
    # Leap year?
    if year % 4 == 0 and (year % 100 != 0 or year % 400 == 0):
        dsum = [0,31,60,91,121,152,182,213,244,274,305,335]
    n = dsum[month-1]+day
    return(n)

def ydrift(year):
    d = [0,0.75,0.5,0.25][year%4]
    return(d)

def eot(y,m,d):
    x = dayofyear(y,m,d)+ydrift(y)
    e = 10.2713*sin(0.0337*x-2.6848)+7.8075*sin(0.0179*x+2.9616)
    return(e)

def eot2(y,m,d):
    jd = julianday(y,m,d)
    T = (jd-J2000)/JC
    ra,dec = sunradec(T)
    L0 = sunmeanlongitude(T)
    e = L0-0.0057183-15*ra
    if abs(e) > 20: e = e-360
    eot = e*4
    return(eot)

# Astronomical algorithms (Meeus)

J2000 = 2451545.0   # mean equinox J2000.0
JC = 36525          # Julian century
e2000 = 23.4392911  # obliquity of ecliptic J2000.0

def sunmeanlongitude(T):
    L0 = 280.46646+36000.76983*T+0.0003032*T**2 
    L0 = L0%360
    return(L0)

def suncentre(T,M):
    C = (1.914602-0.004817*T-0.000014*T**2)*Sin(M) \
        +(0.019993-0.000101*T)*Sin(2*M) \
        +0.000289*Sin(3*M)
    return(C)

def obliquityecliptic(T):
    e0 = e2000-0.01300416667*T-0.00000001639*T**2+0.00000005036*T**3
    e0 = e0+0.00256*Cos(125.04-1934.136*T)
    return(e0)

def sunmeananomaly(T):
    M = 357.52911+35999.05029*T-0.0001537*T**2
    return(M)

def sunapparentlongitude(T,sl):
    om = 125.04-1934.136*T
    l = sl-0.00569-0.00478*Sin(om)
    return(l)

def sunradec(T):
    L0 = sunmeanlongitude(T)
    M = sunmeananomaly(T)%360
    C = suncentre(T,M)
    sl = L0+C
    l = sunapparentlongitude(T,sl)
    oc = obliquityecliptic(T)
    ra = (aTan2(Cos(oc)*Sin(l),Cos(l))%360)/15
    dec = aSin(Sin(oc)*Sin(l))
    return(ra,dec)

def ghaaries(y,m,d):
    jd = julianday(y,m,d)
    T = (jd-J2000)/JC
    g = 280.46061837+360.98564736629*(jd-2451545.0)+0.000387933*T*T-T*T*T/38710000
    gha = g%360
    return(gha)

def ghasun(y,m,d):
    jd = julianday(y,m,d)
    T = (jd-J2000)/JC
    ra,dec = sunradec(T)
    garies = ghaaries(y,m,d)
    sha = 360-ra*15
    gha = (garies+sha)%360
    return(gha)

# Matrix operations

# Dot product of two vectors

def dotp(v1,v2):
    return(sum([v1[i]*v2[i] for i in range(len(v1))]))

# Extract a matrix column

def mcol(M,n):
    return([M[j][n] for j in range(len(M))])

# Scalar multiplication of a vector by a scalar

def smult(s,M):
    return([s*x for x in M])

# Multiplication of two matrices

def mmult(M1,M2):
    M = []
    for i in range(len(M1)):
        r = []
        for j in range(len(M2[0])):
            r.append(dotp(M1[i],mcol(M2,j)))
        M.append(r)
    return(M)

# Multiply a vector by a matrix

def mvmult(M1,v):
    M = []
    for i in range(len(M1)):
        M.append(dotp(M1[i],v))
    return(M)

# Subraction of one vector from another

def vsub(v1,v2):
    return([v1[i]-v2[i] for i in range(len(v1))])

# Augment a matrix with a column

def maugment(M,v):
    MA = [x[:] for x in M]
    for i in range(len(MA)):
        MA[i].append(v[i])
    return(MA)

# Calculate the transpose of a matrix

def mtranspose(M):
    MT = []
    for i in range(len(M[0])):
        MT.append(mcol(M,i))
    return(MT)

# Solve an augmented matrix by row reduction

def msolve(M):
    MRR = [x[:] for x in M]
    for i in range(len(M)):
        for j in range(len(M)):
            if (i == j): continue
            MRR[j] = vsub(smult(MRR[i][i],MRR[j]),smult(MRR[j][i],MRR[i]))
    for i in range(len(M)):
        MRR[i] = smult(1/MRR[i][i],MRR[i])
    return(MRR)