surface.scheme

; A surface is a function from 6 numbers to either #f or a list of 7 numbers.
; The 6 numbers are:
;   x,y,z:             the coordinates of a starting point;
;   dx,dy,dz:          the direction of a ray.
; The 7 numbers are:
;   root:              the factor to multiply dx,dy,dz by;
;   new-x,new-y,new-z: the intersection of the ray with the surface
;                      (i.e., x+dx*root, y+dy*root, z+dz*root, respectively);
;   out-x,out-y,out-z: the normal vector pointing away from the surface.
; If the ray does not intersect the surface, then the function returns #f.
; The dx,dy,dz input must not be all zero.  The out-x,out-y,out-z output must
; not all be zero.  The root output must be positive.

(define (xy z0)
  ; The x-y plane at z = z0
  (lambda (x y z dx dy dz)
    (if (zero? dz) #f
      (let ([root (/ (- z0 z) dz)])
        (if (<= root 0) #f
          (list root (+ x (* dx root)) (+ y (* dy root)) z0 0 0 (- dz)))))))

(define (sphere x0 y0 z0 r)
  ; The sphere centered at (x0,y0,z0) with radius r
  (let ([sqr-r (square r)])
    (lambda (x y z dx dy dz)
      (let ([x1 (- x x0)] [y1 (- y y0)] [z1 (- z z0)])
        (let ([a (square-norm dx dy dz)] ; positive
              [neg-half-b (- 0 (* dx x1) (* dy y1) (* dz z1))]
              [c (- (square-norm x1 y1 z1) sqr-r)])
          (let ([discriminant (- (square neg-half-b) (* a c))])
            (if (negative? discriminant) #f
              (let ([sqrt-discriminant (sqrt discriminant)])
                (let ([root-more (/ (+ neg-half-b sqrt-discriminant) a)])
                  (if (<= root-more 0) #f
                    (let ([root-less (/ (- neg-half-b sqrt-discriminant) a)])
                      (if (<= root-less 0)
                        (let ([new-x (+ x (* dx root-more))]
                              [new-y (+ y (* dy root-more))]
                              [new-z (+ z (* dz root-more))])
                          (list root-more new-x new-y new-z
                            (- x0 new-x) (- y0 new-y) (- z0 new-z)))
                        (let ([new-x (+ x (* dx root-less))]
                              [new-y (+ y (* dy root-less))]
                              [new-z (+ z (* dz root-less))])
                          (list root-less new-x new-y new-z
                            (- new-x x0) (- new-y y0) (- new-z z0)))))))))))))))