{-# LANGUAGE MultiParamTypeClasses #-} 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