# satellite program

Discussion in 'Home Power and Microgeneration' started by David Williams, Sep 19, 2007.

1. ### David WilliamsGuest

Off-grid people often depend on satellite communications, so I'll post
the following BASIC (this version is QBasic) program that does some
relevant calculations. It finds:

A) The position of the satellite in the sky, as compass bearing and
angle of elevation, so you know where to point your dish.

B) The rotation of the plane of polarization of the satellite's
signals, so you can set up the receiver properly,

and

C) The dates and times when the satellite passes in front of, or very
close to, the sun in the sky. When this happens, radio waves from the
sun get received by the dish, degrading or even blocking out the signal
feom the satellite. These events are called "solar outages".

I hope someone finds it useful.

dow

-------------------------------------------------------

' SAT_CALC.BAS
' TV Satellites, Commodore PET version, David Williams, 1982

' David Williams
' P.O. Box 48512
' 3605 Lakeshore Blvd. West
' Toronto, Ontario, M8W 4Y6

' Updated for other computers 1995, 2000, 2002
' Solar outage code added March 2003
' Polarization code added February 2007

' This version dated 2007 Feb 12

DECLARE SUB InNum (Pt\$, V!, MX%)
DECLARE FUNCTION YesNo\$ ()
DECLARE FUNCTION ROff\$ (J!, P%)
DECLARE FUNCTION ATan (A!, B!)
DECLARE FUNCTION CompBear\$ (BE!)
DECLARE SUB Instructions ()
DECLARE SUB SpaceBar ()
DECLARE SUB SlrOuts (LD!, LR!, LG!, TZ!)
DECLARE SUB AzEl (LD!, LR!, B!, A!)
DECLARE FUNCTION ZNum\$ (X%)
DECLARE FUNCTION Force (X!, L!)
DECLARE FUNCTION Dirn\$ (X!, D\$, E\$)

CONST PI = 3.14159265#
CONST CF = PI / 180 ' Degree/radian factor

ON ERROR GOTO Etrap
F% = FREEFILE
N\$ = "TVSATDAT.DAT"
E% = 0
OPEN N\$ FOR INPUT AS F%
' Note: File does not have to exist for program to run. If no file
' is found, default values are used. File is created or updated
' at end of program.

IF E% THEN

FA = 45 ')
FG = 80 ') Default latitude & longitudes, used if file not found
FS = 90 ')
FZ = -5 ' Default time zone
FM = 4 ' Default magnetic deviation

ELSE

INPUT #F%, FA, FG, FS, FZ, FM

END IF

CLOSE F%

ON ERROR GOTO 0

LA = FA
LG = FG
LS = FS
TZ = FZ
MD = FM

CLS
PRINT "Do you want instructions";
IF YesNo\$ = "y" THEN CALL Instructions

Newcalc:

' Input section

CLS
PRINT "For your current position on the ground:"
PRINT "Latitude is"; Dirn\$(LA, "south", "north")
PRINT "Longitude is"; Dirn\$(LG, "east", "west")

PRINT
PRINT "*** Use negative numbers for angles in ";
PRINT "opposite directions to those shown. ***"
PRINT

InNum "Your latitude (deg. north)", LA, 90
InNum "Your longitude (deg. west)", LG, 180

PRINT "For current satellite: ";
PRINT "Longitude is"; Dirn\$(LS, "east", "west")
InNum "Satellite's longitude (deg. west)", LS, 180

PRINT "Current magnetic deviation is"; Dirn\$(MD, "east", "west")
PRINT "Input zero to get true bearings."
InNum "Magnetic deviation (deg. west)", MD, 20

Z\$ = ROff\$(TZ, -1) + "hour"
IF VAL(Z\$) > 0 THEN Z\$ = " +" + MID\$(Z\$, 2)
IF ABS(VAL(Z\$)) <> 1 THEN Z\$ = Z\$ + "s"
PRINT "Current time zone offset is"; Z\$; " from GMT / UT."
InNum "Time zone (+/- hours from GMT / UT)", TZ, 14

PRINT "Confirm these entries";
IF YesNo\$ = "n" GOTO Newcalc

CLS

LD = (LG - LS) * CF ' longitude difference in radians
LR = CF * LA ' latitude in radians

' Azimuth-Elevation calculation
CALL AzEl(LD, LR, B1, A1)

IF A1 < 0 THEN

PRINT "Satellite invisible, below horizon"

ELSE

