{-|
Module: Data.Astro.CelestialObject.RiseSet
Description: Computations rise and set of selestial objects
Copyright: Alexander Ignatyev, 2016

Computations rise and set of selestial objects.

= Examples

== /Stars/

See "Data.Astro.Star" module for example.

== /Planets/

See "Data.Astro.Planet" module for example.
-}

module Data.Astro.CelestialObject.RiseSet
(
  RiseSet(..)
  , RSInfo(..)
  , RiseSetLST(..)
  , RiseSetLCT(..)
  , RiseSetMB(..)
  , riseAndSet
  , riseAndSet2
  , riseAndSetLCT
  , toRiseSetLCT
)

where

import Data.Astro.Types (DecimalDegrees, DecimalHours(..)
                        , GeographicCoordinates(..)
                        , toRadians, fromRadians
                        , toDecimalHours)
import Data.Astro.Utils (reduceToZeroRange)
import Data.Astro.Time (lstToLCT)
import Data.Astro.Time.JulianDate (JulianDate(..), LocalCivilTime(..), LocalCivilDate(..), addHours)
import Data.Astro.Time.Sidereal (LocalSiderealTime, dhToLST)
import Data.Astro.Coordinate (EquatorialCoordinates1(..))

-- | Some Info of Rise and Set of a celestial object

data RiseSet a
  -- | Some Info of Rise and Set of the celestial object

  = RiseSet a a
  -- | The celestial object is always above the horizon

  | Circumpolar
  -- | The celestial object is always below the horizon

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


-- | Rise or Set time and azimuth

type RSInfo a = (a, DecimalDegrees)


-- | LST (Local Sidereal Time) and Azimuth of Rise and Set

type RiseSetLST = RiseSet (RSInfo LocalSiderealTime)


-- | Local Civil Time and Azimuth of Rise and Set

type RiseSetLCT = RiseSet (RSInfo LocalCivilTime)


-- | The optional Rise And optinal Set Information (LocalCivilTime and Azimuth)

type RiseSetMB = RiseSet (Maybe (RSInfo LocalCivilTime))


-- | Calculate rise and set local sidereal time of a celestial object.

-- It takes the equatorial coordinates of the celestial object,

-- vertical shift and the latitude of the observation.

-- To calculate /vertical shift/ for stars use function 'refract' from "Data.Astro.Effects".

-- In most cases you can assume that /vertical shift/ equals 0.566569 (34 arcmins ~ 'refract (DD 0) 12 1012').

riseAndSet :: EquatorialCoordinates1 -> DecimalDegrees -> DecimalDegrees -> RiseSetLST
riseAndSet :: EquatorialCoordinates1
-> DecimalDegrees -> DecimalDegrees -> RiseSetLST
riseAndSet (EC1 DecimalDegrees
delta DecimalHours
alpha) DecimalDegrees
shift DecimalDegrees
lat =
  let delta' :: Double
delta' = DecimalDegrees -> Double
toRadians DecimalDegrees
delta
      shift' :: Double
shift' = DecimalDegrees -> Double
toRadians DecimalDegrees
shift
      lat' :: Double
lat' = DecimalDegrees -> Double
toRadians DecimalDegrees
lat
      cosH :: Double
cosH = Double -> Double -> Double -> Double
cosOfHourAngle Double
delta' Double
shift' Double
lat'
  in Double -> Double -> Double -> Double -> RiseSetLST
sortRiseSet Double
cosH Double
delta' Double
shift' Double
lat'

  where sortRiseSet :: Double -> Double -> Double -> Double -> RiseSetLST
        sortRiseSet :: Double -> Double -> Double -> Double -> RiseSetLST
sortRiseSet Double
cosH Double
delta Double
shift Double
latitude
          | Double
cosH Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< -Double
1 = RiseSetLST
forall a. RiseSet a
Circumpolar
          | Double
cosH Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1 = RiseSetLST
forall a. RiseSet a
NeverRises
          | Bool
