TIP: Click on subject to list as thread! ANSI
echo: quik_bas
to: ALL
from: DAVID WILLIAMS
date: 2003-04-02 22:26:38
subject: satellites again

I managed to check that program I posted here recently that calculates 
the dates and times when geostatinary satellites pass in front of the 
sun, causing signal outages. The program is very accurate. 
  
For anyone who might be interested, here's a version that combines the 
outage calculation with a calculation of the position of the satellite 
in the sky, expressed as compass bearing (azimuth) and angle of 
elevation. 
  
Have fun! 
  
                                       dow 
  
---------------------------------------------------- 
  
' TV_SATS.BAS 
' TV Satellites, Commodore PET version, David Williams, 1982 
  
' david.williams{at}ablelink.org 
  
' Updated for other computers 1995, 2000, 2002 
' Solar outage code added March 2003 
  
' This version dated 2003 March 29 
  
DECLARE SUB InNum (Pt$, V, MX%) 
DECLARE FUNCTION YesNo$ () 
DECLARE FUNCTION ROff$ (J, P%) 
DECLARE FUNCTION ATan (A, B) 
DECLARE SUB Instructions () 
DECLARE SUB SpaceBar () 
  
DIM SHARED PI 
PI = 4 * ATN(1) 
CF = PI / 180  ' Degree/radian factor 
CF2 = PI / 182.62 
R = 6.615  ' Radius of geosynchronous orbit in units of earth 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 
  
ELSE 
  
  INPUT #F%, FA, FG, FS, FZ 
  
END IF 
  
CLOSE F% 
  
ON ERROR GOTO 0 
  
LA = FA 
LG = FG 
LS = FS 
TZ = FZ 
  
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$(LG, -2); "degrees west." 
  
PRINT 
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 
InNum "Your longitude (deg. west)", LG, 180 
  
PRINT "For current satellite:" 
PRINT "Longitude is"; ROff$(LS, -2); "degrees west." 
PRINT 
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 
  
PRINT "Current time zone offset is"; 
PRINT ROff$(TZ, -1); "hours from GMT / UT." 
InNum "Time zone (+/- hours from G.M.T. / U.T.)", TZ, 14 
  
PRINT "Confirm these entries"; 
IF YesNo$ = "n" GOTO Newcalc 
  
LD = (LG - LS) * CF ' longitude difference in radians 
LR = CF * LA        ' latitude in radians for outages 
LT = LR + PI / 2    ' latitude from s. pole 
  
' 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 
  
CLS 
  
IF Y > -1 THEN 
  
  PRINT "Satellite is invisible (below horizon)." 
  
ELSE 
  
  ' calculate altitude 
  
  AL = ATan(-1 - Y, SQR(X * X + Z * Z)) 
  
  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 
  
    BE = ATan(X, Z) 
  
    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 
  
  ' Solar outage calculation 
  
  ' coordinate differences 
  X1 = -R * SIN(LD) 
  Y1 = -SIN(LR) 
  Z1 = R * COS(LD) - COS(LR) 
  
  ' calculate apparent declination 
  AD = ATN(Y1 / SQR(X1 * X1 + Z1 * Z1)) 
  
  ' calculate apparent longitude 
  AG = ATan(X1, Z1) 
  
  V = 4 * (AG / CF + LG) + 60 * TZ + 720.5 
  
  PRINT 
  PRINT "Nominal solar outage dates and times:" 
  
  ' find dates (days after New Year) when sun is at same declination 
  FOR J = 0 TO 1 
  
    U = 171 
    L = 342 * J 
    FOR K = 1 TO 20 
      D = (U + L) / 2 
      A = CF2 * (D + 11) + .033 * SIN(CF2 * (D - 2)) 
      IF INT(L) = INT(U) THEN EXIT FOR 
      S = -.3987 * COS(A) 
      S = ATN(S / SQR(1 - S * S)) 
      IF S > AD THEN U = D ELSE L = D 
    NEXT 
    D = INT(D) 
  
    ' convert date to month and day 
    M = 0 
    RESTORE 
    DO 
      READ L 
      IF D < L THEN EXIT DO 
      D = D - L 
      M = M + 1 
    LOOP 
    DATA 31,28,31,30,31,30,31,31,30,31,30,31 
    D = D + 1 
    M$ = MID$("JanFebMarAprMayJunJulAugSepOctNovDec", 3 * M + 1, 3) 
    PRINT M$; D, 
  
    ' calculate equation of time 
    ET = 9.7 * SIN(A + A) + 7.9 * SIN(A - .224) 
  
    ' calculate outage time 
    T = INT(V + ET) 
    T = T - 1440 * INT(T / 1440) 
    H = T \ 60 
    M = T MOD 60 
    PRINT H; "hrs,"; M; "mins" 
  
  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 IF 
  
PRINT 
  
PRINT "Another calculation"; 
  
IF YesNo$ = "y" GOTO Newcalc 
  
IF LA  FA OR LG  FG OR LS  FS OR TZ  FZ THEN 
  PRINT "Keep current latitude, longitudes & time zone 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) 
    END IF 
    CLOSE F% 
    ON ERROR GOTO 0 
  END IF 
END IF 
  
END 
  
Etrap: 
  E% = 1 
  RESUME NEXT 
  
  
FUNCTION ATan (A, B) 
  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 InNum (Pt$, V, MX%) 
  
 DO 
  
   PRINT Pt$; " (or ENTER for no change)? "; 
   LINE INPUT in$ 
   IF in$ = "" THEN 
     PRINT "Value unchanged ("; LTRIM$(RTRIM$(ROff$(V, -3)));
")" 
     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 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 "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." 
   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. You will also be asked for your" 
   PRINT "time zone, which is used in the solar outage calculation." 
  
  
   CALL SpaceBar 
  
   CLS 
   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 "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." 
  
   CALL SpaceBar 
  
   CLS 
   PRINT "Your entries of latitude, longitudes and time zone can be" 
   PRINT "kept on disk and used in subsequent runs. Initially," 
   PRINT "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 LEN(INKEY$) 
  LOOP 
  
  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 
  
------------------------------------------------------ 
--- Platinum Xpress/Win/WINServer v3.0pr5
* Origin: The Bayman BBS,Toronto, (416)698-6573 - 1:250/514 (1:250/514)
SEEN-BY: 633/267 270
@PATH: 250/514 140/1 106/2000 633/267

SOURCE: echomail via fidonet.ozzmosis.com

Email questions or comments to sysop@ipingthereforeiam.com
All parts of this website painstakingly hand-crafted in the U.S.A.!
IPTIA BBS/MUD/Terminal/Game Server List, © 2025 IPTIA Consulting™.