open Graphics;; class point (a: float) (b: float) (m: float) = object val mutable x = a val mutable y = b val mutable mass = m val mutable dt = 1. method sdt t = dt <- t method x = x method y = y method mass = mass method smass m = mass <- m method xy = (x, y) method moveto a b = x <- x +. (a *. dt); y <- y +. (b *. dt) method teleport a b = x <- a; y <- b val mutable fx = 0. val mutable fy = 0. method fx = fx method fy = fy method fxy = (fx, fy) method addforce a b = fx <- fx +. (a /. mass); fy <- fy +. (b /. mass) method addfriction q = fx <- fx /. q; fy <- fy /. q method move () = x <- x +. fx; y <- y +. fy end;; let point_num = 100;; let point_size = 3.;; let start_zone = 1;; let gconst = 1.;; let force_increase = 20.;; Random.self_init ();; let ( >> ) x f = f x;; let sqr x = x *. x;; let vector_length x y = sqr x +. sqr y >> sqrt;; let cut_vector x y l z = if z = 0. then (0., 0.) else (x *. l /. z, y *. l /. z);; let gforce m1 m2 r = if r <= 0. then 0. else (gconst *. m1 *. m2) /. (sqr r);; open_graph "";; let points = Array.init point_num (fun _ -> new point (size_x () / start_zone >> Random.int >> ((+) ((start_zone / 2) * (size_x () / start_zone))) >> float_of_int) (size_y () / start_zone >> Random.int >> ((+) ((start_zone / 2) * (size_y () / start_zone))) >> float_of_int) (Random.float 3.5 +. 0.5));; let rec full_mass i = if i < point_num then points.(i)#mass +. (full_mass (i + 1)) else 0.;; let full_mass () = full_mass 0;; let rec mass_centre_x i = if i < point_num then (points.(i)#x *. points.(i)#mass) +. (mass_centre_x (i + 1)) else 0.;; let rec mass_centre_y i = if i < point_num then (points.(i)#y *. points.(i)#mass) +. (mass_centre_y (i + 1)) else 0.;; let mass_centre () = ((*int_of_float *)((mass_centre_x 0) /. full_mass ()), (*int_of_float *)((mass_centre_y 0) /. full_mass ()));; (*points.(0)#smass 50.;; points.(0)#teleport ((mass_centre_x 0) /. full_mass ()) ((mass_centre_y 0) /. full_mass ());;*) let time = ref (Unix.time ());; while (not (key_pressed ())) do auto_synchronize false; clear_graph (); for i = 0 to point_num - 1 do for j = i + 1 to point_num - 1 do let (offx, offy) = (points.(j)#x -. points.(i)#x, points.(j)#y -. points.(i)#y) in let dist = vector_length offx offy in let nforce = gforce (points.(i)#mass) (points.(j)#mass) dist in let (nfx, nfy) = cut_vector offx offy nforce dist in points.(i)#addforce nfx nfy; points.(j)#addforce (-.nfx) (-.nfy); (* if dist < 20. then Printf.printf "|> %f %f %f <| " nfx nfy dist;*) done; let (x, y) = (points.(i)#x >> int_of_float, points.(i)#y >> int_of_float) in draw_circle x y (int_of_float (point_size *. points.(i)#mass)); (* moveto x y; draw_string (Printf.sprintf "(%F, %F)" points.(i)#x points.(i)#y);*) moveto x y; let (x, y) = (points.(i)#fx, points.(i)#fy) in let len = vector_length x y in let (x, y) = cut_vector x y (len *. force_increase) len in rlineto (int_of_float x) (int_of_float y); points.(i)#sdt (Unix.time () -. !time); points.(i)#move (); (* if ((int_of_float points.(i)#x) > (size_x ())) then points.(i)#teleport 0. points.(i)#y else if (points.(i)#x < 0.) then points.(i)#teleport (size_x () >> float_of_int) points.(i)#y; if ((int_of_float points.(i)#y) > (size_y ())) then points.(i)#teleport points.(i)#x 0. else if (points.(i)#y < 0.) then points.(i)#teleport points.(i)#x (size_y () >> float_of_int);*) done; time := Unix.time (); let (x, y) = mass_centre () in fill_circle (int_of_float x) (int_of_float y) 5; moveto 0 0; draw_string (Printf.sprintf "(%F, %F)" x y); synchronize (); done;; close_graph ();;