Celestial Nav Calculator for Scilab

//hs is hight read off sextant

//aries correction sha and dec are values read off a nauitcal almanic such as https://thenauticalalmanac.com/TNARegular/2021_Nautical_Almanac.pdf

//aplat and aplong are Assumed Points for lat and long. Basically this is an guess of where you think you are

//Basically this is what happens. It finds out where the star is in the sky using the SHA aries correction and DEC. Then it calculates what the sextant reading would have been if you were actually on the assumed point. It compares that answer with the actual sextant reading and adjust accordingly creating a line that you are either on or near. If you do this with two or more stars your position should be on or near the intersection of the two lines. Probably accurate to 50 NM.. This script may not work in the Eastern Hemisphere. I have not checked yet


//I broke something lol

hs1 =
aries1 =
correction1 =
sha1 =
dec1 =
aplat1 =
aplong1 =
star1 = aries1 + correction1 + sha1
LHA1 = star1 + aplong1

hs2 =
aries2 =
correction2 =
sha2 =
dec2 =
aplat2 =
aplong2 =
star2 = aries2 + correction2 + sha2
LHA2 = star2 + aplong2


if LHA1 < 360 then
    LHA1 = LHA1 + 360
end
if LHA1 > 360 then
    LHA1 = LHA1 - 360
end
if LHA2 < 360 then
    LHA2 = LHA2 + 360
end
if LHA2 > 360 then
    LHA2 = LHA2 - 360
end

if aplat1 > 0 then
    aplatL1 = "N"
else
    aplatL1 = "S"
end
if aplat2 > 0 then
    aplatL2 = "N"
else
    aplatL2 = "S"
end

//Find the Alt of the star at the AP
if aplat1/abs(aplat1) == dec1/abs(dec1) then
    hc1 = asind(cosd(LHA1)*cosd(aplat1)*cosd(dec1) + sind(aplat1)*sind(dec1))
else
    hc1 = asind(cosd(LHA1)*cosd(aplat1)*cosd(dec1) - sind(aplat1)*sind(dec1))
end
hd1=hs1-hc1

if aplat2/abs(aplat2) == dec2/abs(dec2) then
    hc2 = asind(cosd(LHA2)*cosd(aplat2)*cosd(dec2) + sind(aplat2)*sind(dec2))
else
    hc2 = asind(cosd(LHA2)*cosd(aplat2)*cosd(dec2) - sind(aplat2)*sind(dec2))
end
hd2=hs2-hc2


//Find the direction to the star at the AP
a1=abs(tand(aplat1)/tand(LHA1))
if LHA1 > 90 & LHA1 < 270 then
    aL1 = aplatL1
elseif aplatL1 == "S" then
    aL1 = "N"
elseif aplatL1 == "N" then
    aL1 = "S"
end
a2=abs(tand(aplat2)/tand(LHA2))
if LHA2 > 90 & LHA2 < 270 then
    aL2 = aplatL2
elseif aplatL2 == "S" then
    aL2 = "N"
elseif aplatL2 == "N" then
    aL2 = "S"
end

b1=abs(tand(dec1)/sind(LHA1))
if dec1 > 0 then
    bL1 = "N"
else
    bL1 = "S"
end
b2=abs(tand(dec2)/sind(LHA2))
if dec2 > 0 then
    bL2 = "N"
else
    bL2 = "S"
end



if bL1 == aL1 then
    c1=a1+b1
elseif bL1 ~= aL1 then
    c1=abs(a1-b1)
end
if a1 > b1 then
    cL1 = aL1
elseif a1 < b1
    cL1 = bL1
end
if bL2 == aL2 then
    c2=a2+b2
elseif bL2 ~= aL2 then
    c2=abs(a2-b2)
end
if a2 > b2 then
    cL2 = aL2
elseif a2 < b2
    cL2 = bL2
end




az1=atand(1/c1/cosd(aplat1))
azVL1 = cL1
if LHA1 < 180 then
    azHL1 = "W"
elseif LHA1 > 180
    azHL1 = "E"
end

if azVL1 == "N" & azHL1 == "E" then
    azT1 = az1
elseif azVL1 == "S" & azHL1 == "E" then
    azT1 = 180 - az1
elseif azVL1 == "S" & azHL1 == "W" then
    azT1 = 180 + az1
elseif azVL1 == "N" & azHL1 == "W" then
    azT1 = 270 + (90 - az1)
end

az2=atand(1/c2/cosd(aplat2))
azVL2 = cL2
if LHA2 < 180 then
    azHL2 = "W"
elseif LHA2 > 180
    azHL2 = "E"
end
if azVL2 == "N" & azHL2 == "E" then
    azT2 = az2
