'===========================================================================
' Subject: TV Satellite Location              Date: 12-15-02 (  :  )       
'  Author: David Williams                     Code: Qbasic,QB,PDS          
'  Origin: david.williams@ablelink.org      Packet: ALGOR.ABC
'===========================================================================
' TV_SATS.BAS 
' TV Satellites, Commodore PET version, David Williams, 1982 
  
' david.williams@ablelink.org 
  
' Updated for other computers 1995, 2000, 2002 
' This version dated 2002 Nov. 01 
  
DECLARE SUB InNum (Pt$, V, MX%, CFlag%) 
DECLARE FUNCTION YesNo$ () 
DECLARE FUNCTION ROff$ (J, P%) 
DECLARE SUB Instructions () 
DECLARE SUB Spacebar () 
  
CONST R = 6.615 
' Radius of geosynchronous orbit in units of earth radius 
CONST PI = 3.141593 
CONST CF = PI / 180  ' Degree/radian factor 
  
ON ERROR GOTO Etrap 
F% = FREEFILE 
N$ = "TVSATDAT.DAT" 
E% = 0 
OPEN N$ FOR INPUT AS F% 
  
IF E% THEN 
  CFlag% = 1 
  
  LA = 45   ') 
  LO = 80   ') initial default latitude and longitudes 
  LS = 90   ') 
  
ELSE 
  INPUT #F%, LA, LO, LS 
  CFlag% = 0 
  
END IF 
  
CLOSE F% 
  
ON ERROR GOTO 0 
  
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"; ROff$(LA, -2); "degrees north." 
PRINT "Longitude is"; ROff$(LO, -2); "degrees west." 
PRINT 
PRINT "Keep these values"; 
  
IF YesNo$ = "n" THEN 
  
    PRINT "Input latitude and longitude in degrees." 
    PRINT "Use negative numbers for angles in" 
    PRINT "opposite directions to those shown." 
    PRINT 
  
    InNum "Your latitude (deg. north)", LA, 90, CFlag% 
  
    InNum "Your longitude (deg. west)", LO, 180, CFlag% 
  
END IF 
  
LT = LA 
LG = LO 
  
PRINT "For current satellite:" 
PRINT "Longitude is"; ROff$(LS, -2); "degrees west." 
PRINT 
PRINT "Keep this value"; 
  
IF YesNo$ = "n" THEN 
  
  PRINT "Input longitude of satellite in" 
  PRINT "degrees west. (Use negative number" 
  PRINT "if you want to input degrees east.)" 
  PRINT 
  
  InNum "Satellite's longitude", LS, 180, CFlag% 
  
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 = SGN(Z) * PI / 2 
' 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 
  
IF Y > -1 THEN 
  
  PRINT "Satellite is invisible (below horizon)." 
  
ELSE 
  
  ' calculate altitude 
  
  IF X = 0 AND Z = 0 THEN 
    AL = PI / 2 
  ELSE 
    AL = ATN((-1 - Y) / SQR(X * X + Z * Z)) 
  END IF 
  
  AL = AL / CF ' alt. in degrees 
  A$ = ROff$(AL, 1) 
  
  IF VAL(A$) = 90 THEN 
  
     PRINT "Satellite is vertically overhead." 
     PRINT "(Altitude is 90.0 degrees.)" 
  
  ELSE 
  
    ' calculate new bearing 
  
    IF Z = 0 THEN 
      BE = SGN(X) * PI / 2 
    ELSE 
      BE = ATN(X / Z) 
      IF Z < 0 THEN BE = BE + PI 
    END IF 
  
    BE = BE / CF ' bearing in degrees 
  
    ' roundoff and printout 
  
    BE = VAL(ROff$(BE, 1)) 
    BE = BE - 360 * INT(BE / 360) 
  
    PRINT "Satellite's position:" 
    PRINT 
  
    B$ = ROff$(BE, 1) 
    BE = VAL(B$) 
  
    PRINT "Bearing is"; B$; "degrees." 
  
    X = BE / 90 
    IF X = INT(X) THEN 
      Y1$ = RTRIM$(MID$("NorthEast SouthWest", 5 * X + 1, 5)) 
      PRINT "(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 
      PRINT Y1$; ROff$(Z, 1); "degrees "; Y2$ 
    END IF 
  
    PRINT 
    PRINT "Altitude is"; A$; "degrees." 
    IF AL < 5 THEN 
      PRINT 
      PRINT "Very low altitude - Communication unreliable" 
    END IF 
  
  END IF 
  
END IF 
  
PRINT 
  
PRINT "Another calculation"; 
  
IF YesNo$ = "y" GOTO Newcalc 
  
IF CFlag% THEN 
  PRINT "Keep current latitude & longitudes 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%, LA; ","; LO; ","; LS 
    CLOSE F% 
    ON ERROR GOTO 0 
  END IF 
END IF 
  
END 
  
Etrap: 
  E% = 1 
  RESUME NEXT 
  
  
SUB InNum (Pt$, V, MX%, CFlag%) 
  
 DO 
  
   PRINT Pt$; " (or ENTER for no change)? "; 
   LINE INPUT IN$ 
   PRINT 
   IF IN$ = "" THEN 
     PRINT "Value unchanged" 
     PRINT 
     EXIT DO 
   END IF 
   W = VAL(IN$) 
   IF ABS(W) <= MX% AND INSTR(IN$, ",") = 0 THEN 
     IF W <> V THEN 
       V = W 
       CFlag% = 1 
     END IF 
     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 true" 
   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 "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 "Note that the bearings are true. If you want a magnetic" 
   PRINT "bearing, look up the local magnetic deviation. To find" 
   PRINT "the magnetic bearing, add the deviation to the true" 
   PRINT "bearing if magnetic north is west of true north. Subtract" 
   PRINT "the deviation if magnetic north is east of true north." 
   PRINT 
   PRINT "Your entries of latitude and longitudes can be kept on disk" 
   PRINT "and used in subsequent runs. Initially, arbitary defaults" 
   PRINT "are used." 
  
   CALL Spacebar 
  
END SUB 
  
FUNCTION ROff$ (J, P%) 
  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 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