{-|
Module: Data.Astro.Effects.Precession
Description: Luni-solar precession
Copyright: Alexander Ignatyev, 2016

Luni-solar precession.
-}

module Data.Astro.Effects.Precession
(
  AstronomyEpoch(..)
  , precession1
  , precession2
)

where

import Data.Matrix

import qualified Data.Astro.Utils as U
import Data.Astro.Types (DecimalDegrees(..), DecimalHours(..), toDecimalHours, fromDecimalHours, toRadians, fromRadians)
import Data.Astro.Time.JulianDate (JulianDate(..), numberOfYears, numberOfCenturies)
import Data.Astro.Time.Epoch (b1900, b1950, j2000, j2050)
import Data.Astro.Coordinate (EquatorialCoordinates1(..))


-------------------------------------------------------------------------------

-- Low-precision Precession


-- | Epoch Enumeration. See also "Data.Astro.Time.JulianDate" module.

data AstronomyEpoch = B1900  -- ^ Epoch B1900.0

                    | B1950  -- ^ Epoch B1950.0

                    | J2000  -- ^ Epoch J2000.0

                    | J2050  -- ^ Epoch J2050.0

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


-- | Get the start date of the specified Epoch.

epochToJD :: AstronomyEpoch -> JulianDate
epochToJD :: AstronomyEpoch -> JulianDate
epochToJD AstronomyEpoch
B1900 = JulianDate
b1900
epochToJD AstronomyEpoch
B1950 = JulianDate
b1950
epochToJD AstronomyEpoch
J2000 = JulianDate
j2000
epochToJD AstronomyEpoch
J2050 = JulianDate
j2050


-- | Precisional Constants

data PrecessionalConstants = PrecessionalConstants {
  PrecessionalConstants -> Double
pcM :: Double     -- ^ seconds

  , PrecessionalConstants -> Double
pcN :: Double   -- ^ seconds

  , PrecessionalConstants -> Double
pcN' :: Double  -- ^ arcsec

  }


-- | Get Precision Constants of the Epoch

precessionalConstants :: AstronomyEpoch -> PrecessionalConstants
precessionalConstants :: AstronomyEpoch -> PrecessionalConstants
precessionalConstants AstronomyEpoch
B1900 = Double -> Double -> Double -> PrecessionalConstants
PrecessionalConstants Double
3.07234 Double
1.33645 Double
20.0468
precessionalConstants AstronomyEpoch
B1950 = Double -> Double -> Double -> PrecessionalConstants
PrecessionalConstants Double
3.07327 Double
1.33617 Double
20.0426
precessionalConstants AstronomyEpoch
J2000 = Double -> Double -> Double -> PrecessionalConstants
PrecessionalConstants Double
3.07420 Double
1.33589 Double
20.0383
precessionalConstants AstronomyEpoch
J2050 = Double -> Double -> Double -> PrecessionalConstants
PrecessionalConstants Double
3.07513 Double
1.33560 Double
20.0340


-- | Low-precision method to calculate luni-solar precession.

-- It takes Epoch, Equatorial Coordinates those correct at the given epoch, Julian Date of the observation.

-- It returns corrected Equatorial Coordinates.

precession1 :: AstronomyEpoch -> EquatorialCoordinates1 -> JulianDate -> EquatorialCoordinates1
precession1 :: AstronomyEpoch
-> EquatorialCoordinates1 -> JulianDate -> EquatorialCoordinates1
precession1 AstronomyEpoch
epoch (EC1 DecimalDegrees
delta DecimalHours
alpha) JulianDate
jd =
  let delta' :: Double
delta' = DecimalDegrees -> Double
toRadians DecimalDegrees
delta
      alpha' :: Double
alpha' = DecimalDegrees -> Double
toRadians (DecimalDegrees -> Double) -> DecimalDegrees -> Double
forall a b. (a -> b) -> a -> b
$ DecimalHours -> DecimalDegrees
fromDecimalHours DecimalHours
alpha
      years :: Double
years = JulianDate -> JulianDate -> Double
numberOfYears (AstronomyEpoch -> JulianDate
epochToJD AstronomyEpoch
epoch) JulianDate
jd
      PrecessionalConstants Double
m Double
n Double
n' = AstronomyEpoch -> PrecessionalConstants
precessionalConstants AstronomyEpoch
epoch
      s1 :: DecimalHours