elseif azVL1 == "S" & azHL2 == "E" then
    azT2 = 180 - az2
elseif azVL2 == "S" & azHL2 == "W" then
    azT2 = 180 + az2
elseif azVL2 == "N" & azHL2 == "W" then
    azT2 = 270 + (90 - az2)
end






//Drawing corrections
//Caplat = Corrected aplat

Caplat1 = aplat1 + hd1*cosd(azT1)
Caplong1 = aplong1 + hd1*sind(azT1)


Caplat2 = aplat2 + hd2*cosd(azT2)
Caplong2 = aplong2 + hd2*sind(azT2)


//The Intersection Point
x=Caplong1-Caplong2
y=Caplat1-Caplat2
c=sqrt(x^2+y^2)
theta1=atand(y/x)
theta2=atand(x/y)
phi1=90 - (azT1-90) - (theta1)
phi2= 180 - theta2 - (360-azT2-90)
//phii2 = 90 - (azT2-90) - theta2
angleC = 180 - phi1 - phi2
L1=c*sind(phi1)/sind(angleC)
L2=c*sind(phi2)/sind(angleC)
EPY=L2*sind(90-(azT1-90)) + Caplat1
EPX=L2*cosd(90-(azT1-90)) + Caplong1






//Chart Plotting
Appx1 = Caplong1 + (L2+1)*sind(azT1+90)
Appy1 = Caplat1 + (L2+1)*cosd(azT1+90)
Bppx1 = Caplong1 - (L2+1)*sind(azT1+90)
Bppy1 = Caplat1 - (L2+1)*cosd(azT1+90)
xpts1 = [Appx1 Bppx1];
ypts1 = [Appy1 Bppy1];
Clinex1 = [aplong1 Caplong1];
Cliney1 = [aplat1 Caplat1];

Appx2 = Caplong2 + (L1+1)*sind(azT2+90)
Appy2 = Caplat2 + (L1+1)*cosd(azT2+90)
Bppx2 = Caplong2 - (L1+1)*sind(azT2+90)
Bppy2 = Caplat2 - (L1+1)*cosd(azT2+90)
xpts2 = [Appx2 Bppx2];
ypts2 = [Appy2 Bppy2];
Clinex2 = [aplong2 Caplong2];
Cliney2 = [aplat2 Caplat2];








//Create a new AP on the EPX and EPY and stuff
aplat3 = EPY
aplong3 = EPX
LHA3 = star1 + aplong3

aplat4 = EPY
aplong4 = EPX
LHA4 = star2 + aplong4




if LHA3 < 360 then
    LHA3 = LHA3 + 360
end
if LHA3 > 360 then
    LHA3 = LHA3 - 360
end
if aplat3 > 0 then
    aplatL3 = "N"
else
    aplatL3 = "S"
end

if LHA4 < 360 then
    LHA4 = LHA4 + 360
end
if LHA4 > 360 then
    LHA4 = LHA4 - 360
end
if aplat4 > 0 then
    aplatL4 = "N"
else
    aplatL4 = "S"
end




if aplat3/abs(aplat3) == dec1/abs(dec1) then
    hc3 = asind(cosd(LHA3)*cosd(aplat3)*cosd(dec1) + sind(aplat3)*sind(dec1))
else
    hc3 = asind(cosd(LHA3)*cosd(aplat3)*cosd(dec1) - sind(aplat3)*sind(dec1))
end
hd3=hs1-hc3
if aplat4/abs(aplat4) == dec2/abs(dec2) then
    hc4 = asind(cosd(LHA4)*cosd(aplat4)*cosd(dec2) + sind(aplat4)*sind(dec2))
else
    hc4 = asind(cosd(LHA4)*cosd(aplat4)*cosd(dec2) - sind(aplat4)*sind(dec2))
end
hd4=hs2-hc4



//Gpood
a3=abs(tand(aplat3)/tand(LHA3))
if LHA3 > 90 & LHA3 < 270 then
    aL3 = aplatL3
elseif aplatL3 == "S" then
    aL3 = "N"
elseif aplatL3 == "N" then
    aL3 = "S"
end
a4=abs(tand(aplat4)/tand(LHA4))
if LHA4 > 90 & LHA4 < 270 then
    aL4 = aplatL4
elseif aplatL4 == "S" then
    aL4 = "N"
elseif aplatL4 == "N" then
    aL4 = "S"
end


//Gpood
b3=abs(tand(dec1)/sind(LHA3))
if dec1 > 0 then
    bL3 = "N"
else
    bL3 = "S"
end
b4=abs(tand(dec2)/sind(LHA4))
if dec2 > 0 then
    bL4 = "N"
else
    bL4 = "S"
end



//Gpood

if bL3 == aL3 then
    c3=a3+b3
