be molecules3
	max_	=	400
	min_	=	1
	maxb	=	6

	dopt	=	25
	dbind	=	trunc 2*dopt

	ffein	=	100
	Temperature=300
	cfein	=	160000/Temperature
;	dopf	=	dopt*ffein
	dopf	=	dopt*ffein*0.7		;calibration
	dopf2	=	dopt*ffein*1.2
	maxf	=	ffein*dbind

	tE		=	1.1
	anfE	=	0.4

	expo	=	5
	ep		=	1
	ep2		=	0.2

	fac	=	0.25	;0.5
	fac2=	fac*0.1
	gravV=-0.01*fac
	mov	=	10*fac
	vfac=100

	cmin	=	20
	deltaTFac	=	1.5
	sqrDeltaTFac	=	deltaTFac ;*deltaTFac

	radx=4
	rady=3

	ox=FloatArray max_
	oy=FloatArray max_
	c=IntArray max_
	x=FloatArray max_
	y=FloatArray max_
	rx=IntArray max_
	ry=IntArray max_
	vx=FloatArray max_
	vy=FloatArray max_
	ax=FloatArray max_
	ay=FloatArray max_
	banz=IntArray max_
	b=Array max_
	for [i 1 max_]
	[	b.i=[]
	]
	
	f=(FloatArray maxf+1 0)
	f2=(FloatArray maxf+1 0)
	
	sizehx=Int 400/dopt+1
	sizex=2*sizehx+1
	sizehy=Int 300/dopt+1
	sizey=2*sizehy+1
	m=(Array sizex -sizehx)
	for [mi -sizehx sizehx]
	[	mx=(Array  sizey -sizehy)
		for [mj -sizehy sizehy]
		[	mx.mj=[]
		]
		m.mi=mx
	]

	onePoint=true
	gravity	=false

	tooSlow	= false
	tooFast	= false

	topteil = 0
	disposalY=	0

	col=[]
	lineColor=RGB 1 0 1
	bindColor=RGB 1 1 1
	setScreenColor 0
	norefresh
	setUpdateGraph false
	disableRoundLineEnds
	setPenSize [0 0]
	hideTurtle
	PenUp

	be init

		be initforcetable
