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 (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)
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
data PrecessionalConstants = PrecessionalConstants {
PrecessionalConstants -> Double
pcM :: Double
, PrecessionalConstants -> Double
pcN :: Double
, PrecessionalConstants -> Double
pcN' :: Double
}
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
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))
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