--- /dev/null
+(* Source of this demo:
+ * http://www.ffconsultancy.com/ocaml/index.html
+ * Non-free, not for redistribution in Fedora.
+ *)
+
+open Printf
+
+let pi = 4. *. atan 1.
+
+module Array = struct
+ include Array
+
+ let for_all p a =
+ fold_left (fun b x -> p x && b) true a
+end
+
+type vec2 = { x: float; y: float }
+
+let vec2 x y = {x=x; y=y}
+
+let zero = vec2 0. 0.
+
+let of_list = function
+ [x; y] -> vec2 x y
+ | _ -> invalid_arg "Vec2.of_list"
+
+let ( +| ) a b = {x = a.x +. b.x; y = a.y +. b.y}
+let ( -| ) a b = {x = a.x -. b.x; y = a.y -. b.y}
+let ( ~| ) a = {x = -. a.x; y = -. a.y}
+let ( *| ) s r = {x = s *. r.x; y = s *. r.y}
+
+let normal a = { x = a.y; y = -. a.x }
+let dot a b = a.x *. b.x +. a.y *. b.y
+let length2 r = dot r r
+let length r = sqrt(length2 r)
+let unitise r = 1. /. length r *| r
+
+(* Get the time since the program started. *)
+let time =
+ let t = Unix.gettimeofday() in
+ fun () -> Unix.gettimeofday() -. t
+
+let width = ref 1 and height = ref 1
+
+(* Reshape the viewport and store the width and height. *)
+let reshape ~w ~h =
+ GlDraw.viewport ~x:0 ~y:0 ~w ~h;
+ width := max 1 w;
+ height := max 1 h
+
+(* Pass a single vertex to the OpenGL. *)
+let vertex {x=x; y=y} =
+ GlDraw.vertex ~x ~y ()
+
+let translate r =
+ GlMat.translate ~x:r.x ~y:r.y ()
+
+let rotate angle =
+ let angle = angle *. 180. /. pi in
+ GlMat.rotate ~angle ~z:1. ()
+
+let scale s =
+ GlMat.scale ~x:s ~y:s ()
+
+let protect f g =
+ try
+ f();
+ g()
+ with e ->
+ g();
+ raise e
+
+let mat f =
+ GlMat.push();
+ protect f GlMat.pop
+
+let render prim f =
+ GlDraw.begins prim;
+ protect f GlDraw.ends
+
+let circle x = vec2 (sin x) (cos x)
+
+(* Memoize OpenGL calls in a display list. *)
+let gl_memoize f =
+ let dl = ref None in
+ fun () -> match !dl with
+ | Some dl -> GlList.call dl
+ | None ->
+ let dl' = GlList.create `compile in
+ f ();
+ GlList.ends ();
+ dl := Some dl'
+
+let rec iter f l u = if l<u then (f l; iter f (l+1) u)
+
+(* Render a ball at the origin. *)
+let render_ball =
+ gl_memoize
+ (fun () ->
+ let n = 36 in
+ let aux i = circle (2. *. pi *. float i /. float n) in
+ render `triangle_fan
+ (fun () ->
+ vertex zero;
+ iter (fun i -> vertex (aux i)) 0 (n/4+1));
+ render `triangle_fan
+ (fun () ->
+ vertex zero;
+ iter (fun i -> vertex (aux i)) (n/2) (3*n/4+1));
+ render `triangle_strip
+ (fun () -> iter (fun i ->
+ vertex (0.9 *| aux i);
+ vertex (aux i)) 0 (n+1)))
+
+type circle = { center: vec2; radius: float }
+
+type ball = { circle: circle;
+ velocity: vec2;
+ angle: float;
+ angular_velocity: float }
+
+let make_ball ?(v=vec2 0.1 0.) ?(r=0.05) x y =
+ { circle = { center = vec2 x y; radius = r };
+ velocity = v;
+ angle = 0.;
+ angular_velocity = 0. }
+
+type surface =
+ Line of vec2 * vec2
+ | Circle of vec2 * float
+
+let surfaces = ref [Line (vec2 0. 0., vec2 0. 1.);
+ Line (vec2 0. 1., vec2 1. 1.);
+ Line (vec2 1. 1., vec2 1. 0.);
+ Line (vec2 1. 0.2, vec2 0. 0.)]
+
+(* Split balls. *)
+let balls =
+ ref [|{(make_ball 0.5 0.9) with velocity = vec2 0. (-1.)};
+ {(make_ball 0.45 0.05) with velocity = vec2 0. 0.};
+ {(make_ball 0.55 0.05) with velocity = vec2 0. 0.}|]
+let surfaces = ref [Line (vec2 0. 0., vec2 0. 1.);
+ Line (vec2 0. 1., vec2 1. 1.);
+ Line (vec2 1. 1., vec2 1. 0.);
+ Line (vec2 1. 0., vec2 0. 0.)]
+let g = vec2 0. 0.1
+
+let n_balls = try min 100 (max 0 (int_of_string(Sys.argv.(1)))) with _ -> 3
+
+(* Funnel *)
+let balls =
+ let n = 10 in
+ ref (Array.init n_balls
+ (fun x -> make_ball ~v:(vec2 0. 0.) ~r:0.02
+ (0.2 +. 0.03 *. float ((x/n) mod 2) +. 0.06 *. float (x mod n))
+ (0.9 -. 0.04 *. float (x/n))))
+let surfaces =
+ let n = 128 in
+ let aux i =
+ let r = circle (2. *. pi *. float (i + n/4) /. float n) in
+ vec2 (0.5 +. 0.5 *. r.x) (1. +. r.y) in
+ ref (Array.to_list (Array.init n (fun i -> Line (aux (i-1), aux i))))
+let surfaces = ref (Circle (vec2 0.37 0.3, 0.11) ::
+ Circle (vec2 0.63 0.3, 0.11) ::
+ !surfaces)
+let g = vec2 0. 0.5
+
+let display_ball ball =
+ mat (fun () ->
+ translate ball.circle.center;
+ rotate ball.angle;
+ scale ball.circle.radius;
+ render_ball())
+
+let display_balls () =
+ Array.iter display_ball !balls
+
+let display_surface = function
+ Line (v1, v2) ->
+ GlDraw.begins `lines;
+ List.iter vertex [v1; v2];
+ GlDraw.ends ()
+ | Circle (c, r) ->
+ GlMat.push();
+ mat (fun () ->
+ translate c;
+ scale r;
+ let n = 360 in
+ GlDraw.begins `line_loop;
+ for i=0 to n-1 do
+ vertex (circle (2. *. pi *. float i /. float n))
+ done;
+ GlDraw.ends ())
+
+let display_surfaces =
+ gl_memoize (fun () -> List.iter display_surface !surfaces)
+
+let display () =
+ GlClear.clear [ `color; `depth ];
+ Gl.enable `depth_test;
+ GlFunc.depth_func `lequal;
+
+ (* Initialise a orthogonal projection of the 2D scene. *)
+ GlMat.mode `projection;
+ GlMat.load_identity ();
+ GlMat.ortho ~x:(0., 1.) ~y:(0., 1.) ~z:(0., 1.);
+ GlMat.mode `modelview;
+ GlMat.load_identity ();
+
+ GlDraw.color (1., 1., 1.) ~alpha:1.;
+
+ display_balls();
+ display_surfaces();
+
+ Gl.finish ();
+
+ Unix.sleep 0;
+
+ Glut.swapBuffers ()
+
+let clamp l u (x : float) =
+ if x<l then l else if x>u then u else x
+
+let circle_circle c1 r1 c2 r2 =
+ let s = c2 -| c1 in
+ c1 +| 0.5 *. (length s +. r1 -. r2) *| unitise s
+
+let circle_circle c1 r1 c2 r2 =
+ let s = c2 -| c1 in
+ let l = length s in
+ let s = (1. /. l) *| s in
+ c1 +| 0.5 *. (l +. r1 -. r2) *| s
+
+(* Find the point of impact between a circle and a surface. *)
+let point_of_impact circle = function
+ Line (p1, p2) ->
+ (* Find the closest approach of the finite line v1 -> v2 to the point
+ r. *)
+ let p21 = p2 -| p1 in
+ let x = clamp 0. 1. (dot (circle.center -| p1) p21 /. length2 p21) in
+ p1 +| x *| p21
+ | Circle (c2, r2) -> circle_circle circle.center circle.radius c2 r2
+
+let circle_circle_intersects c1 c2 =
+ let p = circle_circle c1.center c1.radius c2.center c2.radius in
+ length(c1.center -| p) < c1.radius
+
+let rec update_ball dt (ball, impacts as accu) surface =
+ let p = point_of_impact ball.circle surface in
+
+ let d = ball.circle.center -| p in
+
+ if length d > ball.circle.radius then accu else begin
+ (* Local basis through center and perpendicular. *)
+ let a = unitise d in
+ let b = normal a in
+
+ (* Resolve the velocity at the point of impact into components through the
+ center and perpendicular. *)
+ let dva = dot ball.velocity a in
+ let dvb = dot ball.velocity b in
+ if dva >= 0. then accu else
+ { ball with
+ velocity = ball.velocity +| 2. *. abs_float dva *| a -|
+ 0. *. ball.angular_velocity *. ball.circle.radius *| b;
+ angular_velocity = -. dvb /. ball.circle.radius }, p::impacts
+ end
+
+(* Bisect the time slice until there is only one impact. *)
+let rec update t t' balls =
+ let dt = t' -. t in
+
+ (* Ball-ball impacts. *)
+ let balls', impacts =
+ let balls = Array.copy balls in
+ let impacts = Array.map (fun _ -> []) balls in
+ let n = Array.length balls in
+
+ for i=0 to n-1 do
+ (* Ball-surface impacts *)
+ (fun (b, is) ->
+ balls.(i) <- b;
+ impacts.(i) <- is)
+ (List.fold_left (update_ball dt) (balls.(i), impacts.(i)) !surfaces);
+
+ (* Ball-ball impacts. *)
+ if i>0 then
+ let b1 = balls.(i) in
+ for j=0 to i-1 do
+ let b2 = balls.(j) in
+ let c1 = b1.circle and c2 = b2.circle in
+ let p = circle_circle c1.center c1.radius c2.center c2.radius in
+ if length(c1.center -| p) < c1.radius then begin
+ let u1 = b1.velocity and u2 = b2.velocity in
+ let da1 = b1.angular_velocity and da2 = b2.angular_velocity in
+ let r1 = c1.radius and r2 = c2.radius in
+
+ (* Find the velocity difference to the center-of-momentum frame. *)
+ let com = 0.5 *| (u1 +| u2) in
+
+ (* Move to COM frame. *)
+ let u1 = u1 -| com and u2 = u2 -| com in
+ let u = unitise (c2.center -| c1.center) in
+ let v = normal u in
+ let impulse = u2 -| u1 +| (da1 *. r1 -. da2 *. r2) *| v in
+ let impulse_u = dot u impulse *| u in
+ let impulse_v = dot v impulse in
+ let v1 = u1 +| impulse_u and v2 = u2 -| impulse_u in
+ let da1' = da1 +. impulse_v *. r1 in
+ let da2' = da2 -. impulse_v *. r2 in
+
+ (* Move from COM frame. *)
+ let v1 = v1 +| com and v2 = v2 +| com in
+
+ balls.(i) <- { balls.(i) with
+ velocity = v1;
+ angular_velocity = da1' };
+ balls.(j) <- { balls.(j) with
+ velocity = v2;
+ angular_velocity = da2' };
+ impacts.(i) <- p :: impacts.(i);
+ impacts.(j) <- p :: impacts.(j);
+ end
+ done;
+ done;
+
+ for i=0 to n-1 do
+ let b = balls.(i) in
+ let r = b.circle.center and dr = b.velocity in
+ let a = b.angle and da = b.angular_velocity in
+ let rec aux dr = function
+ [] ->
+ { balls.(i) with
+ circle = { balls.(i).circle with center = r +| dt *| dr };
+ angle = mod_float (a +. dt *. da) (2. *. pi);
+ velocity = dr -| (t' -. t) *| g }
+ | p::t ->
+ (* Make sure the velocity is not pointing towards the impact
+ point. *)
+ let d = unitise(r -| p) in
+ aux (dr -| min 0. (dot dr d) *| d) t in
+ balls.(i) <- aux dr impacts.(i);
+ match impacts.(i) with
+ [] | [_] -> ()
+ | _ ->
+ balls.(i) <- { balls.(i) with angular_velocity = 0. }
+ done;
+ balls, impacts in
+
+ (* Bisect if there was at least one impact and the time slice was large
+ enough. *)
+ let aux = function [] | [_] -> true | _ -> false in
+ if dt<0.01 && (Array.for_all aux impacts || dt < 1e-3) then balls' else
+ let t2 = (t +. t') *. 0.5 in
+ update t2 t' (update t t2 balls)
+
+let old_time = ref 0.
+
+let idle () =
+ let time' = time() in
+ balls := update !old_time time' !balls;
+ old_time := time';
+ Glut.postRedisplay ()
+
+let () =
+ ignore (Glut.init Sys.argv);
+ Glut.initDisplayMode ~alpha:true ~double_buffer:true ~depth:true ();
+ Glut.initWindowSize ~w:1024 ~h:1024;
+ ignore (Glut.createWindow ~title:"Bouncing balls");
+ GlClear.color (0., 0., 0.) ~alpha:0.;
+
+ Gl.enable `blend;
+ List.iter Gl.enable [`line_smooth; `polygon_smooth];
+ List.iter (fun x -> GlMisc.hint x `nicest) [`line_smooth; `polygon_smooth];
+ GlDraw.line_width 4.;
+ GlFunc.blend_func ~src:`src_alpha ~dst:`one;
+
+ Glut.reshapeFunc ~cb:reshape;
+ Glut.displayFunc ~cb:display;
+ Glut.idleFunc ~cb:(Some idle);
+ Glut.keyboardFunc ~cb:(fun ~key ~x ~y -> if key=27 then exit 0);
+ Glut.mainLoop ()
git.annexia.org / fedora-mingw.git / commitdiff