SyntaxBomb  Indie Coders
Languages & Coding => SmallBASIC => Topic started by: lettersquash on March 27, 2020, 12:55:22 AM

Here's a fun little demo of the smallBASIC graphics  a simulation of stars / planets in mutual orbit around each other. You can mess with the formulae and get some crazy results. Bodies are drawn as circles with size proportional to mass. I love the way you get tiny ones whizzing about, and now and then something gets ejected out of the system.
If the swarm drifts off towards the edge of the window, you can pan about in any direction by clicking on the window towards the place you want to see more of, and it will scroll that way proportionally to how far from the centre the mouse is (click and hold to keep moving that way). The + and  buttons zoom in and out, but have odd effects on the physics and speed. Tested on Windows 7, 64bit.
REM SmallBASIC
REM created: 23/03/2020
'By lettersquash. A sort of gravity simulation. Have fun with the numbers.
'pen on 'Doesn't seem to need this.
option predef antialias off
n=50
dim x (n), y (n), dx (n), dy (n), m (n)
border=xmax50
mix=xmax\2:miy=ymax\2
rect border,0,xmax,ymax,3 filled
rect border+2,miy150,xmax2,miy100,2 filled
rect border+12,miy126,xmax12,miy124,7 filled
rect border+24,miy138,border+26,miy112,7 filled
rect border+2,miy98,xmax2,miy48,2 filled
rect border+12,miy74,xmax12,miy72,7 filled
for i=1 to n
' initial positions
x(i)=rnd*400+mix200
y(i)=rnd*400+miy200
'initial velocity (perhaps rotation)
' dx(i)=(y(i)miy)*0.005
' dy(i)=(mixx(i))*0.005
m(i)=rnd*4+0.1
next
repeat
for i=1 to n
'Blank out body
if x(i)<border3 then circle x(i),y(i),m(i),1,0 filled
for j=1 to n
'Gravitational forces caused by each body are summed as x,y components
if i<>j then
'x,y components of distance between the two
a=x(j)x(i):b=y(j)y(i)
'and distance squared (don't need the root)
d2=(a^2+b^2)
'approximation of gravitational equation for dx,dy
'Fgrav = m1*m2/d2, but using product too much
dx(i)+=m(j)*0.02*a/d2
dy(i)+=m(j)*0.02*b/d2
fi
next j
'change position (experimented with a bit of 'drag' with weird results)
x(i)+=dx(i)/m(i)':dx(i)*=0.998
y(i)+=dy(i)/m(i)':dy(i)*=0.998
if pen(3) then click ' allow pan/zoom via LMB
if x(i)<border3 then circle x(i),y(i),m(i),1,i%16 filled 'i%16 'draw bodies
next i
showpage
delay 1
' if inkey=" " then delay 100:pause
until 0
sub click
p4=pen(4):p5=pen(5)
if p4>border+2 and p4<xmax2 and p5>miy150 and p5<miy100 then
'zoom in with "+" button
x(i)=mix(mixx(i))*1.05
y(i)=miy(miyy(i))*1.05
elseif p4>border+2 and p4<xmax2 and p5>miy98 and p5<miy48 then
'zoom out with "" button
x(i)=mix(mixx(i))*0.95
y(i)=miy(miyy(i))*0.95
else
'"Pan"  move screen proportional to mouse distance from centre with LMB down
'Would be good with normal drag to move, but couldn't figure it out.
px=p4mix:py=p5miy
x(i)=(px/100):y(i)=(py/100)
fi
end sub

Hey LS, that's pretty cool!
Kinda unexpected reaction to mouse location.
Is there some kind of center of gravity you have to calculate for all the bodies to orbit (or not)?

