home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Source Code / C / Applications / SML⁄NJ 93+ / Documentation / examples / ray.sml < prev    next >
Encoding:
Text File  |  1995-12-30  |  2.2 KB  |  61 lines  |  [TEXT/R*ch]
  1. (* ray.sml *)
  2. (* Find the intersection of a ray with a cylinder in three-space.  If
  3.    there is more than one intersection, return the parameter of
  4.    intersection closest to the start of the ray.  If there is no
  5.    intersection, return NONE.
  6.  
  7.    The ray is specified by its starting point (rayStart) and direction
  8.    vector (rayVec).  The cylinder is specified by a starting point (cStart)
  9.    a length vector (cVec), and a radius (cRad).  The cylinder is finite;
  10.    it does not extend beyond the direction vector at either end.
  11.  
  12. *)
  13.  
  14. local open Real
  15.       val EP = 0.0000001
  16.  
  17.   fun newton f x =
  18.     let val y = f x 
  19.      in if abs(y)>EP then newton f (x+EP*y/(y-f(x+EP))) else x 
  20.     end
  21.  
  22.   infix 7 @   fun (a,b,c) @ (d,e,f) = a*d+b*e+c*f     (* dot product *)
  23.   infix 7 *!  fun (a,b,c) *! t = (a*t,b*t,c*t)        (* vector * scalar *)
  24.   infix 6 +!  fun (a,b,c) +! (d,e,f) = (a+d,b+e,c+f)  (* vector addition *)
  25.   infix 6 -!  fun (a,b,c) -! (d,e,f) = (a-d,b-e,c-f)  (* vector subtraction *)
  26.  
  27.   fun along (C,D) t = C +! (D *! t)
  28.  
  29. in 
  30.  
  31. fun cylhit (ray as (rayStart, rayVec), cyl as (cStart, cVec), cRad) =
  32. let val rayStart = rayStart -! cStart
  33.     fun cylpos p = (cVec@p) / (cVec@cVec)
  34.         (* parameter of closest point on cVec to the point p *)
  35.  
  36.     fun cyldist p = let val v = p -! (along cyl (cylpos p)) in v@v end
  37.         (* distance of p from closest point on cVec *)
  38.  
  39.     fun cyldist' p = p*!(cVec@cVec) -! cVec*!(p@cVec)
  40.     (* gradient of the cyldist function *)
  41.  
  42.     val tf = newton(fn t => cyldist' (along ray t) @ rayVec) 0.0
  43.     (* parameter along the ray of the closest point to cVec *)
  44.  
  45.  in if cyldist (along ray tf) > cRad*cRad then NONE else
  46.    let val t0 = newton(fn t => cyldist(along ray t)-(cRad*cRad)) 0.0
  47.     (* one of the two intersection points *)
  48.        val t1 = newton(fn t => cyldist(along ray t)-(cRad*cRad)) (tf+tf-t0)
  49.     (* the other intersection point *)
  50.        fun inRange v = v >= 0.0 andalso v <= 1.0
  51.     in case (t0 > 0.0 andalso inRange (cylpos (along ray t0)),
  52.          t1 > 0.0 andalso inRange (cylpos (along ray t1)))
  53.        of (false,false) => NONE
  54.         | (true,false) => SOME t0
  55.         | (false,true) => SOME t1
  56.         | (true,true) => SOME (if (t0<t1) then t0 else t1)
  57.    end
  58. end
  59.  
  60. end
  61.