; create_uranus_ringframes_french_et_al_1988_v1.pro ; ; Create a uranus ringframes tf kernel file from orbital elements ; in Table XIV of French et al. 1988 Icarus 73, 349. ; Test implementation by reproducing Table II of Gresh et al. 1989. ; ; Revisions: ; 2019 Feb 13 - rfrench - original version ; 2019 Dec 21 - rfrench - use local directory for portability ; 2020 Jan 10 - rfrench - change name of output file ; 2020 Jul 21 - rfrench - minor format changes program = 'create_uranus_ringframes_french_et_al_1988_v1.pro' out_tf_kernelfile = 'uranus_ringframes_french_et_al_1988_v1.tf' tls = 'naif0012.tls' ; END OF USER INPUT cspice_kclear cspice_furnsh,tls frame_ring_parameters = 'B1950' epoch_ring_parameters = 'UTC Mar 10, 1977 20:00:00' cspice_UTC2ET,epoch_ring_parameters,ETepoch_ring_parameters cspice_deltet,ETepoch_ring_parameters,"ET",delta cspice_et2UTC,ETepoch_ring_parameters+delta,'C',6,cal_ring_parameters c = cal_ring_parameters ; shorthand frame_ring_epoch = '@' + strmid(c,0,4) + '-' + strmid(c,5,3) + '-' + $ strmid(c,9,2) + '/' + strmid(c,12,999) RPdeg = 76.5969d0 ; B1950 From Table XIV French et al 1988 DPdeg = 15.1117d0 RP_dot = 0.d0 DP_dot = 0.d0 format = '(" TKFRAME_1799000_ANGLES = ( ",F16.10,", ",F16.10,", 0 )")' pole_entry = string(-90.d0 - RPdeg,DPdeg-90.d0,format=format) ringnames = ['6','5','4','Alpha','Beta','Eta','Gamma','Delta','Lambda','Epsilon'] ringlabels = ['Ring 6','Ring 5','Ring 4','Alpha','Beta','Eta','Gamma','Delta','Lambda','Epsilon'] parmlabel = ' Ring peri(deg) peri_dot inc(deg) node(deg) node_dot' hashline = ' -----------------------------------------------------------' nrings = n_elements(ringnames) ; Table XIV French et al. 1988 GM = 5.793939d6 J2 = 3.34343d-3 J4 = -2.885d-5 J6 = 0.d0 J8 = 0.d0 J10 = 0.d0 J12 = 0.d0 J14 = 0.d0 Rref_km = 26200.d0 str_elements = [ $ {ringname:'6',ringlabel:'Ring 6', akm:41837.15d0, e:1.103d-3, perilon_deg:242.80d0, peri_dot_degsec:0.d0, incl_deg:0.0616d0, nodelon_deg:012.12d0, node_dot_degsec:0.d0} $ ,{ringname:'5',ringlabel:'Ring 5', akm:42234.82d0, e:1.899d-3, perilon_deg:170.31d0, peri_dot_degsec:0.d0, incl_deg:0.0536d0, nodelon_deg:286.57d0, node_dot_degsec:0.d0} $ ,{ringname:'4',ringlabel:'Ring 4', akm:42570.91d0, e:1.059d-3, perilon_deg:127.28d0, peri_dot_degsec:0.d0, incl_deg:0.0323d0, nodelon_deg:089.26d0, node_dot_degsec:0.d0} $ ,{ringname:'Alpha',ringlabel:'Alpha', akm:44718.45d0, e:0.761d-3, perilon_deg:333.24d0, peri_dot_degsec:0.d0, incl_deg:0.0152d0, nodelon_deg:063.08d0, node_dot_degsec:0.d0} $ ,{ringname:'Beta',ringlabel:'Beta', akm:45661.03d0, e:0.442d-3, perilon_deg:224.88d0, peri_dot_degsec:0.d0, incl_deg:0.0051d0, nodelon_deg:310.05d0, node_dot_degsec:0.d0} $ ,{ringname:'Eta',ringlabel:'Eta', akm:47175.91d0, e:0.004d-3, perilon_deg:228.10d0, peri_dot_degsec:0.d0, incl_deg:0.0011d0, nodelon_deg:188.73d0, node_dot_degsec:0.d0} $ ; FROM TABLE XV ,{ringname:'Gamma',ringlabel:'Gamma', akm:47626.87d0, e:0.109d-3, perilon_deg:132.10d0, peri_dot_degsec:0.d0, incl_deg:0.0015d0, nodelon_deg:251.30d0, node_dot_degsec:0.d0} $ ,{ringname:'Delta',ringlabel:'Delta', akm:48300.12d0, e:0.004d-3, perilon_deg:216.70d0, peri_dot_degsec:0.d0, incl_deg:0.0011d0, nodelon_deg:260.70d0, node_dot_degsec:0.d0} $ ; FROM TABLE XIV ,{ringname:'Lambda',ringlabel:'Lambda', akm:50023.94d0, e:0.000d-3, perilon_deg:000.00d0, peri_dot_degsec:0.d0, incl_deg:0.0000d0, nodelon_deg:000.00d0, node_dot_degsec:0.d0} $ ,{ringname:'Epsilon',ringlabel:'Epsilon',akm:51149.32d0, e:7.936d-3, perilon_deg:214.97d0, peri_dot_degsec:0.d0, incl_deg:0.0002d0, nodelon_deg:246.60d0, node_dot_degsec:0.d0} $ ] nrings = n_elements(str_elements) ring_entries = [] for n=0,nrings-1 do begin str = str_elements[n] a = str.akm e = str.e incl_rad = str.incl_deg * cspice_rpd() str.peri_dot_degsec = ringfit_domdt_uranus(a,e,incl_rad,J2,J4,J6,J8,J10,J12,J14,GM, Rref_km,0) * cspice_dpr() str.node_dot_degsec = ringfit_domdt_uranus(a,e,incl_rad,J2,J4,J6,J8,J10,J12,J14,GM, Rref_km,1) * cspice_dpr() ring_entry = string(str.ringlabel,str.perilon_deg,str.peri_dot_degsec*cspice_spd(),str.incl_deg,str.nodelon_deg,$ str.node_dot_degsec*cspice_spd(),format='(5x,A-8,2x,F10.5,1x,F8.5,2x,F9.5,1x,F9.5,1x,F9.5)') ring_entries=[ring_entries,ring_entry] str_elements[n] = str endfor ringframenames = 'URING_' + strupcase(ringnames) frameIDs = string(1799001L+lindgen(nrings),format='(I7)') frame_defs = [' '] for n=0,nrings-1 do begin str = str_elements[n] node_val = str.nodelon_deg node_dot_val = str.node_dot_degsec peri_dot_val = str.peri_dot_degsec ; BUG fixed 2019 Dec 26 ; arg_peri_val = str.perilon_deg arg_peri_val = str.nodelon_deg - str.perilon_deg inc_val = str.incl_deg rfn = ringframenames[n] fID = frameIDs[n] pre = ' FRAME_'+fID+'_' frame_defs = [frame_defs,$ ' FRAME_'+rfn+ ' = '+frameIDs[n],$ pre + 'NAME = '+" '"+rfn+"'",$ pre + 'CLASS = 5',$ pre + 'CLASS_ID = '+fID,$ pre + 'CENTER = 799',$ pre + "RELATIVE = 'URANUS_EQUATORIAL'",$ pre + "DEF_STYLE = 'PARAMETERIZED'",$ pre + "FAMILY = 'EULER'",$ pre + 'EPOCH = '+frame_ring_epoch,$ pre + 'AXES = ( 3 1 3 )',$ pre + "UNITS = 'DEGREES'",$ pre + 'ANGLE_1_COEFFS = (' + string(-node_val,format='(F13.8)') + ' ' + $ string(-node_dot_val,format='(E21.14," )")'),$ pre + 'ANGLE_2_COEFFS = (' + string(-inc_val,format='(F13.8," 0.0 )")'), $ pre + 'ANGLE_3_COEFFS = (' + string(arg_peri_val,format='(F13.8)') + ' '+ $ string(node_dot_val-peri_dot_val,format='(E21.14," )")'),$ ' '] endfor openw,/get_lun,olun,out_tf_kernelfile cd,curr=cwd contents= ['\begintext',out_tf_kernelfile,'produced on '+systime(),$ getenv('USER')+'@'+hostname()+':'+cwd+'/'+program,$ 'Orbital elements are from French et al. (1988) Icarus 73, 349 Table XIV',$ 'This file contains definitions of Uranus equatorial and ring reference frames',$ '','The Uranus equatorial frame is defined by:',' ',$ ' - The +Z axis is parallel to the spin angular velocity vector of',$ ' Uranus at the specified epoch (in this case the B1950 epoch).','',$ ' - The +X axis is parallel to the cross product of the',$ ' +Z axis of the base frame (in this case the B1950 frame)',$ ' and the +Z axis of the equatorial frame.','',$ ' This frame is implemented as a "TK" (type 4) frame. The Euler',$ ' angle sequence used here is the inverse of that used to define',$ ' the orientation of the equatorial frame relative to the base frame.',' ',$ '\begindata','',$ ' FRAME_URANUS_EQUATORIAL = 1799000',$ " FRAME_1799000_NAME = 'URANUS_EQUATORIAL'",$ ' FRAME_1799000_CLASS = 4',$ ' FRAME_1799000_CLASS_ID = 1799000',$ ' FRAME_1799000_CENTER = 799',$ " TKFRAME_1799000_RELATIVE = 'B1950'",$ " TKFRAME_1799000_SPEC = 'ANGLES'",$ " TKFRAME_1799000_UNITS = 'DEGREES'",$ ' TKFRAME_1799000_AXES = ( 3, 1, 3 )',$ pole_entry,'','','\begintext','',$ ' The ring frames are implemented using the Euler family of the',$ ' parameterized dynamic frame class.',' ',$ ' The epoch of the elements is',$ ' '+epoch_ring_parameters,' ',$ ' The corresponding TDB date is',$ ' '+cal_ring_parameters+ ' TDB','',$ ' The Euler angle sequence for the ring frames is','',$ ' [ arg of periapse ] [ inclination ] [ longitude of node ]',$ ' 3 1 3',' ',$ ' The specific constants for the orientation of the rings are:',$ ' ',$ parmlabel,hashline,ring_entries,$ ' ',$ ' Euler frames require the inverse of this sequence, since the',$ " mapping from the Euler frame to its base frame is what's defined.",$ ' The required sequence is:',' ',$ ' [ -longitude of node ] [ -inclination ] [ -arg of periapse ]',$ ' 3 1 3',' ',$ ' Note that in the frame definitions below, the angular rates are',$ ' scaled from degrees/day to degrees/second, as required by the',$ ' SPICE frames system.',' ','\begindata',$ frame_defs] printf,olun,transpose(contents) free_lun,olun spawn,'cat '+ out_tf_kernelfile print print,'Created '+out_tf_kernelfile print,'**************' cspice_furnsh,out_tf_kernelfile UTC_RIP = '1986 Jan 24 ' + [ $ '19:43:36.672' $ ; ei ,'19:48:41.478' $ ; di ,'19:50:04.813' $ ; gi ,'19:50:59.318' $ ; ti ,'19:54:07.459' $ ; bi ,'19:55:58.385' $ ; ai ,'20:00:26.589' $ ; 4i ,'20:00:57.680' $ ; 5i ,'20:01:51.081' $ ; 6i ,'22:33:33.858' $ ; 6e ,'22:34:17.110' $ ; 5e ,'22:35:06.301' $ ; 4e ,'22:39:26.294' $ ; ae ,'22:41:26.104' $ ; be ,'22:44:28.220' $ ; te ,'22:45:22.825' $ ; ge ,'22:46:43.642' $ ; de ,'22:52:50.317' $ ; ee ] cspice_str2et,'UTC '+UTC_RIP,ETsec_RIP rings = [ $ 'EPSILON' $ ,'DELTA' $ ,'GAMMA' $ ,'ETA' $ ,'BETA' $ ,'ALPHA' $ ,'4' $ ,'5' $ ,'6' $ ] rings = [rings,reverse(rings)] ; egress, too frames = 'URING_' + rings print,'Comparison with Gresh et al. 1989 Icarus 78, 131 -- Table II' print print,'Ring Frame I/E Ring Intercept Time (UTC) pi0(deg) Omega0(deg) RA(deg) Dec(deg)' print,'-----------------------------------------------------------------------------------' for m=0,n_elements(UTC_RIP)-1 do begin if m le 8 then n=8-m else n=m frame = frames[n] ETsec = ETsec_RIP[n] cspice_pxform,frame,'B1950',ETsec,rotate ring_pole = [0.d0,0.d0,1.d0] cspice_mxv,rotate,ring_pole,pole_B1950 cspice_recrad,pole_B1950,len,RA,DE RAdeg = RA * cspice_dpr() DEdeg = DE * cspice_dpr() LL = where(strupcase(str_elements.ringname) eq rings[n],nL) perilon_deg = str_elements[LL].perilon_deg peri_dot_degsec = str_elements[LL].peri_dot_degsec nodelon_deg = str_elements[LL].nodelon_deg node_dot_degsec = str_elements[LL].node_dot_degsec dtsec = ETsec - ETepoch_ring_parameters pi0 = perilon_deg + peri_dot_degsec * dtsec Omega0 = nodelon_deg + node_dot_degsec * dtsec pi0 = ((pi0 mod 360.d0) + 360.d0) mod 360.d0 Omega0 = ((Omega0 mod 360.d0) + 360.d0) mod 360.d0 if n lt n_elements(UTC_RIP)/2 then IE = 'I' else IE='E' print,frame,IE,UTC_RIP[n],pi0,Omega0,RAdeg,DEdeg,$ format = '(A-15,2x,A1,2x,A,2x,F8.3,2x,F8.3,2x,F8.4,2x,F7.4)' endfor print,'-----------------------------------------------------------------------------------' cspice_kclear ; start fresh ; here is a code fragment to illustrate the use of the frame kernel cspice_furnsh,['naif0012.tls','uranus_ringframes_French_et_al_1988.tf'] cspice_str2et,'UTC 1986 Jan 24 20:01:51.081',ETsec_ring6_ingress ring_pole = [0.d0,0.d0,1.d0] cspice_pxform,'URING_6','B1950',ETsec_ring6_ingress,rotate cspice_mxv,rotate,ring_pole,pole_B1950 cspice_recrad,pole_B1950,len,RA,DE RAdeg = RA * cspice_dpr() DEdeg = DE * cspice_dpr() print,'B1950 pole direction for Voyager 2 RSS ring 6 ingress:' print,RAdeg,DEdeg,format='("RA = ",F7.4," Dec = ",F7.4," deg")' cspice_pxform,'URING_6','J2000',ETsec_ring6_ingress,rotate cspice_mxv,rotate,ring_pole,pole_J2000 cspice_recrad,pole_J2000,len,RA,DE RAdeg = RA * cspice_dpr() DEdeg = DE * cspice_dpr() print,'J2000 pole direction for Voyager 2 RSS ring 6 ingress:' print,RAdeg,DEdeg,format='("RA = ",F7.4," Dec = ",F7.4," deg")' end