Space Probe : Simple Gravitational Force Simulation

[Tosk]

Back in 1986 I coded a game for the ZX Spectrum that was published in Popular Computing Weekly. Remember that magazine, where you could type in Basic programs and perhaps get them to work?

In this thread I'll be working on converting this to BBC Basic as part of my relearning of this language.

The aim of the game is this:-

For input, the player puts in an initial X position, velocity and starting angle to try and set the Space Probe off towards the exit goal, which will be randomly placed on the far side of the screen (a simple line). If the space probe is plotted into the exit goal area then you've won the game. The graduated axis on the left is to aid the player on deciding that starting position - though I've yet to add some number against the divisions to help with that. The planet will also be placed randomly on the screen with a random size which in turn will affect it's mass. The game demonstrates the interaction of the gravitational force between two bodies.

For anyone else that may have published their BASIC programs during that time, there's two useful archives here:-

1. Index for Spectrum/QL typed in scripts : http://www.users.globalnet.co.uk/~jg27p ... _names.htm
2. This archive is from the "World of Spectrum" but the mag covers all scripts from different machines of the time. https://worldofspectrum.org/archive/mag ... ing-weekly

Here's the original mag and published page:-






I've been using the excellent Owlet editor on my Mac (I can type faster on a modern keyboard) and my first question is on the use of RND, which is very much prevalent in the code above which I'll need in the BBC script.

I've noticed in Owlet that when using RND(X) it ends up giving you the same number each time. On my actual BBC B it does give random numbers each time (pseudo random anyway). So how to test the code in Owlet? Is there something I can do about this or is there another editor I can use that doesn't have this limitation.

Here's my script so far (still needs work) which doesn't use RND at the moment because of the restriction above.

Code: Select all

MODE 1
CLS

PROC_Initialise

:REM Drawing Saturn (planet body) as a circle in the center
saturnX% = 640
saturnY% = 512
saturnRadius% = 240

PROC_DrawPlanet(saturnX%, saturnY%, saturnRadius%)

:REM Drawing Saturn's Rings using filled ellipses
GCOL 0, 3
FOR ringNum% = 1 TO 3
  ringWidth% = saturnRadius% + ringNum% * 70
  ringHeight% = 20 + ringNum% * 12
  
  PROC_DrawEllipse(saturnX%, saturnY%, ringWidth%, ringHeight%)
NEXT

:REM Initialize spacecraft's position and velocity for slingshot trajectory
spacecraftX% = 200
spacecraftY% = 750
vx% = 2
vy% = -3

:REM Simulate and draw spacecraft's path
GCOL 0, 7  :REM Color for spacecraft
FOR t% = 1 TO 2000
  PROC_DrawSpacecraft(spacecraftX%, spacecraftY%)
  
  :REM Compute gravitational force
  dx% = saturnX% - spacecraftX%
  dy% = saturnY% - spacecraftY%
  distance% = SQR(dx%^2 + dy%^2)
  
  :REM For simplicity, we'll use a scaled gravitational constant
  G% = 1E-3
  force% = G% / distance%^2
  
  :REM Adjust velocity based on gravitational force
  accelerationX% = (force% * dx% / distance%)
  accelerationY% = (force% * dy% / distance%)
  
  vx% = vx% + accelerationX%
  vy% = vy% + accelerationY%
  
  :REM Update spacecraft's position
  spacecraftX% = spacecraftX% + vx%
  spacecraftY% = spacecraftY% + vy%
  
  :REM IF distance% < saturnRadius% THEN EXIT
NEXT

STOP

DEF PROC_Initialise
  VDU 23,1,0;0;0;0;
  COLOUR 0
  COLOUR 1
ENDPROC

DEF PROC_DrawPlanet(h%, k%, r%)
  GCOL 0, 5
  MOVE h%, k%
  PLOT 145, r%, 0
ENDPROC

DEF PROC_DrawEllipse(h%, k%, hr%, vr%)
  MOVE h%, k%
  PLOT 0, hr%, 0
  PLOT 193, 0, vr%
ENDPROC

DEF PROC_DrawSpacecraft(x%, y%)
  PLOT 69, x%, y%
ENDPROC

[BigEd]

To set the random number generator to a random starting point you need to give it a negative number. To make that unpredictable, in Owlet, I think you need some kind of user interaction. For example

