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(..))
data AstronomyEpoch = B1900
| B1950
| J2000
| J2050
deriving (Show, Eq)
epochToJD :: AstronomyEpoch -> JulianDate
epochToJD B1900 = b1900
epochToJD B1950 = b1950
epochToJD J2000 = j2000
epochToJD J2050 = j2050
data PrecessionalConstants = PrecessionalConstants {
pcM :: Double
, pcN :: Double
, pcN' :: Double
}
precessionalConstants :: AstronomyEpoch -> PrecessionalConstants
precessionalConstants B1900 = PrecessionalConstants 3.07234 1.33645 20.0468
precessionalConstants B1950 = PrecessionalConstants 3.07327 1.33617 20.0426
precessionalConstants J2000 = PrecessionalConstants 3.07420 1.33589 20.0383
precessionalConstants J2050 = PrecessionalConstants 3.07513 1.33560 20.0340
precession1 :: AstronomyEpoch -> EquatorialCoordinates1 -> JulianDate -> EquatorialCoordinates1
precession1 epoch (EC1 delta alpha) jd =
let delta' = toRadians delta
alpha' = toRadians $ fromDecimalHours alpha
years = numberOfYears (epochToJD epoch) jd
PrecessionalConstants m n n' = precessionalConstants epoch
s1 = DH $ (m + n*(sin alpha')*(tan delta'))*years / 3600
s2 = DD $ (n'*(cos alpha')) * years / 3600
in (EC1 (delta + s2) (alpha + s1))
precession2 :: JulianDate -> EquatorialCoordinates1 -> JulianDate -> EquatorialCoordinates1
precession2 epoch ec jd =
let p' = prepareMatrixP' $ numberOfCenturies j2000 epoch
v = prepareColumnVectorV ec
p = transpose $ prepareMatrixP' $ numberOfCenturies j2000 jd
[m, n, k] = toList $ p*(p'*v)
alpha = atan2 n m
delta = asin k
in EC1 (fromRadians delta) (toDecimalHours $ fromRadians alpha)
prepareMatrixP' n =
let x = U.toRadians $ 0.6406161*n + 0.0000839*n*n + 0.0000050*n*n*n
z = U.toRadians $ 0.6406161*n + 0.0003041*n*n + 0.0000051*n*n*n
t = U.toRadians $ 0.5567530*n - 0.0001185*n*n - 0.0000116*n*n*n
cx = cos x
sx = sin x
cz = cos z
sz = sin z
ct = cos t
st = sin t
matrix = [ [cx*ct*cz-sx*sz, cx*ct*sz+sx*cz, cx*st]
, [(-sx)*ct*cz-cx*sz, (-sx)*ct*sz+cx*cz, (-sx)*st]
, [(-st)*cz, (-st)*sz, ct] ]
in fromLists matrix
prepareColumnVectorV (EC1 delta alpha) =
let d = toRadians delta
a = toRadians $ fromDecimalHours alpha
cd = cos d
sd = sin d
ca = cos a
sa = sin a
v = [ca*cd, sa*cd, sd]
in fromList 3 1 v