{-|
Module: Data.Astro.Coordinate
Description: Celestial Coordinate Systems
Copyright: Alexander Ignatyev, 2016

See "Data.Astro.Types" module for Georgraphic Coordinates.

= Celestial Coordinate Systems

== /Horizon coordinates/

* __altitude, α__ - /'how far up'/ angle from the horizontal plane in degrees
* __azimuth,  Α__ - /'how far round'/ agle from the north direction in degrees to the east


== /Equatorial coordinates/

Accoring to the equatorial coordinates system stars move westwards along the circles centered in the north selestial pole,
making the full cicrle in 24 hours of sidereal time (see "Data.Astro.Time.Sidereal").

* __declination, δ__ - /'how far up'/ angle from the quatorial plane;
* __right ascension, α__  - /'how far round'/ angle from the /vernal equinox/ to the east; __/or/__
* __hour angle__ - /'how far round'/ angle from the meridian to the west


== /Ecliptic Coordinate/

Accoring to the ecliptic coordinates system the Sun moves eastwards along the trace of th ecliptic. The Sun's ecplitic latitude is always 0.

* __ecliptic latitude, β__ - /'how far up'/ angle from the ecliptic
* __ecliptic longitude, λ__ - /'how far round'/ angle from the /vernal equinox/ to the east


== /Galactic Coordinates/

* __galactic latitute, b__ - /'how far up'/ angle from the plane of the Galaxy
* __galactiv longitude, l__ - - /'how far round'/ angle from the direction the Sun - the centre of the Galaxy


== /Terms/

* __ecliptic__ - the plane containing the Earth's orbit around the Sun
* __vernal equinox__, ♈ - fixed direction lies along the line of the intersection of the equatorial plane and the ecliptic
* __obliquity of the ecliptic, β__ - the angle between the plane of the Earth's equator and the ecliptic
* __north selestial pole, P__ - the point on the selestial sphere, right above the Earth's North Pole


= Examples

== /Horizontal Coordinate System/
@
import Data.Astro.Coordinate
import Data.Astro.Types

hc :: HorizonCoordinates
hc = HC (DD 30.5) (DD 180)
-- HC {hAltitude = DD 30.0, hAzimuth = DD 180.0}
@

== /Equatorial Coordinate System/
@
import Data.Astro.Coordinate
import Data.Astro.Types

ec1 :: EquatorialCoordinates1
ec1 = EC1 (DD 71.7) (DH 8)
-- EC1 {e1Declination = DD 71.7, e1RightAscension = DH 8.0}

ec2 :: EquatorialCoordinates2
ec2 = EC1 (DD 77.7) (DH 11)
-- EC2 {e2Declination = DD 77.7, e2HoursAngle = DH 11.0}
@

== /Transformations/
@
import Data.Astro.Time.JulianDate
import Data.Astro.Coordinate
import Data.Astro.Types

ro :: GeographicCoordinates
ro = GeoC (fromDMS 51 28 40) (-(fromDMS 0 0 5))

dt :: LocalCivilTime
dt = lctFromYMDHMS (DH 1) 2017 6 25 10 29 0

sunHC :: HorizonCoordinates
sunHC = HC (fromDMS 49 18 21.77) (fromDMS 118 55 19.53)
-- HC {hAltitude = DD 49.30604722222222, hAzimuth = DD 118.92209166666666}

sunEC2 :: EquatorialCoordinates2
sunEC2 = horizonToEquatorial (geoLatitude ro) sunHC
-- EC2 {e2Declination = DD 23.378295912623855, e2HoursAngle = DH 21.437117068873537}

sunEC1 :: EquatorialCoordinates1
sunEC1 = EC1 (e2Declination sunEC2) (haToRA (e2HoursAngle sunEC2) (geoLongitude ro) (lctUniversalTime dt))
-- EC1 {e1Declination = DD 23.378295912623855, e1RightAscension = DH 6.29383725890224}


sunEC2' :: EquatorialCoordinates2
sunEC2' = EC2 (e1Declination sunEC1) (raToHA (e1RightAscension sunEC1) (geoLongitude ro) (lctUniversalTime dt))
-- EC2 {e2Declination = DD 23.378295912623855, e2HoursAngle = DH 21.437117068873537}

sunHC' :: HorizonCoordinates
sunHC' = equatorialToHorizon (geoLatitude ro) sunEC2'
-- HC {hAltitude = DD 49.30604722222222, hAzimuth = DD 118.92209166666666}
@

=== /Function-shortcuts/

@
import Data.Astro.Time.JulianDate
import Data.Astro.Coordinate
import Data.Astro.Types

ro :: GeographicCoordinates
ro = GeoC (fromDMS 51 28 40) (-(fromDMS 0 0 5))

dt :: LocalCivilTime
dt = lctFromYMDHMS (DH 1) 2017 6 25 10 29 0

sunHC :: HorizonCoordinates
sunHC = HC (fromDMS 49 18 21.77) (fromDMS 118 55 19.53)
-- HC {hAltitude = DD 49.30604722222222, hAzimuth = DD 118.92209166666666}

sunEC1 :: EquatorialCoordinates1
sunEC1 = hcToEC1 ro (lctUniversalTime dt) sunHC
-- EC1 {e1Declination = DD 23.378295912623855, e1RightAscension = DH 6.29383725890224}

sunHC' :: HorizonCoordinates
sunHC' = ec1ToHC ro (lctUniversalTime dt) sunEC1
-- HC {hAltitude = DD 49.30604722222222, hAzimuth = DD 118.92209166666666}
@
-}

