falling sand

[tsh73]

Inspired by this RC task
rosettacode.org/wiki/Abelian_sandpile_model
Thoght really I can't made sence out of it :(
Apparently numbers there is NOT sand level, because it is set never be more then 4. So wht is it then?

So I just drop a sandpiece at 0,0
and see if it towers over one of four neighbours over one unit height.
If it does it falls that way (if several ways, it falls randomly)
So pile is steadily grows.
I visualise it as it go (only top level particles, not whole pile undeneath)
in isometric or dimetric projection (two lines in inmost loop, uncomment any two of them).

The problem is, it looks pretty rational, stright and pretty boring.
(and I tried turn off that "it falls randomly" bit - no visible difference)

Any suggestions?


'todo a sandpile
'isometric
'then sand makes 2-pile it falls to random empty side
   'probably start with a top view
   '01 - setup colors
   '02 - make square field (btw, may be do not show level 0?)
   '03 - switch to isometry or dimetry ( dimetry is easier, start with it)
   '04 - add sand drops

'for i = -3 to 3
'print i, ifPos(i)
'next
'end

global w2, w
w2=11
w=2*w2
maxCol=10
dim e(w,w)  'sand level
dim dx(4):dx(1)=1:dx(2)=-1:dx(3)=0:dx(4)=0
dim dy(4):dy(1)=0:dy(2)= 0:dy(3)=1:dy(4)=-1


'nomainwin
open "test" for graphics_nsb_nf as #gr
#gr "trapclose [quit]"
#gr "down"
#gr "fill white; flush"
#gr "home; posxy cx cy"

pi=acs(-1)
sq32=sqr(3)/2
sz=7
#gr "size ";sz
#gr "size ";sz-2

goto [skeepDrawSpeedTest]
goto [skipColors]
#gr "size 5"
h=0.02
'for x = 0 to 1 step h    'red to red
'for x = 0 to 2/3 step h     'RED TO BLUE
for x = 0 to 1 step h    'now blue to red
   '0->2/3, 1->0
   c=(1-x)*2/3
   #gr "color ";rainbow$(c)
   dy=0: if x mod (1/6) <h then dy = 5
   #gr "set ";x*100*2+50;" ";cy+dy
next

[skipColors]
[fill]
s=0
for i = 0-w2 to w2
   for j= 0-w2 to w2
       'e(i+w2,j+w2)=int(maxCol*rnd(0))+1
       'e(i+w2,j+w2)=int((1-(abs(i)+abs(j))/w)*maxCol)
       'e(i+w2,j+w2)=(abs(i)+abs(j))
       'e(i+w2,j+w2)=int(cos((abs(i)+abs(j))/w/2*pi)*maxCol)
       'e(i+w2,j+w2)=int(sqr(abs(1-(i/w2)^2-(j/w2)^2))*maxCol)
       'e(i+w2,j+w2)=int(sqr(ifPos(1-(i/w2)^2-(j/w2)^2))*maxCol)
       e(i+w2,j+w2)=int((ifPos(1-(i/w2)^2-(j/w2)^2))*maxCol)
       s=s+e(i+w2,j+w2)
       print using("##",e(i+w2,j+w2));
   next
   print
next
print "Pt num: ";s

t0=time$("ms")
N=10
for ii = 1 to N
   gosub [draw]
next
t1=time$("ms")
print "frames per sec", N/(t1-t0)*1000
print "ms per frame", (t1-t0)/N

[skeepDrawSpeedTest]
t0=time$("ms")
N=1000
while 1
   frame=frame+1
   'i=int(rnd(0)*w)
   'j=int(rnd(0)*w)
   'e(i,j)=e(i,j)+1
   call add w2,w2
   gosub[draw]
   if frame>=N then exit while
wend
t1=time$("ms")

s=0
for i = 0-w2 to w2
   for j= 0-w2 to w2
       print using("##",e(i+w2,j+w2));
   next
   print
next

print "frames per sec", N/(t1-t0)*1000
print "ms per frame", (t1-t0)/N

#gr "flush"
wait