Yeah the mouse thing is a bit weird at first and might be better with a more normal clickdragrelease thing you do with documents, maps etc., although when you get used to it and remember that it scrolls at a speed and direction dictated by the distance between the mouse cursor (held down of course) and the centre of the screen, it means you can follow objects (especially  if you have it  using a click lock  I just doubleclickhold on my laptop's trackpad and then I can move the mouse around without holding down the physical left mouse button all the time).
The gravity equation is very much approximated and I'm not sure it's anything like correct, but it seems to be giving believable results. I've often tried to get this to work, but never had the speed of hardware or software (again, sB is surprisingly fast for calculations and plotting to the screen), and then it's a lot of guesswork, intuition and tweaking. I might have a go at writing the correct formula at some point and see how it compares. The correct method, I think, would be to work out for each object the vector of the force (direction and strength) of gravity caused by each of the other objects, which would combine to be one single vector, then apply that to the dx() and dy() value. Those are velocities in the x and y directions, of course, and just added to the x() and y() array values each time round the loop. Newton's  erm  second? law says that just keeps happening  object moving (in a vacuum) keeps going the same velocity unless acted on by a force.
However, because combining all those vectors is slow (and I often get stuck with it when using cosines and sines or whatever to work out the direction), and because vectors can be divided up into their x and y components, I just work out those components (or a 'quick and cheap' approximation) for each interaction between the (i) index object and the (j) object as separate x and y components and alter their dx and dy component speeds. I did that first ignoring the masses, and was surprised to find I'd cracked it at all, then added the mass involvement  which is even more guesswork and trial and error than the first part. The gravitational force between two objects is M1*M2/distance^2, where M1 and M2 are their masses, but putting both Ms in seemed too strong. Then, since the actual acceleration of the object is inversely proportional to its own mass (with no effect from the others, of course)  F=ma, so a=F/m, that's got to go into the equations somewhere. Happily, dividing the velocity components as they're updated (rather than the acceleration) seemed to give pleasing results...I've no idea why!
Small (lowmass) bodies do indeed get whizzed around fast and slingshot more easily into "outer space". Then, because they're less massive, they often return reasonably easily because the force of the combined gravity of the bunch can have a big effect on them (a=F/smallm). Heavy bodies take longer to be dragged around, but are also generally more stable to those highspeed ejections from a system (although that can still happen if they're very close as they pass another object). Of course, for the purposes of the calculations, they're just points, despite being drawn as circles (and it's faster if you change circle to pset). Another hope with this idea is/was to allow lumps of stuff to actually combine into a new object, combining their mass. Then one could model the formation of heavenly bodies from dust and debris.
Thanks for asking  sorry about the essay response!

Wow LS!
Thanks for taking the time for "lecture". You have me thinking... :)
BTW there is a collision of galaxies in SB code library, you might like to compare, I might also ;)

Wow LS!
Thanks for taking the time for "lecture". You have me thinking... :)
Then my work here is done. ;)
BTW there is a collision of galaxies in SB code library, you might like to compare, I might also ;)
Oh that's very interesting  I can't find it though  what's it called please?

