-- start 20/03/2014 14:38 - modified from nbody1 -- nearbyPairs :: (Float, Spatial v) -> Map (Vec3,Vec3) (v,v) -- nearbyJoin (Float, Spatial v, Spatial w) -> Map (Vec3,Vec3) (v,w) -- spatialFromMap :: Map Vec3 v -> Spatial v -- mapFromSpatial :: Spatial v -> Map Vec3 v -- helper functions let G = 0.00000000000667384 in -- gravitational constant let vec0 = (0.0, 0.0, 0.0) in let app3 = \(f,(x,y,z)) -> (f x, f y, f z) in -- apply unary op elementwise to vec3 let app23 = \(f,(x,y,z),(x',y',z')) -> (f (x,x'), f (y,y'), f (z,z')) in -- apply bin op to two vec3s let red3 = \(f,(x,y,z)) -> f (f (x,y), z) in -- reduce vec3 using f let distsq = \dd -> red3 (addf, app3 (sqr dd)) in -- distance squared let fmag = \(m,m',d2) -> divf (mulf (G, mulf (m,m')), d2) in -- magnitude of force -- calculate force between two objects let force = \((p,m),(p',m')) -> -- gravitational force between two objects let dd = app23 (subf,p',p) in -- diff position let d2 = distsq (dx,dy,dz) in -- distance squared let a = divf (fmag (m,m',d2), sqrt d2) in -- fmag / dist app3 (\k -> mulf (a,k), dd) -- 3D components of force app3 (\k -> mulf (a,k), dd) -- 3D components of force -- move an object according to it's force let move = \(p,v,f,m,dt) -> -- p position, v velocity, f force, m mass, dt timestep let (dt_m, dt_2) = (divf (dt,m), divf (dt, 2)) in let dv = app3 (\a -> mulf (a, dt_m), f) in -- velocity change let dp = app3 (\a -> mulf (a, dt_2), app23 (addf, v, dv)) in -- position change let v' = app23 (addf, v, dv) in -- new velocity let p' = app23 (addf, p, dp) in -- new position (p',v') -- update all objects for one timestep -- takes array of masses, array of positions and velocities, and timestep let sim = \(ml,pl,dt) :: (Arr Int Float, Arr Int (Vec3,Vec3), Float) -> let pl1 = allPairsArr (nullF, pl) in -- all combinations of particles, Map (id1,id2) (p1,p2) let pl2 = eqJoinArr (id, snd.fst, ol, eqJoinArr (id, fst, ol, pl1)) in -- join to get masses, Map (id2,(id1,(id1,id2))) (m2,(m1,(p1,p2))) let pl3 = mapArrInv (\(i1,(i2,_)) -> (i1,i2), \(i1,i2) -> (i1,(i2,(i1,i2))), -- calculate forces between pairs \(_,(m2,(m1,((p1,v1),(p2,v2))))) -> let f = force ((p1,m1),(p2,m2)) in (f, subf (0.0, f)), pl2) in let fl = mapArrInv (id, id, fst) in -- forces for left hand particle let fr = mapArrInv (swap, swap, snd) in -- forces for right hand particle let fsum = groupReduceArr (fst, snd, addf, disUnionArr (fl, fr)) in -- total forces for particles Map Int Vec3 let data = eqJoinArr (id, fst, ml, eqJoinArr (id, id, pl, fsum)) in -- join mass, position, velocity, and total forces mapArrInv (fst, \i -> (i,(i,i)), \(_,(m,((p,v),f))) -> move (p,v,f,m,dt), data) -- move particles using forces -- generate random particles let N = 1000000 in let ml = mapArrInv (id, id, \_ -> randf(), rangeArr (1,N,1)) in let pl0 = mapArrInv (id, id, \_ -> (app3 (\_ -> randf (), vec0), app3 (\_ -> randf (), vec0)), rangeArr (1,N,1)) in -- simulation loop (while t < 1000, steps of 1.0) let (maxT, stepT) = (1000.0, 1.0) in while (\(pl,t) -> ((sim (ml,pl,stepT), addf (t, stepT)), ltf (t,maxT)), (pl0,0.0))