sub add i,j 'recursive  'w w2 global
   for d = 1 to 4'direction
       'valid?
       if i+dx(d)>=0 and i+dx(d)<=w and j+dy(d)>=0 and j+dy(d)<=w then
           'diff>1?
           if e(i,j)-e(i+dx(d),j+dy(d))>0 then dirs$=dirs$;d
       end if
   next

   if len(dirs$)=0 then e(i,j)=e(i,j)+1: exit sub 'no way to spill
   randN=int(rnd(0)*len(dirs$)+1)
   'randN=1    'always first
   d=val(mid$(dirs$,randN,1))
   call add i+dx(d),j+dy(d)    'recursive
end sub

[draw]
#gr "cls"
#gr "place 10 30"
#gr "\";frame
for i = 0-w2 to w2
   SCAN
   for j= 0-w2 to w2
       c=(1-e(i+w2,j+w2)/maxCol)*2/3
       if c<>lastC then   'extra frame per sec
           #gr "color ";rainbow$(c)
           lastC=c
       end if
       x=j
       y=i
       z=e(i+w2,j+w2)
       'print using("##",e(i+w2,j+w2));
       'dimetry
       xx=cx+sz*(x+y/sqr(2))
       yy=cy+sz*(0-z+y/sqr(2))
       'izometry
       'xx=cx+sz*(x/sq32-y/sq32)
       'yy=cy+sz*(0-z+y/sq32/2+x/sq32/2)
       #gr "set ";int(xx);" ";int(yy)
   next
   'print
next
return

[quit]
timer 0
close #gr
end

function ifPos(x)
   ifPos=x*(x>0)
end function
'---------------------------------------------
' 0..1 into red-green-blue-red continuous colors
function rainbow$(x)
   hi = int((x*6) mod 6)+ 5*(x<0)  'fixed to 0..5
   f = (x*6) mod 1 + (x<0)         'frac, 0..1
   q = (1-f)
   select case hi
   case 0
       r = 1: g = f: b = 0
   case 1
       r = q: g = 1: b = 0
   case 2
       r = 0: g = 1: b = f
   case 3
       r = 0: g = q: b = 1
   case 4
       r = f: g = 0: b = 1
   case 5
       r = 1: g = 0: b = q
   end select
   R = int(r*255)
   G = int(g*255)
   B = int(b*255)
   rainbow$= R;" ";G;" ";B
end function


[Rod]

After a very quick read I think every grain of sand needs reassessed every cycle and any grain of sand that now has four neighbours or perhaps four grains above will move or they will. So the whole pile moves not just grains piling on top with a static centre.

Need to read more.

[tenochtitlanuk]

It's a fun one! See also my posting on LB forum linked to sandpile

I like the 3D representation, Anatoly!

[tsh73]

a bit unfinished
(I really should calculate height of piramid from N and set maxCol,
may be some over stuff)

You can pause things - reset controls - hit Redraw