module Data.Astro.Coordinate
(
  DecimalDegrees(..)
  , DecimalHours(..)
  , HorizonCoordinates(..)
  , EquatorialCoordinates1(..)
  , EquatorialCoordinates2(..)
  , EclipticCoordinates(..)
  , GalacticCoordinates(..)
  , raToHA
  , haToRA
  , equatorialToHorizon
  , horizonToEquatorial
  , ec1ToHC
  , hcToEC1
  , ecHCConv
  , obliquity
  , eclipticToEquatorial
  , equatorialToEcliptic
  , galacticToEquatorial
  , equatorialToGalactic
)

where

import Data.Astro.Time (utToLST)
import Data.Astro.Time.JulianDate (JulianDate(..), numberOfCenturies, splitToDayAndTime)
import Data.Astro.Time.Epoch (j2000)
import Data.Astro.Time.Sidereal (LocalSiderealTime(..), lstToDH)
import Data.Astro.Types (DecimalDegrees(..), DecimalHours(..)
                        , fromDecimalHours, toDecimalHours
                        , toRadians, fromRadians, fromDMS
                        , GeographicCoordinates(..))
import Data.Astro.Utils (fromFixed)
import Data.Astro.Effects.Nutation (nutationObliquity)


-- | Horizon Coordinates, for details see the module's description