There are at least 2 one is g_col for palm but that might be corrupted, it did not work for me.
The other is this:
https://raw.githubusercontent.com/smallbasic/smallbasic.samples/master/graphics%201/g_col_09.bas
'
'';Galaxy Collision Program ver 0.9
'';Keijo Koskinen 20/4/2005 for smallbasic
'';Displays collision, views from XY and XZ direction
'';core algorithm modified from old Galaxy Collision
'';Astronomy magazine Dec 88 page 92.
'';by M.C.Schroeder and Neil F. Comins,
''
DIM Xs(2,1000), Ys(2,1000), Zs(2,1000) ''; target gal=(1,) intruder=(2,) coordinates
dim VXs(2,1000), VYs(2,1000), VZs(2,1000) ''; speed
dim NRR(2),NRS(2),NS(2),DR(2)
CLS
color 10,0
''; TARGET GALAXY
NRR(2) = 23 ''; orbits in target gal
NRS(2) = 10 ''; stars / orbit in target galaxy
NS(2) = NRR(2) * NRS(2) ''; Total of stars in target galaxy
DR(2) = 20 / (NRR(2)  1)
M2 = 6 ''; Target galaxy mass
X2 = 200 ''; initial Target galaxy coordinates
Y2 = 90
Z2 = 10
''; *** Try here different speed factor to greate different action
''; VX2=0.085 ; distroys totally galaxies
''; VY2=0
''; VX2=0.085 ''; really hot relationship in the beginning
''; VY2=0.065
''; VX2= .085 ; galazies keep fighting over 1000 million years
''; VY2=0.29
''; VX2= 0.4 ''; galazies keep dancing over and over
''; VY2=0.4
''; VX2=0.116 ''; What a tension
''; VY2=0.41
VX2=0.116
VY2=0.52
'';
VZ2 = 0
''; INTRUDER GALAXY
NRR(1) = 23 ''; rings in intruder galaxy
NRS(1) = 10 ''; stars/ring in intruder galaxy
NS(1) = NRR(1) * NRS(1) ''; Total of stars in intruder galaxy
DR(1) = 20 / (NRR(1)  1)
M1 = 6 ''; Intruder galaxy mass
X1 = 150
Y1 = 100 ''; Initial coordinates of intruder galaxy
Z1 = 0
VX1 = 0 ''; Intruder galaxy initial speed
VY1 = 0
VZ1 = 0
''; common factors
SF2 = 2 ''; softening factor
k=0.1 ''; size factor
gsize=k*xmax/((NRR(1)+NRR(2))/2)
NTSPR = 2000 ''; Time / millions of years
TIIME = 0
''; place stars in orbits in target galaxy
''; from r=4 and p at intervals of dr
I=1
FOR IR = 1 TO NRR(1)
R = 4 + (IR/1  1) * DR(1)
V = SQR(((M1*DR(1))/1) / R)
TH = (.5 * V / R) * (180 / 3.14159)
IF R <10 THEN V = .99 * V
if R<7 then V=.98*V
if R <5 then V=0.95*V
FOR IT = 1 TO NRS(1)
T = rnd * 360
I = I + 1
''; initialize stars positions
Xs(1,I) = R * COS(T / 57.2958) + X1
Ys(1,I) = R * SIN(T / 57.2958) + Y1
Zs(1,I) = Xz1+0
VZs(1,I) = VZ1+0
''; init star velocities in orbits
VXs(1,I) = V * SIN(T/57.2958) +VX1
VYs(1,I) = V * COS(T/57.2959) +VY1
NEXT IT
NEXT IR
''; place stars in orbits in intruder galaxy
''; from r=4 and up at intervals of dr
I=1
FOR IR = 1 TO NRR(2)
R = 4 + (IR/1  1) * DR(2)
V = SQR(((M2*DR(2))/1) / R)
TH = (.5 * V / R) * (180 / 3.14159)
IF R <9 THEN V = .93 * V
if R<7 then V=.92*V
if R <5 then V=0.9*V
FOR IT = 1 TO NRS(2)
T=rnd*360
I = I + 1
''; initialize stars positions
Xs(2,I) = R * COS(T / 57.2958) + X2
Ys(2,I) = R * SIN(T / 57.2958) + Y2
Zs(2,I) = Z2+0
VZs(2,I) = VZ2+0
'';init star velocities in orbits
VXs(2,I) = V * SIN(T/57.2958)+VX2
VYs(2,I) = V * COS(T/57.2958) +VY2
NEXT IT
NEXT IR
scr_update
''; particle pusher routine
label timesteps
FOR K = 1 TO NTSPR ''; ok lets go
''; determine the forces pulling stars in intruder gal
J = 2
FOR I = 1 TO NS(J)
R1 = M1 / ((Xs(J,I)  X1) ^ 2 + (Ys(J,I)  Y1) ^ 2 + (Zs(J,I)  Z1) ^ 2 + SF2) ^ 1.5
R2 = M2 / ((Xs(J,I)  X2) ^ 2 + (Ys(J,I)  Y2) ^ 2 + (Zs(J,I)  Z2) ^ 2 + SF2) ^ 1.5
''; calculate stars x,y,z accelerations
AX = R1 * (X1  Xs(J,I)) + R2 * (X2  Xs(J,I))
AY = R1 * (Y1  Ys(J,I)) + R2 * (Y2  Ys(J,I))
AZ = R1 * (Z1  Zs(J,I)) + R2 * (Z2  Zs(J,I))
''; update star positions and velocities
VXs(J,I) = VXs(J,I) + AX
VYs(J,I) = VYs(J,I) + AY
VZs(J,I) = VZs(J,I) + AZ
Xs(J,I) = Xs(J,I) + VXs(J,I)
Ys(J,I) = Ys(J,I) + VYs(J,I)
Zs(J,I) = Zs(J,I) + VZs(J,I)
NEXT I
''; determine the forces pulling stars in target gal
J = 1
FOR I = 1 TO NS(J)
R1 = M1 / ((Xs(J,I)  X1) ^ 2 + (Ys(J,I)  Y1) ^ 2 + (Zs(J,I)  Z1) ^ 2 + SF2) ^ 1.5
R2 = M2 / ((Xs(J,I)  X2) ^ 2 + (Ys(J,I)  Y2) ^ 2 + (Zs(J,I)  Z2) ^ 2 + SF2) ^ 1.5
''; calculate stars x,y,z accelerations
AX = R1 * (X1  Xs(J,I)) + R2 * (X2  Xs(J,I)) ''kokeilu nurin x1 ja x2
AY = R1 * (Y1  Ys(J,I)) + R2 * (Y2  Ys(J,I)) '' samoin
AZ = R1 * (Z1  Zs(J,I)) + R2 * (Z2  Zs(J,I)) '' kuin myos
''; update star positions and velocities
VXs(J,I) = VXs(J,I) + AX
VYs(J,I) = VYs(J,I) + AY
VZs(J,I) = VZs(J,I) + AZ
Xs(J,I) = Xs(J,I) + VXs(J,I)
Ys(J,I) = Ys(J,I) + VYs(J,I)
Zs(J,I) = Zs(J,I) + VZs(J,I)
NEXT I
''; update postions and velocities of the galactic centers
RR = ((X1  X2) ^ 2 + (Y1  Y2) ^ 2 + (Z1  Z2) ^ 2 + SF2) ^ 1.5
AX = (X2  X1) / RR: AY = (Y2  Y1) / RR: AZ = (Z2  Z1) / RR
VX1 = VX1 + M2 * AX: VY1 = VY1 + M2 * AY: VZ1 = VZ1 + M1 * AZ ''; +M1*AZ
VX2 = VX2  M1 * AX: VY2 = VY2  M1 * AY: VZ2 = VZ2  M1 * AZ
X1 = X1 + VX1
Y1 = Y1 + VY1
Z1 = Z1 + VZ1
X2 = X2 + VX2
Y2 = Y2 + VY2
Z2 = Z2 + VZ2
TIIME = TIIME + 1
'';FOR III = 0 TO 100: NEXT ''; use this if your program is running too fast
scr_update
NEXT K
end
'' ;update screen display
sub scr_update
''; calculate center of mass of whole system
XC = (M1 * X1 + M2 * X2) / (M1 + M2)
YC = (M1 * Y1 + M2 * Y2) / (M1 + M2)
ZC = (M1 * Z1 + M2 * Z2) / (M1 + M2)
''; calculate position of galactic centers and display on screem
XX1 = (X1  XC)
YY1 = (Y1  YC)
ZZ1 = (Z1  ZC)
XX2 = (X2  XC)
YY2 = (Y2  YC)
ZZ2 = (Z2  ZC)
''; make black backgrounds
rect 1,1,xmax,ymax/25,0 filled
rect 1,ymax/2+5,xmax,ymax,0 filled
''; galaxy centers on XY view
CIRCLE (XX1*gsize) + (xmax/2), (YY1*gsize) + (ymax*0.75), 2,1, 9 '';11
CIRCLE (XX2*gsize) + (xmax/2), (YY2*gsize) + (ymax*0.75), 2,1, 13 '';15
''; galaxy centers on XYZ view
CIRCLE (XX1*gsize) + (xmax/2), (ZZ1*gsize) + (ymax*0.25), 2,1,9 '';11
CIRCLE (XX2*gsize) + (xmax/2),(ZZ2*gsize)+(ymax*0.25), 2,1, 13 '';15
''; place stars on screen
J=1
FOR I = 1 TO NS(J)
XP = (Xs(J,I)  XC): YP = (Ys(J,I)  YC): ZP = (Zs(J,I)  ZC)
''; stars of intruder galaxy on XY view
PSET (XP*gsize) + (xmax/2), (YP*gsize) + (ymax*0.75), 11
''; stars of intruder galaxy on XZ view
PSET (XP*gsize) + (xmax/2),(ZP*gsize) + (ymax*0.25), 11
NEXT I
J=2
FOR I = 1 TO NS(J)
XP = (Xs(J,I)  XC): YP = (Ys(J,I)  YC): ZP = (Zs(J,I)  ZC)
''; stars of target galaxy on XY view
PSET (XP*gsize) + (xmax/2), (YP*gsize) + (ymax*0.75), 15
''; stars of target galaxy on XZ view
PSET (XP*gsize) + (xmax/2),(ZP*gsize) + (ymax*0.25), 15
NEXT I
''; print out display headings
LOCATE 1,1: PRINT "Upper View = XZ ....Time > ", TIIME ,"Million years"
LOCATE 2,1: PRINT "Lower View = XY"
end sub
You will have to slow it down at least, made for older, slower versions of SB.