Code: Select all

10 PRINT "press any key":IF GET
20 IF RND(-TIME)
30 PRINT RND(1000000)

It's possible I've missed a trick, but I think the Owlet always starts in exactly the same state.

Nice idea to port this program!

[MartinB]

I think in Owlet, it’s just RND(0) that gives the same number each time whereas any unsigned integer in the brackets gives a different 0-integer random number each time. (Variables do work of course so if your example was literal and X is 0, you will get the same number but if X has a value, it’ll work as expected.)

[Tosk]

To set the random number generator to a random starting point you need to give it a negative number. To make that unpredictable, in Owlet, I think you need some kind of user interaction. For example

Code: Select all

10 PRINT "press any key":IF GET
20 IF RND(-TIME)
30 PRINT RND(1000000)

It's possible I've missed a trick, but I think the Owlet always starts in exactly the same state.

Nice idea to port this program!

This makes sense - had a feeling that Owlet was kind of cemented in an initial state each time you ran it. Thanks.

[Tosk]

In the meantime, I've recreated the other graphics set up for the game. As I'm coding I'm sort of appalled by all those GOSUB and GOTO in the ZX Spectrum Basic I used!

New developments tonight includes setting the target zone of where we are hoping to guide the space probe, the X-axis divisions on the left to give the user a sense of what number to input for the starting position of the space probe and even some random stars.

I'm using MODE 1 (320 x 200), but confused by Owlet's screen which seems to allow many more pixels. So at the moment I'm making the graphics larger to fit them in the Owlet screen.

Code: Select all

MODE 1
CLS

PROC_Initialise

:REM Set up target line variables

LET targetlinelen% =200
LET axis% = RND(1) 
LET targetlinestart% = RND(256-targetlinelen%)

: REM Draw target line - X or Y axis on the far side

IF axis% <.5 THEN MOVE  1200,targetlinestart% : DRAW 1200,targetlinestart%+targetlinelen% ELSE MOVE  targetlinestart%,0 : DRAW targetlinestart%+targetlinelen%,0

:REM Set up marked X axis and draw some stars for effect!

MOVE 0,0 : DRAW 0, 1000
FOR n = 1000 TO 0 STEP -100
    MOVE 0,n : DRAW 20,n
    PLOT 69,RND(1000),RND(1000)
NEXT n

:REM Drawing Saturn (planet body) as a circle in the center
saturnX% = 640
saturnY% = 512
saturnRadius% = 240

:REM Draw the planet
PROC_DrawPlanet(saturnX%, saturnY%, saturnRadius%)

:REM Drawing Saturn's Rings using filled ellipses
GCOL 0, 3
FOR ringNum% = 1 TO 3
  ringWidth% = saturnRadius% + ringNum% * 70
  ringHeight% = 20 + ringNum% * 12
  
  PROC_DrawEllipse(saturnX%, saturnY%, ringWidth%, ringHeight%)
NEXT

:REM Initialize spacecraft's position and velocity for slingshot trajectory
spacecraftX% = 200
spacecraftY% = 750
vx% = 2
vy% = -3

:REM Simulate and draw spacecraft's path
GCOL 0, 7  :REM Color for spacecraft
FOR t% = 1 TO 2000
  PROC_DrawSpacecraft(spacecraftX%, spacecraftY%)
  
  :REM Compute gravitational force
  dx% = saturnX% - spacecraftX%
  dy% = saturnY% - spacecraftY%
  distance% = SQR(dx%^2 + dy%^2)
  
  :REM For simplicity, we'll use a scaled gravitational constant
  G% = 1E-3
  force% = G% / distance%^2
  
  :REM Adjust velocity based on gravitational force
  accelerationX% = (force% * dx% / distance%)
  accelerationY% = (force% * dy% / distance%)
  
  vx% = vx% + accelerationX%
  vy% = vy% + accelerationY%
  
  :REM Update spacecraft's position
  spacecraftX% = spacecraftX% + vx%
  spacecraftY% = spacecraftY% + vy%
  
  :REM IF distance% < saturnRadius% THEN EXIT
NEXT

STOP

DEF PROC_Initialise
  VDU 23,1,0;0;0;0;
  COLOUR 0
  COLOUR 1