A\$ = ROff\$(A1 / CF, 1)
AL = VAL(A\$)

PRINT

IF AL = 90 THEN

PRINT
PRINT "Polarization:"
PRINT "Horizontal: East - West"
PRINT "Vertical: North - South"

ELSE

BE = B1 / CF + MD
BE = BE - 360 * INT(BE / 360)
B\$ = ROff\$(BE, 1)
BE = VAL(B\$)
IF BE = 360 THEN
BE = 0
B\$ = " 0.0 "
END IF
PRINT "Satellite's Position:"
PRINT "Altitude (elevation):"; A\$; "degrees"
PRINT "Bearing (azimuth):"; B\$; "degrees "; CompBear\$(BE)

PRINT

IF AL > 85 THEN
P = B1
ELSE
CALL AzEl(LD + .001, LR, B2, A2)
P = ATan(A1 - A2, COS(A1) * Force(B2 - B1, PI + PI))
END IF
P = INT(Force(P, PI) / CF + .5)
PRINT "Polarization Rotation:"; ABS(P); "degree";
IF ABS(P) <> 1 THEN PRINT "s";
SELECT CASE P
CASE 0, 90, -90: PRINT
CASE IS > 0: PRINT " clockwise"
CASE ELSE: PRINT " anticlockwise"
END SELECT

END IF

' Solar outage calculation
CALL SlrOuts(LD, LR, LG, TZ)

END IF

PRINT

PRINT "Another calculation";

IF YesNo\$ = "y" GOTO Newcalc

IF LA <> FA OR LG <> FG OR LS <> FS OR TZ <> FZ OR MD <> FM THEN
PRINT "Keep current parameters for next run";
IF YesNo\$ = "y" THEN
ON ERROR GOTO Etrap
E% = 0
OPEN N\$ FOR OUTPUT AS F%
IF E% = 0 THEN
PRINT #F%, ROff\$(LA, -3); ","; ROff\$(LG, -3); ",";
PRINT #F%, ROff\$(LS, -3); ","; ROff\$(TZ, -1); ",";
PRINT #F%, ROff\$(MD, -2)
END IF
CLOSE F%
ON ERROR GOTO 0
END IF
END IF

END

Etrap:
E% = 1
RESUME NEXT

FUNCTION ATan (A, B)
' Calculates angle from X axis to vector from origin to point (B,A)
IF B = 0 THEN
X = SGN(A) * PI / 2
ELSE
X = ATN(A / B)
IF B < 0 THEN X = X + PI
END IF
ATan = X
END FUNCTION

SUB AzEl (LD, LR, BE, AL)
' Calculates azimuth ("bearing") and elevation ("altitude") of
' geostationary satellite. Variables: LD = longitude difference in
' radians between observer and satellite. Positive means observer
' is west of satellite. LR = observer's latitude in radians. BE and
' AL are output variables, satellite's bearing and altitude, in

LT = LR + PI / 2 ' latitude from s. pole

' satellite's x,y,z coordinates
X = R * SIN(LD)
'Y = 0 ' equatorial orbit
Z = R * COS(LD)

' rotate system to put observer at s. pole, i.e. at (0,-1,0)
D = ABS(Z) ' satellite's distance from x-axis
AN = SGN(Z) * PI / 2 'azimuth angle onto y-z plane (y-axis as zero)
AR = AN + LT ' rotate system
Y = D * COS(AR) ' satellite's new y
Z = D * SIN(AR) ' satellite's new z

' calculate altitude

AL = ATan(-1 - Y, SQR(X * X + Z * Z))

' calculate new bearing

BE = ATan(X, Z)

END SUB

FUNCTION CompBear\$ (BE)
' Expresses azimuth bearing in compass format

X = BE / 90
IF X = INT(X) THEN
Y1\$ = RTRIM\$(MID\$("NorthEast SouthWest", 5 * X + 1, 5))
CompBear\$ = "(Due " + Y1\$ + ")"
ELSE
Z = ABS(180 - BE)
IF Z < 90 THEN
Y1\$ = "(South,"
ELSE
Y1\$ = "(North,"
Z = 180 - Z
END IF
IF X < 2 THEN
Y2\$ = "East)"
ELSE
Y2\$ = "West)"
END IF
CompBear\$ = Y1\$ + ROff\$(Z, 1) + "degrees " + Y2\$
END IF
END FUNCTION

