'===========================================================================
' Subject: CALCULATE TV SATELLITE POSITION   Date: 08-08-00 (16:45)       
' Author:  David Williams                    Code: QB, QBasic, PDS        
' Origin:  david.williams@capcanada.com    Packet: ALGOR.ABC
'===========================================================================
' TV Satellites, Commodore PET version, David Williams, 1982
' Updated for other computers, without major changes, 1995, 2000

DECLARE FUNCTION InNum! (MX%, ER%)
DECLARE FUNCTION YesNo$ ()
DECLARE SUB Instructions ()
DECLARE SUB Spacebar ()
                       
DEFSNG A-Z

CONST R = 6.6169 ' Radius of geosynchronous orbit in units of earth radius

PI = 4 * ATN(1)
CF = PI / 180  ' Degree/radian factor

LA = 43.65
LO = 79.33  ' Latitude & longitude of Toronto (default values)

DEF FNRO (X) = INT(100 * X + .5) / 100

' Read satellite data from table at end of program.

N% = -1
DO
   READ SN$
   IF SN$ = "ENDDATA" THEN EXIT DO
   READ I
   N% = N% + 1
LOOP

IF N% > -1 THEN
   DIM SN$(N%), LS(N%)
   RESTORE
   FOR I% = 0 TO N%
      READ SN$(I%), LS(I%)
   NEXT
END IF

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

' Input section

PRINT
PRINT "Default latitude is "; LA; " degrees north"
PRINT "Default longitude is "; LO; " degrees west"
PRINT "Keep these";

IF YesNo$ = "n" THEN
  Inflag% = -1
ELSE
  Inflag% = 0
END IF

Newcoords:

IF Inflag% = -1 THEN
 
    PRINT
    PRINT "Input latitudes & longitudes in degrees."
    PRINT "Use negative numbers for angles in"
    PRINT "opposite directions to those shown."
   
    DO
      PRINT
      PRINT "Your latitude (deg. north)";
      LA = InNum(90, ER%)
    LOOP UNTIL ER% = 0

    DO
      PRINT
      PRINT "Your longitude (deg. west)";
      LO = InNum(360, ER%)
    LOOP UNTIL ER% = 0

END IF

Newsat:

LT = LA
LG = LO

Flag% = 0

IF N% > -1 THEN

  PRINT
  PRINT "List of pre-programmed satellites:"
  PRINT

  FOR I% = 0 TO N%
    PRINT CHR$(97 + I%); ".  "; SN$(I%)
  NEXT

  PRINT
  PRINT "Press letter to left of name of"
  PRINT "desired satellite, or ENTER if"
  PRINT "it is not on list."
  PRINT

  DO
    DO
      G$ = LCASE$(INKEY$)
    LOOP UNTIL G$ <> ""
    G% = ASC(G$) - 97
  LOOP UNTIL G$ = CHR$(13) OR (G% >= 0 AND G% <= N%)

  IF G$ <> CHR$(13) THEN
    
     PRINT "Satellite selected: "; SN$(G%);
     LS = LS(G%)
     Flag% = -1
 
  END IF

END IF
    
IF Flag% = 0 THEN

     DO
       PRINT
       PRINT "Input longitude of satellite in"
       PRINT "degrees west. (Use negative number"
       PRINT "if you want to input degrees east.)"
       PRINT "Satellite's longitude ";
       LS = InNum(360, ER%)
     LOOP UNTIL ER% = 0
    
     PRINT
     PRINT "New satellites may be included in"
     PRINT "program by adding them to data table"
     PRINT "at end."
 
END IF
   
LD = LG - LS ' longitude difference
LT = 90 + LT ' latitude from s. pole
LT = CF * LT ' latitude in radians
LD = CF * LD ' longitude diff. in radians

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

' rotate system to put observer at s. pole
D = ABS(Z) ' distance from x-axis
AN = PI * ((Z < 0) + .5) ' azimuth angle onto y-z plane (y-axis as zero)
AN = AN + LT ' rotate system
Y = D * COS(AN) ' new y
Z = D * SIN(AN) ' new z

' calculate new bearings
IF Z = 0 THEN
  BE = PI * ((X < 0) + .5)
ELSE
  BE = ATN(X / Z) - PI * (Z < 0) ' bearing
END IF
 
IF X = 0 AND Z = 0 THEN
  AL = (PI / 2) * SGN(-1 - Y)
