code
code
time solver
time solver
In[]:=
ReflectTime[x1_,0,bounds_]=Infinity;
In[]:=
ReflectTime[x1_,v1_,bounds_]:=Max[Divide[2bounds[[1]]+1-2x1,2v1],Divide[2bounds[[2]]-1-2x1,2v1]]
In[]:=
RefreshReflectionTimes[state_,inds_]:=Module[{nextState=state},MapThread[Set[nextState["ReflectionTimes",#1],MapThread[ReflectTime,{#2["Position"],#2["Velocity"],state["BoundingBox"]}]]&,{inds,Lookup[state["PhaseSpace"],First/@inds]}];nextState]
In[]:=
CollideTime[_,vx1_,_,vx1_,_,_,_,_]:=Infinity
In[]:=
CollideTime[x1_,vx1_,x2_,vx2_,y1_,vy1_,y2_,vy2_]:=With[{dt=If[Or[Abs[x1-x2]<1,Sign[x1-x2]==Sign[vx1-vx2]],-1,If[x1>x2,Divide[(x1-x2)-1,vx2-vx1],Divide[(x2-x1)-1,vx1-vx2]]]},If[And[dt>=0,EuclideanDistance[y1+vy1*dt,y2+vy2*dt]<=1],dt,Infinity]]
In[]:=
RefreshCollisionTimes[state_,inds_]:=Module[{nextState=state},Map[Set[nextState["CollisionTimes",#],MapThread[CollideTime,Join[#,Reverse/@#]&[Catenate[Map[Lookup[#,{"Position","Velocity"}]&,Lookup[state["PhaseSpace"],#]]]]]]&,inds];nextState]
initialization
initialization
In[]:=
InitializeCHS[phaseSpace_,bounds_]:=Module[{inds,state=Association["Time"->0,"ExternalMomentum"->{0,0},"BoundingBox"->bounds,"ReflectionTimes"->Association[],"CollisionTimes"->Association[],"PhaseSpace"->Association[MapIndexed[Rule[#2[[1]],#1]&,phaseSpace]]]},inds=Range[Length[phaseSpace]];state=RefreshReflectionTimes[state,List/@inds];state=RefreshCollisionTimes[state,Subsets[inds,{2}]];state]
dynamics
dynamics
In[]:=
PropagateMass[element_Association,time_]:=Module[{nextElement=element},nextElement["Position"]=Plus[nextElement["Position"],Times[time,nextElement["Velocity"]]];nextElement]
In[]:=
PropagateState[state_,nextTime_]:=Module[{nextState=state},Map[Set[nextState["PhaseSpace",#],PropagateMass[nextState["PhaseSpace",#],nextTime]]&,Keys[nextState["PhaseSpace"]]];Set[nextState["ReflectionTimes"],nextState["ReflectionTimes"]-nextTime];Set[nextState["CollisionTimes"],nextState["CollisionTimes"]-nextTime];nextState["Time"]=state["Time"]+nextTime;nextState]
In[]:=
FindZeroTimes[eventTimes_]:=Association[KeyValueMap[Function[{key,val},If[SameQ[#,{}],Nothing,key->#]&[MapIndexed[If[PossibleZeroQ[#1],#2[[1]],Nothing]&,val]]],eventTimes]]
In[]:=
ResolveReflection[element_,inds_]:=Module[{nextElement=element},Set[nextElement["Velocity"],ReplacePart[element["Velocity"],#->-element["Velocity"][[#]]&/@inds]];nextElement]
In[]:=
ResolveReflections[state_,events_]:=Module[{nextState=state,dp},dp=Total[KeyValueMap[Function[{key,value},ReplacePart[ConstantArray[0,2],Rule[#,Times[2,state["PhaseSpace",First[key],"Mass"],state["PhaseSpace",First[key],"Velocity"][[#]]]]&/@value]],events]];Set[nextState["ExternalMomentum"],Plus[dp,state["ExternalMomentum"]]];KeyValueMap[Set[nextState["PhaseSpace",First[#1]],ResolveReflection[state["PhaseSpace",First[#1]],#2]]&,events];nextState]
In[]:=
ElasticCollision[m1_,m2_,u1_,u2_]:=Plus[Times[Divide[m1-m2,m1+m2],u1],Times[Divide[2m2,m1+m2],u2]]
In[]:=
ResolveCollision[element1_,element2_,inds_]:=Module[{nextElement=element1},nextElement["Velocity"]=ReplacePart[element1["Velocity"],Map[Rule[#,ElasticCollision[element1["Mass"],element2["Mass"],element1["Velocity"][[#]],element2["Velocity"][[#]]]]&,inds]];nextElement]
iterator and control
iterator and control
utilities
utilities
Basic Testing
Basic Testing
Critical Testing
Critical Testing
fail testing
fail testing