1 -- 
    2 --      Patricia Fasel
    3 --      Los Alamos National Laboratory
    4 --      1990 August
    5 --
    6 module Distance (distSurf) where
    7 
    8 import GamtebType
    9 import Consts
   10 
   11 -- find distance of a particle to a surface
   12 
   13 distSurf :: Point -> Point -> (Value, Indx)
   14 distSurf (p1@(x,y,z)) (p2@(u,v,w)) =
   15         (dSurf+small, surf)
   16         where
   17             (d1, s1) = ((distPlane y v 0), 1)
   18             (d2, s2) = ((distCyl p1 p2), 2)
   19             (d3, s3) = ((distPlane y v cylLen), 3)
   20             (d4, s4) = ((distPlane y v cylLen2), 4)
   21             (dSurf, surf) = minP (minP (d1,s1) (d2,s2)) (minP (d3,s3) (d4,s4))
   22             minP (d, j) (d', j') =
   23                 if (d < d')
   24                     then (d, j)
   25                     else (d', j')
   26 
   27 
   28 -- find distance to a cylinder
   29 
   30 distCyl :: Point -> Point -> Value
   31 distCyl (x,y,z) (u,v,w)
   32         | (u*u + w*w) == 0     = big             -- w*w used to be v*v  LA
   33         | (u /= 0)         = 
   34                 let
   35                     m = w/u
   36                     b = z - m*x
   37                     s = m*m + 1
   38                     r = sqrt (s - b*b)
   39                     x' = if (u > 0)
   40                             then (-m*b + r) / s
   41                             else (-m*b - r) / s
   42                 in
   43                 (x'-x) / u
   44         | (u == 0 && v /= 0)   = 
   45                 let
   46                     m = w/v
   47                     b = z - m*y
   48                     r = sqrt (1 - x*x)
   49                     y' = if (v > 0)
   50                             then (r-b) / m
   51                             else (-r-b) / m
   52                 in
   53                 (y'-y) / v
   54         | w > 0                     = (sqrt (1 - x*x)) - z
   55         | otherwise       = -(sqrt (1 - x*x)) - z
   56 
   57 
   58 -- find distance of a particle to a plane
   59 
   60 distPlane :: Coord -> Coord -> Value -> Value
   61 distPlane y v yPlane 
   62         | v == 0       = big
   63         | y >= yPlane  = big
   64         | otherwise    = (yPlane-y) / v