'PLOTSTAR.bas has the goal of propagating the orbit using direct newtonian
integration.
'It is of questionable accuracy.
'Keypad arrow keys control x, y velocity, type q to quit.
'Screen and background setup
30 SCREEN 2
'big fat and white, poor aspect not proportional
'Numbers for full scale default =640x=across, 180 down:
360 'WINDOW (0, 0)-(400, 300)
370 'LINE (10, 3)-(620, 190), 1, B
'1 is color, B is box
380 CIRCLE (320, 100), 20, 4
'location radius color
'Background done Units are display units not engineering
units.
'Force center is at origin
'Position coordinates and velocity for circular orbit
'x = 330 Comet
x = 90 'Planet
y = 0
vx = 0
'vy=.05
vy = .95
xp = 0
yp = 0
'Type z to break loop
WHILE a$ <> "z"
'Paint the point at x, y
PSET (sx, sy)
' Delay loop will be removed when
calculations get longer.
FOR i = 1 TO 30: k = SIN(2):
NEXT i
'Keyboard input processing adds or subtracts 0.1 unit of acceleration
a$ = INKEY$
IF a$ = "6" THEN vx = vx
+ .1
IF a$ = "4" THEN vx = vx
- .1
IF a$ = "8" THEN vy = vy
- .1
IF a$ = "2" THEN vy = vy
+ .1
IF a$ = "q" THEN END
'Save previous value to
enable erasure.
sxp = sx
syp = sy
'Determine central force
using r with GM=1
rr = x * x + y * y
r = SQR(rr)
sr = x * x + 4 * y * y
IF sr < 600 THEN PRINT
"Crashed": END
f = 100 / rr
fx = f * x / r
fy = f * y / r
'Integrate central force
acceleration
vx = vx - fx
vy = vy - fy
'integrate velocity
x = x + vx
y = y + vy
'Screen coordinate generation
with wrap limits
sx = (ABS(x) + 320) MOD
640
sx = 320 - 320 * SGN(x)
+ sx * SGN(x)
sy = (ABS(y) + 100) MOD
200
sy = 100 - 100 * SGN(y)
+ sy * SGN(y)
'Erasure is done after the maximum
time to brighten display
PRESET (sxp, syp)
WEND
END
'LOSTinSPACE.bas demonstrates that using too much thrust will
exceed the escape velocity.
'Keypad arrow keys control x, y velocity, type q to quit.
'Screen and background setup
30 SCREEN 2
'big fat and white
'Numbers for full scale default =640x=across, 180 down:
360 'WINDOW (0, 0)-(400, 300)
370 LINE (10, 3)-(620, 190), 1, B
'1 is color, B is box
380 CIRCLE (320, 100), 20 'location radius
color
'Background done Units are display units not engineering
units.
'Position coordinates, force center is at initial x,y
x = 330
y = 0
vx = 0
vy = .05
xp = 0
yp = 0
'Type z to break loop
WHILE a$ <> "z"
'Paint the point at x, y
PSET (sx, sy)
' Delay loop will be removed when
calculations get longer.
FOR i = 1 TO 30: k = SIN(2):
NEXT i
'Keyboard input processing adds or subtracts one unit of acceleration
a$ = INKEY$
IF a$ = "6" THEN vx = vx
+ .1
IF a$ = "4" THEN vx = vx
- .1
IF a$ = "8" THEN vy = vy
- .1
IF a$ = "2" THEN vy = vy
+ .1
IF a$ = "q" THEN END
'Save previous value to
enable erasure.
sxp = sx
syp = sy
'Determine central force
using r with GM=1
rr = x * x + y * y
r = SQR(rr)
IF r < 15 THEN PRINT
"Crashed": END
f = 100 / rr
fx = f * x / r
fy = f * y / r
'Integrate central force
acceleration
vx = vx - fx
vy = vy - fy
'integrate velocity
x = x + vx
y = y + vy
'Screen coordinate generation
with wrap limits
sx = (ABS(x) + 320) MOD
640
sx = 320 - 320 * SGN(x)
+ sx * SGN(x)
sy = (ABS(y) + 100) MOD
200
sy = 100 - 100 * SGN(y)
+ sy * SGN(y)
'Erasure is done after the maximum
time to brighten display
PRESET (sxp, syp)
WEND
END
5 'ORBIT5.BAS translated 11/27/96 from the Commodore 64. Display is
not fixed.
10 REM ORBITB CALCULATES ORBITS IN 2DIM AND ALLOWS MODIFICATION
11/8/87
20 REM STARTING SET UP
30 E = 0: NU = 0: W = 0
40 DT = 10: PI = 3.141596
50 MU = 3.98601E+14: REM MU=GM FOR EARTH
60 P = 6000000!: REM P IS SEMIMINOR AXIS
70 H = 4.8E+10: REM H IS SPECIFIC ANGULAR MOMENTUM = RXV
80 S = 1E+07: REM S IS SCREEN SCALE FACTOR
90 REM Commodore 64 leftovers BIT MAP GRAPHICS SET UP
100 'PIX=8192
110 REM PUT BIT MAP AT 8192
120 'POKE53272,PEEK(53272)OR8
130 REM ENTER BIT MAP MODE
140 'POKE53265,PEEK(53265)OR32
150 REM CLEAR BIT MAP
160 FOR I = PIX TO PIX + 7999: NEXT I'POKEI,0
170 REM SET COLOR TO CYAN AND BLACK
180 FOR I = 1024 TO 2023: NEXT I: 'POKEI,3
190 REM END SET UP
200 REM BEGIN LOOP
210 REM ORBITAL MOVEMENT ITERATION
220 R = P / (1 + E * COS(NU))
230 DN = H / R / R
240 NU = NU + DN * DT: IF NU > 2 * PI THEN NU = NU - 2 * PI
250 REM SETING UP OUTPUT
260 TH = W + NU: REM ANGLE FROM X AXIS TO R
270 REM X IS TO RIGHT, Y IS DOWN, TH IS MEASURED FROM X
GOING CW
280 X = R * COS(TH): Y = R * SIN(TH)
400 REM-----------------PLOT SUBROUTINE
410 XP% = 159 + INT(159 * X / S)
420 YP% = 99 + INT(99 * Y / S)
430 IF XP% < 0 OR XP% > 319 THEN S = S * 1.1: GOTO 400
440 IF YP% < 0 OR YP% > 199 THEN S = S * 1.1: GOTO 400
450 REM CHARACTER POSITION IN HORIZ LIN
460 CH = INT(XP% / 8)
470 RO = INT(YP% / 8): REM ROW
480 REM HORIZONTAL LINE IN CHARACTER
500 LN = YP% AND 7
510 REM BIT MAP ADDRESS
520 BY = PIX + RO * 320 + 8 * CH + LN
530 REM ONES COMPLIMENT OF X ADDRESS FLIPS LINE AROUND
540 BI = 7 - XP% AND 7
550 REM ONES ALREADY PRESENT ARE KEPT
560 POKE BY, PEEK(BY) OR 2 ^ BI
600 REM----INPUT ORBIT VELOCITY CHANGES
605 'A$ = "" failed.
610 A$ = INKEY$: IF A$ = "Q" THEN GOTO 980
620 IF A$ <> "C" THEN GOTO 200
640 VN = H / R: VR = E * SIN(NU) * SQR(MU / P) * (1 - E * E) ^ .25:
REM ANGULAR AND RADIAL VELOCITY
660 A$ = INKEY$: IF A$ = "" THEN GOTO 660
670 IF A$ = "+" THEN VN = VN * 1.01
680 IF A$ = "-" THEN VN = VN * .99
690 REM INPUT IS DONE. UPDATE ELEMENTS
740 H = R * VN
750 P = H * H / MU
760 REM MAKING E
770 V2 = VN * VN + VR * VR
780 D0 = R * VR
800 ER = V2 * R / MU - 1 - D0 * VR / MU
810 EN = -D0 * VN / MU
820 E = SQR(ER * ER + EN * EN)
840 REM E IS DONE. NOW MAKE W THE ANGLE OF PERIGEE WITH X AXIS
850 IF E = 0 THEN W = 0: GOTO 880
860 IF ER = 0 THEN W = TH - (EN > 0) * PI / 2: GOTO 880
865 AT = ATN(EN / ER): IF ER < 0 THEN AT = PI + AT
870 W = TH + AT: IF W > 2 * PI THEN W = W - 2 * PI
880 REM MAKING THE NEW NU
890 NU = TH - W
910 REM PRINT"R P TH NU W E H VX VY AN",R;P;TH;NU;W;E;H;VX;VY;AN
950 REM A$ =inkey$:IF A$="" THEN GOTO 950
969 IF A$ = "Q" THEN GOTO 980
970 GOTO 200
980 REM THIS UNDOES THE BIT MAP
990 POKE 53265, PEEK(53265) AND 223
1000 POKE 53272, PEEK(53272) AND 247
1010 REM NORMAL DISPLAY
1020 END