CNWJ - Calculate a least squares position using Watkins & Janiczek's method
Description: This least squares method uses matrices to produce a position from two or more sights (4 or more are recommended for increased accuracy). Note that the matrix operations are ideally carried out using a numerical library package, but simple row reduction is included here for use on calculators.
Inputs:
| Number of sights | The number of sights to be reduced |
| GHA i | GHA of the celestial object from the ith sight |
| Dec i | Declination of the celestial object from the ith sight |
| H i | Height/Altitude of the celestial object from the ith sight |
| ... repeated for other sights |
Outputs:
| Lat | The latitude of the calculated position |
| Lon | The longitude of the calculated position |
Sample execution:
Number of sights? 4 GHA 1? 126,13.1 Dec 1? 46,01.3 H 1? 18,41.5 GHA 2? 28,15.3 Dec 2? 14,26.2 H 2? 49,45.2 GHA 3? 4,12.6 Dec 3? -11,17.4 H 3? 28,44.4 GHA 4? 286,23.1 Dec 4? 38,48.2 H 4? 35,24.5 lat = 50° 00.5' lon = -5° 00.4'
(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/.
# CNWJ - Multiple sights least squares solution using Watkins & Janiczek
from math import *
from CN_LIB import *
A = []
b = []
n = int(input("Number of sights? "))
for i in range(1,n+1):
gha = stod("GHA "+str(i)+"? ")
theta = -gha if gha < 180 else 360-gha
dec = stod("Dec "+str(i)+"? ")
h = stod("H "+str(i)+"? ")
A.append([Cos(theta)*Cos(dec),Sin(theta)*Cos(dec),Sin(dec)])
b.append(Sin(h))
# If n = 2, the third row represents the DR
if n == 2:
drlat = stod("DR Lat? ")
drlon = stod("DR Lon? ")
A.append([Cos(drlon)*Cos(drlat),Sin(drlon)*Cos(drlat),Sin(drlat)])
b.append(1)
# If n > 3, use the normal equations for an over-determined system
if n > 3:
AT = mtranspose(A)
A = mmult(AT,A)
b = mvmult(AT,b)
A = maugment(A,b)
# Solve Ax = b by Gaussian elimination (row reduction)
M = msolve(A)
x = M[0][3]
y = M[1][3]
z = M[2][3]
r = sqrt(x*x+y*y+z*z)
lat = aSin(z/r)
lon = aTan2(y,x)
print("")
print("lat = "+dtos(lat))
print("lon = "+dtos(lon))