Program: Vsun.bas

DEFDBL A-Z
'                       VSUN.BAS
'       Revision 1.0    Original                         10 Dec 1996
'
DECLARE SUB CALCULATE (t$)
DECLARE FUNCTION FUN# (a#, t#)

DIM SHARED ttime(400), dt(4)
DIM SHARED v, accuracy, time
DIM SHARED aearth, ajupiter, acallisto

CONST day = 24# * 60 * 60
CONST year = 365.25# * 24 * 60 * 60
CONST pi = 3.141592653589793#
CONST pi2 = 2 * pi / 360
CONST wearth = 2 * pi / year       ' omega earth
CONST rearth = 149600000
CONST tjupiter = 11.86#               ' 11.86
CONST wjupiter = 2 * pi / (tjupiter * year)   ' omega jupiter
CONST rjupiter = 778300000
CONST tcallisto = 16.69#              '16.69
CONST wcallisto = 2 * pi / (tcallisto * day)   ' omega callisto
CONST rcallisto = 1888000
CONST c = 300000             ' km/sec
vbase = 300                  ' km/sec
CONST xcenter = 640 / 2: ycenter = 480 / 2

CLS : SCREEN 12
accuracy = 25

LOCATE 1, 20: PRINT "Initialization Display"
INPUT "v in km/sec "; v1
INPUT "phi in degrees "; phi

IF v1 = 0 THEN vbase = 300 ELSE vbase = v1
  v = vbase
  PRINT "v = "; v; "km/sec   phi = "; phi; " degrees"
  acallisto = 0: ajupiter = 0
  aearth = 0
  CALCULATE "t1"
  t1 = time
  PRINT
  aearth = pi
  CALCULATE "t2"
  t2 = time
  PRINT USING "t2-t1 ######.##"; TAB(50); t2 - t1
  acallisto = pi: ajupiter = pi
  aearth = 0
  CALCULATE "t3"
  PRINT
  t3 = time
  aearth = pi
  CALCULATE "t4"
  t4 = time
  PRINT USING "t3-t4 ######.##"; TAB(50); t3 - t4
  PRINT USING "delta ######.##"; TAB(50); t3 - t4 - (t2 - t1)
  INPUT "next "; a$

CLS
COLOR 15
LOCATE 1, 20: PRINT "Technical Explanation Display"
COLOR 12
LOCATE 2, 1: PRINT " Inner circle = Earth"
COLOR 7
LOCATE 3, 1: PRINT " Outer circle ="
LOCATE 4, 3: PRINT " Jupiter & Callisto"
COLOR 14: LOCATE 16, 40: PRINT "S"
PSET (xcenter, ycenter)
LOCATE 5, 1: PRINT " Center = Sun"

r1 = 200# * rearth * 1.5 / rjupiter: r2 = 200# * rjupiter / rjupiter
CIRCLE (xcenter, ycenter), r1, 12
CIRCLE (xcenter, ycenter), r2, 7

a1 = 0 * pi / 180                               ' start, earth and jupiter

x1 = r1 * COS(a1): y1 = r1 * SIN(a1): x2 = r2 * COS(a1): y2 = r2 * SIN(a1)
LINE (xcenter + x1, ycenter - y1)-(xcenter + x2, ycenter - y2), 15
a31 = a1 + 15 * pi / 180
x31 = r1 * COS(a31): y31 = r1 * SIN(a31):
LINE (xcenter + x31, ycenter - y31)-(xcenter + x2, ycenter - y2), 15, , &HCCCC
COLOR 15: LOCATE 15, 45: PRINT "t2"

a2 = a1 + 20 * pi / 180                         ' jupiter
a3 = a1 + 30 * pi / 180                         ' earth

x3 = r1 * COS(a3): y3 = r1 * SIN(a3): x4 = r2 * COS(a2): y4 = r2 * SIN(a2)
LINE (xcenter + x3, ycenter - y3)-(xcenter + x4, ycenter - y4), 7
x5 = x4 + 2 * r1 * COS(a1): y5 = y4 + 2 * r1 * SIN(a1)
LINE (xcenter + x4, ycenter - y4)-(xcenter + x5, ycenter - y5), 15   't0
x6 = x4 + 2 * r1 * COS(a2): y6 = y4 + 2 * r1 * SIN(a2)
LINE (xcenter + x4, ycenter - y4)-(xcenter + x6, ycenter - y6), 7   't1
a31 = a3 + 15 * pi / 180
x31 = r1 * COS(a31): y31 = r1 * SIN(a31):
LINE (xcenter + x31, ycenter - y31)-(xcenter + x4, ycenter - y4), 7, , &HCCCC
COLOR 15: LOCATE 12, 72: PRINT "t0"
COLOR 7: LOCATE 9, 72: PRINT "t1"
COLOR 7: LOCATE 14, 44: PRINT "t2"
COLOR 7
LOCATE 16, 67: PRINT "1"
LOCATE 12, 66: PRINT "2"
LOCATE 8, 63: PRINT "3"
LOCATE 4, 55: PRINT "4"

