- From: david.williams@xxxxxxxxxx (David Williams)
- Date: Tue, 09 Sep 2008 22:51:42 -0500
This program calculates various things concerned with geostationary
satellites such as the ones that are used for TV transmission. It
calculates the position of any satellite in the sky, as seen from any
place on the earth, the rotation of the polarization of the satellite's
signals, and the dates and times when the satellite passes in front of
the sun in the sky, when radio waves from the sun can interfere with
signals from the satellite.
I wrote the "position" part of the program in the early 1980s on a
Commodore PET. It became quite widely used by people who were working
in the then-new field of home satellite receivers. The polarization and
solar outage routines have been added more recently.
It should still be useful to someone...
' 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
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
INPUT #F%, FA, FG, FS, FZ, FM
ON ERROR GOTO 0
LA = FA
LG = FG
LS = FS
TZ = FZ
MD = FM
PRINT "Do you want instructions";
IF YesNo$ = "y" THEN CALL Instructions
' Input section
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
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"
A$ = ROff$(A1 / CF, 1)
AL = VAL(A$)
IF AL = 90 THEN
PRINT "Satellite vertically overhead"
PRINT "Horizontal: East - West"
PRINT "Vertical: North - South"
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 "
PRINT "Satellite's Position:"
PRINT "Altitude (elevation):"; A$; "degrees"
PRINT "Bearing (azimuth):"; B$; "degrees "; CompBear$(BE)
IF AL > 85 THEN
P = B1
CALL AzEl(LD + .001, LR, B2, A2)
P = ATan(A1 - A2, COS(A1) * Force(B2 - B1, PI + PI))
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"
' Solar outage calculation
CALL SlrOuts(LD, LR, LG, TZ)
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)
ON ERROR GOTO 0
E% = 1
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
X = ATN(A / B)
IF B < 0 THEN X = X + PI
ATan = X
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)
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$ + ")"
Z = ABS(180 - BE)
IF Z < 90 THEN
Y1$ = "(South,"
Y1$ = "(North,"
Z = 180 - Z
IF X < 2 THEN
Y2$ = "East)"
Y2$ = "West)"
CompBear$ = Y1$ + ROff$(Z, 1) + "degrees " + Y2$
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$
S$ = S$ + D$
Dirn = S$ + "."
FUNCTION Force (X, L)
'Subtracts nearest multiple of L, so Force is in range +/- L/2
Force = X - L * INT(X / L + .5)
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).
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)"
W = VAL(In$)
IF ABS(W) <= MX% AND INSTR(In$, ",") = 0 THEN
V = W
M$ = LTRIM$(STR$(MX%))
PRINT "Input illegal or out of range! (-"; M$; " to "; M$; ")"
PRINT "Try again..."
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."
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."
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."
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
T$ = "." + RIGHT$(S$, F%)
IF P% <= 0 THEN
DO WHILE RIGHT$(T$, 1) = "0"
T$ = LEFT$(T$, LEN(T$) - 1)
IF T$ = "." THEN T$ = ""
ROff$ = M$ + LEFT$(S$, L% - F%) + T$ + " "
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
' 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
L = VAL(MID$(" 28 31122 31 30999", P, 3))
IF D% <= L THEN EXIT DO
D% = D% - L
P = P + 3
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."
PRINT MN$; STR$(D%); SX$; TAB(20); ZNum$(H%); ":"; ZNum$(M%)
PRINT "Outages occur for a few days before and after nominal dates,"
PRINT "and for a few minutes before and after nominal times."
' Waits for spacebar to be pressed, and clears screen
PRINT "Press Space Bar to continue";
DO WHILE LEN(INKEY$)
DO WHILE INKEY$ <> " "
' Waits for "y" or "n" key to be pressed, and puts lowercase value
' into output.
PRINT "? (y/n) ";
DO WHILE INKEY$ <> ""
G$ = LCASE$(INKEY$)
LOOP UNTIL G$ = "y" OR G$ = "n"
YesNo$ = G$
FUNCTION ZNum$ (X%)
ZNum$ = RIGHT$(STR$(100 + X%), 2)
- Re: satellite calculator
- From: Warren Gay
- Re: satellite calculator
- Prev by Date: great circle program
- Next by Date: Re: heliostat and sun tracker program in qbasic
- Previous by thread: great circle program
- Next by thread: Re: satellite calculator