packages feed

astro (empty) → 0.4.1.0

raw patch · 29 files changed

+3422/−0 lines, 29 filesdep +HUnitdep +QuickCheckdep +aesonsetup-changed

Dependencies added: HUnit, QuickCheck, aeson, astro, base, bytestring, matrix, optparse-applicative, test-framework, test-framework-hunit, test-framework-quickcheck2, time

Files

+ LICENSE view
@@ -0,0 +1,30 @@+Copyright Alexander Ignatyev (c) 2016++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++    * Redistributions of source code must retain the above copyright+      notice, this list of conditions and the following disclaimer.++    * Redistributions in binary form must reproduce the above+      copyright notice, this list of conditions and the following+      disclaimer in the documentation and/or other materials provided+      with the distribution.++    * Neither the name of Author name here nor the names of other+      contributors may be used to endorse or promote products derived+      from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ app/Main.hs view
@@ -0,0 +1,314 @@+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE OverloadedStrings #-}+module Main where++import GHC.Generics+import Data.Aeson+import qualified Data.ByteString.Lazy.Char8 as B+import Data.Maybe (fromMaybe)+import Data.Time.LocalTime (ZonedTime, getZonedTime)+import Options.Applicative+import Data.Monoid((<>))++-- Astro Imports+import Data.Astro.Time.JulianDate+import Data.Astro.Time.Conv (zonedTimeToLCT, zonedTimeToLCD, lctToZonedTime)++import Data.Astro.Effects (refract)+import Data.Astro.CelestialObject.RiseSet(riseAndSetLCT, riseAndSet2, RiseSetMB(..), RiseSetLCT(..))++import Data.Astro.Sun++import Data.Astro.Star++import Data.Astro.Types+import Data.Astro.Coordinate++import Data.Astro.Moon (moonPosition1, moonDistance1, moonAngularSize)+import Data.Astro.Moon.MoonDetails (j2010MoonDetails, mduToKm)++import Data.Astro.Planet (Planet(..), planetPosition, planetTrueAnomaly1, planetDistance1, planetAngularDiameter)+import Data.Astro.Planet.PlanetDetails (j2010PlanetDetails)+++main :: IO ()+main = execParser opts >>= run+  where opts = info (cmdOptions <**> helper)+          ( progDesc "Amateur astronomical computations"+            <> header "Astro" )+++run :: CmdOptions -> IO ()+run cmdOptions = do+  defParams <- defaultParams+  let params = fromMaybe defParams $ fromMaybe defParams <$> decode <$> B.pack <$> cmdJson cmdOptions+      res = processQuery params+  B.putStrLn $ encode res+++-- Calcs+calculateSunResult :: Params -> PlanetaiResult+calculateSunResult params = PR {+  riseSet = riseSet+  , distance = DR distance "km"+  , angularSize = angularSize+  , position = hcPosition+  }+  where coords = paramsCoordinates params+        date = paramsDate params+        lct = paramsDateTime params+        jd = lctUniversalTime lct+        rs = sunRiseAndSet coords 0.833333 date+        riseSet = toRiseSetResult rs+        distance = sunDistance jd+        DD angularSize = sunAngularSize jd+        ec1 = sunPosition2 jd+        hcPosition = toHorizonCoordinatesResult coords jd ec1+++calculateMoonResult :: Params -> PlanetaiResult+calculateMoonResult params = PR {+  riseSet = riseSet+  , distance = DR distance "km"+  , angularSize = angularSize+  , position = hcPosition+  }+  where position = moonPosition1 j2010MoonDetails+        coords = paramsCoordinates params+        verticalShift = refract (DD 0) 12 1012+        date = paramsDate params+        lct = paramsDateTime params+        jd = lctUniversalTime lct+        rs = riseAndSet2 0.000001 position coords verticalShift date+        riseSet = toRiseSetResult rs+        mdu = moonDistance1 j2010MoonDetails jd+        distance = mduToKm mdu+        DD angularSize = moonAngularSize mdu+        ec1 = position jd+        hcPosition = toHorizonCoordinatesResult coords jd ec1+++calculatePlanetResult :: Params -> Planet -> PlanetaiResult+calculatePlanetResult params planet = PR {+  riseSet = riseSet+  , distance = DR distance "AU"+  , angularSize = angularSize+  , position = hcPosition+  }+  where coords = paramsCoordinates params+        verticalShift = refract (DD 0) 12 1012+        date = paramsDate params+        lct = paramsDateTime params+        jd = lctUniversalTime lct+        planetDetails = j2010PlanetDetails planet+        earthDetails = j2010PlanetDetails Earth+        position = planetPosition planetTrueAnomaly1 planetDetails earthDetails+        rs = riseAndSet2 0.000001 position coords verticalShift date+        riseSet = toRiseSetResult rs+        au = planetDistance1 planetDetails earthDetails jd+        AU distance = au+        DD angularSize = planetAngularDiameter planetDetails au+        ec1 = position jd+        hcPosition = toHorizonCoordinatesResult coords jd ec1+++calculateStarResult :: Params -> Star -> StarResult+calculateStarResult params star = SR {+  starRiseSet = riseSet+  , starPosition = hcPosition+  }+  where coords = paramsCoordinates params+        verticalShift = refract (DD 0) 12 1012+        date = paramsDate params+        lct = paramsDateTime params+        jd = lctUniversalTime lct+        ec1 = starCoordinates star+        rs = riseAndSetLCT coords date verticalShift ec1+        riseSet = fromRiseSetLCT rs+        hcPosition = toHorizonCoordinatesResult coords jd ec1+++toRiseSetResult :: RiseSetMB -> RiseSetResult+toRiseSetResult rs = case rs of+  RiseSet rise set -> RSR { rise = lctToZonedTime <$> fst <$> rise+                          , riseAzimuth = ddValue <$> snd <$> rise+                          , set = lctToZonedTime <$> fst <$> set+                          , setAzimuth = ddValue <$> snd <$> set+                          , state = "Rise and/or set"+                          }+  Circumpolar -> RSR Nothing Nothing Nothing Nothing "Circumpolar"+  NeverRises -> RSR Nothing Nothing Nothing Nothing "NeverRises"+++fromRiseSetLCT :: RiseSetLCT -> RiseSetResult+fromRiseSetLCT rs = case rs of+  RiseSet rise set -> RSR { rise = Just $ lctToZonedTime $ fst rise+                          , riseAzimuth = Just $ ddValue $ snd $ rise+                          , set = Just $ lctToZonedTime $ fst set+                          , setAzimuth = Just $ ddValue $ snd $ set+                          , state = "Rise and Set"+                          }+  Circumpolar -> RSR Nothing Nothing Nothing Nothing "Circumpolar"+  NeverRises -> RSR Nothing Nothing Nothing Nothing "NeverRises"+++ddValue :: DecimalDegrees -> Double+ddValue (DD value) = value++toHorizonCoordinatesResult :: GeographicCoordinates+                           -> JulianDate+                           -> EquatorialCoordinates1+                           -> HorizonCoordinatesResult+toHorizonCoordinatesResult (GeoC lat long) jd (EC1 delta alpha) = HCR altitude azimuth+  where ec2 = EC2 delta (raToHA alpha long jd)+        hc = equatorialToHorizon lat ec2+        HC (DD altitude) (DD azimuth) = hc+        +        +        ++processQuery :: Params -> AstroResult+processQuery params = AstroResult {+  request = params+  , sun = calculateSunResult params+  , moon = calculateMoonResult params+  , mercury = calculatePlanetResult params Mercury+  , venus = calculatePlanetResult params Venus+  , mars = calculatePlanetResult params Mars+  , jupiter = calculatePlanetResult params Jupiter+  , saturn = calculatePlanetResult params Saturn+  , uranus = calculatePlanetResult params Uranus+  , neptune = calculatePlanetResult params Neptune+  , polaris = calculateStarResult params Polaris+  , alphaCrucis = calculateStarResult params AlphaCrucis+  , sirius = calculateStarResult params Sirius+  , betelgeuse = calculateStarResult params Betelgeuse+  , rigel = calculateStarResult params Rigel+  , vega = calculateStarResult params Vega+  , antares = calculateStarResult params Antares+  , canopus = calculateStarResult params Canopus+  , pleiades = calculateStarResult params Pleiades+  }+++-- Command Line Options+data CmdOptions = CmdOptions {+  cmdJson :: Maybe String+  }+++cmdOptions :: Parser CmdOptions+cmdOptions = CmdOptions+  <$> (optional $ strOption ( long "json" <> short 'j' <> help "JSON-encoded params") )+++-- Params+data CoordinatesParam = CoordinatesParam {+    latitude :: Double+  , longitude  :: Double+  } deriving (Generic, Show)++instance ToJSON CoordinatesParam+instance FromJSON CoordinatesParam+++data Params = Params {+  coordinates :: CoordinatesParam+  , datetime :: ZonedTime+  } deriving (Generic, Show)++instance ToJSON Params+instance FromJSON Params+++paramsCoordinates :: Params -> GeographicCoordinates+paramsCoordinates params = GeoC (DD $ latitude coords) (DD $ longitude coords)+  where coords = coordinates params+++paramsDateTime :: Params -> LocalCivilTime+paramsDateTime = zonedTimeToLCT . datetime+++paramsDate :: Params -> LocalCivilDate+paramsDate = zonedTimeToLCD . datetime+++greenwichCoordinates :: CoordinatesParam+greenwichCoordinates = CoordinatesParam 51.4768 0+++defaultParams :: IO (Params)+defaultParams = do+  time <- getZonedTime+  return Params {+    coordinates = greenwichCoordinates+    , datetime = time+    }+++-- Result+data HorizonCoordinatesResult = HCR {+  altitude :: Double+  , azimuth :: Double+  } deriving (Generic, Show)++instance ToJSON HorizonCoordinatesResult++data RiseSetResult = RSR {+  rise :: Maybe ZonedTime+  , riseAzimuth :: Maybe Double+  , set :: Maybe ZonedTime+  , setAzimuth :: Maybe Double+  , state :: String+  } deriving (Generic, Show)++instance ToJSON RiseSetResult+++data DistanceResult = DR {+  value :: Double+  , units :: String+  } deriving (Generic, Show)++instance ToJSON DistanceResult++data PlanetaiResult = PR {+  riseSet :: RiseSetResult+  , distance :: DistanceResult+  , angularSize:: Double+  , position :: HorizonCoordinatesResult             +  } deriving (Generic, Show)++instance ToJSON PlanetaiResult++data StarResult = SR {+  starRiseSet :: RiseSetResult+  , starPosition :: HorizonCoordinatesResult+  } deriving (Generic, Show)++instance ToJSON StarResult++data AstroResult = AstroResult {+  request :: Params+  , sun :: PlanetaiResult+  , moon :: PlanetaiResult+  , mercury :: PlanetaiResult+  , venus :: PlanetaiResult+  , mars :: PlanetaiResult+  , jupiter :: PlanetaiResult+  , saturn :: PlanetaiResult+  , uranus :: PlanetaiResult+  , neptune :: PlanetaiResult+  , polaris :: StarResult+  , alphaCrucis :: StarResult+  , sirius :: StarResult+  , betelgeuse :: StarResult+  , rigel :: StarResult+  , vega :: StarResult+  , antares :: StarResult+  , canopus :: StarResult+  , pleiades :: StarResult+  } deriving (Generic, Show)++instance ToJSON AstroResult
+ astro.cabal view
@@ -0,0 +1,76 @@+name:                astro+version:             0.4.1.0+synopsis:            Astro+description:         Please see README.md+homepage:            https://github.com/alexander-ignatyev/astro+license:             BSD3+license-file:        LICENSE+author:              Alexander Ignatyev+maintainer:          Alexander Ignatyev+copyright:           2016-2017 Alexander Ignatyev+category:            Science+build-type:          Simple+-- extra-source-files:+cabal-version:       >=1.10++library+  hs-source-dirs:      src+  exposed-modules:     Data.Astro.Time+                     , Data.Astro.Time.GregorianCalendar+                     , Data.Astro.Time.JulianDate+                     , Data.Astro.Time.Sidereal+                     , Data.Astro.Time.Epoch+                     , Data.Astro.Time.Conv+                     , Data.Astro.Coordinate+                     , Data.Astro.Types+                     , Data.Astro.Utils+                     , Data.Astro.CelestialObject+                     , Data.Astro.CelestialObject.RiseSet+                     , Data.Astro.Effects+                     , Data.Astro.Effects.Parallax+                     , Data.Astro.Star+                     , Data.Astro.Sun+                     , Data.Astro.Sun.SunInternals+                     , Data.Astro.Planet+                     , Data.Astro.Planet.PlanetDetails+                     , Data.Astro.Planet.PlanetMechanics+                     , Data.Astro.Moon+                     , Data.Astro.Moon.MoonDetails+  other-modules:       Data.Astro.Effects.Precession+                     , Data.Astro.Effects.Nutation+                     , Data.Astro.Effects.Aberration+  build-depends:       base >= 4.7 && < 5+                     , time+                     , matrix+  default-language:    Haskell2010++executable astro-app+  hs-source-dirs:      app+  main-is:             Main.hs+  ghc-options:         -threaded -rtsopts -with-rtsopts=-N+  build-depends:       base+                     , bytestring+                     , time+                     , aeson+                     , optparse-applicative+                     , astro+  default-language:    Haskell2010++test-suite astro-test+  type:                exitcode-stdio-1.0+  hs-source-dirs:      test+  main-is:             Main.hs+  build-depends:       base+                     , astro+                     , time+                     , test-framework+                     , test-framework-hunit+                     , test-framework-quickcheck2+                     , HUnit+                     , QuickCheck > 2.0+  ghc-options:         -threaded -rtsopts -with-rtsopts=-N+  default-language:    Haskell2010++source-repository head+  type:     git+  location: https://github.com/alexander-ignatyev/astro
+ src/Data/Astro/CelestialObject.hs view
@@ -0,0 +1,42 @@+{-|+Module: Data.Astro.CelestialObject+Description: Computations characteristics of selestial objects+Copyright: Alexander Ignatyev, 2016++Computations characteristics of selestial objects.+-}++module Data.Astro.CelestialObject+(+   angleEquatorial+  , angleEcliptic+)++where++import Data.Astro.Types (DecimalDegrees, toRadians, fromRadians, fromDecimalHours)+import Data.Astro.Coordinate (EquatorialCoordinates1(..), EclipticCoordinates(..))+++-- | Calculate angle between two celestial objects+-- whose coordinates specified in Equatorial Coordinate System.+angleEquatorial :: EquatorialCoordinates1 -> EquatorialCoordinates1 -> DecimalDegrees+angleEquatorial (EC1 delta1 alpha1) (EC1 delta2 alpha2) =+  calcAngle (delta1, fromDecimalHours alpha1) (delta2, fromDecimalHours alpha2)+++-- | Calculate angle between two celestial objects+-- whose coordinates specified in Ecliptic Coordinate System.+angleEcliptic :: EclipticCoordinates -> EclipticCoordinates -> DecimalDegrees+angleEcliptic (EcC beta1 lambda1) (EcC beta2 lambda2) =+  calcAngle (beta1, lambda1) (beta2, lambda2)+++calcAngle :: (DecimalDegrees, DecimalDegrees) -> (DecimalDegrees, DecimalDegrees) -> DecimalDegrees+calcAngle (up1, round1) (up2, round2) =+  let up1' = toRadians up1+      round1' = toRadians round1+      up2' = toRadians up2+      round2' = toRadians round2+      d = acos $ (sin up1')*(sin up2') + (cos up1')*(cos up2')*cos(round1'-round2')+  in fromRadians d
+ src/Data/Astro/CelestialObject/RiseSet.hs view
@@ -0,0 +1,191 @@+{-|+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)
+ src/Data/Astro/Coordinate.hs view
@@ -0,0 +1,362 @@+{-|+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, &#x3B1;__ - /'how far up'/ angle from the horizontal plane in degrees+* __azimuth,  &#x391;__ - /'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, &#x3B4;__ - /'how far up'/ angle from the quatorial plane;+* __right ascension, &#x3B1;__  - /'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, &#x3B2;__ - /'how far up'/ angle from the ecliptic+* __ecliptic longitude, &#x3BB;__ - /'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__, &#x2648; - fixed direction lies along the line of the intersection of the equatorial plane and the ecliptic+* __obliquity of the ecliptic, &#x3B2;__ - 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 {+  hAltitude :: DecimalDegrees   -- ^ alpha+  , hAzimuth :: DecimalDegrees  -- ^ big alpha+  } deriving (Show, Eq)+++-- | Equatorial Coordinates, defines fixed position in the sky+data EquatorialCoordinates1 = EC1 {+  e1Declination :: DecimalDegrees     -- ^ delta+  , e1RightAscension :: DecimalHours  -- ^ alpha+  } deriving (Show, Eq)+++-- | Equatorial Coordinates+data EquatorialCoordinates2 = EC2 {+  e2Declination :: DecimalDegrees    -- ^ delta+  , e2HoursAngle :: DecimalHours     -- ^ H+  } deriving (Show, Eq)+++-- | Ecliptic Coordinates+data EclipticCoordinates = EcC {+  ecLatitude :: DecimalDegrees      -- ^ beta+  , ecLongitude :: DecimalDegrees   -- ^ lambda+  } deriving (Show, Eq)+++-- | Galactic Coordinates+data GalacticCoordinates = GC {+  gLatitude :: DecimalDegrees       -- ^ b+  , gLongitude :: DecimalDegrees    -- ^ l+  } deriving (Show, Eq)+++-- | Convert Right Ascension to Hour Angle for specified longitude and Universal Time+raToHA :: DecimalHours -> DecimalDegrees -> JulianDate -> DecimalHours+raToHA = haRAConv+++-- | Convert Hour Angle to Right Ascension for specified longitude and Universal Time+haToRA :: DecimalHours -> DecimalDegrees -> JulianDate -> DecimalHours+haToRA = haRAConv+++-- | HA <-> RA Conversions+haRAConv :: DecimalHours -> DecimalDegrees -> JulianDate -> DecimalHours+haRAConv dh longitude ut =+  let lst = utToLST longitude ut  -- Local Sidereal Time+      DH hourAngle = (lstToDH lst) - dh+  in if hourAngle < 0 then (DH $ hourAngle+24) else (DH 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 latitude (EC2 dec hourAngle) =+  let hourAngle' = fromDecimalHours hourAngle+      (altitude, azimuth) = ecHCConv latitude (dec, hourAngle')+  in HC altitude 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 latitude (HC altitude azimuth) =+  let (dec, hourAngle) = ecHCConv latitude (altitude, azimuth)+  in EC2 dec $ toDecimalHours 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 (GeoC latitude longitude) jd (EC1 delta alpha) =+  let ec2 = EC2 delta (raToHA alpha longitude jd)+  in equatorialToHorizon latitude 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 (GeoC latitude longitude) jd hc =+  let (EC2 dec hourAngle) = horizonToEquatorial latitude hc+  in EC1 dec (haToRA hourAngle longitude 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 latitude (up, round) =+  let latitude' = toRadians latitude+      up' = toRadians up+      round' = toRadians round+      sinUpResult = (sin up')*(sin latitude') + (cos up')*(cos latitude')*(cos round')+      upResult = asin sinUpResult+      roundResult = acos $ ((sin up') - (sin latitude')*sinUpResult) / ((cos latitude') * (cos upResult))+      roundResult' = if (sin round') < 0 then roundResult else (2*pi - roundResult)+  in ((fromRadians upResult), (fromRadians roundResult'))+++-- | Calculate the obliquity of the ecpliptic on JulianDate+obliquity :: JulianDate -> DecimalDegrees+obliquity jd =+  let DD baseObliquity = fromDMS 23 26 21.45+      t = numberOfCenturies j2000 jd+      de = (46.815*t + 0.0006*t*t - 0.00181*t*t*t) / 3600  -- 3600 number of seconds in 1 degree+  in (DD $ baseObliquity - de) + (nutationObliquity jd)+++-- | Converts Ecliptic Coordinates on specified Julian Date to Equatorial Coordinates+eclipticToEquatorial :: EclipticCoordinates -> JulianDate -> EquatorialCoordinates1+eclipticToEquatorial (EcC beta gamma) jd =+  let epsilon' = toRadians $ obliquity jd+      beta' = toRadians beta+      gamma' = toRadians gamma+      delta = asin $ (sin beta')*(cos epsilon') + (cos beta')*(sin epsilon')*(sin gamma')+      y = (sin gamma')*(cos epsilon') - (tan beta')*(sin epsilon')+      x = cos gamma'+      alpha = reduceToZero2PI $ atan2 y x+  in EC1 (fromRadians delta) (toDecimalHours $ fromRadians alpha)+++-- | Converts Equatorial Coordinates to Ecliptic Coordinates on specified Julian Date+equatorialToEcliptic :: EquatorialCoordinates1 -> JulianDate -> EclipticCoordinates+equatorialToEcliptic (EC1 delta alpha) jd =+  let epsilon' = toRadians $ obliquity jd+      delta' = toRadians delta+      alpha' = toRadians $ fromDecimalHours alpha+      beta = asin $ (sin delta')*(cos epsilon') - (cos delta')*(sin epsilon')*(sin alpha')+      y = (sin alpha')*(cos epsilon') + (tan delta')*(sin epsilon')+      x = cos alpha'+      gamma = reduceToZero2PI $ atan2 y x+  in EcC (fromRadians beta) (fromRadians gamma)+++-- | Galactic Pole Coordinates+galacticPole :: EquatorialCoordinates1+galacticPole = EC1 (DD 27.4) (toDecimalHours $ DD 192.25)++galacticPoleInRadians = (delta, alpha)+  where delta = toRadians $ e1Declination galacticPole+        alpha = toRadians $ fromDecimalHours $ e1RightAscension galacticPole+++-- | Ascending node of the galactic place on equator+ascendingNode :: DecimalDegrees+ascendingNode = DD 33+++-- | Convert Galactic Coordinates Equatorial Coordinates+galacticToEquatorial :: GalacticCoordinates -> EquatorialCoordinates1+galacticToEquatorial (GC b l) =+  let b' = toRadians b+      l' = toRadians l+      (poleDelta, poleAlpha) = galacticPoleInRadians+      an = toRadians ascendingNode+      delta = asin $ (cos b')*(cos poleDelta)*(sin (l'-an)) + (sin b')*(sin poleDelta)+      y = (cos b')*(cos (l'-an))+      x = (sin b')*(cos poleDelta) - (cos b')*(sin poleDelta)*(sin (l'-an))+      alpha = reduceToZero2PI $ (atan2 y x) + poleAlpha+  in EC1 (fromRadians delta) (toDecimalHours $ fromRadians alpha)+++-- | Convert Equatorial Coordinates to Galactic Coordinates+equatorialToGalactic :: EquatorialCoordinates1 -> GalacticCoordinates+equatorialToGalactic (EC1 delta alpha) =+  let delta' = toRadians delta+      alpha' = toRadians $ fromDecimalHours alpha+      (poleDelta, poleAlpha) = galacticPoleInRadians+      sinb = (cos delta')*(cos poleDelta)*(cos (alpha'-poleAlpha)) + (sin delta') * (sin poleDelta)+      y = (sin delta') - sinb*(sin poleDelta)+      x = (cos delta')*(sin (alpha'-poleAlpha))*(cos poleDelta)+      b = asin sinb+      l = reduceToZero2PI $ (atan2 y x) + (toRadians ascendingNode)+  in GC (fromRadians b) (fromRadians 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 rad = if rad < 0 then rad + 2*pi else rad
+ src/Data/Astro/Effects.hs view
@@ -0,0 +1,50 @@+{-|+Module: Data.Astro.Effects+Description: Physical effects+Copyright: Alexander Ignatyev, 2016++Physical effects which influence on accuracy of astronomical calculations.+-}++module Data.Astro.Effects+(+  refract+  , Precession.AstronomyEpoch(..)+  , Precession.precession1+  , Precession.precession2+  , Nutation.nutationLongitude+  , Nutation.nutationObliquity+  , Aberration.includeAberration+  , Parallax.parallax+)++where++import Data.Astro.Types (DecimalDegrees(..), toRadians)++import qualified Data.Astro.Effects.Precession as Precession+import qualified Data.Astro.Effects.Nutation as Nutation+import qualified Data.Astro.Effects.Aberration as Aberration+import qualified Data.Astro.Effects.Parallax as Parallax++-- | Calculate the atmospheric refraction angle.+-- It takes the observed altitude (of Horizon Coordinates), temperature in degrees centigrade and barometric pressure in millibars.+-- The average sea level atmospheric pressure is 1013 millibars.+refract :: DecimalDegrees -> Double -> Double -> DecimalDegrees+refract altitude temperature pressure =+  let f = if altitude > (DD 15) then refractBigAlpha else refractSmallAlpha+  in f altitude temperature pressure+++-- | Calculate the atmospheric refraction angle for big values of alpha (altitude) (> 15 decimal degrees)+refractBigAlpha :: DecimalDegrees -> Double -> Double -> DecimalDegrees+refractBigAlpha altitude temperature pressure =+  let z = toRadians $ 90 - altitude  -- zenith angle+  in DD $ 0.00452*pressure*(tan z) /(273+temperature) +++-- | Calculate the atmospheric refraction angle for small values of alpha (altitude)+refractSmallAlpha :: DecimalDegrees -> Double -> Double -> DecimalDegrees+refractSmallAlpha altitude temperature pressure =+  let a = toRadians altitude+  in DD $ pressure*(0.1594+0.0196*a+0.00002*a*a)/((273+temperature)*(1+0.505*a+0.0845*a*a))
+ src/Data/Astro/Effects/Aberration.hs view
@@ -0,0 +1,32 @@+{-|+Module: Data.Astro.Effects.Aberration+Description: Calculation effects of aberration.+Copyright: Alexander Ignatyev, 2016++Calculation effects of aberration.+-}++module Data.Astro.Effects.Aberration+(+  includeAberration+)++where++import Data.Astro.Types (DecimalDegrees, toRadians, fromDMS)+import Data.Astro.Time.JulianDate (JulianDate)+import Data.Astro.Coordinate (EclipticCoordinates(..))+++-- | Includes aberration effect.+-- It takes true Ecliptic Coordinates,+-- the Sun's longitude at the given Julian Day (the third parameter).+-- Returns apparent ecliptic coordinates.+-- The Sun's longitude can be calculated using 'sunEclipticLongitude1' or 'sunEclipticLongitude2' of "Data.Astro.Sun" module.+includeAberration :: EclipticCoordinates -> JulianDate -> DecimalDegrees -> EclipticCoordinates+includeAberration (EcC beta lambda) jd sunLambda =+  let lambdaDiff = toRadians $ sunLambda - lambda+      beta' = toRadians beta+      dLambda = -20.5 * (cos lambdaDiff) / (cos beta')+      dBeta = -20.5 * (sin lambdaDiff) * (sin beta')+  in EcC (beta + fromDMS 0 0 dBeta) (lambda + fromDMS 0 0 dLambda)
+ src/Data/Astro/Effects/Nutation.hs view
@@ -0,0 +1,59 @@+{-|+Module: Data.Astro.Effects.Nutation+Description: Calculation effects of nutation+Copyright: Alexander Ignatyev, 2016++Calculation effects of nutation.+-}++module Data.Astro.Effects.Nutation+(+  nutationLongitude+  , nutationObliquity+)++where++import qualified Data.Astro.Utils as U+import Data.Astro.Types (DecimalDegrees(..), toRadians, fromDMS)+import Data.Astro.Time.JulianDate (JulianDate, numberOfCenturies)+import Data.Astro.Time.Epoch (j1900)+++-- | Calculates the nutation on the ecliptic longitude at the given JulianDate+nutationLongitude :: JulianDate -> DecimalDegrees+nutationLongitude jd =+  let t = numberOfCenturies j1900 jd+      l = sunMeanLongutude t+      omega = moonNode t+      dPsi = -17.2*(sin omega) - 1.3*(sin $ 2*l)+  in fromDMS 0 0 dPsi+++-- | Calculates the nutation on the obliquity of the ecliptic at the given JulianDate+nutationObliquity :: JulianDate -> DecimalDegrees+nutationObliquity jd =+  let t = numberOfCenturies j1900 jd+      l = sunMeanLongutude t+      omega = moonNode t+      dEps = 9.2*(cos omega) + 0.5*(cos $ 2*l)+  in fromDMS 0 0 dEps+++-- | It takes a number of centuries and returns the Sun's mean longitude in radians+sunMeanLongutude :: Double -> Double+sunMeanLongutude t =+  let a = 100.002136 * t+  in U.toRadians $ U.reduceToZeroRange 360 $ 279.6967 + 360 * (a - int a)+++-- | It takes a number of centuries and returns the Moon's node in radians+moonNode :: Double -> Double+moonNode t =+  let b = 5.372617 * t+  in U.toRadians $ U.reduceToZeroRange 360 $ 259.1833 - 360*(b - int b)+++-- | 'round' function that returns Double+int :: Double -> Double+int = fromIntegral . round
+ src/Data/Astro/Effects/Parallax.hs view
@@ -0,0 +1,70 @@+{-|+Module: Data.Astro.Effects.Parallax+Description: Calculation effects of geocentric parallax+Copyright: Alexander Ignatyev, 2016+++Calculation effects of geocentric parallax.+-}++module Data.Astro.Effects.Parallax+(+  parallaxQuantities+  , parallax+)++where++import Data.Astro.Types (DecimalDegrees(..)+                        , DecimalHours(..)+                        , AstronomicalUnits(..)+                        , GeographicCoordinates(..)+                        , toRadians, fromRadians+                        , fromDMS+                        , toDecimalHours, fromDecimalHours)+import Data.Astro.Time (utToLST)+import Data.Astro.Time.JulianDate (JulianDate(..))+import Data.Astro.Time.Sidereal (LocalSiderealTime(..), utToGST, gstToLST)+import Data.Astro.Coordinate (EquatorialCoordinates1(..), raToHA)+++-- | It takes latitude of the observer+-- and height above sea-level of the observer measured in metres+-- Returns palallax quantities (p*(sin phi'), p*(cos phi')),+-- where phi' is the geocentric latitude+-- and p is the distance of the obserbve from the centre of the Earth.+parallaxQuantities :: DecimalDegrees -> Double -> (Double, Double)+parallaxQuantities latitude height =+  let c = 0.996647+      phi = toRadians latitude+      h = earthRadiusUnits height+      u = atan (c*(tan phi))+      pSin = c * (sin u) + h*(sin phi)+      pCos = (cos u) + h*(cos phi)+  in (pSin, pCos)+++-- | Calculate the apparent position of the celestial object (the Sun or a planet).+-- It takes geocraphic coordinates of the observer and height above sea-level of the observer measured in metres,+-- distance from the celestial object to the Earth measured in AU, the Universal Time and geocentric equatorial coordinates.+-- It returns adjusted equatorial coordinates.+parallax :: GeographicCoordinates -> Double -> AstronomicalUnits -> JulianDate -> EquatorialCoordinates1 -> EquatorialCoordinates1+parallax (GeoC latitude longitude) height distance ut (EC1 delta alpha) =+  let piD = earthRadiusUnitsAU distance+      lst = utToLST longitude ut+      (pSin, pCos) = parallaxQuantities latitude height+      ha = toRadians $ fromDecimalHours $ raToHA alpha longitude ut+      delta' = toRadians delta+      dAlpha = (toDecimalHours piD) * (DH $ (sin ha)*pCos/(cos delta'))+      dDelta = piD * (DD $ pSin*(cos delta') - pCos*(cos ha)*(sin delta'))+  in EC1 (delta-dDelta) (alpha-dAlpha)+++-- | It takes the distance in metres and+-- returns the distance measured in units of qquatorial Earth radius+earthRadiusUnits :: Double -> Double+earthRadiusUnits d = d / 6378140+++--earthRadiusUnitsAU :: AstronomicalUnits -> DecimalDegrees+earthRadiusUnitsAU (AU d) = fromDMS 0 0 (8.794/d)
+ src/Data/Astro/Effects/Precession.hs view
@@ -0,0 +1,117 @@+{-|+Module: Data.Astro.Effects.Precession+Description: Luni-solar precession+Copyright: Alexander Ignatyev, 2016++Luni-solar precession.+-}++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(..))+++-------------------------------------------------------------------------------+-- Low-precision Precession++-- | Epoch Enumeration. See also "Data.Astro.Time.JulianDate" module.+data AstronomyEpoch = B1900  -- ^ Epoch B1900.0+                    | B1950  -- ^ Epoch B1950.0+                    | J2000  -- ^ Epoch J2000.0+                    | J2050  -- ^ Epoch J2050.0+                    deriving (Show, Eq)+++-- | Get the start date of the specified Epoch.+epochToJD :: AstronomyEpoch -> JulianDate+epochToJD B1900 = b1900+epochToJD B1950 = b1950+epochToJD J2000 = j2000+epochToJD J2050 = j2050+++-- | Precisional Constants+data PrecessionalConstants = PrecessionalConstants {+  pcM :: Double     -- ^ seconds+  , pcN :: Double   -- ^ seconds+  , pcN' :: Double  -- ^ arcsec+  }+++-- | Get Precision Constants of the Epoch+precessionalConstants :: AstronomyEpoch -> PrecessionalConstants+precessionalConstants B1900 = PrecessionalConstants 3.07234 1.33645 20.0468+precessionalConstants B1950 = PrecessionalConstants 3.07327 1.33617 20.0426+precessionalConstants J2000 = PrecessionalConstants 3.07420 1.33589 20.0383+precessionalConstants J2050 = PrecessionalConstants 3.07513 1.33560 20.0340+++-- | Low-precision method to calculate luni-solar precession.+-- It takes Epoch, Equatorial Coordinates those correct at the given epoch, Julian Date of the observation.+-- It returns corrected Equatorial Coordinates.+precession1 :: AstronomyEpoch -> EquatorialCoordinates1 -> JulianDate -> EquatorialCoordinates1+precession1 epoch (EC1 delta alpha) jd =+  let delta' = toRadians delta+      alpha' = toRadians $ fromDecimalHours alpha+      years = numberOfYears (epochToJD epoch) jd+      PrecessionalConstants m n n' = precessionalConstants epoch+      s1 = DH $ (m + n*(sin alpha')*(tan delta'))*years / 3600+      s2 = DD $ (n'*(cos alpha')) * years / 3600+  in (EC1 (delta + s2) (alpha + s1))+++-------------------------------------------------------------------------------+-- Rigorous Method+++-- | Rigorous method to calculate luni-solar precession.+-- It takes julian date at whose the coordinates are correct, Equatorial Coordinates, Julian Date of the observation.+-- It returns corrected Equatorial Coordinates.+precession2 :: JulianDate -> EquatorialCoordinates1 -> JulianDate -> EquatorialCoordinates1+precession2 epoch ec jd =+  let p' = prepareMatrixP' $ numberOfCenturies j2000 epoch+      v = prepareColumnVectorV ec+      p = transpose $ prepareMatrixP' $ numberOfCenturies j2000 jd+      [m, n, k] = toList $ p*(p'*v)+      alpha = atan2 n m+      delta = asin k+  in EC1 (fromRadians delta) (toDecimalHours $ fromRadians alpha)+++prepareMatrixP' n =+  let x = U.toRadians $ 0.6406161*n + 0.0000839*n*n + 0.0000050*n*n*n+      z = U.toRadians $ 0.6406161*n + 0.0003041*n*n + 0.0000051*n*n*n+      t = U.toRadians $ 0.5567530*n - 0.0001185*n*n - 0.0000116*n*n*n+      cx = cos x+      sx = sin x+      cz = cos z+      sz = sin z+      ct = cos t+      st = sin t+      matrix = [ [cx*ct*cz-sx*sz,    cx*ct*sz+sx*cz,    cx*st]+               , [(-sx)*ct*cz-cx*sz, (-sx)*ct*sz+cx*cz, (-sx)*st]+               , [(-st)*cz,          (-st)*sz,          ct] ]+  in fromLists matrix++prepareColumnVectorV (EC1 delta alpha) =+  let d = toRadians delta+      a = toRadians $ fromDecimalHours alpha+      cd = cos d+      sd = sin d+      ca = cos a+      sa = sin a+      v = [ca*cd, sa*cd, sd]+  in fromList 3 1 v
+ src/Data/Astro/Moon.hs view
@@ -0,0 +1,217 @@+{-|+Module: Data.Astro.Moon+Description: Calculation characteristics of the Moon+Copyright: Alexander Ignatyev, 2016++Calculation characteristics of the Moon.++= Example++@+import Data.Astro.Time.JulianDate+import Data.Astro.Coordinate+import Data.Astro.Types+import Data.Astro.Effects+import Data.Astro.CelestialObject.RiseSet+import Data.Astro.Moon++ro :: GeographicCoordinates+ro = GeoC (fromDMS 51 28 40) (-(fromDMS 0 0 5))++dt :: LocalCivilTime+dt = lctFromYMDHMS (DH 1) 2017 6 25 10 29 0++today :: LocalCivilDate+today = lcdFromYMD (DH 1) 2017 6 25++jd :: JulianDate+jd = lctUniversalTime dt++-- distance from the Earth to the Moon in kilometres+mdu :: MoonDistanceUnits+mdu = moonDistance1 j2010MoonDetails jd+-- MDU 0.9550170577020396++distance :: Double+distance = mduToKm mdu+-- 367109.51199772174++-- Angular Size+angularSize :: DecimalDegrees+angularSize = moonAngularSize mdu+-- DD 0.5425033990980761++-- The Moon's coordinates+position :: JulianDate -> EquatorialCoordinates1+position = moonPosition1 j2010MoonDetails++ec1 :: EquatorialCoordinates1+ec1 = position jd+-- EC1 {e1Declination = DD 18.706180658927323, e1RightAscension = DH 7.56710547682055}++hc :: HorizonCoordinates+hc = ec1ToHC ro jd ec1+-- HC {hAltitude = DD 34.57694951316064, hAzimuth = DD 103.91119101451832}++-- Rise and Set+riseSet :: RiseSetMB+riseSet = riseAndSet2 0.000001 position ro verticalShift today+-- RiseSet+--    (Just (2017-06-25 06:22:51.4858 +1.0,DD 57.81458864497365))+--    (Just (2017-06-25 22:28:20.3023 +1.0,DD 300.4168238905249))++-- Phase+phase :: Double+phase = moonPhase j2010MoonDetails jd+-- 2.4716141948212922e-2+++sunEC1 :: EquatorialCoordinates1+sunEC1 = sunPosition2 jd+-- EC1 {e1Declination = DD 23.37339098989099, e1RightAscension = DH 6.29262026252748}++limbAngle :: DecimalDegrees+limbAngle = moonBrightLimbPositionAngle ec1 sunEC1+-- DD 287.9869373767473+@+-}++module Data.Astro.Moon+(+  moonPosition1+  , moonDistance1+  , moonAngularSize+  , moonHorizontalParallax+  , moonPhase+  , moonBrightLimbPositionAngle+)++where++import qualified Data.Astro.Utils as U+import Data.Astro.Types (DecimalDegrees(..), toRadians, fromRadians)+import Data.Astro.Time.JulianDate (JulianDate(..), numberOfDays)+import Data.Astro.Coordinate (EquatorialCoordinates1(..), EclipticCoordinates(..), eclipticToEquatorial)+import Data.Astro.Planet (planetBrightLimbPositionAngle)+import Data.Astro.Sun (sunDetails, sunMeanAnomaly2, sunEclipticLongitude2)+import Data.Astro.Moon.MoonDetails (MoonDetails(..), MoonDistanceUnits(..), j2010MoonDetails)+++-- | Reduce the value to the range [0, 360)+reduceDegrees :: DecimalDegrees -> DecimalDegrees+reduceDegrees = U.reduceToZeroRange 360+++-- | Calculate Equatorial Coordinates of the Moon with the given MoonDetails and at the given JulianDate.+-- It is recommended to use 'j2010MoonDetails' as a first parameter.+moonPosition1 :: MoonDetails -> JulianDate -> EquatorialCoordinates1+moonPosition1 md ut =+  let sd = sunDetails ut+      lambdaS = sunEclipticLongitude2 sd+      ms = sunMeanAnomaly2 sd+      mmq = meanMoonQuantities md ut+      MQ lm'' _ nm' = correctedMoonQuantities lambdaS ms mmq+      a = toRadians $ lm''-nm'+      i = toRadians $ mdI md+      y = (sin a) * (cos i)+      x = cos a+      at = reduceDegrees $ fromRadians $ atan2 y x+      lambdaM = at + nm'+      betaM = fromRadians $ asin $ (sin a) * (sin i)+  in eclipticToEquatorial (EcC betaM lambdaM) ut+++-- | Calculates the Moon's Distance at the given julian date.+-- Returns distance to the Moon+-- moonDistance1 :: JulianDate -> MoonDistanceUnits+-- you can use 'mduToKm' (defined in "Data.Astro.Moon.MoonDetails") to convert result to kilometers+moonDistance1 :: MoonDetails -> JulianDate -> MoonDistanceUnits+moonDistance1 md ut =+  let sd = sunDetails ut+      lambdaS = sunEclipticLongitude2 sd+      ms = sunMeanAnomaly2 sd+      mmq = meanMoonQuantities md ut+      cmq = correctedMoonQuantities lambdaS ms mmq+      mm' = toRadians $ mqAnomaly cmq+      ec = toRadians $ centreEquation mm'+      e = mdE md+  in MDU $ (1 - e*e)/(1+e*(cos(mm'+ec)))+++-- | Calculate the Moon's angular size at the given distance.+moonAngularSize :: MoonDistanceUnits -> DecimalDegrees+moonAngularSize (MDU p) = (mdBigTheta j2010MoonDetails) / (DD p)+++-- | Calculates the Moon's horizontal parallax at the given distance.+moonHorizontalParallax :: MoonDistanceUnits -> DecimalDegrees+moonHorizontalParallax (MDU p) = (mdPi j2010MoonDetails) / (DD p)+++-- | Calculates the Moon's phase (the area of the visible segment expressed as a fraction of the whole disk)+-- at the given universal time.+moonPhase :: MoonDetails -> JulianDate -> Double+moonPhase md ut =+  let sd = sunDetails ut+      lambdaS = sunEclipticLongitude2 sd+      ms = sunMeanAnomaly2 sd+      mmq = meanMoonQuantities md ut+      MQ ml _ _ = correctedMoonQuantities lambdaS ms mmq+      d = toRadians $ ml - lambdaS+      f = 0.5 * (1 - cos d)+  in f++++-- | Calculate the Moon's position-angle of the bright limb.+-- It takes the Moon's coordinates and the Sun's coordinates.+-- Position-angle is the angle of the midpoint of the illuminated limb+-- measured eastwards from the north point of the disk.+moonBrightLimbPositionAngle :: EquatorialCoordinates1 -> EquatorialCoordinates1 -> DecimalDegrees+moonBrightLimbPositionAngle = planetBrightLimbPositionAngle+++-- | The Moon's quantities+-- Used to store intermidiate results+data MoonQuantities = MQ {+  mqLongitude :: DecimalDegrees        -- ^ the Moon's longitude+  , mqAnomaly :: DecimalDegrees        -- ^ the Moon's anomaly+  , mqAscendingNode :: DecimalDegrees  -- ^ the Moon's ascending node's longitude+  }+++-- | Calculates the Moon's mean quantities on the given date.+-- It takes the Moon's orbita details and julian date+meanMoonQuantities :: MoonDetails -> JulianDate -> MoonQuantities+meanMoonQuantities md ut =+  let d = DD $ numberOfDays (mdEpoch md) ut+      lm = reduceDegrees $ (mdL md) + 13.1763966*d  -- Moon's mean longitude+      mm = reduceDegrees $ lm - 0.1114041*d - (mdP md)  -- Moon's mean anomaly+      nm = reduceDegrees $ (mdN md) - 0.0529539*d  -- ascending node's mean longitude+  in MQ lm mm nm+++-- | Calculates correction for the equation of the centre+-- It takes the Moon's corrected anomaly in radians+centreEquation :: Double -> DecimalDegrees+centreEquation mm = DD $ 6.2886 * (sin mm)+++-- | Calculates the Moon's corrected longitude, anomaly and asceding node's longitude+-- It takes the Sun's longitude, the Sun's mean anomaly and the Moon's mean quantities+correctedMoonQuantities :: DecimalDegrees -> DecimalDegrees -> MoonQuantities -> MoonQuantities+correctedMoonQuantities lambdaS ms (MQ lm mm nm) =+  let ms' = toRadians ms+      c = lm - lambdaS+      ev = DD $ 1.2739 * (sin $ toRadians $ 2*c - mm)  -- correction for evection+      ae = DD $ 0.1858 * (sin ms')  -- correction for annual equation+      a3 = DD $ 0.37 * (sin ms')  -- third correction+      mm' = mm + (ev - ae - a3) -- Moon's corrected anomaly+      mm'' = toRadians mm'+      ec = centreEquation mm''  -- correction for the equation of the centre+      a4 = DD $ 0.214 * (sin $ 2*mm'') -- fourth correction term+      lm' = lm + (ev + ec -ae + a4) -- Moon's corrected longitude+      v = DD $ 0.6583 * (sin $ toRadians $ 2*(lm' - lambdaS))-- correction for variation+      lm'' = lm' + v -- Moon's true orbital longitude+      nm' = nm - (DD $ 0.16 * (sin ms')) -- ascending node's corrected longitude+  in MQ lm'' mm' nm'
+ src/Data/Astro/Moon/MoonDetails.hs view
@@ -0,0 +1,47 @@+{-|+Module: Data.Astro.Moon.MoonDetails+Description: Planet Details+Copyright: Alexander Ignatyev, 2016++Moon Details.+-}++module Data.Astro.Moon.MoonDetails+(+  MoonDetails(..)+  , MoonDistanceUnits(..)+  , j2010MoonDetails+  , mduToKm+)++where++import Data.Astro.Types (DecimalDegrees)+import Data.Astro.Time.Epoch (j2010)+import Data.Astro.Time.JulianDate (JulianDate(..))+++-- | Details of the Moon's orbit at the epoch+data MoonDetails = MoonDetails {+  mdEpoch :: JulianDate     -- ^ the epoch+  , mdL :: DecimalDegrees   -- ^ mean longitude at the epoch+  , mdP :: DecimalDegrees   -- ^ mean longitude of the perigee at the epoch+  , mdN :: DecimalDegrees   -- ^ mean longitude of the node at the epoch+  , mdI :: DecimalDegrees   -- ^ inclination of the orbit+  , mdE :: Double           -- ^ eccentricity of the orbit+  , mdA :: Double           -- ^ semi-major axis of the orbit+  , mdBigTheta :: DecimalDegrees  -- ^ angular diameter at the distance `mdA` from the Earth+  , mdPi :: DecimalDegrees        -- ^ parallax at distance `mdA` from the Earth+  } deriving (Show)+++-- | Moon distance units, 1 MDU = semi-major axis of the Moon's orbit+newtype MoonDistanceUnits = MDU Double deriving (Show)+++j2010MoonDetails = MoonDetails j2010 91.929336 130.143076 291.682547 5.145396 0.0549 384401 0.5181 0.9507+++-- | Convert MoonDistanceUnits to km+mduToKm :: MoonDistanceUnits -> Double+mduToKm (MDU p) = p * (mdA j2010MoonDetails)
+ src/Data/Astro/Planet.hs view
@@ -0,0 +1,100 @@+{-|+Module: Data.Astro.Planet+Description: Planet calculations+Copyright: Alexander Ignatyev, 2016++Planet calculations.++= Example++=== /Initialisation/++@+import Data.Astro.Time.JulianDate+import Data.Astro.Coordinate+import Data.Astro.Types+import Data.Astro.Effects+import Data.Astro.CelestialObject.RiseSet+import Data.Astro.Planet++ro :: GeographicCoordinates+ro = GeoC (fromDMS 51 28 40) (-(fromDMS 0 0 5))++dt :: LocalCivilTime+dt = lctFromYMDHMS (DH 1) 2017 6 25 10 29 0++today :: LocalCivilDate+today = lcdFromYMD (DH 1) 2017 6 25++jupiterDetails :: PlanetDetails+jupiterDetails = j2010PlanetDetails Jupiter++earthDetails :: PlanetDetails+earthDetails = j2010PlanetDetails Earth++jupiterPosition :: JulianDate -> EquatorialCoordinates1+jupiterPosition = planetPosition planetTrueAnomaly1 jupiterDetails earthDetails+@++=== /Calcaulate Coordinates/+@+jupiterEC1 :: EquatorialCoordinates1+jupiterEC1 = jupiterPosition (lctUniversalTime dt)+-- EC1 {e1Declination = DD (-4.104626810672402), e1RightAscension = DH 12.863365504382228}++jupiterHC :: HorizonCoordinates+jupiterHC = ec1ToHC ro (lctUniversalTime dt) jupiterEC1+-- HC {hAltitude = DD (-30.67914598469227), hAzimuth = DD 52.29376845044007}+@++=== /Calculate Distance/+@+jupiterDistance :: AstronomicalUnits+jupiterDistance = planetDistance1 jupiterDetails earthDetails (lctUniversalTime dt)+-- AU 5.193435872521039+@++=== /Calculate Angular Size/+@+jupiterAngularSize :: DecimalDegrees+jupiterAngularSize = planetAngularDiameter jupiterDetails jupiterDistance+-- DD 1.052289877865987e-2++toDMS jupiterAngularSize+-- (0,0,37.88243560317554)+@++=== /Calculate Rise and Set/++@+verticalShift :: DecimalDegrees+verticalShift = refract (DD 0) 12 1012+-- DD 0.5660098245614035++jupiterRiseSet :: RiseSetMB+jupiterRiseSet = riseAndSet2 0.000001 jupiterPosition ro verticalShift today+-- RiseSet+--    (Just (2017-06-25 13:53:27.3109 +1.0,DD 95.88943953535569))+--    (Just (2017-06-25 01:21:23.5835 +1.0,DD 264.1289033612776))+@+-}++module Data.Astro.Planet+(+  Details.Planet(..)+  , Details.PlanetDetails(..)+  , Details.j2010PlanetDetails+  , Mechanics.planetTrueAnomaly1+  , Mechanics.planetPosition+  , Mechanics.planetPosition1+  , Mechanics.planetDistance1+  , Mechanics.planetAngularDiameter+  , Mechanics.planetPhase1+  , Mechanics.planetBrightLimbPositionAngle+)++where+++import qualified Data.Astro.Planet.PlanetDetails as Details+import qualified Data.Astro.Planet.PlanetMechanics as Mechanics
+ src/Data/Astro/Planet/PlanetDetails.hs view
@@ -0,0 +1,71 @@+{-|+Module: Data.Astro.Planet.PlanetDetails+Description: Planet Details+Copyright: Alexander Ignatyev, 2016++Planet Details.+-}++module Data.Astro.Planet.PlanetDetails+(+  Planet(..)+  , PlanetDetails(..)+  , j2010PlanetDetails+  , isInnerPlanet+)++where++import Data.Astro.Types (DecimalDegrees(..), AstronomicalUnits, fromDMS)+import Data.Astro.Time.JulianDate (JulianDate)+import Data.Astro.Time.Epoch (j2010)+++-- | Planets of the Solar System+data Planet = Mercury+             | Venus+             | Earth +             | Mars+             | Jupiter+             | Saturn+             | Uranus+             | Neptune+               deriving (Show, Eq)+++-- | Details of the planetary orbit at the epoch+data PlanetDetails = PlanetDetails {+  pdPlanet :: Planet+  , pdEpoch :: JulianDate+  , pdTp :: Double               -- ^ Orbital period in tropical years+  , pdEpsilon :: DecimalDegrees  -- ^ Longitude at the Epoch+  , pdOmegaBar :: DecimalDegrees -- ^ Longitude of the perihelion+  , pdE :: Double                -- ^ Eccentricity of the orbit+  , pdAlpha :: AstronomicalUnits -- ^ Semi-major axis of the orbit in AU+  , pdI :: DecimalDegrees        -- ^ Orbital inclination+  , pdBigOmega :: DecimalDegrees -- ^ Longitude of the ascending node+  , pdBigTheta :: DecimalDegrees -- ^ Angular diameter at 1 AU+  } deriving (Show, Eq)+++-- | Return True if the planet is inner (its orbit lies inside the Earth's orbit)+isInnerPlanet :: PlanetDetails -> Bool+isInnerPlanet pd+  | pdPlanet pd == Mercury = True+  | pdPlanet pd == Venus = True+  | otherwise = False+++-- | PlanetDetails at the reference Epoch J2010.0+j2010PlanetDetails :: Planet -> PlanetDetails+j2010PlanetDetails Mercury = PlanetDetails Mercury j2010 0.24085    75.5671    77.612     0.205627 0.387098 7.0051   48.449    (arcsecs 6.74)+j2010PlanetDetails Venus   = PlanetDetails Venus   j2010 0.615207   272.30044  131.54     0.006812 0.723329 3.3947   76.769    (arcsecs 16.92)+j2010PlanetDetails Earth   = PlanetDetails Earth   j2010 0.999996   99.556772  103.2055   0.016671 0.999985 0        0         (arcsecs 0)+j2010PlanetDetails Mars    = PlanetDetails Mars    j2010 1.880765   109.09646  336.217    0.093348 1.523689 1.8497   49.632    (arcsecs 9.36)+j2010PlanetDetails Jupiter = PlanetDetails Jupiter j2010 11.857911  337.917132 14.6633    0.048907 5.20278  1.3035   100.595   (arcsecs 196.74)+j2010PlanetDetails Saturn  = PlanetDetails Saturn  j2010 29.310579  172.398316 89.567     0.053853 9.51134  2.4873   113.752   (arcsecs 165.6)+j2010PlanetDetails Uranus  = PlanetDetails Uranus  j2010 84.039492  271.063148 172.884833 0.046321 19.21814 0.773059 73.926961 (arcsecs 65.8)+j2010PlanetDetails Neptune = PlanetDetails Neptune j2010 165.845392 326.895127 23.07      0.010483 30.1985  1.7673   131.879   (arcsecs 62.2)++-- | arcseconds to DecimalHours+arcsecs = fromDMS 0 0
+ src/Data/Astro/Planet/PlanetMechanics.hs view
@@ -0,0 +1,286 @@+{-|+Module: Data.Astro.Planet.PlanetMechanics+Description: Planet mechanics+Copyright: Alexander Ignatyev, 2016++Planet mechanics.+-}++module Data.Astro.Planet.PlanetMechanics+(+  planetMeanAnomaly+  , planetTrueAnomaly1+  , planetHeliocentricRadiusVector+  , planetHeliocentricLongitude+  , planetHeliocentricLatitude+  , planetProjectedRadiusVector+  , planetProjectedLongitude+  , planetEclipticLongitude+  , planetEclipticLatitude+  , planetPosition+  , planetPosition1+  , planetDistance1+  , planetAngularDiameter+  , planetPhase1+  , planetPertubations+  , planetBrightLimbPositionAngle+)++where++import qualified Data.Astro.Utils as U+import Data.Astro.Types (DecimalDegrees(..), AstronomicalUnits(..), toRadians, fromRadians, fromDecimalHours)+import Data.Astro.Time.Epoch (j1900)+import Data.Astro.Time.JulianDate (JulianDate, numberOfDays, numberOfCenturies)+import Data.Astro.Coordinate (EquatorialCoordinates1(..), EclipticCoordinates(..), eclipticToEquatorial)+import Data.Astro.Planet.PlanetDetails (Planet(..), PlanetDetails(..), isInnerPlanet)+import Data.Astro.Sun.SunInternals (solveKeplerEquation)++{-+1. Calculate the planet position on its own orbital plane+2. Convert the planet's position to planetHeliocentric coordinates.+3. Convert from planetHeliocentric coordinates to ecliptic coordinates.+-}+++-- | reduce DecimalDegrees to the range [0, 360)+reduceDegrees :: DecimalDegrees -> DecimalDegrees+reduceDegrees = U.reduceToZeroRange 360+++-- | Calculate the planet mean anomaly.+planetMeanAnomaly pd jd =+  let d =  numberOfDays (pdEpoch pd) jd+      n = reduceDegrees $ DD $ (360/U.tropicalYearLen) * (d/(pdTp pd))+  in reduceDegrees $ n + (pdEpsilon pd) - (pdOmegaBar pd)+++-- | Calculate the planet true anomaly using approximate method+planetTrueAnomaly1 pd jd =+  let meanAnomaly = toRadians $ planetMeanAnomaly pd jd+      e = pdE pd+  in reduceDegrees $ fromRadians $ meanAnomaly + 2*e*(sin meanAnomaly)+++-- | Calculate Heliocentric Longitude.+-- It takes Planet Details and true anomaly.+planetHeliocentricLongitude :: PlanetDetails -> DecimalDegrees -> DecimalDegrees+planetHeliocentricLongitude pd trueAnomaly = reduceDegrees $ (pdOmegaBar pd) + trueAnomaly+++-- | Calculate Heliocentric Latitude.+-- It takes Planet Details and heliocentric longitude.+planetHeliocentricLatitude :: PlanetDetails -> DecimalDegrees -> DecimalDegrees+planetHeliocentricLatitude pd hcl =+  let l' = toRadians hcl+      i' = toRadians $ pdI pd+      bigOmega' = toRadians $ pdBigOmega pd+  in fromRadians $ asin $ (sin $ l' - bigOmega')*(sin i')+++-- | Calculate Heliocentric Radius Vector.+-- It takes Planet Details and true anomaly.+planetHeliocentricRadiusVector :: PlanetDetails -> DecimalDegrees -> AstronomicalUnits+planetHeliocentricRadiusVector pd trueAnomaly =+  let nu = toRadians trueAnomaly+      AU alpha = pdAlpha pd+      e = pdE pd+  in AU $ alpha*(1 - e*e)/(1+e*(cos nu))+++-- | Calculate Heliocentric Longitude projected to the ecliptic.+-- It takes Planet Details and Heliocentric Longitude+planetProjectedLongitude :: PlanetDetails -> DecimalDegrees -> DecimalDegrees+planetProjectedLongitude pd hcl =+  let hcl' = toRadians hcl+      bigOmega = pdBigOmega pd+      bigOmega' = toRadians $ bigOmega+      i' = toRadians $ pdI pd+      y = (sin $ hcl'-bigOmega')*(cos i')+      x = (cos $ hcl'-bigOmega')+      n = fromRadians $ atan2 y x+  in n + bigOmega+++-- | Calculate Heliocentric Radius Vector projected to the ecliptic.+-- It takes Planet Details, planetHeliocentric latitude and Radius Vector+planetProjectedRadiusVector :: PlanetDetails -> DecimalDegrees -> AstronomicalUnits -> AstronomicalUnits+planetProjectedRadiusVector pd psi (AU hcr) = AU $ hcr*cos(toRadians psi)+++-- | Calculate ecliptic longitude for outer planets.+-- It takes planet projected longitude, planet projected radius vector+-- the Earth's longitude and radius vector.+outerPlanetEclipticLongitude :: DecimalDegrees -> AstronomicalUnits -> DecimalDegrees -> AstronomicalUnits -> DecimalDegrees+outerPlanetEclipticLongitude lp (AU rp) le (AU re) =+  let lp' = toRadians lp+      le' = toRadians le+      x = atan $ re * (sin $ lp'-le')/(rp - re*(cos $ lp'-le'))+  in reduceDegrees $ (fromRadians x) + lp+++-- | Calculate ecliptic longitude for inner planets.+-- It takes planet projected longitude, planet projected radius vector+-- the Earth's longitude and radius vector.+innerPlanetEclipticLongitude :: DecimalDegrees -> AstronomicalUnits -> DecimalDegrees -> AstronomicalUnits -> DecimalDegrees+innerPlanetEclipticLongitude lp (AU rp) le (AU re) =+  let lp' = toRadians lp+      le' = toRadians le+      x = atan $ rp * (sin $ le'-lp')/(re - rp*(cos $ le'-lp'))+  in reduceDegrees $ (fromRadians x) + le + 180+++-- | Calculate Ecliptic Longitude.+-- It takes planet projected longitude, planet projected radius vector+-- the Earth's longitude and radius vector.+planetEclipticLongitude :: PlanetDetails -> DecimalDegrees -> AstronomicalUnits -> DecimalDegrees -> AstronomicalUnits -> DecimalDegrees+planetEclipticLongitude pd+  | isInnerPlanet pd = innerPlanetEclipticLongitude+  | otherwise = outerPlanetEclipticLongitude+++-- | Calculate ecliptic Latitude.+-- It takes the planet's: heliocentric latitude, projected heliocentric longutide,+-- projected heliocentric longitude;+-- the Earth's: heliocentric longitede and heliocentric radius vector.+-- Also it takes the planet's ecliptic longitude.+planetEclipticLatitude :: DecimalDegrees+                          -> DecimalDegrees+                          -> AstronomicalUnits+                          -> DecimalDegrees+                          -> AstronomicalUnits+                          -> DecimalDegrees+                          -> DecimalDegrees+planetEclipticLatitude psi lp (AU rp) le (AU re) lambda =+  let psi' = toRadians psi+      lp' = toRadians lp+      le' = toRadians le+      lambda' = toRadians lambda+      y = rp*(tan psi')*(sin $ lambda' - lp')+      x = re * (sin $ lp' -le')+  in fromRadians $ atan (y/x)+++-- | Calculate the planet's postion at the given date.+-- It takes a function to calculate true anomaly,+-- planet details of the planet, planet details of the Earth+-- and JulianDate.+planetPosition :: (PlanetDetails -> JulianDate -> DecimalDegrees)+                  -> PlanetDetails -> PlanetDetails -> JulianDate+                  -> EquatorialCoordinates1+planetPosition trueAnomaly pd ed jd =+      -- planet+  let nup = trueAnomaly pd jd+      lp = planetHeliocentricLongitude pd nup+      rp = planetHeliocentricRadiusVector pd nup+      psi = planetHeliocentricLatitude pd lp+      lp' = planetProjectedLongitude pd lp+      rp' = planetProjectedRadiusVector pd psi rp+      -- earth+      nue = trueAnomaly ed jd+      le = planetHeliocentricLongitude ed nue+      re = planetHeliocentricRadiusVector ed nue+      -- position+      lambda = planetEclipticLongitude pd lp' rp' le re+      beta = planetEclipticLatitude psi lp' rp' le re lambda+      ec = eclipticToEquatorial (EcC beta lambda) jd+    in ec+++-- | Calculates the distance betweeth the planet and the Earth at the given date.+-- It takes the planet's detail, the Earth's details and the julian date.+planetDistance1 :: PlanetDetails -> PlanetDetails -> JulianDate -> AstronomicalUnits+planetDistance1 pd ed jd =+  let nup = planetTrueAnomaly1 pd jd+      lp = planetHeliocentricLongitude pd nup+      AU rp = planetHeliocentricRadiusVector pd nup+      psi = planetHeliocentricLatitude pd lp+      -- earth+      nue = planetTrueAnomaly1 ed jd+      le = planetHeliocentricLongitude ed nue+      AU re = planetHeliocentricRadiusVector ed nue+      -- distance+      ro = sqrt $ re*re + rp*rp - 2*re*rp*(cos . toRadians $ lp - le)*(cos $ toRadians psi)+    in AU ro+++-- | Calculates the planet's angular diameter for the given distance.+planetAngularDiameter :: PlanetDetails -> AstronomicalUnits -> DecimalDegrees+planetAngularDiameter pd (AU ro) = (pdBigTheta pd)/(DD ro)+++-- | Calculate the planet's phase at the given phase.+-- Phase is a fraction of the visible disc that is illuminated.+-- It takes the planet's details, the Earth's details and the julian date.+-- Returns fraction values from 0 to 1.+planetPhase1 :: PlanetDetails -> PlanetDetails -> JulianDate -> Double+planetPhase1 pd ed jd =+      -- planet+  let nup = planetTrueAnomaly1 pd jd+      lp = planetHeliocentricLongitude pd nup+      rp = planetHeliocentricRadiusVector pd nup+      psi = planetHeliocentricLatitude pd lp+      lp' = planetProjectedLongitude pd lp+      rp' = planetProjectedRadiusVector pd psi rp+      -- earth+      nue = planetTrueAnomaly1 ed jd+      le = planetHeliocentricLongitude ed nue+      re = planetHeliocentricRadiusVector ed nue++      lambda = planetEclipticLongitude pd lp' rp' le re+      d = toRadians $ lambda - lp+    in (1+ (cos d)) * 0.5+++-- | Calculate the planet's postion at the given date using the approximate algoruthm.+-- It takes a function to calculate true anomaly,+-- planet details of the planet, planet details of the Earth+-- and JulianDate.+planetPosition1 :: PlanetDetails -> PlanetDetails -> JulianDate+                  -> EquatorialCoordinates1+planetPosition1 = planetPosition planetTrueAnomaly1+++-- | Calculates pertubations for the planet at the given julian date.+-- Returns a value that should be added to the mean longitude (planet heliocentric longitude).+planetPertubations :: Planet -> JulianDate -> DecimalDegrees+planetPertubations Jupiter jd =+  let (a, _, v, _) = pertubationsQuantities jd+      v' = toRadians v+      dl = (0.3314-0.0103*a)*(sin v') - 0.0644*a*(cos v')+  in DD dl+planetPertubations Saturn jd =+  let (a, q, v, b) = pertubationsQuantities jd+      q' = toRadians q+      v' = toRadians v+      b' = toRadians b+      dl = (0.1609*a-0.0105)*(cos v') + (0.0182*a-0.8142)*(sin v') - 0.1488*(sin b')+        - 0.0408*(sin $ 2*b') + 0.0856*(sin b')*(cos q') + 0.0813*(cos b')*(sin q')+  in DD dl+planetPertubations _ _ = 0+++-- pertrubationsQuantities :: JulianDate+pertubationsQuantities jd =+  let t = numberOfCenturies j1900 jd+      a = t*0.2 + 0.1+      p = DD $ 237.47555 + 3034.9061*t+      q = DD $ 265.91650 + 1222.1139*t+      v = 5*q - 2*p+      b = q - p+  in (a, q, v, b)+++-- | Calculate the planet's position-angle of the bright limb.+-- It takes the planet's coordinates and the Sun's coordinates.+-- Position-angle is the angle of the midpoint of the illuminated limb+-- measured eastwards from the north point of the disk.+planetBrightLimbPositionAngle :: EquatorialCoordinates1 -> EquatorialCoordinates1 -> DecimalDegrees+planetBrightLimbPositionAngle (EC1 deltaP alphaP) (EC1 deltaS alphaS) =+  let dAlpha = toRadians $ fromDecimalHours $ alphaS - alphaP+      deltaP' = toRadians deltaP+      deltaS' = toRadians deltaS+      y = (cos deltaS')*(sin dAlpha)+      x = (cos deltaP')*(sin deltaS') - (sin deltaP')*(cos deltaS')*(cos dAlpha)+      chi = reduceDegrees $ fromRadians $ atan2 y x+  in chi
+ src/Data/Astro/Star.hs view
@@ -0,0 +1,104 @@+{-|+Module: Data.Astro.Star+Description: Stars+Copyright: Alexander Ignatyev, 2017++Stars.++= Examples++== /Location/++@+import Data.Astro.Time.JulianDate+import Data.Astro.Coordinate+import Data.Astro.Types+import Data.Astro.Star+++ro :: GeographicCoordinates+ro = GeoC (fromDMS 51 28 40) (-(fromDMS 0 0 5))++dt :: LocalCivilTime+dt = lctFromYMDHMS (DH 1) 2017 6 25 10 29 0++-- Calculate location of Betelgeuse++betelgeuseEC1 :: EquatorialCoordinates1+betelgeuseEC1 = starCoordinates Betelgeuse+-- EC1 {e1Declination = DD 7.407064, e1RightAscension = DH 5.919529}++betelgeuseHC :: HorizonCoordinates+betelgeuseHC = ec1ToHC ro (lctUniversalTime dt) betelgeuseEC1+-- HC {hAltitude = DD 38.30483892505852, hAzimuth = DD 136.75755644642248}+@++== /Rise and Set/++@+import Data.Astro.Time.JulianDate+import Data.Astro.Coordinate+import Data.Astro.Types+import Data.Astro.Effects+import Data.Astro.CelestialObject.RiseSet+import Data.Astro.Star+++ro :: GeographicCoordinates+ro = GeoC (fromDMS 51 28 40) (-(fromDMS 0 0 5))++today :: LocalCivilDate+today = lcdFromYMD (DH 1) 2017 6 25++-- Calculate location of Betelgeuse++rigelEC1 :: EquatorialCoordinates1+rigelEC1 = starCoordinates Rigel++verticalShift :: DecimalDegrees+verticalShift = refract (DD 0) 12 1012+-- DD 0.5660098245614035++rigelRiseSet :: RiseSetLCT+rigelRiseSet = riseAndSetLCT ro today verticalShift rigelEC1+-- RiseSet (2017-06-25 06:38:18.4713 +1.0,DD 102.51249855335433) (2017-06-25 17:20:33.4902 +1.0,DD 257.48750144664564)+@+-}+++module Data.Astro.Star+(+  Star(..)+  , starCoordinates+)++where++import Data.Astro.Coordinate (EquatorialCoordinates1(..))+import Data.Astro.Types (fromDMS, fromHMS)+++-- | Some of the stars+data Star = Polaris+            | AlphaCrucis+            | Sirius+            | Betelgeuse+            | Rigel+            | Vega+            | Antares+            | Canopus+            | Pleiades+              deriving (Show, Eq)+++-- | Returns Equatorial Coordinates for the given star+starCoordinates :: Star -> EquatorialCoordinates1+starCoordinates Polaris = EC1 (fromDMS 89 15 51) (fromHMS 2 31 48.7)+starCoordinates AlphaCrucis = EC1 (-(fromDMS 63 5 56.73)) (fromHMS 12 26 35.9)+starCoordinates Sirius = EC1 (-(fromDMS 16 42 58.02)) (fromHMS 6 45 8.92)+starCoordinates Betelgeuse = EC1 (fromDMS 07 24 25.4304) (fromHMS 5 55 10.30536)+starCoordinates Rigel = EC1 (-(fromDMS 8 12 05.8981)) (fromHMS 5 14 32.27210)+starCoordinates Vega = EC1 (fromDMS 38 47 01.2802) (fromHMS 18 36 56.33635)+starCoordinates Antares = EC1 (-(fromDMS 26 25 55.2094)) (fromHMS 16 29 24.45970)+starCoordinates Canopus = EC1 (-(fromDMS 52 41 44.3810)) (fromHMS 6 23 57.10988)+starCoordinates Pleiades = EC1 (fromDMS 24 7 00) (fromHMS 3 47 24)
+ src/Data/Astro/Sun.hs view
@@ -0,0 +1,252 @@+{-|+Module: Data.Astro.Sun+Description: Calculation characteristics of the Sun+Copyright: Alexander Ignatyev, 2016++= Calculation characteristics of the Sun.++== /Terms/++* __perihelion__ - minimal distance from the Sun to the planet+* __aphelion__ - maximal distance from the Sun to the planet++* __perigee__ - minimal distance from the Sun to the Earth+* __apogee__ - maximal distance from the Sun to the Earth+++= Example++@+import Data.Astro.Time.JulianDate+import Data.Astro.Coordinate+import Data.Astro.Types+import Data.Astro.Sun++ro :: GeographicCoordinates+ro = GeoC (fromDMS 51 28 40) (-(fromDMS 0 0 5))++dt :: LocalCivilTime+dt = lctFromYMDHMS (DH 1) 2017 6 25 10 29 0++today :: LocalCivilDate+today = lcdFromYMD (DH 1) 2017 6 25++jd :: JulianDate+jd = lctUniversalTime dt++verticalShift :: DecimalDegrees+verticalShift = refract (DD 0) 12 1012++-- distance from the Earth to the Sun in kilometres+distance :: Double+distance = sunDistance jd+-- 1.5206375976421073e8++-- Angular Size+angularSize :: DecimalDegrees+angularSize = sunAngularSize jd+-- DD 0.5244849215333616++-- The Sun's coordinates+ec1 :: EquatorialCoordinates1+ec1 = sunPosition2 jd+-- EC1 {e1Declination = DD 23.37339098989099, e1RightAscension = DH 6.29262026252748}++hc :: HorizonCoordinates+hc = ec1ToHC ro jd ec1+-- HC {hAltitude = DD 49.312050979507404, hAzimuth = DD 118.94723825710143}+++-- Rise and Set+riseSet :: RiseSetMB+riseSet = sunRiseAndSet ro 0.833333 today+-- RiseSet+--    (Just (2017-06-25 04:44:04.3304 +1.0,DD 49.043237261724215))+--    (Just (2017-06-25 21:21:14.4565 +1.0,DD 310.91655607595595))+@+-}++module Data.Astro.Sun+(+  SunDetails(..)+  , RiseSet(..)+  , sunDetails+  , j2010SunDetails+  , sunMeanAnomaly2+  , sunEclipticLongitude1+  , sunEclipticLongitude2+  , sunPosition1+  , sunPosition2+  , sunDistance+  , sunAngularSize+  , sunRiseAndSet+  , equationOfTime+  , solarElongation+)++where++import qualified Data.Astro.Utils as U+import Data.Astro.Types (DecimalDegrees(..), DecimalHours(..)+                        , toDecimalHours, fromDecimalHours+                        , toRadians, fromRadians+                        , GeographicCoordinates(..) )+import Data.Astro.Time.JulianDate (JulianDate(..), LocalCivilTime(..), LocalCivilDate(..), numberOfDays, numberOfCenturies, splitToDayAndTime, addHours)+import Data.Astro.Time.Sidereal (gstToUT, dhToGST)+import Data.Astro.Time.Epoch (j1900, j2010)+import Data.Astro.Coordinate (EquatorialCoordinates1(..), EclipticCoordinates(..), eclipticToEquatorial)+import Data.Astro.Effects.Nutation (nutationLongitude)+import Data.Astro.CelestialObject.RiseSet (RiseSet(..), RiseSetMB, RSInfo(..), riseAndSet2)+import Data.Astro.Sun.SunInternals (solveKeplerEquation)+++-- | Details of the Sun's apparent orbit at the given epoch+data SunDetails = SunDetails {+  sdEpoch :: JulianDate             -- ^ Epoch+  , sdEpsilon :: DecimalDegrees     -- ^ Ecliptic longitude at the Epoch+  , sdOmega :: DecimalDegrees       -- ^ Ecliptic longitude of perigee at the Epoch+  , sdE :: Double                   -- ^ Eccentricity of the orbit at the Epoch+  } deriving (Show)++-- | SunDetails at the Sun's reference Epoch J2010.0+j2010SunDetails :: SunDetails+j2010SunDetails = SunDetails j2010 (DD 279.557208) (DD 283.112438) 0.016705+++-- | Semi-major axis+r0 :: Double+r0 = 1.495985e8+++-- | Angular diameter at r = r0+theta0 :: DecimalDegrees+theta0 = DD 0.533128+++-- | Reduce the value to the range [0, 360)+reduceTo360 :: Double -> Double+reduceTo360 = U.reduceToZeroRange 360+++-- | Reduce the value to the range [0, 360)+reduceDegrees :: DecimalDegrees -> DecimalDegrees+reduceDegrees = U.reduceToZeroRange 360+++-- | Calculate SunDetails for the given JulianDate.+sunDetails :: JulianDate -> SunDetails+sunDetails jd =+  let t = numberOfCenturies j1900 jd+      epsilon = reduceTo360 $ 279.6966778 + 36000.76892*t + 0.0003025*t*t+      omega = reduceTo360 $ 281.2208444 + 1.719175*t + 0.000452778*t*t+      e = 0.01675104 - 0.0000418*t - 0.000000126*t*t+  in SunDetails jd (DD epsilon) (DD omega) e+++-- | Calculate the ecliptic longitude of the Sun with the given SunDetails at the given JulianDate+sunEclipticLongitude1 :: SunDetails -> JulianDate -> DecimalDegrees+sunEclipticLongitude1 sd@(SunDetails epoch (DD eps) (DD omega) e) jd =+  let d = numberOfDays epoch jd+      n = reduceTo360 $ (360/U.tropicalYearLen) * d+      meanAnomaly = reduceTo360 $ n + eps - omega+      ec = (360/pi)*e*(sin $ U.toRadians meanAnomaly)+      DD nutation = nutationLongitude jd+  in DD $ reduceTo360 $ n + ec + eps + nutation+++-- | Calculate Equatorial Coordinates of the Sun with the given SunDetails at the given JulianDate.+-- It is recommended to use 'j2010SunDetails' as a first parameter.+sunPosition1 :: SunDetails -> JulianDate -> EquatorialCoordinates1+sunPosition1 sd jd =+  let lambda = sunEclipticLongitude1 sd jd+      beta = DD 0+  in eclipticToEquatorial (EcC beta lambda) jd+++-- | Calculate mean anomaly using the second 'more accurate' method+sunMeanAnomaly2 :: SunDetails -> DecimalDegrees+sunMeanAnomaly2 sd = reduceDegrees $ (sdEpsilon sd) - (sdOmega sd)+++-- | Calculate true anomaly using the second 'more accurate' method+trueAnomaly2 :: SunDetails -> DecimalDegrees+trueAnomaly2 sd =+  let m = toRadians $ sunMeanAnomaly2 sd+      e = sdE sd+      bigE = solveKeplerEquation e m 0.000000001+      tanHalfNu = sqrt((1+e)/(1-e)) * tan (0.5 * bigE)+      nu = reduceTo360 $ U.fromRadians $ 2 * (atan tanHalfNu)+  in DD nu+++-- | Calculate the ecliptic longitude of the Sun+sunEclipticLongitude2 :: SunDetails -> DecimalDegrees+sunEclipticLongitude2 sd =+  let DD omega = sdOmega sd+      DD nu = trueAnomaly2 sd+      DD nutation = nutationLongitude $ sdEpoch sd+  in DD $ reduceTo360 $ nu + omega + nutation+++-- | More accurate method to calculate position of the Sun+sunPosition2 :: JulianDate -> EquatorialCoordinates1+sunPosition2 jd =+  let sd = sunDetails jd+      lambda = sunEclipticLongitude2 sd+      beta = DD 0+  in eclipticToEquatorial (EcC beta lambda) jd+++-- Distance and Angular Size helper function+dasf sd =+  let e = sdE sd+      nu = toRadians $ trueAnomaly2 sd+  in (1 + e*(cos nu)) / (1 - e*e)+++-- | Calculate Sun-Earth distance.+sunDistance :: JulianDate -> Double+sunDistance jd = r0 / (dasf $ sunDetails jd)+++-- | Calculate the Sun's angular size (i.e. its angular diameter).+sunAngularSize :: JulianDate -> DecimalDegrees+sunAngularSize jd = theta0 * (DD $ dasf $ sunDetails jd)+++-- | Calculatesthe Sun's rise and set+-- It takes coordinates of the observer,+-- local civil date,+-- vertical shift (good value is 0.833333).+-- It returns Nothing if fails to calculate rise and/or set.+-- It should be accurate to within a minute of time.+sunRiseAndSet :: GeographicCoordinates+                 -> DecimalDegrees+                 -> LocalCivilDate+                 -> RiseSetMB+sunRiseAndSet = riseAndSet2 0.000001 (sunPosition1 j2010SunDetails)+++-- | Calculates discrepancy between the mean solar time and real solar time+-- at the given date.+equationOfTime :: JulianDate -> DecimalHours+equationOfTime jd =+  let (day, _) = splitToDayAndTime jd+      midday = addHours (DH 12) day  -- mean solar time+      EC1 _ ra = sunPosition1 j2010SunDetails midday+      ut = gstToUT day $ dhToGST ra+      JD time = midday - ut+  in DH $ time*24+++-- | Calculates the angle between the lines of sight to the Sun and to a celestial object+-- specified by the given coordinates at the given Universal Time.+solarElongation :: EquatorialCoordinates1 -> JulianDate -> DecimalDegrees+solarElongation (EC1 deltaP alphaP) jd =+  let (EC1 deltaS alphaS) = sunPosition1 j2010SunDetails jd+      deltaP' = toRadians deltaP+      alphaP' = toRadians $ fromDecimalHours alphaP+      deltaS' = toRadians deltaS+      alphaS' = toRadians $ fromDecimalHours alphaS+      eps = acos $ (sin deltaP')*(sin deltaS') + (cos $ alphaP' - alphaS')*(cos deltaP')*(cos deltaS')+  in fromRadians eps
+ src/Data/Astro/Sun/SunInternals.hs view
@@ -0,0 +1,28 @@+{-|+Module: Data.Astro.Sun.SunInternals+Description: Internal functions of Sun module.+Copyright: Alexander Ignatyev, 2016++Internal functions of Sun module. Exposed only for Unit Tests+-}++module Data.Astro.Sun.SunInternals+(+  solveKeplerEquation+)++where+++-- | Solve Kepler's Equation: E - e * (sin E) = M+-- It takes eccentricity,+-- mean anomaly in radians equals epsilon - omega (see 'SunDetails').+-- It returns E in radians.+solveKeplerEquation :: Double -> Double -> Double -> Double+solveKeplerEquation e m eps = iter m+  where iter x =+          let delta = x - e*(sin x) - m+              dx = delta / (1 - e*(cos x))+          in if abs delta < eps+             then x+             else iter (x-dx)
+ src/Data/Astro/Time.hs view
@@ -0,0 +1,60 @@+{-|+Module: Data.Astro.Time+Description: Time+Copyright: Alexander Ignatyev, 2016++Root Time module+-}+++module Data.Astro.Time+(+  utToLST+  , lctToLST+  , lstToLCT+)++where++import Data.Astro.Types (DecimalDegrees)+import Data.Astro.Time.JulianDate (JulianDate(..), LocalCivilTime(..), LocalCivilDate(..), splitToDayAndTime, addHours)+import Data.Astro.Time.Sidereal (LocalSiderealTime, utToGST, gstToUT, gstToLST, lstToGST, lstToGSTwDC)+++-- | Universal Time to Local Sidereal Time.+-- It takes longitude in decimal degrees and local civil time+utToLST :: DecimalDegrees -> JulianDate -> LocalSiderealTime+utToLST longitude ut = gstToLST longitude $ utToGST ut+++-- | Local Civil Time to Local Sidereal Time.+-- It takes longitude in decimal degrees and local civil time+lctToLST :: DecimalDegrees -> LocalCivilTime -> LocalSiderealTime+lctToLST longitude lct = utToLST longitude $ lctUniversalTime lct+++-- | Local Sidereal Time to Local Civil Time.+-- It takes longitude in decimal degrees, local civil date and local sidereal time+lstToLCT :: DecimalDegrees -> LocalCivilDate -> LocalSiderealTime -> LocalCivilTime+lstToLCT longitude lcd lst =+  let gst = lstToGST longitude lst+      ut = gstToUT (lcdDate lcd) gst+      lct = LCT (lcdTimeZone lcd) ut+  in if sameDay lcd lct+     then lct -- lstToLCTwDC longitude timeZone jd lst+     else lstToLCTwDC longitude lcd lst+++lstToLCTwDC :: DecimalDegrees -> LocalCivilDate -> LocalSiderealTime -> LocalCivilTime+lstToLCTwDC longitude lcd lst =+  let gst = lstToGSTwDC longitude lst+      ut = gstToUT (lcdDate lcd) gst+      lct = LCT (lcdTimeZone lcd) ut+  in lct+++-- | Returns True if both JulianDates hve the same day+sameDay :: LocalCivilDate -> LocalCivilTime -> Bool+sameDay (LCD _ (JD d1)) (LCT tz jd2) =+  let (JD d2, _) = splitToDayAndTime $ addHours tz jd2+  in abs (d1 - d2) < 0.000001
+ src/Data/Astro/Time/Conv.hs view
@@ -0,0 +1,74 @@+{-|+Module: Data.Astro.Time.Conv+Description: Julian Date+Copyright: Alexander Ignatyev, 2017+++Conversion functions between datetime types defined in Data.Time and Data.Astro.Time modules.+-}+module Data.Astro.Time.Conv+(+  zonedTimeToLCT+  , zonedTimeToLCD+  , lctToZonedTime+)++where+++import Data.Time.LocalTime (ZonedTime(..), LocalTime(..)+                           , TimeOfDay(..), TimeZone(..)+                           , minutesToTimeZone)+import Data.Time.Calendar (toGregorian, fromGregorian)++import Data.Astro.Types(DecimalHours(..))+import Data.Astro.Utils (fromFixed)+import Data.Astro.Time.JulianDate (JulianDate(..)+                                  , LocalCivilTime(..)+                                  , LocalCivilDate(..)+                                  , fromYMDHMS, toYMDHMS+                                  , lctFromYMDHMS, lcdFromYMD+                                  , lctToYMDHMS)+++-----------------------------------------------------------+-- Data.Time types -> Data.Astro types+timeZoneToDH :: TimeZone -> DecimalHours+timeZoneToDH  tz = DH hours+  where toMinutes = fromIntegral . timeZoneMinutes+        hours = (toMinutes tz) / 60.0+++-- | Convert ZonedTime to LocalCivilTime+zonedTimeToLCT :: ZonedTime -> LocalCivilTime+zonedTimeToLCT zonedTime = lctFromYMDHMS tz y m d hours mins (fromFixed secs)+  where tz = timeZoneToDH (zonedTimeZone zonedTime)+        lt = zonedTimeToLocalTime zonedTime+        (y, m, d) = toGregorian (localDay lt)+        TimeOfDay hours mins secs = localTimeOfDay lt+++-- | Convert ZonedTime to LocalCivilDate+zonedTimeToLCD :: ZonedTime -> LocalCivilDate+zonedTimeToLCD zonedTime = lcdFromYMD tz y m d+  where tz = timeZoneToDH (zonedTimeZone zonedTime)+        lt = zonedTimeToLocalTime zonedTime+        (y, m, d) = toGregorian (localDay lt)+++-----------------------------------------------------------+-- Data.Astro Types -> Data.Time types++dhToTimeZone :: DecimalHours -> TimeZone+dhToTimeZone (DH hours) = minutesToTimeZone minutes+  where minutes = round (60*hours)+++-- | Convert LocalCivilTime to ZonedTime+lctToZonedTime :: LocalCivilTime -> ZonedTime+lctToZonedTime lct = ZonedTime { zonedTimeToLocalTime = lt, zonedTimeZone = tz }+  where tz = dhToTimeZone $ lctTimeZone lct+        (y, m, d, hours, mins, secs) = lctToYMDHMS lct+        day = fromGregorian y m d+        time = TimeOfDay hours mins (realToFrac secs)+        lt = LocalTime { localDay = day, localTimeOfDay = time }
+ src/Data/Astro/Time/Epoch.hs view
@@ -0,0 +1,51 @@+{-|+Module: Data.Astro.Time.Epoch+Description: Astronomical Epochs+Copyright: Alexander Ignatyev, 2016++Definitions of well-known astronomical epochs.+-}+module Data.Astro.Time.Epoch+(+    -- * Epochs+    -- ** Besselian Epochs+  b1900+  , b1950+    -- ** New Epochs+  , j1900+  , j2000+  , j2050+    -- ** Well-known epochs+  , j2010+)++where+++import Data.Astro.Time.JulianDate (JulianDate(..))++-- | Epoch B1900.0, 1900 January 0.8135+b1900 :: JulianDate+b1900 = JD 2415020.3135++-- | Epoch B1950.0, January 0.9235+b1950 :: JulianDate+b1950 = JD 2433282.4235+++-- | Epoch J1900.0 1900 January 0.5+j1900 :: JulianDate+j1900 = JD 2415020.0++-- | Epoch J2000.0, 12h on 1 January 2000+j2000 :: JulianDate+j2000 = JD 2451545.0++-- | Epoch J2050.0, 12h on 1 January 2000+j2050 :: JulianDate+j2050 = JD 2469807.50+++-- | The Sun's and planets reference Epoch J2010.0 (2010 January 0.0)+j2010 :: JulianDate+j2010 = JD 2455196.5
+ src/Data/Astro/Time/GregorianCalendar.hs view
@@ -0,0 +1,86 @@+{-|+Module: Data.Astro.Time.GregorianCalendar+Description: Gregorian Calendar+Copyright: Alexander Ignatyev, 2016+++Gregorian Calendar was introduced by Pope Gregory XIII.+He abolished the days 1582-10-05 to 1582-10-14 inclusive to bring back civil and tropical years back to line.+-}++module Data.Astro.Time.GregorianCalendar+(+  isLeapYear+  , dayNumber+  , easterDayInYear+  , gregorianDateAdjustment+)++where++import Data.Time.Calendar (Day(..), fromGregorian, toGregorian)++-- Date after 15 October 1582 belongs to Gregorian Calendar+-- Before this date - to Julian Calendar+isGregorianDate :: Integer -> Int -> Int -> Bool+isGregorianDate y m d = y > gyear+  || (y == gyear && m > gmonth)+  || (y == gyear && m == gmonth && d >= gday)+  where gyear = 1582+        gmonth = 10+        gday = 15+++gregorianDateAdjustment :: Integer -> Int ->Int -> Int+gregorianDateAdjustment year month day =+  if isGregorianDate year month day+  then let y = if month < 3 then year - 1 else year+           y' = fromIntegral y+           a = truncate (y' / 100)+       in 2 - a + truncate(fromIntegral a/4)+  else 0+++-- | Check Gregorian calendar leap year+isLeapYear :: Integer -> Bool+isLeapYear year =+  year `mod` 4 == 0+  && (year `mod` 100 /= 0 || year `mod` 400 == 0)+++-- | Day Number in a year+dayNumber :: Day -> Int+dayNumber date =+  (daysBeforeMonth year month) + day+  where (year, month, day) = toGregorian date+++-- | Get Easter date+-- function uses absolutely crazy Butcher's algorithm+easterDayInYear :: Int -> Day+easterDayInYear year =+  let  a = year `mod` 19+       b = year `div` 100+       c = year `mod` 100+       d = b `div` 4+       e = b `mod` 4+       f = (b+8) `div` 25+       g = (b-f+1) `div` 3+       h = (19*a+b-d-g+15) `mod` 30+       i = c `div` 4+       k = c `mod` 4+       l = (32+2*e+2*i-h-k) `mod` 7+       m = (a+11*h+22*l) `div` 451+       n' = (h+l-7*m+114)+       n = n' `div` 31+       p = n' `mod` 31+  in fromGregorian (fromIntegral year) n (p+1)+++daysBeforeMonth :: Integer -> Int -> Int+daysBeforeMonth year month =+  let a = if isLeapYear year then 62 else 63+      month' = (fromIntegral month) :: Double+  in if month > 2 then+    truncate $ ((month' + 1.0) * 30.6) - a+  else truncate $ (month' - 1.0)*a*0.5
+ src/Data/Astro/Time/JulianDate.hs view
@@ -0,0 +1,240 @@+{-|+Module: Data.Astro.Time.JulianDate+Description: Julian Date+Copyright: Alexander Ignatyev, 2016+++Julian date is the continuous count of days since noon on January 1, 4713 BC,+the beginning of the Julian Period.++= Examples++== /JulianDate/+@+import Data.Astro.Time.JulianDate++-- 2017-06-25 9:29:00 (GMT)+jd :: JulianDate+jd = fromYMDHMS 2017 6 25 9 29 0+-- JD 2457929.895138889+@++== /LocalCiviTime and LocalCivilDate/++@+import Data.Astro.Time.JulianDate+import Data.Astro.Types++-- 2017-06-25 10:29:00 +0100 (BST)+lct :: LocalCivilTime+lct = lctFromYMDHMS (DH 1) 2017 6 25 10 29 0+-- 2017-06-25 10:29:00.0000 +1.0++lctJD :: JulianDate+lctJD = lctUniversalTime lct+-- JD 2457929.895138889++lctTZ :: DecimalHours+lctTZ = lctTimeZone lct+-- DH 1.0++lcd :: LocalCivilDate+lcd = lcdFromYMD (DH 1) 2017 6 25++lcdJD :: JulianDate+lcdJD = lcdDate lcd+-- JD 2457929.5++lcdTZ :: DecimalHours+lcdTZ = lcdTimeZone lcd+-- DH 1.0+@+-}++module Data.Astro.Time.JulianDate+(+  JulianDate(..)+  , julianStartDateTime+  , LocalCivilTime(..)+  , LocalCivilDate(..)+  , TimeBaseType+  , numberOfDays+  , numberOfYears+  , numberOfCenturies+  , addHours+  , fromYMD+  , fromYMDHMS+  , toYMDHMS+  , dayOfWeek+  , splitToDayAndTime+  , lctFromYMDHMS+  , lctToYMDHMS+  , lcdFromYMD+  , printLctHs+)++where++import Text.Printf (printf)++import Data.Astro.Types(DecimalHours(..), fromHMS, toHMS)+import Data.Astro.Time.GregorianCalendar (gregorianDateAdjustment)+import Data.Astro.Utils (trunc, fraction)+++type TimeBaseType = Double++-- | A number of days since noon of 1 January 4713 BC+newtype JulianDate = JD TimeBaseType+                     deriving (Show, Eq)+++-- | Represents Local Civil Time+data LocalCivilTime = LCT {+  lctTimeZone :: DecimalHours   -- Time Zone correction+  , lctUniversalTime :: JulianDate+  } deriving (Eq)+++instance Show LocalCivilTime where+  show = printLct+++-- | Local Civil Date, used for time conversions when base date is needed+data LocalCivilDate = LCD {+  lcdTimeZone :: DecimalHours+  , lcdDate :: JulianDate+  } deriving (Eq)+++-- | Beginning of the Julian Period+julianStartDateTime = fromYMDHMS (-4712) 1 1 12 0 0+++instance Num JulianDate where+  (+) (JD d1) (JD d2) = JD (d1+d2)+  (-) (JD d1) (JD d2) = JD (d1-d2)+  (*) (JD d1) (JD d2) = JD (d1*d2)+  negate (JD d) = JD (negate d)+  abs (JD d) = JD (abs d)+  signum (JD d) = JD (signum d)+  fromInteger int = JD (fromInteger int)+++-- | Return number of days since the first argument till the second one+numberOfDays :: JulianDate -> JulianDate -> TimeBaseType+numberOfDays (JD jd1) (JD jd2) = jd2 - jd1+++-- | Return number of years since the first argument till the second one+numberOfYears :: JulianDate -> JulianDate -> TimeBaseType+numberOfYears (JD jd1) (JD jd2) = (jd2-jd1) / 365.25+++-- | Return number of centuries since the first argument till the second one+numberOfCenturies :: JulianDate -> JulianDate -> TimeBaseType+numberOfCenturies (JD jd1) (JD jd2) = (jd2-jd1) / 36525+++-- | add Decimal Hours+addHours :: DecimalHours -> JulianDate -> JulianDate+addHours (DH hours) jd = jd + (JD $ hours/24)+++-- | Create Julian Date.+-- It takes year, month [1..12], Day [1..31].+fromYMD :: Integer -> Int -> Int -> JulianDate+fromYMD year month day =+  let (y, m) = if month < 3 then (year-1, month+12) else (year, month)+      y' = fromIntegral y+      m' = fromIntegral m+      b = gregorianDateAdjustment year month day+      c = if y < 0+          then truncate (365.25*y' - 0.75)  -- 365.25 - number of solar days in a year+          else truncate (365.25*y')+      d = truncate (30.6001 * (m'+1))+      jd = fromIntegral (b + c + d + day) + 1720994.5  -- add 1720994.5 to process BC/AC border+  in JD jd+++-- | Create Julian Date.+-- It takes year, month [1..12], Day [1..31], hours, minutes, seconds.+fromYMDHMS :: Integer -> Int -> Int -> Int -> Int -> TimeBaseType -> JulianDate+fromYMDHMS year month day hs ms ss = addHours (fromHMS hs ms ss) (fromYMD year month day)+++-- | It returns year, month [1..12], Day [1..31], hours, minutes, seconds.+toYMDHMS :: JulianDate -> (Integer, Int, Int, Int, Int, TimeBaseType)+toYMDHMS (JD jd) =+  let (i, time) = fraction (jd + 0.5)+      b = if i > 2299160  -- 2299161 - first day of Georgian Calendar+          then let a = trunc $ (i-1867216.25)/36524.25+               in i + a - trunc (a*0.25) + 1+          else i+      c = b + 1524+      d = trunc $ (c-122.1)/365.25+      e = trunc $ d * 365.25+      g = trunc $ (c-e)/30.6001+      day = truncate $ c - e - trunc (30.6001*g)+      month = truncate $ if g < 13.5 then g - 1 else g - 13+      year = truncate $ if month > 2 then d-4716 else d-4715+      (h, m, s) = toHMS $ DH $ 24*time+   in (year, month, day, h, m, s)+++-- | Get Day of the Week+-- 0 is for Sunday, 1 for manday and 6 for Saturday+dayOfWeek :: JulianDate -> Int+dayOfWeek jd =+  let JD d = removeHours jd+      (_, f) = properFraction $ (d+1.5) / 7+  in round (7*f)+++-- | Extract Day and Time parts of Date+splitToDayAndTime :: JulianDate -> (JulianDate, JulianDate)+splitToDayAndTime jd@(JD n) =+  let day = JD $ 0.5 + trunc (n - 0.5)+      time = jd - day+  in (day, time)+++-- | Get Julian date corresponding to midnight+removeHours :: JulianDate -> JulianDate+removeHours jd =+  let (d, _) = splitToDayAndTime jd+  in d+++-- | Create LocalCivilTime from tize zone, local year, local month, local day, local hours, local minutes and local secunds.+lctFromYMDHMS :: DecimalHours ->Integer -> Int -> Int -> Int -> Int -> TimeBaseType -> LocalCivilTime+lctFromYMDHMS tz y m d hs ms ss =+  let jd = fromYMDHMS y m d hs ms ss+      jd' = addHours (-tz) jd+  in LCT tz jd'+++-- | Get from LocalCivilTime local year, local month, local day, local hours, local minutes and local secunds.+lctToYMDHMS :: LocalCivilTime -> (Integer, Int, Int, Int, Int, TimeBaseType)+lctToYMDHMS (LCT tz jd)= toYMDHMS (addHours tz jd)+++-- Create LocalCivilDate from time zone, local year, local month, local day+lcdFromYMD :: DecimalHours -> Integer -> Int -> Int -> LocalCivilDate+lcdFromYMD tz y m d = LCD tz (fromYMD y m d)+++-- | Print Local Civil Time in human-readable format+printLct :: LocalCivilTime -> String+printLct lct =+  printf "%d-%02d-%02d %02d:%02d:%07.4f %+03.1f" y m d hs ms ss tz+  where (y, m, d, hs, ms, ss) = lctToYMDHMS lct+        DH tz = lctTimeZone lct+++-- | Print local civil time in machine readable format+printLctHs :: LocalCivilTime -> String+printLctHs lct =+  printf "lctFromYMDHMS (%1.0f) %d %d %d %d %d %.4f" tz y m d hs ms ss+  where (y, m, d, hs, ms, ss) = lctToYMDHMS lct+        DH tz = lctTimeZone lct
+ src/Data/Astro/Time/Sidereal.hs view
@@ -0,0 +1,143 @@+{-|+Module: Data.Astro.Time.Sidereal+Description: Sidereal Time+Copyright: Alexander Ignatyev, 2016++According to the Sidereal Clock any observed star returns to the same position+in the sky every 24 hours.++Each sidereal day is shorter than the solar day, 24 hours of sidereal time+corresponding to 23:56:04.0916 of solar time.+-}++module Data.Astro.Time.Sidereal+(+  GreenwichSiderealTime+  , LocalSiderealTime+  , dhToGST+  , dhToLST+  , gstToDH+  , lstToDH+  , hmsToGST+  , hmsToLST+  , utToGST+  , gstToUT+  , gstToLST+  , lstToGST+  , lstToGSTwDC+)+where++import Data.Astro.Types (DecimalHours(..), fromHMS)+import Data.Astro.Time.JulianDate (JulianDate(..), TimeBaseType, numberOfCenturies, splitToDayAndTime)+import Data.Astro.Time.Epoch (j2000)+import Data.Astro.Utils (reduceToZeroRange)+import qualified Data.Astro.Types as C+++-- | Greenwich Sidereal Time+-- GST can be in range [-12h, 36h] carrying out a day correction+newtype GreenwichSiderealTime = GST TimeBaseType deriving (Show, Eq)+++-- | Local Sidereal Time+newtype LocalSiderealTime = LST TimeBaseType deriving (Show, Eq)+++-- | Convert Decimal Hours to Greenwich Sidereal Time+dhToGST :: DecimalHours -> GreenwichSiderealTime+dhToGST (DH t) = GST t+++-- | Convert Decimal Hours to Local Sidereal Time+dhToLST :: DecimalHours -> LocalSiderealTime+dhToLST (DH t) = LST t+++-- | Convert Greenwich Sidereal Time to Decimal Hours+gstToDH :: GreenwichSiderealTime -> DecimalHours+gstToDH (GST t) = DH t+++-- | Convert Local Sidereal Time to Decimal Hours+lstToDH :: LocalSiderealTime -> DecimalHours+lstToDH (LST t) = DH t+++-- | Comvert Hours, Minutes, Seconds to Greenwich Sidereal Time+hmsToGST :: Int -> Int -> TimeBaseType -> GreenwichSiderealTime+hmsToGST h m s = dhToGST $ fromHMS h m s+++-- | Comvert Hours, Minutes, Seconds to Local Sidereal Time+hmsToLST :: Int -> Int -> TimeBaseType -> LocalSiderealTime+hmsToLST h m s = dhToLST $ fromHMS h m s+++-- | Convert from Universal Time (UT) to Greenwich Sidereal Time (GST)+utToGST :: JulianDate -> GreenwichSiderealTime+utToGST jd =+  let (JD day, JD time) = splitToDayAndTime jd+      t = solarSiderealTimesDiff day+      time' = reduceToZeroRange 24 $ time*24/siderealDayLength + t+  in GST $ time'+++-- | Convert from Greenwich Sidereal Time (GST) to Universal Time (UT).+-- It takes GST and Greenwich Date, returns JulianDate.+-- Because the sidereal day is shorter than the solar day (see comment to the module).+-- In case of such ambiguity the early time will be returned.+-- You can easily check the ambiguity: if time is equal or less 00:03:56+-- you can get the second time by adding 23:56:04+gstToUT :: JulianDate -> GreenwichSiderealTime -> JulianDate+gstToUT jd gst =+  let (day, time) = dayTime jd gst+      t = solarSiderealTimesDiff day+      time' = (reduceToZeroRange 24 (time-t)) * siderealDayLength+  in JD $ day + time'/24+  where dayTime jd (GST gst)+          | gst < 0   = (day-1, gst+24)+          | gst >= 24 = (day+1, gst-24)+          | otherwise = (day,   gst)+            where (JD day, _) = splitToDayAndTime jd+++-- | Convert Greenwich Sidereal Time to Local Sidereal Time.+-- It takes GST and longitude in decimal degrees+gstToLST :: C.DecimalDegrees -> GreenwichSiderealTime -> LocalSiderealTime+gstToLST longitude (GST gst) =+  let C.DH dhours = C.toDecimalHours longitude+      lst = reduceToZeroRange 24 $ gst + dhours+  in LST lst+++-- | Convert Local Sidereal Time to Greenwich Sidereal Time+-- It takes LST and longitude in decimal degrees+lstToGST :: C.DecimalDegrees -> LocalSiderealTime -> GreenwichSiderealTime+lstToGST longitude (LST lst) =+  let C.DH dhours = C.toDecimalHours longitude+      gst = reduceToZeroRange 24 $ lst - dhours+  in GST gst+++-- | Convert Local Sidereal Time to Greenwich Sidereal Time with Day Correction.+-- It takes LST and longitude in decimal degrees+lstToGSTwDC :: C.DecimalDegrees -> LocalSiderealTime -> GreenwichSiderealTime+lstToGSTwDC longitude (LST lst) =+  let C.DH dhours = C.toDecimalHours longitude+      gst = lst - dhours+  in GST gst+++-- Sidereal time internal functions++-- sidereal 24h correspond to 23:56:04 of solar time+siderealDayLength :: TimeBaseType+siderealDayLength = hours/24+  where C.DH hours = fromHMS 23 56 04.0916+++solarSiderealTimesDiff :: TimeBaseType -> TimeBaseType+solarSiderealTimesDiff d =+  let t = numberOfCenturies j2000 (JD d)+  in reduceToZeroRange 24 $ 6.697374558 + 2400.051336*t + 0.000025862*t*t
+ src/Data/Astro/Types.hs view
@@ -0,0 +1,207 @@+{-|+Module: Data.Astro.Types+Description: Common Types+Copyright: Alexander Ignatyev, 2016++Common Types are usfull across all subsystems like Time and Coordinate.++= Examples++== /Decimal hours and Decimal degrees/++@+import Data.Astro.Types++-- 10h 15m 19.7s+dh :: DecimalHours+dh = fromHMS 10 15 19.7+-- DH 10.255472222222222++(h, m, s) = toHMS dh+-- (10,15,19.699999999999562)+++-- 51°28′40″+dd :: DecimalDegrees+dd = fromDMS 51 28 40+-- DD 51.477777777777774++(d, m, s) = toDMS dd+-- (51,28,39.999999999987494)+@++== /Geographic Coordinates/+@+import Data.Astro.Types++-- the Royal Observatory, Greenwich+ro :: GeographicCoordinates+ro = GeoC (fromDMS 51 28 40) (-(fromDMS 0 0 5))+-- GeoC {geoLatitude = DD 51.4778, geoLongitude = DD (-0.0014)}+@+-}++module Data.Astro.Types+(+  DecimalDegrees(..)+  , DecimalHours (..)+  , GeographicCoordinates(..)+  , AstronomicalUnits(..)+  , lightTravelTime+  , toDecimalHours+  , fromDecimalHours+  , toRadians+  , fromRadians+  , fromDMS+  , toDMS+  , fromHMS+  , toHMS+)++where++import qualified Data.Astro.Utils as U+++newtype DecimalDegrees = DD Double+                         deriving (Show, Eq, Ord)+++instance Num DecimalDegrees where+  (+) (DD d1) (DD d2) = DD (d1+d2)+  (-) (DD d1) (DD d2) = DD (d1-d2)+  (*) (DD d1) (DD d2) = DD (d1*d2)+  negate (DD d) = DD (negate d)+  abs (DD d) = DD (abs d)+  signum (DD d) = DD (signum d)+  fromInteger int = DD (fromInteger int)++instance Real DecimalDegrees where+  toRational (DD d) = toRational d++instance Fractional DecimalDegrees where+  (/) (DD d1) (DD d2) = DD (d1/d2)+  recip (DD d) = DD (recip d)+  fromRational r = DD (fromRational r)++instance RealFrac DecimalDegrees where+  properFraction (DD d) =+    let (i, f) = properFraction d+    in (i, DD f)+++newtype DecimalHours = DH Double+                       deriving (Show, Eq, Ord)+++instance Num DecimalHours where+  (+) (DH d1) (DH d2) = DH (d1+d2)+  (-) (DH d1) (DH d2) = DH (d1-d2)+  (*) (DH d1) (DH d2) = DH (d1*d2)+  negate (DH d) = DH (negate d)+  abs (DH d) = DH (abs d)+  signum (DH d) = DH (signum d)+  fromInteger int = DH (fromInteger int)++instance Real DecimalHours where+  toRational (DH d) = toRational d++instance Fractional DecimalHours where+  (/) (DH d1) (DH d2) = DH (d1/d2)+  recip (DH d) = DH (recip d)+  fromRational r = DH (fromRational r)++instance RealFrac DecimalHours where+  properFraction (DH d) =+    let (i, f) = properFraction d+    in (i, DH f)+++-- | Convert decimal degrees to decimal hours+toDecimalHours :: DecimalDegrees -> DecimalHours+toDecimalHours (DD d) = DH $ d/15  -- 360 / 24 = 15++-- | Convert decimal hours to decimal degrees+fromDecimalHours :: DecimalHours -> DecimalDegrees+fromDecimalHours (DH h) = DD $ h*15+++-- | Geographic Coordinates+data GeographicCoordinates = GeoC {+  geoLatitude :: DecimalDegrees+  , geoLongitude :: DecimalDegrees+  } deriving (Show, Eq)+++-- | Astronomical Units, 1AU = 1.4960×1011 m+-- (originally, the average distance of Earth's aphelion and perihelion).+newtype AstronomicalUnits = AU Double deriving (Show, Eq, Ord)+++instance Num AstronomicalUnits where+  (+) (AU d1) (AU d2) = AU (d1+d2)+  (-) (AU d1) (AU d2) = AU (d1-d2)+  (*) (AU d1) (AU d2) = AU (d1*d2)+  negate (AU d) = AU (negate d)+  abs (AU d) = AU (abs d)+  signum (AU d) = AU (signum d)+  fromInteger int = AU (fromInteger int)++instance Real AstronomicalUnits where+  toRational (AU d) = toRational d++instance Fractional AstronomicalUnits where+  (/) (AU d1) (AU d2) = AU (d1/d2)+  recip (AU d) = AU (recip d)+  fromRational r = AU (fromRational r)++instance RealFrac AstronomicalUnits where+  properFraction (AU d) =+    let (i, f) = properFraction d+    in (i, AU f)+++-- | Light travel time of the distance in Astronomical Units+lightTravelTime :: AstronomicalUnits -> DecimalHours+lightTravelTime (AU ro) = DH $ 0.1386*ro++-- | Convert from DecimalDegrees to Radians+toRadians (DD deg) = U.toRadians deg+++-- | Convert from Radians to DecimalDegrees+fromRadians rad = DD $ U.fromRadians rad+++-- | Convert Degrees, Minutes, Seconds to DecimalDegrees+fromDMS :: RealFrac a => Int -> Int -> a -> DecimalDegrees+fromDMS d m s =+  let d' = fromIntegral d+      m' = fromIntegral m+      s' = realToFrac s+  in DD $ d'+(m'+(s'/60))/60+++-- | Convert DecimalDegrees to Degrees, Minutes, Seconds+toDMS (DD dd) =+  let (d, rm) = properFraction dd+      (m, rs) = properFraction $ 60 * rm+      s = 60 * rs+  in (d, m, s)+++-- | Comvert Hours, Minutes, Seconds to DecimalHours+fromHMS :: RealFrac a => Int -> Int -> a -> DecimalHours+fromHMS h m s =+  let h' = fromIntegral h+      m' = fromIntegral m+      s' = realToFrac s+  in DH $ h'+(m'+(s'/60))/60+++-- | Convert DecimalDegrees to Degrees, Minutes, Seconds+toHMS (DH dh) =+  let (h, rm) = properFraction dh+      (m, rs) = properFraction $ 60 * rm+      s = 60 * rs+  in (h, m, s)
+ src/Data/Astro/Utils.hs view
@@ -0,0 +1,68 @@+{-|+Module: Data.Astro.Utils+Description: Utility functions+Copyright: Alexander Ignatyev, 2016++Utility functions.+-}+++module Data.Astro.Utils+(+  fromFixed+  , trunc+  , fraction+  , reduceToZeroRange+  , toRadians+  , fromRadians+  , roundToN+  , tropicalYearLen+)++where++import Data.Fixed(Fixed(MkFixed), HasResolution(resolution))++-- | Convert From Fixed to Fractional+fromFixed :: (Fractional a, HasResolution b) => Fixed b -> a+fromFixed fv@(MkFixed v) = (fromIntegral v) / (fromIntegral $ resolution fv)+++-- | return the integral part of a number+-- almost the same as truncate but result type is Real+trunc :: RealFrac a => a -> a+trunc = fromIntegral . truncate+++-- | Almost the same the properFraction function but result type+fraction :: (RealFrac a, Num b) => a -> (b, a)+fraction v = let (i, f) = (properFraction v)+             in (fromIntegral i, f)+++-- | Reduce to range from 0 to n+reduceToZeroRange :: RealFrac a => a -> a -> a+reduceToZeroRange r n =+  let b = n - (trunc (n / r)) * r+  in if b < 0 then b + r else b+++-- | Convert from degrees to radians+toRadians :: Floating a => a -> a+toRadians deg = deg*pi/180+++-- | Convert from radians to degrees+fromRadians :: Floating a => a -> a+fromRadians rad = rad*180/pi+++-- | Round to a specified number of digits+roundToN :: RealFrac a => Int -> a -> a+roundToN n f = (fromInteger $ round $ f * factor) / factor+  where factor = 10.0^^n+++-- | Length of a tropical year in days+tropicalYearLen :: Double+tropicalYearLen = 365.242191
+ test/Main.hs view
@@ -0,0 +1,43 @@+import Test.Framework (defaultMain, testGroup)+++import qualified Data.Astro.TimeTest as Time+import qualified Data.Astro.Time.GregorianCalendarTest as Time.GregorianCalendar+import qualified Data.Astro.Time.JulianDateTest as Time.JulianDate+import qualified Data.Astro.Time.SiderealTest as Time.Sidereal+import qualified Data.Astro.Time.ConvTest as Time.Conv+import qualified Data.Astro.CoordinateTest as Coordinate+import qualified Data.Astro.TypesTest as Types+import qualified Data.Astro.UtilsTest as Utils+import qualified Data.Astro.CelestialObjectTest as CelestialObject+import qualified Data.Astro.CelestialObject.RiseSetTest as CelestialObject.RiseSet+import qualified Data.Astro.EffectsTest as Effects+import qualified Data.Astro.Effects.ParallaxTest as Effects.Parallax+import qualified Data.Astro.SunTest as Sun+import qualified Data.Astro.Sun.SunInternalsTest as SunInternals+import qualified Data.Astro.Planet.PlanetDetailsTest as PlanetDetails+import qualified Data.Astro.Planet.PlanetMechanicsTest as PlanetMechanics+import qualified Data.Astro.MoonTest as Moon+++main = defaultMain tests++tests = [+  testGroup "Data.Astro.Time" Time.tests+  , testGroup "Data.Astro.Time.GregorianCalendar" Time.GregorianCalendar.tests+  , testGroup "Data.Astro.Time.JulianDate" Time.JulianDate.tests+  , testGroup "Data.Astro.Time.Sidereal" Time.Sidereal.tests+  , testGroup "Data.Astro.Time.Conv" Time.Conv.tests+  , testGroup "Data.Astro.Coordinate" Coordinate.tests+  , testGroup "Data.Astro.Types" Types.tests+  , testGroup "Data.Astro.Utils" Utils.tests+  , testGroup "Data.Astro.CelestialObject" CelestialObject.tests+  , testGroup "Data.Astro.CelestialObject.RiseSet" CelestialObject.RiseSet.tests+  , testGroup "Data.Astro.Effects" Effects.tests+  , testGroup "Data.Astro.Effects.Parallax" Effects.Parallax.tests+  , testGroup "Data.Astro.Sun" Sun.tests+  , testGroup "Data.Astro.Sun.SunInternals" SunInternals.tests+  , testGroup "Data.Astro.Planet.PlanetDetails" PlanetDetails.tests+  , testGroup "Data.Astro.Planet.PlanetMechanics" PlanetMechanics.tests+  , testGroup "Data.Astro.Moon" Moon.tests+  ]