-- |

-- Module:      Data.Geo.Jord.Triangle

-- Copyright:   (c) 2020 Cedric Liegeois

-- License:     BSD3

-- Maintainer:  Cedric Liegeois <ofmooseandmen@yahoo.fr>

-- Stability:   experimental

-- Portability: portable

--

-- Types and functions for working with triangles at the surface of a __spherical__ celestial body.

--

-- In order to use this module you should start with the following imports:

--

-- @

-- import qualified Data.Geo.Jord.Geodetic as Geodetic

-- import qualified Data.Geo.Jord.Triangle as Triangle

-- @

module Data.Geo.Jord.Triangle
    ( Triangle
    , vertex0
    , vertex1
    , vertex2
    , make
    , unsafeMake
    , centroid
    , circumcentre
    , contains
    ) where

import Control.Monad (join)
import Data.Maybe (fromJust)

import Data.Geo.Jord.Geodetic (HorizontalPosition)
import qualified Data.Geo.Jord.Geodetic as Geodetic
import qualified Data.Geo.Jord.GreatCircle as GreatCircle
import qualified Data.Geo.Jord.Math3d as Math3d
import Data.Geo.Jord.Model (Spherical)

-- | A triangle whose vertices are horizontal geodetic positions.

data Triangle a =
    Triangle (HorizontalPosition a) (HorizontalPosition a) (HorizontalPosition a)
    deriving (Triangle a -> Triangle a -> Bool
(Triangle a -> Triangle a -> Bool)
-> (Triangle a -> Triangle a -> Bool) -> Eq (Triangle a)
forall a. Model a => Triangle a -> Triangle a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Triangle a -> Triangle a -> Bool
$c/= :: forall a. Model a => Triangle a -> Triangle a -> Bool
== :: Triangle a -> Triangle a -> Bool
$c== :: forall a. Model a => Triangle a -> Triangle a -> Bool
Eq, Int -> Triangle a -> ShowS
[Triangle a] -> ShowS
Triangle a -> String
(Int -> Triangle a -> ShowS)
-> (Triangle a -> String)
-> ([Triangle a] -> ShowS)
-> Show (Triangle a)
forall a. Model a => Int -> Triangle a -> ShowS
forall a. Model a => [Triangle a] -> ShowS
forall a. Model a => Triangle a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Triangle a] -> ShowS
$cshowList :: forall a. Model a => [Triangle a] -> ShowS
show :: Triangle a -> String
$cshow :: forall a. Model a => Triangle a -> String
showsPrec :: Int -> Triangle a -> ShowS
$cshowsPrec :: forall a. Model a => Int -> Triangle a -> ShowS
Show)

-- | First vertex of given triangle.

vertex0 :: (Spherical a) => Triangle a -> HorizontalPosition a
vertex0 :: Triangle a -> HorizontalPosition a
vertex0 (Triangle HorizontalPosition a
v HorizontalPosition a
_ HorizontalPosition a
_) = HorizontalPosition a
v

-- | Second vertex of given triangle.

vertex1 :: (Spherical a) => Triangle a -> HorizontalPosition a
vertex1 :: Triangle a -> HorizontalPosition a
vertex1 (Triangle HorizontalPosition a
_ HorizontalPosition a
v HorizontalPosition a
_) = HorizontalPosition a
v

-- | Third vertex of given triangle.

vertex2 :: (Spherical a) => Triangle a -> HorizontalPosition a
vertex2 :: Triangle a -> HorizontalPosition a
vertex2 (Triangle HorizontalPosition a
_ HorizontalPosition a
_ HorizontalPosition a
v) = HorizontalPosition a
v

-- | Triangle from given vertices. Returns 'Nothing' if some vertices are equal or some are antipodes of others.

make ::
       (Spherical a)
    => HorizontalPosition a
    -> HorizontalPosition a
    -> HorizontalPosition a
    -> Maybe (Triangle a)
make :: HorizontalPosition a
-> HorizontalPosition a
-> HorizontalPosition a
-> Maybe (Triangle a)
make HorizontalPosition a
v0 HorizontalPosition a
v1 HorizontalPosition a
v2
    | HorizontalPosition a
