(* Planets: A celestial simulator Copyright (C) 2001-2003 Yaron M. Minsky This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA*)openStdLabelsopenMoreLabelsopenPrintfopenStateopenConstantsopenCommon(* Faster physics implementation that is array-based *)(* Only implements the act_all_on_all function. *)(* Another useful thing to implement faster would be collision detection *)(* all-float record, to let the compiler unbox in a loop *)typefbody={mutablex_pos:float;mutabley_pos:float;mutablex_vel:float;mutabley_vel:float;fb_radius:float;fb_mass:float;}letcubex=x*.x*.xletsquarex=x*.x(* let sqrtcube x = sqrt (x *. x *. x) (* Faster but less flexible implementaqtion of sqrtcube *) *)letbody_to_fbodybody=let(x_pos,y_pos)=body.posand(x_vel,y_vel)=body.velocity(*and (x_accel,y_accel) = match body.i with None -> (0.,0.) | Some accel -> accel *)in{x_pos=x_pos;y_pos=y_pos;fb_radius=body.radius;fb_mass=body.mass;x_vel=x_vel;y_vel=y_vel;}letbuild_fbody_arraybodies=Array.of_list(List.map~f:body_to_fbodybodies)(* convert a body and a b to an updated body *)letfbody_and_body_to_bodybodyb={bodywithvelocity=(b.x_vel,b.y_vel);pos=(b.x_pos,b.y_pos);}(******** List iterator w/index ********)letreclist_iteri_rec~flisti=matchlistwith[]->()|hd::tl->f~ihd;list_iteri_rec~ftl(i+1)letlist_iteri~flist=list_iteri_rec~flist0(**************************************************)letarray_to_bodiesbodiesarray=letb_list=Array.to_listarrayinList.map2~f:fbody_and_body_to_bodybodiesb_list(**********************************************************************)lettimer=MTimer.create()letsqrtcubex=sqrt(x*.x*.x)letact_all_on_all_rk~bouncebodies=letgexp=grav_exp#vinletexp=((1.0+.gexp)/.2.0)inletsqrtcube=ifgexp=2.0thensqrtcubeelse(funx->x**exp)inletconst=gconst#vin(* t is the time. This function has no time dependence. s is the position is state space dsdt is the array of derivatives *)letderiv_tsdsdt=(* initialize derivatives to 0 *)fori=0toArray.lengthdsdt-1dodsdt.(i)<-0.done;fori=0toArray.lengthbodies-1do(* x and y pos derivative *)dsdt.(i*4)<-s.(i*4+2);dsdt.(i*4+1)<-s.(i*4+3);done;fori=0toArray.lengthbodies-1doforj=i+1toArray.lengthbodies-1do(* compute i's action on j and vice-versa. That way you only need to * compute the force once. It's nearly twice as fast that way. *)letx_pos_i=s.(i*4)andy_pos_i=s.(i*4+1)inletx_pos_j=s.(j*4)andy_pos_j=s.(j*4+1)inletxdiff=x_pos_i-.x_pos_jandydiff=y_pos_i-.y_pos_jinletdist_sq=xdiff*.xdiff+.ydiff*.ydiffinletmult=const/.sqrtcubedist_sqinletmult_i=-.mult*.bodies.(j).fb_massandmult_j=mult*.bodies.(i).fb_massin(* x vel derivative *)dsdt.(i*4+2)<-dsdt.(i*4+2)+.xdiff*.mult_i;dsdt.(j*4+2)<-dsdt.(j*4+2)+.xdiff*.mult_j;(* y vel derivative *)dsdt.(i*4+3)<-dsdt.(i*4+3)+.ydiff*.mult_i;dsdt.(j*4+3)<-dsdt.(j*4+3)+.ydiff*.mult_j;ifbounce&&(lettotal_radius=bodies.(j).fb_radius+.bodies.(i).fb_radiusindist_sq<total_radius*.total_radius)then(lettotal_radius=bodies.(j).fb_radius+.bodies.(i).fb_radiusinletmult=-.mult*.bodies.(i).fb_mass*.bodies.(j).fb_mass*.(total_radius*.total_radius/.dist_sq)inletmult_i=-.mult/.bodies.(i).fb_massandmult_j=mult/.bodies.(j).fb_massin(* x vel derivative *)dsdt.(i*4+2)<-dsdt.(i*4+2)+.xdiff*.mult_i;dsdt.(j*4+2)<-dsdt.(j*4+2)+.xdiff*.mult_j;(* y vel derivative *)dsdt.(i*4+3)<-dsdt.(i*4+3)+.ydiff*.mult_i;dsdt.(j*4+3)<-dsdt.(j*4+3)+.ydiff*.mult_j;)donedoneinlets=Array.init(4*Array.lengthbodies)~f:(funi->matchimod4with|0->bodies.(i/4).x_pos|1->bodies.(i/4).y_pos|2->bodies.(i/4).x_vel|_->bodies.(i/4).y_vel)in(* dumb_solver s state.delta#v deriv; *)Rk4.step~y:s~x:0.0~h:state.delta#v~yout:s~derivs:deriv;fori=0toArray.lengthbodies-1dobodies.(i).x_pos<-s.(4*i);bodies.(i).y_pos<-s.(4*i+1);bodies.(i).x_vel<-s.(4*i+2);bodies.(i).y_vel<-s.(4*i+3);doneletact_all_on_all~bouncebodies=letarray=build_fbody_arraybodiesinact_all_on_all_rk~bouncearray;array_to_bodiesbodiesarray