otherwise = DecimalHours
-> DecimalHours -> Double -> Double -> Double -> RiseSetLST
calcTimesAndAzimuths DecimalHours
alpha (Double -> DecimalHours
toHours (Double -> DecimalHours) -> Double -> DecimalHours
forall a b. (a -> b) -> a -> b
$ Double -> Double
forall a. Floating a => a -> a
acos Double
cosH) Double
delta Double
shift Double
latitude

        toHours :: Double -> DecimalHours
        toHours :: Double -> DecimalHours
toHours = DecimalDegrees -> DecimalHours
toDecimalHours (DecimalDegrees -> DecimalHours)
-> (Double -> DecimalDegrees) -> Double -> DecimalHours
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> DecimalDegrees
fromRadians

        cosOfHourAngle :: Double -> Double -> Double -> Double
        cosOfHourAngle :: Double -> Double -> Double -> Double
cosOfHourAngle Double
delta Double
shift Double
latitude = -((Double -> Double
forall a. Floating a => a -> a
sin Double
shift) 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
sin Double
delta)) 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
delta))

        calcTimesAndAzimuths :: DecimalHours -> DecimalHours -> Double -> Double -> Double -> RiseSetLST
        calcTimesAndAzimuths :: DecimalHours
-> DecimalHours -> Double -> Double -> Double -> RiseSetLST
calcTimesAndAzimuths DecimalHours
alpha DecimalHours
hourAngle Double
delta Double
shift Double
latitude =
          let lstRise :: LocalSiderealTime
lstRise = DecimalHours -> LocalSiderealTime
dhToLST (DecimalHours -> LocalSiderealTime)
-> DecimalHours -> LocalSiderealTime
forall a b. (a -> b) -> a -> b
$ DecimalHours -> DecimalHours -> DecimalHours
forall a. RealFrac a => a -> a -> a
reduceToZeroRange DecimalHours
24 (DecimalHours -> DecimalHours) -> DecimalHours -> DecimalHours
forall a b. (a -> b) -> a -> b
$ DecimalHours
alpha DecimalHours -> DecimalHours -> DecimalHours
forall a. Num a => a -> a -> a
- DecimalHours
hourAngle
              lstSet :: LocalSiderealTime
lstSet = DecimalHours -> LocalSiderealTime
dhToLST (DecimalHours -> LocalSiderealTime)
-> DecimalHours -> LocalSiderealTime
forall a b. (a -> b) -> a -> b
$ DecimalHours -> DecimalHours -> DecimalHours
forall a. RealFrac a => a -> a -> a
reduceToZeroRange DecimalHours
24 (DecimalHours -> DecimalHours) -> DecimalHours -> DecimalHours
forall a b. (a -> b) -> a -> b
$ DecimalHours
alpha DecimalHours -> DecimalHours -> DecimalHours
forall a. Num a => a -> a -> a
+ DecimalHours
hourAngle
              azimuthRise :: Double
azimuthRise = Double -> Double -> Double
forall a. RealFrac a => a -> a -> a
reduceToZeroRange (Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall a. Floating a => a
pi) (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ 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
delta) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double -> Double
forall a. Floating a => a -> a
sin Double
shift)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. Fractional a => a -> a -> a
/ ((Double -> Double
forall a. Floating a => a -> a
cos Double
shift)Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
cos Double
latitude))
              azimuthSet :: Double
azimuthSet = 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
azimuthRise
          in (LocalSiderealTime, DecimalDegrees)
-> (LocalSiderealTime, DecimalDegrees) -> RiseSetLST
forall a. a -> a -> RiseSet a
RiseSet (LocalSiderealTime
lstRise, Double -> DecimalDegrees
fromRadians Double
azimuthRise) (LocalSiderealTime
lstSet, Double -> DecimalDegrees
fromRadians Double
azimuthSet)


-- | Calculate rise and set local sidereal time of a celestial object

-- that changes its equatorial coordinates during the day (the Sun, the Moon, planets).