Redraw button does better recizing (of 2d/3d), but you can switch 2d/3d/etc on the go.
It really looks better with 2d :( and probably smaller half-size value.
Sand trail is a path of last sand particle - or marks where Abelian sand avalanches
Rundom jump refers to condition where particle could go several ways. It either goes at random one, or in order
(does not used for Abelian sand)
Top level only works in 3d view - shows only top of pile (vs. each particle from bottom. Significally slows down).

Things that do not work without starting anew are
* field size (half-width) - because center on screen is 0,0 but in array it is (half_width, half_width)
* N (it just does not get read after start, so as for now you cannot add another 300 sandparticles)
* switching Abelian on/off (it does technically work, it just does not make any sence)

'   Sandpile (Abelian too)
' tsh73 Nov 2021

   nomainwin

global w2, w
global cx, cy, sz
global isShowGrid, isTopLev, isRandJump
global isShowTrail, trail$
global isDimetry, is2d, isAbelian

   WindowWidth = 608
   WindowHeight = 625

   UpperLeftX=int((DisplayWidth-WindowWidth)/2)
   UpperLeftY=int((DisplayHeight-WindowHeight)/2)

   texteditor #main.log, 374, 361, 216, 185
   graphicbox #main.gr, 5, 5, 350, 350
   button #main.btnStart, "Start Anew", [StartClick], UL, 366, 11, 122, 25
   button #main.btnPause, "Pause/Cont", [PauseClick], UL, 366, 46, 122, 25
   statictext #main.statictext5, "N", 366, 86, 25, 20
   textbox #main.txtN, 406, 76, 64, 25
   statictext #main.statictext7, "curr", 366, 111, 40, 20
textbox #main.txtCurr, 406, 106, 64, 25
   radiobutton #main.rb2d, "2d", [2dSet], [dummy], 366, 216, 40, 20
   radiobutton #main.rb3d, "3d", [3dSet], [dummy], 366, 241, 40, 20
   groupbox #main.groupbox12, "3d projection", 414, 241, 152, 100
   radiobutton #main.rbDimetry, "Dimetric", [dimetrySet], [dummy], 430, 261, 80, 20
   radiobutton #main.rbIsomery, "Isometric", [isometrySet], [dummy], 430, 281, 88, 20
   checkbox #main.chkAbelian, "Abelian sand", [abelianSet], [abelianReset], 366, 136, 112, 20
   checkbox #main.chkTopLev, "Top level only ", [topLevSet], [topLevReset], 430, 306, 136, 20
   statictext #main.statictext18, "Log", 366, 336, 144, 20
   checkbox #main.chkTrail, "Sand trail", [trailSet], [trailReset], 366, 161, 96, 20
   checkbox #main.chkGrid, "Floor grid", [gridSet], [gridReset], 366, 186, 96, 20
   statictext #main.statictext21, "half size", 510, 81, 64, 20
   textbox #main.txtW2, 510, 106, 64, 25
   button #main.btnRedraw, "Redraw", [RedrawClick], UL, 494, 46, 88, 25
   checkbox #main.chkRand, "Random jump", [randSet], [randReset], 462, 161, 120, 20  
   menu #main, "Edit"    '<--- Texteditor Menu can be moved but not removed.

   open "Sandpile" for window as #main
   print #main.gr,"down; fill white; flush"
#main.gr "home; posxy cx cy"
   print #main, "trapclose [quit.main]"

   print #main, "font ms_sans_serif 10"
#main.txtCurr "!disable"
gosub [init]


   wait

[quit.main]
   Close #main
   END

'chkBoxes/radioButtons ==========

[2dSet]
is2d=1
   wait

[3dSet]
is2d=0
   wait

[dimetrySet]
isDimetry=1
   wait

[isometrySet]
isDimetry=0
   wait

[abelianSet]
isAbelian=1
   wait

[abelianReset]
isAbelian=0
   wait

[gridSet]
isShowGrid=1
   wait

[gridReset]    
isShowGrid=0
   wait

[trailSet]  
isShowTrail=1
   wait

[trailReset]
isShowTrail=0
   wait

[randSet]
isRandJump=1
   wait

[randReset]
isRandJump=0
   wait

[topLevSet]
isTopLev=1
   wait

[topLevReset]
isTopLev=0
   wait

[dummy]    'non-used radio reset handler
   wait

'buttons =====================

[RedrawClick]
gosub [initCalc]
gosub [doDraw]
   wait

[StartClick]
doWork=1
doPause=0
frame=0
#main.txtCurr frame
#main.log "started ";time$()
gosub [initCalc]
redim e(w,w) 'clear
goto [doFrame] 'will call itself with a timer
   wait

[PauseClick]
if doPause then
doPause=0
#main.log "continued ";time$()
goto [doFrame]
else
if doWork=1 then
doPause=1
#main.log "paused ";time$()
end if
end if
   wait

'returnable SUBs ======================
[init]
w2=11
#main.txtW2 w2
w=2*w2

N=300
#main.txtN N

is2d=1
#main.rb2d iif$(is2d,"set","reset")
#main.rb3d iif$(not(is2d),"set","reset")

isDimetry=1
#main.rbDimetry iif$(isDimetry,"set","reset")
#main.rbIsomery iif$(not(isDimetry),"set","reset")

isShowGrid=1
#main.chkGrid  iif$(isShowGrid,"set","reset")

isTopLev=1
#main.chkTopLev iif$(isTopLev,"set","reset")

isAbelian=0
#main.chkAbelian iif$(isAbelian,"set","reset")

isShowTrail=1
#main.chkTrail iif$(isShowTrail,"set","reset")

isRandJump=1
#main.chkRand iif$(isRandJump,"set","reset")

if isAbelian then
maxCol=3
else
maxCol=10 'should be counted from N
end if


dim e(w,w)  'sand level
'directions array
dim dx(4):dx(1)=1:dx(2)= 0:dx(3)=-1:dx(4)=0
dim dy(4):dy(1)=0:dy(2)= 1:dy(3)= 0:dy(4)=-1
return

[initCalc]
#main.txtN "!contents? N"

#main.txtW2 "!contents? w2"
w=2*w2

if is2d then 'fill the place
sz=int(2*cx/(w+2))
ptSize=sz-2
else
if isDimetry then 'trial and error factor
sz=int(2*cx/(w+2)/1.5)
ptSize=sz-4
else
sz=int(2*cx/(w+2)/2)
ptSize=sz
end if
end if
return

[doDraw]
#main.gr "cls"
if isShowGrid then
#main.gr "size 1; color black"
z=0
for i = 0-w2 to w2
' x=j
y=i
xx1=sx(0-w2,y,z)
yy1=sy(0-w2,y,z)
xx2=sx(w2,y,z)
yy2=sy(w2,y,z)
#main.gr "line ";xx1;" ";yy1;" ";xx2;" ";yy2
next

for j= 0-w2 to w2
x=j
'y=i
xx1=sx(x,0-w2,z)
yy1=sy(x,0-w2,z)
xx2=sx(x,w2,z)
yy2=sy(x,w2,z)
#main.gr "line ";xx1;" ";yy1;" ";xx2;" ";yy2
next
end if

#main.gr "size ";ptSize
for i = 0-w2 to w2
   'SCAN
for j= 0-w2 to w2
x=j
y=i
zz=e(i+w2,j+w2)
if zz then 'and we skip 0
if isTopLev or is2d then minZZ=zz else minZZ=1
for z = minZZ to zz
'print using("##",e(i+w2,j+w2));
c=(1-(z-1)/maxCol)*2/3
#main.gr "color ";rainbow$(c)
xx=sx(x,y,z)
yy=sy(x,y,z)
#main.gr "set ";xx;" ";yy
next
end if
next
'print
next

return

' timer sub ==========================
[doFrame]
timer 0 'just for a case
frame=frame+1
#main.txtCurr frame
if frame>=N then
doWork=0
#main.log "finished ";time$()
wait
end if
trail$="color red"
if isAbelian then
call AbelianAdd w2,w2
else
call Add w2,w2
end if
gosub [doDraw]
if isShowTrail then #main.gr trail$

if doWork and not(doPause) then timer 150, [doFrame] 'recursively calls itself
   wait

' Specific subs/functions ====================
'------------------------------------------------
sub Add i,j 'recursive  'w w2 global
'stores trail as graphic commands in global trail$
x=j-w2
y=i-w2
z=e(i,j)+1
xx=sx(x,y,z)
yy=sy(x,y,z)
trail$=trail$;";set ";int(xx);" ";int(yy)

'print i,j,x,y,cx,cy,sz,xx,yy

   for d = 1 to 4'direction
       'valid?
       if i+dx(d)>=0 and i+dx(d)<=w and j+dy(d)>=0 and j+dy(d)<=w then
           'diff>1?
           if e(i,j)-e(i+dx(d),j+dy(d))>0 then dirs$=dirs$;d
       end if
   next

   if len(dirs$)=0 then e(i,j)=e(i,j)+1: exit sub 'no way to spill
'actually it should fall off borders - not realised
   if isRandJump then
randN=int(rnd(0)*len(dirs$)+1)
else
randN=1    'always first
end if
   d=val(mid$(dirs$,randN,1))
   call Add i+dx(d),j+dy(d)    'recursive
end sub

sub AbelianAdd i,j 'recursive  'w w2 global
   e(i,j)=e(i,j)+1
   if e(i,j)<4 then exit sub
x=j-w2
y=i-w2
z=e(i,j)+1
xx=sx(x,y,z)
yy=sy(x,y,z)
trail$=trail$;";set ";int(xx);" ";int(yy)

'print i,j,x,y,cx,cy,sz,xx,yy

   e(i,j)=e(i,j)-4
   for d = 1 to 4'direction
       'valid?
       if i+dx(d)>=0 and i+dx(d)<=w and j+dy(d)>=0 and j+dy(d)<=w then
           call  AbelianAdd i+dx(d),j+dy(d)    'recursive
       end if
   next
end sub

'3d to screen
'---------------------------------
function sx(x,y,z)  'screen x
if is2d then sx=int(cx+sz*x):exit function
if isDimetry then
sx=int(cx+sz*(x+y/sqr(2)))
else
sx=int(cx+sz*((x+y)/sqr(3)*2))
end if
end function

function sy(x,y,z)  'screen y
if is2d then sy=int(cx+sz*y):exit function
if isDimetry then
sy=int(cy+sz*(0-z+y/sqr(2)))
else
sy=int(cy+sz*(0-z+(y-x)/sqr(3)))
end if
end function

' General subs/functions ====================
'------------------------------------------------
function iif(test, valYes, valNo)
   iif =  valNo
   if test then iif =  valYes
end function

function iif$(test, valYes$, valNo$)
   iif$ =  valNo$
   if test then iif$ =  valYes$
end function

'---------------------------------------------
' 0..1 into red-green-blue-red continuous colors
function rainbow$(x)
   hi = int((x*6) mod 6)+ 5*(x<0)  'fixed to 0..5
   f = (x*6) mod 1 + (x<0)         'frac, 0..1
   q = (1-f)
   select case hi
   case 0
       r = 1: g = f: b = 0
   case 1
       r = q: g = 1: b = 0
   case 2
       r = 0: g = 1: b = f
   case 3
       r = 0: g = q: b = 1
   case 4
       r = f: g = 0: b = 1
   case 5
       r = 1: g = 0: b = q
   end select
   R = int(r*255)
   G = int(g*255)
   B = int(b*255)
   rainbow$= R;" ";G;" ";B
end function

[carlg]

Use this code to produce a rendering of the Just BASIC logo icon.

[tenochtitlanuk]

Easier to produce multiple frames by displaying the JB logo then erasing it in random place, saving every so often. Then create an animated GIF either with photoshop/GIMP or one of the free online converters,

[tsh73]

Ok I get it to some finished point ("no bugs I am aware of", "no ideas as of how make better")
It's a toy - go ahead and play with it
If you set No Draw it will calculate fast but after finishing it will draw final picture.
(if you push Redraw you will get instant snapshot)
If you want add more, increase N and press Pause/Cont
To change picture while in finished/paused mode, change parameters and hit Redraw.
Checkboxes under groupbox radiobuttons relate only to one corresponding selection
Abelian sand mode works recursively, then I set number of points N more then 2000, it died with out of memory.

'   Sandpile (Abelian too)
' tsh73 Nov 2021
' v.1.1

   nomainwin

global w2, w
global cx, cy, sz
global isAbelian
global isShowGrid, nxtDraw
global isShowTrail, trail$
global is2d, isDimetry
global isTopLev 'only 3d
global isRandJump, isSpill 'non-Abelian

   WindowWidth = 608
   WindowHeight = 650

   UpperLeftX=int((DisplayWidth-WindowWidth)/2)
   UpperLeftY=int((DisplayHeight-WindowHeight)/2)

   graphicbox #main.gr, 5, 5, 350, 350
   statictext #main.statictext18, "Log", 5, 355, 144, 20
   texteditor #main.log, 5, 380, 350, 190

   button #main.btnStart, "Start Anew", [StartClick], UL, 366, 11, 122, 25
   button #main.btnPause, "Pause/Cont", [PauseClick], UL, 366, 46, 122, 25
   button #main.btnRedraw, "Single step", [SingleStepClick], UL, 494, 11, 88, 25
   button #main.btnRedraw, "Redraw", [RedrawClick], UL, 494, 46, 88, 25

   statictext #main.statictext5, "N", 366, 86, 25, 20
   textbox #main.txtN, 406, 76, 64, 25
   statictext #main.statictext7, "curr", 366, 111, 40, 20
textbox #main.txtCurr, 406, 106, 64, 25
   statictext #main.statictext21, "half size", 510, 81, 64, 20
   textbox #main.txtW2, 510, 106, 64, 25

   checkbox #main.chkNoDraw, "No Draw", [noDrawSet], [noDrawReset], 366, 135, 112, 20
   checkbox #main.chkGrid, "Floor grid", [gridSet], [gridReset], 366, 160, 96, 20
   checkbox #main.chkTrail, "Sand trail", [trailSet], [trailReset], 366, 185, 96, 20

   groupbox #main.groupbox13, "Sand style", 366, 215, 210, 50
   radiobutton #main.rbSandNormal, "Normal", [sandNormalSet], [dummy], 375, 235, 85, 20
   radiobutton #main.rbSandAbelian, "Abelian", [sandAbelianSet], [dummy], 460, 235, 85, 20
   checkbox #main.chkSpill "Spill over", [spillSet], [spillReset], 375, 270, 120, 20  
   checkbox #main.chkRand, "Random jump", [randSet], [randReset], 375, 295, 120, 20  

   groupbox #main.groupbox12, "View", 366, 330, 210, 60
   radiobutton #main.rb2d, "2d", [2dSet], [dummy], 375, 355, 40, 20
   radiobutton #main.rb3d, "3d", [3dSet], [dummy], 460, 355, 40, 20

   groupbox #main.groupbox12, "3d projection", 460, 400, 115, 80
   radiobutton #main.rbDimetry, "Dimetric", [dimetrySet], [dummy], 470, 425, 80, 20
   radiobutton #main.rbIsomery, "Isometric", [isometrySet], [dummy], 470, 450, 88, 20
   checkbox #main.chkTopLev, "Top level only ", [topLevSet], [topLevReset], 460, 490, 136, 20
   
menu #main, "Edit"    '<--- Texteditor Menu can be moved but not removed.

   open "Sandpile" for window as #main
   print #main, "trapclose [quit.main]"
   print #main.gr,"down; fill white; flush"
#main.gr "home; posxy cx cy"

   print #main, "font ms_sans_serif 10"
#main.txtCurr "!disable"
gosub [init]


   wait

[quit.main]
   Close #main
   END

'chkBoxes/radioButtons ==========

[2dSet]
is2d=1
gosub [initCalc] 'recalc coords
   wait

[3dSet]
is2d=0
gosub [initCalc]
   wait

[dimetrySet]
isDimetry=1
gosub [initCalc]
   wait

[isometrySet]
isDimetry=0
gosub [initCalc]
   wait

[sandAbelianSet]
isAbelian=1
   wait

[sandNormalSet]
isAbelian=0
   wait

[gridSet]
isShowGrid=1
   wait

[gridReset]    
isShowGrid=0
   wait

[trailSet]  
isShowTrail=1
   wait

[trailReset]
isShowTrail=0
   wait

[randSet]
isRandJump=1
   wait

[randReset]
isRandJump=0
   wait

[topLevSet]
isTopLev=1
   wait

[topLevReset]
isTopLev=0
   wait

[spillSet]
isSpill=1
   wait

[spillReset]
isSpill=0
   wait

[noDrawSet]
nxtDraw=0
   wait

[noDrawReset]
nxtDraw=1000000 'big enough
   wait

[dummy]    'non-used radio reset handler
   wait

'buttons =====================

[RedrawClick]
gosub [initCalc]
if nxtDraw<=0 then nxtDraw=1
gosub [doDraw]
   wait

[StartClick]
doWork=1
doPause=0
frame=0
#main.txtCurr frame
#main.log "started ";time$()
gosub [initCalc]
redim e(maxW,maxW) 'clear
goto [doFrame] 'will call itself with a timer
   wait

[PauseClick]
gosub [initCalc]
if doWork=0 and frame<=N then doWork=1 :doPause=1
if doPause then
doPause=0
#main.log "continued ";time$()
goto [doFrame]
else
if doWork=1 then
doPause=1
#main.log "paused ";time$()
end if
end if
   wait

[SingleStepClick]
gosub [initCalc]
goto [doFrame] 'in pause mode that does single frame
wait

'returnable SUBs ======================
[init]
w2=8
#main.txtW2 w2
w=2*w2

N=300
#main.txtN N

is2d=0
#main.rb2d iif$(is2d,"set","reset")
#main.rb3d iif$(not(is2d),"set","reset")

isDimetry=1
#main.rbDimetry iif$(isDimetry,"set","reset")
#main.rbIsomery iif$(not(isDimetry),"set","reset")

isShowGrid=1
#main.chkGrid  iif$(isShowGrid,"set","reset")

isTopLev=1
#main.chkTopLev iif$(isTopLev,"set","reset")

isAbelian=0
#main.rbSandAbelian iif$(isAbelian,"set","reset")
#main.rbSandNormal iif$(isAbelian,"reset","set")

isShowTrail=1
#main.chkTrail iif$(isShowTrail,"set","reset")

isRandJump=1
#main.chkRand iif$(isRandJump,"set","reset")

nxtDraw=1000000 'big or 0
#main.chkNoDraw iif$(nxtDraw<=0,"set","reset")

maxW=100 'so we can increase halfW and did not break
dim e(maxW,maxW)  'sand level
'directions array
dim dx(4):dx(1)=1:dx(2)= 0:dx(3)=-1:dx(4)=0
dim dy(4):dy(1)=0:dy(2)= 1:dy(3)= 0:dy(4)=-1
return

[initCalc]
#main.txtN "!contents? N"

if isAbelian then
maxCol=3
else
i = 1
ttl =1
while ttl <N
s=(i-1)^2+i^2
ttl=ttl+s
i=i+1
wend
maxCol=i-3 'counted from N, but -3 is trial and error
if maxCol< 1 then maxCol =1
'notice maxCol;" ";N;" ";ttl

end if

#main.txtW2 "!contents? w2"
w=2*w2

if is2d then 'fill the place
sz=int(2*cx/(w+2))
ptSize=sz-2
else
if isDimetry then 'trial and error factor
sz=int(2*cx/(w+2)/1.5)
ptSize=sz-4
else
sz=int(2*cx/(w+2)/2)
ptSize=sz
end if
end if
r=int(ptSize/2): if r<1 then r=1
return

[doDraw]
if nxtDraw<=0 then return
nxtDraw=nxtDraw-1
#main.gr "cls"
if isShowGrid then
#main.gr "size 1; color black"
z=0.5
for i = 0-w2 to w2
' x=j
y=i
xx1=sx(0-w2,y,z)
yy1=sy(0-w2,y,z)
xx2=sx(w2,y,z)
yy2=sy(w2,y,z)
#main.gr "line ";xx1;" ";yy1;" ";xx2;" ";yy2
next

for j= 0-w2 to w2
x=j
'y=i
xx1=sx(x,0-w2,z)
yy1=sy(x,0-w2,z)
xx2=sx(x,w2,z)
yy2=sy(x,w2,z)
#main.gr "line ";xx1;" ";yy1;" ";xx2;" ";yy2
next
end if

' #main.gr "size ";ptSize
#main.gr "size 1; color black"
for i = 0-w2 to w2
   'SCAN
for j= 0-w2 to w2
x=j
y=i
zz=e(i+w2,j+w2)
if zz then 'and we skip 0
if isTopLev or is2d then minZZ=zz else minZZ=1
for z = minZZ to zz
'print using("##",e(i+w2,j+w2));
if isAbelian then
colr$=word$("blue cyan pink",z)
else
c=(1-(z-1)/maxCol)*2/3
colr$=rainbow$(c)
end if
#main.gr "backcolor ";colr$
xx=sx(x,y,z)
yy=sy(x,y,z)
#main.gr "place ";xx;" ";yy;";circlefilled ";r


next
end if
next
'print
next

return

' timer sub ==========================
[doFrame]
timer 0 'just for a case
#main.txtCurr frame
if frame>N then
doWork=0
#main.log "finished ";time$()
GOTO [RedrawClick]
wait
end if
trail$="color red"
if isAbelian then
call AbelianAdd w2,w2
else
call Add w2,w2
end if
gosub [doDraw]
if isShowTrail and (nxtDraw>0) then #main.gr "size ";ptSize;";";trail$

if doWork and not(doPause) then
if nxtDraw>0 then
timer 100, [doFrame] 'recursively calls itself
else
timer 1, [doFrame] 'recursively calls itself
end if
end if
frame=frame+1
   wait

' Specific subs/functions ====================
'------------------------------------------------
sub Add i,j 'recursive  'w w2 global
'stores trail as graphic commands in global trail$
x=j-w2
y=i-w2
z=e(i,j)+1
xx=sx(x,y,z)
yy=sy(x,y,z)
trail$=trail$;";set ";int(xx);" ";int(yy)

'print i,j,x,y,cx,cy,sz,xx,yy

   for d = 1 to 4'direction
       'valid?
       if i+dx(d)>=0 and i+dx(d)<=w and j+dy(d)>=0 and j+dy(d)<=w then
           'diff>1?
           if e(i,j)-e(i+dx(d),j+dy(d))>0 then dirs$=dirs$;d
else
if isSpill and e(i,j)>0 then dirs$=dirs$;d
       end if
   next

   if len(dirs$)=0 then e(i,j)=e(i,j)+1: exit sub 'no way to spill
'actually it should fall off borders - not realised
   if isRandJump then
randN=int(rnd(0)*len(dirs$)+1)
else
randN=1    'always first
end if
   d=val(mid$(dirs$,randN,1))
'valid?
if i+dx(d)>=0 and i+dx(d)<=w and j+dy(d)>=0 and j+dy(d)<=w then
call Add i+dx(d),j+dy(d)    'recursive
'else spill - just forget that particle
end if
end sub

sub AbelianAdd i,j 'recursive  'w w2 global
   e(i,j)=e(i,j)+1
   if e(i,j)<4 then exit sub
x=j-w2
y=i-w2
z=e(i,j)+1
xx=sx(x,y,z)
yy=sy(x,y,z)
trail$=trail$;";set ";int(xx);" ";int(yy)

'print i,j,x,y,cx,cy,sz,xx,yy

   e(i,j)=e(i,j)-4
   for d = 1 to 4'direction
       'valid?
       if i+dx(d)>=0 and i+dx(d)<=w and j+dy(d)>=0 and j+dy(d)<=w then
           call  AbelianAdd i+dx(d),j+dy(d)    'recursive
       end if
   next
end sub

'3d to screen
'---------------------------------
function sx(x,y,z)  'screen x
if is2d then sx=int(cx+sz*x):exit function
if isDimetry then
sx=int(cx+sz*(x+y/sqr(2)))
else
sx=int(cx+sz*((x+y)/sqr(3)*2))
end if
end function

function sy(x,y,z)  'screen y
if is2d then sy=int(cx+sz*y):exit function
if isDimetry then
sy=int(cy+sz*(0-z+y/sqr(2)))
else
sy=int(cy+sz*(0-z+(y-x)/sqr(3)))
end if
end function

' General subs/functions ====================
'------------------------------------------------
function iif(test, valYes, valNo)
   iif =  valNo
   if test then iif =  valYes
end function

function iif$(test, valYes$, valNo$)
   iif$ =  valNo$
   if test then iif$ =  valYes$
end function

'---------------------------------------------
' 0..1 into red-green-blue-red continuous colors
function rainbow$(x)
   hi = int((x*6) mod 6)+ 5*(x<0)  'fixed to 0..5
   f = (x*6) mod 1 + (x<0)         'frac, 0..1
   q = (1-f)
   select case hi
   case 0
       r = 1: g = f: b = 0
   case 1
       r = q: g = 1: b = 0
   case 2
       r = 0: g = 1: b = f
   case 3
       r = 0: g = q: b = 1
   case 4
       r = f: g = 0: b = 1
   case 5
       r = 1: g = 0: b = q
   end select
   R = int(r*255)
   G = int(g*255)
   B = int(b*255)
   rainbow$= R;" ";G;" ";B
end function