let equations move =
    move.(3* !root+0) <- Some 0. ;
    move.(3* !root+1) <- Some 0. ;
    move.(3* !root+2) <- Some 0. ;
    let m =
      let m = ref 0 in
        Array.iter (fun a -> if a <> None then incr m) move ;
        !m
    in
    let p = Array.length pos in
    let b = nb_bones in
    let l = 2*b+3*p in (* vector's length *)
    let eq = Array.make (m+3*b) [||] in
    let cst = Array.make (m+3*b) 0. in
    let eq_i = ref 0 in

      (* Moving constraints *)
      Array.iteri
        (fun i -> function
           | None -> ()
           | Some d ->
               let e = Array.make l 0. in
                 e.(2*b+i) <- 1. ;
                 eq.(!eq_i) <- e ; cst.(!eq_i) <- d ; incr eq_i )
        move ;

      (* Skeleton constraints *)
      for k = 0 to b-1 do
        let (i,j,d) = bones.(k) in
        let {theta=theta;phi=phi} = angles.(k) in
        let ex = Array.make l 0. in
        let ey = Array.make l 0. in
        let ez = Array.make l 0. in
        let o = 2*k in

          (* Here we encode
           * dx_j = dx_i - dtheta sin phi sin theta + dphi cos phi cos theta *)
          ex.(o+0) <- -. d *. (sin phi) *. (sin theta) ;
          ex.(o+1) <-    d *. (cos phi) *. (cos theta) ;
          ex.(2*b+3*j+0) <- -1. ;
          ex.(2*b+3*i+0) <- 1. ;
          eq.(!eq_i) <- ex ; incr eq_i ;

          (* Similar stuff for dy ... *)
          ey.(o+0) <- d *. (sin phi) *. (cos theta) ;
          ey.(o+1) <- d *. (cos phi) *. (sin theta) ;
          ey.(2*b+3*j+1) <- -1. ;
          ey.(2*b+3*i+1) <- 1. ;
          eq.(!eq_i) <- ey ; incr eq_i ;

          (* ... and dz. *)
          ez.(o+1) <- -. d *. (sin phi) ;
          ez.(2*b+3*j+2) <- -1. ;
          ez.(2*b+3*i+2) <- 1. ;
          eq.(!eq_i) <- ez ; incr eq_i ;

      done ;

      assert (!eq_i = m+3*b) ;
      (eq, cst)