module Geometry.Hyperbolic.TwoD.HalfSpace
( Point(..)
, Geodesic(..)
) where
import Data.Fixed (mod')
import Geometry.Class
import qualified Geometry.Flat.TwoD.Space as E
data Point = Point{ pX, pY :: !Double } deriving Show
data Geodesic
= Line{ gDir :: !Bool, gCX :: !Double }
| Arc { gDir :: !Bool, gCX :: !Double, gR :: !Double }
deriving Show
instance Geometry Point Geodesic where
distance p q = case geodesic p q of
Line{} -> abs (log (pY p / pY q))
g@Arc{} ->
let a = gCX g gR g
b = gCX g + gR g
pa = sqrt ((pX p a) * (pX p a) + pY p * pY p)
pb = sqrt ((pX p b) * (pX p b) + pY p * pY p)
qa = sqrt ((pX q a) * (pX q a) + pY q * pY q)
qb = sqrt ((pX q b) * (pX q b) + pY q * pY q)
in abs (log ((pa / pb) / (qa / qb)))
angle _ _ _ = error "Geometry.Hyperbolic.TwoD.HalfPlane.angle"
geodesic p q
| abs (pX p pX q) < eps =
Line{ gCX = 0.5 * (pX p + pX q), gDir = pY p < pY q }
| otherwise =
let cx = 0.5 * (pX q * pX q + pY q * pY q pX p * pX p pY p * pY p) / (pX q pX p)
r = sqrt ((pX p cx) * (pX p cx) + pY p * pY p)
in Arc{ gCX = cx, gDir = pX p < pX q, gR = r }
rotate a p q =
let g = geodesic p q
b = angleAt g p
h = fromPointAngle p (b + a)
d = distance p q
in translate d h p
translate d g p = atDist g p (abs d) (d >= 0)
instance Embedding Point Point where
embed = id
model = Just
instance Embedding E.Point Point where
embed (Point x y) = E.Point x y
model (E.Point x y)
| y > 0 = Just (Point x y)
| otherwise = Nothing
eps :: Double
eps = 1e-12
atParam :: Geodesic -> Double -> Point
atParam g@Line{} t = Point (gCX g) t
atParam g@Arc {} t = Point (gR g * cos (pi t) + gCX g) (gR g * sin t)
paramAt :: Geodesic -> Point -> Double
paramAt Line{} p = pY p
paramAt g@Arc {} p = pi atan2 (pY p) (pX p gCX g)
angleAt :: Geodesic -> Point -> Double
angleAt g@Line{} _ = (if gDir g then (subtract pi) else id) (pi/2)
angleAt g@Arc {} p = (if gDir g then (pi +) else id) (atan2 (pY p) (pX p gCX g) + pi/2)
fromPointAngle :: Point -> Double -> Geodesic
fromPointAngle p a
| abs ((a `mod'` pi) pi/2) < eps =
Line{ gCX = pX p, gDir = cos a >= 0 }
| otherwise =
let b = a + pi/2
co = cos b
si = sin b
r = abs (pY p / si)
cx = (pY p * co pX p * si) / si
in Arc{ gCX = cx, gDir = cos a >= 0, gR = r }
atDist :: Geodesic -> Point -> Double -> Bool -> Point
atDist g@Line{} p d dir
| dir /= gDir g = binarySearch g p d (0 + eps) pp (<)
| otherwise = binarySearch g p d pp (pp * 100) (>)
where pp = paramAt g p
atDist g@Arc {} p d dir
| dir /= gDir g = binarySearch g p d (0 + eps) pp (<)
| otherwise = binarySearch g p d pp (pi eps) (>)
where pp = paramAt g p
binarySearch :: Geodesic -> Point -> Double -> Double -> Double -> (Double -> Double -> Bool) -> Point
binarySearch g p d lo0 hi0 cmp = go lo0 hi0
where
go lo hi =
let mid = 0.5 * (lo + hi)
q = atParam g mid
m = distance p q
in if hi lo < eps
then q
else if m `cmp` d
then go lo mid
else go mid hi