Gravity - a fun demo

[lettersquash]

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.

Code: [Select]

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=xmax-50
mix=xmax\2:miy=ymax\2
rect border,0,xmax,ymax,3 filled

rect border+2,miy-150,xmax-2,miy-100,2 filled
rect border+12,miy-126,xmax-12,miy-124,7 filled
rect border+24,miy-138,border+26,miy-112,7 filled

rect border+2,miy-98,xmax-2,miy-48,2 filled
rect border+12,miy-74,xmax-12,miy-72,7 filled

for i=1 to n
  ' initial positions
  x(i)=rnd*400+mix-200
  y(i)=rnd*400+miy-200
  'initial velocity (perhaps rotation)
'  dx(i)=(y(i)-miy)*0.005
'  dy(i)=(mix-x(i))*0.005
  m(i)=rnd*4+0.1
next

repeat
  for i=1 to n
    'Blank out body
    if x(i)<border-3 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)<border-3 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<xmax-2 and p5>miy-150 and p5<miy-100 then
   'zoom in with "+" button
    x(i)=mix-(mix-x(i))*1.05
    y(i)=miy-(miy-y(i))*1.05
  elseif p4>border+2 and p4<xmax-2 and p5>miy-98 and p5<miy-48 then
    'zoom out with "-" button
    x(i)=mix-(mix-x(i))*0.95
    y(i)=miy-(miy-y(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=p4-mix:py=p5-miy
    x(i)-=(px/100):y(i)-=(py/100)
  fi
end sub


[bplus]

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)?

[lettersquash]

Yeah the mouse thing is a bit weird at first and might be better with a more normal click-drag-release 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 double-click-hold 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 (low-mass) 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 high-speed 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!

[bplus]

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 ;-)

[lettersquash]

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?

[bplus]

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

Code: [Select]

'
'';Galaxy Collision Program  ver 0.9
'';Keijo Koskinen 20/4/2005  for smallbasic
'';Displays collision, views from X-Y and X-Z 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/2-5,0 filled
 rect 1,ymax/2+5,xmax,ymax,0 filled
 
 ''; galaxy centers on X-Y 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 X-YZ 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 X-Y view
  PSET (XP*gsize) + (xmax/2), (YP*gsize) + (ymax*0.75), 11
 ''; stars of intruder galaxy on X-Z 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 X-Y view
  PSET (XP*gsize) + (xmax/2), (YP*gsize) + (ymax*0.75), 15
 ''; stars of target galaxy on X-Z view
   PSET (XP*gsize) + (xmax/2),(ZP*gsize) + (ymax*0.25), 15
NEXT I

 ''; print out display headings
LOCATE 1,1: PRINT "Upper View = X-Z                             ....Time > ", TIIME ,"Million years"
LOCATE 2,1: PRINT "Lower View = X-Y"

 end sub


You will have to slow it down at least, made for older, slower versions of SB.

[chrisws]

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

[lettersquash]

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 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.

[lettersquash]

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.

[bplus]

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!

[lettersquash]

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 Ctrl-B out of it, of course if it's going on too long. One or two other tweaks. Here's my slightly modified version:

Code: [Select]

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 X-Y and X-Z 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 X-Y 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 X-Z 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 X-Y 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 X-Z view
  if (ZP*gsize) + (ymax*0.25) < ymax/2-5 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 X-Y 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 X-Z view
    if (ZP*gsize) + (ymax*0.25) < ymax/2-5 then PSET (XP*gsize) + (xmax/2),(ZP*gsize) + (ymax*0.25), 15
NEXT I

 ''; print out display headings
LOCATE 1,1: PRINT "Upper View = X-Z                             ....Time > ", TIIME ,"Million years"
LOCATE 2,1: PRINT "Lower View = X-Y"
showpage
 end sub



[bplus]

Oh wow! what an improvement! Nice detective work.

And it runs fast now so not agonizing at all, beautiful.

[lettersquash]

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.

[lettersquash]

I managed to write that collapsing galaxy/solar-system 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.

Code: [Select]

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 click-hold 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=(xmax-border)\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*(xmax-border)+border
  y(i)=rnd*ymax
  'initial velocity (rotation by this formula or rem out)
  dx(i)=(y(i)-miy)*0.003
  dy(i)=(mix-x(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<border-4 and p5>100 and p5<148 then
    'zoom in with "+" button
    x(i)=mix-(mix-x(i))*1.05
    y(i)=miy-(miy-y(i))*1.05
  elseif p4<border-4 and p5>150 and p5<198 then
    'zoom out with "-" button
    x(i)=mix-(mix-x(i))*0.95
    y(i)=miy-(miy-y(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=p4-mix:py=p5-miy
    x(i)-=(px/100):y(i)-=(py/100)
  fi
end sub


[bplus]

Hey LS,

Check out the gravity in this little Galaxy called Aurel's toy 

Code: [Select]

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)-(ly-lx*lx)*sin(a)
ly=lx*sin(a)+(ly-lx*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