data HorizonCoordinates = HC {
  HorizonCoordinates -> DecimalDegrees
hAltitude :: DecimalDegrees   -- ^ alpha

  , HorizonCoordinates -> DecimalDegrees
hAzimuth :: DecimalDegrees  -- ^ big alpha

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


-- | Equatorial Coordinates, defines fixed position in the sky

data EquatorialCoordinates1 = EC1 {
  EquatorialCoordinates1 -> DecimalDegrees
e1Declination :: DecimalDegrees     -- ^ delta

  , EquatorialCoordinates1 -> DecimalHours
e1RightAscension :: DecimalHours  -- ^ alpha

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


-- | Equatorial Coordinates

data EquatorialCoordinates2 = EC2 {
  EquatorialCoordinates2 -> DecimalDegrees
e2Declination :: DecimalDegrees    -- ^ delta

  , EquatorialCoordinates2 -> DecimalHours
e2HoursAngle :: DecimalHours     -- ^ H

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


-- | Ecliptic Coordinates

data EclipticCoordinates = EcC {
  EclipticCoordinates -> DecimalDegrees
ecLatitude :: DecimalDegrees      -- ^ beta

  , EclipticCoordinates -> DecimalDegrees
ecLongitude :: DecimalDegrees   -- ^ lambda

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


-- | Galactic Coordinates

data GalacticCoordinates = GC {
  GalacticCoordinates -> DecimalDegrees
gLatitude :: DecimalDegrees       -- ^ b

  , GalacticCoordinates -> DecimalDegrees
gLongitude :: DecimalDegrees    -- ^ l

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


-- | Convert Right Ascension to Hour Angle for specified longitude and Universal Time

raToHA :: DecimalHours -> DecimalDegrees -> JulianDate -> DecimalHours
raToHA :: DecimalHours -> DecimalDegrees -> JulianDate -> DecimalHours
raToHA = DecimalHours -> DecimalDegrees -> JulianDate -> DecimalHours
haRAConv


-- | Convert Hour Angle to Right Ascension for specified longitude and Universal Time

haToRA :: DecimalHours -> DecimalDegrees -> JulianDate -> DecimalHours
haToRA :: DecimalHours -> DecimalDegrees -> JulianDate -> DecimalHours
haToRA = DecimalHours -> DecimalDegrees -> JulianDate -> DecimalHours
haRAConv


-- | HA <-> RA Conversions

haRAConv :: DecimalHours -> DecimalDegrees -> JulianDate -> DecimalHours
haRAConv :: DecimalHours -> DecimalDegrees -> JulianDate -> DecimalHours
haRAConv DecimalHours
dh DecimalDegrees
longitude JulianDate
ut =
  let lst :: LocalSiderealTime
lst = DecimalDegrees -> JulianDate -> LocalSiderealTime
utToLST DecimalDegrees
longitude JulianDate
ut  -- Local Sidereal Time

      DH Double
hourAngle = (LocalSiderealTime -> DecimalHours
lstToDH LocalSiderealTime
lst) DecimalHours -> DecimalHours -> DecimalHours
forall a. Num a => a -> a -> a
- DecimalHours
dh
  in if Double
hourAngle Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 then (Double -> DecimalHours
DH (Double -> DecimalHours) -> Double -> DecimalHours
forall a b. (a -> b) -> a -> b
$ Double
hourAngleDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
24) else (Double -> DecimalHours
DH Double
hourAngle)


-- | Convert Equatorial Coordinates to Horizon Coordinates.

-- It takes a latitude of the observer and 'EquatorialCoordinates2'.

-- If you need to convert 'EquatorialCoordinates1'

-- you may use 'raToHa' function to obtain 'EquatorialCoordinates2'

-- or just use function-shortcut 'ec1ToHC' straightaway.

-- The functions returns 'HorizonCoordinates'.

equatorialToHorizon :: DecimalDegrees -> EquatorialCoordinates2 -> HorizonCoordinates
equatorialToHorizon :: DecimalDegrees -> EquatorialCoordinates2 -> HorizonCoordinates
equatorialToHorizon DecimalDegrees
latitude (EC2 DecimalDegrees
dec DecimalHours
hourAngle) =
  let hourAngle' :: DecimalDegrees
hourAngle' = DecimalHours -> DecimalDegrees
fromDecimalHours DecimalHours
hourAngle
      (DecimalDegrees
altitude, DecimalDegrees
azimuth) = DecimalDegrees
-> (DecimalDegrees, DecimalDegrees)
-> (DecimalDegrees, DecimalDegrees)
ecHCConv DecimalDegrees
latitude (DecimalDegrees
dec, DecimalDegrees
hourAngle')
  in DecimalDegrees -> DecimalDegrees -> HorizonCoordinates
HC DecimalDegrees
altitude DecimalDegrees
azimuth


-- | Convert Horizon Coordinates to Equatorial Coordinates.

-- It takes a latitude of the observer and 'HorizonCoordinates'.

-- The functions returns 'EquatorialCoordinates2'.

-- If you need to obtain 'EquatorialCoordinates1' you may use 'haToRa' function,

-- or function-shortcut `hcToEC1`.

horizonToEquatorial :: DecimalDegrees -> HorizonCoordinates -> EquatorialCoordinates2
horizonToEquatorial :: DecimalDegrees -> HorizonCoordinates -> EquatorialCoordinates2
horizonToEquatorial DecimalDegrees
latitude (HC DecimalDegrees
altitude DecimalDegrees
azimuth) =
  let (DecimalDegrees
dec, DecimalDegrees
hourAngle) = DecimalDegrees
-> (DecimalDegrees, DecimalDegrees)
-> (DecimalDegrees, DecimalDegrees)
ecHCConv DecimalDegrees
latitude (DecimalDegrees
altitude, DecimalDegrees
azimuth)
  in DecimalDegrees -> DecimalHours -> EquatorialCoordinates2
EC2 DecimalDegrees
dec (DecimalHours -> EquatorialCoordinates2)
-> DecimalHours -> EquatorialCoordinates2
forall a b. (a -> b) -> a -> b
$ DecimalDegrees -> DecimalHours
toDecimalHours DecimalDegrees
hourAngle


-- | Convert Equatorial Coordinates (Type 1) to Horizon Coordinates.

-- This is function shortcut - tt combines `equatorialToHorizon` and `raToHA`.

-- It takes geographic coordinates of the observer, universal time and equatorial coordinates.

ec1ToHC :: GeographicCoordinates -> JulianDate -> EquatorialCoordinates1 -> HorizonCoordinates
ec1ToHC :: GeographicCoordinates
-> JulianDate -> EquatorialCoordinates1 -> HorizonCoordinates
ec1ToHC (GeoC DecimalDegrees
latitude DecimalDegrees
longitude) JulianDate
jd (EC1 DecimalDegrees
delta DecimalHours
alpha) =
  let ec2 :: EquatorialCoordinates2
ec2 = DecimalDegrees -> DecimalHours -> EquatorialCoordinates2
EC2 DecimalDegrees
delta (DecimalHours -> DecimalDegrees -> JulianDate -> DecimalHours
raToHA DecimalHours
alpha DecimalDegrees
longitude JulianDate
jd)
  in DecimalDegrees -> EquatorialCoordinates2 -> HorizonCoordinates
equatorialToHorizon DecimalDegrees
latitude EquatorialCoordinates2
ec2


-- | Convert Horizon Coordinates to Equatorial Coordinates (Type 1).

-- This is function shortcut - tt combines `horizonToEquatorial` and `haToRA`.

-- It takes geographic coordinates of the observer, universal time and horizon coordinates.

hcToEC1 :: GeographicCoordinates -> JulianDate -> HorizonCoordinates -> EquatorialCoordinates1
hcToEC1 :: GeographicCoordinates
-> JulianDate -> HorizonCoordinates -> EquatorialCoordinates1
hcToEC1 (GeoC DecimalDegrees
latitude DecimalDegrees
longitude) JulianDate
jd HorizonCoordinates
hc =
  let (EC2 DecimalDegrees
dec DecimalHours
hourAngle) = DecimalDegrees -> HorizonCoordinates -> EquatorialCoordinates2
horizonToEquatorial DecimalDegrees
latitude HorizonCoordinates
hc
  in DecimalDegrees -> DecimalHours -> EquatorialCoordinates1
EC1 DecimalDegrees
dec (DecimalHours -> DecimalDegrees -> JulianDate -> DecimalHours
haToRA DecimalHours
hourAngle DecimalDegrees
longitude JulianDate
jd)


-- | Function converts Equatorial Coordinates To Horizon Coordinates and vice versa

-- It takes a latitide of the observer as a first parameter and a pair of 'how far up' and 'how far round' coordinates

-- as a second parameter. It returns a pair of 'how far up' and 'how far round' coordinates.

ecHCConv :: DecimalDegrees -> (DecimalDegrees, DecimalDegrees) -> (DecimalDegrees, DecimalDegrees)
ecHCConv :: DecimalDegrees
-> (DecimalDegrees, DecimalDegrees)
-> (DecimalDegrees, DecimalDegrees)
ecHCConv DecimalDegrees
latitude (DecimalDegrees
up, DecimalDegrees
round) =
  let latitude' :: Double
latitude' = DecimalDegrees -> Double
toRadians DecimalDegrees
latitude
      up' :: Double
up' = DecimalDegrees -> Double
toRadians DecimalDegrees
up
      round' :: Double
round' = DecimalDegrees -> Double
toRadians DecimalDegrees
round
      sinUpResult :: Double
sinUpResult = (Double -> Double
forall a. Floating a => a -> a
sin Double
up')Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
sin Double
latitude') Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double -> Double
forall a. Floating a => a -> a
cos Double
up')Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
cos Double
latitude')Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
cos Double
round')
      upResult :: Double
upResult = Double -> Double
forall a. Floating a => a -> a
asin Double
sinUpResult
      roundResult :: Double
roundResult = Double -> Double
forall a. Floating a => a -> a
acos (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ ((Double -> Double
forall a. Floating a => a -> a
sin Double
up') Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double -> Double
forall a. Floating a => a -> a
sin Double
latitude')Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
sinUpResult) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ ((Double -> Double
forall a. Floating a => a -> a
cos Double
latitude') Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double -> Double
forall a. Floating a => a -> a
cos Double
upResult))
      roundResult' :: Double
roundResult' = if (Double -> Double
forall a. Floating a => a -> a
sin Double
round') Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 then Double
roundResult else (Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
roundResult)
  in ((Double -> DecimalDegrees
fromRadians Double
upResult), (Double -> DecimalDegrees
fromRadians Double
roundResult'))


-- | Calculate the obliquity of the ecpliptic on JulianDate

obliquity :: JulianDate -> DecimalDegrees
obliquity :: JulianDate -> DecimalDegrees
obliquity JulianDate
jd =
  let DD Double
baseObliquity = Int -> Int -> Double -> DecimalDegrees
forall a. RealFrac a => Int -> Int -> a -> DecimalDegrees
fromDMS Int
23 Int
26 Double
21.45
      t :: Double
t = JulianDate -> JulianDate -> Double
numberOfCenturies JulianDate
j2000 JulianDate
jd
      de :: Double
de = (Double
46.815Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
0.0006Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
tDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
0.00181Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
tDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
tDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
t) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3600  -- 3600 number of seconds in 1 degree

  in (Double -> DecimalDegrees
DD (Double -> DecimalDegrees) -> Double -> DecimalDegrees
forall a b. (a -> b) -> a -> b
$ Double
baseObliquity Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
de) DecimalDegrees -> DecimalDegrees -> DecimalDegrees
forall a. Num a => a -> a -> a
+ (JulianDate -> DecimalDegrees
nutationObliquity JulianDate
jd)


-- | Converts Ecliptic Coordinates on specified Julian Date to Equatorial Coordinates

eclipticToEquatorial :: EclipticCoordinates -> JulianDate -> EquatorialCoordinates1
eclipticToEquatorial :: EclipticCoordinates -> JulianDate -> EquatorialCoordinates1
eclipticToEquatorial (EcC DecimalDegrees
beta DecimalDegrees
gamma) JulianDate
jd =
  let epsilon' :: Double
epsilon' = DecimalDegrees -> Double
toRadians (DecimalDegrees -> Double) -> DecimalDegrees -> Double
forall a b. (a -> b) -> a -> b
$ JulianDate -> DecimalDegrees
obliquity JulianDate
jd
      beta' :: Double
beta' = DecimalDegrees -> Double
toRadians DecimalDegrees
beta
      gamma' :: Double
gamma' = DecimalDegrees -> Double
toRadians DecimalDegrees
gamma
      delta :: Double
delta = Double -> Double
forall a. Floating a => a -> a
asin (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ (Double -> Double
forall a. Floating a => a -> a
sin Double
beta')Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
cos Double
epsilon') Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double -> Double
forall a. Floating a => a -> a
cos Double
beta')Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
sin Double
epsilon')Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
sin Double
gamma')
      y :: Double
y = (Double -> Double
forall a. Floating a => a -> a
sin Double
gamma')Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
cos Double
epsilon') Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double -> Double
forall a. Floating a => a -> a
tan Double
beta')Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
sin Double
epsilon')
      x :: Double
x = Double -> Double
forall a. Floating a => a -> a
cos Double
gamma'
      alpha :: Double
alpha = Double -> Double
forall a. (Floating a, Ord a) => a -> a
reduceToZero2PI (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 Double
y Double
x
  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)


-- | Converts Equatorial Coordinates to Ecliptic Coordinates on specified Julian Date

equatorialToEcliptic :: EquatorialCoordinates1 -> JulianDate -> EclipticCoordinates
equatorialToEcliptic :: EquatorialCoordinates1 -> JulianDate -> EclipticCoordinates
equatorialToEcliptic (EC1 DecimalDegrees
delta DecimalHours
alpha) JulianDate
jd =
  let epsilon' :: Double
epsilon' = DecimalDegrees -> Double
toRadians (DecimalDegrees -> Double) -> DecimalDegrees -> Double
forall a b. (a -> b) -> a -> b
$ JulianDate -> DecimalDegrees
obliquity JulianDate
jd
      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
      beta :: Double
beta = Double -> Double
forall a. Floating a => a -> a
asin (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ (Double -> Double
forall a. Floating a => a -> a
sin Double
delta')Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
cos Double
epsilon') Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double -> Double
forall a. Floating a => a -> a
cos Double
delta')Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
sin Double
epsilon')Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
sin Double
alpha')
      y :: Double
y = (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
cos Double
epsilon') 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 -> Double
forall a. Floating a => a -> a
sin Double
epsilon')
      x :: Double
x = Double -> Double
forall a. Floating a => a -> a
cos Double
alpha'
      gamma :: Double
gamma = Double -> Double
forall a. (Floating a, Ord a) => a -> a
reduceToZero2PI (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 Double
y Double
x
  in DecimalDegrees -> DecimalDegrees -> EclipticCoordinates
EcC (Double -> DecimalDegrees
fromRadians Double
beta) (Double -> DecimalDegrees
fromRadians Double
gamma)


-- | Galactic Pole Coordinates

galacticPole :: EquatorialCoordinates1
galacticPole :: EquatorialCoordinates1
galacticPole = DecimalDegrees -> DecimalHours -> EquatorialCoordinates1
EC1 (Double -> DecimalDegrees
DD Double
27.4) (DecimalDegrees -> DecimalHours
toDecimalHours (DecimalDegrees -> DecimalHours) -> DecimalDegrees -> DecimalHours
forall a b. (a -> b) -> a -> b
$ Double -> DecimalDegrees
DD Double
192.25)

galacticPoleInRadians :: (Double, Double)
galacticPoleInRadians = (Double
delta, Double
alpha)
  where delta :: Double
delta = DecimalDegrees -> Double
toRadians (DecimalDegrees -> Double) -> DecimalDegrees -> Double
forall a b. (a -> b) -> a -> b
$ EquatorialCoordinates1 -> DecimalDegrees
e1Declination EquatorialCoordinates1
galacticPole
        alpha :: Double
alpha = DecimalDegrees -> Double
toRadians (DecimalDegrees -> Double) -> DecimalDegrees -> Double
forall a b. (a -> b) -> a -> b
$ DecimalHours -> DecimalDegrees
fromDecimalHours (DecimalHours -> DecimalDegrees) -> DecimalHours -> DecimalDegrees
forall a b. (a -> b) -> a -> b
$ EquatorialCoordinates1 -> DecimalHours
e1RightAscension EquatorialCoordinates1
galacticPole


-- | Ascending node of the galactic place on equator

ascendingNode :: DecimalDegrees
ascendingNode :: DecimalDegrees
ascendingNode = Double -> DecimalDegrees
DD Double
33


-- | Convert Galactic Coordinates Equatorial Coordinates

galacticToEquatorial :: GalacticCoordinates -> EquatorialCoordinates1
galacticToEquatorial :: GalacticCoordinates -> EquatorialCoordinates1
galacticToEquatorial (GC DecimalDegrees
b DecimalDegrees
l) =
  let b' :: Double
b' = DecimalDegrees -> Double
toRadians DecimalDegrees
b
      l' :: Double
l' = DecimalDegrees -> Double
toRadians DecimalDegrees
l
      (Double
poleDelta, Double
poleAlpha) = (Double, Double)
galacticPoleInRadians
      an :: Double
an = DecimalDegrees -> Double
toRadians DecimalDegrees
ascendingNode
      delta :: Double
delta = Double -> Double
forall a. Floating a => a -> a
asin (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ (Double -> Double
forall a. Floating a => a -> a
cos Double
b')Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
cos Double
poleDelta)Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
sin (Double
l'Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
an)) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double -> Double
forall a. Floating a => a -> a
sin Double
b')Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
sin Double
poleDelta)
      y :: Double
y = (Double -> Double
forall a. Floating a => a -> a
cos Double
b')Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
cos (Double
l'Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
an))
      x :: Double
x = (Double -> Double
forall a. Floating a => a -> a
sin Double
b')Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
cos Double
poleDelta) Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double -> Double
forall a. Floating a => a -> a
cos Double
b')Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
sin Double
poleDelta)Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
sin (Double
l'Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
an))
      alpha :: Double
alpha = Double -> Double
forall a. (Floating a, Ord a) => a -> a
reduceToZero2PI (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ (Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 Double
y Double
x) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
poleAlpha
  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)


-- | Convert Equatorial Coordinates to Galactic Coordinates

equatorialToGalactic :: EquatorialCoordinates1 -> GalacticCoordinates
equatorialToGalactic :: EquatorialCoordinates1 -> GalacticCoordinates
equatorialToGalactic (EC1 DecimalDegrees
delta DecimalHours
alpha) =
  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
      (Double
poleDelta, Double
poleAlpha) = (Double, Double)
galacticPoleInRadians
      sinb :: Double
sinb = (Double -> Double
forall a. Floating a => a -> a
cos Double
delta')Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
cos Double
poleDelta)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
poleAlpha)) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double -> Double
forall a. Floating a => a -> a
sin Double
delta') Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double -> Double
forall a. Floating a => a -> a
sin Double
poleDelta)
      y :: Double
y = (Double -> Double
forall a. Floating a => a -> a
sin Double
delta') Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
sinbDouble -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
sin Double
poleDelta)
      x :: Double
x = (Double -> Double
forall a. Floating a => a -> a
cos Double
delta')Double -> 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
poleAlpha))Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
cos Double
poleDelta)
      b :: Double
b = Double -> Double
forall a. Floating a => a -> a
asin Double
sinb
      l :: Double
l = Double -> Double
forall a. (Floating a, Ord a) => a -> a
reduceToZero2PI (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ (Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 Double
y Double
x) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (DecimalDegrees -> Double
toRadians DecimalDegrees
ascendingNode)
  in DecimalDegrees -> DecimalDegrees -> GalacticCoordinates
GC (Double -> DecimalDegrees
fromRadians Double
b) (Double -> DecimalDegrees
fromRadians Double
l)


-- | Reduce angle from [-pi, pi] to [0, 2*pi]

-- Usefull to correct results of atan2 for 'how far round' coordinates

reduceToZero2PI :: (Floating a, Ord a) => a -> a
reduceToZero2PI :: a -> a
reduceToZero2PI a
rad = if a
rad a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0 then a
rad a -> a -> a
forall a. Num a => a -> a -> a
+ a
2a -> a -> a
forall a. Num a => a -> a -> a
*a
forall a. Floating a => a
pi else a
rad