{-|
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 {
  SunDetails -> JulianDate
sdEpoch :: JulianDate             -- ^ Epoch

  , SunDetails -> DecimalDegrees
sdEpsilon :: DecimalDegrees     -- ^ Ecliptic longitude at the Epoch

  , SunDetails -> DecimalDegrees
sdOmega :: DecimalDegrees       -- ^ Ecliptic longitude of perigee at the Epoch

  , SunDetails -> Double
sdE :: Double                   -- ^ Eccentricity of the orbit at the Epoch

  } deriving (Int -> SunDetails -> ShowS
[SunDetails] -> ShowS
SunDetails -> String
(Int -> SunDetails -> ShowS)
-> (SunDetails -> String)
-> ([SunDetails] -> ShowS)
-> Show SunDetails
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [SunDetails] -> ShowS
$cshowList :: [SunDetails] -> ShowS
show :: SunDetails -> String
$cshow :: SunDetails -> String
showsPrec :: Int -> SunDetails -> ShowS
$cshowsPrec :: Int -> SunDetails -> ShowS
Show)

-- | SunDetails at the Sun's reference Epoch J2010.0

j2010SunDetails :: SunDetails
j2010SunDetails :: SunDetails
j2010SunDetails = JulianDate
-> DecimalDegrees -> DecimalDegrees -> Double -> SunDetails
SunDetails JulianDate
j2010 (Double -> DecimalDegrees
DD Double
279.557208) (Double -> DecimalDegrees
DD Double
283.112438) Double
0.016705


-- | Semi-major axis

r0 :: Double
r0 :: Double
r0 = Double
1.495985e8


-- | Angular diameter at r = r0

theta0 :: DecimalDegrees
theta0 :: DecimalDegrees
theta0 = Double -> DecimalDegrees
DD Double
0.533128


-- | Reduce the value to the range [0, 360)

reduceTo360 :: Double -> Double
reduceTo360 :: Double -> Double
reduceTo360 = Double -> Double -> Double
forall a. RealFrac a => a -> a -> a
U.reduceToZeroRange Double
360


-- | Reduce the value to the range [0, 360)

reduceDegrees :: DecimalDegrees -> DecimalDegrees
reduceDegrees :: DecimalDegrees -> DecimalDegrees
reduceDegrees = DecimalDegrees -> DecimalDegrees -> DecimalDegrees
forall a. RealFrac a => a -> a -> a
U.reduceToZeroRange DecimalDegrees
360


-- | Calculate SunDetails for the given JulianDate.

sunDetails :: JulianDate -> SunDetails
sunDetails :: JulianDate -> SunDetails
sunDetails JulianDate
jd =
  let t :: Double
t = JulianDate -> JulianDate -> Double
numberOfCenturies JulianDate
j1900 JulianDate
jd
      epsilon :: Double
epsilon = Double -> Double
reduceTo360 (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double
279.6966778 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
36000.76892Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
0.0003025Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
tDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
t
      omega :: Double
omega = Double -> Double
reduceTo360 (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double
281.2208444 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1.719175Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
0.000452778Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
tDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
t
      e :: Double
e = Double
0.01675104 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
0.0000418Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
0.000000126Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
tDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
t
  in JulianDate
-> DecimalDegrees -> DecimalDegrees -> Double -> SunDetails
SunDetails JulianDate
jd (Double -> DecimalDegrees
DD Double
epsilon) (Double -> DecimalDegrees
DD Double
omega) Double
e


-- | Calculate the ecliptic longitude of the Sun with the given SunDetails at the given JulianDate

sunEclipticLongitude1 :: SunDetails -> JulianDate -> DecimalDegrees
sunEclipticLongitude1 :: SunDetails -> JulianDate -> DecimalDegrees
sunEclipticLongitude1 sd :: SunDetails
sd@(SunDetails JulianDate
epoch (DD Double
eps) (DD Double
omega) Double
e) JulianDate
jd =
  let d :: Double
d = JulianDate -> JulianDate -> Double
numberOfDays JulianDate
epoch JulianDate
jd
      n :: Double
n = Double -> Double
reduceTo360 (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ (Double
360Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
U.tropicalYearLen) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d
      meanAnomaly :: Double
meanAnomaly = Double -> Double
reduceTo360 (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
eps Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
omega
      ec :: Double
ec = (Double
360Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
forall a. Floating a => a
pi)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
eDouble -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
sin (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double -> Double
forall a. Floating a => a -> a
U.toRadians Double
meanAnomaly)
      DD Double
nutation = JulianDate -> DecimalDegrees
nutationLongitude JulianDate
jd
  in Double -> DecimalDegrees
DD (Double -> DecimalDegrees) -> Double -> DecimalDegrees
forall a b. (a -> b) -> a -> b
$ Double -> Double
reduceTo360 (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
ec Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
eps Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
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 :: SunDetails -> JulianDate -> EquatorialCoordinates1
sunPosition1 SunDetails
sd JulianDate
jd =
  let lambda :: DecimalDegrees
lambda = SunDetails -> JulianDate -> DecimalDegrees
sunEclipticLongitude1 SunDetails
sd JulianDate
jd
      beta :: DecimalDegrees
beta = Double -> DecimalDegrees
DD Double
0
  in EclipticCoordinates -> JulianDate -> EquatorialCoordinates1
eclipticToEquatorial (DecimalDegrees -> DecimalDegrees -> EclipticCoordinates
EcC DecimalDegrees
beta DecimalDegrees
lambda) JulianDate
jd


-- | Calculate mean anomaly using the second 'more accurate' method

sunMeanAnomaly2 :: SunDetails -> DecimalDegrees
sunMeanAnomaly2 :: SunDetails -> DecimalDegrees
sunMeanAnomaly2 SunDetails
sd = DecimalDegrees -> DecimalDegrees
reduceDegrees (DecimalDegrees -> DecimalDegrees)
-> DecimalDegrees -> DecimalDegrees
forall a b. (a -> b) -> a -> b
$ (SunDetails -> DecimalDegrees
sdEpsilon SunDetails
sd) DecimalDegrees -> DecimalDegrees -> DecimalDegrees
forall a. Num a => a -> a -> a
- (SunDetails -> DecimalDegrees
sdOmega SunDetails
sd)


-- | Calculate true anomaly using the second 'more accurate' method

trueAnomaly2 :: SunDetails -> DecimalDegrees
trueAnomaly2 :: SunDetails -> DecimalDegrees
trueAnomaly2 SunDetails
sd =
  let m :: Double
m = DecimalDegrees -> Double
toRadians (DecimalDegrees -> Double) -> DecimalDegrees -> Double
forall a b. (a -> b) -> a -> b
$ SunDetails -> DecimalDegrees
sunMeanAnomaly2 SunDetails
sd
      e :: Double
e = SunDetails -> Double
sdE SunDetails
sd
      bigE :: Double
bigE = Double -> Double -> Double -> Double
solveKeplerEquation Double
e Double
m Double
0.000000001
      tanHalfNu :: Double
tanHalfNu = Double -> Double
forall a. Floating a => a -> a
sqrt((Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
e)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
e)) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
tan (Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
bigE)
      nu :: Double
nu = Double -> Double
reduceTo360 (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double -> Double
forall a. Floating a => a -> a
U.fromRadians (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double -> Double
forall a. Floating a => a -> a
atan Double
tanHalfNu)
  in Double -> DecimalDegrees
DD Double
nu


-- | Calculate the ecliptic longitude of the Sun

sunEclipticLongitude2 :: SunDetails -> DecimalDegrees
sunEclipticLongitude2 :: SunDetails -> DecimalDegrees
sunEclipticLongitude2 SunDetails
sd =
  let DD Double
omega = SunDetails -> DecimalDegrees
sdOmega SunDetails
sd
      DD Double
nu = SunDetails -> DecimalDegrees
trueAnomaly2 SunDetails
sd
      DD Double
nutation = JulianDate -> DecimalDegrees
nutationLongitude (JulianDate -> DecimalDegrees) -> JulianDate -> DecimalDegrees
forall a b. (a -> b) -> a -> b
$ SunDetails -> JulianDate
sdEpoch SunDetails
sd
  in Double -> DecimalDegrees
DD (Double -> DecimalDegrees) -> Double -> DecimalDegrees
forall a b. (a -> b) -> a -> b
$ Double -> Double
reduceTo360 (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double
nu Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
omega Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
nutation


-- | More accurate method to calculate position of the Sun

sunPosition2 :: JulianDate -> EquatorialCoordinates1
sunPosition2 :: JulianDate -> EquatorialCoordinates1
sunPosition2 JulianDate
jd =
  let sd :: SunDetails
sd = JulianDate -> SunDetails
sunDetails JulianDate
jd
      lambda :: DecimalDegrees
lambda = SunDetails -> DecimalDegrees
sunEclipticLongitude2 SunDetails
sd
      beta :: DecimalDegrees
beta = Double -> DecimalDegrees
DD Double
0
  in EclipticCoordinates -> JulianDate -> EquatorialCoordinates1
eclipticToEquatorial (DecimalDegrees -> DecimalDegrees -> EclipticCoordinates
EcC DecimalDegrees
beta DecimalDegrees
lambda) JulianDate
jd


-- Distance and Angular Size helper function

dasf :: SunDetails -> Double
dasf SunDetails
sd =
  let e :: Double
e = SunDetails -> Double
sdE SunDetails
sd
      nu :: Double
nu = DecimalDegrees -> Double
toRadians (DecimalDegrees -> Double) -> DecimalDegrees -> Double
forall a b. (a -> b) -> a -> b
$ SunDetails -> DecimalDegrees
trueAnomaly2 SunDetails
sd
  in (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
eDouble -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
cos Double
nu)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
eDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
e)


-- | Calculate Sun-Earth distance.

sunDistance :: JulianDate -> Double
sunDistance :: JulianDate -> Double
sunDistance JulianDate
jd = Double
r0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (SunDetails -> Double
dasf (SunDetails -> Double) -> SunDetails -> Double
forall a b. (a -> b) -> a -> b
$ JulianDate -> SunDetails
sunDetails JulianDate
jd)


-- | Calculate the Sun's angular size (i.e. its angular diameter).

sunAngularSize :: JulianDate -> DecimalDegrees
sunAngularSize :: JulianDate -> DecimalDegrees
sunAngularSize JulianDate
jd = DecimalDegrees
theta0 DecimalDegrees -> DecimalDegrees -> DecimalDegrees
forall a. Num a => a -> a -> a
* (Double -> DecimalDegrees
DD (Double -> DecimalDegrees) -> Double -> DecimalDegrees
forall a b. (a -> b) -> a -> b
$ SunDetails -> Double
dasf (SunDetails -> Double) -> SunDetails -> Double
forall a b. (a -> b) -> a -> b
$ JulianDate -> SunDetails
sunDetails JulianDate
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 :: GeographicCoordinates
-> DecimalDegrees -> LocalCivilDate -> RiseSetMB
sunRiseAndSet = DecimalHours
-> (JulianDate -> EquatorialCoordinates1)
-> GeographicCoordinates
-> DecimalDegrees
-> LocalCivilDate
-> RiseSetMB
riseAndSet2 DecimalHours
0.000001 (SunDetails -> JulianDate -> EquatorialCoordinates1
sunPosition1 SunDetails
j2010SunDetails)


-- | Calculates discrepancy between the mean solar time and real solar time

-- at the given date.

equationOfTime :: JulianDate -> DecimalHours
equationOfTime :: JulianDate -> DecimalHours
equationOfTime JulianDate
jd =
  let (JulianDate
day, JulianDate
_) = JulianDate -> (JulianDate, JulianDate)
splitToDayAndTime JulianDate
jd
      midday :: JulianDate
midday = DecimalHours -> JulianDate -> JulianDate
addHours (Double -> DecimalHours
DH Double
12) JulianDate
day  -- mean solar time

      EC1 DecimalDegrees
_ DecimalHours
ra = SunDetails -> JulianDate -> EquatorialCoordinates1
sunPosition1 SunDetails
j2010SunDetails JulianDate
midday
      ut :: JulianDate
ut = JulianDate -> GreenwichSiderealTime -> JulianDate
gstToUT JulianDate
day (GreenwichSiderealTime -> JulianDate)
-> GreenwichSiderealTime -> JulianDate
forall a b. (a -> b) -> a -> b
$ DecimalHours -> GreenwichSiderealTime
dhToGST DecimalHours
ra
      JD Double
time = JulianDate
midday JulianDate -> JulianDate -> JulianDate
forall a. Num a => a -> a -> a
- JulianDate
ut
  in Double -> DecimalHours
DH (Double -> DecimalHours) -> Double -> DecimalHours
forall a b. (a -> b) -> a -> b
$ Double
timeDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
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 :: EquatorialCoordinates1 -> JulianDate -> DecimalDegrees
solarElongation (EC1 DecimalDegrees
deltaP DecimalHours
alphaP) JulianDate
jd =
  let (EC1 DecimalDegrees
deltaS DecimalHours
alphaS) = SunDetails -> JulianDate -> EquatorialCoordinates1
sunPosition1 SunDetails
j2010SunDetails JulianDate
jd
      deltaP' :: Double
deltaP' = DecimalDegrees -> Double
toRadians DecimalDegrees
deltaP
      alphaP' :: Double
alphaP' = DecimalDegrees -> Double
toRadians (DecimalDegrees -> Double) -> DecimalDegrees -> Double
forall a b. (a -> b) -> a -> b
$ DecimalHours -> DecimalDegrees
fromDecimalHours DecimalHours
alphaP
      deltaS' :: Double
deltaS' = DecimalDegrees -> Double
toRadians DecimalDegrees
deltaS
      alphaS' :: Double
alphaS' = DecimalDegrees -> Double
toRadians (DecimalDegrees -> Double) -> DecimalDegrees -> Double
forall a b. (a -> b) -> a -> b
$ DecimalHours -> DecimalDegrees
fromDecimalHours DecimalHours
alphaS
      eps :: Double
eps = Double -> Double
forall a. Floating a => a -> a
acos (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ (Double -> Double
forall a. Floating a => a -> a
sin Double
deltaP')Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
sin Double
deltaS') Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double -> Double
forall a. Floating a => a -> a
cos (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double
alphaP' Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
alphaS')Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
cos Double
deltaP')Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
cos Double
deltaS')
  in Double -> DecimalDegrees
fromRadians Double
eps