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(..))
data RiseSet a
= RiseSet a a
| Circumpolar
| NeverRises
deriving (Show, Eq)
type RSInfo a = (a, DecimalDegrees)
type RiseSetLST = RiseSet (RSInfo LocalSiderealTime)
type RiseSetLCT = RiseSet (RSInfo LocalCivilTime)
type RiseSetMB = RiseSet (Maybe (RSInfo LocalCivilTime))
riseAndSet :: EquatorialCoordinates1 -> DecimalDegrees -> DecimalDegrees -> RiseSetLST
riseAndSet (EC1 delta alpha) shift lat =
let delta' = toRadians delta
shift' = toRadians shift
lat' = toRadians lat
cosH = cosOfHourAngle delta' shift' lat'
in sortRiseSet cosH delta' shift' lat'
where sortRiseSet :: Double -> Double -> Double -> Double -> RiseSetLST
sortRiseSet cosH delta shift latitude
| cosH < -1 = Circumpolar
| cosH > 1 = NeverRises
| otherwise = calcTimesAndAzimuths alpha (toHours $ acos cosH) delta shift latitude
toHours :: Double -> DecimalHours
toHours = toDecimalHours . fromRadians
cosOfHourAngle :: Double -> Double -> Double -> Double
cosOfHourAngle delta shift latitude = -((sin shift) + (sin latitude)*(sin delta)) / ((cos latitude)*(cos delta))
calcTimesAndAzimuths :: DecimalHours -> DecimalHours -> Double -> Double -> Double -> RiseSetLST
calcTimesAndAzimuths alpha hourAngle delta shift latitude =
let lstRise = dhToLST $ reduceToZeroRange 24 $ alpha - hourAngle
lstSet = dhToLST $ reduceToZeroRange 24 $ alpha + hourAngle
azimuthRise = reduceToZeroRange (2*pi) $ acos $ ((sin delta) + (sin shift)*(sin latitude)) / ((cos shift)*(cos latitude))
azimuthSet = 2*pi - azimuthRise
in RiseSet (lstRise, fromRadians azimuthRise) (lstSet, fromRadians azimuthSet)
riseAndSet2 :: DecimalHours
-> (JulianDate -> EquatorialCoordinates1)
-> GeographicCoordinates
-> DecimalDegrees
-> LocalCivilDate
-> RiseSetMB
riseAndSet2 eps getPosition geoc shift lcd =
let day = lcdDate lcd
pos = getPosition (addHours 12 day)
rs = riseAndSetLCT geoc lcd shift pos
rise = calc getRiseTime (getRiseTime rs) 0
set = calc getSetTime (getSetTime rs) 0
in case rs of
Circumpolar -> Circumpolar
NeverRises -> NeverRises
_ -> buildResult rise set
where calc :: (RiseSetLCT -> RSInfo LocalCivilTime) -> RSInfo LocalCivilTime -> Int -> RiseSetLCT
calc getRSInfo rsi@(time, _) iterNo =
let pos = getPosition $ lctUniversalTime time
rs = riseAndSetLCT geoc lcd shift pos
rsi' = getRSInfo rs
in case rs of
Circumpolar -> Circumpolar
NeverRises -> NeverRises
_ -> if isOK rsi rsi' || iterNo >= maxIters
then rs
else calc getRSInfo rsi' (iterNo+1)
isOK :: RSInfo LocalCivilTime -> RSInfo LocalCivilTime -> Bool
isOK (t1, _) (t2, _) = (abs d) < (h/24)
where JD d = (lctUniversalTime t1) - (lctUniversalTime t2)
DH h = eps
maxIters = 3
getRiseTime :: RiseSetLCT -> RSInfo LocalCivilTime
getRiseTime (RiseSet r _) = r
getSetTime :: RiseSetLCT -> RSInfo LocalCivilTime
getSetTime (RiseSet _ s) = s
buildResult (RiseSet r _) (RiseSet _ s) = RiseSet (Just r) (Just s)
buildResult (RiseSet r _) _ = RiseSet (Just r) Nothing
buildResult _ (RiseSet _ s) = RiseSet Nothing (Just s)
riseAndSetLCT :: GeographicCoordinates
-> LocalCivilDate
-> DecimalDegrees
-> EquatorialCoordinates1
-> RiseSetLCT
riseAndSetLCT (GeoC latitude longitude) lcd shift ec
= toRiseSetLCT longitude lcd $ riseAndSet ec shift latitude
toRiseSetLCT :: DecimalDegrees
-> LocalCivilDate
-> RiseSetLST
-> RiseSetLCT
toRiseSetLCT longitude lcd (RiseSet (rise, azRise) (set, azSet)) =
let toLCT lst = lstToLCT longitude lcd lst
rise' = toLCT rise
set' = toLCT set
in RiseSet (rise', azRise) (set', azSet)
toRiseSetLCT _ _ Circumpolar = Circumpolar
toRiseSetLCT _ _ NeverRises = NeverRises
dhToJD :: DecimalHours -> JulianDate -> JulianDate
dhToJD (DH hours) day = day + (JD $ hours/24)