be molecules3dc mol end be mol max_ = 1000 min_ = 1 maxb = 6 dopt = 35 dbind = trunc 2*dopt radius = dopt/8 ffein = 100 Temperature=300 cfein = 4000000/Temperature dopf = dopt*ffein*1.07 ;calibration dopf2 = dopt*ffein*0.8 maxf = ffein*dbind tE = 1.2 anfE = 0.4 expo = 5 ep = 1 fac = 1 fac2 = fac*0.001 gravV = -0.1*fac mov = 1*fac vfac = 50 cmin = 20 deltaTFac = 1.5 sqrDeltaTFac= deltaTFac*deltaTFac gravity=false threeDimensional=1 radx=dopt/2 rady=dopt/2 radz=dopt/2 sizeX=2 sizeY=2 c=IntArray max_ r=Array 3 ri=Array 3 v=Array 3 a=Array 3 for [k 1 3] [ r.k=FloatArray max_ ri.k=IntArray max_ v.k=FloatArray max_ a.k=FloatArray max_ ] banz=IntArray max_ b=Array max_ for [i 1 max_] [ b.i=IntArray maxb ] f=(FloatArray maxf+1 0) f2=(FloatArray maxf+1 0) sizehx=Int 400/dbind+1 sizex=2*sizehx+1 sizehy=Int 300/dbind+1 sizey=2*sizehy+1 sizehz=Int 300/dbind+1 sizez=2*sizehz+1 sizeh=Array 3 sizeh.1=sizehx sizeh.2=sizehy sizeh.3=sizehz m=(Array sizex -sizehx) for [mi -sizehx sizehx] [ mx=(Array sizey -sizehy) for [mj -sizehy sizehy] [ my=(Array sizez -sizehz) for [mk -sizehz sizehz] [ my.mk=[] ] mx.mj=my ] m.mi=mx ] onePoint=true tooSlow = false tooFast = false topteil = 0 disposalY= 0 clicked=false col=[] lineColor=RGB 1 0 1 bindColor=RGB 1 1 1 eye=Array 3 phi=0 theta=30 dtheta=1 center={0 0 0} upvector={0 1 0} dphi=1 ddphi=dphi/3 rotatescene_r=800 dr=1.1 setScreenColor 0 refresh setUpdateGraph false perspective 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=fac*( -((i/dopf)^(-expo))+(i/dopf)^(-expo-ep)) +fac*0.07*exp -sqr (i/maxf-0.6)*3 ] for [i 0 3] [ f2.i=0 ] end be squarehex 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 r.(1).i= x_+vxx*kx+vyx*ky r.(2).i= y_+vxy*kx+vyy*ky r.(3).i= 0 for [k 1 3] [ ri.k.i=sizeh.k+round r.k.i/dbind ] v.(1).i=vx_ v.(2).i=vy_ v.(3).i=0 banz.i=0 ifElse i <= max_ [ i=i+1 ][ print [Too many parts!] ] ] ] topteil=i end be squaresquare 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 ky=yi r.(1).i= x_+vxx*kx+vyx*ky r.(2).i= y_+vxy*kx+vyy*ky r.(3).i= 0 for [k 1 3] [ ri.k.i=sizeh.k+round r.k.i/dbind ] v.(1).i=vx_ v.(2).i=vy_ v.(3).i=0 banz.i=0 ifElse i <= max_ [ i=i+1 ][ print [Too many parts!] ] ] ] topteil=i end be cubesc side x_ y_ z_ 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) z_= z_-dopt local [i] i=1 for [yi 1 side] [ for [xi 1 side] [ for [zi 1 side] [ kx=xi ky=yi kz=zi r.(1).i= x_+vxx*kx+vyx*ky r.(2).i= y_+vxy*kx+vyy*ky r.(3).i= z_+dopt*kz for [k 1 3] [ ri.k.i=sizeh.k+round r.k.i/dbind ] v.(1).i=vx_ v.(2).i=vy_ v.(3).i=0 banz.i=0 ifElse i <= max_ [ i=i+1 ][ print [Too many parts!] ] ] ] ] topteil=i end be cubehcp side x_ y_ z_ 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) z_= z_-dopt local [i] i=1 for [yi 1 side] [ for [xi 1 side] [ for [zi 1 side] [ kx=xi+(mod yi 2)/2-(mod zi 2)/2 ky=yi*(Sqrt 3)/2-(mod zi 2)*(Sqrt 3)/4 kz=zi r.(1).i= x_+vxx*kx+vyx*ky r.(2).i= y_+vxy*kx+vyy*ky r.(3).i= z_+dopt*kz for [k 1 3] [ ri.k.i=sizeh.k+round r.k.i/dbind ] v.(1).i=vx_ v.(2).i=vy_ v.(3).i=0 banz.i=0 ifElse i <= max_ [ i=i+1 ][ print [Too many parts!] ] ] ] ] topteil=i end ifElse threeDimensional==1 [ ;cubehcp 5 0 0 0 30 0 90 cubesc 5 0 0 0 30 0 90 ][ ;squarehex 10 0 0 30 0 90 squaresquare 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 repeat count col [ i=repcount cc=reRGB col.i col.i=RGBA cc.1 cc.2 cc.3 0.5 ] end be movethem be faster local [k l] for [k 1 topteil] [ for [l 1 3] [ v.l.k=v.l.k*deltaTFac ] ] for [k 0 maxf] [ f.k=f.k*sqrDeltaTFac ] gravV:= gravV*sqrDeltaTFac end be slower local [k] for [k 1 topteil] [ for [l 1 3] [ v.l.k=v.l.k/deltaTFac ] ] for [k 0 maxf] [ f.k=f.k/sqrDeltaTFac ] gravV:= gravV/sqrDeltaTFac end be preparevars local [i k] if tooSlow [faster] if tooFast [slower] tooSlow=true tooFast=false for [i 1 topteil] [ for [k 1 3] [ a.k.i=0 ] c.i=0 ] end be energyloss local [h k] h=FloatArray 3 for [k 1 3] [ h.k=(v.k.i+v.k.j)/2 ] if i >= min_ [ for [k 1 3] [ v.k.i=h.k ] ] banz.i=banz.i+1 b.i=fput j b.i for [k 1 3] [ v.k.j=h.k ] banz.j=banz.j+1 b.j=fput i b.j 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 ] end be unboundf i j output not member? j b.i end local [i j bi di d da fx fy f0 force _c nomml] unbound=true da=FloatArray 3 Tag "nomml preparevars for [i 1 topteil] [ for [ix ri.(1).i-1 ri.(1).i+1] [ for [iy ri.(2).i-1 ri.(2).i+1] [ for [iz ri.(3).i-1 ri.(3).i+1] [ l=m.ix.iy.iz while [not empty? l] [ j=first l l=butFirst l da.1=r.(1).i-r.(1).j if (abs da.1) > dbind [continueLoop] da.2=r.(2).i-r.(2).j if (abs da.2) > dbind [continueLoop] da.3=r.(3).i-r.(3).j if (abs da.3) > dbind [continueLoop] d=Sqrt (Sqr da.1)+ (Sqr da.2)+ (Sqr da.3) if d > dbind [ ionize ] if d < dbind [ unbound=unboundf i j if unbound and2 (banz.i < maxb) and2 (banz.j < maxb) [ ;if not yet bound & free if (abs d-dopt) < 0.01 [ ;and d around dopt energyloss ;then "emitt a Photon" ] ] ] if unbound and2 (d > dopt) [ continueLoop ] d=d*ffein di=Int d if di >= maxf-1 [di=maxf-1] f0=f.di force=f0 ;+(d-Int d)*(f.(di+1)-f0) for [k 1 3] [ a.k.i=a.k.i+da.k*force ] ] ] ] ] ] for [i min_ topteil] [ c.i=Int (sqrt (sqr a.(1).i)+ (sqr a.(2).i)+ (sqr a.(3).i))*cfein if c.i > 250 [ tooFast=true tooSlow=false goto "nomml ] if c.i > cmin [ tooSlow=false ] ] for [i 1 topteil] [ for [k 1 3] [ v.k.i=v.k.i+a.k.i ] if gravity [ v.(2).i=v.(2).i+gravV ] rxi=ri.(1).i ryi=ri.(2).i rzi=ri.(3).i m.rxi.ryi.rzi=remove i m.rxi.ryi.rzi for [k 1 3] [ r.k.i=r.k.i+v.k.i ] if r.(1).i < -400+radx or2 r.(1).i > 400-radx [ v.(1).i=-v.(1).i r.(1).i=r.(1).i+v.(1).i ] if r.(2).i < -300+rady or2 r.(2).i > 300-rady [ v.(2).i=-v.(2).i r.(2).i=r.(2).i+v.(2).i ] if r.(3).i < -300+radz or2 r.(3).i > 300-radz [ v.(3).i=-v.(3).i r.(3).i=r.(3).i+v.(3).i ] for [k 1 3] [ ri.k.i=round r.k.i/dopt ] rxi=ri.(1).i ryi=ri.(2).i rzi=ri.(3).i m.rxi.ryi.rzi=fPut i m.rxi.ryi.rzi ] end be drawthem be draw i setXYZ r.(1).i r.(2).i r.(3).i setPC col.(c.i+2) Sphere radius setPC lineColor PenDown setXYZ r.(1).i+v.(1).i*vfac r.(2).i+v.(2).i*vfac r.(3).i+v.(3).i*vfac PenUp end cs for [i 1 topteil] [ draw i repeat maxb [ j=b.i.repcount if j > 0 [ setXYZ r.(1).i r.(2).i r.(3).i setPC bindColor PenDown setXYZ r.(1).j r.(2).j r.(3).j PenUp ] ] ] end be cooling local [i] for [i 1 topteil] [ for [k 1 3] [ v.k.i=v.k.i/tE ] ] end be heating local [i] for [i 1 topteil] [ for [k 1 3] [ v.k.i=v.k.i*tE ] ] end be stopthem local [i] for [i 1 topteil] [ for [k 1 3] [ v.k.i=0 ] ] end be findnearest hx hy local [i j dmin d] dmin=IntMax for [i 1 topteil] [ d=trunc Sqrt (Sqr hx-r.(1).i)+(Sqr hy-r.(2).i) if d < dmin [ dmin=d j=i ] ] output j end be showmark x y setPC 12 setXYZ x y z circle dopt/4 setPixelXYZ x y z 0 updateGraph setPC 0 setXYZ x y z circle dopt/4 end be mousepulling mou=MousePos if not clicked [ clicki=findnearest mou.1 mou.2 clicked=true ] i=clicki ; showmark(ox,oy); d=((Sqr mou.1-r.(1).i)+Sqr mou.2-r.(2).i)^0.3 for [k 1 3] [ v.k.i=v.k.i+mov*(mou.k-r.k.i)/d/dbind/te ; r.k.i=r.k.i+mov*(mou.k-r.k.i)/d ] end be mousespecials local [i mx my] mx=MouseX my=MouseY i=findnearest mx my showmark r.1.i r.2.i r.3.i for [i 1 3] [ v.k.i=0 ] r.1.i=radx r.2.i=rady+dopt*disposalY r.3.i=0 disposalY= Mod (disposalY+1) 6 while [MouseButtons!=0] [ dispatchMessages ] end be onCharHandler ch=KeyboardValue if ch==wxk_escape [OnChar [] running=false] if ch==wxk_return [onePoint=not onePoint] if ch==ASCII "- [cooling] if ch==ASCII "+ [heating] if ch==ASCII 0 [stopthem] if (ch==ASCII "g) or2 (ch==ASCII "G) [gravity=not gravity] if ch==ASCII "\ [clearScreen] 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 setEye eye center upvector setLightPos {1000 1000 1000} redraw updateGraph end init pd Line [[-400 0][400 0]] 1 setPixelXY rSeqFA -400 400 maxf+1 f*1000 2 setPixelXY rSeqFA -400 400 maxf+1 f2*1000 4 onCharHandler ignore readChar OnChar [onCharHandler] t0=0 running=true while [running] [ ;ignore[ Molecules3dMoveThem topteil maxf maxb r v a ri f f2 c banz b m gravV dbind dopt deltaTFac ffein cfein cmin gravity ;] movethem drawthem setXYZ -350 -250 0 setPC RGB 1 1 1 PenDown setHeading 90 t=TimeMilli Label list 1000/(t-t0) "fps t0=TimeMilli PenUp updateGraph ;updateVars ConsoleSetFocus dispatchMessages ifElse MouseButtons==1 [ mousepulling ][ ifElse MouseButtons==2 [ mousespecials ][ clicked=false ] ] ] end