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(..))
data RiseSet a
= RiseSet a a
| Circumpolar
| NeverRises
deriving (Int -> RiseSet a -> ShowS
[RiseSet a] -> ShowS
RiseSet a -> String
(Int -> RiseSet a -> ShowS)
-> (RiseSet a -> String)
-> ([RiseSet a] -> ShowS)
-> Show (RiseSet a)
forall a. Show a => Int -> RiseSet a -> ShowS
forall a. Show a => [RiseSet a] -> ShowS
forall a. Show a => RiseSet a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [RiseSet a] -> ShowS
$cshowList :: forall a. Show a => [RiseSet a] -> ShowS
show :: RiseSet a -> String
$cshow :: forall a. Show a => RiseSet a -> String
showsPrec :: Int -> RiseSet a -> ShowS
$cshowsPrec :: forall a. Show a => Int -> RiseSet a -> ShowS
Show, RiseSet a -> RiseSet a -> Bool
(RiseSet a -> RiseSet a -> Bool)
-> (RiseSet a -> RiseSet a -> Bool) -> Eq (RiseSet a)
forall a. Eq a => RiseSet a -> RiseSet a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: RiseSet a -> RiseSet a -> Bool
$c/= :: forall a. Eq a => RiseSet a -> RiseSet a -> Bool
== :: RiseSet a -> RiseSet a -> Bool
$c== :: forall a. Eq a => RiseSet a -> RiseSet a -> Bool
Eq)
type RSInfo a = (a, DecimalDegrees)
type RiseSetLST = RiseSet (RSInfo LocalSiderealTime)
type RiseSetLCT = RiseSet (RSInfo LocalCivilTime)
type RiseSetMB = RiseSet (Maybe (RSInfo LocalCivilTime))
riseAndSet :: EquatorialCoordinates1 -> DecimalDegrees -> DecimalDegrees -> RiseSetLST
riseAndSet :: EquatorialCoordinates1
-> DecimalDegrees -> DecimalDegrees -> RiseSetLST
riseAndSet (EC1 DecimalDegrees
delta DecimalHours
alpha) DecimalDegrees
shift DecimalDegrees
lat =
let delta' :: Double
delta' = DecimalDegrees -> Double
toRadians DecimalDegrees
delta
shift' :: Double
shift' = DecimalDegrees -> Double
toRadians DecimalDegrees
shift
lat' :: Double
lat' = DecimalDegrees -> Double
toRadians DecimalDegrees
lat
cosH :: Double
cosH = Double -> Double -> Double -> Double
cosOfHourAngle Double
delta' Double
shift' Double
lat'
in Double -> Double -> Double -> Double -> RiseSetLST
sortRiseSet Double
cosH Double
delta' Double
shift' Double
lat'
where sortRiseSet :: Double -> Double -> Double -> Double -> RiseSetLST
sortRiseSet :: Double -> Double -> Double -> Double -> RiseSetLST
sortRiseSet Double
cosH Double
delta Double
shift Double
latitude
| Double
cosH Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< -Double
1 = RiseSetLST
forall a. RiseSet a
Circumpolar
| Double
cosH Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1 = RiseSetLST
forall a. RiseSet a
NeverRises
| Bool
otherwise = DecimalHours
-> DecimalHours -> Double -> Double -> Double -> RiseSetLST
calcTimesAndAzimuths DecimalHours
alpha (Double -> DecimalHours
toHours (Double -> DecimalHours) -> Double -> DecimalHours
forall a b. (a -> b) -> a -> b
$ Double -> Double
forall a. Floating a => a -> a
acos Double
cosH) Double
delta Double
shift Double
latitude
toHours :: Double -> DecimalHours
toHours :: Double -> DecimalHours
toHours = DecimalDegrees -> DecimalHours
toDecimalHours (DecimalDegrees -> DecimalHours)
-> (Double -> DecimalDegrees) -> Double -> DecimalHours
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> DecimalDegrees
fromRadians
cosOfHourAngle :: Double -> Double -> Double -> Double
cosOfHourAngle :: Double -> Double -> Double -> Double
cosOfHourAngle Double
delta Double
shift Double
latitude = -((Double -> Double
forall a. Floating a => a -> a
sin Double
shift) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double -> Double
forall a. Floating a => a -> a
sin Double
latitude)Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
sin Double
delta)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ ((Double -> Double
forall a. Floating a => a -> a
cos Double
latitude)Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
cos Double
delta))
calcTimesAndAzimuths :: DecimalHours -> DecimalHours -> Double -> Double -> Double -> RiseSetLST
calcTimesAndAzimuths :: DecimalHours
-> DecimalHours -> Double -> Double -> Double -> RiseSetLST
calcTimesAndAzimuths DecimalHours
alpha DecimalHours
hourAngle Double
delta Double
shift Double
latitude =
let lstRise :: LocalSiderealTime
lstRise = DecimalHours -> LocalSiderealTime
dhToLST (DecimalHours -> LocalSiderealTime)
-> DecimalHours -> LocalSiderealTime
forall a b. (a -> b) -> a -> b
$ DecimalHours -> DecimalHours -> DecimalHours
forall a. RealFrac a => a -> a -> a
reduceToZeroRange DecimalHours
24 (DecimalHours -> DecimalHours) -> DecimalHours -> DecimalHours
forall a b. (a -> b) -> a -> b
$ DecimalHours
alpha DecimalHours -> DecimalHours -> DecimalHours
forall a. Num a => a -> a -> a
- DecimalHours
hourAngle
lstSet :: LocalSiderealTime
lstSet = DecimalHours -> LocalSiderealTime
dhToLST (DecimalHours -> LocalSiderealTime)
-> DecimalHours -> LocalSiderealTime
forall a b. (a -> b) -> a -> b
$ DecimalHours -> DecimalHours -> DecimalHours
forall a. RealFrac a => a -> a -> a
reduceToZeroRange DecimalHours
24 (DecimalHours -> DecimalHours) -> DecimalHours -> DecimalHours
forall a b. (a -> b) -> a -> b
$ DecimalHours
alpha DecimalHours -> DecimalHours -> DecimalHours
forall a. Num a => a -> a -> a
+ DecimalHours
hourAngle
azimuthRise :: Double
azimuthRise = Double -> Double -> Double
forall a. RealFrac a => a -> a -> a
reduceToZeroRange (Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall a. Floating a => a
pi) (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ 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
delta) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double -> Double
forall a. Floating a => a -> a
sin Double
shift)Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
sin Double
latitude)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ ((Double -> Double
forall a. Floating a => a -> a
cos Double
shift)Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double -> Double
forall a. Floating a => a -> a
cos Double
latitude))
azimuthSet :: Double
azimuthSet = Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
azimuthRise
in (LocalSiderealTime, DecimalDegrees)
-> (LocalSiderealTime, DecimalDegrees) -> RiseSetLST
forall a. a -> a -> RiseSet a
RiseSet (LocalSiderealTime
lstRise, Double -> DecimalDegrees
fromRadians Double
azimuthRise) (LocalSiderealTime
lstSet, Double -> DecimalDegrees
fromRadians Double
azimuthSet)
riseAndSet2 :: DecimalHours
-> (JulianDate -> EquatorialCoordinates1)
-> GeographicCoordinates
-> DecimalDegrees
-> LocalCivilDate
-> RiseSetMB
riseAndSet2 :: DecimalHours
-> (JulianDate -> EquatorialCoordinates1)
-> GeographicCoordinates
-> DecimalDegrees
-> LocalCivilDate
-> RiseSetMB
riseAndSet2 DecimalHours
eps JulianDate -> EquatorialCoordinates1
getPosition GeographicCoordinates
geoc DecimalDegrees
shift LocalCivilDate
lcd =
let day :: JulianDate
day = LocalCivilDate -> JulianDate
lcdDate LocalCivilDate
lcd
pos :: EquatorialCoordinates1
pos = JulianDate -> EquatorialCoordinates1
getPosition (DecimalHours -> JulianDate -> JulianDate
addHours DecimalHours
12 JulianDate
day)
rs :: RiseSetLCT
rs = GeographicCoordinates
-> LocalCivilDate
-> DecimalDegrees
-> EquatorialCoordinates1
-> RiseSetLCT
riseAndSetLCT GeographicCoordinates
geoc LocalCivilDate
lcd DecimalDegrees
shift EquatorialCoordinates1
pos
rise :: RiseSetLCT
rise = (RiseSetLCT -> RSInfo LocalCivilTime)
-> RSInfo LocalCivilTime -> Int -> RiseSetLCT
calc RiseSetLCT -> RSInfo LocalCivilTime
getRiseTime (RiseSetLCT -> RSInfo LocalCivilTime
getRiseTime RiseSetLCT
rs) Int
0
set :: RiseSetLCT
set = (RiseSetLCT -> RSInfo LocalCivilTime)
-> RSInfo LocalCivilTime -> Int -> RiseSetLCT
calc RiseSetLCT -> RSInfo LocalCivilTime
getSetTime (RiseSetLCT -> RSInfo LocalCivilTime
getSetTime RiseSetLCT
rs) Int
0
in case RiseSetLCT
rs of
RiseSetLCT
Circumpolar -> RiseSetMB
forall a. RiseSet a
Circumpolar
RiseSetLCT
NeverRises -> RiseSetMB
forall a. RiseSet a
NeverRises
RiseSetLCT
_ -> RiseSetLCT -> RiseSetLCT -> RiseSetMB
forall a. RiseSet a -> RiseSet a -> RiseSet (Maybe a)
buildResult RiseSetLCT
rise RiseSetLCT
set
where calc :: (RiseSetLCT -> RSInfo LocalCivilTime) -> RSInfo LocalCivilTime -> Int -> RiseSetLCT
calc :: (RiseSetLCT -> RSInfo LocalCivilTime)
-> RSInfo LocalCivilTime -> Int -> RiseSetLCT
calc RiseSetLCT -> RSInfo LocalCivilTime
getRSInfo rsi :: RSInfo LocalCivilTime
rsi@(LocalCivilTime
time, DecimalDegrees
_) Int
iterNo =
let pos :: EquatorialCoordinates1
pos = JulianDate -> EquatorialCoordinates1
getPosition (JulianDate -> EquatorialCoordinates1)
-> JulianDate -> EquatorialCoordinates1
forall a b. (a -> b) -> a -> b
$ LocalCivilTime -> JulianDate
lctUniversalTime LocalCivilTime
time
rs :: RiseSetLCT
rs = GeographicCoordinates
-> LocalCivilDate
-> DecimalDegrees
-> EquatorialCoordinates1
-> RiseSetLCT
riseAndSetLCT GeographicCoordinates
geoc LocalCivilDate
lcd DecimalDegrees
shift EquatorialCoordinates1
pos
rsi' :: RSInfo LocalCivilTime
rsi' = RiseSetLCT -> RSInfo LocalCivilTime
getRSInfo RiseSetLCT
rs
in case RiseSetLCT
rs of
RiseSetLCT
Circumpolar -> RiseSetLCT
forall a. RiseSet a
Circumpolar
RiseSetLCT
NeverRises -> RiseSetLCT
forall a. RiseSet a
NeverRises
RiseSetLCT
_ -> if RSInfo LocalCivilTime -> RSInfo LocalCivilTime -> Bool
isOK RSInfo LocalCivilTime
rsi RSInfo LocalCivilTime
rsi' Bool -> Bool -> Bool
|| Int
iterNo Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
maxIters
then RiseSetLCT
rs
else (RiseSetLCT -> RSInfo LocalCivilTime)
-> RSInfo LocalCivilTime -> Int -> RiseSetLCT
calc RiseSetLCT -> RSInfo LocalCivilTime
getRSInfo RSInfo LocalCivilTime
rsi' (Int
iterNoInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
isOK :: RSInfo LocalCivilTime -> RSInfo LocalCivilTime -> Bool
isOK :: RSInfo LocalCivilTime -> RSInfo LocalCivilTime -> Bool
isOK (LocalCivilTime
t1, DecimalDegrees
_) (LocalCivilTime
t2, DecimalDegrees
_) = (Double -> Double
forall a. Num a => a -> a
abs Double
d) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< (Double
hDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
24)
where JD Double
d = (LocalCivilTime -> JulianDate
lctUniversalTime LocalCivilTime
t1) JulianDate -> JulianDate -> JulianDate
forall a. Num a => a -> a -> a
- (LocalCivilTime -> JulianDate
lctUniversalTime LocalCivilTime
t2)
DH Double
h = DecimalHours
eps
maxIters :: Int
maxIters = Int
3
getRiseTime :: RiseSetLCT -> RSInfo LocalCivilTime
getRiseTime :: RiseSetLCT -> RSInfo LocalCivilTime
getRiseTime (RiseSet RSInfo LocalCivilTime
r RSInfo LocalCivilTime
_) = RSInfo LocalCivilTime
r
getSetTime :: RiseSetLCT -> RSInfo LocalCivilTime
getSetTime :: RiseSetLCT -> RSInfo LocalCivilTime
getSetTime (RiseSet RSInfo LocalCivilTime
_ RSInfo LocalCivilTime
s) = RSInfo LocalCivilTime
s
buildResult :: RiseSet a -> RiseSet a -> RiseSet (Maybe a)
buildResult (RiseSet a
r a
_) (RiseSet a
_ a
s) = Maybe a -> Maybe a -> RiseSet (Maybe a)
forall a. a -> a -> RiseSet a
RiseSet (a -> Maybe a
forall a. a -> Maybe a
Just a
r) (a -> Maybe a
forall a. a -> Maybe a
Just a
s)
buildResult (RiseSet a
r a
_) RiseSet a
_ = Maybe a -> Maybe a -> RiseSet (Maybe a)
forall a. a -> a -> RiseSet a
RiseSet (a -> Maybe a
forall a. a -> Maybe a
Just a
r) Maybe a
forall a. Maybe a
Nothing
buildResult RiseSet a
_ (RiseSet a
_ a
s) = Maybe a -> Maybe a -> RiseSet (Maybe a)
forall a. a -> a -> RiseSet a
RiseSet Maybe a
forall a. Maybe a
Nothing (a -> Maybe a
forall a. a -> Maybe a
Just a
s)
riseAndSetLCT :: GeographicCoordinates
-> LocalCivilDate
-> DecimalDegrees
-> EquatorialCoordinates1
-> RiseSetLCT
riseAndSetLCT :: GeographicCoordinates
-> LocalCivilDate
-> DecimalDegrees
-> EquatorialCoordinates1
-> RiseSetLCT
riseAndSetLCT (GeoC DecimalDegrees
latitude DecimalDegrees
longitude) LocalCivilDate
lcd DecimalDegrees
shift EquatorialCoordinates1
ec
= DecimalDegrees -> LocalCivilDate -> RiseSetLST -> RiseSetLCT
toRiseSetLCT DecimalDegrees
longitude LocalCivilDate
lcd (RiseSetLST -> RiseSetLCT) -> RiseSetLST -> RiseSetLCT
forall a b. (a -> b) -> a -> b
$ EquatorialCoordinates1
-> DecimalDegrees -> DecimalDegrees -> RiseSetLST
riseAndSet EquatorialCoordinates1
ec DecimalDegrees
shift DecimalDegrees
latitude
toRiseSetLCT :: DecimalDegrees
-> LocalCivilDate
-> RiseSetLST
-> RiseSetLCT
toRiseSetLCT :: DecimalDegrees -> LocalCivilDate -> RiseSetLST -> RiseSetLCT
toRiseSetLCT DecimalDegrees
longitude LocalCivilDate
lcd (RiseSet (LocalSiderealTime
rise, DecimalDegrees
azRise) (LocalSiderealTime
set, DecimalDegrees
azSet)) =
let toLCT :: LocalSiderealTime -> LocalCivilTime
toLCT LocalSiderealTime
lst = DecimalDegrees
-> LocalCivilDate -> LocalSiderealTime -> LocalCivilTime
lstToLCT DecimalDegrees
longitude LocalCivilDate
lcd LocalSiderealTime
lst
rise' :: LocalCivilTime
rise' = LocalSiderealTime -> LocalCivilTime
toLCT LocalSiderealTime
rise
set' :: LocalCivilTime
set' = LocalSiderealTime -> LocalCivilTime
toLCT LocalSiderealTime
set
in RSInfo LocalCivilTime -> RSInfo LocalCivilTime -> RiseSetLCT
forall a. a -> a -> RiseSet a
RiseSet (LocalCivilTime
rise', DecimalDegrees
azRise) (LocalCivilTime
set', DecimalDegrees
azSet)
toRiseSetLCT DecimalDegrees
_ LocalCivilDate
_ RiseSetLST
Circumpolar = RiseSetLCT
forall a. RiseSet a
Circumpolar
toRiseSetLCT DecimalDegrees
_ LocalCivilDate
_ RiseSetLST
NeverRises = RiseSetLCT
forall a. RiseSet a
NeverRises
dhToJD :: DecimalHours -> JulianDate -> JulianDate
dhToJD :: DecimalHours -> JulianDate -> JulianDate
dhToJD (DH Double
hours) JulianDate
day = JulianDate
day JulianDate -> JulianDate -> JulianDate
forall a. Num a => a -> a -> a
+ (Double -> JulianDate
JD (Double -> JulianDate) -> Double -> JulianDate
forall a b. (a -> b) -> a -> b
$ Double
hoursDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
24)