a1 = a2                                         ' jupiter old
a2 = a2 + 20 * pi / 180                         ' jupiter
a3 = a3 + 30 * pi / 180                         ' earth
x3 = r1 * COS(a3): y3 = r1 * SIN(a3): x4 = r2 * COS(a2): y4 = r2 * SIN(a2)
LINE (xcenter + x3, ycenter - y3)-(xcenter + x4, ycenter - y4), 14
x5 = x4 + 2 * r1 * COS(a1): y5 = y4 + 2 * r1 * SIN(a1)
LINE (xcenter + x4, ycenter - y4)-(xcenter + x5, ycenter - y5), 7   't0
x6 = x4 + 2 * r1 * COS(a2): y6 = y4 + 2 * r1 * SIN(a2)
LINE (xcenter + x4, ycenter - y4)-(xcenter + x6, ycenter - y6), 14  't1
a31 = a3 + 15 * pi / 180
x31 = r1 * COS(a31): y31 = r1 * SIN(a31):
LINE (xcenter + x31, ycenter - y31)-(xcenter + x4, ycenter - y4), 14, , &HCCCC
COLOR 7: LOCATE 7, 68: PRINT "t0"
COLOR 14: LOCATE 4, 65: PRINT "t1"
COLOR 14: LOCATE 13, 42: PRINT "t2"

a1 = a2                                         ' jupiter old
a2 = a2 + 20 * pi / 180                         ' jupiter
a3 = a3 + 30 * pi / 180                         ' earth
x3 = r1 * COS(a3): y3 = r1 * SIN(a3): x4 = r2 * COS(a2): y4 = r2 * SIN(a2)
LINE (xcenter + x3, ycenter - y3)-(xcenter + x4, ycenter - y4), 13
a31 = a3 + 15 * pi / 180
x31 = r1 * COS(a31): y31 = r1 * SIN(a31):
LINE (xcenter + x31, ycenter - y31)-(xcenter + x4, ycenter - y4), 13, , &HCCCC

  COLOR 15
  LOCATE 30, 1
  INPUT "next "; a$

CLS

pass = 1

s00:
  t = 0                                   ' time allignment
  tmin = 0: tmax = 0: min1 = 0: min2 = 0
  totalx = 0: totalmin = 0: totalmax = 0
  printline = 1
  IF pass = 0 THEN v = 0:  ELSE v = vbase
  IF pass < 2 THEN
    PRINT "v = "; v; "km/sec   phi = "; phi; " degrees"
    COLOR 14
    PRINT TAB(4); "t0"; TAB(16); "t1"; TAB(24); "t1 - t0"; TAB(33); "time";
    PRINT TAB(41); "dist"; TAB(52); "t2=t1+time"; TAB(68); "delta t2"
    COLOR 15
  ELSE
    INPUT "next "; a$
    CLS
    PSET (xcenter, ycenter)
    COLOR 15: PRINT "Jupiter "
    COLOR 14: PRINT "Callisto "
    COLOR 13: PRINT "Earth "
    COLOR 14: LOCATE 16, 40: PRINT "S"
    PSET (xcenter, ycenter)
  END IF