Really great! Thank for sharing. :)
Btw, I have something on my tablet that I will share later. All done on flights back in the good ol' days lol

There are at least 2 one is g_col for palm but that might be corrupted, it did not work for me.
<snip>
You will have to slow it down at least, made for older, slower versions of SB.
Thanks bplus. I tried Keijo Koskinen's one you posted as code above, which doesn't quite work for me out of the box  it's drawing on the line that I assume is meant to be blank in between the two views  does it work ok on yours? It might be a screen size issue or a case of tweaking an IF statement or two. Anyway, yes that looks like a proper job, judging by the sines and cosines in there, and I spotted a ^0.5 to complete the Pythagoras. It also deals with proper physical space with 3 dimensions. I'm amazed it still draws so fast with all that calculation, although it looks like the page updating is left to the OS to do when it's got time  there is no SHOWPAGE in there at all, and it's pretty jerky on my system...
Ah, I'm starting to get it. I had to change k to something a bit smaller to see most of the interaction (although it's still drawing outlying starts over the middle strip), and I stuck a showpage down at the bottom of the sub after the text printing. Now I can try some of those different speed factors remmed out near the top.... Nice. BTW, that is on the page of examples (http://"https://smallbasic.github.io/pages/samples.html") at github too, g_col_09.bas
My program is only 2D, so maybe it's a question for a maths/physics genius what the correct 'gravitational' equation would be, since 2D reality doesn't exist, or at least presumably doesn't have mass! I'm more of a bodger. :))

