```{-|
Module: Data.Astro.Planet.PlanetMechanics
Description: Planet mechanics

Planet mechanics.
-}

module Data.Astro.Planet.PlanetMechanics
(
planetMeanAnomaly
, planetTrueAnomaly1
, planetHeliocentricLongitude
, planetHeliocentricLatitude
, planetProjectedLongitude
, planetEclipticLongitude
, planetEclipticLatitude
, planetPosition
, planetPosition1
, planetDistance1
, planetAngularDiameter
, planetPhase1
, planetPertubations
, planetBrightLimbPositionAngle
)

where

import qualified Data.Astro.Utils as U
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 =
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
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 =
bigOmega = pdBigOmega pd
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

-- | 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) =
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) =
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 =
y = rp*(tan psi')*(sin \$ lambda' - lp')
x = re * (sin \$ lp' -le')

-- | 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
psi = planetHeliocentricLatitude pd lp
lp' = planetProjectedLongitude pd lp
rp' = planetProjectedRadiusVector pd psi rp
-- earth
nue = trueAnomaly ed jd
le = planetHeliocentricLongitude 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
psi = planetHeliocentricLatitude pd lp
lp' = planetProjectedLongitude pd lp
rp' = planetProjectedRadiusVector pd psi rp
-- earth
nue = planetTrueAnomaly1 ed jd
le = planetHeliocentricLongitude 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
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
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