{-# 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