# ring_longitude_example.py # # Revisions: # 2020 Jul 21 - rfrench@wellesley.edu - original version # 2020 Dec 05 - rfrench@wellesley.edu - newer frame kernel # # Illustrate several methods of computing the longitude of # periapse of a ring, both by explicit computations and # using matrix transformations based on two different # uranus_ring_frames *tf frame kernels: # # 1) uranus_ringframes_french_et_al_1988_v1.tf # French et al. (1988) Icarus 73, 349 Table XIV # Orbital elements in the B1950 frame. # 2) uranus_ringframes_french20201201_v1.tf # based on a ring orbit fit performed on # 2020 Dec 1 and used as the basis for # computing the geometry of Uranus ring occultations # submitted to NASA's Planetary Data System # in 2020 by Richard G. French, Wellesley College. # This kernel is in the J2000 frame # # The frame kernels are text files and they contain sufficiently # detailed information to enable users to construct their own for # any ring orbit model and Uranus pole direction. # # In this example program: # The periapse longitude is calculated in the J2000 frame # using each kernel for the ten classical Uranian rings, for # a single event time. ######################################################## import numpy as np import spiceypy as spice VERBOSE = True # set to True for additional printed output, else False # !!! Modify this line to point to your local kernels directory !!! kernels_dir = '/Volumes/dione_raid2/Research/kernels/' tls_kernel = 'naif0012.tls' # leap-seconds kernel Uranus_frame_kernels = [ 'uranus_ringframes_french_et_al_1988_v1.tf', 'uranus_ringframes_rfrench20201201_v1.tf'] # the framenames for the rings: UEQframename = 'URANUS_EQUATORIAL' # Uranus equatorial frame framenames = ['URING_6','URING_5','URING_4', 'URING_ALPHA','URING_BETA','URING_ETA', 'URING_GAMMA','URING_DELTA','URING_LAMBDA', 'URING_EPSILON'] # specify the event time at which to calculate the longitude of periapse event_UTC = '1997 March 10 20:00:00' # **************** END OF USER INPUT **************************** nUframe_kernels = len(Uranus_frame_kernels) nframenames = len(framenames) # MAIN LOOP OVER URANUS FRAME KERNELS for nUframe_kernel in range(0,nUframe_kernels): Uranus_frame_kernel = Uranus_frame_kernels[nUframe_kernel] print('Uranus frame kernel: ',Uranus_frame_kernel) spice.kclear # start off fresh, only one Uranus frame kernel is valid at a time spice.furnsh(kernels_dir + tls_kernel) spice.furnsh(kernels_dir + Uranus_frame_kernel) ETsec = spice.str2et(event_UTC) # get the event time, now that the tls kernel is available # obtain the reference frame used for the orbital elements UEQframeID = spice.namfrm(UEQframename) sUEQframeID = str(UEQframeID) TKsearch = 'TKFRAME_' + sUEQframeID + '_RELATIVE' Frame_elements = spice.gcpool(TKsearch,0,1,6) Frame_elements = Frame_elements[0] if VERBOSE: print('This is in the '+Frame_elements+' inertial frame') # obtain the Uranus pole direction TKsearch = 'TKFRAME_' + sUEQframeID + '_ANGLES' pole_dir = spice.gdpool(TKsearch,0,2) RApole = -pole_dir[0] - 90.0 DEpole = pole_dir[1] + 90.0 print('Pole direction (' + Frame_elements + ') : ',RApole,DEpole) # if not in J2000 coordinates, convert to J2000 if Frame_elements != 'J2000': MMtoJ2000 = spice.pxform(Frame_elements,'J2000',ETsec) nhat_pole = spice.radrec(1.0,RApole*spice.rpd(),DEpole*spice.rpd()) nhat_pole_J2000 = spice.mxv(MMtoJ2000,nhat_pole) len,RApoleJ2000,DEpoleJ2000 = spice.recrad(nhat_pole_J2000) print('Pole direction (J2000) : ',RApoleJ2000*spice.dpr(), DEpoleJ2000*spice.dpr()) # loop over all the ring frames.... for nframe in range (0,nframenames): framename = framenames[nframe] print(framename) # obtain the corresponding frameID and construct # prefix string for variables to be extracted from the kernel pool frameID = spice.namfrm(framename) prefix = 'FRAME_'+str(frameID)+'_' # obtain the epoch of the ring elements, and print the result in UTC epoch_ETsec = spice.gdpool(prefix+'EPOCH',0,1) epoch_UTC = spice.et2utc(epoch_ETsec,'C',5) if VERBOSE: print('Epoch of orbital elements (UTC) : '+epoch_UTC[0]) # obtain node and node rate at epoch from kernel pool node_vals = spice.gdpool(prefix+'ANGLE_1_COEFFS',0,2) node_deg_epoch = -node_vals[0] node_rate_deg_day = -node_vals[1] * spice.spd() if VERBOSE: print('Node and node rate: ',node_deg_epoch,node_rate_deg_day) # obtain ring inclination (just for information, not used below) incl_vals = spice.gdpool(prefix+'ANGLE_2_COEFFS',0,1) inclination_deg = -incl_vals[0] if VERBOSE: print('Inclination (deg) :',inclination_deg) # obtain argument of periapse and rate at epoch from the kernel pool arg_peri_vals = spice.gdpool(prefix+'ANGLE_3_COEFFS',0,2) arg_peri_deg_epoch = -arg_peri_vals[0] arg_peri_rate_deg_day = -arg_peri_vals[1] * spice.spd() # form periapse longitude and rate at epoch apse_deg_epoch = (arg_peri_deg_epoch + node_deg_epoch) % 360.00 apse_rate_deg_day = arg_peri_rate_deg_day + node_rate_deg_day if VERBOSE: print('Apse and apse rate: ',apse_deg_epoch,apse_rate_deg_day) # explicitly calculate node and apse longitudes at event time node_deg_event = node_deg_epoch + node_rate_deg_day * (ETsec - epoch_ETsec)/spice.spd() node_deg_event = ((node_deg_event % 360.0) + 360.0) % 360.0 node_deg_event = node_deg_event[0] apse_deg_event = apse_deg_epoch + apse_rate_deg_day * (ETsec - epoch_ETsec)/spice.spd() apse_deg_event = ((apse_deg_event % 360.0) + 360.0) % 360.0 apse_deg_event = apse_deg_event[0] # compute the node from cross product of pole of ring plane and pole of # Uranus pole direction, all in Uranus equatorial coordinates # mm_RPL2UR is the matrix for transforming from the Ring Plane # to the URanus equatorial frame mm_RPL2UR= spice.pxform(framename,'URANUS_EQUATORIAL',ETsec) # use the transformation matrix to determine epsilon ring pole direction # in the Uranus equatorial frame nhat_pole_ringplane = [0.0,0.0,1.0] nhat_pole_ring_ur = spice.mxv(mm_RPL2UR, nhat_pole_ringplane) # an equivalent way to form the pole direction by extracting a the final # column of the transformation matrix: # nhat_pole_ring_ur = [mm_RPL2UR[0,2],mm_RPL2UR[1,2],mm_RPL2UR[2,2]] nhat_pole_planet_ur = [0.0,0.0,1.0] node_line_UR = spice.vcrss(nhat_pole_planet_ur,nhat_pole_ring_ur) # calculate node longitude in degrees from the node line... node_deg_event_calc1 = np.arctan2(node_line_UR[1],node_line_UR[0]) * spice.dpr() node_deg_event_calc1 = (node_deg_event_calc1 + 360.0) % 360.0 # or use the appropriate elements of the transformation matrix... # (note that array indices are in different order from IDL) node_deg_event_calc2 = np.arctan2(mm_RPL2UR[0,2],-mm_RPL2UR[1,2]) * spice.dpr() node_deg_event_calc2 = (node_deg_event_calc2 + 360.0) % 360.0 print('Compare the three different calculations of the node longitude at the event time:') print(node_deg_event) print(node_deg_event_calc1) print(node_deg_event_calc2) # The longitude of periapse in the ring plane relative to the ascending node # is the argument of periapse: arg_peri_event = ((apse_deg_event - node_deg_event) + 360.0) % 360.0 if VERBOSE: print('argument of periapse at the event time = ',arg_peri_event) #compute node line in ring plane node_line_rp = spice.mtxv(mm_RPL2UR,node_line_UR) arg_peri1 = (-np.arctan2(node_line_rp[1],node_line_rp[0]) *spice.dpr() + 360.0) % 360. # combine to form periapse longitude as a broken angle apse_deg_event_calc1 = (arg_peri_event + node_deg_event) % 360.0 # For low-inclination rings, the longitude of periapse can be approximated from the # longitude of the projection of the periapse line in the ring plane into the equatorial # plane: peri_ur_plane = spice.mxv(mm_RPL2UR,[1.0,0.0,0.0]) apse_deg_event_projection = (np.arctan2(peri_ur_plane[1],peri_ur_plane[0]) * spice.dpr() +360.0) % 360.0 print('Compare the two different calculations of the periapse longitude at the event time,') print('and the approximate result obtained by projecting the periapse line in the ring plane') print('into the equatorial plane.') print(apse_deg_event_calc1) print(apse_deg_event) print(apse_deg_event_projection) print(' ')