'===========================================================================
' Subject: Find TV sattalite position         Date: 09-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 August 29 
  
DECLARE SUB InNum (Pt$, V, MX%, CFlag%) 
DECLARE FUNCTION YesNo$ () 
DECLARE FUNCTION ROff$ (J, P%) 
DECLARE SUB Instructions () 
DECLARE SUB Spacebar () 
  
CONST R = 6.61614 
' 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 
  
' calculate altitude 
  
IF X = 0 AND Z = 0 THEN 
  AL = SGN(-1 - Y) * PI / 2 
ELSE 
  AL = ATN((-1 - Y) / SQR(X * X + Z * Z)) 
END IF 
  
IF AL < 0 THEN 
  
  PRINT "Satellite invisible (below horizon)" 
  
ELSE 
  
  ' calculate new bearings 
  
  IF Z = 0 THEN 
    BE = SGN(X) * PI / 2 
  ELSE 
    BE = ATN(X / Z) - PI * (Z < 0) ' bearing 
  END IF 
  
  BE = BE / CF ' bearing in degrees 
  AL = AL / CF ' alt. 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"; ROff$(AL, 1); "degrees." 
  IF AL < 5 THEN 
    PRINT 
    PRINT "Very low altitude -" 
    PRINT "Communication unreliable" 
  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" 
   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 "Note that the bearings are true. If" 
   PRINT "you want a magnetic bearing, you must" 
   PRINT "find the local magnetic deviation from" 
   PRINT "a current map. The magnetic bearing is" 
   PRINT "found by adding the deviation to the true" 
   PRINT "bearing if magnetic north is west of true" 
   PRINT "north, or subtracting the deviation if" 
   PRINT "magnetic north is east of true north." 
   PRINT 
   PRINT "Your entries of latitude and longitudes" 
   PRINT "are kept on disk and used subsequently." 
   PRINT "Initially, arbitary defaults 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