module Data.Astro.Planet.PlanetMechanics
(
planetMeanAnomaly
, planetTrueAnomaly1
, planetTrueAnomaly2
, planetHeliocentricRadiusVector
, planetHeliocentricLongitude
, planetHeliocentricLatitude
, planetProjectedRadiusVector
, planetProjectedLongitude
, planetEclipticLongitude
, planetEclipticLatitude
, planetPosition
, planetPosition1
, planetDistance
, planetDistance1
, planetAngularDiameter
, planetPhase1
, planetPertubations
, planetBrightLimbPositionAngle
)
where
import qualified Data.Astro.Utils as U
import Data.Astro.Types (DecimalDegrees(..), AstronomicalUnits(..), toRadians, fromRadians, fromDecimalHours)
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)
reduceDegrees :: DecimalDegrees -> DecimalDegrees
reduceDegrees = U.reduceToZeroRange 360
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)
planetTrueAnomaly1 pd jd =
let meanAnomaly = toRadians $ planetMeanAnomaly pd jd
e = pdE pd
in reduceDegrees $ fromRadians $ meanAnomaly + 2*e*(sin meanAnomaly)
planetTrueAnomaly2 pd jd =
let meanAnomaly = toRadians $ planetMeanAnomaly pd jd
e = pdE pd
eccentricAnomaly = solveKeplerEquation e meanAnomaly 0.000000001
trueAnomaly = 2 * atan (sqrt((1 + e) / (1 - e)) * tan (eccentricAnomaly / 2))
in reduceDegrees $ fromRadians trueAnomaly
planetHeliocentricLongitude :: PlanetDetails -> DecimalDegrees -> DecimalDegrees
planetHeliocentricLongitude pd trueAnomaly = reduceDegrees $ (pdOmegaBar pd) + trueAnomaly
planetHeliocentricLatitude :: PlanetDetails -> DecimalDegrees -> DecimalDegrees
planetHeliocentricLatitude pd hcl =
let l' = toRadians hcl
i' = toRadians $ pdI pd
bigOmega' = toRadians $ pdBigOmega pd
in fromRadians $ asin $ (sin $ l' - bigOmega')*(sin i')
planetHeliocentricRadiusVector :: PlanetDetails -> DecimalDegrees -> AstronomicalUnits
planetHeliocentricRadiusVector pd trueAnomaly =
let nu = toRadians trueAnomaly
AU alpha = pdAlpha pd
e = pdE pd
in AU $ alpha*(1 - e*e)/(1+e*(cos nu))
planetProjectedLongitude :: PlanetDetails -> DecimalDegrees -> DecimalDegrees
planetProjectedLongitude pd hcl =
let hcl' = toRadians hcl
bigOmega = pdBigOmega pd
bigOmega' = toRadians $ bigOmega
i' = toRadians $ pdI pd
y = (sin $ hcl'-bigOmega')*(cos i')
x = (cos $ hcl'-bigOmega')
n = fromRadians $ atan2 y x
in n + bigOmega
planetProjectedRadiusVector :: PlanetDetails -> DecimalDegrees -> AstronomicalUnits -> AstronomicalUnits
planetProjectedRadiusVector pd psi (AU hcr) = AU $ hcr*cos(toRadians psi)
outerPlanetEclipticLongitude :: DecimalDegrees -> AstronomicalUnits -> DecimalDegrees -> AstronomicalUnits -> DecimalDegrees
outerPlanetEclipticLongitude lp (AU rp) le (AU re) =
let lp' = toRadians lp
le' = toRadians le
x = atan $ re * (sin $ lp'-le')/(rp - re*(cos $ lp'-le'))
in reduceDegrees $ (fromRadians x) + lp
innerPlanetEclipticLongitude :: DecimalDegrees -> AstronomicalUnits -> DecimalDegrees -> AstronomicalUnits -> DecimalDegrees
innerPlanetEclipticLongitude lp (AU rp) le (AU re) =
let lp' = toRadians lp
le' = toRadians le
x = atan $ rp * (sin $ le'-lp')/(re - rp*(cos $ le'-lp'))
in reduceDegrees $ (fromRadians x) + le + 180
planetEclipticLongitude :: PlanetDetails -> DecimalDegrees -> AstronomicalUnits -> DecimalDegrees -> AstronomicalUnits -> DecimalDegrees
planetEclipticLongitude pd
| isInnerPlanet pd = innerPlanetEclipticLongitude
| otherwise = outerPlanetEclipticLongitude
planetEclipticLatitude :: DecimalDegrees
-> DecimalDegrees
-> AstronomicalUnits
-> DecimalDegrees
-> AstronomicalUnits
-> DecimalDegrees
-> DecimalDegrees
planetEclipticLatitude psi lp (AU rp) le (AU re) lambda =
let psi' = toRadians psi
lp' = toRadians lp
le' = toRadians le
lambda' = toRadians lambda
y = rp*(tan psi')*(sin $ lambda' - lp')
x = re * (sin $ lp' -le')
in fromRadians $ atan (y/x)
planetPosition :: (PlanetDetails -> JulianDate -> DecimalDegrees)
-> PlanetDetails -> PlanetDetails -> JulianDate
-> EquatorialCoordinates1
planetPosition trueAnomaly pd ed jd =
let nup = trueAnomaly pd jd
lp = planetHeliocentricLongitude pd nup
rp = planetHeliocentricRadiusVector pd nup
psi = planetHeliocentricLatitude pd lp
lp' = planetProjectedLongitude pd lp
rp' = planetProjectedRadiusVector pd psi rp
nue = trueAnomaly ed jd
le = planetHeliocentricLongitude ed nue
re = planetHeliocentricRadiusVector ed nue
lambda = planetEclipticLongitude pd lp' rp' le re
beta = planetEclipticLatitude psi lp' rp' le re lambda
ec = eclipticToEquatorial (EcC beta lambda) jd
in ec
planetDistance :: (PlanetDetails -> JulianDate -> DecimalDegrees)
-> PlanetDetails -> PlanetDetails -> JulianDate
-> AstronomicalUnits
planetDistance trueAnomaly pd ed jd =
let nup = trueAnomaly pd jd
lp = planetHeliocentricLongitude pd nup
AU rp = planetHeliocentricRadiusVector pd nup
psi = planetHeliocentricLatitude pd lp
nue = trueAnomaly ed jd
le = planetHeliocentricLongitude ed nue
AU re = planetHeliocentricRadiusVector ed nue
ro = sqrt $ re*re + rp*rp - 2*re*rp*(cos . toRadians $ lp - le)*(cos $ toRadians psi)
in AU ro
planetAngularDiameter :: PlanetDetails -> AstronomicalUnits -> DecimalDegrees
planetAngularDiameter pd (AU ro) = (pdBigTheta pd)/(DD ro)
planetPhase1 :: PlanetDetails -> PlanetDetails -> JulianDate -> Double
planetPhase1 pd ed jd =
let nup = planetTrueAnomaly1 pd jd
lp = planetHeliocentricLongitude pd nup
rp = planetHeliocentricRadiusVector pd nup
psi = planetHeliocentricLatitude pd lp
lp' = planetProjectedLongitude pd lp
rp' = planetProjectedRadiusVector pd psi rp
nue = planetTrueAnomaly1 ed jd
le = planetHeliocentricLongitude ed nue
re = planetHeliocentricRadiusVector ed nue
lambda = planetEclipticLongitude pd lp' rp' le re
d = toRadians $ lambda - lp
in (1+ (cos d)) * 0.5
planetPosition1 :: PlanetDetails -> PlanetDetails -> JulianDate
-> EquatorialCoordinates1
planetPosition1 = planetPosition planetTrueAnomaly1
planetDistance1 :: PlanetDetails -> PlanetDetails -> JulianDate
-> AstronomicalUnits
planetDistance1 = planetDistance planetTrueAnomaly1
planetPertubations :: Planet -> JulianDate -> DecimalDegrees
planetPertubations Jupiter jd =
let (a, _, v, _) = pertubationsQuantities jd
v' = toRadians v
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
q' = toRadians q
v' = toRadians v
b' = toRadians b
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
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)
planetBrightLimbPositionAngle :: EquatorialCoordinates1 -> EquatorialCoordinates1 -> DecimalDegrees
planetBrightLimbPositionAngle (EC1 deltaP alphaP) (EC1 deltaS alphaS) =
let dAlpha = toRadians $ fromDecimalHours $ alphaS - alphaP
deltaP' = toRadians deltaP
deltaS' = toRadians deltaS
y = (cos deltaS')*(sin dAlpha)
x = (cos deltaP')*(sin deltaS') - (sin deltaP')*(cos deltaS')*(cos dAlpha)
chi = reduceDegrees $ fromRadians $ atan2 y x
in chi