FUNCTION Dirn\$ (X, D\$, E\$)
T\$ = ROff\$(ABS(X), -2)
S\$ = T\$ + "degree"
IF T\$ <> " 1 " THEN S\$ = S\$ + "s"
IF T\$ <> " 0 " THEN
S\$ = S\$ + " "
IF X > 0 THEN
S\$ = S\$ + E\$
ELSE
S\$ = S\$ + D\$
END IF
END IF
Dirn = S\$ + "."
END FUNCTION

FUNCTION Force (X, L)
'Subtracts nearest multiple of L, so Force is in range +/- L/2

Force = X - L * INT(X / L + .5)

END FUNCTION

SUB InNum (Pt\$, V, MX%)
' Inputs a number, with various bells and whistles. Input variables:
' Pt\$ is prompt. MX% is maximum allowed absolute value of number. V
' is an input variable, holding the previous value of the quantity,
' and an output variable, holding the new value (if different).

DO

PRINT Pt\$; " (or ENTER to keep value)? ";
C = POS(0)
LINE INPUT In\$
IF In\$ = "" THEN
LOCATE CSRLIN - 1, C - 1
PRINT ROff\$(V, -3); "(kept)"
PRINT
EXIT DO
END IF
W = VAL(In\$)
IF ABS(W) <= MX% AND INSTR(In\$, ",") = 0 THEN
V = W
PRINT
EXIT DO
END IF
BEEP
M\$ = LTRIM\$(STR\$(MX%))
PRINT "Input illegal or out of range! (-"; M\$; " to "; M\$; ")"
PRINT "Try again..."
PRINT

LOOP

END SUB

SUB Instructions

CLS
PRINT "This program calculates the position in the sky, as"
PRINT "compass bearing and altitude (or angle of elevation), of"
PRINT "any satellite that is in geostationary orbit. (Almost all"
PRINT "T.V. broadcasting and relay satellites are geostationary.)"
PRINT
PRINT "The rotation of the plane of polarization of the"
PRINT "satellite's signal is also shown. The sense of the"
PRINT "rotation is calculated looking toward the satellite"
PRINT "from the ground."
PRINT
PRINT "The program also calculates the dates and times of 'solar"
PRINT "outages' which occur when the satellite passes in front"
PRINT "of, or very near, the sun in the sky. Radio waves from the"
PRINT "sun then interfere with the satellite's communications."

CALL SpaceBar

PRINT "You will be asked for your latitude and longitude, and for"
PRINT "the longitude of the satellite. Enter these quantities, in"
PRINT "degrees, accurate to at least one place of decimals (0.1"
PRINT "degree) if possible. Errors greater than 0.1 degree will"
PRINT "cause significantly inaccurate calculated results. It is"
PRINT "a good idea to use a G.P.S. receiver to find your own"
PRINT "latitude and longitude."
PRINT
PRINT "You will also be asked for your time zone, which is used"
PRINT "in the solar outage calculation. Enter the number of hours"
PRINT "by which it is ahead of (+) or behind (-) G.M.T. / U.T."
PRINT
PRINT "You will be asked for the local magnetic deviation, which"
PRINT "will be used in calculating the compass bearing of the"
PRINT "satellite. If you want the true (i.e. not magnetic) bearing,"
PRINT "input a value of zero."

CALL SpaceBar

PRINT "The solar outage dates and times shown by the program"
PRINT "are those when the satellite is closest to the centre of"
PRINT "the sun's disk, as seen from your location on the earth."
PRINT "Outages can occur when the satellite is a couple of degrees"
PRINT "from the centre, depending on sunspot activity and the"
PRINT "geometry of your dish. This means that outages can be"
PRINT "several minutes in length, centred on the nominal time, and"
PRINT "can occur for several successive days (at the same time of"
PRINT "day), centred on the nominal date."
PRINT
PRINT "Your entries of latitude, longitudes, time zone and"
PRINT "magnetic deviation can be kept on disk and used in subsequent"
PRINT "runs. Initially, arbitary defaults are used."

CALL SpaceBar

END SUB

FUNCTION ROff\$ (J, P%)
' String showing number J rounded to ABS(P%) places of decimals.
' Trailing decimal zeroes are truncated if P% is negative.
' Not if not. Leading and trailing spaces are included.

