' a star and planetary system governed by newton's laws.
' when the escape velocity of a planet is exceded, it's trajectory turns red
' the range of x values on the screen is -130 to +130
' the range of y values on the screen is -100 to +100

WINDOW -640, -512,640, 512

func gryb (yb, yp, ys, yt) = mp * (yp - yb) / rbp + ms * (ys - yb) / rbs + mt * (yt - yb) / rbt
func grxb (xb, xp, xs, xt) = mp * (xp - xb) / rbp + ms * (xs - xb) / rbs + mt * (xt - xb) / rbt
func gryp (yb, yp, ys, yt) = ms * (ys - yp) / rps + mb * (yb - yp) / rbp + mt * (yt - yp) / rpt
func grxp (xb, xp, xs, xt) = ms * (xs - xp) / rps + mb * (xb - xp) / rbp + mt * (xt - xp) / rpt
func grxs (xb, xp, xs, xt) = mp * (xp - xs) / rps + mb * (xb - xs) / rbs + mt * (xt - xs) / rst
func grys (yb, yp, ys, yt) = mp * (yp - ys) / rps + mb * (yb - ys) / rbs + mt * (yt - ys) / rst
func grxt (xb, xp, xs, xt) = mb * (xb - xt) / rbt + mp * (xp - xt) / rpt + ms * (xs - xt) / rst
func gryt (yb, yp, ys, yt) = mb * (yb - yt) / rbt + mp * (yp - yt) / rpt + ms * (ys - yt) / rst
PRINT "automatic momentum zero"
INPUT "central mass (50000)"; mb
xb = 0: yb = 0
INPUT "primary mass (2000)"; mp
INPUT "initial primary x (400)"; xp
INPUT "initial primary y (0)"; yp
zxp = -yp * (mb / (xp ^ 2 + yp ^ 2) ^ 1.5) ^ .5
zyp = xp * (mb / (xp ^ 2 + yp ^ 2) ^ 1.5) ^ .5
INPUT "circular orbit ? y/n"; acc$: IF acc$ = "y" THEN GOTO 10
PRINT "primary x-circular velo.="; zxp
INPUT "initial dx/dt"; zxp
PRINT "primary y-circular velo.="; zyp
INPUT "initial dy/dt"; zyp
10 INPUT "secondary mass (1000)"; ms
INPUT "initial secondary x (0)"; xs
INPUT "initial secondary y (300)"; ys
zxs = -ys * (mb / (xs ^ 2 + ys ^ 2) ^ 1.5) ^ .5
zys = xs * (mb / (xs ^ 2 + ys ^ 2) ^ 1.5) ^ .5
INPUT "circular orbit ? y/n"; acc$: IF acc$ = "y" THEN GOTO 20
PRINT "secondary x-circular velo.="; zxs
INPUT "initial secondary dx/dt"; zxs
PRINT "secondary y-circular velo.="; zys
INPUT "initial secondary dy/dt"; zys
20 INPUT "tertiary mass (500)"; mt
INPUT "initial tertiary x (-300)"; xt
INPUT "initial tertiary y (-400)"; yt
zxt = -yt * (mb / (xt ^ 2 + yt ^ 2) ^ 1.5) ^ .5
zyt = xt * (mb / (xt ^ 2 + yt ^ 2) ^ 1.5) ^ .5
INPUT "circular orbit ? y/n"; acc$: IF acc$ = "y" THEN GOTO 30
PRINT "tertiary x-circular velo.="; zxt
INPUT "initial tertiary dx/dt"; zxt
PRINT "tertiary y-circular velo.="; zyt
INPUT "initial tertiary dy/dt"; zyt
30 INPUT "final t (1600)"; tf
INPUT "number of intervals (800000)"; n
zxb = -(mp * zxp + ms * zxs + mt * zxt) / mb
zyb = -(mp * zyp + ms * zys + mt * zyt) / mb
h = tf / n
CLS
FOR q = 1 TO n
rbp = ((xb - xp) ^ 2 + (yb - yp) ^ 2) ^ 1.5
rbs = ((xb - xs) ^ 2 + (yb - ys) ^ 2) ^ 1.5
rps = ((xp - xs) ^ 2 + (yp - ys) ^ 2) ^ 1.5
rbt = ((xb - xt) ^ 2 + (yb - yt) ^ 2) ^ 1.5
rpt = ((xp - xt) ^ 2 + (yp - yt) ^ 2) ^ 1.5
rst = ((xs - xt) ^ 2 + (ys - yt) ^ 2) ^ 1.5
kxb1 = h * zxb: lxb = h * grxb(xb, xp, xs, xt)
kyb1 = h * zyb: lyb = h * gryb(yb, yp, ys, yt)
kxp1 = h * zxp: lxp = h * grxp(xb, xp, xs, xt)
kyp1 = h * zyp: lyp = h * gryp(yb, yp, ys, yt)
kxs1 = h * zxs: lxs = h * grxs(xb, xp, xs, xt)
kys1 = h * zys: lys = h * grys(yb, yp, ys, yt)
kxt1 = h * zxt: lxt = h * grxt(xb, xp, xs, xt)
kyt1 = h * zyt: lyt = h * gryt(yb, yp, ys, yt)
kxb2 = h * (zxb + .5 * lxb)
kyb2 = h * (zyb + .5 * lyb)
kxp2 = h * (zxp + .5 * lxp)
kyp2 = h * (zyp + .5 * lyp)
kxs2 = h * (zxs + .5 * lxs)
kys2 = h * (zys + .5 * lys)
kxt2 = h * (zxt + .5 * lxt)
kyt2 = h * (zyt + .5 * lyt)
nxb = xb + .5 * kxb1: nyb = yb + .5 * kyb1
nxp = xp + .5 * kxp1: nyp = yp + .5 * kyp1
nxs = xs + .5 * kxs1: nys = ys + .5 * kys1
nxt = xt + .5 * kxt1: nyt = yt + .5 * kyt1
rbp = ((nxb - nxp) ^ 2 + (nyb - nyp) ^ 2) ^ 1.5
rbs = ((nxb - nxs) ^ 2 + (nyb - nys) ^ 2) ^ 1.5
rps = ((nxp - nxs) ^ 2 + (nyp - nys) ^ 2) ^ 1.5
rbt = ((nxb - nxt) ^ 2 + (nyb - nyt) ^ 2) ^ 1.5
rpt = ((nxp - nxt) ^ 2 + (nyp - nyt) ^ 2) ^ 1.5
rst = ((nxs - nxt) ^ 2 + (nys - nyt) ^ 2) ^ 1.5
lxb = h * grxb(nxb, nxp, nxs, nxt)
lyb = h * gryb(nyb, nyp, nys, nyt)
lxp = h * grxp(nxb, nxp, nxs, nxt)
lyp = h * gryp(nyb, nyp, nys, nyt)
lxs = h * grxs(nxb, nxp, nxs, nxt)
lys = h * grys(nyb, nyp, nys, nyt)
lxt = h * grxt(nxb, nxp, nxs, nxt)
lyt = h * gryt(nyb, nyp, nys, nyt)
kxb3 = h * (zxb + .5 * lxb)
kyb3 = h * (zyb + .5 * lyb)
kxp3 = h * (zxp + .5 * lxp)
kyp3 = h * (zyp + .5 * lyp)
kxs3 = h * (zxs + .5 * lxs)
kys3 = h * (zys + .5 * lys)
kxt3 = h * (zxt + .5 * lxt)
kyt3 = h * (zyt + .5 * lyt)
nxb = xb + .5 * kxb2: nyb = yb + .5 * kyb2
nxp = xp + .5 * kxp2: nyp = yp + .5 * kyp2
nxs = xs + .5 * kxs2: nys = ys + .5 * kys2
nxt = xt + .5 * kxt2: nyt = yt + .5 * kyt2
rbp = ((nxb - nxp) ^ 2 + (nyb - nyp) ^ 2) ^ 1.5
rbs = ((nxb - nxs) ^ 2 + (nyb - nys) ^ 2) ^ 1.5
rps = ((nxp - nxs) ^ 2 + (nyp - nys) ^ 2) ^ 1.5
rbt = ((nxb - nxt) ^ 2 + (nyb - nyt) ^ 2) ^ 1.5
rpt = ((nxp - nxt) ^ 2 + (nyp - nyt) ^ 2) ^ 1.5
rst = ((nxs - nxt) ^ 2 + (nys - nyt) ^ 2) ^ 1.5
lxb = h * grxb(nxb, nxp, nxs, nxt)
lyb = h * gryb(nyb, nyp, nys, nyt)
lxp = h * grxp(nxb, nxp, nxs, nxt)
lyp = h * gryp(nyb, nyp, nys, nyt)
lxs = h * grxs(nxb, nxp, nxs, nxt)
lys = h * grys(nyb, nyp, nys, nyt)
lxt = h * grxt(nxb, nxp, nxs, nxt)
lyt = h * gryt(nyb, nyp, nys, nyt)
kxb4 = h * (zxb + lxb): kyb4 = h * (zyb + lyb)
kxp4 = h * (zxp + lxp): kyp4 = h * (zyp + lyp)
kxs4 = h * (zxs + lxs): kys4 = h * (zys + lys)
kxt4 = h * (zxt + lxt): kyt4 = h * (zyt + lyt)
xb = xb + (kxb1 + 2 * kxb2 + 2 * kxb3 + kxb4) / 6
yb = yb + (kyb1 + 2 * kyb2 + 2 * kyb3 + kyb4) / 6
xp = xp + (kxp1 + 2 * kxp2 + 2 * kxp3 + kxp4) / 6
yp = yp + (kyp1 + 2 * kyp2 + 2 * kyp3 + kyp4) / 6
xs = xs + (kxs1 + 2 * kxs2 + 2 * kxs3 + kxs4) / 6
ys = ys + (kys1 + 2 * kys2 + 2 * kys3 + kys4) / 6
xt = xt + (kxt1 + 2 * kxt2 + 2 * kxt3 + kxt4) / 6
yt = yt + (kyt1 + 2 * kyt2 + 2 * kyt3 + kyt4) / 6
zxb = kxb4 / h: zyb = kyb4 / h
zxp = kxp4 / h: zyp = kyp4 / h
zxs = kxs4 / h: zys = kys4 / h
zxt = kxt4 / h: zyt = kyt4 / h
IF q = n THEN GOTO 200
IF (INT(q / 400000) - q / 400000) <> 0 THEN GOTO 200
CLS
200 IF zxp * zxp + zyp * zyp > 2 * mb / (rbp ^ (1 / 3)) THEN cp = 4 ELSE cp = 9
IF zxs * zxs + zys * zys > 2 * mb / (rbs ^ (1 / 3)) THEN cs = 4 ELSE cs = 3
IF zxt * zxt + zyt * zyt > 2 * mb / (rbt ^ (1 / 3)) THEN ct = 4 ELSE ct = 6
PSET xb, yb, 5: PSET xp, yp, cp: PSET xs, ys, cs: PSET xt, yt, ct
NEXT q
END