Program: Hubble.bas
Implementation Select: Implementation
Operation Select: Program: HUBBLE.BAS
DEFSNG A-H, M-Z
DEFINT I-K
' HUBBLE.BAS
' Revision 1.0 Original 11 Feb 1996
' Revision 2 Modified input parameters 21 Feb 2016
'
DECLARE SUB COORD (m0!, m1!, dist!, x0!, x1!, vy0!, vy1!)
DECLARE SUB GETSCREEN ()
DECLARE SUB INITIALISE ()
CONST ESC = 27, ENTER = 13
CONST pi = 3.141592653589#
DIM SHARED x0(10), x1(10), y0(10), y1(10), x2(10), y2(10)
DIM SHARED vx0(10), vx1(10), vy0(10), vy1(10)
DIM SHARED x00(10), y00(10) ' initial position
DIM SHARED vx00(10), vy00(10) ' initial velocity
DIM SHARED xx0(10, 2), yy0(10, 2), rx(10, 10)
DIM SHARED ax(10), ay(10)
DIM SHARED Sax(10), Say(10) ' sum of ax integrated
DIM SHARED SSax(10), SSay(10) ' sum of sum ax integrated
DIM SHARED Sva(10), Svy(10) ' sum of vx integrated
DIM SHARED rr(10, 10), dx(10, 10), dy(10, 10)
DIM SHARED m(10)
DIM SHARED dist1 ' distance 1
DIM SHARED dist2 ' distance 2
DIM SHARED t, tcnt, TT ' t + tdelta
DIM SHARED G, initcond, state
DIM SHARED stars
DIM SHARED tdelta
DIM SHARED scren%
DIM SHARED xmax, ymax, xchrmax, ychrmax, fract
DIM SHARED title$
DIM SHARED x0disp, y0disp, scalex, scaley
DIM SHARED dispcond ' - 1 once each rev
DIM SHARED endc ' max number of endc
DIM SHARED agrad, alpha0, alpha, alphastr, alphaend, dalpha
title$ = "3 Objects m0, m1 and m2 Demonstration"
s0:
CLS
COLOR 15: PRINT "Alpha start angle": LOCATE 1, 32: INPUT ""; a$
alphastr = VAL(a$)
COLOR 15: PRINT "Alpha end angle ": LOCATE 2, 32: INPUT ""; a$
alphaend = VAL(a$)
COLOR 15: PRINT "delta t between 0.001 and 0.01": LOCATE 3, 32: INPUT ""; a$
deltat = VAL(a$)
IF deltat < .001 THEN deltat = .001
IF deltat > .01 THEN deltat = .01
tdelta = deltat
alpha0 = alphastr
dalpha = (alphaend - alphastr) / 10
IF dalpha <= 0 THEN dalpha = 0: GOTO s01
IF dalpha <= .1 THEN dalpha = .1: GOTO s01
IF dalpha <= 1 THEN dalpha = 1: GOTO s01
IF dalpha <= 18 THEN dalpha = 10: GOTO s01
dalpha = 20
s01:
PRINT "alphastr "; alphastr; "alphaend"; alphaend; "dalpha"; dalpha; deltat; t; deltat
COLOR 15: INPUT "Input parameters okay (Y or N) ", a$
IF UCASE$(a$) = "N" THEN GOTO s0
s1:
INITIALISE ' Part 2
IF color15 <> 15 THEN color15 = 15 ELSE color15 = 7
IF COLOR12 <> 12 THEN COLOR12 = 12 ELSE COLOR12 = 12
IF COLOR10 <> 10 THEN COLOR10 = 10 ELSE COLOR10 = 10
COLOR 15
IF alpha0 = alphastr THEN CLS : SCREEN scren%
LOCATE 1, 25: PRINT title$
LOCATE 2, 1: PRINT "Delta time "; tdelta; dispcond; "# stars "; stars;
COLOR color15: LOCATE 3, 1: PRINT "alpha "; alpha;
stars = stars - 1
t = 0: tcnt = 0
end1 = 0
time1 = TIMER
FOR i = 0 TO stars
FOR j = 0 TO stars
IF i <> j AND m(j) > 0 THEN
dx = (x00(j) - x00(i))
dy = (y00(j) - y00(i))
rr = dx * dx + dy * dy:
rr(i, j) = rr: dx(i, j) = dx: dy(i, j) = dy
END IF
NEXT j
x1(i) = x0(i): y1(i) = y0(i) ' all equal
x2(i) = x1(i): y2(i) = y1(i)
NEXT i
FOR i = 0 TO stars
Sax(i) = 0: SSax(i) = 0
Say(i) = 0: SSay(i) = 0
NEXT i
' calculate for t = -1
tcnt = -1
FOR i = 0 TO stars
ax(i) = 0: ay(i) = 0
IF m(i) >= 0 THEN
FOR j = 0 TO stars
IF i <> j THEN
rr = rr(i, j): dx = dx(i, j): dy = dy(i, j)
a = m(j) * G / rr
ax(i) = ax(i) + dx * a / SQR(rr)
ay(i) = ay(i) + dy * a / SQR(rr)
END IF
NEXT j
Sax(i) = Sax(i) + ax(i)
Say(i) = Say(i) + ay(i)
SSax(i) = SSax(i) + Sax(i)
SSay(i) = SSay(i) + Say(i)
END IF
NEXT i
FOR i = 0 TO stars
x1(i) = x00(i) + tdelta * (tcnt * vx00(i) + tdelta * SSax(i))
y1(i) = y00(i) + tdelta * (tcnt * vy00(i) + tdelta * SSay(i))
FOR j = 0 TO stars
IF i <> j THEN
dx = (x00(j) - x00(i)) + tdelta * (tcnt * (vx00(j) - vx00(i)) + tdelta * (SSax(j) - SSax(i)))
dy = (y00(j) - y00(i)) + tdelta * (tcnt * (vy00(j) - vy00(i)) + tdelta * (SSay(j) - SSay(i)))
rr = dx * dx + dy * dy
rr(i, j) = rr: dx(i, j) = dx: dy(i, j) = dy
IF rr < 1 THEN end1 = 1
END IF
NEXT j
NEXT i
' initialise rr(i,j)
FOR i = 0 TO stars
FOR j = 0 TO stars
IF i <> j THEN
dx = (x00(j) - x00(i))
dy = (y00(j) - y00(i))
rr = dx * dx + dy * dy
rr(i, j) = rr: dx(i, j) = dx: dy(i, j) = dy
END IF
NEXT j
NEXT i
FOR i = 0 TO stars
Sax(i) = 0: SSax(i) = 0
Say(i) = 0: SSay(i) = 0
NEXT i
t = 0: tcnt = 0: endc = 0
str1:
tcnt = tcnt + 1: t = tcnt * tdelta
IF end1 = 0 THEN
FOR i = 0 TO stars
ax(i) = 0: ay(i) = 0
IF m(i) >= 0 THEN
FOR j = 0 TO stars
IF i <> j THEN
rr = rr(i, j)
dx = dx(i, j): dy = dy(i, j)
a = m(j) * G / rr
ax(i) = ax(i) + dx * a / SQR(rr)
ay(i) = ay(i) + dy * a / SQR(rr)
END IF
NEXT j
' vnt = v0 + (a0 + a1 ...an)t = v0 + Sax * t
' xn = x0 + nv0 + (n*a0 + n-1*a1 ... 1*an) tư = x0 + n*v0 + SSax*tư
Sax(i) = Sax(i) + ax(i)
Say(i) = Say(i) + ay(i)
SSax(i) = SSax(i) + Sax(i)
SSay(i) = SSay(i) + Say(i)
' vx1(i) = vx0(i) + ax(i) * tdelta
' vy1(i) = vy0(i) + ay(i) * tdelta
' x0(i) = x00(i) + tdelta * (tcnt * vx00(i) + tdelta * SSax(i))
' y0(i) = y00(i) + tdelta * (tcnt * vy00(i) + tdelta * SSay(i))
END IF
NEXT i
'x0(1) = dist1 / 2 * COS(2 * pi * t / TT + agrad)
'y0(1) = dist1 / 2 * SIN(2 * pi * t / TT + agrad)
'x0(0) = -x0(1)
'y0(0) = -y0(1)
FOR i = 0 TO stars
x2(i) = x1(i): y2(i) = y1(i)
x1(i) = x0(i): y1(i) = y0(i)
x0(i) = x00(i) + tdelta * (tcnt * vx00(i) + tdelta * SSax(i))
y0(i) = y00(i) + tdelta * (tcnt * vy00(i) + tdelta * SSay(i))
FOR j = 0 TO stars
IF i <> j THEN
' dx = (x00(j) - x00(i)) + tdelta * (tcnt * (vx00(j) - vx00(i)) + tdelta * (Ssax(j) - Ssax(i)))
' dy = (y00(j) - y00(i)) + tdelta * (tcnt * (vy00(j) - vy00(i)) + tdelta * (Ssay(j) - Ssay(i)))
dx = x0(j) - x0(i) ' low accuracy
dy = y0(j) - y0(i) ' low accuracy
rr = dx * dx + dy * dy
rr(i, j) = rr
dx(i, j) = dx: dy(i, j) = dy
END IF
NEXT j
NEXT i
rrmin = rr(2, 0)
IF rr(2, 1) < rrmin THEN rrmin = rr(2, 1)
IF state = 0 THEN
IF rrmin < .5 * dist2 * dist2 THEN state = 1
ELSE
IF rrmin > 2.5 * dist2 * dist2 THEN end1 = 1
IF rrmin < 10 THEN end1 = 1
END IF
END IF
IF (tcnt MOD dispcond = 0 AND dispcond > 0) OR end1 = 1 OR (dispcond < 0 AND prtest% = 1) THEN
vx1(2) = vx0(2) + Sax(2) * tdelta
vy1(2) = vy0(2) + Say(2) * tdelta
vnor = SQR(vx1(2) * vx1(2) + vy1(2) * vy1(2))
vesc = (SQR(2 * (m(0) / SQR(rr(2, 0)) + m(1) / SQR(rr(2, 1)))))
COLOR color15: LOCATE 2, 40: PRINT "cycle "; tcnt; TAB(56); "r #3"; INT(SQR((rr(2, 0) + rr(2, 1)) / 2))
LOCATE 3, 20: PRINT USING "v escape ##.###"; vesc
IF vnor > vesc THEN COLOR 12
LOCATE 3, 40: PRINT "v #3 "; : PRINT USING "##.###"; vnor
PSET (x0disp + x0(0) * scalex * fract, y0disp - y0(0) * scaley), COLOR12
PSET (x0disp + x0(1) * scalex * fract, y0disp - y0(1) * scaley), COLOR10
COLOR color15
IF vnor > vesc THEN COLOR 12
IF stars > 1 THEN PSET (x0disp + x0(2) * scalex * fract, y0disp - y0(2) * scaley)
END IF
s9:
a$ = INKEY$: IF a$ = CHR$(ESC) THEN GOTO s10:
IF end1 = 1 THEN GOTO s10:
GOTO str1
s10:
IF subtest% = 2 AND alpha0 < alphaend THEN
alpha0 = alpha0 + dalpha:
IF end1 = 1 THEN GOTO s1
END IF
COLOR 14
LOCATE ychrmax, 35: PRINT "NEXT TEST ? ";
DO: a$ = INKEY$: LOOP UNTIL a$ <> ""
IF a$ = CHR$(ESC) THEN END
COLOR 15
GOTO s0
SUB COORD (m0, m1, dist, x0, x1, vy0, vy1)
' COORD
' x0 and vy0 are base position and speed
mm01 = m0 + m1
r0 = dist * m1 / mm01
r1 = dist - r0
vr0 = SQR(G * r0 * m1) / dist
vr1 = SQR(G * r1 * m0) / dist
x1 = x0 ' same base
vy1 = vy0 ' same base
x0 = x0 - r0
x1 = x1 + r1
vy0 = vy0 - vr0
vy1 = vy1 + vr1
END SUB
SUB GETSCREEN
' GETSCREEN
SELECT CASE scren%
CASE 7
xmax = 320: ymax = 200: fract = 1.2
xchrmax = 40: ychrmax = 25
CASE 8
xmax = 640: ymax = 200: fract = 2.4
xchrmax = 80: ychrmax = 25
CASE 9
xmax = 640: ymax = 350: fract = 1.368
xchrmax = 80: ychrmax = 25
CASE 12
xmax = 640: ymax = 480: fract = 1
xchrmax = 80: ychrmax = 30
CASE ELSE
PRINT "SCREEN ERROR, not 7, 8, 9 or 12"
END SELECT
END SUB
SUB INITIALISE
G = 1
scren% = 12 ' screen mode
alpha = alpha0 ' alpha for test 13
stars = 3
' tdelta = .001 ' delta time 20/2/2016
dispcond = 100 ' display condition
initcond = 1 ' 1 = towards -1 = away
dist1 = 200
dist2 = 500
m(0) = 1000
m(1) = 1000
m(2) = 1
FOR i = 0 TO stars - 1
x0(i) = 0 ' initial positions
y0(i) = 0 ' initial positions
vx0(i) = 0 ' initial velocity
vy0(i) = 0 ' initial velocity
NEXT i
GETSCREEN ' Get screen parameters
alpha = alpha0 ' alpha for test 13
COORD m(0), m(1), dist1, x0(0), x0(1), vy0(0), vy0(1)
agrad = alpha * pi / 180
r1 = ABS(x0(1)): V1 = ABS(vy0(1))
x0(1) = r1 * COS(agrad): y0(1) = r1 * SIN(agrad)
vx0(1) = V1 * COS(agrad + pi / 2): vy0(1) = V1 * SIN(agrad + pi / 2)
r0 = ABS(x0(0)): v0 = ABS(vy0(0))
x0(0) = r0 * COS(agrad + pi): y0(0) = r0 * SIN(agrad + pi)
vx0(0) = v0 * COS(agrad + 1.5 * pi): vy0(0) = v0 * SIN(agrad + 1.5 * pi)
x0(2) = dist2 ' far away
dist20 = x0(2) - x0(0): dist21 = x0(2) - x0(1)
vx0(2) = -initcond * .9 * (SQR(2 * G * (m(0) / dist20 + m(1) / dist21)))
x0disp = xmax / 2: y0disp = ymax / 2
scalex = xmax / dist2 / 4: scaley = scalex
endc = 3 ' maximum numbers of revolutions
TT = pi * dist1 / V1
FOR i = 0 TO stars - 1
x00(i) = x0(i) ' initial positions
y00(i) = y0(i) ' initial positions
vx00(i) = vx0(i) ' initial velocity
vy00(i) = vy0(i) ' initial velocity
NEXT i
state = 0
END SUB
Back to my home page: Contents of This Document