CNVA - Calculate the Observer's Position Using Van Allen's Method
Description: Van Allen's Method uses two sights to calculate the observer's position from two circles of equal altitude. It gives two solutions, the choice of which is made by inspection (or a third sight). The solutions are exact if the sights are simultaneous. For a running fix, some sort of estimate of the original DR and the observer's displacement is required (\(DMG > 0\)). Set \(DMG = 0\) for simultaneous sights.
Inputs:
| GHA1 | GHA of the celestial object from sight 1 |
| Dec1 | Declination of the celestial object from sight 1 |
| H1 | Height/Altitude of the celestial object from sight 1 |
| DMG | Distance made good between sights 1 and 2 (running fix) |
| CMG | Course made good between sights (\(DMG > 0\)) |
| Lat | The original DR latitude of the observer (\(DMG > 0\)) |
| Lon | The original DR longitude of the observer (\(DMG > 0\)) |
| GHA2 | GHA of the celestial object from sight 2 |
| Dec2 | Declination of the celestial object from sight 2 |
| H2 | Height/Altitude of the celestial object from sight 2 |
Outputs:
| Lat1 | The latitude of the first point of intersection of the two circles |
| Lon1 | The longitude of the first point of intersection of the two circles |
| Lat2 | The latitude of the second point of intersection of the two circles |
| Lon2 | The longitude of the second point of intersection of the two circles |
Sample execution:
GHA1? 29,30.4 Dec1? -11,17 H1? 31,19 DMG? 0 GHA2? 129,34.6 Dec2? -16,45 H2? 22,48.4 Lat1: 25° 54.7' Lon1: -76° 10.7' Lat2: -66° 38.4' Lon2: -58° 29.2'
(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/.
# CNVA - Calculate a fix from 2 Sights using Van Allen's method
from math import *
from CN_LIB import *
# First sight
gha1 = stod("GHA1? ")
dec1 = stod("Dec1? ")
h1 = stod("H1? ")
# Change GP of 1st sight if running fix
dmg = float(input("DMG? "))
if (dmg > 0):
cmg = stod("CMG? ")
lat = stod("DR Lat? ")
lon = stod("DR Lon? ")
lha = (gha1+lon+360)%360
hc = aSin(Sin(lat)*Sin(dec1)+Cos(lat)*Cos(dec1)*Cos(lha))
z = aCos((Sin(dec1)-Sin(hc)*Sin(lat))/(Cos(hc)*Cos(lat)))
zn = 360-z if lha <= 180 else z
dlat = dmg*Cos(cmg)/60
dlon = dmg*Sin(cmg)/Cos(lat)/60
lat = lat+dlat
lon = lon+dlon
sind = Sin(lat)*Sin(hc)+Cos(lat)*Cos(hc)*Cos(zn)
cosd = sqrt(1-sind*sind)
dec1 = aTan(sind/cosd)
t = aCos((Sin(hc)-Sin(lat)*sind)/(Cos(lat)*cosd))
gha1 = -lon+t if zn >= 180 else -lon-t
# Second sight
gha2 = stod("GHA2? ")
dec2 = stod("Dec2? ")
h2 = stod("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)
print("")
print("Lat1: "+dtos(lat1))
print("Lon1: "+dtos(lon1))
print("")
print("Lat2: "+dtos(lat2))
print("Lon2: "+dtos(lon2))