s1 = Double -> DecimalHours
DH (Double -> DecimalHours) -> Double -> DecimalHours
forall a b. (a -> b) -> a -> b
$ (Double
m Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
nDouble -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
sin Double
alpha')Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
tan Double
delta'))Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
years Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3600
      s2 :: DecimalDegrees
s2 = Double -> DecimalDegrees
DD (Double -> DecimalDegrees) -> Double -> DecimalDegrees
forall a b. (a -> b) -> a -> b
$ (Double
n'Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
cos Double
alpha')) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
years Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3600
  in (DecimalDegrees -> DecimalHours -> EquatorialCoordinates1
EC1 (DecimalDegrees
delta DecimalDegrees -> DecimalDegrees -> DecimalDegrees
forall a. Num a => a -> a -> a
+ DecimalDegrees
s2) (DecimalHours
alpha DecimalHours -> DecimalHours -> DecimalHours
forall a. Num a => a -> a -> a
+ DecimalHours
s1))


-------------------------------------------------------------------------------

-- Rigorous Method



-- | Rigorous method to calculate luni-solar precession.

-- It takes julian date at whose the coordinates are correct, Equatorial Coordinates, Julian Date of the observation.

-- It returns corrected Equatorial Coordinates.

precession2 :: JulianDate -> EquatorialCoordinates1 -> JulianDate -> EquatorialCoordinates1
precession2 :: JulianDate
-> EquatorialCoordinates1 -> JulianDate -> EquatorialCoordinates1
precession2 JulianDate
epoch EquatorialCoordinates1
ec JulianDate
jd =
  let p' :: Matrix Double
p' = Double -> Matrix Double
forall a. Floating a => a -> Matrix a
prepareMatrixP' (Double -> Matrix Double) -> Double -> Matrix Double
forall a b. (a -> b) -> a -> b
$ JulianDate -> JulianDate -> Double
numberOfCenturies JulianDate
j2000 JulianDate
epoch
      v :: Matrix Double
v = EquatorialCoordinates1 -> Matrix Double
prepareColumnVectorV EquatorialCoordinates1
ec
      p :: Matrix Double
p = Matrix Double -> Matrix Double
forall a. Matrix a -> Matrix a
transpose (Matrix Double -> Matrix Double) -> Matrix Double -> Matrix Double
forall a b. (a -> b) -> a -> b
$ Double -> Matrix Double
forall a. Floating a => a -> Matrix a
prepareMatrixP' (Double -> Matrix Double) -> Double -> Matrix Double
forall a b. (a -> b) -> a -> b
$ JulianDate -> JulianDate -> Double
numberOfCenturies JulianDate
j2000 JulianDate
jd
      [Double
m, Double
n, Double
k] = Matrix Double -> [Double]
forall a. Matrix a -> [a]
toList (Matrix Double -> [Double]) -> Matrix Double -> [Double]
forall a b. (a -> b) -> a -> b
$ Matrix Double
pMatrix Double -> Matrix Double -> Matrix Double
forall a. Num a => a -> a -> a
*(Matrix Double
p'Matrix Double -> Matrix Double -> Matrix Double
forall a. Num a => a -> a -> a
*Matrix Double
v)
      alpha :: Double
alpha = Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 Double
n Double
m
      delta :: Double
delta = Double -> Double
forall a. Floating a => a -> a
asin Double
k
  in DecimalDegrees -> DecimalHours -> EquatorialCoordinates1
EC1 (Double -> DecimalDegrees
fromRadians Double
delta) (DecimalDegrees -> DecimalHours
toDecimalHours (DecimalDegrees -> DecimalHours) -> DecimalDegrees -> DecimalHours
forall a b. (a -> b) -> a -> b
$ Double -> DecimalDegrees
fromRadians Double
alpha)


prepareMatrixP' :: a -> Matrix a
prepareMatrixP' a
n =
  let x :: a
x = a -> a
forall a. Floating a => a -> a
U.toRadians (a -> a) -> a -> a
forall a b. (a -> b) -> a -> b
$ a
0.6406161a -> a -> a
forall a. Num a => a -> a -> a
*a
n a -> a -> a
forall a. Num a => a -> a -> a
+ a
0.0000839a -> a -> a
forall a. Num a => a -> a -> a
*a
na -> a -> a
forall a. Num a => a -> a -> a
*a
n a -> a -> a
forall a. Num a => a -> a -> a
+ a
0.0000050a -> a -> a
forall a. Num a => a -> a -> a
*a
na -> a -> a
forall a. Num a => a -> a -> a
*a
na -> a -> a
forall a. Num a => a -> a -> a
*a
n
      z :: a
z = a -> a
forall a. Floating a => a -> a
U.toRadians (a -> a) -> a -> a
forall a b. (a -> b) -> a -> b
$ a
0.6406161a -> a -> a
forall a. Num a => a -> a -> a
*a
n a -> a -> a
forall a. Num a => a -> a -> a
+ a
0.0003041a -> a -> a
forall a. Num a => a -> a -> a
*a
na -> a -> a
forall a. Num a => a -> a -> a
*a
n a -> a -> a
forall a. Num a => a -> a -> a
+ a
0.0000051a -> a -> a
forall a. Num a => a -> a -> a
*a
na -> a -> a
forall a. Num a => a -> a -> a
*a
na -> a -> a
forall a. Num a => a -> a -> a
*a
n
      t :: a
t = a -> a
forall a. Floating a => a -> a
U.toRadians (a -> a) -> a -> a
forall a b. (a -> b) -> a -> b
$ a
0.5567530a -> a -> a
forall a. Num a => a -> a -> a
*a
n a -> a -> a
forall a. Num a => a -> a -> a
- a
0.0001185a -> a -> a
forall a. Num a => a -> a -> a
*a
na -> a -> a
forall a. Num a => a -> a -> a
*a
n a -> a -> a
forall a. Num a => a -> a -> a
- a
0.0000116a -> a -> a
forall a. Num a => a -> a -> a
*a
na -> a -> a
forall a. Num a => a -> a -> a
*a
na -> a -> a
forall a. Num a => a -> a -> a
*a
n
      cx :: a
cx = a -> a
forall a. Floating a => a -> a
cos a
x
      sx :: a
sx = a -> a
forall a. Floating a => a -> a
sin a
x
      cz :: a
cz = a -> a
forall a. Floating a => a -> a
cos a
z
      sz :: a
sz = a -> a
forall a. Floating a => a -> a
sin a
z
      ct :: a
ct = a -> a
forall a. Floating a => a -> a
cos a
t
      st :: a
st = a -> a
forall a. Floating a => a -> a
sin a
t
      matrix :: [[a]]
matrix = [ [a
cxa -> a -> a
forall a. Num a => a -> a -> a
*a
cta -> a -> a
forall a. Num a => a -> a -> a
*a
cza -> a -> a
forall a. Num a => a -> a -> a
-a
sxa -> a -> a
forall a. Num a => a -> a -> a
*a
sz,    a
cxa -> a -> a
forall a. Num a => a -> a -> a
*a
cta -> a -> a
forall a. Num a => a -> a -> a
*a
sza -> a -> a
forall a. Num a => a -> a -> a
+a
sxa -> a -> a
forall a. Num a => a -> a -> a
*a
cz,    a
cxa -> a -> a
forall a. Num a => a -> a -> a
*a
st]
               , [(-a
sx)a -> a -> a
forall a. Num a => a -> a -> a
*a
cta -> a -> a
forall a. Num a => a -> a -> a
*a
cza -> a -> a
forall a. Num a => a -> a -> a
-a
cxa -> a -> a
forall a. Num a => a -> a -> a
*a
sz, (-a
sx)a -> a -> a
forall a. Num a => a -> a -> a
*a
cta -> a -> a
forall a. Num a => a -> a -> a
*a
sza -> a -> a
forall a. Num a => a -> a -> a
+a
cxa -> a -> a
forall a. Num a => a -> a -> a
*a
cz, (-a
sx)a -> a -> a
forall a. Num a => a -> a -> a
*a
st]
               , [(-a
st)a -> a -> a
forall a. Num a => a -> a -> a
*a
cz,          (-a
st)a -> a -> a
forall a. Num a => a -> a -> a
*a
sz,          a
ct] ]
  in [[a]] -> Matrix a
forall a. [[a]] -> Matrix a
fromLists [[a]]
matrix

prepareColumnVectorV :: EquatorialCoordinates1 -> Matrix Double
prepareColumnVectorV (EC1 DecimalDegrees
delta DecimalHours
alpha) =
  let d :: Double
d = DecimalDegrees -> Double
toRadians DecimalDegrees
delta
      a :: Double
a = DecimalDegrees -> Double
toRadians (DecimalDegrees -> Double) -> DecimalDegrees -> Double
forall a b. (a -> b) -> a -> b
$ DecimalHours -> DecimalDegrees
fromDecimalHours DecimalHours
alpha
      cd :: Double
cd = Double -> Double
forall a. Floating a => a -> a
cos Double
d
      sd :: Double
sd = Double -> Double
forall a. Floating a => a -> a
sin Double
d
      ca :: Double
ca = Double -> Double
forall a. Floating a => a -> a
cos Double
a
      sa :: Double
sa = Double -> Double
forall a. Floating a => a -> a
sin Double
a
      v :: [Double]
v = [Double
caDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
cd, Double
saDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
cd, Double
sd]
  in Int -> Int -> [Double] -> Matrix Double
forall a. Int -> Int -> [a] -> Matrix a
fromList Int
3 Int
1 [Double]
v