In [98]:

    

from __future__ import (absolute_import, division, print_function,
                        unicode_literals)
import numpy as np
from astropy.time import Time
from astropy.coordinates import (SkyCoord, ICRS, GCRS, EarthLocation, 
                                 cartesian_to_spherical, get_sun, 
                                 CartesianRepresentation)
import astropy.units as u
import sys
sys.path.insert(0, '/Users/bmorris/git/astroplan')

def get_moon(time, location):
    '''
    Use JPLephem to get approximate ICRS coord for moon at ``time``. Only accepts
    ``times`` between years 1550–2650.
    '''
    # Get position of moon relative to Earth-moon barycenter with JPLephem
    import jplephem
    from jplephem.spk import SPK
    # Note, this table is only good for years 1550–2650. 
    # Table source: 
    #     http://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de430.bsp
    # About: 
    #     http://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/aareadme_de430-de431.txt
    kernel = SPK.open('de430.bsp')
    #t = Time(time)
    # Position of moon relative to Earth-moon barycenter
    cartesian_position = kernel[3,301].compute(t.jd)
    x, y, z = cartesian_position*u.km

    # Convert to GCRS coordinates
    #d, dec, ra = cartesian_to_spherical(*cartesian_position)
    #return GCRS(ra=ra, dec=dec, distance=d, obsgeoloc=location)
    cartrep = CartesianRepresentation(x=x, y=y, z=z)
    return SkyCoord(cartrep, frame=GCRS(obstime=time, obsgeoloc=location))

def calc_moon_phase_angle(moon, sun):
    '''
    Calculate lunar orbital phase [radians]
    '''
    distance_moon_sun = np.sqrt(sun.distance**2 + moon.distance**2 - 
                                2*sun.distance*moon.distance*np.cos(sun.separation(moon)))
    i = np.arccos((distance_moon_sun**2 + moon.distance**2 - 
                   sun.distance**2)/(2*distance_moon_sun*moon.distance))
    return i

def calc_moon_illumination(moon, sun):
    '''
    Calculate fraction of Moon illuminated when moon is at phase angle ``i``
    '''
    i = calc_moon_phase_angle(moon, sun)
    k = (1 + np.cos(i))/2.0
    return k

# Set up observer
lat = '00:00:00'
lon = '00:00:00'
location = EarthLocation.from_geodetic(lon, lat, 0)
t = Time('2015-01-02 00:00:00')

# Get the position of the moon for this observer
moon = get_moon(t, location)
sun = get_sun(t)
print('Moon RA, dec (JPLephem+astropy):', moon.ra.degree, moon.dec.degree)

# Compute alt/az of moon at time with astropy+JPLephem
from astroplan import Observer
obs = Observer(location=location, pressure=0)
altaz = obs.altaz(t, moon)
print('Moon alt, az (JPLephem+astropy):', altaz.alt.degree, altaz.az.degree)
print('Illumination (JPLephem+astropy):', calc_moon_illumination(moon, sun))

# Compute comparisions with PyEphem
import ephem
pyephem_obs = ephem.Observer()
pyephem_obs.lat = lat
pyephem_obs.lon = lon
pyephem_obs.elevation = 0
pyephem_obs.date = t.iso
pyephem_obs.pressure = 0
pyephem_moon = ephem.Moon()
pyephem_moon.compute(pyephem_obs)
print('\nMoon RA, dec (PyEphem):', np.degrees(pyephem_moon.a_ra), np.degrees(pyephem_moon.a_dec))
print('Moon alt, az (PyEphem):', np.degrees(pyephem_moon.alt), np.degrees(pyephem_moon.az))
print('Illumination (PyEphem): ', pyephem_moon.moon_phase)


    

    




Moon RA, dec (JPLephem+astropy): 62.229557107 17.1665467007
Moon alt, az (JPLephem+astropy): [-15.80035594] [ 261.29933774]
Illumination (JPLephem+astropy): [ 0.89074998]

Moon RA, dec (PyEphem): 62.2399250651 17.1678454694
Moon alt, az (PyEphem): 47.4206681054 296.266913816
Illumination (PyEphem):  0.896286911345

In [101]:

    

print(moon, '\n', sun)


    

    




<SkyCoord (GCRS: obstime=2015-01-02 00:00:00.000, obsgeoloc=[ 6378137.        0.        0.] m, obsgeovel=[ 0.  0.  0.] m / s):00:00.000, obsgeoloc=[ 6378137.        0.        0.] m, obsgeovel=[ 0.  0.  0.] m / s): (ra, dec, distance) in (deg, deg, km)
    (62.22955711, 17.1665467, 382245.99875975)> 
 <SkyCoord (GCRS: obstime=2015-01-02 00:00:00.000, obsgeoloc=[ 0.  0.  0.] m, obsgeovel=[ 0.  0.  0.] m / s):00:00.000, obsgeoloc=[ 0.  0.  0.] m, obsgeovel=[ 0.  0.  0.] m / s): (ra, dec, distance) in (deg, deg, AU)
    (283.09239612, -22.8918508, 0.98329387)>

In [ ]: