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)
data SunDetails = SunDetails {
sdEpoch :: JulianDate
, sdEpsilon :: DecimalDegrees
, sdOmega :: DecimalDegrees
, sdE :: Double
} deriving (Show)
j2010SunDetails :: SunDetails
j2010SunDetails = SunDetails j2010 (DD 279.557208) (DD 283.112438) 0.016705
r0 :: Double
r0 = 1.495985e8
theta0 :: DecimalDegrees
theta0 = DD 0.533128
reduceTo360 :: Double -> Double
reduceTo360 = U.reduceToZeroRange 360
reduceDegrees :: DecimalDegrees -> DecimalDegrees
reduceDegrees = U.reduceToZeroRange 360
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
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
sunPosition1 :: SunDetails -> JulianDate -> EquatorialCoordinates1
sunPosition1 sd jd =
let lambda = sunEclipticLongitude1 sd jd
beta = DD 0
in eclipticToEquatorial (EcC beta lambda) jd
sunMeanAnomaly2 :: SunDetails -> DecimalDegrees
sunMeanAnomaly2 sd = reduceDegrees $ (sdEpsilon sd) - (sdOmega sd)
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
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
sunPosition2 :: JulianDate -> EquatorialCoordinates1
sunPosition2 jd =
let sd = sunDetails jd
lambda = sunEclipticLongitude2 sd
beta = DD 0
in eclipticToEquatorial (EcC beta lambda) jd
dasf sd =
let e = sdE sd
nu = toRadians $ trueAnomaly2 sd
in (1 + e*(cos nu)) / (1 - e*e)
sunDistance :: JulianDate -> Double
sunDistance jd = r0 / (dasf $ sunDetails jd)
sunAngularSize :: JulianDate -> DecimalDegrees
sunAngularSize jd = theta0 * (DD $ dasf $ sunDetails jd)
sunRiseAndSet :: GeographicCoordinates
-> DecimalDegrees
-> LocalCivilDate
-> RiseSetMB
sunRiseAndSet = riseAndSet2 0.000001 (sunPosition1 j2010SunDetails)
equationOfTime :: JulianDate -> DecimalHours
equationOfTime jd =
let (day, _) = splitToDayAndTime jd
midday = addHours (DH 12) day
EC1 _ ra = sunPosition1 j2010SunDetails midday
ut = gstToUT day $ dhToGST ra
JD time = midday - ut
in DH $ time*24
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