ENDPROC

DEF PROC_DrawPlanet(h%, k%, r%)
  GCOL 0, 5
  MOVE h%, k%
  PLOT 145, r%, 0
ENDPROC

DEF PROC_DrawEllipse(h%, k%, hr%, vr%)
  MOVE h%, k%
  PLOT 0, hr%, 0
  PLOT 193, 0, vr%
ENDPROC

DEF PROC_DrawSpacecraft(x%, y%)
  PLOT 69, x%, y%
ENDPROC

[SteveF]

I'm using MODE 1 (320 x 200), but confused by Owlet's screen which seems to allow many more pixels. So at the moment I'm making the graphics larger to fit them in the Owlet screen.

I haven't used Owlet, but just in case this is the source of confusion: on a real BBC, the OS presents a "virtual" screen resolution of 1280x1024 in all modes, so in mode 1 you need to multiply "physical" pixel coordinates by 4 when using them with commands like MOVE/DRAW/PLOT. The idea is that programs using graphics commands can be run in any screen mode rather than being tied to a specific mode's physical size.

Great to see your progress, I always enjoy a dev blog thread!

Edit: And apologies if this is just pedantry on my part, but FWIW, at least on a real UK machine, mode 1 is 320x256, not 320x200.

[Tosk]

I'm using MODE 1 (320 x 200), but confused by Owlet's screen which seems to allow many more pixels. So at the moment I'm making the graphics larger to fit them in the Owlet screen.

I haven't used Owlet, but just in case this is the source of confusion: on a real BBC, the OS presents a "virtual" screen resolution of 1280x1024 in all modes, so in mode 1 you need to multiply "physical" pixel coordinates by 4 when using them with commands like MOVE/DRAW/PLOT. The idea is that programs using graphics commands can be run in any screen mode rather than being tied to a specific mode's physical size.

Great to see your progress, I always enjoy a dev blog thread!

Edit: And apologies if this is just pedantry on my part, but FWIW, at least on a real UK machine, mode 1 is 320x256, not 320x200.

Great thanks for your help on that - the virtual screen resolution was something I completely had forgotten about.

[lovebug]

I love programs that simulate gravitational forces
great progress so far, well done
I cant wait to see more

[Tosk]

Now the code is in a much better position.

I had forgotten that variables with a % at the end were integers - which was screwing up some of the maths.

Another big clunker was the mass of the planet and spacecraft was left out of the math as well.

The only bit left is checking for success by using POINT and checking when the plotted spacecraft goes off screen.

Here's the link on Owlet:-

https://bbcmic.ro/#%7B%22v%22%3A1%2C%22 ... PROC%22%7D

Code: Select all

REM Space probe by Paul Bussey (C)2023

REM G% is the gravitational constant
G% = 10 : REM was scaled at 1E-3

saturnX% = 640
saturnY% = 512
saturnRadius% = 140
saturnMass% = saturnRadius% * 50

spaceprobeMass% = 1
spacecraftX% = 10
spacecraftY% = 600
dt = 1

v = 10 : REM Spacecraft velocity
a = -40 : REM Initial angle of space probe

MODE 1
CLS

VDU 23,1,0;0;0;0;
COLOUR 0
COLOUR 1

REM Set up the target line

LET targetlinelen% =200
LET axis = RND(1) 
LET targetlinestart% = RND(256-targetlinelen%)

REM Draw target line - either X or Y axis on the far side

IF axis <.5 THEN MOVE  1200,targetlinestart% : DRAW 1200,targetlinestart%+targetlinelen% ELSE MOVE  targetlinestart%,0 : DRAW targetlinestart%+targetlinelen%,0

REM Set up marked X axis and draw some stars for effect!

MOVE 0,0 : DRAW 0, 1000
FOR n = 1000 TO 0 STEP -100
MOVE 0,n : DRAW 20,n
PLOT 69,RND(1000),RND(1000)
NEXT n

REM Drawing a Saturn-like planet body as a circle in the center

:REM Draw the planet
PROC_DrawPlanet(saturnX%, saturnY%, saturnRadius%)