Really great! Thank for sharing. :)
Btw, I have something on my tablet that I will share later. All done on flights back in the good ol' days lol
Thanks Chris, I look forward to seeing those.

Galaxies Collide might be trying to show two different views, from above and profile, I know there was a version of this I've seen in past with 2 views and yes! this is running way too fast, from what I remember it ran agonizingly slow.
I think LS, yours more interesting!

You're too kind!
Yes, that's what it's "trying" to do  on my PC the two views are spilling over, and I have no idea from the code why they wouldn't when it was written. I've put some IF statements checking the y values first and fixed that, and also bumped up the number of millions of years it runs  you can still CtrlB out of it, of course if it's going on too long. One or two other tweaks. Here's my slightly modified version:
REM SmallBASIC
REM adapted by lettersquash: 28/03/2020
'
'';Galaxy Collision Program ver 0.9
'';Keijo Koskinen 20/4/2005 for smallbasic
'';Displays collision, views from XY and XZ direction
'';core algorithm modified from old Galaxy Collision
'';Astronomy magazine Dec 88 page 92.
'';by M.C.Schroeder and Neil F. Comins,
''
DIM Xs(2,1000), Ys(2,1000), Zs(2,1000) ''; target gal=(1,) intruder=(2,) coordinates
dim VXs(2,1000), VYs(2,1000), VZs(2,1000) ''; speed
dim NRR(2),NRS(2),NS(2),DR(2)
CLS
color 10,0
''; TARGET GALAXY
NRR(2) = 23 ''; orbits in target gal
NRS(2) = 10 ''; stars / orbit in target galaxy
NS(2) = NRR(2) * NRS(2) ''; Total of stars in target galaxy
DR(2) = 20 / (NRR(2)  1)
M2 = 6 ''; Target galaxy mass
X2 = 200 ''; initial Target galaxy coordinates
Y2 = 90
Z2 = 10
''; *** Try here different speed factor to greate different action
''; VX2=0.085 ; distroys totally galaxies
''; VY2=0
''; VX2=0.085 ''; really hot relationship in the beginning
''; VY2=0.065
' VX2= .085 '; galazies keep fighting over 1000 million years
' VY2=0.29
' VX2= 0.4 ''; galazies keep dancing over and over
' VY2=0.4
VX2=0.116 ''; What a tension
VY2=0.41
'VX2=0.116
'VY2=0.52
'';
VZ2 = 0
''; INTRUDER GALAXY
NRR(1) = 23 ''; rings in intruder galaxy
NRS(1) = 10 ''; stars/ring in intruder galaxy
NS(1) = NRR(1) * NRS(1) ''; Total of stars in intruder galaxy
DR(1) = 20 / (NRR(1)  1)
M1 = 6 ''; Intruder galaxy mass
X1 = 150
Y1 = 100 ''; Initial coordinates of intruder galaxy
Z1 = 0
VX1 = 0 ''; Intruder galaxy initial speed
VY1 = 0
VZ1 = 0
''; common factors
SF2 = 2 ''; softening factor
k=0.05 ''; size factor
gsize=k*xmax/((NRR(1)+NRR(2))/2)
NTSPR = 12000 ''; Time / millions of years
TIIME = 0
''; place stars in orbits in target galaxy
''; from r=4 and p at intervals of dr
I=1
FOR IR = 1 TO NRR(1)
R = 4 + (IR/1  1) * DR(1)
V = SQR(((M1*DR(1))/1) / R)
TH = (.5 * V / R) * (180 / 3.14159)
IF R <10 THEN V = .99 * V
if R<7 then V=.98*V
if R <5 then V=0.95*V
FOR IT = 1 TO NRS(1)
T = rnd * 360
I = I + 1
''; initialize stars positions
Xs(1,I) = R * COS(T / 57.2958) + X1
Ys(1,I) = R * SIN(T / 57.2958) + Y1
Zs(1,I) = Xz1+0
VZs(1,I) = VZ1+0
''; init star velocities in orbits
VXs(1,I) = V * SIN(T/57.2958) +VX1
VYs(1,I) = V * COS(T/57.2959) +VY1
NEXT IT
NEXT IR
''; place stars in orbits in intruder galaxy
''; from r=4 and up at intervals of dr
I=1
FOR IR = 1 TO NRR(2)
R = 4 + (IR/1  1) * DR(2)
V = SQR(((M2*DR(2))/1) / R)
TH = (.5 * V / R) * (180 / 3.14159)
IF R <9 THEN V = .93 * V
if R<7 then V=.92*V
if R <5 then V=0.9*V
FOR IT = 1 TO NRS(2)
T=rnd*360
I = I + 1
''; initialize stars positions
Xs(2,I) = R * COS(T / 57.2958) + X2
Ys(2,I) = R * SIN(T / 57.2958) + Y2
Zs(2,I) = Z2+0
VZs(2,I) = VZ2+0
'';init star velocities in orbits
VXs(2,I) = V * SIN(T/57.2958)+VX2
VYs(2,I) = V * COS(T/57.2958) +VY2
NEXT IT
NEXT IR
scr_update
''; particle pusher routine
label timesteps
FOR K = 1 TO NTSPR ''; ok lets go
''; determine the forces pulling stars in intruder gal
J = 2
FOR I = 1 TO NS(J)
R1 = M1 / ((Xs(J,I)  X1) ^ 2 + (Ys(J,I)  Y1) ^ 2 + (Zs(J,I)  Z1) ^ 2 + SF2) ^ 1.5
R2 = M2 / ((Xs(J,I)  X2) ^ 2 + (Ys(J,I)  Y2) ^ 2 + (Zs(J,I)  Z2) ^ 2 + SF2) ^ 1.5
''; calculate stars x,y,z accelerations
AX = R1 * (X1  Xs(J,I)) + R2 * (X2  Xs(J,I))
AY = R1 * (Y1  Ys(J,I)) + R2 * (Y2  Ys(J,I))
AZ = R1 * (Z1  Zs(J,I)) + R2 * (Z2  Zs(J,I))
''; update star positions and velocities
VXs(J,I) = VXs(J,I) + AX
VYs(J,I) = VYs(J,I) + AY
VZs(J,I) = VZs(J,I) + AZ
Xs(J,I) = Xs(J,I) + VXs(J,I)
Ys(J,I) = Ys(J,I) + VYs(J,I)
Zs(J,I) = Zs(J,I) + VZs(J,I)
NEXT I
''; determine the forces pulling stars in target gal
J = 1
FOR I = 1 TO NS(J)
R1 = M1 / ((Xs(J,I)  X1) ^ 2 + (Ys(J,I)  Y1) ^ 2 + (Zs(J,I)  Z1) ^ 2 + SF2) ^ 1.5
R2 = M2 / ((Xs(J,I)  X2) ^ 2 + (Ys(J,I)  Y2) ^ 2 + (Zs(J,I)  Z2) ^ 2 + SF2) ^ 1.5
''; calculate stars x,y,z accelerations
AX = R1 * (X1  Xs(J,I)) + R2 * (X2  Xs(J,I)) ''kokeilu nurin x1 ja x2
AY = R1 * (Y1  Ys(J,I)) + R2 * (Y2  Ys(J,I)) '' samoin
AZ = R1 * (Z1  Zs(J,I)) + R2 * (Z2  Zs(J,I)) '' kuin myos
''; update star positions and velocities
VXs(J,I) = VXs(J,I) + AX
VYs(J,I) = VYs(J,I) + AY
VZs(J,I) = VZs(J,I) + AZ
Xs(J,I) = Xs(J,I) + VXs(J,I)
Ys(J,I) = Ys(J,I) + VYs(J,I)
Zs(J,I) = Zs(J,I) + VZs(J,I)
NEXT I
''; update postions and velocities of the galactic centers
RR = ((X1  X2) ^ 2 + (Y1  Y2) ^ 2 + (Z1  Z2) ^ 2 + SF2) ^ 1.5
AX = (X2  X1) / RR: AY = (Y2  Y1) / RR: AZ = (Z2  Z1) / RR
VX1 = VX1 + M2 * AX: VY1 = VY1 + M2 * AY: VZ1 = VZ1 + M1 * AZ ''; +M1*AZ
VX2 = VX2  M1 * AX: VY2 = VY2  M1 * AY: VZ2 = VZ2  M1 * AZ
X1 = X1 + VX1
Y1 = Y1 + VY1
Z1 = Z1 + VZ1
X2 = X2 + VX2
Y2 = Y2 + VY2
Z2 = Z2 + VZ2
TIIME = TIIME + 1
' FOR III = 0 TO 100: NEXT ''; use this if your program is running too fast
scr_update
NEXT K
end
'' ;update screen display
sub scr_update
''; calculate center of mass of whole system
XC = (M1 * X1 + M2 * X2) / (M1 + M2)
YC = (M1 * Y1 + M2 * Y2) / (M1 + M2)
ZC = (M1 * Z1 + M2 * Z2) / (M1 + M2)
''; calculate position of galactic centers and display on screem
XX1 = (X1  XC)
YY1 = (Y1  YC)
ZZ1 = (Z1  ZC)
XX2 = (X2  XC)
YY2 = (Y2  YC)
ZZ2 = (Z2  ZC)
''; make black backgrounds
rect 0,0,xmax+1,ymax+1,0 filled
' rect 1,ymax/2+5,xmax,ymax,0 filled
''; galaxy centers on XY view
if (YY1*gsize) + (ymax*0.75)>ymax/2+7 then CIRCLE (XX1*gsize) + (xmax/2), (YY1*gsize) + (ymax*0.75), 2,1, 9 '';11
if (YY2*gsize) + (ymax*0.75)>ymax/2+7 then CIRCLE (XX2*gsize) + (xmax/2), (YY2*gsize) + (ymax*0.75), 2,1, 13 '';15
''; galaxy centers on XZ view
CIRCLE (XX1*gsize) + (xmax/2), (ZZ1*gsize) + (ymax*0.25), 2,1,9 '';11
CIRCLE (XX2*gsize) + (xmax/2),(ZZ2*gsize)+(ymax*0.25), 2,1, 13 '';15
''; place stars on screen
J=1
FOR I = 1 TO NS(J)
XP = (Xs(J,I)  XC): YP = (Ys(J,I)  YC): ZP = (Zs(J,I)  ZC)
''; stars of intruder galaxy on XY view
if (YP*gsize) + (ymax*0.75)>ymax/2+5 then PSET (XP*gsize) + (xmax/2), (YP*gsize) + (ymax*0.75), 14
''; stars of intruder galaxy on XZ view
if (ZP*gsize) + (ymax*0.25) < ymax/25 then PSET (XP*gsize) + (xmax/2),(ZP*gsize) + (ymax*0.25), 14
NEXT I
J=2
FOR I = 1 TO NS(J)
XP = (Xs(J,I)  XC): YP = (Ys(J,I)  YC): ZP = (Zs(J,I)  ZC)
''; stars of target galaxy on XY view
if (YP*gsize) + (ymax*0.75)>ymax/2+5 then PSET (XP*gsize) + (xmax/2), (YP*gsize) + (ymax*0.75), 15
''; stars of target galaxy on XZ view
if (ZP*gsize) + (ymax*0.25) < ymax/25 then PSET (XP*gsize) + (xmax/2),(ZP*gsize) + (ymax*0.25), 15
NEXT I
''; print out display headings
LOCATE 1,1: PRINT "Upper View = XZ ....Time > ", TIIME ,"Million years"
LOCATE 2,1: PRINT "Lower View = XY"
showpage
end sub