-- It takes epsilon, the function that returns equatorial coordinates of the celestial object for a given julian date,

-- vertical shift and the latitude of the observation.

-- To calculate /vertical shift/ for stars use function 'refract' from "Data.Astro.Effects".

-- In most cases you can assume that /vertical shift/ equals 0.566569 (34 arcmins ~ 'refract (DD 0) 12 1012').

riseAndSet2 :: DecimalHours
               -> (JulianDate -> EquatorialCoordinates1)
               -> GeographicCoordinates
               -> DecimalDegrees
               -> LocalCivilDate
               -> RiseSetMB
riseAndSet2 :: DecimalHours
-> (JulianDate -> EquatorialCoordinates1)
-> GeographicCoordinates
-> DecimalDegrees
-> LocalCivilDate
-> RiseSetMB
riseAndSet2 DecimalHours
eps JulianDate -> EquatorialCoordinates1
getPosition GeographicCoordinates
geoc DecimalDegrees
shift LocalCivilDate
lcd =
  let day :: JulianDate
day = LocalCivilDate -> JulianDate
lcdDate LocalCivilDate
lcd
      pos :: EquatorialCoordinates1
pos = JulianDate -> EquatorialCoordinates1
getPosition (DecimalHours -> JulianDate -> JulianDate
addHours DecimalHours
12 JulianDate
day)
      rs :: RiseSetLCT
rs = GeographicCoordinates
-> LocalCivilDate
-> DecimalDegrees
-> EquatorialCoordinates1
-> RiseSetLCT
riseAndSetLCT GeographicCoordinates
geoc LocalCivilDate
lcd DecimalDegrees
shift EquatorialCoordinates1
pos
      rise :: RiseSetLCT
rise = (RiseSetLCT -> RSInfo LocalCivilTime)
-> RSInfo LocalCivilTime -> Int -> RiseSetLCT
calc RiseSetLCT -> RSInfo LocalCivilTime
getRiseTime (RiseSetLCT -> RSInfo LocalCivilTime
getRiseTime RiseSetLCT
rs) Int
0
      set :: RiseSetLCT
set = (RiseSetLCT -> RSInfo LocalCivilTime)
-> RSInfo LocalCivilTime -> Int -> RiseSetLCT
calc RiseSetLCT -> RSInfo LocalCivilTime
getSetTime (RiseSetLCT -> RSInfo LocalCivilTime
getSetTime RiseSetLCT
rs) Int
0
  in case RiseSetLCT
rs of
    RiseSetLCT
Circumpolar -> RiseSetMB
forall a. RiseSet a
Circumpolar
    RiseSetLCT
NeverRises -> RiseSetMB
forall a. RiseSet a
NeverRises
    RiseSetLCT
_ -> RiseSetLCT -> RiseSetLCT -> RiseSetMB
forall a. RiseSet a -> RiseSet a -> RiseSet (Maybe a)
buildResult RiseSetLCT
rise RiseSetLCT
set

  where calc :: (RiseSetLCT -> RSInfo LocalCivilTime) -> RSInfo LocalCivilTime -> Int -> RiseSetLCT
        calc :: (RiseSetLCT -> RSInfo LocalCivilTime)
-> RSInfo LocalCivilTime -> Int -> RiseSetLCT
calc RiseSetLCT -> RSInfo LocalCivilTime
getRSInfo rsi :: RSInfo LocalCivilTime
rsi@(LocalCivilTime
time, DecimalDegrees
_) Int
iterNo =
          let pos :: EquatorialCoordinates1
pos = JulianDate -> EquatorialCoordinates1
getPosition (JulianDate -> EquatorialCoordinates1)
-> JulianDate -> EquatorialCoordinates1
forall a b. (a -> b) -> a -> b
$ LocalCivilTime -> JulianDate
lctUniversalTime LocalCivilTime
time
              rs :: RiseSetLCT