:REM Drawing Saturn's Rings using  ellipses
GCOL 0, 3
FOR ringNum% = 1 TO 3

    ringWidth% = saturnRadius% + ringNum% * 40
    ringHeight% = 20 + ringNum% * 12

    PROC_DrawEllipse(saturnX%, saturnY%, ringWidth%, ringHeight%)

NEXT

:REM Start calculating - initial velocity of space probe

a=a*(PI/180)

vx = v*COS(a)
vy = v*SIN(a)

REM Simulate and draw spacecraft's path

GCOL 0, 2  :REM Color for spacecraft
FOR t% = 1 TO 2000

    REM Compute gravitational force
    dx = saturnX% - spacecraftX%
    dy = saturnY% - spacecraftY%
    distance = SQR(dx^2 + dy^2)

    force = G% * saturnMass% * spaceprobeMass% / distance^2 

    Phi = ATN(ABS(dy/dx))

    fx=force*COS(Phi*SGN(dx))
    fy=force*SIN(Phi*SGN(dy))

    ax = fx/spaceprobeMass%
    ay = fy/saturnMass%

    sx=vx*dt+(ax*0.5)*dt^2
    vx=vx+ax*dt : 
    sy=vy*dt+(ay*0.5)*dt^2

    REM Plot the spacecraft

    PLOT 69, spacecraftX%, spacecraftY%

    REM Update spacecraft's position
    spacecraftX% = spacecraftX% + sx
    spacecraftY% = spacecraftY% + sy

    DRAW spacecraftX%,spacecraftY%

    REM Check for success

    REM PROC_CheckSuccess(axis,spacecraftX%,spacecraftY%)

    :REM IF distance% < saturnRadius% THEN EXIT

NEXT

STOP


DEF PROC_DrawPlanet(h%, k%, r%)
GCOL 0, 5
MOVE h%, k%
PLOT 145, r%, 0
ENDPROC

DEF PROC_DrawEllipse(h%, k%, hr%, vr%)
MOVE h%, k%
PLOT 0, hr%, 0
PLOT 193, 0, vr%
ENDPROC

REM DEF PROC_CheckSuccess(axis,x%,y%)
REM IF axis <0.5 THEN startloop% = targetlinestart% : endloop%= targetlinestart%+targetlinelen% ELSE startloop% = 1200+targetlinestart% : endloop%= 1200+targetlinestart%+targetlinelen%



ENDPROC

[BigEd]

Looks nice! I was expecting something more symmetrical though - unless you're in 3D or if you have a moving planet and are modelling a slingshot, wouldn't the exit path be like a mirror of the entry path?

Edit: just a thought, in case it's an artefact of the time step, could you run it twice with a half-size timestep and see how sensitive it is to that?

[Tosk]

Looks nice! I was expecting something more symmetrical though - unless you're in 3D or if you have a moving planet and are modelling a slingshot, wouldn't the exit path be like a mirror of the entry path?

Edit: just a thought, in case it's an artefact of the time step, could you run it twice with a half-size timestep and see how sensitive it is to that?

Thanks! I'll add an explanation of the game which will help explain the design and I'll copy this in to the main post. At the moment I haven't added inputs for the player - they are set as constants at the beginning.

For input, the player puts in an initial X position, velocity and starting angle to try and set the Space Probe off towards the exit goal, which will be randomly placed on the far side of the screen. The graduated axis on the left is to aid the player on deciding that starting position - though I've yet to add some number against the divisions to help with that.

I think from memory (I was 17 at the time!) the time step is the variable DT, which at the moment is set 1, when you increase that you see the plots are calculated further apart. The reason I'm drawing a line between the plots is because it allows for using detection (POINT) to check if the probe is in the goal exit area.

[Tosk]

Now with detection if space probe is out of limits or has reached the goal area. Because of these tests the code runs significantly slower, so I've increased DT to 4 (time step) to speed things up.

If you have ideas on how to speed things up - let me know!

Just need to add an introduction screen and invite the user to input the Space Probe Y Axis start position, start angle direction and velocity.

https://bbcmic.ro/#%7B%22v%22%3A1%2C%22 ... PROC%22%7D

[joachim]

I feel like there are some bugs in the physics here, which likely explains BigEd's concern (can't tell whether they were also in the Spectrum original):

Code: Select all

970 ax = fx/spaceprobeMass%
980 ay = fy/saturnMass%

