# Copyright 2014 - 2015 Anna-Lea Lesage for MASCARA
#
# This file is part of the mascara package.
# mascara is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# mascara is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with mascara. If not, see <http://www.gnu.org/licenses/>.
'''
.. module: coordinates
.. Author: Anna-Lea Lesage
The coordinates module defines a set of function which are exclusively related
to coordinates conversions and coordinates updates.
:mod:`mascara.observer.coordinates` contains the functions to compute nutation and
precession parameters. Furthermore, it also contains the functions to compute many
planetary elements (mean longitude, mean anomaly ...) for the Sun, the Earth, the
Moon and some other planets.
'''
__version__ = '15.10.07'
### *** Latest Changes Record *** ###
### 17.11.2014 Corrected the sun coordinates function to deliver correct altitudes:
### Added the JD1900 entry.
### 08.07.2015 Updated the SunAlt function to use now JD2000 instead of JD1900
### 07.10.2015 Moved the functions greenwich_mean_siderial_time, and Earth_rotation_angle
### from Masfunct to here.
### ******************************* ###
import numpy as np
from mascara.constants import deg2as, ciras, as2rad, JD2000, JDyr
[docs]def Nutation(jd): # TESTED, WORKING ... but carefull when comparing with other results...
''' Nutation calculates for a given date the mean obliquity, the nutation
longitude and the nutation obliquity according to the IAU 2000B model.
The Delaunay coefficients are from Simon and al. (1994).
The model is simplified to the coefficients higher than 0.05".
Obliquity is calculated for J2000.
Input:
- scalar julian date (only one date per call)
Outputs :
- Obl, the obliquity of the ecliptic
- psi, the nutation in the longitude of the ecliptic
- eps, the nutation in the obliquity of the ecliptic
All outputs are in RADIAN
Comments:
- The Delaunay angles are extracted from the paper of Simon and al (1994)
- The nutation angles are calculated according to the formulas of the
IAU2000B model which can be found in the paper of McCarthy and Luzum (2002)
'''
from mascara.constants import NutationCoeff as NC # Nutation coefficients are tested and OK
from numpy.polynomial.polynomial import polyval
t = (jd - JD2000)/JDyr
# Calculating the mean obliquity for the epoch 2000, using the coefficients
# from N.Capitaine and al (2003), IAU2000 Precession Table 8, PO3.
epsa = (84381.406, -46.836769, -0.0001831, 0.00200340, -0.000000576)
Obl = (polyval(t, epsa))*as2rad
# Defining the Delaunay angles (from Simon and al. 1994)
# Mean Anomaly of the Moon (l)
maM = MoonMeanAnomaly(t)
# Mean anomaly of the Sun (l')
maS = SunMeanAnomaly(t)
# Mean elongation of the Moon (D)
meM = MoonMeanElong(t)[0]
# Mean argument of the latitude of the Moon (F)
mlM = MoonArgLatitude(t)
# Longitude of the ascending node of the Moon (Omega)
Om = MoonAscendingNode(t)
args = np.outer(NC.ll, maM) + np.outer(NC.lpl, maS) + np.outer(NC.Dl, meM) + \
np.outer(NC.Fl, mlM) + np.outer(NC.oml, Om)
sarg = np.sin(args)
carg = np.cos(args)
uas2rad = 1e-7*as2rad
mas2rad = 1e-3*as2rad
sls = np.meshgrid(t, NC.sls)[1]
cls = np.meshgrid(t, NC.cls)[1]
csls = np.meshgrid(t, NC.csls)[1]
ccls = np.meshgrid(t, NC.ccls)[1]
DP = np.sum((sls + np.outer(NC.dsls,t))*sarg + csls*carg, axis=0) * uas2rad
DE = np.sum((cls + np.outer(NC.dcls,t))*carg + ccls*sarg, axis=0) * uas2rad
CDP = -0.135*mas2rad
CDE = 0.388*mas2rad
Dpsi = DP + CDP
Deps = DE + CDE
return Obl, Dpsi, Deps # in radians-> to be used with rotation matrix rotmat
[docs]def precession_matrix(jd): # TESTED WORKING very fine!
''' Precession calculates the RA and Dec coordinates changes from a reference
epoch to a new epoch.
Inputs:
- the date, in Julian Date format
Output:
- the precession rotation matrix
Comments: the formulas used here respected the IAU 2000 resolution.
They were extracted from the paper by Capitaine and al (2003),
"Expressions for IAU 2000 precession quantities"
'''
from numpy.polynomial.polynomial import polyval
from mascara.funcs.utils import rotmat
# Convert the date to Julian century:
tt = (jd - JD2000)/JDyr
# The precession matrix coefficients extracted from
# Capitaine and al (2003) "Precession quantities" :
pzta = [2.650545, 2306.083227, 0.2988499, 0.01801828, 0.000005971, 0.0000003173]
pza = [-2.650545, 2306.077181, 1.0927348, 0.01826837, -0.000028596, -0.0000002904]
ptheta = [0, 2004.191903, -0.4294934, -0.04182264, -0.000007089, -0.0000001274]
zta = polyval(tt,pzta)/3600
za = polyval(tt, pza)/3600
theta = polyval(tt, ptheta)/3600
Pmatrix = rotmat(-za,'z',degree=True)*rotmat(theta,'y',degree=True)*rotmat(-zta,'z',degree=True)
return Pmatrix
[docs]def nutation_matrix(*args):
''' Nutation_matrix returns the rotation matrix required for the correction
of the Earth's nutation effects. The angles should have been calculated via
Nutation before input. The matrix is calculated according to the IAU 2000
standard, using the formulas from the USNO circular n179.
Inputs:
- a date (in jd or datetime format) OR
- obli, the obliquity of the ecliptic
- dpsi, the nutation longitude
- deps, the nutation in obliquity
Keywords: ?
Output : Nmatrix, the nutation rotation matrix.
'''
from datetime import datetime
from mascara.funcs.utils import rotmat
from mascara.funcs.timing import calendar2jd
# If only one argument is supplied then it is probably a date,
# but it could be a datetime object....
if type(args) == datetime:
jd = calendar2jd(args)
obli, dpsi, deps = Nutation(jd)
if len(args)==1:
jd = np.float64(args)
obli,dpsi,deps = Nutation(jd)
if len(args) == 3:
obli, dpsi, deps = args
Nmatrix = rotmat(-(obli+deps),'x')*rotmat(-dpsi,'z')*rotmat(obli,'x')
return Nmatrix
[docs]def bias_matrix(FK5='False'):
''' A Bias matrix is required to convert from ICRS data to the dynamical mean
equation and equinox J2000.
Keywords:
- if FK5 is true , then the bias matrix is calculated for ICRS to
FK5 system, else it generates a general matrix
Output:
- The Bias rotation matrix between the ICRS and the J2000 frame, or
the bias matrix between the GCRS and the CIO.
'''
from mascara.funcs.utils import rotmat
from mascara.constants import as2rad
if FK5:
da0 = -22.9e-3 * as2rad
et0 = 9.1e-3 * as2rad
nt0 = 19.9e-3 * as2rad
else:
da0 = -14.6e-3 * as2rad
et0 = -16.6170e-3 * as2rad
nt0 = -6.8192e-3 * as2rad
Bmatrix = rotmat(-nt0, 'x')*rotmat(et0, 'y')*rotmat(da0,'z')
return Bmatrix
def Aberration(ra, dec, jd):
''' Calculates the changes in ra and dec due to annual aberration
Inputs:
- ra, dec the equatorial coordinates
- date, the date in Julian Date format
Outputs:
- sra, sdec, the corrected coordinates.
'''
from mascara.constants import aber, deg2rad, rad2as
jt = (jd - JD2000)/JDyr
## Calculates the Sun's coordinates at the given date.
sunra, sundec, sunlon = SunCoord(jd, longitude=True)
## Calculate obliquity of the ecliptic:
obli, dpsi, deps = Nutation(jd)
eps0 = (23*deg2as+26.*60+21.448)-46.815*jt-0.00059*jt**2 +0.001813*jt**3
eps = (eps0+deps*rad2as)/3600*deg2rad
## Earth's orbital excentricity:
Ee = EarthExcentricity(jt)
## Longitude of Perihelion:
epi = EarthLongPerihelion(jt)
cd = np.cos(dec * deg2rad)
sd = np.sin(dec * deg2rad)
ce = np.cos(eps)
te = np.tan(eps)
cp = np.cos(epi)
sp = np.sin(epi)
cs = np.cos(sunlon*deg2rad)
ss = np.sin(sunlon*deg2rad)
ca = np.cos(ra*deg2rad)
sa = np.sin(ra*deg2rad)
term1 = (ca*cs*ce+sa*ss)/cd
term2 = (ca*cp*ce+sa*sp)/cd
term3 = (cs*ce*(te*cd-sa*sd)+ca*sd*ss)
term4 = (cp*ce*(te*cd-sa*sd)+ca*sd*sp)
sra = - aber*term1+Ee*aber*term2
sdec = -aber *term3 + Ee*aber*term4
return sra/3600, sdec/3600
[docs]def greenwich_mean_sidereal_time(date, apparent=True): #TESTED WORKING
''' Calculates the GMST at the current observation date.
This is the first step towards the calculation of the Local Mean
Sidereal Time (LMST). The coefficients are extracted from the circular
179 of USNO.
Input:
- date, the date in Julian days, or in datetime format
Keyword:
- Apparent, if set, the program will calculate the apparement sidereal
time at the given date, taking into account the Earth's nutation.
Output:
- the GMST in fractional hours
'''
from numpy.polynomial.polynomial import polyval
from mascara.constants import as2rad
from mascara.funcs.timing import calendar2jd
from datetime import datetime
# Check what is the input format, and convert the date to JD
if type(date)==datetime:
jd = calendar2jd(date, UTC=True)
elif type(date) == np.float64:
jd = date
else:
raise SyntaxError('The date format is not recognized')
t = (jd - JD2000)/JDyr
pgmst = np.array([0.014506, 4612.156534, 1.3915817,-0.00000044, -0.000029956, -0.0000000368])
pgmst = pgmst*as2rad
era = earth_rotation_angle(jd, degree=False)
gmst = era + polyval(t, pgmst)
if apparent:
eps0, dpsi, deps = Nutation(jd)
gmst = gmst + dpsi*np.cos(eps0)
return (gmst*12/ np.pi)%24
def earth_rotation_angle(date, UTC=True, degree= None): # TESTED WORKING
''' Calculates the Earth Rotation Angle (ERA) at the input date.
The coefficients and formula are from the circular 179 of the USNO.
Input:
- date, the date chosen to calculate the ERA. Either a datetime object or JD date
Output:
- era the earth rotation angle in term of rotation units (rad or degree)
'''
from mascara.funcs.timing import calendar2jd
from datetime import datetime
if type(date)== datetime:
jd = calendar2jd(date, UTC=UTC)
else: jd=date
Du = jd - JD2000
era = (0.7790572732640 + 0.00273781191135448 * Du + (Du%1.0)) %1.0
if degree == None:
return era
if degree == True:
return era/360
if degree == False:
return era*2*np.pi
[docs]def EarthExcentricity(jdc):
''' Calculates the Earth's orbital excentricity at the given date;
Input:
- t, the date in Julian format
Ouput:
- the Earth's Excentricity (no dimension)
'''
if (jdc > 2e6).all():
jdc = (jdc - JD2000)/JDyr
ee = 0.0167086342 - 0.00004203654*jdc - 0.000000126734*jdc**2
return ee
def EarthLongPerihelion(jdc):
''' Calculates the longitude of perihelion for the Earth at the given date
Input: jdc the date in Julian format
Output: elp in radian
'''
if (jdc > 2e6).all():
jdc = (jdc - JD2000)/JDyr
elp = (102.93734808*deg2as + 11612.35290*jdc + 53.27577*jdc**2)%ciras * as2rad
return elp
def EarthMeanLong(jdc):
''' Calculates the Mean Longitude of the Earth
Input: jdc the date in Julian format
Output: Earth Mean Longitude in radian
'''
if (jdc > 2e6).all():
jdc = (jdc - JD2000)/JDyr
eml = (100.46645683*deg2as + 129597742.283429*jdc - 2.04411*jdc**2)%ciras * as2rad
return eml
def EarthMeanAnomaly(jdc):
''' Forms the mean anomaly term for the Earth at the given Julian date t.
Input:
- t, the date in Julian days
Outputs:
- ema, the Earth's mean anomaly in radian
- ecorr, the ellipticity of the orbit in radian!!
**Comments**: the correction term is calculated from the equation of centre in
the Moon-Earth system
'''
if (jdc > 2e6).all():
jdc = (jdc - JD2000)/JDyr
return SunMeanAnomaly(jdc)
def MoonDistance(jdc, perturbate = True):
''' Calculates the Moon-Earth distance (in km!!!) at the given date.
Input:
- jd, the date in Julian format
Keyword:
- perturbate, if True include perturbations terms in the computatino
Output:
- d, the Moon-Earth distance in km!!!
'''
if (jdc > 2e6).all():
jdc = (jdc - JD2000)/JDyr
dm = 383397.7725 + 0.0040*jdc
if perturbate:
D = MoonMeanElong(jdc)[0]
l = MoonMeanAnomaly(jdc)
lp = SunMeanAnomaly(jdc)
pert = 3400.4*np.cos(2*D) - 635.6*np.cos(2*D-l) - 235.6*np.cos(l) + \
218.1*np.cos(2*D-lp) + 181.0*np.cos(2*D+l)
dm = dm + pert
return dm
def MoonLongPerihelion(jdc, perturbate =True):
''' Calculates the longitude of perihelion for the Moon at the given date
Input:
-jd, the date in Julian format
Keyword:
- perturbate, to compute the additional perturbation terms of the Moon's
excentricity
'''
if (jdc > 2e6).all():
jdc = (jdc - JD2000)/JDyr
mlp = (83.35324312*deg2as + 14643420.2669*jdc - 38.2702*jdc**2)%ciras * as2rad
if perturbate:
D = MoonMeanElong(jdc)[0]
l = MoonMeanAnomaly(jdc)
lp = SunMeanAnomaly(jdc)
ad1 = (-55609*np.sin(2*D-l) - 34711*np.sin(2*D-2*l) - 9792*np.sin(l))*as2rad
ad2 = (9385*np.sin(4*D-3*l) + 7505*np.sin(4*D-2*l) + 5318*np.sin(2*D+l))*as2rad
ad3 = (3484*np.sin(4*D-4*l) - 3417*np.sin(2*D-lp-l) - 2530*np.sin(6*D-4*l))*as2rad
ad4 = (-2376*np.sin(2*D) - 2075*np.sin(2*D-3*l) - 1883*np.sin(2*l))*as2rad
ad5 = (-1736*np.sin(6*D-5*l) + 1626*np.sin(lp) - 1370*np.sin(6*D-3*l))*as2rad
mlp = mlp + ad1 + ad2 + ad3 + ad4 + ad5
return mlp
def MoonExcentricity(jdc, perturbate=False):
''' Computes the excentricity of the Moon at the given date
Input:
-jd, the date in Julian format
Keyword:
- perturbate, to compute the additional perturbation terms of the Moon's
excentricity
Output:
- Mex, the Moon's excentricity
'''
if (jdc > 2e6).all():
jdc = (jdc - JD2000)/JDyr
mex = 0.055545526-0.000000016*jdc
if perturbate:
D = MoonMeanElong(jdc)[0]
l = MoonMeanAnomaly(jdc)
lp = SunMeanAnomaly(jdc)
F = MoonArgLatitude(jdc)
ad1 = (0.014216*np.cos(2*D-l) + 0.008551*np.cos(2*D - 2*l) - 0.001383*np.cos(l))
ad2 = (0.001356*np.cos(2*D+l) - 0.001147*np.cos(4*D-3*l) - 0.000914*np.cos(4*D - 2*l))
ad3 = (0.000869*np.cos(2*D-lp-l)-0.000627*np.cos(2*D) - 0.000394*np.cos(4*D-4*l))
ad4 = (0.000282*np.cos(2*D-lp-2*l)-0.000279*np.cos(D-l) - 0.000236*np.cos(2*l))
ad5 = (0.000231*np.cos(4*D) + 0.000229*np.cos(6*D-4*l) - 0.000201*np.cos(2*l-2*F))
mex = mex + ad1 + ad2 + ad3 + ad4 + ad5
return mex
def MoonInclination(jdc, perturbate=True):
''' computes the inclination angle of the Moon at the given date
'''
if (jdc > 2e6).all():
jdc = (jdc - JD2000)/JDyr
mincl = (5.15668983*deg2as - 0.00008*jdc)*as2rad
if perturbate:
D = MoonMeanElong(jdc)[0]
l = MoonMeanAnomaly(jdc)
lp = SunMeanAnomaly(jdc)
F = MoonArgLatitude(jdc)
pert = (486.26*np.cos(2*D-2*F) - 40.13*np.cos(2*D) + 37.51*np.cos(2*F) + \
25.73*np.cos(2*l - 2*F) + 19.97*np.cos(2*D-lp-2*F))*as2rad
mincl = mincl+pert
return mincl
def MoonArgLatitude(jdc): ###WORKING
''' Forms the Moon's argument of latitude at the given Julian date t.
Input:
- jdc, the Julian date (or in format (JD - JD2000)/Jdyr )
Output:
- the Mean argument of latitude of the Moon
Comment: the formula is from the Simon and al (1994) paper
" Precession formulae and mean elements for the Moon and the planets"
'''
if (jdc > 2e6).all():
jdc = (jdc - JD2000)/JDyr
mal = (93.27209062*deg2as + 1739527262.8478*jdc -12.7512*jdc**2)%ciras * as2rad
return mal
def MoonAscendingNode(jdc, perturbate = True): ###WORKING
''' Calculates the longitude of the ascending node of the Moon at the given date.
Input:
- jdc, the Julian date perturbat
Keyword:
- perturbate, if True, include the perturbation terms
Ouput:
- the ascending node Omega
Comment: From Simon and al (1994)
'''
if (jdc > 2e6).all():
jdc = (jdc - JD2000)/JDyr
Om = (125.04455501*deg2as - 6962890.5431*jdc + 7.4722*jdc**2)%ciras * as2rad
if perturbate:
D = MoonMeanElong(jdc)[0]
l = MoonMeanAnomaly(jdc)
lp = SunMeanAnomaly(jdc)
F = MoonArgLatitude(jdc)
per1 = (-5392*np.sin(2*D-2*l) - 540*np.sin(lp) - 441*np.sin(2*D)+423*np.sin(2*F) - \
288*np.sin(2*l-2*F))*as2rad
Om = Om + per1
return Om
[docs]def MoonMeanAnomaly(jdc): ### WORKING
''' Forms the mean anomaly term for the Moon at the given Julian date t.
Input:
- t, the date dafined as (Jd - JD2000)/JDyr
Ouput:
- mma, the Moon's mean anomaly
Comment: The formula is extracted from Simon and al (1994).
'''
if (jdc > 2e6).all():
jdc = (jdc - JD2000)/JDyr
mam = (134.96340251*deg2as + 1717915923.2178*jdc + 31.8792*jdc**2)%ciras * as2rad
return mam
[docs]def MoonMeanLong(jdc, perturbate = True):
''' Defines the Mean longitude of the Moon referred to the mean equinox of the date
Input:
- t, the date expressed in Julian days
Keyword:
- perturbate, if True then include perturbation terms to the calculation of the
Moon's Mean Longitude
Outputs:
- mml, the Moon's mean longitude
'''
if (jdc > 2e6).all():
jdc = (jdc - JD2000)/JDyr
mml = (218.31664563*deg2as + 1732564372.30470*jdc - 5.2790*jdc**2)%ciras *as2rad
if perturbate:
print 'perturabting'
D = MoonMeanElong(jdc)[0]
l = MoonMeanAnomaly(jdc)
lp = SunMeanAnomaly(jdc)
pert = (-3332.9*np.sin(2*D) + 1197.4*np.sin(2*D-l) - 662.5*np.sin(lp) + \
396.3*np.sin(l) - 218.0*np.sin(2*D-lp))*as2rad
mml = mml + pert
return mml
[docs]def MoonMeanElong(jdc): # OK
''' Gives the mean elongation of the Moon from the Sun, and the
corresponding correction term at the given date t.
Input:
- jdc, the date in Julian format (JD - JD2000)/JDyr
Outputs:
- mme, the Moon's mean elongation term
- mocorr, the correction term
Comment: THe formula is from the Simon and al (1994) paper
"Precession formulae and mean elements for the Moon and the planets"
'''
if (jdc > 2e6).all():
jdc = (jdc - JD2000)/JDyr
mme = (297.85019547*deg2as + 1602961601.2090*jdc - 6.33706*jdc**2 + 0.006593*jdc**3)%ciras * as2rad
mocorr = 6.5 * np.sin(mme) * as2rad
return mme, mocorr
def SunMeanLng(jdc):
''' Defines the Sun's mean longitude for a given time t.
Input:
- t, the date expressed in Julian days
Output:
- sml, the Sun's mean longitude
Comment: using the formulas in the Simon and al (1994) paper
'''
from mascara.constants import deg2rad
if (jdc > 2e6).all():
jdc = (jdc - JD2000)/JDyr
sml = ((280.46645 +36000.768925*jdc+0.0003032*jdc**2)%360)*deg2rad
return sml
[docs]def SunMeanAnomaly(jdc):
''' Forms the mean anomaly for the Sun at teh given Julian date t.
Input:
- t, the date defined as (JD - JD2000)/JDyr
Output:
- sma, the Sun's mean anomaly
'''
if (jdc > 2e6).all():
jdc = (jdc - JD2000)/JDyr
sma = (357.52910918*deg2as + 129596581.0481*jdc -0.5532*jdc**2+0.000136*jdc**3)%ciras * as2rad
return sma
[docs]def EquationCenter(jdc, sma):
''' Computes the equation of centre of the Sun at the given date.
Input:
- t the Julian date
- sma the mean Anomaly
Output:
- stl, the equation of centre for the Sun, in RADIAN
'''
if (jdc > 2e6).all():
jdc = (jdc - JD2000)/JDyr
# Sun's equation of center in RADIAN!!!:
ee = EarthExcentricity(jdc)
sec = (2*ee-0.25*ee**3)*np.sin(sma) + 5./4*ee**2*np.sin(2*sma) + 13./12 *ee**3 * np.sin(3*sma)
return sec
def SunTrueLongitude(*args):
''' Calculates the True Longitude of the Sun/Earth at the given date...
Input:
- t, the date in Julian format
- sml, the mean longitude of the Sun
Output:
- the True Longitude of the Sun
'''
### Number of arguments passed:
if len(args) == 1:
jdc = args[0]
if (jdc > 2e6).all():
jdc = (jdc - JD2000)/JDyr
sml = SunMeanLng(jdc)
elif len(args) == 2:
jdc, sml = args
if (jdc > 2e6).all():
jdc = (jdc - JD2000)/JDyr
### Get the excentric anomaly
sec = EquationCenter(jdc)
stl = sml + sec
return stl
def SunTrueAnomaly(*args):
''' Calculates the True Anomaly of the Sun/Earth at the given date...
Input:
- t, the date in Julian format
- sml, the mean longitude of the Sun
Output:
- the True Anomaly of the Sun
'''
### Number of arguments passed:
if len(args) == 1:
jdc = args[0]
if (jdc > 2e6).all():
jdc = (jdc - JD2000)/JDyr
sma = SunMeanAnomaly(jdc)
elif len(args) == 2:
jdc, sml = args
if (jdc > 2e6).all():
jdc = (jdc - JD2000)/JDyr
### Get the excentric anomaly
sec = EquationCenter(jdc)
sta = sma + sec
return sta
def VenusMeanLong(jdc):
''' Forms the mean Longitude term for Venus at the given Julian date jdc
Input:
- jdc, in Julian days
Output:
- vml, Venus' mean longitude
'''
if (jdc > 2e6).all():
jdc = (jdc - JD2000)/JDyr
vml = (181.97980085*deg2as + 2106641364*jdc + 0.59381*jdc**2)%ciras *as2rad
return vml
def VenusPerihelion(jdc):
''' Compute the longitude of the perihelion of Venus for the given date
'''
if (jdc > 2e6).all():
jdc = (jdc - JD2000)/JDyr
vpe = (131.56370300*deg2as + 175.48640*jdc - 498.48184*jdc**2 - \
20.50042*jdc**3 - 0.72432*jdc**4)%ciras * as2rad
return vpe
def VenusMeanAnomaly(jdc):
''' Forms the mean anomaly term for Venus at the given julian date t.
Input:
- jdc, in Julian days (JD - JD1900)/JDyr
Output:
- vma, venus's mean anomaly term in radian
- vcorr, the correction term for Venus's perturbations in radian
'''
from mascara.constants import deg2rad
if (jdc > 2e6).all():
jdc = (jdc - JD2000)/JDyr
vma = VenusMeanLong(jdc) - VenusPerihelion(jdc)
sma = SunMeanAnomaly(jdc)
vcorr = 4.8*np.cos(299.1017*deg2rad + 1.0*vma - 1.0*sma) + \
5.5*np.cos(148.3133*deg2rad + 2.0*vma - 2.0*sma) + \
2.5*np.cos(315.9433*deg2rad + 2.0*vma - 3.0*sma) + \
1.6*np.cos(345.2533*deg2rad + 3.0*vma - 4.0*sma) + \
1.0*np.cos(318.1500*deg2rad + 3.0*vma - 5.0*sma)
return vma, vcorr *as2rad
def JupiterPerihelion(jdc):
''' Computes the Longitude of perihelion for Jupiter at the given time
'''
if (jdc > 2e6).all():
jdc = (jdc - JD2000)/JDyr
jpe = (14.33120687*deg2as + 7758.75163*jdc + 259.95938*jdc**2 - 16.14731*jdc**3)%ciras *as2rad
return jpe
def JupiterMeanLong(jdc):
''' Computes the Mean Longitude of Jupiter at the given Julian Date
'''
if (jdc > 2e6).all():
jdc = (jdc - JD2000)/JDyr
jml = (34.35151874*deg2as + 109256603.77991*jdc - 30.60378*jdc**2)%ciras *as2rad
return jml
def JupiterMeanAnomaly(jdc):
''' Forms the mean anomaly term for Jupiter, and the corresponding perturbations
at the given date.
Input:
- jdc, the date in Julian format (JD - JD1900)/JDyr
Ouputs:
- jma, Jupiter's mean anomaly term
- jcorr, the correction term
'''
from mascara.constants import deg2rad
if (jdc > 2e6).all():
jdc = (jdc - JD2000)/JDyr
jma = JupiterMeanLong(jdc) - JupiterPerihelion(jdc)
ema = EarthMeanAnomaly(jdc)
jcorr = 7.2*np.cos(179.5317*deg2rad - jma + ema) + \
2.6*np.cos(263.2167*deg2rad - jma) + \
2.7*np.cos( 87.1450*deg2rad - 2.0*jma + 2.0*ema) + \
1.6*np.cos(109.4933*deg2rad - 2.0*jma + 1.0*ema)
return jma, jcorr * as2rad
def MarsPerihelion(jdc):
''' Computes the longitude of Perihelion of Mars for the given Julian date
'''
if (jdc > 2e6).all():
jdc = (jdc - JD2000)/JDyr
mpe = (336.06023395*deg2as + 15980.45908*jdc - 62.3280*jdc**2 + 1.86464*jdc**3)%ciras *as2rad
return mpe
def MarsMeanLong(jdc):
''' Computes the mean longitude of Mars at the given Julian date
'''
if (jdc > 2e6).all():
jdc = (jdc - JD2000)/JDyr
maml = (355.43299958*deg2as + 689050774.93988*jdc+0.94264*jdc**2)%ciras * as2rad
return maml
def MarsMeanAnomaly(jdc):
'''Same as above but for Mars...
Input:
- jdc, the date in Julian format (JD - JD1900)/JDyr
'''
from mascara.constants import deg2rad
if (jdc > 2e6).all():
jdc = (jdc - JD2000)/JDyr
#mma = 319.529425 + ((19139.858500*jdc) %360.0)
mma = MarsMeanLong(jdc) - MarsPerihelion(jdc)
ema = EarthMeanAnomaly(jdc)
macorr = 2.0*np.cos(343.8883*deg2rad - 2.0*mma + 2.0*ema) + \
1.8*np.cos(200.4017*deg2rad - 2.0*mma + ema)
return mma, macorr * as2rad
def ComputeExcentricAnomaly(ex, ma):
''' Computes iteratively the Excentric Anomaly for the given excentricity
and mean anomaly.
Inputs:
- ex, the excentricity
- ma, the mean anomaly in RADIAN
Output:
- exa, the excentric anomaly
'''
### Define the ending criterion: tolerance for the iterations:
endtol = 1e-3
### First computation:
E0 = ma + ex*np.sin(ma) * ( 1.+ex*np.cos(ma))
diff = np.abs(E0 - 1.)
while diff > endtol:
E1 = E0 - (E0 - ex * np.sin(E0)-ma)/(1.-ex*np.cos(E0))
diff = np.abs(E1-E0)
E0 = E1
return E1
[docs]def SunCoord(*args, **kwargs):
''' SunCorrLng computes the corrected Sun's Mean Longitude. It takes into account
the perturbation terms of Venus, Mars, Jupiter and the Moon, the mean ellipticity
of the Earth, the long term variations and the abberations.
Input:
- t, in Julian format
Outputs:
- ra, dec, the Sun's coordinates at the given time.
There remain some long term annual/semi-annual variations in the coordinates.
These affect the ra and dec values at the most at the solstice. However the errors
induced are less than one degree.
The code uses the latest formulas for the mean elements of the planets.
'''
from mascara.constants import rad2deg
from datetime import datetime
from mascara.funcs.timing import calendar2jd
from mascara.funcs.utils import pol2rect, rect2pol
try:
t = args[0]
except TypeError:
print 'No input given. Input is: date in Julian Date format'
return
### What is t? datetime or julian date?
if isinstance(t, datetime):
jd = calendar2jd(t)
elif isinstance(t, np.float64):
jd = t
else:
print 'No valid input'
return
jdc = (jd - JD2000)/(JDyr)
longitude = kwargs.get('longitude', False)
### All these quantities are in radiants
sml = np.pi + EarthMeanLong(jdc)
sma = SunMeanAnomaly(jdc)
sec = EquationCenter(jdc, sma)
### Sun's true longitude:
stl = sml + sec
# Longitude term converted to degree:
if longitude:
elong = stl * rad2deg
## Get the true Obliquity
Obl, dpsi, deps = Nutation(t)
### And convert in ra dec
ra = np.arctan2(np.cos(Obl)*np.sin(stl), np.cos(stl))
dec = np.arcsin(np.sin(Obl)*np.sin(stl))
### Some slight correction for Precession?
P = precession_matrix(jd)
B = bias_matrix()
X = pol2rect(ra, dec, degree=False)
Y = P*B*X.T
ra, dec = rect2pol(Y)
if longitude:
return ra, dec, elong
else:
return ra, dec
def MoonCoord(t, longitude=False):
''' MoonCoord computes the Moon's coordinates (right ascension and declination)
at the given date.
Inputs:
- date, in Julian date format or datetime format
Outputs:
- ra, dec, the Moon's celestial coordinates at that date
'''
from mascara.constants import rad2deg
from datetime import datetime
from mascara.funcs.timing import calendar2jd
### What is t? datetime or julian date?
if isinstance(t, datetime):
jd = calendar2jd(t)
elif isinstance(t, np.float64):
jd = t
else:
print 'No valid input'
return
jdc = (jd - JD2000)/JDyr
### Get Delaunay elements of the Moon:
Dm = MoonMeanElong(jdc)[0]
Fm = MoonArgLatitude(jdc)
mphi = (310.17137918*deg2as - 6967051.4360*jdc + 6.2068*jdc**2)*as2rad
### Get Orbital mean elements of the Moon:
mml = MoonMeanLong(jdc, perturbate=True)
mincl = MoonInclination(jdc, perturbate=True)
mex = MoonExcentricity(jdc, perturbate=True)
mperih = MoonLongPerihelion(jdc, perturbate = True)
mom = MoonAscendingNode(jdc, perturbate=True)
### Add complementary terms:
mml = mml + ((2.116*np.sin(mphi)-0.111*np.sin(2*Dm-2*Fm-mphi))*jdc - \
0.0015*np.sin(mphi)*jdc**2)*as2rad
mperih = mperih + ((2.116*np.sin(mphi)-0.111*np.sin(2*Dm-2*Fm-mphi))*jdc - \
0.0015*np.sin(mphi)*jdc**2)*as2rad
mom = mom + ((-520.77*np.sin(mphi)+13.66*np.sin(2*Dm-2*Fm+mphi))*jdc*as2rad)
mincl = mincl + 46.997*np.cos(mphi)*jdc*as2rad
### Compute the Mean Anomaly with the perturbation:
MMA = MoonMeanAnomaly(jdc) #mml - mperih
### Compute the excenttric anomaly:
MEA = ComputeExcentricAnomaly(mex, MMA)
### Compute true anomaly and distance:
MTA = 2* np.arctan(np.sqrt((1-mex)/(1+mex)) * np.tan(MEA/2))
### Compute the heliocentric ecliptic longitude and latitude
hlong = np.arctan(np.tan(mperih - mom + MTA) * np.cos(mincl)) + mom
hlat = np.arctan(np.sin(hlong - mom) * np.tan(mincl))
### get orbit obliquity
obl, npsi, neps = Nutation(jd)
### Convert in Ra and Dec:
mra = np.arctan2(np.sin(hlong)*np.cos(obl) - np.tan(hlat)*np.sin(obl), np.cos(hlong))
mdec = np.arcsin(np.sin(hlat)*np.cos(obl) + np.cos(hlat)*np.sin(obl)*np.sin(hlong))
if longitude:
return mra*rad2deg %360, mdec*rad2deg, hlong, hlat
else:
return mra*rad2deg %360, mdec*rad2deg
