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

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

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

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.

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

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

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

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