Vers Détermination de la trajectoire de comètes (Présentation)

include "A:/Tracé de trajectoires.ml";;

let rec divise =function
    |[]->([],[])
    |[e]->([e],[])
    |a::b::r->let (m1,m2)=divise r in (a::m1,b::m2);;

let rec fusion_pts lista listb=match lista,listb with
    |l,[]->lista
    |[],l->listb
    |(a::r),(b::s)->if (a.x)<>b.x then if (a.x)<(b.x) then a::(fusion_pts r listb)
                                                                                  else b::(fusion_pts lista s)
         else if (a.y)<(b.y) then a::(fusion_pts r listb)
                                                                           else b::(fusion_pts lista s);;
let rec tri_pts=fun
    |[]->[]
    |[e]->[e]
    |l->let (m1,m2)=divise l in
                     fusion_pts (tri_pts m1) (tri_pts m2);;

let distance pa pb=sqrt((((pa.x)-.(pb.x))*.((pa.x)-.(pb.x)))+.(((pa.y)-.(pb.y))*.((pa.y)-.(pb.y))));;

let intersection lista listb=
    if lista=[] or listb=[] then
      failwith "liste vide"
    else
      let ha=ref(hd lista) and hb=ref(hd listb)
      and resulta=ref(hd lista) and resultb=ref(hd listb)
      and qa=ref(tl lista) and qb=ref(tl listb) in
      let epsilon=ref(distance (!ha) (!hb)) in
      while ((!qa)<>[] && (!qb)<>[]) do
        if ((!ha).x)<>((!hb).x) then
        
          if ((!ha).x)<((!hb).x) then
            
            if (distance (!ha) (!hb))<(!epsilon) then
              begin
                epsilon:=(distance (!ha) (!hb));
                resulta:=(!ha);
                resultb:=(!hb);
                ha:=(hd(!qa));
                qa:=(tl(!qa));
              end
            else
              begin
                ha:=(hd(!qa));
                qa:=(tl(!qa));
              end
            
          else
            
            if (distance (!ha) (!hb))<(!epsilon) then
              begin
                epsilon:=(distance (!ha) (!hb));
                resulta:=(!ha);
                resultb:=(!hb);
                hb:=(hd(!qb));
                qb:=(tl(!qb));
              end
            else
              begin
                hb:=(hd(!qb));
                qb:=(tl(!qb));
              end
        
        else
        
          if ((!ha).y)<=((!hb).y) then
            
            if (distance (!ha) (!hb))<(!epsilon) then
              begin
                epsilon:=(distance (!ha) (!hb));
                resulta:=(!ha);
                resultb:=(!hb);
                ha:=(hd(!qa));
                qa:=(tl(!qa));
              end
            else
              begin
                ha:=(hd(!qa));
                qa:=(tl(!qa));
              end
              
          else
            
            if (distance (!ha) (!hb))<(!epsilon) then
              begin
                epsilon:=(distance (!ha) (!hb));
                resulta:=(!ha);
                resultb:=(!hb);
                hb:=(hd(!qb));
                qb:=(tl(!qb));
              end
            else
              begin
                hb:=(hd(!qb));
                qb:=(tl(!qb));
              end;
      done;
    (!epsilon,(!resulta),(!resultb));;

let fusion_intersection liste1 liste2=
    let (epsilon,p1,p2)=intersection (tri_pts liste1) (tri_pts liste2) in
                   (epsilon,{x=(((p1.x)+.(p2.x))/.2.);y=(((p1.y)+.(p2.y))/.2.)});;

let equ_hyperbole f1 f2 m=
      let c=((distance f1 f2)/.2.)
      and a=(abs_float((distance f1 m)-.(distance f2 m))/.2. )
      and aap=(atan( ((f2.y)-.(f1.y))/.((f2.x)-.(f1.x)) )) in
{e=(c/.a);p=(abs_float( ((c*.c)-.(a*.a))/.a ));aap=aap};;

let equ_ellipse f1 f2 m =
      let c = ((distance f1 f2)/.2.)
      and a = (((distance f1 m)+.(distance f2 m))/.2.)
      and aap = ( atan( (f2.y)/.(f2.x) )) in
{e=(c/.a);p=(abs_float( ((c*.c)-.(a*.a))/.a ));aap=aap};;

let branche_hyperbole hyperbole foyer n =
     let p=hyperbole.p
     and e=hyperbole.e
     and aap=hyperbole.aap in
     let liste=ref[]
     and t=ref( (acos(-.1./.(e)))+.aap )
     and dt=( (2.*.(acos(-.1./.(e))))/.(float_of_int n) ) in
       for i=n downto 1 do
          let r = p/.(1. +. e*.(cos((!t)-.aap))) in
          let x = (r)*.(cos(!t))
          and y = (r)*.(sin(!t)) in
          liste := {x=(x+.(foyer.x));y=(y+.(foyer.y))}::(!liste);
          t := (!t)-.dt
       done;
(!liste);;