CNDEWIT - Calculate a least squares position using DeWit's method
Description: This least squares method uses intercepts to produce a position from three or more sights (4 or more are recommended for increased accuracy). Each sight can be assigned a weight according to confidence in the observation. If the weights are equal, a vanilla least squares solution results. For running fixes, the distance made good (DMG) and course made good (CMG) can be specified for each leg between sights.
Inputs:
| Number of sights | The number of sights to be reduced |
| DR Lat | The DR latitude of the observer |
| DR Lon | The DR longitude of the observer |
| 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 |
| SD i | Standard deviation assigned to the ith observation |
| DMG i->i+1 | Distance made good between sights i and i+1 |
| CMG i->i+1 | Course made good between sights i and i+1 (DMG > 0) |
| ... 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 DR Lat? 49,55 DR Lon? -4,41 GHA 1? 126,13.1 Dec 1? 46,01.3 H 1? 18,41.5 SD 1? 1 DMG 1->2? 0 GHA 2? 28,15.3 Dec 2? 14,26.2 H 2? 49,45.2 SD 2? 1 DMG 2->3? 0 GHA 3? 4,12.6 Dec 3? -11,17.4 H 3? 28,44.4 SD 3? 1 DMG 3->4? 0 GHA 4? 286,23.1 Dec 4? 38,48.2 H 4? 35,24.5 SD 4? 1 Lat = 49° 59.5'(SD = 0.7 NM) Lon = -5° 00.4'(SD = 0.8 NM)
(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/.
# CNDEWIT - Multiple sights weighted least squares solution using DeWit's method
from math import *
from CN_LIB import *
n = int(input("Number of sights? "))
lat = drlat = stod("DR Lat? ")
lon = drlon = stod("DR Lon? ")
Ai = [] # Azimuths
Zi = [] # Intercepts
Si = [] # Weights (standard deviations)
for i in range(1,n+1):
gha = stod("GHA "+str(i)+"? ")
dec = stod("Dec "+str(i)+"? ")
ho = stod("H "+str(i)+"? ")
s = float(input("SD "+str(i)+"? "))
lha = (gha+lon+360)%360
hc = aSin(Sin(lat)*Sin(dec)+Cos(lat)*Cos(dec)*Cos(lha))
z = aCos((Sin(dec)-Sin(hc)*Sin(lat))/(Cos(hc)*Cos(lat)))
zn = 360-z if lha <= 180 else z
Ai.append(zn)
Zi.append(ho-hc)
Si.append(s)
if (i < n):
dmg = float(input("DMG "+str(i)+"->"+str(i+1)+"? "))
if (dmg > 0):
cmg = stod("CMG "+str(i)+"->"+str(i+1)+"? ")
lat = lat + dmg*Cos(cmg)/60
lon = lon + dmg*Sin(cmg)/Cos(lat)/60
a11 = 0
a12 = 0
a22 = 0
b1 = 0
b2 = 0
for i in range(len(Ai)):
a11+= Cos(Ai[i])*Cos(Ai[i])/(Si[i]*Si[i])
a12+= Cos(Ai[i])*Sin(Ai[i])/(Si[i]*Si[i])
a22+= Sin(Ai[i])*Sin(Ai[i])/(Si[i]*Si[i])
b1+= Zi[i]*Cos(Ai[i])/(Si[i]*Si[i])
b2+= Zi[i]*Sin(Ai[i])/(Si[i]*Si[i])
A = a11*a22-a12*a12
x1 = (b1*a22-b2*a12)/A
x2 = (b2*a11-b1*a12)/A
lat = lat+x1
lon = lon+x2/Cos(lat)
sd1 = sqrt(a22/A)
sd2 = sqrt(a11/A)
print("")
print("Lat = "+dtos(lat)+"(SD = "+"%.1f NM)" % (sd1))
print("Lon = "+dtos(lon)+"(SD = "+"%.1f NM)" % (sd2))