FOR n = 1 TO 320
  ' calculate moment t that sun, jupiter and callisto are in one line
  t0 = t
  FOR i = 1 TO accuracy
    ajupiter = FUN(wjupiter, t)
    acallisto = FUN(wcallisto, t)
    IF acallisto > 6 AND ajupiter < 1 THEN acallisto = acallisto - 2 * pi
    deltat = (ajupiter - acallisto) / wcallisto
    t = t + deltat
  NEXT i
  t1 = t
  ' position Callisto and Jupiter
  xc = rcallisto * COS(acallisto): yc = rcallisto * SIN(acallisto):
  xj = rjupiter * COS(ajupiter): yj = rjupiter * SIN(ajupiter):
  xc = xc + xj: yc = yc + yj:

  ' Calculate time time when signal reaches earth.
  time = 2000#
  FOR i = 1 TO 2 * accuracy
    aearth = FUN(wearth, t + time)
    ' position earth
    xe = rearth * COS(aearth) + v * time * COS(phi * pi / 180): ye = rearth * SIN(aearth) + v * time * SIN(phi * pi / 180)
    dist2 = (xc - xe) * (xc - xe) + (yc - ye) * (yc - ye)
    dist = SQR(dist2)
    time = dist / c                                 ' time
  NEXT i
  t2 = t + time
 
  IF pass = 2 THEN
      year1 = INT(t2 / year)
      IF year1 MOD 2 = 0 THEN COLOR 15:  ELSE COLOR 7
      PSET (xcenter * (1 + .9 * xj / (rjupiter * 1.36)), ycenter * (1 - .9 * yj / rjupiter))
      PSET (xcenter * (1 + .97 * xc / (rjupiter * 1.36)), ycenter * (1 - .97 * yc / rjupiter)), 14
      PSET (xcenter * (1 + .9 * xe / (rjupiter * 1.36)), ycenter * (1 - .9 * ye / rjupiter)), 13
  END IF
 
  ttime(n) = t + time:
  dt(2) = ttime(n) - ttime(n - 1)
  IF printline = 1 AND pass <= 1 THEN
    PRINT INT(t0); TAB(12); INT(t1);
    PRINT USING "######.#"; TAB(23); t1 - t0;
    PRINT TAB(32); INT(time); TAB(40); INT(dist); TAB(54); INT(t + time);
    PRINT TAB(67); INT(dt(2))
  END IF
  IF n >= 4 THEN
    dt(2) = ttime(n) - ttime(n - 1)       'now
    dt(1) = ttime(n - 1) - ttime(n - 2)
    dt(0) = ttime(n - 2) - ttime(n - 3)
    IF dt(2) < dt(1) AND dt(1) > dt(0) THEN
      a1 = dt(0)
      b1 = (-3 * dt(0) + 4 * dt(1) - dt(2)) / 2
      c1 = (dt(0) - 2 * dt(1) + dt(2)) / 2
      t1 = -b1 / (2 * c1)
      max1 = a1 + b1 * t1 + c1 * t1 * t1
    END IF
    IF dt(2) > dt(1) AND dt(1) < dt(0) THEN
      a1 = dt(0)
      b1 = (-3 * dt(0) + 4 * dt(1) - dt(2)) / 2
      c1 = (dt(0) - 2 * dt(1) + dt(2)) / 2
      t1 = -b1 / (2 * c1)
      min1 = a1 + b1 * t1 + c1 * t1 * t1
      IF pass < 2 THEN
        IF printline = 1 THEN
          COLOR 14
          PRINT TAB(38); "max"; TAB(49); "min"; TAB(58); "max-min"; TAB(69); "int dt"
          COLOR 15
        END IF
        PRINT "min"; INT(dt(2)); INT(dt(1)); INT(dt(0));
        PRINT USING "#######.#"; TAB(35); max1;
        PRINT USING "#######.#"; TAB(46); min1;
        PRINT USING "#####.#"; TAB(57); max1 - min1;
        PRINT USING "#####.#"; TAB(68); totalx / 2
        IF printline = 1 THEN
          lpos1 = CSRLIN
          INPUT "next "; a$
          LOCATE lpos1, 1
        END IF
        printline = 0
      END IF
      IF totalmax = 0 THEN totalmax = totalx: totalmin = totalx
      IF totalx > totalmax THEN totalmax = totalx
      IF totalx < totalmin THEN totalmin = totalx
      totalx = 0
    END IF
    IF min1 > 0 THEN
      totalx = totalx + dt(2) - min1
    END IF
  END IF
  t = t + tcallisto * day         ' update time for next allignment
NEXT n
  pass = pass + 1
  IF pass < 2 THEN
    PRINT "Time start "; INT(ttime(1)); "Time end"; INT(ttime(n - 1)); "Average time "; INT((ttime(n - 1) - ttime(1)) / (n - 2))
  END IF
  IF pass <= 2 THEN
    dtt = (totalmax - totalmin) / 2
    PRINT USING "min #####.#"; totalmin / 2;
    PRINT USING " max #####.#"; totalmax / 2;
    PRINT USING "   max-min ####.#"; dtt
    PRINT USING "v1 ######.#"; dtt * c * c / (4 * rearth);
    PRINT USING " v2 ######.#"; dtt * c * c / (2 * pi * rearth);
    PRINT USING "   (v) ######.#"; v
  END IF
IF pass < 3 THEN GOTO s00

END

SUB CALCULATE (t$)
'                                                                CALCULATE
  ' position Callisto and Jupiter
  xc = rcallisto * COS(acallisto): yc = rcallisto * SIN(acallisto):
  xj = rjupiter * COS(ajupiter): yj = rjupiter * SIN(ajupiter):
  xc = xc + xj: yc = yc + yj:

  ' Calculate time time when signal reaches earth.
  time = 2000#
  FOR i = 1 TO 2 * accuracy
    ' position earth
    xe = rearth * COS(aearth) + v * time: ye = rearth * SIN(aearth)
    dist2 = (xc - xe) * (xc - xe) + (yc - ye) * (yc - ye)
    dist = SQR(dist2)
    time = dist / c                                 ' time
  NEXT i
  PRINT "angle"; ajupiter * 180 / pi; TAB(12); "dist "; INT(dist); TAB(30); t$; " time ";
  PRINT USING " ######.##"; time;

END SUB

FUNCTION FUN (a, t)
fun1 = a * t
f1:
IF fun1 >= 2 * pi THEN fun1 = fun1 - 2 * pi: GOTO f1
FUN = fun1
END FUNCTION