Connect with us

satellite program

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

  1. 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

    ---------------------------------------------------------
     
Ask a Question
Want to reply to this thread or ask your own question?
You'll need to choose a username for the site, which only take a couple of moments (here). After that, you can post your question and our members will help you out.
Electronics Point Logo
Continue to site
Quote of the day

-