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 {
SunDetails -> JulianDate
sdEpoch :: JulianDate
, SunDetails -> DecimalDegrees
sdEpsilon :: DecimalDegrees
, SunDetails -> DecimalDegrees
sdOmega :: DecimalDegrees
, SunDetails -> Double
sdE :: Double
} 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)
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
r0 :: Double
r0 :: Double
r0 = Double
1.495985e8
theta0 :: DecimalDegrees
theta0 :: DecimalDegrees
theta0 = Double -> DecimalDegrees
DD Double
0.533128
reduceTo360 :: Double -> Double
reduceTo360 :: Double -> Double
reduceTo360 = Double -> Double -> Double
forall a. RealFrac a => a -> a -> a
U.reduceToZeroRange Double
360
reduceDegrees :: DecimalDegrees -> DecimalDegrees
reduceDegrees :: DecimalDegrees -> DecimalDegrees
reduceDegrees = DecimalDegrees -> DecimalDegrees -> DecimalDegrees
forall a. RealFrac a => a -> a -> a
U.reduceToZeroRange DecimalDegrees
360
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
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
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
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)
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
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
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
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)
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)
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)
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)
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
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
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