rs = GeographicCoordinates
-> LocalCivilDate
-> DecimalDegrees
-> EquatorialCoordinates1
-> RiseSetLCT
riseAndSetLCT GeographicCoordinates
geoc LocalCivilDate
lcd DecimalDegrees
shift EquatorialCoordinates1
pos
              rsi' :: RSInfo LocalCivilTime
rsi' = RiseSetLCT -> RSInfo LocalCivilTime
getRSInfo RiseSetLCT
rs
          in case RiseSetLCT
rs of
            RiseSetLCT
Circumpolar -> RiseSetLCT
forall a. RiseSet a
Circumpolar
            RiseSetLCT
NeverRises -> RiseSetLCT
forall a. RiseSet a
NeverRises
            RiseSetLCT
_ -> if RSInfo LocalCivilTime -> RSInfo LocalCivilTime -> Bool
isOK RSInfo LocalCivilTime
rsi RSInfo LocalCivilTime
rsi' Bool -> Bool -> Bool
|| Int
iterNo Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
maxIters
                 then RiseSetLCT
rs
                 else (RiseSetLCT -> RSInfo LocalCivilTime)
-> RSInfo LocalCivilTime -> Int -> RiseSetLCT
calc RiseSetLCT -> RSInfo LocalCivilTime
getRSInfo RSInfo LocalCivilTime
rsi' (Int
iterNoInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)

        isOK :: RSInfo LocalCivilTime -> RSInfo LocalCivilTime -> Bool
        isOK :: RSInfo LocalCivilTime -> RSInfo LocalCivilTime -> Bool
isOK (LocalCivilTime
t1, DecimalDegrees
_) (LocalCivilTime
t2, DecimalDegrees
_) = (Double -> Double
forall a. Num a => a -> a
abs Double
d) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< (Double
hDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
24)
          where JD Double
d = (LocalCivilTime -> JulianDate
lctUniversalTime LocalCivilTime
t1) JulianDate -> JulianDate -> JulianDate
forall a. Num a => a -> a -> a
- (LocalCivilTime -> JulianDate
lctUniversalTime LocalCivilTime
t2)
                DH Double
h = DecimalHours
eps

        maxIters :: Int
maxIters = Int
3

        getRiseTime :: RiseSetLCT -> RSInfo LocalCivilTime
        getRiseTime :: RiseSetLCT -> RSInfo LocalCivilTime
getRiseTime (RiseSet RSInfo LocalCivilTime
r RSInfo LocalCivilTime
_) = RSInfo LocalCivilTime
r

        getSetTime :: RiseSetLCT -> RSInfo LocalCivilTime
        getSetTime :: RiseSetLCT -> RSInfo LocalCivilTime
getSetTime (RiseSet RSInfo LocalCivilTime
_ RSInfo LocalCivilTime
s) = RSInfo LocalCivilTime
s

        buildResult :: RiseSet a -> RiseSet a -> RiseSet (Maybe a)
buildResult (RiseSet a
r a
_) (RiseSet a
_ a
s) = Maybe a -> Maybe a -> RiseSet (Maybe a)
forall a. a -> a -> RiseSet a
RiseSet (a -> Maybe a
forall a. a -> Maybe a
Just a
r) (a -> Maybe a
forall a. a -> Maybe a
Just a
s)
        buildResult (RiseSet a
r a
_) RiseSet a
_ = Maybe a -> Maybe a -> RiseSet (Maybe a)
forall a. a -> a -> RiseSet a
RiseSet (a -> Maybe a
forall a. a -> Maybe a
Just a
r) Maybe a
forall a. Maybe a
Nothing
        buildResult RiseSet a
_ (RiseSet a
_ a
s) = Maybe a -> Maybe a -> RiseSet (Maybe a)
forall a. a -> a -> RiseSet a
RiseSet Maybe a
forall a. Maybe a
Nothing (a -> Maybe a
forall a. a -> Maybe a
Just a
s)



-- | Calculates set and rise of the celestial object

