; ring_longitude_example.pro ; ; Revisions: ; 2020 Jul 21 - rfrench@wellesley.edu - original version ; 2020 Dec 05 - rfrench@wellesley.edu - updated tf 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 B1950 or J2000 frame ; using each kernel for the ten classical Uranian rings, for ; a single event time. ; ###################################################################### VERBOSE = 1 ; set to 1 for additional printed output ; !!! 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 = n_elements(Uranus_frame_kernels) ; MAIN LOOP OVER URANUS FRAME KERNELS for nUframe_kernel = 0,nUframe_kernels-1 do begin Uranus_frame_kernel = Uranus_frame_kernels[nUframe_kernel] print,'Uranus frame kernel: '+Uranus_frame_kernel cspice_kclear ; start off fresh, only one Uranus frame kernel is valid at a time cspice_furnsh,kernels_dir + tls_kernel cspice_furnsh,kernels_dir + Uranus_frame_kernel cspice_str2ET,event_UTC,ETsec ; get the event time, now that the tls kernel is available ; obtain the reference frame used for the orbital elements cspice_namfrm,UEQframename,UEQframeID sUEQframeID = strtrim(string(UEQframeID),2) TKsearch = 'TKFRAME_' + sUEQframeID + '_RELATIVE' cspice_gcpool,TKsearch,0,1,6,Frame_elements,found if VERBOSE then print,'This is in the '+Frame_elements+' inertial frame' ; obtain the Uranus pole direction TKsearch = 'TKFRAME_' + sUEQframeID + '_ANGLES' cspice_gdpool,TKsearch,0,2,pole_dir,found RApole = -pole_dir[0] - 90.d0 DEpole = pole_dir[1] + 90.d0 print,'Pole direction (' + Frame_elements[0] + ') : ',RApole,DEpole ; if not in J2000 coordinates, convert to J2000 if Frame_elements ne 'J2000' then begin cspice_pxform,Frame_elements[0],'J2000',ETsec,MMtoJ2000 cspice_radrec,1.d0,RApole*cspice_rpd(),DEpole*cspice_rpd(),nhat_pole cspice_mxv,MMtoJ2000,nhat_pole,nhat_pole_J2000 cspice_recrad,nhat_pole_J2000,len,RApoleJ2000,DEpoleJ2000 print,'Pole direction (J2000) : ',RApoleJ2000*cspice_dpr(),$ DEpoleJ2000*cspice_dpr() endif nframes = n_elements(framenames) ; loop over all the ring frames.... for nframe = 0,nframes-1 do begin framename = framenames[nframe] print,framename ; obtain the corresponding frameID and construct ; prefix string for variables to be extracted from the kernel pool cspice_namfrm,framename,frameID prefix = 'FRAME_'+strtrim(string(frameID),2)+'_' ; obtain the epoch of the ring elements, and print the result in UTC cspice_gdpool,prefix+'EPOCH',0,1,epoch_ETsec,found cspice_et2utc,epoch_ETsec,'C',5,epoch_UTC if VERBOSE then print,'Epoch of orbital elements (UTC) : '+epoch_UTC ; obtain node and node rate at epoch from kernel pool cspice_gdpool,prefix+'ANGLE_1_COEFFS',0,2,node_vals,found node_deg_epoch = -node_vals[0] node_rate_deg_day = -node_vals[1] * cspice_spd() if VERBOSE then print,'Node and node rate: ',node_deg_epoch,node_rate_deg_day ; obtain ring inclination (just for information, not used below) cspice_gdpool,prefix+'ANGLE_2_COEFFS',0,2,incl_vals,found inclination_deg = -incl_vals[0] if VERBOSE then print,'Inclination (deg) :',inclination_deg ; obtain argument of periapse and rate at epoch from the kernel pool cspice_gdpool,prefix+'ANGLE_3_COEFFS',0,2,arg_peri_vals,found arg_peri_deg_epoch = -arg_peri_vals[0] arg_peri_rate_deg_day = -arg_peri_vals[1] * cspice_spd() ; form periapse longitude and rate at epoch apse_deg_epoch = (arg_peri_deg_epoch + node_deg_epoch) mod 360.d0 apse_rate_deg_day = arg_peri_rate_deg_day + node_rate_deg_day if VERBOSE then 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)/cspice_spd() node_deg_event = ((node_deg_event mod 360.d0) + 360.d0) mod 360.d0 apse_deg_event = apse_deg_epoch + apse_rate_deg_day * (ETsec - epoch_ETsec)/cspice_spd() apse_deg_event = ((apse_deg_event mod 360.d0) + 360.d0) mod 360.d0 ; 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 cspice_pxform,framename,'URANUS_EQUATORIAL',ETsec,mm_RPL2UR ; use the transformation matrix to determine epsilon ring pole direction ; in the Uranus equatorial frame nhat_pole_ringplane = [0.d0,0.d0,1.d0] cspice_mxv,mm_RPL2UR, nhat_pole_ringplane,nhat_pole_ring_ur ; an equivalent way to form the pole direction by extracting a the final ; column of the transformation matrix: ; nhat_pole_ring_ur = reform(mm_RPL2UR[2,*]) nhat_pole_planet_ur = [0.d0,0.d0,1.d0] node_line_UR = crossp(nhat_pole_planet_ur,nhat_pole_ring_ur) ; calculate node longitude in degrees from the node line... node_deg_event_calc1 = atan(node_line_UR[1],node_line_UR[0]) * cspice_dpr() node_deg_event_calc1 = (node_deg_event_calc1 + 360.d0) mod 360.d0 ; or use the appropriate elements of the transformation matrix... node_deg_event_calc2 = atan(mm_RPL2UR[2,0],-mm_RPL2UR[2,1]) * cspice_dpr() node_deg_event_calc2 = (node_deg_event_calc2 + 360.d0) mod 360.d0 print,'Compare the three different calculations of the node longitude at the event time:' format = '(F18.13)' print,node_deg_event,format=format print,node_deg_event_calc1,format=format print,node_deg_event_calc2,format=format ; 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.d0) mod 360.d0 if VERBOSE then print,'argument of periapse at the event time = ',arg_peri_event ;compute node line in ring plane cspice_mtxv,mm_RPL2UR,node_line_UR,node_line_rp arg_peri1 = (-atan(node_line_rp[1],node_line_rp[0]) *cspice_dpr() + 360.d0) mod 360 ; combine to form periapse longitude as a broken angle apse_deg_event_calc1 = (arg_peri_event + node_deg_event) mod 360.d0 ; 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: cspice_mxv,mm_RPL2UR,[1.d0,0.d0,0.d0],peri_ur_plane apse_deg_event_projection = (atan(peri_ur_plane[1],peri_ur_plane[0]) * cspice_dpr() +360.d0) mod 360.d0 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,format=format print,apse_deg_event,format=format print,apse_deg_event_projection,format=format print endfor ; end of loop over ring frames endfor ; end of loop over frame kernels end