MODULE R3;

(* Functions on 3D real vectors. *)

CONST Origin = [0, 0, 0];
(* "Origin" is the origin in 3-space. *)

(* \section{Arithmetic Functions} *)

FUNC r = Plus(p, q) IS 
  (E px, py, pz, qx, qy, qz :: 
    p = [px, py, pz] AND 
    q = [qx, qy, qz] AND 
    r = [px + qx, py + qy, pz + qz]) 
END;

(* "r" is the vector sum "p" + "q". *)

FUNC r = Minus(p, q) IS 
  (E px, py, pz, qx, qy, qz :: 
    p = [px, py, pz] AND 
    q = [qx, qy, qz] AND 
    r = [px - qx, py - qy, pz - qz]) 
END;

(* "r" is the vector difference "p" - 
   "q". *)

FUNC p = Times(t, q) IS 
  (E qx, qy, qz :: 
    q = [qx, qy, qz] AND 
    p = [t * qx, t * qy, t * qz]) 
END;

(* "p" is the product of the scalar "t" 
   and the vector "q". *)

FUNC len2 = Length2(v) IS 
  (E vx, vy, vz :: 
    v = [vx, vy, vz] AND 
    len2 = 
      (vx * vx) + (vy * vy) + (vz * vz)) 
END;

(* "len2" is the square of the L2 norm 
   of the vector "v". *)

FUNC len = Length(v) IS 
  len = Math.Sqrt(Length2(v)) 
END;

(* "len" is the L2 norm of the vector 
   "v", namely, its Euclidean length. *)

FUNC u = Normalize(v) IS
  Times(Length(v), u) = v
END;

(* "u" is the vector "v" normalized
   to have unit length. *)

FUNC t = Dot(p, q) IS 
  (E px, py, pz, qx, qy, qz :: 
    p = [px, py, pz] AND 
    q = [qx, qy, qz] AND 
    t = 
      (px * qx) + (py * qy) + (pz * qz)) 
END;

(* The scalar "t" is the dot-product of 
   "p" and "q". *)

FUNC r = Cross(p, q) IS 
  (E px, py, pz, qx, qy, qz :: 
    p = [px, py, pz] AND 
    q = [qx, qy, qz] AND 
    r = 
      [(py * qz) - (qy * pz), 
       (pz * qx) - (qz * px), 
       (px * qy) - (qx * py)]) 
END;

(* "r" is the cross-product of "p" and 
   "q". *)

(* \section{Projection Functions} *)

FUNC x = X(p) IS 
  (E y, z :: p = [x, y, z]) 
END;

(* "x" is the x-coordinate of the 3D 
   point "p". *)

FUNC y = Y(p) IS 
  (E x, z :: p = [x, y, z]) 
END;

(* "y" is the y-coordinate of the 3D 
   point "p". *)

FUNC z = Z(p) IS 
  (E x, y :: p = [x, y, z]) 
END;

(* "z" is the y-coordinate of the 3D 
   point "p". *)

FUNC r = FromR2(p) IS 
  (E px, py :: 
    p = (px, py) AND r = [px, py, 0]) 
END;

(* Convert the R2 point "p" to an R3 
   point whose "z" coordinate is 0. *)

FUNC r = ToR2(p) IS 
  (E px, py, pz :: 
    p = [px, py, pz] AND r = (px, py)) 
END;

(* "r" is the R2 projection of the R3 
   point "p" onto the XY-plane. *)