module DSMC.Traceables
(
Body
, plane
, sphere
, cylinder
, cylinderFrustum
, cone
, coneFrustum
, intersect
, unite
, complement
, HitPoint(..)
, hitPoint
, HitSegment
, Trace
, trace
, inside
)
where
import Prelude hiding (Just, Nothing, Maybe, fst, snd)
import Data.Functor
import Data.Strict.Maybe
import Data.Strict.Tuple
import DSMC.Particles
import DSMC.Util
import DSMC.Util.Vector
data HitPoint = HitPoint !Double (Maybe Vec3)
deriving (Eq, Ord, Show)
type HitSegment = (Pair HitPoint HitPoint)
type Trace = [HitSegment]
infinityP :: Double
infinityP = (/) 1 0
infinityN :: Double
infinityN = infinityP
hitN :: HitPoint
hitN = (HitPoint infinityN Nothing)
hitP :: HitPoint
hitP = (HitPoint infinityP Nothing)
data Body = Plane !Vec3 !Double
| Sphere !Vec3 !Double
| Cylinder !Vec3 !Point !Double
| Cone !Vec3 !Point !Double !Matrix !Double !Double
| Union !Body !Body
| Intersection !Body !Body
| Complement !Body
deriving Show
plane :: Point -> Vec3 -> Body
plane p n = Plane nn (p .* nn)
where
nn = normalize n
sphere :: Vec3 -> Double -> Body
sphere o r = Sphere o r
cylinder :: Point -> Point -> Double -> Body
cylinder p1 p2 r = Cylinder (normalize $ p2 <-> p1) p1 r
cylinderFrustum :: Point -> Point -> Double -> Body
cylinderFrustum pb pt r =
intersect (plane pt axis)
(intersect (plane pb $ invert axis)
(cylinder pb pt r))
where
axis = pt <-> pb
cone :: Vec3 -> Point -> Double -> Body
cone a o h =
let
h' = cos $ (h * pi / 180)
n = normalize $ invert a
gamma = diag (h' * h')
m = addM (n `vxv` n) gamma
ta = tan $ h
odelta = n .* o
in
Cone n o h' m ta odelta
coneFrustum :: (Point, Double) -> (Point, Double) -> Body
coneFrustum (p1, r1) (p2, r2) =
let
(pb, rb, pt, rt) = case (r1 < r2) of
True -> (p2, r2, p1, r1)
False -> (p1, r1, p2, r2)
gap = pt <-> pb
height = norm gap
axis = normalize gap
dist = if rt == 0
then 0
else height / (rb / rt 1)
apex = pt <+> (axis .^ dist)
degs = atan (rb / (dist + norm (pt <-> pb))) * (180 / pi)
in
intersect (plane pt axis)
(intersect (plane pb $ invert axis)
(cone axis apex degs))
intersect :: Body -> Body -> Body
intersect !b1 !b2 = Intersection b1 b2
unite :: Body -> Body -> Body
unite !b1 !b2 = Union b1 b2
complement :: Body -> Body
complement !b = Complement b
trace :: Body -> Particle -> Trace
trace !b@(Plane n d) !p@(pos, v) =
let
!f = (n .* v)
in
if f == 0
then
if inside b p
then [(HitPoint infinityN Nothing) :!: (HitPoint infinityP Nothing)]
else []
else
let
!t = (pos .* n d) / f
in
if f > 0
then [(HitPoint t (Just n)) :!: (HitPoint infinityP Nothing)]
else [(HitPoint infinityN Nothing) :!: (HitPoint t (Just n))]
trace !(Sphere c r) !(pos, v) =
let
!d = pos <-> c
!roots = solveq (v .* v) (v .* d * 2) (d .* d r * r)
normal !u = normalize (u <-> c)
in
case roots of
Nothing -> []
Just (t1 :!: t2) ->
[HitPoint t1 (Just $ normal $ moveBy pos v t1) :!:
HitPoint t2 (Just $ normal $ moveBy pos v t2)]
trace !(Cylinder n c r) !(pos, v) =
let
d = (pos <-> c) >< n
e = v >< n
roots = solveq (e .* e) (d .* e * 2) (d .* d r * r)
normal u = normalize $ h <-> (n .^ (h .* n))
where h = u <-> c
in
case roots of
Nothing -> []
Just (t1 :!: t2) ->
[HitPoint t1 (Just $ normal $ moveBy pos v t1) :!:
HitPoint t2 (Just $ normal $ moveBy pos v t2)]
trace !(Cone n c _ m ta odelta) !(pos, v) =
let
delta = pos <-> c
c2 = dotM v v m
c1 = dotM v delta m
c0 = dotM delta delta m
roots = solveq c2 (2 * c1) c0
normal !u = normalize $ nx .^ (1 / ta) <-> ny .^ ta
where h = u <-> c
ny' = n .^ (n .* h)
ny = normalize ny'
nx = normalize $ h <-> ny'
in
case roots of
Nothing -> []
Just (t1 :!: t2) ->
let
pos1 = moveBy pos v t1
pos2 = moveBy pos v t2
in
case ((pos1 .* n odelta) > 0, (pos2 .* n odelta) > 0) of
(True, True) -> [HitPoint t1 (Just $ normal pos1) :!:
HitPoint t2 (Just $ normal pos2)]
(True, False) -> [HitPoint infinityN Nothing :!:
HitPoint t1 (Just $ normal pos1)]
(False, True) -> [HitPoint t2 (Just $ normal pos2) :!:
HitPoint infinityP Nothing]
(False, False) -> []
trace !(Intersection b1 b2) !p =
intersectTraces tr1 tr2
where
tr1 = trace b1 p
tr2 = trace b2 p
trace !(Union b1 b2) !p =
uniteTraces tr1 tr2
where
tr1 = trace b1 p
tr2 = trace b2 p
trace !(Complement b) !p =
complementTrace $ trace b p
uniteTraces :: Trace -> Trace -> Trace
uniteTraces u [] = u
uniteTraces u (v:t2) =
uniteTraces (unite1 u v) t2
where
merge :: HitSegment -> HitSegment -> HitSegment
merge !(a1 :!: b1) !(a2 :!: b2) = (min a1 a2) :!: (max b1 b2)
unite1 :: Trace -> HitSegment -> Trace
unite1 [] hs = [hs]
unite1 t@(hs1@(a1 :!: b1):tr') hs2@(a2 :!: b2)
| b1 < a2 = hs1:(unite1 tr' hs2)
| a1 > b2 = hs2:t
| otherwise = unite1 tr' (merge hs1 hs2)
intersectTraces :: Trace -> Trace -> Trace
intersectTraces tr1 tr2 =
let
overlap :: HitSegment -> HitSegment -> HitSegment
overlap !(a1 :!: b1) !(a2 :!: b2) = (max a1 a2) :!: (min b1 b2)
in
case tr2 of
[] -> []
(hs2@(a2 :!: b2):tr2') ->
case tr1 of
[] -> []
(hs1@(a1 :!: b1):tr1') ->
case (b1 < a2) of
True -> (intersectTraces tr1' tr2)
False ->
case (b2 < a1) of
True -> intersectTraces tr1 tr2'
False -> (overlap hs1 hs2):(intersectTraces tr1' tr2)
complementTrace :: Trace -> Trace
complementTrace ((sp@(HitPoint ts _) :!: ep):xs) =
start ++ (complementTrace' ep xs)
where
flipNormals :: HitSegment -> HitSegment
flipNormals !((HitPoint t1 n1) :!: (HitPoint t2 n2)) =
(HitPoint t1 (invert <$> n1)) :!: (HitPoint t2 (invert <$> n2))
start = if (isInfinite ts)
then []
else [flipNormals $ hitN :!: sp]
complementTrace' :: HitPoint -> Trace -> Trace
complementTrace' c ((a :!: b):tr) =
(flipNormals (c :!: a)):(complementTrace' b tr)
complementTrace' a@(HitPoint t _) [] =
if (isInfinite t)
then []
else [flipNormals (a :!: hitP)]
complementTrace [] = [hitN :!: hitP]
hitPoint :: Time -> Body -> Particle -> Maybe HitPoint
hitPoint !dt !b !p =
let
lastHit = [(HitPoint (dt) Nothing) :!: (HitPoint 0 Nothing)]
in
case (intersectTraces lastHit $ trace b p) of
[] -> Nothing
(hs:_) -> Just $ fst hs
inside :: Body -> Particle -> Bool
inside !(Plane n d) !(pos, _) = (pos .* n d) < 0
inside !(Sphere c r) !(pos, _) = (norm $ pos <-> c) < r
inside !(Cylinder n c r) !(pos, _) =
(norm $ h <-> (n .^ (h .* n))) < r
where
h = pos <-> c
inside !(Cone n c a _ _ _) !(pos, _) =
(n .* (normalize $ pos <-> c)) > a
inside !(Intersection b1 b2) !p = inside b1 p && inside b2 p
inside !(Union b1 b2) !p = inside b1 p || inside b2 p
inside !(Complement b) !p = not $ inside b p