 May 25, 2020, 08:59:05 AM

### Author Topic: Gravity - a fun demo  (Read 573 times)

#### lettersquash ##### Gravity - a fun demo
« 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.
Code: [Select]
`REM SmallBASICREM 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 offn=50dim x (n), y (n), dx (n), dy (n), m (n)border=xmax-50mix=xmax\2:miy=ymax\2rect border,0,xmax,ymax,3 filledrect border+2,miy-150,xmax-2,miy-100,2 filledrect border+12,miy-126,xmax-12,miy-124,7 filledrect border+24,miy-138,border+26,miy-112,7 filledrect border+2,miy-98,xmax-2,miy-48,2 filledrect border+12,miy-74,xmax-12,miy-72,7 filledfor 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.1nextrepeat  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:pauseuntil 0sub 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)  fiend sub`

#### bplus ##### Re: Gravity - a fun demo
« Reply #1 on: March 27, 2020, 01:13:38 AM »
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 ##### Re: Gravity - a fun demo
« Reply #2 on: March 27, 2020, 11:22:36 AM »
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.

#### bplus ##### Re: Gravity - a fun demo
« Reply #3 on: March 27, 2020, 03:02:21 PM »
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 ##### Re: Gravity - a fun demo
« Reply #4 on: March 27, 2020, 05:30:32 PM »
Wow LS!

Thanks for taking the time for "lecture". You have me thinking... Then my work here is done. Quote
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 ##### Re: Gravity - a fun demo
« Reply #5 on: March 27, 2020, 08:44:51 PM »
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,) coordinatesdim  VXs(2,1000), VYs(2,1000), VZs(2,1000)  ''; speeddim  NRR(2),NRS(2),NS(2),DR(2)CLScolor 10,0''; TARGET GALAXYNRR(2) = 23                ''; orbits in target galNRS(2) = 10                ''; stars / orbit in target galaxyNS(2) = NRR(2) * NRS(2)    ''; Total of stars in target galaxyDR(2) = 20 / (NRR(2) - 1)M2 = 6                     ''; Target galaxy massX2 = 200                   ''; initial Target galaxy coordinatesY2 = 90Z2 = 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.41VX2=-0.116VY2=0.52                  '';VZ2 = 0''; INTRUDER GALAXYNRR(1) = 23                 ''; rings in intruder galaxyNRS(1) = 10                 ''; stars/ring in  intruder galaxyNS(1) = NRR(1) * NRS(1)     ''; Total of stars in intruder galaxyDR(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 drI=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 drI=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 routinelabel 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 displaysub 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=1FOR 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), 11NEXT I J=2FOR 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), 15NEXT I ''; print out display headingsLOCATE 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

• Jr. Member
•  • Posts: 43 ##### Re: Gravity - a fun demo
« Reply #6 on: March 28, 2020, 12:12:46 AM »
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 ##### Re: Gravity - a fun demo
« Reply #7 on: March 28, 2020, 11:26:18 AM »
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 ##### Re: Gravity - a fun demo
« Reply #8 on: March 28, 2020, 11:28:52 AM »
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 ##### Re: Gravity - a fun demo
« Reply #9 on: March 28, 2020, 04:17:20 PM »
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 ##### Re: Gravity - a fun demo
« Reply #10 on: March 28, 2020, 06:35:03 PM »
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 SmallBASICREM 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,) coordinatesdim  VXs(2,1000), VYs(2,1000), VZs(2,1000)  ''; speeddim  NRR(2),NRS(2),NS(2),DR(2)CLScolor 10,0''; TARGET GALAXYNRR(2) = 23                ''; orbits in target galNRS(2) = 10                ''; stars / orbit in target galaxyNS(2) = NRR(2) * NRS(2)    ''; Total of stars in target galaxyDR(2) = 20 / (NRR(2) - 1)M2 = 6                     ''; Target galaxy massX2 = 200                   ''; initial Target galaxy coordinatesY2 = 90Z2 = 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 GALAXYNRR(1) = 23                 ''; rings in intruder galaxyNRS(1) = 10                 ''; stars/ring in  intruder galaxyNS(1) = NRR(1) * NRS(1)     ''; Total of stars in intruder galaxyDR(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 drI=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 drI=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 routinelabel 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 displaysub 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=1FOR 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), 14NEXT I J=2FOR 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), 15NEXT I ''; print out display headingsLOCATE 1,1: PRINT "Upper View = X-Z                             ....Time > ", TIIME ,"Million years"LOCATE 2,1: PRINT "Lower View = X-Y"showpage end sub`

#### bplus ##### Re: Gravity - a fun demo
« Reply #11 on: March 28, 2020, 07:24:51 PM »
Oh wow! what an improvement! Nice detective work.

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

#### lettersquash ##### Re: Gravity - a fun demo
« Reply #12 on: March 28, 2020, 11:03:39 PM »
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 ##### Re: Gravity - Part 2, The Collapse!
« Reply #13 on: March 30, 2020, 01:31:34 AM »
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 SmallBASICREM 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 offn=100 'number of initial objectsG=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=56mix=(xmax-border)\2:miy=ymax\2'rect 0,0,border,ymax,3 filled'zoom in buttonrect 4,100,step 48,48,3 filledrect 12,122,step 32,4,7 filledrect 26,108,step 4,32,7 filled'zoom out buttonrect 4,150,step 48,48,3 filledrect 12,172,step 32,4,7 filledlocate 0,0:?" n="'SETUPfor 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.75next'MAIN ROUTINErepeat  gravitate  showpage  delay 1 ' if inkey=" " then delay 100:pauseuntil 0END '===============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 iend subsub 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 bodiesend subsub 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)  fiend sub`

#### bplus ##### Re: Gravity - a fun demo
« Reply #14 on: April 08, 2020, 12:28:25 AM »
Hey LS,

Check out the gravity in this little Galaxy called Aurel's toy Code: [Select]
`a=-10.1 : sc=175 : ox=400 : oy=240while a < 1cls 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 nextnextshowpage a=a+0.005wend`