Notes on Celestial Navigation

CNLDIFF - Lunar differentials using Chichester's method

Description: Estimates UT and the observer's longitude from Moon and 2nd object sights.

Inputs:

Outputs:

Sample execution:

Lat? 50
H Moon? 33,9.7
H 2nd object? 47,14.9
Time guess (hh:mm:ss)? 15:00
Time at T1 (hh:mm:ss)? 14:30
Time at T2 (hh:mm:ss)? 15:30
GHA Moon at T1? 320,44.8
Dec Moon at T1? 0,55.2
GHA 2nd object at T1? 36,30.4
Dec 2nd object at T1? 23,3.6
Lat1,Lon1:  50° 17.0' 6° 03.6'
GHA Moon at T2? 335,20.3
Dec Moon at T2? 0,40.6
GHA 2nd object at T2? 51,30.3
Dec 2nd object at T2? 23,3.5
Lat2,Lon2:  49° 54.3' -8° 41.2'

Lon:  -4° 59.0'
Time:  15h 14m 56s

(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/.

# CNLDIFF - Calculate Longitude and Time by Chichester's Lunar Difference Method

from math import *
from CN_LIB import *

# Van Allen's method for intersecting circles of equal altitude

def va(gha1,dec1,h1,gha2,dec2,h2):
    theta1 = -gha1 if gha1 < 180 else 360-gha1
    theta2 = -gha2 if gha2 < 180 else 360-gha2

    # Rectangular parameters of planes of 2 circles

    a1 = Cos(theta1)*Cos(dec1)
    a2 = Cos(theta2)*Cos(dec2)
    b1 = Sin(theta1)*Cos(dec1)
    b2 = Sin(theta2)*Cos(dec2)
    c1 = Sin(dec1)
    c2 = Sin(dec2)
    p1 = Sin(h1)
    p2 = Sin(h2)

    # Parameters of line of intersection of 2 planes

    A = a1*b2-a2*b1
    B = b2*c1-b1*c2
    C = b2*p1-b1*p2
    D = a1*c2-a2*c1
    E = b1*c2-b2*c1
    F = c2*p1-c1*p2

    # Quadratic equation (intersection of line with unit sphere)

    a = 1+D*D/(E*E)+A*A/(B*B)
    b = -2*D*F/(E*E)-2*A*C/(B*B)
    c = F*F/(E*E)+C*C/(B*B)-1

    sqr = sqrt(b*b-(4*a*c))

    # Roots of quadratic and corresponding y,z values

    x1 = (-b+sqr)/(2*a)
    x2 = (-b-sqr)/(2*a)
    y1 = (F-D*x1)/E
    y2 = (F-D*x2)/E
    z1 = (C-A*x1)/B
    z2 = (C-A*x2)/B

    # Convert from rectangular to polar coordinates

    lat1 = aSin(z1)
    lon1 = aTan2(y1,x1)

    lat2 = aSin(z2)
    lon2 = aTan2(y2,x2)

    if abs(lat1-lat) < abs(lat2-lat):
        return(lat1,lon1)
    else:
        return(lat2,lon2)

# Latitude
lat = stod("Lat? ")

# Moon
hm = stod("H Moon? ")

# 2nd object
hs = stod("H 2nd object? ")

# Estimated time
utguess = stot("Time guess (hh:mm:ss)? ")
t1 = stot("Time at T1 (hh:mm:ss)? ")
t2 = stot("Time at T2 (hh:mm:ss)? ")

# T1 intersection
gham = stod("GHA Moon at T1? ")
decm = stod("Dec Moon at T1? ")
ghas = stod("GHA 2nd object at T1? ")
decs = stod("Dec 2nd object at T1? ")

lat1,lon1 = va(gham,decm,hm,ghas,decs,hs)
print("Lat1,Lon1: ",dtos(lat1),dtos(lon1))

# T2 intersection
gham = stod("GHA Moon at T2? ")
decm = stod("Dec Moon at T2? ")
ghas = stod("GHA 2nd object at T2? ")
decs = stod("Dec 2nd object at T2? ")

lat2,lon2 = va(gham,decm,hm,ghas,decs,hs)
print("Lat2,Lon2: ",dtos(lat2),dtos(lon2))

# Interpolation to known Latitude
lon = lon1+((lat-lat1)/(lat2-lat1))*(lon2-lon1)
t = t1+((lat-lat1)/(lat2-lat1))*(t2-t1)

print("")
print("Lon: ",dtos(lon))
print("Time: ",ttos(t))