(* X-MP-i.ml : code Pascal
 *  pour me compiler et m'exécuter (sous Linux) :
 *  
 *  camlc -o essai X-MP-i.ml
 *  ./essai >essai.ps
 *  gs essai.ps
 *)

#open "printf";;


                 (* +---------------------------------------+
                    |  Déclarations figurant dans l'énoncé  |
                    +---------------------------------------+ *)


type vecteur = {x: float; y: float};;
let vecteur_nul = {x=0.0; y=0.0};;

type corps = {
  mass: float;
  mutable pos: vecteur;
  mutable vel: vecteur;
  mutable acc: vecteur;
  couleur: string;         (* ajouté pour l'animation *)
  mutable oldpos: vecteur; (* idem *)
  mutable oldvel: vecteur  (* idem *)
};;

(* initialisation de l'univers *)
let univers = [|
  {mass    = 3000.0;
   pos     = {x=0.0; y=0.0};
   vel     = {x=5.0; y=0.0};
   couleur = "1 0 0 s";      (* rouge *)
   acc     = vecteur_nul;
   oldpos  = vecteur_nul;
   oldvel  = vecteur_nul};

  {mass    = 1500.0;
   pos     = {x=100.0; y=110.0};
   vel     = {x=5.0; y=0.0};
   couleur = "0 0 1 s";      (* bleu *)
   acc     = vecteur_nul;
   oldpos  = vecteur_nul;
   oldvel  = vecteur_nul};

  {mass    = 1500.0;
   pos     = {x=25.0; y=150.0};
   vel     = {x= -9.0; y=0.0};
   couleur = "0 1 0 s";      (* vert *)
   acc     = vecteur_nul;
   oldpos  = vecteur_nul;
   oldvel  = vecteur_nul}
|];;

let N = vect_length univers  (* nb effectif de corps *)
and taille_univers = 500.0   (* dimension maximum    *)
and theta = 0.1         (* paramètre d'approximation *)
;;

type arbre = Noeud of cellule | Feuille of corps | Vide
and  cellule = {
  mutable cm_mass: float;
  mutable cm_pos:  vecteur;
  filles: arbre vect
};;
let arbre_univers = ref(Vide);;

                            (* +-----------------+
                               |  Question 2abc  |
                               +-----------------+ *)

let add v1 v2 = {x=v1.x+.v2.x; y=v1.y+.v2.y};;
let sub v1 v2 = {x=v1.x-.v2.x; y=v1.y-.v2.y};;
let scal m u  = {x=m*.u.x; y=m*.u.y};;
let carre u   = u.x*.u.x +. u.y*.u.y;;

                            (* +----------------+
                               |  Question 3ab  |
                               +----------------+ *)

let indice_fille p_c p =
  (if p_c.y > p.y then 1 else 0) + 2*(if p_c.x > p.x then 1 else 0);;

let position_fille p_c taille i =
  let t = taille/.4.0 in
  match i with
  | 0 -> add p_c {x =   t; y =   t}
  | 1 -> add p_c {x =   t; y = -.t}
  | 2 -> add p_c {x = -.t; y =   t}
  | _ -> add p_c {x = -.t; y = -.t}
;;

                             (* +---------------+
                                |  Question 4a  |
                                +---------------+ *)

(* insère un corps dans un quadtree *)
let rec insere_corps corps arbre p_c taille = match arbre with
| Vide -> Feuille(corps)
| Noeud(cell) ->
   insere_dans_filles corps cell.filles p_c taille;
   arbre
| Feuille(c) ->
   let f = make_vect 4 Vide in
   insere_dans_filles c     f p_c taille;
   insere_dans_filles corps f p_c taille;
   Noeud{cm_mass=0.0; cm_pos=vecteur_nul; filles=f}

(* insère un corps dans un tableau de filles *)
and insere_dans_filles corps filles p_c taille =
  let i = indice_fille   p_c corps.pos in
  let p = position_fille p_c taille i  in
  filles.(i) <- insere_corps corps filles.(i) p (taille/.2.0)
;;

                             (* +--------------+
                                |  Question 5  |
                                +--------------+ *)

(* ajuste les champs$ cm_mass $et$ cm_pos $et retourne *)
(* le couple (masse,somme pondérée des positions)  *)
let rec mass_pos(arbre) = match arbre with
| Vide        -> (0.0, vecteur_nul)
| Feuille(c)  -> (c.mass, scal c.mass c.pos)
| Noeud(cell) ->
    let m = ref(0.0) and p = ref(vecteur_nul) in
    for i=0 to 3 do
      let (mi,pi) = mass_pos(cell.filles.(i)) in
      m := !m +. mi;
      p := add !p pi
    done;
    cell.cm_mass <- !m;
    cell.cm_pos  <- scal (1.0/. !m) !p;
    (!m,!p)
;;

let barycentres(arbre) = let _ = mass_pos(arbre) in ();;

                             (* +---------------+
                                |  Question 6a  |
                                +---------------+ *)

(* fonction auxilliaire : distance numérique et vectorielle *)
let distance a b = let d = sub b a in (sqrt(carre(d)),d);;

(* fonction auxilliaire : accélération de la masse$ m $placée en$ b $sur$ a $*)
let accel a b m = let (r,d) = distance a b in scal (m/.r/.r/.r) d;;

(* accélération de l'univers sur un corps *)
let rec grav_approx pos arbre taille = match arbre with
| Vide        -> vecteur_nul
| Feuille(c)  -> if c.pos = pos then vecteur_nul else accel pos c.pos c.mass
| Noeud(cell) ->
    let (r,_) = distance pos cell.cm_pos in if taille < r*.theta
    then begin
      printf "%%%% approx %f %f %f %f %f\n" pos.x pos.y cell.cm_pos.x cell.cm_pos.y taille;
      accel pos cell.cm_pos cell.cm_mass
    end
    else begin
      let acc = ref(vecteur_nul) and t = taille/.2.0 in
      for i=0 to 3 do acc := add !acc (grav_approx pos cell.filles.(i) t) done;
      !acc
    end
;;

                   (* +-----------------------------------+
                      |  Assemblage fourni dans l'énoncé  |
                      +-----------------------------------+ *)

let construire_arbre() =
  arbre_univers := Vide;
  for i=0 to N-1 do
    arbre_univers := insere_corps univers.(i) !arbre_univers vecteur_nul taille_univers
 done
;;

let ajuste_acc_approx() =
  for i=0 to N-1 do
    univers.(i).acc <- grav_approx univers.(i).pos !arbre_univers taille_univers
  done
;;

                          (* +---------------------+
                             |  Dessine l'univers  |
                             +---------------------+ *)

let rayon = 3.0;; (* rayon des points *)
let rec dessine a p taille = match a with
| Vide -> ()
| Noeud(cell) ->
   let t = taille/.2.0 in
   printf "%f %f %f %f l\n" p.x (p.y-.t) p.x (p.y+.t);
   printf "%f %f %f %f l\n" (p.x-.t) p.y (p.x+.t) p.y;
   for i=0 to 3 do dessine cell.filles.(i) (position_fille p taille i) t done
| Feuille(c) ->
  printf "gsave newpath %s %f %f %f r grestore\n" c.couleur rayon c.pos.x c.pos.y
;;

let dessine_univers() =
  printf "%%%% photo\n0 0 0 s\n%f dup -2 div dup moveto " taille_univers;
  printf "dup 0 rlineto dup 0 exch rlineto neg 0 rlineto closepath stroke\n";
  dessine !arbre_univers vecteur_nul taille_univers;
  printf "%%%% endphoto\n"
;;

                         (* +-----------------------+
                            |  Programme principal  |
                            +-----------------------+ *)


let tmax = 50.0;;  (* temps maximum *)
let t0   = 27.3;;  (* dessine l'univers à cet instant *)
let h    = 0.1;;   (* intervalle élémentaire de temps *)
let h2   = h*.h/.2.0;;

let t = ref 0.0;;
exception fin;;

begin
   (* en-tête Postscript *)
   printf "%%!PS-Adobe-2.0\n";
   printf "%%%%BoundingBox: 0 0 %f %f\n" taille_univers taille_univers;
   printf "/l {4 2 roll moveto lineto stroke} def\n";
   printf "/s {setrgbcolor} def\n";
   printf "/r {translate dup 0 moveto 0 0 3 2 roll 0 360 arc fill} def\n";
   printf "%f 2 div dup translate\n" taille_univers;

   t := 0.0;
   begin try while !t < tmax do

   	 printf "%%%% t=%f\n" !t;

   	 (* division de l'univers en cellules *)
   	 construire_arbre();
   	 barycentres(!arbre_univers);
   	 if abs_float(!t-.t0) < h/.2.0 then dessine_univers();
   
   	 (* calcule les accélérations et avance d'un pas
   	  * arrête tout si on sort de l'univers, les cas de collision ne sont
   	  * pas testés, une collision parfaite étant improbable.
   	  *)
   	 ajuste_acc_approx();
   	 for i=0 to N-1 do let c = univers.(i) in
	    c.oldpos <- c.pos;
	    c.oldvel <- c.vel;
	    c.pos    <- {x = c.pos.x +. h*.c.vel.x +. h2*.c.acc.x;
                         y = c.pos.y +. h*.c.vel.y +. h2*.c.acc.y};
	    c.vel    <- {x = c.vel.x +. h*.c.acc.x;
	                 y = c.vel.y +. h*.c.acc.y};
	    if (abs_float(c.pos.x) > taille_univers/.2.0) or
	       (abs_float(c.pos.y) > taille_univers/.2.0) then raise fin
   	 done;
   	 construire_arbre();
   	 barycentres(!arbre_univers);
   
   	 (* avance d'un deuxième pas et fait la moyenne avec la position initiale *)
   	 ajuste_acc_approx();
   	 for i=0 to N-1 do let c = univers.(i) in
	    c.pos <- {x = (c.pos.x +. h*.c.vel.x +. h2*.c.acc.x +. c.oldpos.x)/.2.0;
	              y = (c.pos.y +. h*.c.vel.y +. h2*.c.acc.y +. c.oldpos.y)/.2.0};
	    c.vel <- {x = (c.vel.x +. h*.c.acc.x +. c.oldvel.x)/.2.0;
	              y = (c.vel.y +. h*.c.acc.y +. c.oldvel.y)/.2.0};
	    if (abs_float(c.pos.x) > taille_univers/.2.0) or
	       (abs_float(c.pos.y) > taille_univers/.2.0) then raise fin
   	 done;
   	 
   	 (* trace les segments *)
   	 for i=0 to N-1 do let c = univers.(i) in
	    printf "%s %f %f %f %f l\n" c.couleur c.oldpos.x c.oldpos.y c.pos.x c.pos.y
   	 done;
	    
   	 t := !t+.h
      done
      with fin -> ()
   end;

   printf "showpage\n";
   flush stdout
   
end;;