In[]:=
HardSquareSimulation[initialPositions_,initialVelocities_,boxSize_,stepSize_,squareRad_,steps_]:=Module{collisions,trajectories,counter},(*Arraytokeeptrackofcollisionsastheyoccur*)collisions={};(*Arraytokeeptrackofthetimestep.Willjuststorecountingnumbers.*)(*k-that'sodd,why'sitalist.Iupdatedtovariable*)counter=0;(*NestListwilliterativelyapplyafunctiontotheoriginalinput.Inthiscase,theinputwillbetheinitialconditions,andthefunctionwillupdatethepositionsandvelocities.Storetrajectoriestoreturnattheend.*)trajectories=NestListFunction{states},Module{positions,velocities,colliding,overlap,p1,p2},positions=states[[1]];velocities=states[[2]];(*Incrementthecounter.*)counter++;(*Updatepositions*)positions+=stepSize*velocities;(*Mutualdistancechecks&scatteringforallparticles*)Do(*Findallparticlesintersectingwiththisone*)colliding=Select[Range[Length[positions]],#!=particle&&Max[Abs[Subtract@@positions[[{#,particle}]]]]<2squareRad&];(*Ifanyarecloseenough(alongbothaxes)toscatter,collidewitheachinanarbitraryorder*)Table(*Recordcollision*)AppendTo[collisions,{particle,other,counter}];(*Findthecornerpointswhichareoverlapping(Thiscouldprobablyusefewervariablesbutthisisclearer)*)p1=positions[[particle]]+squareRad*Sign/@Subtract@@positions[[{other,particle}]];p2=positions[[other]]+squareRad*Sign/@Subtract@@positions[[{particle,other}]];With(*Attempttofindwhichaxiscollidedfirst,bytracingthecornerpositionsbackintime,andseeingwhichcomponentbecameequalmostrecently.Ifithasn'tbecomeequal,itmustnotbecollidinginthisaxis,sosetto-∞*)dimension=First@OrderingIf[#>0,-∞,#]&/@Quiet@Check,1&@@@Transpose[{{p2,p1},stepSize*velocities[[{particle,other}]]},{2,3,1}],-1,Less,(*Calculatehowmuchtheparticlesareoverlappingtocorrecttheissue(anecessaryby-productofdiscretetimesteps)*)overlap=2squareRad-Abs[Subtract@@positions[[{other,particle},dimension]]];positions[[{particle,other},dimension]]+=*{-1,1}*If[positions[[particle,dimension]]<positions[[other,dimension]],1,-1];(*Swapthevelocitiesintheaxiswhichtheycollided*)velocities[[{particle,other},dimension]]=velocities[[{other,particle},dimension]];,{other,colliding},{particle,1,Length[positions]};(*Reflectionoffwalls.Ifparticleiswithin"rad"ofedge,reversethatcomponentofitsvelocity.Setitspositiontobeprecisely"rad"fromtheedge.Implementsthisprocedurecomponentwise.*)Do[If[Abs[positions〚particle,dimension〛]>boxSize-squareRad,velocities〚particle,dimension〛*=-1;positions〚particle,dimension〛=N[Sign[positions〚particle,dimension〛](boxSize-squareRad)];];,{dimension,1,Length[positions〚1〛]},{particle,1,Length[positions]}];{positions,velocities}(*Initialconditionsandnumberofsteps,forNestList.*),{initialPositions,initialVelocities},steps;(*Returnlist{trajectories,collisions}.*){trajectories,collisions}
Subtract@@#1
Subtract@@#2
overlap
2
In[]:=
Table[Graphics[{Style[Disk[#[[1]],1],Orange],Arrow[{#[[1]],#[[1]]+4#[[2]]}]}&/@Transpose[With[{num=50,steps=2000,boxSize=15,rad=1,stepSize=0.1},With[{initialPos=Catenate@Table[{-14+2.5i,14-2.5j},{i,4},{j,4}],initialVel=Table[.5{1,-1},16]},HardSquareSimulation[initialPos,initialVel,boxSize,stepSize,rad,steps]]][[1,t]]],PlotRange->15,Frame->True,FrameTicks->None,ImageSize->100],{t,1,2000,100}]
Out[]=
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
In[]:=
Table[Graphics[{Style[RegularPolygon[#[[1]],1,4],Orange],Arrow[{#[[1]],#[[1]]+4#[[2]]}]}&/@Transpose[With[{num=50,steps=2000,boxSize=15,rad=1,stepSize=0.1},With[{initialPos=3.2N[CirclePoints[10]],initialVel=N[CirclePoints[10]]},HardSquareSimulation[initialPos,initialVel,boxSize,stepSize,rad,steps]]][[1,t]]],PlotRange->15,Frame->True,FrameTicks->None,ImageSize->100],{t,1,1000,100}]
Out[]=
,
,
,
,
,
,
,
,
,
