diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/app/Main.hs b/app/Main.hs new file mode 100644 --- /dev/null +++ b/app/Main.hs @@ -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 diff --git a/astro.cabal b/astro.cabal new file mode 100644 --- /dev/null +++ b/astro.cabal @@ -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 diff --git a/src/Data/Astro/CelestialObject.hs b/src/Data/Astro/CelestialObject.hs new file mode 100644 --- /dev/null +++ b/src/Data/Astro/CelestialObject.hs @@ -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 diff --git a/src/Data/Astro/CelestialObject/RiseSet.hs b/src/Data/Astro/CelestialObject/RiseSet.hs new file mode 100644 --- /dev/null +++ b/src/Data/Astro/CelestialObject/RiseSet.hs @@ -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) diff --git a/src/Data/Astro/Coordinate.hs b/src/Data/Astro/Coordinate.hs new file mode 100644 --- /dev/null +++ b/src/Data/Astro/Coordinate.hs @@ -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, α__ - /'how far up'/ angle from the horizontal plane in degrees +* __azimuth, Α__ - /'how far round'/ agle from the north direction in degrees to the east + + +== /Equatorial coordinates/ + +Accoring to the equatorial coordinates system stars move westwards along the circles centered in the north selestial pole, +making the full cicrle in 24 hours of sidereal time (see "Data.Astro.Time.Sidereal"). + +* __declination, δ__ - /'how far up'/ angle from the quatorial plane; +* __right ascension, α__ - /'how far round'/ angle from the /vernal equinox/ to the east; __/or/__ +* __hour angle__ - /'how far round'/ angle from the meridian to the west + + +== /Ecliptic Coordinate/ + +Accoring to the ecliptic coordinates system the Sun moves eastwards along the trace of th ecliptic. The Sun's ecplitic latitude is always 0. + +* __ecliptic latitude, β__ - /'how far up'/ angle from the ecliptic +* __ecliptic longitude, λ__ - /'how far round'/ angle from the /vernal equinox/ to the east + + +== /Galactic Coordinates/ + +* __galactic latitute, b__ - /'how far up'/ angle from the plane of the Galaxy +* __galactiv longitude, l__ - - /'how far round'/ angle from the direction the Sun - the centre of the Galaxy + + +== /Terms/ + +* __ecliptic__ - the plane containing the Earth's orbit around the Sun +* __vernal equinox__, ♈ - fixed direction lies along the line of the intersection of the equatorial plane and the ecliptic +* __obliquity of the ecliptic, β__ - the angle between the plane of the Earth's equator and the ecliptic +* __north selestial pole, P__ - the point on the selestial sphere, right above the Earth's North Pole + + += Examples + +== /Horizontal Coordinate System/ +@ +import Data.Astro.Coordinate +import Data.Astro.Types + +hc :: HorizonCoordinates +hc = HC (DD 30.5) (DD 180) +-- HC {hAltitude = DD 30.0, hAzimuth = DD 180.0} +@ + +== /Equatorial Coordinate System/ +@ +import Data.Astro.Coordinate +import Data.Astro.Types + +ec1 :: EquatorialCoordinates1 +ec1 = EC1 (DD 71.7) (DH 8) +-- EC1 {e1Declination = DD 71.7, e1RightAscension = DH 8.0} + +ec2 :: EquatorialCoordinates2 +ec2 = EC1 (DD 77.7) (DH 11) +-- EC2 {e2Declination = DD 77.7, e2HoursAngle = DH 11.0} +@ + +== /Transformations/ +@ +import Data.Astro.Time.JulianDate +import Data.Astro.Coordinate +import Data.Astro.Types + +ro :: GeographicCoordinates +ro = GeoC (fromDMS 51 28 40) (-(fromDMS 0 0 5)) + +dt :: LocalCivilTime +dt = lctFromYMDHMS (DH 1) 2017 6 25 10 29 0 + +sunHC :: HorizonCoordinates +sunHC = HC (fromDMS 49 18 21.77) (fromDMS 118 55 19.53) +-- HC {hAltitude = DD 49.30604722222222, hAzimuth = DD 118.92209166666666} + +sunEC2 :: EquatorialCoordinates2 +sunEC2 = horizonToEquatorial (geoLatitude ro) sunHC +-- EC2 {e2Declination = DD 23.378295912623855, e2HoursAngle = DH 21.437117068873537} + +sunEC1 :: EquatorialCoordinates1 +sunEC1 = EC1 (e2Declination sunEC2) (haToRA (e2HoursAngle sunEC2) (geoLongitude ro) (lctUniversalTime dt)) +-- EC1 {e1Declination = DD 23.378295912623855, e1RightAscension = DH 6.29383725890224} + + +sunEC2' :: EquatorialCoordinates2 +sunEC2' = EC2 (e1Declination sunEC1) (raToHA (e1RightAscension sunEC1) (geoLongitude ro) (lctUniversalTime dt)) +-- EC2 {e2Declination = DD 23.378295912623855, e2HoursAngle = DH 21.437117068873537} + +sunHC' :: HorizonCoordinates +sunHC' = equatorialToHorizon (geoLatitude ro) sunEC2' +-- HC {hAltitude = DD 49.30604722222222, hAzimuth = DD 118.92209166666666} +@ + +=== /Function-shortcuts/ + +@ +import Data.Astro.Time.JulianDate +import Data.Astro.Coordinate +import Data.Astro.Types + +ro :: GeographicCoordinates +ro = GeoC (fromDMS 51 28 40) (-(fromDMS 0 0 5)) + +dt :: LocalCivilTime +dt = lctFromYMDHMS (DH 1) 2017 6 25 10 29 0 + +sunHC :: HorizonCoordinates +sunHC = HC (fromDMS 49 18 21.77) (fromDMS 118 55 19.53) +-- HC {hAltitude = DD 49.30604722222222, hAzimuth = DD 118.92209166666666} + +sunEC1 :: EquatorialCoordinates1 +sunEC1 = hcToEC1 ro (lctUniversalTime dt) sunHC +-- EC1 {e1Declination = DD 23.378295912623855, e1RightAscension = DH 6.29383725890224} + +sunHC' :: HorizonCoordinates +sunHC' = ec1ToHC ro (lctUniversalTime dt) sunEC1 +-- HC {hAltitude = DD 49.30604722222222, hAzimuth = DD 118.92209166666666} +@ +-} + +module Data.Astro.Coordinate +( + DecimalDegrees(..) + , DecimalHours(..) + , HorizonCoordinates(..) + , EquatorialCoordinates1(..) + , EquatorialCoordinates2(..) + , EclipticCoordinates(..) + , GalacticCoordinates(..) + , raToHA + , haToRA + , equatorialToHorizon + , horizonToEquatorial + , ec1ToHC + , hcToEC1 + , ecHCConv + , obliquity + , eclipticToEquatorial + , equatorialToEcliptic + , galacticToEquatorial + , equatorialToGalactic +) + +where + +import Data.Astro.Time (utToLST) +import Data.Astro.Time.JulianDate (JulianDate(..), numberOfCenturies, splitToDayAndTime) +import Data.Astro.Time.Epoch (j2000) +import Data.Astro.Time.Sidereal (LocalSiderealTime(..), lstToDH) +import Data.Astro.Types (DecimalDegrees(..), DecimalHours(..) + , fromDecimalHours, toDecimalHours + , toRadians, fromRadians, fromDMS + , GeographicCoordinates(..)) +import Data.Astro.Utils (fromFixed) +import Data.Astro.Effects.Nutation (nutationObliquity) + + +-- | Horizon Coordinates, for details see the module's description +data HorizonCoordinates = HC { + 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 diff --git a/src/Data/Astro/Effects.hs b/src/Data/Astro/Effects.hs new file mode 100644 --- /dev/null +++ b/src/Data/Astro/Effects.hs @@ -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)) diff --git a/src/Data/Astro/Effects/Aberration.hs b/src/Data/Astro/Effects/Aberration.hs new file mode 100644 --- /dev/null +++ b/src/Data/Astro/Effects/Aberration.hs @@ -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) diff --git a/src/Data/Astro/Effects/Nutation.hs b/src/Data/Astro/Effects/Nutation.hs new file mode 100644 --- /dev/null +++ b/src/Data/Astro/Effects/Nutation.hs @@ -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 diff --git a/src/Data/Astro/Effects/Parallax.hs b/src/Data/Astro/Effects/Parallax.hs new file mode 100644 --- /dev/null +++ b/src/Data/Astro/Effects/Parallax.hs @@ -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) diff --git a/src/Data/Astro/Effects/Precession.hs b/src/Data/Astro/Effects/Precession.hs new file mode 100644 --- /dev/null +++ b/src/Data/Astro/Effects/Precession.hs @@ -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 diff --git a/src/Data/Astro/Moon.hs b/src/Data/Astro/Moon.hs new file mode 100644 --- /dev/null +++ b/src/Data/Astro/Moon.hs @@ -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' diff --git a/src/Data/Astro/Moon/MoonDetails.hs b/src/Data/Astro/Moon/MoonDetails.hs new file mode 100644 --- /dev/null +++ b/src/Data/Astro/Moon/MoonDetails.hs @@ -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) diff --git a/src/Data/Astro/Planet.hs b/src/Data/Astro/Planet.hs new file mode 100644 --- /dev/null +++ b/src/Data/Astro/Planet.hs @@ -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 diff --git a/src/Data/Astro/Planet/PlanetDetails.hs b/src/Data/Astro/Planet/PlanetDetails.hs new file mode 100644 --- /dev/null +++ b/src/Data/Astro/Planet/PlanetDetails.hs @@ -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 diff --git a/src/Data/Astro/Planet/PlanetMechanics.hs b/src/Data/Astro/Planet/PlanetMechanics.hs new file mode 100644 --- /dev/null +++ b/src/Data/Astro/Planet/PlanetMechanics.hs @@ -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 diff --git a/src/Data/Astro/Star.hs b/src/Data/Astro/Star.hs new file mode 100644 --- /dev/null +++ b/src/Data/Astro/Star.hs @@ -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) diff --git a/src/Data/Astro/Sun.hs b/src/Data/Astro/Sun.hs new file mode 100644 --- /dev/null +++ b/src/Data/Astro/Sun.hs @@ -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 diff --git a/src/Data/Astro/Sun/SunInternals.hs b/src/Data/Astro/Sun/SunInternals.hs new file mode 100644 --- /dev/null +++ b/src/Data/Astro/Sun/SunInternals.hs @@ -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) diff --git a/src/Data/Astro/Time.hs b/src/Data/Astro/Time.hs new file mode 100644 --- /dev/null +++ b/src/Data/Astro/Time.hs @@ -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 diff --git a/src/Data/Astro/Time/Conv.hs b/src/Data/Astro/Time/Conv.hs new file mode 100644 --- /dev/null +++ b/src/Data/Astro/Time/Conv.hs @@ -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 } diff --git a/src/Data/Astro/Time/Epoch.hs b/src/Data/Astro/Time/Epoch.hs new file mode 100644 --- /dev/null +++ b/src/Data/Astro/Time/Epoch.hs @@ -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 diff --git a/src/Data/Astro/Time/GregorianCalendar.hs b/src/Data/Astro/Time/GregorianCalendar.hs new file mode 100644 --- /dev/null +++ b/src/Data/Astro/Time/GregorianCalendar.hs @@ -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 diff --git a/src/Data/Astro/Time/JulianDate.hs b/src/Data/Astro/Time/JulianDate.hs new file mode 100644 --- /dev/null +++ b/src/Data/Astro/Time/JulianDate.hs @@ -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 diff --git a/src/Data/Astro/Time/Sidereal.hs b/src/Data/Astro/Time/Sidereal.hs new file mode 100644 --- /dev/null +++ b/src/Data/Astro/Time/Sidereal.hs @@ -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 diff --git a/src/Data/Astro/Types.hs b/src/Data/Astro/Types.hs new file mode 100644 --- /dev/null +++ b/src/Data/Astro/Types.hs @@ -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) diff --git a/src/Data/Astro/Utils.hs b/src/Data/Astro/Utils.hs new file mode 100644 --- /dev/null +++ b/src/Data/Astro/Utils.hs @@ -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 diff --git a/test/Main.hs b/test/Main.hs new file mode 100644 --- /dev/null +++ b/test/Main.hs @@ -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 + ]