v0 HorizontalPosition a -> HorizontalPosition a -> Bool
forall a. Eq a => a -> a -> Bool
== HorizontalPosition a
v1 Bool -> Bool -> Bool
|| HorizontalPosition a
v1 HorizontalPosition a -> HorizontalPosition a -> Bool
forall a. Eq a => a -> a -> Bool
== HorizontalPosition a
v2 Bool -> Bool -> Bool
|| HorizontalPosition a
v2 HorizontalPosition a -> HorizontalPosition a -> Bool
forall a. Eq a => a -> a -> Bool
== HorizontalPosition a
v0 = Maybe (Triangle a)
forall a. Maybe a
Nothing
    | HorizontalPosition a
a0 HorizontalPosition a -> HorizontalPosition a -> Bool
forall a. Eq a => a -> a -> Bool
== HorizontalPosition a
v1 Bool -> Bool -> Bool
|| HorizontalPosition a
a0 HorizontalPosition a -> HorizontalPosition a -> Bool
forall a. Eq a => a -> a -> Bool
== HorizontalPosition a
v2 Bool -> Bool -> Bool
|| HorizontalPosition a
a1 HorizontalPosition a -> HorizontalPosition a -> Bool
forall a. Eq a => a -> a -> Bool
== HorizontalPosition a
v2 = Maybe (Triangle a)
forall a. Maybe a
Nothing
    | Bool
otherwise = Triangle a -> Maybe (Triangle a)
forall a. a -> Maybe a
Just (HorizontalPosition a
-> HorizontalPosition a -> HorizontalPosition a -> Triangle a
forall a.
HorizontalPosition a
-> HorizontalPosition a -> HorizontalPosition a -> Triangle a
Triangle HorizontalPosition a
v0 HorizontalPosition a
v1 HorizontalPosition a
v2)
  where
    a0 :: HorizontalPosition a
a0 = HorizontalPosition a -> HorizontalPosition a
forall a. Model a => HorizontalPosition a -> HorizontalPosition a
Geodetic.antipode HorizontalPosition a
v0
    a1 :: HorizontalPosition a
a1 = HorizontalPosition a -> HorizontalPosition a
forall a. Model a => HorizontalPosition a -> HorizontalPosition a
Geodetic.antipode HorizontalPosition a
v1

-- | Triangle from given vertices. This is unsafe, if any vertices are equal or some are antipodes of others,

-- the resulting triangle is actually undefined.

unsafeMake ::
       (Spherical a)
    => HorizontalPosition a
    -> HorizontalPosition a
    -> HorizontalPosition a
    -> Triangle a
unsafeMake :: HorizontalPosition a
-> HorizontalPosition a -> HorizontalPosition a -> Triangle a
unsafeMake = HorizontalPosition a
-> HorizontalPosition a -> HorizontalPosition a -> Triangle a
forall a.
HorizontalPosition a
-> HorizontalPosition a -> HorizontalPosition a -> Triangle a
Triangle

-- | @contains t p@ returns 'True' if position @p@ is enclosed by the vertices of triangle

-- @t@ - see 'GreatCircle.enclosedBy'.

contains :: (Spherical a) => Triangle a -> HorizontalPosition a -> Bool
contains :: Triangle a -> HorizontalPosition a -> Bool
contains (Triangle HorizontalPosition a
v0 HorizontalPosition a
v1 HorizontalPosition a
v2) HorizontalPosition a
p = HorizontalPosition a -> [HorizontalPosition a] -> Bool
forall a.
Spherical a =>
HorizontalPosition a -> [HorizontalPosition a] -> Bool
GreatCircle.enclosedBy HorizontalPosition a
p [HorizontalPosition a
v0, HorizontalPosition a
v1, HorizontalPosition a
v2]

-- | Computes the centroid of the given triangle: the position which is the intersection of the three medians of

-- the triangle (each median connecting a vertex with the midpoint of the opposite side).

--

-- The centroid is always within the triangle.