elseif bL3 ~= aL3 then
    c3=abs(a3-b3)
end
if a3 > b3 then
    cL3 = aL3
elseif a3 < b3
    cL3 = bL3
end
if bL4 == aL4 then
    c4=a4+b4
elseif bL4 ~= aL4 then
    c4=abs(a4-b4)
end
if a4 > b4 then
    cL4 = aL4
elseif a4 < b4
    cL4 = bL4
end


//Gpood
az3=atand(1/c3/cosd(aplat3))
azVL3 = cL3
if LHA3 < 180 then
    azHL3 = "W"
elseif LHA3 > 180
    azHL3 = "E"
end
az4=atand(1/c4/cosd(aplat4))
azVL4 = cL4
if LHA4 < 180 then
    azHL4 = "W"
elseif LHA4 > 180
    azHL4 = "E"
end



if azVL3 == "N" & azHL3 == "E" then
    azT3 = az3
elseif azVL3 == "S" & azHL3 == "E" then
    azT3 = 180 - az3
elseif azVL3 == "S" & azHL3 == "W" then
    azT3 = 180 + az3
elseif azVL3 == "N" & azHL3 == "W" then
    azT3 = 270 + (90 - az3)
end
if azVL4 == "N" & azHL4 == "E" then
    azT4 = az4
elseif azVL4 == "S" & azHL4 == "E" then
    azT4 = 180 - az4
elseif azVL4 == "S" & azHL4 == "W" then
    azT4 = 180 + az4
elseif azVL4 == "N" & azHL4 == "W" then
    azT4 = 270 + (90 - az4)
end




Caplat3 = aplat3 + hd3*cosd(azT3)
Caplong3 = aplong3 + hd3*sind(azT3)

Caplat4 = aplat4 + hd4*cosd(azT4)
Caplong4 = aplong4 + hd4*sind(azT4)

//The Intersection Point 2222
x3=Caplong3-Caplong4
y3=Caplat3-Caplat4
c3=sqrt(x3^2+y3^2)
theta3=atand(y3/x3)
theta4=atand(x3/y3)
phi3=90 - (azT3-90) - (theta3)
phi4= 180 - theta4 - (360-azT4-90)
//phii4 = 90 - (azT4-90) - theta4
angleC3 = 180 - phi3 - phi4
L3=c3*sind(phi3)/sind(angleC3)
L4=c3*sind(phi4)/sind(angleC3)
EPY3=L4*sind(90-(azT3-90)) + Caplat3
EPX3=L4*cosd(90-(azT3-90)) + Caplong3






//Chart Plotting 222222222
Appx3 = Caplong3 + (L4+3)*sind(azT3+90)
Appy3 = Caplat3 + (L4+3)*cosd(azT3+90)
Bppx3 = Caplong3 - (L4+3)*sind(azT3+90)
Bppy3 = Caplat3 - (L4+3)*cosd(azT3+90)
xpts3 = [Appx3 Bppx3];
ypts3 = [Appy3 Bppy3];
Clinex3 = [aplong3 Caplong3];
Cliney3 = [aplat3 Caplat3];

Appx4 = Caplong4 + (L3+3)*sind(azT4+90)
Appy4 = Caplat4 + (L3+3)*cosd(azT4+90)
Bppx4 = Caplong4 - (L3+3)*sind(azT4+90)
Bppy4 = Caplat4 - (L3+3)*cosd(azT4+90)
xpts4 = [Appx4 Bppx4];
ypts4 = [Appy4 Bppy4];
Clinex4 = [aplong4 Caplong4];
Cliney4 = [aplat4 Caplat4];












//Plotting
clf
title('Your Estimated Position:  ' + string(EPX3) + "  " + string(EPY3))



plot(-sha1,dec1,"*Y")
plot(aplong1,aplat1,"*G")
plot(Caplong1,Caplat1,"*R")
plot(xpts1, ypts1)
plot(Clinex1,Cliney1,"G")



plot(-sha2,dec2,"Y--+")
plot(aplong2,aplat2,"G--+")
plot(Caplong2,Caplat2,"R--+")
plot(xpts2, ypts2)
plot(Clinex2,Cliney2, "G")



plot(EPX,EPY,"B--+")
plot(Caplong1+EPX-5,Caplat1+EPY-5,"*W")
plot(Caplong1+EPX+5,Caplat1+EPY+5,"*W")



plot(Caplong3,Caplat3,"o--R")
plot(Caplong4,Caplat4,"o--R")
plot(Clinex3,Cliney3,"G")
plot(Clinex4,Cliney4,"G")
plot(xpts3, ypts3,"R")
plot(xpts4, ypts4,"R")
plot(EPX3,EPY3,"B--+")