Oh wow! what an improvement! Nice detective work.
And it runs fast now so not agonizing at all, beautiful.

I just spotted a mistake in mine. Line 50, I do the circles in COLOR i%16, forgetting that this will include 0s  just in case you're puzzled  there'll be a few black holes in there! It might be nicer with colours related to something meaningful like the mass. Say color m(i)+9.
For cycling all the standard colours 1 to 15, as I intended, I think that should be color (i%15)+1.

I managed to write that collapsing galaxy/solarsystem thing. I also fixed the annoying tendency of the system to float off somewhere. Let me know if you find any bugs, or if you tweak the formulae and it gives great results. Of course, it's a chaotic system that starts with loads of random numbers, so it's always gonna do something different.
REM SmallBASIC
REM created: 23/03/2020
'Gravity2. Simulation of star system formation, by lettersquash.
'Have fun with the numbers. Most of the relevant formulae are commented.
'G below will create different effects, but all the variables will affect the others.
'Body with largest mass (in yellow) changes, all are drawn towards centre of screen with it.
'Can override somewhat with clickhold LMB on screen. +/ buttons zoom (work in progress).
'Continues indefinitely  CTRL+B to stop.
option predef antialias off
n=100 'number of initial objects
G=0.05 'gravitational constant (in this universe), bigger makes a fast pace & can handle larger n.
dim x (n), y (n), dx (n), dy (n), m (n) 'position, speed, mass
' BORDER, SPACE FOR OUTPUT AND BUTTONS (intend more controls soon)
border=56
mix=(xmaxborder)\2:miy=ymax\2
'rect 0,0,border,ymax,3 filled
'zoom in button
rect 4,100,step 48,48,3 filled
rect 12,122,step 32,4,7 filled
rect 26,108,step 4,32,7 filled
'zoom out button
rect 4,150,step 48,48,3 filled
rect 12,172,step 32,4,7 filled
locate 0,0:?" n="
'SETUP
for i=1 to n
' initial positions
x(i)=rnd*(xmaxborder)+border
y(i)=rnd*ymax
'initial velocity (rotation by this formula or rem out)
dx(i)=(y(i)miy)*0.003
dy(i)=(mixx(i))*0.003
'initial masses (bodies drawn with radius the root of this value, but they combine & grow)
m(i)=rnd*3+0.75
next
'MAIN ROUTINE
repeat
gravitate
showpage
delay 1
' if inkey=" " then delay 100:pause
until 0
END '===============
sub gravitate
'Clear the drawing area
rect border,0,xmax+1,ymax+1,0 filled
for i=1 to n
for j=1 to n
'Gravitational forces caused by each body are summed as x,y components
if i<>j then
'x,y components of distance between the two
a=x(j)x(i):b=y(j)y(i)
'and distance squared
d2=(a^2+b^2)
'product of masses
m2=m(i)*m(j)
'############ HOW EASILY THEY COMBINE ###########
'Close enough and slow enough for gravitation to combine masses?...
if d2 < m(i)^0.5+m(j)^0.5 and (dx(i)dx(j)) < m2 and (dy(i)dy(j)) < m2 then combine:exit sub
'approximation of gravitational equation for dx,dy
'Force of gravity = G*m1*m2/d^2, but m1*m2 seems too much, just using m(j)
dx(i)+=m(j)*G*a/d2
dy(i)+=m(j)*G*b/d2
fi
next j
'update position. Inertia is introduced here, /m(i)
x(i)+=dx(i)/m(i):y(i)+=dy(i)/m(i)
'a bit of 'drag' gives quicker collapse, needs to be only just < 1...
'0.99 gives fairly quick collapse into one body, or rem out for the long view.
' dx(i)*=0.9999:dy(i)*=0.9999
if pen(3) then click ' allow pan/zoom via LMB
c=7 ' colour of lower masses
if m(i)=max(m) then
c=14 'colour of biggest mass
'the following provides stability, bringing biggest mass towards centre
px=x(i)mix:py=y(i)miy
for k=1 to n
x(k)=px/(n*10):y(k)=py/(n*10)
next
fi
if x(i)>border+20 then circle x(i),y(i),m(i)^0.5,1,c filled 'draw bodies
next i
end sub
sub combine
'give index i combined values from i and j
dx(i)+=dx(j)/m(i)
dy(i)+=dy(j)/m(i)
m(i)+=m(j)
'undraw j object and delete from arrays
if x(j)>border+20 then circle x(j),y(j),m(j)^0.5,1,0 filled
delete x,j:delete y,j:delete dx,j:delete dy,j:delete m,j
n 'now one less object. Exit sub from main loop avoids out of bounds error
locate 1,0:? usg " ####";n 'print remaining number of bodies
end sub
sub click
p4=pen(4):p5=pen(5)
if p4<border4 and p5>100 and p5<148 then
'zoom in with "+" button
x(i)=mix(mixx(i))*1.05
y(i)=miy(miyy(i))*1.05
elseif p4<border4 and p5>150 and p5<198 then
'zoom out with "" button
x(i)=mix(mixx(i))*0.95
y(i)=miy(miyy(i))*0.95
else
'"Pan"  move screen proportional to mouse distance from centre with LMB down
'Would be good with normal drag to move, but couldn't figure it out.
px=p4mix:py=p5miy
x(i)=(px/100):y(i)=(py/100)
fi
end sub

Hey LS,
Check out the gravity in this little Galaxy called Aurel's toy :))
a=10.1 : sc=175 : ox=400 : oy=240
while a < 1
cls
for x= 0.1 to 0.8 step 0.05
for y=0.1 to 0.8 step 0.05
lx=x : ly=y
for n=1 to 40
xx=lx*cos(a)(lylx*lx)*sin(a)
ly=lx*sin(a)+(lylx*lx)*cos(a)
lx=xx
if abs(lx+ly)>3000 then n=41 'end if
pset ly*sc+ox,lx*sc+oy,rgb(n * 5, 255  n * 5, 128 + n * 2)
next
next
next
showpage
a=a+0.005
wend

That's not gravity. That's witchcraft.