-- It takes geographic coordinates of the observer, local civil date, vertical shift

-- and equatorial coordinates of the celestial object.

riseAndSetLCT :: GeographicCoordinates
                -> LocalCivilDate
                -> DecimalDegrees
                -> EquatorialCoordinates1
                -> RiseSetLCT
riseAndSetLCT :: GeographicCoordinates
-> LocalCivilDate
-> DecimalDegrees
-> EquatorialCoordinates1
-> RiseSetLCT
riseAndSetLCT (GeoC DecimalDegrees
latitude DecimalDegrees
longitude) LocalCivilDate
lcd DecimalDegrees
shift EquatorialCoordinates1
ec
  = DecimalDegrees -> LocalCivilDate -> RiseSetLST -> RiseSetLCT
toRiseSetLCT DecimalDegrees
longitude LocalCivilDate
lcd (RiseSetLST -> RiseSetLCT) -> RiseSetLST -> RiseSetLCT
forall a b. (a -> b) -> a -> b
$ EquatorialCoordinates1
-> DecimalDegrees -> DecimalDegrees -> RiseSetLST
riseAndSet EquatorialCoordinates1
ec DecimalDegrees
shift DecimalDegrees
latitude


-- | Converts Rise and Set in Local Sidereal Time to Rise and Set in Local Civil Time.

-- It takes longutude of the observer and local civil date.

-- To calculate /vertical shift/ for stars use function 'refract' from "Data.Astro.Effects".

-- In most cases you can assume that /vertical shift/ equals 0.566569 (34 arcmins ~ 'refract (DD 0) 12 1012').

toRiseSetLCT :: DecimalDegrees
               -> LocalCivilDate
               -> RiseSetLST
               -> RiseSetLCT
toRiseSetLCT :: DecimalDegrees -> LocalCivilDate -> RiseSetLST -> RiseSetLCT
toRiseSetLCT DecimalDegrees
longitude LocalCivilDate
lcd (RiseSet (LocalSiderealTime
rise, DecimalDegrees
azRise) (LocalSiderealTime
set, DecimalDegrees
azSet)) =
  let toLCT :: LocalSiderealTime -> LocalCivilTime
toLCT LocalSiderealTime
lst = DecimalDegrees
-> LocalCivilDate -> LocalSiderealTime -> LocalCivilTime
lstToLCT DecimalDegrees
longitude LocalCivilDate
lcd LocalSiderealTime
lst
      rise' :: LocalCivilTime
rise' = LocalSiderealTime -> LocalCivilTime
toLCT LocalSiderealTime
rise
      set' :: LocalCivilTime
set' = LocalSiderealTime -> LocalCivilTime
toLCT LocalSiderealTime
set
  in RSInfo LocalCivilTime -> RSInfo LocalCivilTime -> RiseSetLCT
forall a. a -> a -> RiseSet a
RiseSet (LocalCivilTime
rise', DecimalDegrees
azRise) (LocalCivilTime
set', DecimalDegrees
azSet)
toRiseSetLCT DecimalDegrees
_ LocalCivilDate
_ RiseSetLST
Circumpolar  = RiseSetLCT
forall a. RiseSet a
Circumpolar
toRiseSetLCT DecimalDegrees
_ LocalCivilDate
_ RiseSetLST
NeverRises = RiseSetLCT
forall a. RiseSet a
NeverRises


-- | Convert LST in decimal hours to the JuliadDate

-- the second parameter must be desired day at midnignt.

dhToJD :: DecimalHours -> JulianDate -> JulianDate
dhToJD :: DecimalHours -> JulianDate -> JulianDate
dhToJD (DH Double
hours) JulianDate
day = JulianDate
day JulianDate -> JulianDate -> JulianDate
forall a. Num a => a -> a -> a
+ (Double -> JulianDate
JD (Double -> JulianDate) -> Double -> JulianDate
forall a b. (a -> b) -> a -> b
$ Double
hoursDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
24)