ELSE
  AL = ATN((-1 - Y) / (SQR(X * X + Z * Z))) ' altitude
END IF

BE = BE / CF ' bearing in degrees
AL = AL / CF ' alt. in degrees
 
DO WHILE BE < 0
  BE = BE + 360
LOOP

DO WHILE BE >= 360
  BE = BE - 360
LOOP

' roundoff
BE = FNRO(BE)
AL = FNRO(AL)

' printout
PRINT
 
IF AL < 0 THEN
   
  PRINT "Satellite invisible (below horizon)"
 
ELSE
   
  PRINT
  PRINT "Satellite's position:"
  PRINT
 
  PRINT "Bearing is"; BE; "degrees."
   
  IF BE = 0 THEN
    PRINT "(due north)"
  ELSEIF BE < 90 THEN
    PRINT "(north"; BE; "degrees east)"
  ELSEIF BE = 90 THEN
    PRINT "(due east)"
  ELSEIF BE < 180 THEN
    PRINT "(south"; 180 - BE; "degrees east)"
  ELSEIF BE = 180 THEN
    PRINT "(due south)"
  ELSEIF BE < 270 THEN
    PRINT "(south"; BE - 180; "degrees west)"
  ELSEIF BE = 270 THEN
    PRINT "(due west)"
  ELSE
    PRINT "(north"; 360 - BE; "degrees west)"
  END IF

  PRINT
  PRINT "Altitude is"; AL; "degrees."
  IF AL < 3 THEN
    PRINT
    PRINT "Very low altitude -"
    PRINT "Communication unreliable"
  END IF
 
END IF

PRINT
 
Calcloop:
 
  PRINT "Another calculation";
 
  IF YesNo$ = "n" THEN END
 
  PRINT
  PRINT "Another satellite (same latitude and longitude)";
 
  IF YesNo$ = "y" GOTO Newsat
 
  PRINT
  PRINT "Another latitude and/or longitude";
 
  IF YesNo$ = "y" THEN
    Inflag% = -1
    GOTO Newcoords
  ELSE
    Inflag% = 0
  END IF

GOTO Calcloop

END

' satellite data table

DATA ANIK 1,104
DATA ANIK B,109
DATA ANIK A2 & A3,114
DATA ANIK D1,104.5
DATA WESTAR 1,99
DATA WESTAR 2,123.5
DATA WESTAR 3,91
DATA SATCOM 2,119
DATA SATCOM 4,83
DATA COMSTAR D1 & D2,95

' Additional satellites, and their longitudes (deg. west) can be added
' to table above. No other changes to program are needed.

DATA ENDDATA

FUNCTION InNum (MX%, ER%)

   ER% = 0
   INPUT In$
   IF In$ = "" THEN ER% = -1
   IF ER% = 0 THEN
     In = VAL(In$)
     ER% = (ABS(In) > MX%)
   END IF
   IF ER% THEN
      PRINT
      BEEP
      PRINT "Input illegal or out of range!"
      PRINT "Try again..."
      PRINT
      In = 0
   END IF
   InNum = In

END FUNCTION

SUB Instructions

   CLS
   PRINT "This program calculates the position"
   PRINT "(expressed as true compass bearing and"
   PRINT "altitude) of any satellite which is in"
   PRINT "geostationary orbit. (Almost all T.V."
   PRINT "broadcasting and relay satellites are"
   PRINT "geostationary.) You will be asked for"
   PRINT "your latitude and longitude, and for"
   PRINT "the longitude of the satellite."
   PRINT
   PRINT "The program has default values for"
   PRINT "your latitude and longitude, which"
   PRINT "are those of Toronto. If you are in"
   PRINT "Toronto, just keep these values."
   PRINT "It also contains longitude data for"
   PRINT "some frequently used satellites."

   CALL Spacebar

   CLS
   PRINT "You can change the default latitude"
   PRINT "and longitude by editing the info"
   PRINT "near the start of the program."
   PRINT "The satellite data can be found at"
   PRINT "the end of the program."
   
   CALL Spacebar

   PRINT
   PRINT

END SUB

SUB Spacebar
 
  PRINT
  PRINT "Press Space Bar to continue";

  DO WHILE INKEY$ <> " "
  LOOP

END SUB

FUNCTION YesNo$
  
   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