centroid :: (Spherical a) => Triangle a -> HorizontalPosition a
centroid :: Triangle a -> HorizontalPosition a
centroid (Triangle HorizontalPosition a
v0 HorizontalPosition a
v1 HorizontalPosition a
v2) = Maybe (HorizontalPosition a) -> HorizontalPosition a
forall a. HasCallStack => Maybe a -> a
fromJust Maybe (HorizontalPosition a)
c -- this is safe unless triangle was created using unsafeMake.

  where
    m1 :: Maybe (MinorArc a)
m1 = [HorizontalPosition a] -> Maybe (HorizontalPosition a)
forall a.
Spherical a =>
[HorizontalPosition a] -> Maybe (HorizontalPosition a)
GreatCircle.mean [HorizontalPosition a
v1, HorizontalPosition a
v2] Maybe (HorizontalPosition a)
-> (HorizontalPosition a -> Maybe (MinorArc a))
-> Maybe (MinorArc a)
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= HorizontalPosition a -> HorizontalPosition a -> Maybe (MinorArc a)
forall a.
Spherical a =>
HorizontalPosition a -> HorizontalPosition a -> Maybe (MinorArc a)
GreatCircle.minorArc HorizontalPosition a
v0
    m2 :: Maybe (MinorArc a)
m2 = [HorizontalPosition a] -> Maybe (HorizontalPosition a)
forall a.
Spherical a =>
[HorizontalPosition a] -> Maybe (HorizontalPosition a)
GreatCircle.mean [HorizontalPosition a
v2, HorizontalPosition a
v0] Maybe (HorizontalPosition a)
-> (HorizontalPosition a -> Maybe (MinorArc a))
-> Maybe (MinorArc a)
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= HorizontalPosition a -> HorizontalPosition a -> Maybe (MinorArc a)
forall a.
Spherical a =>
HorizontalPosition a -> HorizontalPosition a -> Maybe (MinorArc a)
GreatCircle.minorArc HorizontalPosition a
v1
    c :: Maybe (HorizontalPosition a)
c = Maybe (Maybe (HorizontalPosition a))
-> Maybe (HorizontalPosition a)
forall (m :: * -> *) a. Monad m => m (m a) -> m a
join (MinorArc a -> MinorArc a -> Maybe (HorizontalPosition a)
forall a.
Spherical a =>
MinorArc a -> MinorArc a -> Maybe (HorizontalPosition a)
GreatCircle.intersection (MinorArc a -> MinorArc a -> Maybe (HorizontalPosition a))
-> Maybe (MinorArc a)
-> Maybe (MinorArc a -> Maybe (HorizontalPosition a))
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Maybe (MinorArc a)
m1 Maybe (MinorArc a -> Maybe (HorizontalPosition a))
-> Maybe (MinorArc a) -> Maybe (Maybe (HorizontalPosition a))
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> Maybe (MinorArc a)
m2)

-- | The circumcentre of the triangle: the position which is equidistant from all three vertices.

--

-- The circumscribed circle or circumcircle of a triangle is a circle which passes through all

-- the vertices of the triangle; The circumcentre is not necessarily inside the triangle.

--

-- Thanks to STRIPACK: http://orion.math.iastate.edu/burkardt/f_src/stripack/stripack.f90

circumcentre :: (Spherical a) => Triangle a -> HorizontalPosition a
circumcentre :: Triangle a -> HorizontalPosition a
circumcentre (Triangle HorizontalPosition a
v0 HorizontalPosition a
v1 HorizontalPosition a
v2) = V3 -> a -> HorizontalPosition a
forall a. Model a => V3 -> a -> HorizontalPosition a
Geodetic.nvectorPos' V3
cu (HorizontalPosition a -> a
forall a. Model a => HorizontalPosition a -> a
Geodetic.model HorizontalPosition a
v0)
  where
    e0 :: V3
e0 = V3 -> V3 -> V3
Math3d.subtract (HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
v1) (HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
v0)
    e1 :: V3
e1 = V3 -> V3 -> V3
Math3d.subtract (HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
v2) (HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
v0)
    cu :: V3
cu = V3 -> V3 -> V3
Math3d.cross V3
e0 V3
e1