-- |
-- Description: Manipulate coordinates in flat spacetime.
-- Copyright: Copyright (C) 2023 Yoo Chung
-- License: All rights reserved
-- Maintainer: web@chungyc.org
--
-- Defines coordinates for flat spacetime and
-- provides functions to manipulate and transform them.
-- It is intended to be a simple module sufficient enough
-- to assist in drawing spacetime diagrams.
module Physics.Spacetime.Flat
  ( -- * Coordinates
    Coordinate (..),

    -- * Simple operations
    plus,
    minus,
    scale,
    combine,
    approximates,

    -- * Spacetime interval
    interval,

    -- * Lorentz transformation
    transform,
    transformation,

    -- * Working with Diagrams
    vectorize2D,
    axes,
    axesWith,
    axesOptions,
    axesLength,
    axesLightcone,
    AxesOptions,
  )
where

import Diagrams.Backend.SVG
import Diagrams.Prelude hiding (interval, scale, transform)
import Diagrams.Transform.Matrix (fromMat22)
import GHC.Generics (Generic)

-- | The coordinates \((t, x, y, z)\) in flat spacetime.
--
-- The coordinates use Planck units.  In other words, \(c = 1\).
newtype Coordinate = Coordinate (Double, Double, Double, Double)
  deriving (Coordinate -> Coordinate -> Bool
(Coordinate -> Coordinate -> Bool)
-> (Coordinate -> Coordinate -> Bool) -> Eq Coordinate
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Coordinate -> Coordinate -> Bool
== :: Coordinate -> Coordinate -> Bool
$c/= :: Coordinate -> Coordinate -> Bool
/= :: Coordinate -> Coordinate -> Bool
Eq, Int -> Coordinate -> ShowS
[Coordinate] -> ShowS
Coordinate -> String
(Int -> Coordinate -> ShowS)
-> (Coordinate -> String)
-> ([Coordinate] -> ShowS)
-> Show Coordinate
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Coordinate -> ShowS
showsPrec :: Int -> Coordinate -> ShowS
$cshow :: Coordinate -> String
show :: Coordinate -> String
$cshowList :: [Coordinate] -> ShowS
showList :: [Coordinate] -> ShowS
Show, (forall x. Coordinate -> Rep Coordinate x)
-> (forall x. Rep Coordinate x -> Coordinate) -> Generic Coordinate
forall x. Rep Coordinate x -> Coordinate
forall x. Coordinate -> Rep Coordinate x
forall a.
(forall x. a -> Rep a x) -> (forall x. Rep a x -> a) -> Generic a
$cfrom :: forall x. Coordinate -> Rep Coordinate x
from :: forall x. Coordinate -> Rep Coordinate x
$cto :: forall x. Rep Coordinate x -> Coordinate
to :: forall x. Rep Coordinate x -> Coordinate
Generic)

-- | Treat two coordinates like vectors and add them together.
--
-- >>> Coordinate (2, 1, 0, 0) `plus` Coordinate (-1, 1, 0, 0)
-- Coordinate (1.0,2.0,0.0,0.0)
plus :: Coordinate -> Coordinate -> Coordinate
plus :: Coordinate -> Coordinate -> Coordinate
plus = (Double -> Double -> Double)
-> Coordinate -> Coordinate -> Coordinate
combine Double -> Double -> Double
forall a. Num a => a -> a -> a
(+)

-- | Treat two coordinates like vectors and subtract the second one from the first.
--
-- >>> Coordinate (2, 1, 0, 0) `minus` Coordinate (-1, 1, 0, 0)
-- Coordinate (3.0,0.0,0.0,0.0)
minus :: Coordinate -> Coordinate -> Coordinate
minus :: Coordinate -> Coordinate -> Coordinate
minus = (Double -> Double -> Double)
-> Coordinate -> Coordinate -> Coordinate
combine (-)