Surely you'd want to divide both co-ordinates by spaceprobeMass%.

Code: Select all

1000 sx=vx*dt+(ax*0.5)*dt^2
1010 vx=vx+ax*dt : 
1020 sy=vy*dt+(ay*0.5)*dt^2

vy never gets updated.

Code: Select all

940 fx=force*COS(Phi*SGN(dx))
950 fy=force*SIN(Phi*SGN(dy))

This doesn't actually set the sign of fx the way you want, because COS(-Phi) is equal to COS(Phi) whereas you need -COS(Phi). Simplest solution is to move the brackets so that it becomes `fx=force*COS(Phi)*SGN(dx)`. (fy happens to be correct already because SIN(-Phi) = -SIN(Phi).)

[ChrisB]

If you have ideas on how to speed things up - let me know!

Some suggestions:
Current time to perform a check 250ms

In MODE 1 each pixel is 4 OS units wide. At the moment you're checking each pixel 4 times! So increasing checkpoint by 4 gives 64ms

Your check iterates over the axis check. Splitting this out to two separate procedures gives 35ms

On success you set the exit condition and force the loop to exit early with 'PRINT "Mission Success!" : checkpoint%= endloop% : exitspace = -1' but the loop condition is 'exitspace = -1 OR checkpoint% = endloop%'. Given the success is only true once there is not point exiting the loop early. Simplifying these to set exitspace on success and not checking this on the loop condition but status to print success gives 27ms

Now we have simplified the loop we can change it from being a REPEAT UNTIL with an increment within the loop to a FOR NEXT. This give 16ms

The code now looks as follows:
https://bbcmic.ro/#%7B%22v%22%3A1%2C%22 ... PROC%22%7D

However - the POINT command in quite slow and from what I can tell you're only looking to see if the line has crossed another line - so this can be achieved with a couple of IF statements rather than a loop over each point in the line then the check becomes around 3ms and your check no longer takes a significant amount of time.

https://bbcmic.ro/#%7B%22v%22%3A1%2C%22 ... %5Cn%22%7D

[Tosk]

I feel like there are some bugs in the physics here, which likely explains BigEd's concern (can't tell whether they were also in the Spectrum original):
...

Thanks for casting your eye on this joachim. I'll see to get these items sorted. It's mixture of copy and paste mishaps, having trouble reading some of the scanned text from the original magazine scan and for me to practice attention to detail!

[Tosk]

Some suggestions:
Current time to perform a check 250ms.....

Thanks Chris for these ideas to reduce the check times - fantastic.

I'm very interested to learn how you got these timings. What have you used to measure the timings of these chunks of code?

[Tosk]

This latest now includes Chris's performance improvements and joachim's corrections to the maths.

All in all - it makes the game a bit harder!

I also noticed that vx and vy were being calculated after sx and sy, these two calculations should have been done before - that was also a crucial bug fix.

https://bbcmic.ro/#%7B%22v%22%3A1%2C%22 ... %5Cn%22%7D

[ChrisB]

To measure time use the TIME counter that increments 100 times per second. E.g.

Code: Select all

T%=TIME
REM Do some other stuff
PRINT TIME-T%

Looking back I've quoted everything in miliseconds - but actually they are centiseconds. Now - with apologies for the long post. I got carried away...

You've got a mix of integer and float variables that might be causing errors to creep in. For example if you move from dt=2 to dt=1 then the final position of the spacecraft changes quite significantly.

Let's see what else we can do to speed it up.

So we'll put that TIME construct around the main repeat loop and we get 198ms per "frame" - I've actually counted the frames and taken an average so those actually are milliseconds

Continuing on the theme of types "axis" is only a flag. Changing it to a value of 1 or 2 and making it integer reduces frame time to 196ms. The "Mission success" print is also checked every frame and only needs to be displayed at the end so let's move it to after the UNTIL. This moves s down to 195ms. Now these are only small - but you can see how they all add together. Doing the same with exitspace and changing the UNTIL to be <0 brings down to 193ms. But let's go after some bigger gains.

Floating point maths is expensive so doing it where unnecessary will hurt performance. So for these two lines

Code: Select all

distance = SQR(dx^2 + dy^2)
force = G% * saturnMass% * spaceprobeMass% / distance^2 

