ocaml-lablgl/balls.ml

   1 (* Source of this demo:
   2  * http://www.ffconsultancy.com/ocaml/index.html
   3  * Non-free, not for redistribution in Fedora.
   4  *)
   5
   6 open Printf
   7
   8 let pi = 4. *. atan 1.
   9
  10 module Array = struct
  11   include Array
  12
  13   let for_all p a =
  14     fold_left (fun b x -> p x && b) true a
  15 end
  16
  17 type vec2 = { x: float; y: float }
  18
  19 let vec2 x y = {x=x; y=y}
  20
  21 let zero = vec2 0. 0.
  22
  23 let of_list = function
  24     [x; y] -> vec2 x y
  25   | _ -> invalid_arg "Vec2.of_list"
  26
  27 let ( +| ) a b = {x = a.x +. b.x; y = a.y +. b.y}
  28 let ( -| ) a b = {x = a.x -. b.x; y = a.y -. b.y}
  29 let ( ~| ) a = {x = -. a.x; y = -. a.y}
  30 let ( *| ) s r = {x = s *. r.x; y = s *. r.y}
  31
  32 let normal a = { x = a.y; y = -. a.x }
  33 let dot a b = a.x *. b.x +. a.y *. b.y
  34 let length2 r = dot r r
  35 let length r = sqrt(length2 r)
  36 let unitise r = 1. /. length r *| r
  37
  38 (* Get the time since the program started. *)
  39 let time =
  40   let t = Unix.gettimeofday() in
  41   fun () -> Unix.gettimeofday() -. t
  42
  43 let width = ref 1 and height = ref 1
  44
  45 (* Reshape the viewport and store the width and height. *)
  46 let reshape ~w ~h =
  47   GlDraw.viewport ~x:0 ~y:0 ~w ~h;
  48   width := max 1 w;
  49   height := max 1 h
  50
  51 (* Pass a single vertex to the OpenGL. *)
  52 let vertex {x=x; y=y} =
  53   GlDraw.vertex ~x ~y ()
  54
  55 let translate r =
  56   GlMat.translate ~x:r.x ~y:r.y ()
  57
  58 let rotate angle =
  59   let angle = angle *. 180. /. pi in
  60   GlMat.rotate ~angle ~z:1. ()
  61
  62 let scale s =
  63   GlMat.scale ~x:s ~y:s ()
  64
  65 let protect f g =
  66   try
  67     f();
  68     g()
  69   with e ->
  70     g();
  71     raise e
  72
  73 let mat f =
  74   GlMat.push();
  75   protect f GlMat.pop
  76
  77 let render prim f =
  78   GlDraw.begins prim;
  79   protect f GlDraw.ends
  80
  81 let circle x = vec2 (sin x) (cos x)
  82
  83 (* Memoize OpenGL calls in a display list. *)
  84 let gl_memoize f =
  85   let dl = ref None in
  86   fun () -> match !dl with
  87   | Some dl -> GlList.call dl
  88   | None ->
  89       let dl' = GlList.create `compile in
  90       f ();
  91       GlList.ends ();
  92       dl := Some dl'
  93
  94 let rec iter f l u = if l<u then (f l; iter f (l+1) u)
  95
  96 (* Render a ball at the origin. *)
  97 let render_ball =
  98   gl_memoize
  99     (fun () ->
 100        let n = 36 in
 101        let aux i = circle (2. *. pi *. float i /. float n) in
 102        render `triangle_fan
 103          (fun () ->
 104             vertex zero;
 105             iter (fun i -> vertex (aux i)) 0 (n/4+1));
 106        render `triangle_fan
 107          (fun () ->
 108             vertex zero;
 109             iter (fun i -> vertex (aux i)) (n/2) (3*n/4+1));
 110        render `triangle_strip
 111          (fun () -> iter (fun i ->
 112                             vertex (0.9 *| aux i);
 113                             vertex (aux i)) 0 (n+1)))
 114
 115 type circle = { center: vec2; radius: float }
 116
 117 type ball = { circle: circle;
 118               velocity: vec2;
 119               angle: float;
 120               angular_velocity: float }
 121
 122 let make_ball ?(v=vec2 0.1 0.) ?(r=0.05) x y =
 123   { circle = { center = vec2 x y; radius = r };
 124     velocity = v;
 125     angle = 0.;
 126     angular_velocity = 0. }
 127
 128 type surface =
 129     Line of vec2 * vec2
 130   | Circle of vec2 * float
 131
 132 let surfaces = ref [Line (vec2 0. 0., vec2 0. 1.);
 133                     Line (vec2 0. 1., vec2 1. 1.);
 134                     Line (vec2 1. 1., vec2 1. 0.);
 135                     Line (vec2 1. 0.2, vec2 0. 0.)]
 136
 137 (* Split balls. *)
 138 let balls =
 139   ref [|{(make_ball 0.5 0.9) with velocity = vec2 0. (-1.)};
 140         {(make_ball 0.45 0.05) with velocity = vec2 0. 0.};
 141         {(make_ball 0.55 0.05) with velocity = vec2 0. 0.}|]
 142 let surfaces = ref [Line (vec2 0. 0., vec2 0. 1.);
 143                     Line (vec2 0. 1., vec2 1. 1.);
 144                     Line (vec2 1. 1., vec2 1. 0.);
 145                     Line (vec2 1. 0., vec2 0. 0.)]
 146 let g = vec2 0. 0.1
 147
 148 let n_balls = try min 100 (max 0 (int_of_string(Sys.argv.(1)))) with _ -> 3
 149
 150 (* Funnel *)
 151 let balls =
 152   let n = 10 in
 153   ref (Array.init n_balls
 154          (fun x -> make_ball ~v:(vec2 0. 0.) ~r:0.02
 155             (0.2 +. 0.03 *. float ((x/n) mod 2) +. 0.06 *. float (x mod n))
 156             (0.9 -. 0.04 *. float (x/n))))
 157 let surfaces =
 158   let n = 128 in
 159   let aux i =
 160     let r = circle (2. *. pi *. float (i + n/4) /. float n) in
 161     vec2 (0.5 +. 0.5 *. r.x) (1. +. r.y) in
 162   ref (Array.to_list (Array.init n (fun i -> Line (aux (i-1), aux i))))
 163 let surfaces = ref (Circle (vec2 0.37 0.3, 0.11) ::
 164                       Circle (vec2 0.63 0.3, 0.11) ::
 165                       !surfaces)
 166 let g = vec2 0. 0.5
 167
 168 let display_ball ball =
 169   mat (fun () ->
 170          translate ball.circle.center;
 171          rotate ball.angle;
 172          scale ball.circle.radius;
 173          render_ball())
 174
 175 let display_balls () =
 176   Array.iter display_ball !balls
 177
 178 let display_surface = function
 179     Line (v1, v2) ->
 180       GlDraw.begins `lines;
 181       List.iter vertex [v1; v2];
 182       GlDraw.ends ()
 183   | Circle (c, r) ->
 184       GlMat.push();
 185       mat (fun () ->
 186              translate c;
 187              scale r;
 188              let n = 360 in
 189              GlDraw.begins `line_loop;
 190              for i=0 to n-1 do
 191                vertex (circle (2. *. pi *. float i /. float n))
 192              done;
 193              GlDraw.ends ())
 194
 195 let display_surfaces =
 196   gl_memoize (fun () -> List.iter display_surface !surfaces)
 197
 198 let display () =
 199   GlClear.clear [ `color; `depth ];
 200   Gl.enable `depth_test;
 201   GlFunc.depth_func `lequal;
 202
 203   (* Initialise a orthogonal projection of the 2D scene. *)
 204   GlMat.mode `projection;
 205   GlMat.load_identity ();
 206   GlMat.ortho ~x:(0., 1.) ~y:(0., 1.) ~z:(0., 1.);
 207   GlMat.mode `modelview;
 208   GlMat.load_identity ();
 209
 210   GlDraw.color (1., 1., 1.) ~alpha:1.;
 211
 212   display_balls();
 213   display_surfaces();
 214
 215   Gl.finish ();
 216
 217   Unix.sleep 0;
 218
 219   Glut.swapBuffers ()
 220
 221 let clamp l u (x : float) =
 222   if x<l then l else if x>u then u else x
 223
 224 let circle_circle c1 r1 c2 r2 =
 225   let s = c2 -| c1 in
 226   c1 +| 0.5 *. (length s +. r1 -. r2) *| unitise s
 227
 228 let circle_circle c1 r1 c2 r2 =
 229   let s = c2 -| c1 in
 230   let l = length s in
 231   let s = (1. /. l) *| s in
 232   c1 +| 0.5 *. (l +. r1 -. r2) *| s
 233
 234 (* Find the point of impact between a circle and a surface. *)
 235 let point_of_impact circle = function
 236     Line (p1, p2) ->
 237       (* Find the closest approach of the finite line v1 -> v2 to the point
 238          r. *)
 239       let p21 = p2 -| p1 in
 240       let x = clamp 0. 1. (dot (circle.center -| p1) p21 /. length2 p21) in
 241       p1 +| x *| p21
 242   | Circle (c2, r2) -> circle_circle circle.center circle.radius c2 r2
 243
 244 let circle_circle_intersects c1 c2 =
 245   let p = circle_circle c1.center c1.radius c2.center c2.radius in
 246   length(c1.center -| p) < c1.radius
 247
 248 let rec update_ball dt (ball, impacts as accu) surface =
 249   let p = point_of_impact ball.circle surface in
 250
 251   let d = ball.circle.center -| p in
 252
 253   if length d > ball.circle.radius then accu else begin
 254     (* Local basis through center and perpendicular. *)
 255     let a = unitise d in
 256     let b = normal a in
 257
 258     (* Resolve the velocity at the point of impact into components through the
 259        center and perpendicular. *)
 260     let dva = dot ball.velocity a in
 261     let dvb = dot ball.velocity b in
 262     if dva >= 0. then accu else
 263       { ball with
 264           velocity = ball.velocity +| 2. *. abs_float dva *| a -|
 265               0. *. ball.angular_velocity *. ball.circle.radius *| b;
 266           angular_velocity = -. dvb /. ball.circle.radius }, p::impacts
 267   end
 268
 269 (* Bisect the time slice until there is only one impact. *)
 270 let rec update t t' balls =
 271   let dt = t' -. t in
 272
 273   (* Ball-ball impacts. *)
 274   let balls', impacts =
 275     let balls = Array.copy balls in
 276     let impacts = Array.map (fun _ -> []) balls in
 277     let n = Array.length balls in
 278
 279     for i=0 to n-1 do
 280       (* Ball-surface impacts *)
 281       (fun (b, is) ->
 282          balls.(i) <- b;
 283          impacts.(i) <- is)
 284         (List.fold_left (update_ball dt) (balls.(i), impacts.(i)) !surfaces);
 285
 286       (* Ball-ball impacts. *)
 287       if i>0 then
 288         let b1 = balls.(i) in
 289         for j=0 to i-1 do
 290           let b2 = balls.(j) in
 291           let c1 = b1.circle and c2 = b2.circle in
 292           let p = circle_circle c1.center c1.radius c2.center c2.radius in
 293           if length(c1.center -| p) < c1.radius then begin
 294             let u1 = b1.velocity and u2 = b2.velocity in
 295             let da1 = b1.angular_velocity and da2 = b2.angular_velocity in
 296             let r1 = c1.radius and r2 = c2.radius in
 297
 298             (* Find the velocity difference to the center-of-momentum frame. *)
 299             let com = 0.5 *| (u1 +| u2) in
 300
 301             (* Move to COM frame. *)
 302             let u1 = u1 -| com and u2 = u2 -| com in
 303             let u = unitise (c2.center -| c1.center) in
 304             let v = normal u in
 305             let impulse = u2 -| u1 +| (da1 *. r1 -. da2 *. r2) *| v in
 306             let impulse_u = dot u impulse *| u in
 307             let impulse_v = dot v impulse in
 308             let v1 = u1 +| impulse_u and v2 = u2 -| impulse_u in
 309             let da1' = da1 +. impulse_v *. r1 in
 310             let da2' = da2 -. impulse_v *. r2 in
 311
 312             (* Move from COM frame. *)
 313             let v1 = v1 +| com and v2 = v2 +| com in
 314
 315             balls.(i) <- { balls.(i) with
 316                              velocity = v1;
 317                              angular_velocity = da1' };
 318             balls.(j) <- { balls.(j) with
 319                              velocity = v2;
 320                              angular_velocity = da2' };
 321             impacts.(i) <- p :: impacts.(i);
 322             impacts.(j) <- p :: impacts.(j);
 323           end
 324         done;
 325     done;
 326
 327     for i=0 to n-1 do
 328       let b = balls.(i) in
 329       let r = b.circle.center and dr = b.velocity in
 330       let a = b.angle and da = b.angular_velocity in
 331       let rec aux dr = function
 332           [] ->
 333              { balls.(i) with
 334                  circle = { balls.(i).circle with center = r +| dt *| dr };
 335                  angle = mod_float (a +. dt *. da) (2. *. pi);
 336                  velocity = dr -| (t' -. t) *| g }
 337         | p::t ->
 338             (* Make sure the velocity is not pointing towards the impact
 339                point. *)
 340             let d = unitise(r -| p) in
 341             aux (dr -| min 0. (dot dr d) *| d) t in
 342       balls.(i) <- aux dr impacts.(i);
 343       match impacts.(i) with
 344         [] | [_] -> ()
 345       | _ ->
 346           balls.(i) <- { balls.(i) with angular_velocity = 0. }
 347     done;
 348     balls, impacts in
 349
 350   (* Bisect if there was at least one impact and the time slice was large
 351      enough. *)
 352   let aux = function [] | [_] -> true | _ -> false in
 353   if dt<0.01 && (Array.for_all aux impacts || dt < 1e-3) then balls' else
 354     let t2 = (t +. t') *. 0.5 in
 355     update t2 t' (update t t2 balls)
 356
 357 let old_time = ref 0.
 358
 359 let idle () =
 360   let time' = time() in
 361   balls := update !old_time time' !balls;
 362   old_time := time';
 363   Glut.postRedisplay ()
 364
 365 let () =
 366   ignore (Glut.init Sys.argv);
 367   Glut.initDisplayMode ~alpha:true ~double_buffer:true ~depth:true ();
 368   Glut.initWindowSize ~w:1024 ~h:1024;
 369   ignore (Glut.createWindow ~title:"Bouncing balls");
 370   GlClear.color (0., 0., 0.) ~alpha:0.;
 371
 372   Gl.enable `blend;
 373   List.iter Gl.enable [`line_smooth; `polygon_smooth];
 374   List.iter (fun x -> GlMisc.hint x `nicest) [`line_smooth; `polygon_smooth];
 375   GlDraw.line_width 4.;
 376   GlFunc.blend_func ~src:`src_alpha ~dst:`one;
 377
 378   Glut.reshapeFunc ~cb:reshape;
 379   Glut.displayFunc ~cb:display;
 380   Glut.idleFunc ~cb:(Some idle);
 381   Glut.keyboardFunc ~cb:(fun ~key ~x ~y -> if key=27 then exit 0);
 382   Glut.mainLoop ()