F% = ABS(P%)
H& = INT(10 ^ F% * J + .5)
IF H& < 0 THEN M\$ = " -" ELSE M\$ = " "
S\$ = LTRIM\$(STR\$(ABS(H&)))
L% = LEN(S\$)
IF L% < F% + 1 THEN
S\$ = STRING\$(F% + 1 - L%, "0") + S\$
L% = F% + 1
END IF
T\$ = "." + RIGHT\$(S\$, F%)
IF P% <= 0 THEN
DO WHILE RIGHT\$(T\$, 1) = "0"
T\$ = LEFT\$(T\$, LEN(T\$) - 1)
LOOP
IF T\$ = "." THEN T\$ = ""
END IF
ROff\$ = M\$ + LEFT\$(S\$, L% - F%) + T\$ + " "
END FUNCTION

SUB SlrOuts (LD, LR, LG, TZ)
' Calculates and prints dates and times of "solar outages", when a
' (geostationary) satellite passes in front of the sun in the sky.
' Input variables: LD = longitude difference in radians. Positive
' means observer is west of satellite. LR = observer's latitude, in
' radians. LG = observer's longitude west of Greenwich, in degrees.
' TZ = time zone offset in hours from GMT/UT. The global (CONSTant)
' quantities R (orbit radius), PI and CF (PI/180) are also used.

W = PI / 182.5 ' mean orbital angular velocity
W1 = 12 * W

' coordinate differences
X = -R * SIN(LD)
Y = -SIN(LR)
Z = R * COS(LD) - COS(LR)

' calculate -sin(apparent declination)/sin(axial tilt)
C = Y / (.3979 * SQR(X * X + Y * Y + Z * Z))
AC = PI / 2 - ATan(C, SQR(1 - C * C))' arc-cosine for decl'n calc

' calculate apparent longitude
AG = ATan(X, Z)

V = 4 * (AG / CF + LG) + 60 * TZ + 720.5

PRINT
PRINT "Nominal solar outage dates and times: (Times are HH:MM.)"

' find dates (days after Dec solstice) when sun is at same decl'n
FOR J = 0 TO 1

U = 182
L = 364 * J
FOR K = 1 TO 20
D = (U + L) / 2
G = W * D
A = G + .0334 * SIN(G - W1)
IF ABS(U - L) < .1 THEN EXIT FOR
IF ABS(A - PI) < AC THEN U = D ELSE L = D
NEXT

' calculate equation of time
B = (ATN(TAN(A) / .9174) - G) / PI
ET = (B - INT(B + .5)) * 720

' calculate outage time
T = INT(V + ET)
T = T - 1440 * INT(T / 1440)
H% = T \ 60
M% = T MOD 60

' convert date to month and day
D = D - 40 + TZ / 24 ' 40 days from Dec solstice to Jan 31
D% = INT(D)
D% = D% + INT(D - D% - T / 1440 + .5)

P = 1
DO
L = VAL(MID\$(" 28 31122 31 30999", P, 3))
IF D% <= L THEN EXIT DO
D% = D% - L
P = P + 3
LOOP
MN\$ = "February March April August SeptemberOctober"
MN\$ = RTRIM\$(MID\$(MN\$, 3 * P - 2, 9))

SELECT CASE D%
CASE 1, 21, 31: SX\$ = "st."
CASE 2, 22: SX\$ = "nd."
CASE 3, 23: SX\$ = "rd."
CASE ELSE: SX\$ = "th."
END SELECT

PRINT MN\$; STR\$(D%); SX\$; TAB(20); ZNum\$(H%); ":"; ZNum\$(M%)

NEXT

PRINT
PRINT "Outages occur for a few days before and after nominal dates,"
PRINT "and for a few minutes before and after nominal times."

END SUB

SUB SpaceBar
' Waits for spacebar to be pressed, and clears screen

PRINT
PRINT "Press Space Bar to continue";

DO WHILE LEN(INKEY\$)
LOOP

DO WHILE INKEY\$ <> " "
LOOP

CLS

END SUB

FUNCTION YesNo\$
' Waits for "y" or "n" key to be pressed, and puts lowercase value
' into output.

PRINT "? (y/n) ";
DO WHILE INKEY\$ <> ""
LOOP
DO
G\$ = LCASE\$(INKEY\$)
LOOP UNTIL G\$ = "y" OR G\$ = "n"
PRINT G\$
PRINT
YesNo\$ = G\$

END FUNCTION

FUNCTION ZNum\$ (X%)
ZNum\$ = RIGHT\$(STR\$(100 + X%), 2)
END FUNCTION

---------------------------------------------------------  