be wheelsimulation
	sim 30 10
end

to sim rnum pensize_
	width=5
	N=rnum*width
	partx=(array N 1)
	partv=array partx
	parta=array partx
	a=array partx
	partcolor=array partx
	partbounces=array partx
	noforce=Array partx
	trails=false

	perspective
	setUpdateGraph false
	setPointSize pensize_
;	disablePointSmooth
	enablePointSmooth
;	enableDepthTest
	enableShadows
;	disableLighting
	cs ht pu setpc [50 230 20] 
	setx 0 sety -200-pensize_
	pd rt 90 fd 1000 pu home pd
	OnMouseLeftDown [trails=true]
	OnMouseRightDown 
	[	cs ht pu setpc [50 230 20] 
		setx 0 sety -200-pensize_
		pd rt 90 fd 1000 pu home pd 
		trails=false
	]
	dopt=40
	r=dopt*rnum/(2*pi)
	v0=2
	angle=20
	repeat N
	[	j=repcount
		k=Int (j-1)/rnum
		l=(mod (j-1) rnum)/rnum
;(pr j k l)
		partx.j=Float (List 
			-(r*cos 360*l)*(sin angle)+k*dopt
			(r*cos 360*l)*(cos angle)-50
			r*sin 360*l
			)
		partv.j=(List 0.0 0.0 0.0) 
			;(List v0*cos 360*l -v0*sin 360*l 0.0)
		parta.j=(List 0.0 -0.001 0.0)
		a.j=Float (List 0.0 0.0 0.0)
		partbounces.j=Int 0
		partcolor.j=HSBA 360*j/N 1 1 1
		noforce.j=Int 1
	]
ignore[
	noforce1=Array noforce
	noforce1_=Array noforce
	noforce2=Array noforce
	noforce2_=Array noforce
	noforcernum=Array noforce
	noforcernum_=Array noforce
	noforcernum1=Array noforce
	noforcernum1_=Array noforce
	noforce2rnum=Array noforce
	noforce2rnum_=Array noforce
	for [i 0 width-1]
	[	noforce1.i*rnum+1=Int 0
		noforce1_.(i+1)*rnum=Int 0
		noforce2.i*rnum+1=Int 0
		j=mod (i+1)*rnum N
		if j > 0
		[	noforce2_.j=Int 0
		]
		noforce2.i*rnum+2=Int 0
		j=mod (i+1)*rnum-1 N
		noforce2_.j=Int 0
		noforcernum.1+i*rnum=Int 0
		noforcernum_.N-i*rnum=Int 0
		noforce2rnum.1+i*rnum=Int 0
		noforce2rnum_.N-i*rnum=Int 0
	]
	noforcernum1.1=Int 0
	noforcernum1_.N-rnum=Int 0
	noforcernum1.(rnum+1)=Int 0
	noforcernum1_.N=Int 0
]
;	norefresh
	friction=0.9
	fn=10000
	force=(Array fn+1 0)
	fac=0.02
	ffein=10
	dopf=10*dopt
	for [i 4 fn]
	[	force.i=
			saturateBelow -1 
			saturateAbove 1
			fac*( 1*((i/dopf)^(-8))-2*(i/dopf)^(-4))
	]
	for [i 0 3]
	[	force.i=0
	]
	pd setpc 0 setXY rSeq -400 400 fn tolist force*1000 pu 
	updateGraph 
;stop
	grass=loadImage "grass.jpg
	texGrass=Texture grass
	disableTexture

	eye=Array 3
	phi=25
	theta=5
	dtheta=1
	center={0 0 0}
	upvector={0 1 0}
	dphi=1
	ddphi=dphi/3
	rotatescene_r=800
	dr=1.1
	onCharHandler
	OnChar [onCharHandler]
	video=false