-- | Scale each component of a coordinate with the given function.
--
-- >>> scale (*2) $ Coordinate (1, 2, 0, 0)
-- Coordinate (2.0,4.0,0.0,0.0)
scale :: (Double -> Double) -> Coordinate -> Coordinate
scale :: (Double -> Double) -> Coordinate -> Coordinate
scale Double -> Double
f (Coordinate (Double
t, Double
x, Double
y, Double
z)) = (Double, Double, Double, Double) -> Coordinate
Coordinate (Double -> Double
f Double
t, Double -> Double
f Double
x, Double -> Double
f Double
y, Double -> Double
f Double
z)

-- | Combine the corresponding components of two coordinates with the given function.
--
-- >>> combine (*) (Coordinate (2, 3, 1, 2)) (Coordinate (1, -1, 2, 3))
-- Coordinate (2.0,-3.0,2.0,6.0)
combine :: (Double -> Double -> Double) -> Coordinate -> Coordinate -> Coordinate
combine :: (Double -> Double -> Double)
-> Coordinate -> Coordinate -> Coordinate
combine Double -> Double -> Double
f (Coordinate (Double
t, Double
x, Double
y, Double
z)) (Coordinate (Double
t', Double
x', Double
y', Double
z')) =
  (Double, Double, Double, Double) -> Coordinate
Coordinate (Double -> Double -> Double
f Double
t Double
t', Double -> Double -> Double
f Double
x Double
x', Double -> Double -> Double
f Double
y Double
y', Double -> Double -> Double
f Double
z Double
z')

-- | Returns whether two coordinates are approximately equal to each other.
--
-- >>> Coordinate (1, 1, 0, 0) `approximates` Coordinate (1, 1, 0, 0)
-- True
-- >>> Coordinate (1, -1, 0, 0) `approximates` Coordinate (1, 1, 0, 0)
-- False
--
-- The threshold for being approximately equivalent is somewhat arbitrary.
approximates :: Coordinate -> Coordinate -> Bool
approximates :: Coordinate -> Coordinate -> Bool
approximates (Coordinate (Double
t, Double
x, Double
y, Double
z)) (Coordinate (Double
t', Double
x', Double
y', Double
z')) =
  Double
t Double -> Double -> Bool
forall {a}. (Ord a, Fractional a) => a -> a -> Bool
`closeTo` Double
t' Bool -> Bool -> Bool
&& Double
x Double -> Double -> Bool
forall {a}. (Ord a, Fractional a) => a -> a -> Bool
`closeTo` Double
x' Bool -> Bool -> Bool
&& Double
y Double -> Double -> Bool
forall {a}. (Ord a, Fractional a) => a -> a -> Bool
`closeTo` Double
y' Bool -> Bool -> Bool
&& Double
z Double -> Double -> Bool
forall {a}. (Ord a, Fractional a) => a -> a -> Bool
`closeTo` Double
z'
  where
    closeTo :: a -> a -> Bool
closeTo a
a a
b = a -> a
forall a. Num a => a -> a
abs (a
a a -> a -> a
forall a. Num a => a -> a -> a
- a
b) a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
1e-5

-- | The spacetime interval between the origin and the given coordinates.
--
-- This uses the convention where \(ds^2 = dt^2 - dx^2 - dy^2 - dz^2\).
--
-- >>> interval $ Coordinate (2, 1, 0, 0)
-- 3.0
interval :: Coordinate -> Double
interval :: Coordinate -> Double
interval (Coordinate (Double
t, Double
x, Double
y, Double
z)) = Double
t Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
y Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
z Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Double
2

-- | Apply the Lorentz transformation.
--
-- This transforms coordinates in an original inertial frame to those
-- in a new inertial frame.  It is assumed that the origin coincides
-- between the two frames, and that the velocity only has an \(x\) component.
--
-- >>> let expected = Coordinate (0.414213562,0.414213562,2.0,3.0)
-- >>> let actual = transform (sqrt 0.5) $ Coordinate (1, 1, 2, 3)
-- >>> actual `approximates` expected
-- True
--
-- For reference, look up the
-- [Lorentz transformation](https://chungyc.org/article/reference/physics/relativity/#lorentz).
transform ::
  -- | Velocity of the new frame in the original frame.
  Double ->
  -- | Coordinate in the original frame.
  Coordinate ->
  -- | Coordinate in the new frame.
  Coordinate
transform :: Double -> Coordinate -> Coordinate
transform Double
v (Coordinate (Double
t, Double
x, Double
y, Double
z)) = (Double, Double, Double, Double) -> Coordinate
Coordinate (Double
t', Double
x', Double
y', Double
z')
  where
    gamma :: Double
gamma = Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
v Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Double
2)
    t' :: Double
t' = Double
gamma Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x)
    x' :: Double
x' = Double
gamma Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t)
    y' :: Double
y' = Double
y
    z' :: Double
z' = Double
z

-- | Returns a Lorentz transformation.
--
-- It will transform diagrams in an original inertial frame to those
-- in a new inertial frame.  It is assumed that the origin coincides
-- between the two frames, and that the velocity only has an \(x\) component.
--
-- >>> import Diagrams.Prelude qualified as D
-- >>> let c = Coordinate (2, 1, 0, 0)
-- >>> let v = 0.5
-- >>> let expected = vectorize2D $ transform v c
-- >>> let actual = D.transform (transformation v) $ vectorize2D c
-- >>> abs (actual - expected) < 1e-5
-- True
--
-- For reference, look up the
-- [Lorentz transformation](https://chungyc.org/article/reference/physics/relativity/#lorentz).
transformation ::
  -- | Velocity of the new fram in the original frame.
  Double ->
  -- | Transformation for a two-dimensional "Diagrams" object.
  Transformation V2 Double
transformation :: Double -> Transformation V2 Double
transformation Double
v = M22 Double -> V2 Double -> Transformation V2 Double
forall n. Floating n => M22 n -> V2 n -> T2 n
fromMat22 M22 Double
matrix (Double -> Double -> V2 Double
forall a. a -> a -> V2 a
V2 Double
0 Double
0)
  where
    gamma :: Double
gamma = Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
v Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Double
2)
    matrix :: M22 Double
matrix = V2 Double -> V2 Double -> M22 Double
forall a. a -> a -> V2 a
V2 (Double -> Double -> V2 Double
forall a. a -> a -> V2 a
V2 Double
gamma (-Double
gamma Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v)) (Double -> Double -> V2 Double
forall a. a -> a -> V2 a
V2 (-Double
gamma Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v) Double
gamma)

-- | Convert given coordinates into a two-dimensional "Diagrams" vector.
--
-- It will project the \(x\) coordinate to the \(x\) component
-- and the \(t\) coordinate to the \(y\) component of the vector.
-- The other coordinates will be ignored.
--
-- >>> vectorize2D $ Coordinate (1, 1, 0, 0)
-- V2 1.0 1.0
vectorize2D :: Coordinate -> V2 Double
vectorize2D :: Coordinate -> V2 Double
vectorize2D (Coordinate (Double
t, Double
x, Double
_, Double
_)) = (Double, Double) -> V2 Double
forall n. (n, n) -> V2 n
r2 (Double
x, Double
t)

-- | Stores various options for drawing axes.
data AxesOptions = AxesOptions
  { -- | Length of each axis.  The default is 1.
    --
    -- >>> axesOptions { axesLength = 4 }
    -- AxesOptions {axesLength = 4.0, axesLightcone = True}
    AxesOptions -> Double
axesLength :: Double,
    -- | Whether to draw the lightcone.  It is drawn by default.
    --
    -- >>> axesOptions { axesLightcone = False }
    -- AxesOptions {axesLength = 1.0, axesLightcone = False}
    AxesOptions -> Bool
axesLightcone :: Bool
  }
  deriving (Int -> AxesOptions -> ShowS
[AxesOptions] -> ShowS
AxesOptions -> String
(Int -> AxesOptions -> ShowS)
-> (AxesOptions -> String)
-> ([AxesOptions] -> ShowS)
-> Show AxesOptions
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> AxesOptions -> ShowS
showsPrec :: Int -> AxesOptions -> ShowS
$cshow :: AxesOptions -> String
show :: AxesOptions -> String
$cshowList :: [AxesOptions] -> ShowS
showList :: [AxesOptions] -> ShowS
Show)

-- | Default options for drawing axes.
axesOptions :: AxesOptions
axesOptions :: AxesOptions
axesOptions = AxesOptions {axesLength :: Double
axesLength = Double
1.0, axesLightcone :: Bool
axesLightcone = Bool
True}

-- | Draw axes for a spacetime diagram with default options.
--
-- Both axes will extend up to length 1 from the origin,
-- and the lightcone will also be drawn.
-- Units will be in Planck units, i.e., \(c = 1\).
--
-- >>> let _ = axes
axes :: Diagram B
axes :: Diagram B
axes = AxesOptions -> Diagram B
axesWith AxesOptions
axesOptions

-- | Draw axes for a spacetime diagram with given options.
--
-- The horizontal axis will be for the \(x\) space coordinate,
-- while the vertical axis will be for the \(t\) space coordinate.
-- Units will be in Planck units, i.e., \(c = 1\).
--
-- >>> let _ = axesWith axesOptions { axesLength = 8 }
axesWith :: AxesOptions -> Diagram B
axesWith :: AxesOptions -> Diagram B
axesWith AxesOptions {Double
axesLength :: AxesOptions -> Double
axesLength :: Double
axesLength, Bool
axesLightcone :: AxesOptions -> Bool
axesLightcone :: Bool
axesLightcone} =
  QDiagram B V2 Double Any
xAxis QDiagram B V2 Double Any
-> QDiagram B V2 Double Any -> QDiagram B V2 Double Any
forall a. Semigroup a => a -> a -> a
<> QDiagram B V2 Double Any
tAxis QDiagram B V2 Double Any
-> QDiagram B V2 Double Any -> QDiagram B V2 Double Any
forall a. Semigroup a => a -> a -> a
<> QDiagram B V2 Double Any
lightCone
  where
    xAxis :: QDiagram B V2 Double Any
xAxis = (-Double
PrevDim
  (Point
     (V (QDiagram B V2 Double Any)) (N (QDiagram B V2 Double Any)))
axesLength PrevDim
  (Point
     (V (QDiagram B V2 Double Any)) (N (QDiagram B V2 Double Any)))
-> FinalCoord
     (Point
        (V (QDiagram B V2 Double Any)) (N (QDiagram B V2 Double Any)))
-> Point
     (V (QDiagram B V2 Double Any)) (N (QDiagram B V2 Double Any))
forall c. Coordinates c => PrevDim c -> FinalCoord c -> c
^& FinalCoord
  (V (QDiagram B V2 Double Any) (N (QDiagram B V2 Double Any)))
FinalCoord
  (Point
     (V (QDiagram B V2 Double Any)) (N (QDiagram B V2 Double Any)))
0) Point (V (QDiagram B V2 Double Any)) (N (QDiagram B V2 Double Any))
-> Point
     (V (QDiagram B V2 Double Any)) (N (QDiagram B V2 Double Any))
-> QDiagram B V2 Double Any
forall t (v :: * -> *) n.
(V t ~ v, N t ~ n, TrailLike t) =>
Point v n -> Point v n -> t
~~ (Double
PrevDim
  (Point
     (V (QDiagram B V2 Double Any)) (N (QDiagram B V2 Double Any)))
axesLength PrevDim
  (Point
     (V (QDiagram B V2 Double Any)) (N (QDiagram B V2 Double Any)))
-> FinalCoord
     (Point
        (V (QDiagram B V2 Double Any)) (N (QDiagram B V2 Double Any)))
-> Point
     (V (QDiagram B V2 Double Any)) (N (QDiagram B V2 Double Any))
forall c. Coordinates c => PrevDim c -> FinalCoord c -> c
^& FinalCoord
  (V (QDiagram B V2 Double Any) (N (QDiagram B V2 Double Any)))
FinalCoord
  (Point
     (V (QDiagram B V2 Double Any)) (N (QDiagram B V2 Double Any)))
0)
    tAxis :: QDiagram B V2 Double Any
tAxis = (Double
PrevDim (Point V2 Double)
0 PrevDim (Point V2 Double)
-> FinalCoord (Point V2 Double) -> Point V2 Double
forall c. Coordinates c => PrevDim c -> FinalCoord c -> c
^& (-Double
axesLength)) Point V2 Double -> Point V2 Double -> QDiagram B V2 Double Any
forall n b.
(TypeableFloat n, Renderable (Path V2 n) b) =>
Point V2 n -> Point V2 n -> QDiagram b V2 n Any
`arrowBetween` (Double
PrevDim (Point V2 Double)
0 PrevDim (Point V2 Double)
-> FinalCoord (Point V2 Double) -> Point V2 Double
forall c. Coordinates c => PrevDim c -> FinalCoord c -> c
^& Double
FinalCoord (Point V2 Double)
axesLength)
    lightCone :: QDiagram B V2 Double Any
lightCone
      | Bool
axesLightcone =
          (QDiagram B V2 Double Any
forwardEdge QDiagram B V2 Double Any
-> QDiagram B V2 Double Any -> QDiagram B V2 Double Any
forall a. Semigroup a => a -> a -> a
<> QDiagram B V2 Double Any
backwardEdge)
            # opacity 0.75
            # dashingN [0.01, 0.02] 0
            # lw ultraThin
      | Bool
otherwise = QDiagram B V2 Double Any
forall a. Monoid a => a
mempty
      where
        forwardEdge :: QDiagram B V2 Double Any
forwardEdge = Point (V (QDiagram B V2 Double Any)) (N (QDiagram B V2 Double Any))
lowerLeftCorner Point (V (QDiagram B V2 Double Any)) (N (QDiagram B V2 Double Any))
-> Point
     (V (QDiagram B V2 Double Any)) (N (QDiagram B V2 Double Any))
-> QDiagram B V2 Double Any
forall t (v :: * -> *) n.
(V t ~ v, N t ~ n, TrailLike t) =>
Point v n -> Point v n -> t
~~ Point (V (QDiagram B V2 Double Any)) (N (QDiagram B V2 Double Any))
upperRightCorner
        backwardEdge :: QDiagram B V2 Double Any
backwardEdge = QDiagram B V2 Double Any -> QDiagram B V2 Double Any
forall (v :: * -> *) n t.
(InSpace v n t, R1 v, Transformable t) =>
t -> t
reflectX QDiagram B V2 Double Any
forwardEdge
        lowerLeftCorner :: Point (V (QDiagram B V2 Double Any)) (N (QDiagram B V2 Double Any))
lowerLeftCorner = (-Double
axesLength) PrevDim
  (Point
     (V (QDiagram B V2 Double Any)) (N (QDiagram B V2 Double Any)))
-> FinalCoord
     (Point
        (V (QDiagram B V2 Double Any)) (N (QDiagram B V2 Double Any)))
-> Point
     (V (QDiagram B V2 Double Any)) (N (QDiagram B V2 Double Any))
forall c. Coordinates c => PrevDim c -> FinalCoord c -> c
^& (-Double
axesLength)
        upperRightCorner :: Point (V (QDiagram B V2 Double Any)) (N (QDiagram B V2 Double Any))
upperRightCorner = Double
PrevDim
  (Point
     (V (QDiagram B V2 Double Any)) (N (QDiagram B V2 Double Any)))
axesLength PrevDim
  (Point
     (V (QDiagram B V2 Double Any)) (N (QDiagram B V2 Double Any)))
-> FinalCoord
     (Point
        (V (QDiagram B V2 Double Any)) (N (QDiagram B V2 Double Any)))
-> Point
     (V (QDiagram B V2 Double Any)) (N (QDiagram B V2 Double Any))
forall c. Coordinates c => PrevDim c -> FinalCoord c -> c
^& Double
FinalCoord
  (Point
     (V (QDiagram B V2 Double Any)) (N (QDiagram B V2 Double Any)))
axesLength