(* 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(*******************************************************)(*** Physics *****************************************)(*******************************************************)letcubex=x*.x*.xletsquarex=x*.xletmagsq(x,y)=x*.x+.y*.yletmag(x,y)=sqrt(x*.x+.y*.y)letdistsqv1v2=magsq(v1<->v2)letdistv1v2=mag(v1<->v2)(* return unit vector pointing from (x1,y1) to (x2,y2) *)letunit_vectv1v2=distv1v2<|>(v2<->v1)letupdate_posbody={bodywithpos=body.pos<+>(state.delta#v<*>body.velocity)}letupdate_pos_bodiesbodies=List.map~f:update_posbodiesletprint_bodybody=let(x,y)=body.posand(x_v,y_v)=body.velocityinprintf"pos: (%3f,%3f), vel: (%3f,%3f) "xyx_vy_vletprint_bodiesbodies=List.iter~f:print_bodybodies(***********************************************)(** Energy Calculations **********************)(***********************************************)letrecpairfold~flist~init=matchlistwith[]->init|hd::tl->letinit=List.fold_left~f:(funpartialel->fpartialhdel)~inittlinpairfold~ftl~init(* How to compute potential: sum of pair energy for all pairs (don't do pairs twice), where pair energy is G m_1 * m_2 / d Only works with grav_exp = 2.0 *)letpair_energyb1b2=letdist=mag(b1.pos<->b2.pos)in-.gconst#v*.b1.mass*.b2.mass/.dist(* returns potential energy *)letpenergybodies=pairfold~f:(funeb1b2->e+.pair_energyb1b2)~init:0.bodies(* returns kinetic energy *)letkenergybodies=List.fold_left~f:(funeb->e+.0.5*.b.mass*.magsqb.velocity)~init:0.bodiesletenergybodies=(* penergy bodies +. *)kenergybodies(***********************************************)(***********************************************)(***********************************************)letcenter_of_massbodies=matchbodieswith[]->(0.0,0.0)|_->letmpositions=List.map~f:(funbody->body.mass<*>body.pos)bodiesandmasses=List.map~f:(funbody->body.mass)bodiesin(summasses<|>vsummpositions)letcentral_velocitybodies=matchbodieswith[]->(0.0,0.0)|_->letmomenta=List.map~f:(funbody->body.mass<*>body.velocity)bodiesandmasses=List.map~f:(funbody->body.mass)bodiesin(summasses<|>vsummomenta)(***********************************************)letorbital_velocitybodies~posdir=(* first compute some global facts about the system *)letcom=center_of_massbodiesandmasses=List.map~f:(funbody->body.mass)bodiesandcv=central_velocitybodiesinletmass=summassesin(* now we compute the orbital speed *)letradius_vect=com-|posinletr=sqrt(dotradius_vectradius_vect)inletspeed=sqrt(gconst#v*.mass/.r)inletuvect=r/|radius_vectinletuvect=ifdirthenrotleftuvectelserotrightuvectin(speed*|uvect)+|cvletinduced_orbital_velocitybodies~posdir=ifList.lengthbodies=0then(0.0,0.0)elseletcv=central_velocitybodiesinletinduced_accel_list=List.map~f:(funbody->letdvect=body.pos-|posinletd=magdvectinletuvect=d/|dvectin(body.mass/.(d*.d))*|uvect)bodiesinletinduced_accel=List.fold_left~f:(+|)induced_accel_list~init:vzeroinlettotal_mass=sum(List.map~f:(funbody->body.mass)bodies)inletimplied_dist=sqrt(total_mass/.maginduced_accel)inletimplied_uvect=maginduced_accel/|induced_accelinletspeed=sqrt(gconst#v*.total_mass/.implied_dist)inletuvect=(ifdirthenrotleftimplied_uvectelserotrightimplied_uvect)in(speed*|uvect)+|cv(***********************************************)letsub_velocityvelbody={bodywithvelocity=body.velocity<->vel;}letzero_speed_bodiesselected_bodies=letvelocity=central_velocityselected_bodiesinstate.bodies<-List.map~f:(sub_velocityvelocity)state.bodiesletcenter_bodiesselected_bodies=letcenter=center_of_massselected_bodiesinstate.center#setcenter(***********************************************)letzero_speed()=zero_speed_bodiesstate.bodiesletcenter()=center_bodiesstate.bodiesletbodies_from_idsids=List.filter~f:(funbody->List.membody.id~set:ids)state.bodiesletzero_speeds_idsids=zero_speed_bodies(bodies_from_idsids)letcenter_idsids=center_bodies(bodies_from_idsids)(************************************************************************)(** Collision Detection ********************************************)(************************************************************************)lettouch~multb1b2=letmdist=maxb1.radiusb2.radiusindistsqb1.posb2.pos<mdist*.mdist*.mult*.multletjoin_bodiesb1b2={pos=center_of_mass[b1;b2];velocity=(b1.mass+.b2.mass)<|>((b1.mass<*>b1.velocity)<+>(b2.mass<*>b2.velocity));radius=((b1.radius**3.0)+.(b2.radius**3.0))**(1.0/.3.0);color=join_colorsb1.colorb1.massb2.colorb2.mass;mass=b1.mass+.b2.mass;id=Random.bits();i=None;}letfind_single_collision~multb1bodies=letrecloopb1bodiesexamined=matchbodieswith[]->b1::examined|b2::tl->iftouch~multb1b2thenloop(join_bodiesb1b2)tlexaminedelseloopb1tl(b2::examined)inloopb1bodies[](* look for a collision. If you find it, return a body list with those two bodies joined. Otherwise, return the original unchanged *)letrecfind_collisions~multbodies=matchbodieswith[]->[]|b1::tl->(find_single_collision~multb1(find_collisions~multtl))(************************************************************************)(** Simulation ******************************************************)(************************************************************************)letcomposefgx=f(gx)letcompose3fghx=f(g(hx))letidentx=xletrecapplynfx=matchnwith0->x|_->apply(n-1)f(fx)letsimulate?(bounce=false)i=letaction=compose(Fast_physics.act_all_on_all~bounce)(find_collisions~mult:(ifbouncethen0.5else1.0))instate.bodies<-applyiactionstate.bodies