' geostationary satellite pointing
' by TassyJim 27 Nov 2012
' based on Nomads Alamnac www.c-com.com.au/nomad.htm
' formula from an old edition of the ARRL Handbook

satC1Long=156 ' position of the satellite of interest
fromFront=1 ' for skew, we are looking at the front of the dish
dim magdev$(36) ' magnetic variation table storage
' 46 degrees of longitude at 5 bytes each gives a string 230 bytes long.
' one string for each of 36 degrees of latitude
open "magdev.txt" for input as #1
for n = 1 to 36
input #1,magdev$(n)
next n

input "Latitude ";lat1
input "Longitude ";lon1
if lat1=0 then lat1= -41.1258 ' default location
if lon1=0 then lon1= 145.9078 ' we could be getting this from a GPS
cls
if mm.device$ <> "" then font 2,1 ' don't try to set the font for DOS MMBasic
print " Optus C1 bearing from "
print " ";lat1;" ";lon1
magDev=deviation(lat1,lon1)
print
print " Magnetic Variation:";format$(magDev,"% 5.1f")
dishAZ= azimuth(satC1Long,lat1,lon1)
print " Azimuth (True): ";format$(dishAZ,"% 5.1f")

dishAZmag=dishAZ-magDev
if dishAZmag < 0 then dishAZmag=dishAZmag+360
print " Azimuth (Magnetic):"format$(dishAZmag,"% 5.1f")
dishEL= elevation(satC1Long,lat1,lon1)
print " Dish Elevation: ";format$(dishEL,"% 5.1f")
dishSK= tilt(satC1Long,lat1,lon1)
dishSK=dishSK-45 ' Optus satellite skew
if fromFront=0 then
print " LNB Skew: ";format$(dishSK,"% 5.1f")
print
print " Skew is clockwise looking at Sat."
else
print " LNB Skew: ";format$(dishSK* -1,"% 5.1f")
print
print " Skew is clockwise looking at dish"
endif
end

function azimuth(satlong,eslat,eslong)
longdiffr=(eslong-satlong)/57.29578
esazimuth=180+57.29578*atn(tan(longdiffr)/sin((eslat/57.2957 8)))
if eslat<0 then esazimuth=esazimuth-180
if esazimuth<0 then esazimuth=esazimuth+360.0
azimuth = esazimuth
end function

function elevation(satlong,eslat,eslong)
longdiffr=(eslong-satlong)/57.29578
eslatr=eslat/57.29578
r1=1+35786/6378.16
v1=r1*cos(eslatr)*cos(longdiffr)-1
v2=r1*sqr(1 - cos(eslatr)*cos(eslatr)*cos(longdiffr)*cos(longdiffr))
eselevation=57.29578*atn(v1/v2)
if eselevation < 30 then
eselrefracted=(eselevation+sqr(eselevation*eselevation+4.132 ))/2
else
eselrefracted=eselevation
endif
if eselrefracted < 1 then eselrefracted=0
elevation= eselrefracted
end function

function tilt(satlong,eslat,eslong)
longdiffr=(satlong-eslong)/57.29578
eslatr=eslat/57.29578
poltilt=-57.29578*atn(sin(longdiffr)/tan(eslatr))
tilt= poltilt
end function

function deviation(lat,long) ' magnetic variation for Australia
latmin=abs(lat) -9 ' our table starts at 10 deg South
longmin= long-109 ' and 110 deg East
if latmin < 1 or longmin < 1 then ' check that we are in range of the table
deviation = 0 ' we could use a better number to indicate
goto _exitfunc ' out of range. eg 99
endif
if latmin > 36 or longmin > 46 then
deviation = 0
goto _exitfunc
endif
latI=int(latmin)
longI=int(longmin)
dev1=val(mid$(magdev$(latI),longI*5-4,5)) ' data is stored in a string
dev2=val(mid$(magdev$(latI),longI*5+1,5)) ' for each degree of latitude
dev3=val(mid$(magdev$(latI+1),longI*5-4,5)) ' we extract the 4 points bounding
dev4=val(mid$(magdev$(latI+1),longI*5+1,5)) ' our position
inclong=dev2-dev1:inclat=dev3-dev1 ' and interpolate for our position
deviation=dev1+inclong*(longmin-longI)+inclat*(latmin-latI)
_exitfunc:
end function