;			setItems 0 f (rSeqFA 1 0 int maxf/2)^2/2
;			setItems int maxf/2 f (rSeqFA 0 1 int maxf/2)^2/100* -1
;		stop
			for [i 4 maxf]
			[	f.i=fac*( -((i/dopf)^(-expo))+(i/dopf)^(-expo-ep))
			]
			for [i 0 3]
			[	f.i=0
			]

			for [i 4 maxf]
			[	f2.i=fac2*( -((i/dopf2)^(-expo))+(i/dopf2)^(-expo-ep2))
			]
			for [i 0 3]
			[	f2.i=0
			]
		end
	
		be square side x_ y_ angle v vangle
			 angle= angle
			vangle=vangle

			vx_=v*Cos vangle
			vy_=v*Sin vangle

			vxx= dopt*Cos angle 
			vxy= dopt*Sin angle
			vyx= dopt*Sin angle
			vyy=-dopt*Cos angle

			kx=side/2+(mod trunc side/2 2)/2
			ky=side/2*(Sqrt 3)/2
			x_= x_-(vxx*kx+vyx*ky)
			y_= y_-(vxy*kx+vyy*ky)

			local [i]
			i=1
			for [yi 1 side]
			[	for [xi 1 side]
				[	kx=xi+(mod yi 2)/2
					ky=yi*(Sqrt 3)/2
					x.i= x_+vxx*kx+vyx*ky
					y.i= y_+vxy*kx+vyy*ky
					rx.i=round x.i/dopt
					ry.i=round y.i/dopt
					vx.i=vx_
					vy.i=vy_
					banz.i=0
					ifElse i <= max_
					[	i=i+1
					][	print [Too many parts!]
					]
				]
			]
			topteil=i
		end
	
		square 10 0 0 30 0 90
		topteil=topteil-1
	
		initforcetable
	
		setXY 0 -270  setH 90
		Label [[RETURN]=splines  [+]=heat  [-]=cool
		  [G]=gravity  [other Key]=cS  Mouse: L=pull R=del]
		col=loadpalette "TEILE.PAL
	end


	be movethem
	
		be faster
			local [k]
			for [k 1 topteil]
			[	vx.k=vx.k*deltaTFac
				vy.k=vy.k*deltaTFac
			]
			for [k 0 maxf]
			[	f.k=f.k*sqrDeltaTFac
				f2.k=f2.k*sqrDeltaTFac
			]
			gravV:= gravV*sqrDeltaTFac
		end
	
		be slower
			local [k]
			for [k 1 topteil]
			[	vx.k=vx.k/deltaTFac
				vy.k=vy.k/deltaTFac
			]
			for [k 0 maxf]
			[	f.k=f.k/sqrDeltaTFac
				f2.k=f2.k/sqrDeltaTFac
			]
			gravV:= gravV/sqrDeltaTFac
		end

		be preparevars
			local [i]
			if tooSlow [faster]
			if tooFast [slower]
	
			tooSlow=true
			tooFast=false
	
			for [i 1 topteil]
			[	ax.i=0
				ay.i=0
				c.i=0
			]
		end
	
		be energyloss
			local [hx hy]
			setPC RGB 1 1 1
			Line List List x.i y.i  List x.j y.j RGB 0 0 1
			hx=(vx.i+vx.j)/2
			hy=(vy.i+vy.j)/2
	
			if i >= min_
			[	vx.i=hx
				vy.i=hy
			]
			banz.i=banz.i+1
			b.i=fput j b.i
			vx.j=hx
			vy.j=hy
			banz.j=banz.j+1
			b.j=fput i b.j
	
			setPC 0
			Line List List x.i y.i  List x.j y.j RGB 0 0 1
		end
	
		be ionize
			if member? j b.i
			[	b.i=remove j b.i
				banz.i=banz.i-1
				b.j=remove i b.j
				banz.j=banz.j-1
(pr "i i j)
			]
		end
	
		be draw i
			setXY x.i y.i
			setFC col.(c.i+2)
			fillCircle dopt/4
			PenDown
			Line List 
				List x.i y.i 
				List x.i+vx.i*vfac y.i+vy.i*vfac lineColor
			PenUp
		end
	
		be del x y
			setXY x y
			setFC 0
		;	fillCircle dopt/2
		end
	
		be unboundf i j
			output not member? j b.i
		end
	
		local [i j bi di
			d dx dy
			fx fy f0 force _c nomml]

		unbound_=true
		tag "nomml
		preparevars
		for [i 1 topteil]
		[	for [ix rx.i-1 rx.i+1]
			[	for [iy ry.i-1 ry.i+1]
				[	l=m.ix.iy
					while [not empty? l]
					[	j=first l
						l=butFirst l
						dx= x.i-x.j
					;	if (abs dx) > dbind [continueLoop]
						dy= y.i-y.j
					;	if (abs dy) > dbind [continueLoop]
						d= Sqrt (Sqr dx)+(Sqr dy)
	
						if d > dopt*1.3
						[
							ionize
						]
						if d < dopt*1.1
						[
							unbound=unboundf i j
							if unbound 
							and2 (banz.i < maxb) 
							and2 (banz.j < maxb)
							[	;if not yet bound & free
								if (abs d-dopt)/dopt < 0.5
								[	;and d around dopt
									energyloss	;then "emitt a Photon"
(pr "e i j)
								]
							]	
						]
						if d > dbind
						[
							continueLoop
						]
						d=d*ffein
						di=Int d
						if di >= maxf-1 [di=maxf-1]
						ifelse unbound
						[	f0=f2.di
						][	f0=f.di
						]						
						force=f0 ;+(d-Int d)*(f.(di+1)-f0)

						fx=dx*force
						fy=dy*force
						ax.i=ax.i+fx
						ay.i=ay.i+fy
;						ax.j=ax.j-fx
;						ay.j=ay.j-fy
					]
				]
			]
		]

		for [i min_ topteil]
		[	c.i=Int (sqrt (sqr ax.i)+(sqr ay.i))*cfein
			if c.i > 250 
			[
				tooFast=true
				tooSlow=false
				goto "nomml
			]
			if c.i > cmin
			[	tooSlow=false
			]
		]
		clearScreen
		for [i 1 min_-1 1]
		[	draw i
		]	
		for [i min_ topteil]
		[	vx.i=vx.i+ax.i
			vy.i=vy.i+ay.i
	
			if gravity
			[	vy.i=vy.i+gravV
			]
			rxi=rx.i
			ryi=ry.i
			m.rxi.ryi=remove i m.rxi.ryi
			x.i=x.i+vx.i
			y.i=y.i+vy.i
	
			if x.i < -400+radx or2 x.i > 400-radx
			[	vx.i=-vx.i
				x.i=x.i+vx.i
			]
			if y.i < -300+rady or2 y.i > 300-rady
			[	vy.i=-vy.i
				y.i=y.i+vy.i
			]
			rx.i=round x.i/dopt
			ry.i=round y.i/dopt
			rxi=rx.i
			ryi=ry.i
			m.rxi.ryi=fPut i m.rxi.ryi
			
			ifElse onePoint
			[	draw i
				l=b.i
				while [not empty? l]
				[	j=first l
					l=butFirst l
					if j > 0
					[	setXY 
							x.i 
							y.i 
						setPC bindColor
						PenDown
						setXY 
							x.j 
							y.j 
						PenUp
					]
				]
			][	setPixelXY x.i y.i c.i+1
			]
		]
	end

	be cooling
		local [i]
		for [i 1 topteil]
		[	vx.i=vx.i/tE
			vy.i=vy.i/tE
		]
	end
	
	be heating
		local [i]
		for [i 1 topteil]
		[	vx.i=vx.i*tE
			vy.i=vy.i*tE
		]
	end
	
	be findnearest hx hy
		local [i j dmin d]
		dmin=IntMax
		for [i 1 topteil]
		[	d=trunc Sqrt (Sqr hx-x.i)+(Sqr hy-y.i)
			if d < dmin
			[	dmin=d
				j=i
			]
		]
		output j
	end
	
	be showmark x y
		setPC 12
		setXY x y
		circle dopt/4
		setPixelXY x y 0
		updateGraph
		setPC 0
		setXY x y
		circle dopt/4
	end
	
	be mousepulling
		mx=MouseX
		my=MouseY
		if not clicked
		[	clicki=findnearest mx my
			clicked=true
		]
		i=clicki
	;	showmark(ox,oy);
		d=((Sqr mx-x)+Sqr my-y)^0.3
		vx.i=0 ;(vx+mov*(mx-x)/d)/te
		vy.i=0 ;(vy+mov*(my-y)/d)/te
		x.i=x.i+mov*(mx-x.i)/d
		y.i=y.i+mov*(my-y.i)/d
	end
	
	be mousespecials
		local [i mx my]
		mx=MouseX
		my=MouseY
		i=findnearest mx my
		showmark ox.i oy.i
		vx.i=0
		vy.i=0
		x.i=radx
		y.i=rady+dopt*disposalY
		disposalY= Mod (disposalY+1) 6
		while [MouseButtons!=0]
		[	dispatchMessages
		]
	end


	init
	setPixelXY rSeqFA -400 400 maxf+1  f*1000 15
	setPixelXY rSeqFA -400 400 maxf+1  f2*1000 4
;stop
	forever
	[	movethem
;updateVars
		updateGraph
		dispatchMessages
		if Key?
		[	ch=lowerCase readChar
			if ch==Char 27	[break]
			if ch==Char 13 [onePoint=not onePoint]
			if ch=="- [cooling]
			if ch=="+ [heating]
			if ch=="g [gravity=not gravity]
			if ch=="  [clearScreen]
		]
		ifElse MouseButtons==1 
		[	mousepulling
		][	ifElse MouseButtons==2 
			[	mousespecials
			][	clicked=false
			]
		]
	]
	pr [End]
end