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

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

kernels_dir = '../../../../kernels/'
kernels	= kernels_dir + [ $
	'ura111.bsp',$
	'uranus_ringframes_rfrench20201201_v1.tf',$
	'naif0012.tls']
cspice_furnsh,kernels

; specify a time at which to calculate the orientation of the epsilon ring

UTCstr	= '2020 Jan 25 12:00:00'
cspice_str2et,UTCstr,ETsec

Re      = 25559.d0 ; Equatorial radius of Uranus                            

; construct a transformation from Earth equatorial to sky plane frames

cspice_spkezr,'Uranus',ETsec,'J2000','CN','Earth',state,ltime
cspice_recrad,state[0:2],dist,ra,dec

cspice_rotate, dec, 2, mm_temp1                            
cspice_rotate, -ra, 3, mm_temp2                            
cspice_mxm, mm_temp2,mm_temp1, mm_SKY2EEQ                  
cspice_pxform,'URANUS_EQUATORIAL','J2000',ETsec,mm_UR2EEQ

; Uranus pole direction in J2000 coordinates

nhat_pole 	= reform(mm_UR2EEQ[2,*])
nhat_pole	*= Re ; scale by radius of planet

; compute sky plane coordinates of pole

cspice_mtxv,mm_SKY2EEQ,nhat_pole,pole_sky

; compute epsilon ring sky plane coordinates

a_km_ring	= 51149.d0
ae_km_ring	= 406.d0
e_ring		= ae_km_ring/a_km_ring
ntheta  	= 360L*50L                                            
theta   	= dindgen(ntheta)* cspice_rpd()/50.d0                      
rvals		= a_km_ring * (1.d0 - e_ring^2)/(1.d0 + e_ring * cos(theta))
xvals		= rvals * cos(theta)
yvals		= rvals * sin(theta)
ff_km 		= dblarr(ntheta)
gg_km 		= dblarr(ntheta)
cspice_pxform,'URING_EPSILON','J2000',ETsec-ltime,mm_RPL2EEQ
cspice_mtxm,mm_SKY2EEQ,mm_RPL2EEQ, mm_RPL2SKY
for nn=0,ntheta-1 do begin
	vec_ring = [xvals[nn],yvals[nn],0.]
      	cspice_mxv,mm_RPL2SKY, vec_ring,vec_sky
       	ff_km[nn] = vec_sky[1]
       	gg_km[nn] = vec_sky[2]
endfor

; plot the results

set_plot,'PS'
ps_figure 	= 'plot_epsilon_ring_example_IDL.ps'
device,file	= ps_figure,xsize=7,ysize=7,yoffset=0.5,/inch
!P.font = 0

title   = 'Epsilon ring from Earth '+UTCstr
xtitle  = 'East (km)'                                              
ytitle  = 'North (km)'                                             
xrange  = -3.0*[-1,1]*re                                           
yrange  =  3.0*[-1,1]*re                                           

; draw outline of circular planet on sky

xx      = cos(theta)                                               
yy      = sin(theta)                                               
x       = Re * xx
y       = Re * yy
plot,/iso,x,y,xrange=xrange,yrange=yrange,/xstyle,/ystyle,$
        title = title,xtitle=xtitle,ytitle=ytitle

; mark the visible pole location

if pole_sky[0] lt 0.d0 then $
	plots,pole_sky[1],pole_sky[2],psym=4 $
else $
	plots,-pole_sky[1],-pole_sky[2],psym=4

; draw epsilon ring and periapse location

oplot,ff_km,gg_km
plots,ff_km[0],gg_km[0],psym=1,symsize=3 ; mark periapse

device,/close
print,'Saved plot as '+ps_figure
end