(******************************************************************************
 ** RealGeometryPrims3d.sml
 ** sml
 **
 ** Franklin Chen and Guy Blelloch 
 ** An implemetation of GEOMETRY_PRIMS in 3d with Real coordinates
 **   This is specialized to Reals for efficiency reasons
 ******************************************************************************)

structure RealGeometryPrims3d
    :> GEOMETRY_PRIMS
    where Number = RealNumber
    =
    struct
        structure Number = RealNumber

	exception BadDimension
	exception NYI

        structure Vec =
            struct
                structure Number = Number
                structure Field = Number
                type vec = Number.t * Number.t * Number.t
		type t = vec

                val zero = (Number.zero, Number.zero, Number.zero)
                val infinity = (Number.posInf, Number.posInf, Number.posInf)

		fun fromSeq #[x, y, z] : t = (x, y, z)
		  | fromSeq _  = raise BadDimension

                fun toSeq (x,y,z) = #[x, y, z]

                fun (x, y, z) + (x', y', z') : t =
                    let
                        open Number
                    in
                        (x + x', y + y', z + z')
                    end

                fun ~ (x, y, z) : t =
                    let
                        open Number
                    in
                        (~x, ~y, ~z)
                    end

                fun (x, y, z) - (x', y', z')  =
		    let
			open Real
		    in
                        (x - x', y - y', z - z')
		    end

                fun scale ((x, y, z) : t, a : Number.t) : t =
                    let
                        open Number
                    in
                        (a * x, a * y, a * z)
                    end

                fun dot ((x, y, z), (x', y', z')) : Number.t =
                    let
                        open Number
                    in
                        x * x' + y * y' + z * z'
                    end

		fun cross #[ (x, y, z), (x', y', z') ] : t =
		    let
			open Number
		    in
			(y * z' - z * y', z * x' - x * z', x * y' - y * x')
		    end
		  | cross _ = raise BadDimension


                (* Derived form. *)
                fun length (p : t) : Number.t =
                    Number.sqrt (dot (p, p))

		fun norm (p : t) : Number.t =
                    Number.sqrt (dot (p, p))

                fun == ((x, y, z), (x', y', z')) =
                    let
                        open Number
                    in
                        == (x, x') andalso == (y, y') andalso == (z, z')
                    end

                fun != ((x, y, z), (x', y', z')) =
                    let
                        open Number
                    in
                        != (x, x') orelse != (y, y') orelse != (z, z')
                    end

		fun toString(x,y,z) =
		    ("(" ^ Number.toString(x) ^ 
		     ", " ^ Number.toString(y) ^ 
		     ", " ^ Number.toString(z) ^ ")")
            end

	structure Point = PointFromVec(structure Vec = Vec)

	fun minCoords((x1,y1,z1),(x2,y2,z2)) =
	    (Number.min(x1,x2),Number.min(y1,y2),Number.min(z1,z2))

	fun maxCoords((x1,y1,z1),(x2,y2,z2)) =
	    (Number.max(x1,x2),Number.max(y1,y2),Number.max(z1,z2))

	type plane = (Vec.t*Number.t)

	fun planeNormal (normal, offset) = normal

	fun planeFromPoints points = 
	    let
		open Number
		val #[(x0,y0,z0),(x1,y1,z1),(x2,y2,z2)] = points
		val (a0,a1,a2) = (x1-x0,y1-y0,z1-z0)
		val (b0,b1,b2) = (x2-x0,y2-y0,z2-z0)
		val (xv,yv,zv) = ((a1*b2 - a2*b1), 
				  (a2*b0 - a0*b2), 
				  (a0*b1 - a1*b0))
		val offset = (xv*x0+yv*y0+zv*z0)
	    in
		((xv,yv,zv),offset)
	    end

	fun planeDistance ((x,y,z), offset) point = 
	    let 
		open Number
		val (x2,y2,z2) = point
		val det = x*x2 + y*y2 + z*z2
	    in
		det-offset
	    end

	fun planeSide (normal, offset) point = 
	    case Number.compare((planeDistance (normal, offset) point),Number.zero) of
		LESS => Side.In
	      | GREATER => Side.Out
	      | EQUAL => Side.On

	fun planeEqual((n1,o1),(n2,o2)) = 
	    Vec.==(n1,n2) andalso Number.==(o1,o2)

	fun planeInvert(normal,offset) = 
	    (Vec.~ normal, offset)

	fun planeToString(a) = raise NYI
	fun planeFromNormal(a) = raise NYI

	type sphere = unit
	fun sphereFromPoints _ = raise NYI
	fun sphereSide _ _ = raise NYI

    end