openState(***************************************************************)(** Collision detection *************************************)(***************************************************************)(* Computes smallest positive solution to the quadratic equation |p + t*v|^2 - dist = 0. The solution is: - v.p +/- sqrt((v.p)^2 - v.v*(p.p - dist)) ------------------------------------------ v.v if the sqrt term is imaginary, then bail out early (before computing the sqrt). This is a big performance win. *)letnext_collision_valuespvdistdotpp=letdotvp=dotvpinletdotvpsq=dotvp*.dotvpanddotpp=dotppanddotvv=dotvvinletpre_sqrt_term=dotvpsq-.dotvv*.(dotpp-.(dist*.dist))inifpre_sqrt_term<=0.theninfinityelseletsqrt_term=sqrtpre_sqrt_terminletsoln1=(-.dotvp-.sqrt_term)/.dotvvandsoln2=(-.dotvp+.sqrt_term)/.dotvvin(* It's important to ignore planets that are currently touching *)(* note that soln2 is always greater than soln1 *)ifsoln2<=0.theninfinityelseifsoln1>0.thensoln1elsesoln2(* computes time of next collision given two bodies *)letnext_collision_bodiesb1b2=letdist=b1.radius+.b2.radiusinletp=b1.pos<->b2.posinletv=b1.velocity<->b2.velocityinletdotpp=dotppinifdotpp<dist*.disttheninfinityelsenext_collision_valuespvdistdotpp(* computes nearest collision *)letnext_collisionbodies=letlen=Array.lengthbodiesinletbest=ref(infinity,None)infori=0tolen-2doforj=i+1tolen-1dolettime=next_collision_bodiesbodies.(i)bodies.(j)inlet(best_time,_)=!bestiniftime<best_timethenbest:=(time,Some(i,j));donedone;match!bestwith(time,None)->raiseNot_found|(time,Some(i,j))->(time,(i,j))(***************************************************************)(*** Bouncing ***********************************************)(***************************************************************)(* Computes the bounce of two touching bodies. It is assumed that the two bodies are touching. *)(* compute new velocities from 1-dimensional bounce problem *)letbounce_1d~s1~s2~m1~m2=letnew_s1=((m1-.m2)*.s1+.2.*.m2*.s2)/.(m1+.m2)andnew_s2=((m2-.m1)*.s2+.2.*.m1*.s1)/.(m1+.m2)in(new_s1,new_s2)letmagsqv=dotvvletmagv=sqrt(magsqv)letuvecv=magv<|>vletbounce_pairb1b2=letv1=b1.velocityandv2=b2.velocityandm1=b1.massandm2=b2.massandp1=b1.posandp2=b2.posinletu=uvec(p2<->p1)in(* speed in collision direction *)lets1=dotv1uands2=dotv2uin(* perpendicular components of velocity *)letv1p=v1<->(s1<*>u)andv2p=v2<->(s2<*>u)inletnew_s1,new_s2=bounce_1d~s1~s2~m1~m2inletnew_v1=(new_s1<*>u)<+>v1pandnew_v2=(new_s2<*>u)<+>v2pin({b1withvelocity=new_v1},{b2withvelocity=new_v2})(***************************************************************)(**** Moving Forward ****************************************)(***************************************************************)letmove_forward_simplebodiestime=fori=0toArray.lengthbodies-1doletbody=bodies.(i)inbodies.(i)<-{bodywithpos=body.pos<+>(time<*>body.velocity)}doneletrecmove_forward_arraytimebodies=trylet(ctime,(i,j))=next_collisionbodiesinifctime>timethenmove_forward_simplebodiestimeelse(move_forward_simplebodiesctime;(* move to collision time *)letb_i,b_j=bodies.(i),bodies.(j)inletb_i,b_j=bounce_pairb_ib_jinbodies.(i)<-b_i;bodies.(j)<-b_j;move_forward_array(time-.ctime)bodies)withNot_found->move_forward_simplebodiestimeletmove_forwardtimebodies=letbodies=Array.of_listbodiesinmove_forward_arraytimebodies;Array.to_listbodies