D
David Williams
- Jan 1, 1970
- 0
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
' Canada
' 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
CONST R = 6.615' Geosynchronous orbit radius / earth's radius
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 "Satellite vertically overhead"
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
' radians
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
---------------------------------------------------------
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
' Canada
' 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
CONST R = 6.615' Geosynchronous orbit radius / earth's radius
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 "*** Use negative numbers for angles in ";
PRINT "opposite directions to those shown. ***"
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$)
IF AL = 90 THEN
PRINT "Satellite vertically overhead"
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)
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 "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
' radians
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)"
EXIT DO
END IF
W = VAL(In$)
IF ABS(W) <= MX% AND INSTR(In$, ",") = 0 THEN
V = W
EXIT DO
END IF
BEEP
M$ = LTRIM$(STR$(MX%))
PRINT "Input illegal or out of range! (-"; M$; " to "; M$; ")"
PRINT "Try again..."
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 "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 "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 "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 "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 "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 "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 "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 "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$
YesNo$ = G$
END FUNCTION
FUNCTION ZNum$ (X%)
ZNum$ = RIGHT$(STR$(100 + X%), 2)
END FUNCTION
---------------------------------------------------------