;	video=true
	if video [(VideoStart "wheelsimulation 30)]
	running=true
	while [running]
	[	a*=Int 0
		
		dx=(partx-rotate partx 1)
		f=force.Int (saturateAbove fn (Norm dx)*ffein)
		a=a+dx*f ;*noforce1

		dx=(partx-rotate partx 2)
		f=force.Int (saturateAbove fn (Norm dx)*ffein)
		a=a-dx*f ;*noforce2

		dx=(partx-rotate partx rnum)
		f=force.Int (saturateAbove fn (Norm dx)*ffein)
		a=a+dx*f ;*noforcernum
	
		dx=(partx-rotate partx rnum+1)
		f=force.Int (saturateAbove fn (Norm dx)*ffein)
		a=a+dx*f ;*noforcernum1
	
		dx=(partx-rotate partx rnum*2)
		f=force.Int (saturateAbove fn (Norm dx)*ffein)
		a=a-dx*f ;*noforce2rnum
	
		dx=(partx-rotate partx -1)
		f=force.Int (saturateAbove fn (Norm dx)*ffein)
		a=a+dx*f ;*noforce1_

		dx=(partx-rotate partx -2)
		f=force.Int (saturateAbove fn (Norm dx)*ffein)
		a=a-dx*f ;*noforce2_

		dx=(partx-rotate partx -rnum)
		f=force.Int (saturateAbove fn (Norm dx)*ffein)
		a=a+dx*f ;*noforcernum_
	
		dx=(partx-rotate partx -rnum-1)
		f=force.Int (saturateAbove fn (Norm dx)*ffein)
		a=a+dx*f ;*noforcernum1_
	
		dx=(partx-rotate partx -rnum*2)
		f=force.Int (saturateAbove fn (Norm dx)*ffein)
		a=a-dx*f ;*noforce2rnum_
ignore [	
		for [i 0 width-1]
		[	j=1+i*rnum
			k=(i+1)*rnum	;1
			dx=(partx.j-partx.k)
			df=dx*(force.Int (saturateAbove fn (Norm dx)*ffein))
			a.j=a.j+df
			a.k=a.k-df
;ignore[
			j=1+i*rnum
			k=(i+1)*rnum-1	;2
			dx=(partx.j-partx.k)
			df=dx*(force.Int (saturateAbove fn (Norm dx)*ffein))
			a.j=a.j-df
			a.k=a.k+df
;ignore[
			j=2+i*rnum
			k=(i+1)*rnum	;-2
			dx=(partx.j-partx.k)
			df=dx*(force.Int (saturateAbove fn (Norm dx)*ffein))/2
			a.j=a.j-df
			a.k=a.k+df
;];ignore[
			j=1+i*rnum
			k=1+(i+1)*rnum	;rnum
			if k <= N
			[	dx=(partx.j-partx.k)
				df=dx*(force.Int (saturateAbove fn (Norm dx)*ffein))
				a.j=a.j+df
				a.k=a.k-df
			]
;ignore[
			j=1+i*rnum
			k=1+mod (i+width-1)*rnum N	;rnum+1
			dx=(partx.j-partx.k)
			df=dx*(force.Int (saturateAbove fn (Norm dx)*ffein))
			a.j=a.j+df
			a.k=a.k-df
;]
			j=1+i*rnum
			k=1+(i+2)*rnum	;rnum*2
			if k <= N
			[	dx=(partx.j-partx.k)
				df=dx*(force.Int (saturateAbove fn (Norm dx)*ffein))
				a.j=a.j-df
				a.k=a.k+df
			]
;ignore [
			j=1+i*rnum
			k=(i+width-2)*rnum+1	;-rnum*2
			if k <= N 
			[	dx=(partx.j-partx.k)
				df=dx*(force.Int (saturateAbove fn (Norm dx)*ffein))
				a.j=a.j-df
				a.k=a.k+df
			]
;]
		]
;]
;ignore [
		for [i 0 width-2]
		[	j=1+i*rnum
			k=1+(i+1)*rnum
			dx=(partx.j-partx.k)
			df=dx*force.Int (saturateAbove fn (Norm dx)*ffein)
			a.j=a.j+df
			a.k=a.k-df
		]
;]
]
		partv=partv+a
		partv=partv*friction
		partv=partv+parta
		partx+=partv
		repeat N
		[	j=repcount
			if (abs partx.j.2) > 300
			[	partx.j.2=partx.j.2-partv.j.2
				partv.j.2=-partv.j.2*0.95 
				partv.j.1=partv.j.1*0.95
				partv.j.3=partv.j.3*0.95
;				partbounces.j=partbounces.j+1
;				if partbounces.j==200
;				[	partx.j=Float (List -300 rnd*100 0)
;					partv.j=(List rnd/100 (rnd+1)/100 (rnd-0.5)/100)
;					partbounces.j=Int 0
;				]
			] 
		]
		clearScreen 
		clearShadows
		setEye eye center upvector
		draw_plane
		setPixel partx partcolor
;ignore [
		PenDown
		SurfaceStart
		i=1
		for [y 1 rnum]
		[	for [x 1 width]
			[	setPC partcolor.i
				setPosXYZ partx.i
				i=i+1
			]
			SurfaceColumn
		]
		SurfaceEnd
		partx+=(List 0.0 2.0 0.0)
		SurfaceStart
		for [y 1 rnum]
		[	i=y
			for [x 1 width]
			[	setPC partcolor.i
				setPosXYZ partx.i
				i=i+rnum
			]
			SurfaceColumn
		]
		SurfaceEnd
		PenUp
		partx-=(List 0.0 2.0 0.0)
;]
;		if Key? [break]
		castShadows
		updateGraph
		if video [VideoFrame]
		dispatchMessages
		GC
	]
	if video [VideoEnd]
end

be draw_plane
	horizon=10000
	PenUp
	Home
	setY -301 down 90
	fd horizon  rt 90  fd horizon  rt 90
	setPenColor HSB 60 0.3 0.7
	PenDown  
	enableTexture
	PolyStart  
		setTexXY   0 300 fd horizon*2  rt 90
		setTexXY 300 300 fd horizon*2  rt 90
		setTexXY 300   0 fd horizon*2  rt 90
		setTexXY   0   0 fd horizon*2  rt 90
	PolyEnd
	disableTexture
end

to onCharHandler
	ch=KeyboardValue
	if ch==ASCII "a 
	[	repeat N
		[	j=repcount
			k=Int (j-1)/rnum
			l=(1+mod j rnum)/rnum
			partv.j=(List 0.0 v0*sin 360*l -v0*(1+cos 360*l))
			p=partv.j
			s=sin angle
			c=cos angle
			p.2= c*p.2+s*p.3
			p.3=-s*p.2+c*p.3
		]
		friction=0.999
	]
	if ch==wxk_escape	[OnChar [] running=false]
	if ch==wxk_return	[onePoint=not onePoint]
	if ch==wxk_right [phi=phi+dphi]
	if ch==wxk_left  [phi=phi-dphi]
	if ch==wxk_up    [theta=theta+dtheta]
	if ch==wxk_down  [theta=theta-dtheta]
	if ch==wxk_prior [rotatescene_r=rotatescene_r/dr]
	if ch==wxk_next  [rotatescene_r=rotatescene_r*dr]
	eye.1=rotatescene_r*(cos theta)*sin phi
	eye.2=rotatescene_r* sin theta
	eye.3=rotatescene_r*(cos theta)*cos phi
	setLightPos {1000 1000 1000}
	setEye eye center upvector
	redraw
	updateGraph
end