satellite calculator
- 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...
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
----------------------------------------------------------
.
- Follow-Ups:
- 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
- Index(es):
Relevant Pages
|