We see a square root function - only to square the result. We can simplify that to

Code: Select all

distance = dx^2 + dy^2
force = G% * saturnMass% * spaceprobeMass% / distance

This bring the frame time to 180ms.

But what about spaceprobeMass%? Objects that are different masses fall at the same rate in the same gravitational field and all that - and you can see this by increasing spaceproemMass% and the resulting path is the same. So we see this code:

Code: Select all

force = G% * saturnMass% * spaceprobeMass% / distance 
Phi = ATN(ABS(dy/dx))
fx=force*COS(Phi*SGN(dx))
fy=force*SIN(Phi*SGN(dy))
ax = fx/spaceprobeMass%
ay = fy/spaceprobeMass%

And you multiply by spaceprobeMass% only to divide by it later. We can remove these (and the fx and fy variables) and go straight to ax/ay. So now we're down to 169ms.

There's also G% * saturnMass%. That's always the same so we can calculate that at the start of the loop and just use it as a constant. 168ms.

We can also look at the dt^2 in the sx/sy calculation. dt never changes so we can pre-calculate this as well. 162ms.

Now we get into the esoterics of BASIC. It turns out that ^2 is a lot slower then multiplying together so we can change the distance calculation to be
distance = dx*dx + dy*dy - Now we're at 158ms. About 0.04 seconds per frame faster that we started with.

Since you're drawing frame to frame the PLOT 69 is unnecessary and a MOVE at the start of the process can replace it. That brings us to 153ms.

We could now consider this "fast enough" to move to dt of 1 - which means we could drop the dt calculations entirely. This now takes us to 142ms. We can remove sx as it's only used as a temporary variable and we get to 140ms. The next steps I'd be looking at woud be lookup tables for the SIN/COS but here's the code as is
https://bbcmic.ro/#%7B%22v%22%3A1%2C%22 ... %5Cn%22%7D

[joachim]

We could now consider this "fast enough" to move to dt of 1 - which means we could drop the dt calculations entirely. This now takes us to 142ms. We can remove sx as it's only used as a temporary variable and we get to 140ms. The next steps I'd be looking at woud be lookup tables for the SIN/COS but here's the code as is
https://bbcmic.ro/#%7B%22v%22%3A1%2C%22 ... %5Cn%22%7D

There's a simpler speedup still available, as you can get rid of all the trig functions at the cost of bringing back one SQR. The key is that give or take some signs, COS(ATN(dy/dx)) simplifies to dx / SQR(dx*dx + dy*dy) and similarly SIN(ATN(dy/dx)) is dy / SQR(dx*dx + dy*dy).

So the calculation of ax and ay can be replaced with (pre-calculating the square root so that it doesn't get computed twice):

Code: Select all

sd=SQR(distance)
ax=force*dx/sd
ay=force*dy/sd

I claim this is the correct calculation, but if you try this you'll notice it behaves differently from the current code. That's because the signs in the calculation of ax are still wrong in the current code . So here is a version where I have made my changes but put an erroneous ABS back in so that the physics is unchanged (up to floating point inaccuracies) from ChrisB's last post, thus ensuring a like-for-like comparison. It runs in 83ms/frame — close to twice as fast as ChrisB's last version.

Meanwhile, here is the same with my preferred physics in which I have:

• corrected the sign of ax
• fixed a typo in line 990 (vy should be incremented by ay, not ax)
• changed a sign in 1020 and 1030 because you probably want to move by the average of the old v and the new v, not by the average of the new v and some extrapolated future value (this is a comparatively small change, i.e. first-order rather than zeroth-order, and you can revert it if you don't like it)

[Tosk]

Thanks both!

With the corrected maths, typos and other mistakes out of the way, together with some optimisation, I've now added:-

- Inputs for the initial angle of flight, velocity and Y value for spacecraft
- Validation if the inputs are out of range
- An option for the player to keep the position of the planet position, size and goal exit area to have another go.

[EDIT]

Final version.

1) If you failed the validation input, then the validation flag or message was not reset. This bug now fixed.
2) You now get a message of success if you reach the goal exit.
3) The Y start position is also printed at the top in the stats area.

Link to Owlet editor:-

https://bbcmic.ro/#%7B%22v%22%3A1%2C%22 ... PROC%22%7D