# plot_epsilon_ring_example.py
#
# Illustrate use of uranus_ring_frames *tf frame kernel
# to plot epsilon ring as seen from Earth

import matplotlib.pyplot as plt
import numpy as np
import spiceypy as spice

# !!!  Modify this line to point to your local kernels directory !!!

kernels_dir = '../../../../kernels/'

kernels	= ['ura111.bsp',
	'uranus_ringframes_rfrench20201201_v1.tf',
	'naif0012.tls']

kernels_list=[kernels_dir + k for k in kernels]
spice.furnsh(kernels_list)

UTCstr	= '2020 Jan 25 12:00:00'
ETsec = spice.str2et(UTCstr)

Re	= 25559.0 # Equatorial radius of Uranus

# construct a transformation from Earth equatorial to sky plane frames

state,ltime = spice.spkezr('Uranus',ETsec,'J2000','CN','Earth')
dist,ra,dec = spice.recrad(state[0:3])

mm_temp1    = spice.rotate(dec,2)
mm_temp2    = spice.rotate(-ra,3)
mm_SKY2EEQ  = spice.mxm(mm_temp2,mm_temp1)
mm_UR2EEQ   = spice.pxform('URANUS_EQUATORIAL','J2000',ETsec)

# Uranus pole direction in J2000 coordinates

nhat_pole   = mm_UR2EEQ[0:3,2]
nhat_pole  *= Re # scale by radius of planet
nhat_pole   = np.array(nhat_pole)

# compute sky plane coordinates of pole

pole_sky    = spice.mtxv(mm_SKY2EEQ,nhat_pole)

# compute epsilon ring sky plane coordinates

a_km_ring	= 51149.
ae_km_ring	= 406.
e_ring		= ae_km_ring/a_km_ring
ntheta          = 360*50
theta           = np.radians(np.arange(ntheta)/50.)
rvals		= a_km_ring * (1.0 - e_ring**2)/(1.0 + e_ring * np.cos(theta))
xvals		= rvals * np.cos(theta)
yvals		= rvals * np.sin(theta)
mm_RPL2EEQ      = spice.pxform('URING_EPSILON','J2000',ETsec-ltime)
mm_RPL2SKY      = spice.mtxm(mm_SKY2EEQ,mm_RPL2EEQ)

ntheta		
ff_km 		= np.zeros(ntheta)
gg_km 		= np.zeros(ntheta)
zvals 		= np.zeros(ntheta)
for nn in range(ntheta):
        vec_ring  = np.array([xvals[nn],yvals[nn],0])
        vec_sky   = spice.mxv(mm_RPL2SKY,vec_ring)
       	ff_km[nn] = vec_sky[1]
       	gg_km[nn] = vec_sky[2]

# plot the results

plt.xlabel('East (km)')
plt.ylabel('North (km)')
plt.title('Epsilon ring from Earth '+UTCstr)
plt.axis('equal')
plt.axis([3*Re,-3*Re,-3*Re,3*Re])

# draw outline of circular planet on sky

xx      = np.cos(theta)                                               
yy      = np.sin(theta)                                               
x       = Re * xx
y       = Re * yy

plt.plot(x,y)

# mark the visible pole location

if pole_sky[0] < 0:
	plt.plot(pole_sky[1],pole_sky[2],"o")
else:
	plt.plot(-pole_sky[1],-pole_sky[2],"o")

# draw epsilon ring and periapse location

plt.plot(ff_km,gg_km)
plt.plot(ff_km[0],gg_km[0],marker='+',markersize=10)

plt.tight_layout()
pdf_figure = 'plot_epsilon_ring_example_python.pdf'
plt.savefig(pdf_figure)
print('Saved plot as '+pdf_figure)
plt.show()