{-| 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 (Show, 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 (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) -- | 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 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) -- | 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 (GeoC latitude longitude) lcd shift ec = toRiseSetLCT longitude lcd $ riseAndSet ec shift 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 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 -- | Convert LST in decimal hours to the JuliadDate -- the second parameter must be desired day at midnignt. dhToJD :: DecimalHours -> JulianDate -> JulianDate dhToJD (DH hours) day = day + (JD $ hours/24)