# occgeom_example1.py # # Example program to compute the geometry of an earth-based # stellar occultation for a ring event observed during the # Uranus occultation of star u14 on 22 April 1982. # # Revisions: # 2020 Jul 21 - rfrench@wellesley.edu - development version # 2020 Dec 05 - rfrench@wellesley.edu - use updated kernels, input data # Requirements: # NAIF toolkit implementation in IDL (ICY) # kernels files listed below # # ********* MAIN ROUTINE ***************** import array as arr import numpy as np import spiceypy as spice program = 'occgeom_example1.py' # Specify all input information # local directory for kernels # !!! Modify this line to point to your local kernels directory !!! kernels_dir = '../../../../kernels/' site_kernel = 'RAJobs_U111+rgf9.spk' # contains observatory coordinates kernels = [ 'ura111.bsp' # contains planetary ephemerides ,'earthstns_itrf93_040916.bsp' # earth rotation ,'earth_720101_031229.bpc' # earth rotation ,'pck00010.tpc' # planetary contstants ,'uranus_ringframes_rfrench20201201_v1.tf' # Uranus ring frames ,'naif0012.tls' # leap seconds kernel ,site_kernel # described above ] # specify details of observation planet = 'URANUS' event = 'u14' star_catalog = 'Hipparcos' starID = '79085' RAstar_J2000 = 242.1493474e0 DEstar_J2000 = -20.80743248e0 pmRA_mas_yr = -1.16e0 pmDE_max_yr = 0.45e0 parallax_mas = -6.e0 star_epochJD = 'JD 2448349.0625' # corrections to star direction, from ring orbit fit V0204 dRA_mas = 8.8145056956578607e0 dDE_mas = -21.8548569599863640e0 Obsname = 'OPMT' Obs = '586' latitude_deg = 42.e0 + 56.e0/60.e0 + 11.6e0/3600.e0 Elongitude_deg = 0.e0 + 8.e0/60.e0 + 32.3e0/3600.e0 altitude_km = 2.890e0 # specific event to calculate bundleID = 'uranus_ringocc_u14_opmt_200cm' UTC_obs = '1982 APR 22 01:37:00.9786' direction = 'ingress' ringname = 'delta' # values from ring orbit model V0204 produced on 2020 Dec 01 GM = 5.793951322e6 J2 = 3.510561e-03 RrefJ2 = 25559.0e0 RApole = 77.3111427895026821e0 DEpole = 15.1721876765447306e0 epoch_elements = 'UTC Jan 01, 1987 12:00:00' frame_heliocentric= 'J2000' dtoff = 3.6890647775855578e0 # time offset correction in sec R_obs = 48295.9251448240538593e0 R_model = 48297.1535286743601318e0 longitude_deg = 36.1801969754847050e0 sinB = -0.9637606993140888e0 anomaly_deg = 325.0064839021191574e0 # ring orbital elements and normal mode parameters, if present akm = 48300.4466287911563995e0 ae_km = 0.2878785923405290e0 pi0_deg = 125.8816916709482427e0 mode_wavenumber = 2 mode_da = 3.1621875079657320e0 mode_longitude_deg = 245.5225568967657921e0 mode_OmegaP_degday = 562.5159513898452133e0 # ----------- END OF USER INPUT ---------------- # load kernels kernels_list=[kernels_dir + k for k in kernels] spice.furnsh(kernels_list) # get ID numbers for Earth and planet ID_Earth = spice.bodn2c('EARTH') ID_Planet= spice.bodn2c(planet) # get Earth radii, if computing geocentric observer coordinates radii = spice.bodvar(ID_Earth,"RADII",3) a = radii[0] b = radii[2] # compute ET of the observation ETsec = spice.str2et(UTC_obs) # Use site kernel to grab position of Earth observer SPK = kernels_dir + site_kernel handle = spice.dafopr(SPK) void = spice.dafbfs(handle) found = spice.daffna() ND = 2 NI = 1 while found: result = spice.dafgs() dc,ic = spice.dafus(result,ND,NI) site = spice.dafgn() void = spice.boddef(site,ic[0]) found = spice.daffna() spice.dafcls(handle) this_obs = Obs starg,ltime = spice.spkezr(this_obs,ETsec,frame_heliocentric,'NONE','EARTH') R_E0_t0 = starg[0:3] Rdot_E0_t0 = starg[3:6] # For comparison, compute the geocentric body-fixed coordinates of the observatory # two ways: e2 = 1.e0 - (b/a)**2 phi = latitude_deg * spice.rpd() Lambda = Elongitude_deg * spice.rpd() Nphi = a/np.sqrt(1.e0 - e2*np.sin(phi)**2) h = altitude_km Xobs = (Nphi +h) * np.cos(phi) * np.cos(Lambda) Yobs = (Nphi +h) * np.cos(phi) * np.sin(Lambda) Zobs = ((b/a)**2 * Nphi + h) * np.sin(phi) # now use the site kernel starg_obs,ltime = spice.spkezr(this_obs, ETsec,'ITRF93','NONE','EARTH') # (Comparison is printed at the end) R_E_t0_state = spice.spkssb(ID_Earth,ETsec, frame_heliocentric) R_E_t0 = R_E_t0_state[0:3] Rdot_E_t0 = R_E_t0_state[3:6] # heliocentric planet barycenter position at ti ID_obs = spice.bodn2c(this_obs) receive,ltime = spice.ltime(ETsec,ID_obs,"<-",ID_Planet) ETsec_ti0 = ETsec - ltime R_S_ti0_state = spice.spkssb(ID_Planet,ETsec_ti0, frame_heliocentric) R_S_ti0 = R_S_ti0_state[0:3] Rdot_S_ti0 = R_S_ti0_state[3:6] # correct star position for proper motion and parallax km_per_AU = 149597870.691e0 spd = spice.spd() dpy = 365.25e0 spy = dpy*spd dpr = spice.dpr() rpd = spice.rpd() pmRA = pmRA_mas_yr pmDEC = pmDE_max_yr # correct the star position for offsets determined from ring orbit fit DEC_J2000 = DEstar_J2000 + 1.e-3 * dDE_mas/3600.e0 RA_J2000 = RAstar_J2000 + 1.e-3 * dRA_mas/(3600.e0*np.cos(DEstar_J2000 * rpd)) # and correct for proper motion pmRA_deg = pmRA/(np.cos(DEC_J2000*rpd)*3600.e0*1000.e0) pmDEC_deg = pmDEC/(3600.e0*1000.e0) ETsec_epc = spice.str2et(star_epochJD) ETsec_occ = ETsec dt_years = (ETsec_occ - ETsec_epc)/spy RA_2 = RA_J2000 + pmRA_deg *dt_years DEC_2 = DEC_J2000 + pmDEC_deg*dt_years nhat_sun_star_ppm_corr = spice.radrec(1.e0,RA_2*rpd, DEC_2*rpd) R_scT = R_E_t0 + R_E0_t0 PLX = parallax_mas # note that if PLX <0, no parallax correction is applied if PLX > 0.e0: D_Star = km_per_AU/np.sin(PLX*rpd/(3600.e0*1000.e0)) rect_sun_star = nhat_sun_star_ppm_corr*D_Star nhat_star_ppm_plx_corr = (rect_sun_star - R_scT)/spice.vnorm(rect_sun_star - R_scT) else: nhat_star_ppm_plx_corr = nhat_sun_star_ppm_corr # compute Uranus pole direction in J2000 coordinates pole = [0.e0,0.e0,1.e0] mmUEQ2J2000 = spice.pxform('URANUS_EQUATORIAL',frame_heliocentric,ETsec_ti0) nhat_pole_eq = spice.mxv(mmUEQ2J2000,pole) # confirm that equator pole direction agrees with ring orbit result radius,raU,decU = spice.recrad(nhat_pole_eq) # compute pole of delta ring at observed time backdated by light travel time ringframe = 'URING_' + ringname.upper() mm_RP2EEQ = spice.pxform(ringframe,frame_heliocentric,ETsec_ti0) mm_EEQ2RP = spice.xpose(mm_RP2EEQ) nhat_pole_delta_ring = spice.mxv(mm_RP2EEQ,pole) radius,raRP,decRP = spice.recrad(nhat_pole_delta_ring) # compute B and P, needed for GR deflection calculation nhat_star = nhat_star_ppm_plx_corr len,ra,dec = spice.recrad(nhat_star) mm_temp1 = spice.rotate(dec,2) mm_temp2 = spice.rotate(-ra,3) mm_SKY2EEQ = spice.mxm(mm_temp2,mm_temp1) mm_EEQ2SKY = spice.xpose(mm_SKY2EEQ) mm_PL2EEQ = spice.pxform('URANUS_EQUATORIAL',frame_heliocentric,ETsec) mm_EEQ2PL = spice.xpose(mm_PL2EEQ) nhat_pole = spice.mxv(mm_PL2EEQ,[0.e0,0.e0,1.e0]) nhat_pole_sky = spice.mxv(mm_EEQ2SKY,nhat_pole) Pdeg = (np.arctan2(-nhat_pole_sky[1],-nhat_pole_sky[2]) * dpr + 360.e0) % 360.e0 Bdeg = np.arcsin(nhat_pole_sky[0]) * dpr input_esec_t0 = ETsec input_esec_ti0 = ETsec_ti0 input_R_S_ti0 = R_S_ti0 input_Rdot_S_ti0= Rdot_S_ti0 nhat_star = nhat_star_ppm_plx_corr nhat_pole = nhat_pole_delta_ring dR_S = [0.e0,0.e0,0.e0] input_R_E_t0 = R_E_t0 input_Rdot_E_t0 = Rdot_E_t0 input_R_E0_t0 = R_E0_t0 input_Rdot_E0_t0= Rdot_E0_t0 input_dtoff = dtoff c = spice.clight() # apply time offset dtoff = input_dtoff if dtoff == 0: esec_t0 = input_esec_t0 esec_ti0 = input_esec_ti0 R_S_ti0 = input_R_S_ti0 Rdot_S_ti0 = input_Rdot_S_ti0 R_E_t0 = input_R_E_t0 R_E0_t0 = input_R_E0_t0 else: esec_t0 = input_esec_t0 + dtoff esec_ti0 = input_esec_ti0 + dtoff R_S_ti0 = input_R_S_ti0 + dtoff*input_Rdot_S_ti0 Rdot_S_ti0 = input_Rdot_S_ti0 R_E_t0 = input_R_E_t0 + dtoff*input_Rdot_E_t0 R_E0_t0 = input_R_E0_t0 + dtoff*input_Rdot_E0_t0 n_iter = 3 R_EE0_t0 = spice.vadd(R_E_t0,R_E0_t0) mR_EE0_t0 = spice.vminus(R_EE0_t0) Delta = (esec_t0 - esec_ti0) * c # initial value nhat_stard = nhat_star # initial value D = 0.e0 rho = [0.e0,0.e0,0.e0] rho_mag = 0.e0 if GM == 0.e0: GR_correc=0 else: GR_correc=1 # compute info needed for J2 term for GR correction if GR_correc ==1: GRfactor = 4.e0*GM/c**2 if J2 != 0.e0: J2fac = J2*(RrefJ2*np.cos(Bdeg*rpd))**2 radius,ra_star,dec_star = spice.recrad(nhat_star) mm_PSK2EEQ = spice.eul2m(-ra_star,dec_star,Pdeg*rpd,3,2,1) mm_EEQ2PSK = spice.xpose(mm_PSK2EEQ) for n in range(1,n_iter+1): nhats_dot = spice.vdot(nhat_stard,nhat_pole) esec_ti = esec_t0 - Delta/c dR_S_ti0 = spice.vscl(esec_ti - esec_ti0,Rdot_S_ti0) R_S_ti = spice.vadd(R_S_ti0,dR_S_ti0) R_S_ti = spice.vadd(R_S_ti,dR_S) # ephemeris offset -- zero for this example vecsum = spice.vadd(R_S_ti,mR_EE0_t0) Delta = spice.vdot(vecsum,nhat_pole)/nhats_dot mR_S_ti = spice.vminus(R_S_ti) vecsum = spice.vadd(R_EE0_t0,mR_S_ti) Delta_nhat_stard = spice.vscl(Delta,nhat_stard) R_SI = spice.vadd(vecsum,Delta_nhat_stard) factor_rho = 0.e0 factor_rho_sphere = 0.e0 if GR_correc == 1 and n == 1: mvecsum = spice.vminus(vecsum) D = spice.vdot(mvecsum,nhat_stard) D_nhat_stard = spice.vscl(D,nhat_stard) rho = spice.vadd(D_nhat_stard,vecsum) rho_mag = spice.vnorm(rho) factor0 = GRfactor/rho_mag**2 if J2 == 0.e0: # spherical bending: factor_rho = spice.vscl(factor0,rho) else: # J2 bending: vecxy = spice.mxv(mm_EEQ2PSK,rho) x = vecxy[1] y = vecxy[2] xyterm1 = J2fac*((3.e0*y**2-x**2)/rho_mag**4) eps_x = factor0 * (1.e0 - xyterm1)*x xyterm2 = J2fac*((3.e0*x**2-y**2)/rho_mag**4) eps_y = factor0 * (1.e0 + xyterm2)*y eps_vec = spice.vpack(0.e0,eps_x,eps_y) factor_rho = spice.mxv(mm_PSK2EEQ,eps_vec) eps_norm = spice.vnorm(eps_vec) nhat_stard = spice.vadd(nhat_star,factor_rho) nhat_stard_norm = spice.vnorm(nhat_stard) nhat_stard = spice.vscl(1.e0/nhat_stard_norm,nhat_stard) Delta1_out = Delta esec_ti_out = esec_ti R_SI_out = R_SI nhat_stard_out = nhat_stard D_out = D rho_out = rho rho_mag_out = rho_mag nhat_pole = spice.mxv(mm_RP2EEQ,[0.e0,0.e0,1.e0]) nhat_pole_sky = spice.mxv(mm_EEQ2SKY,nhat_pole) BdegRP = np.arcsin(nhat_pole_sky[0]) * dpr R_SI_ringplane = spice.mxv(mm_EEQ2RP,R_SI_out) # calculate ring longitude and anomaly R_SI_planet = spice.mxv(mm_EEQ2PL,R_SI_out) rip_londeg = (np.arctan2(R_SI_planet[1],R_SI_planet[0]) * dpr + 360.e0) % 360.e0 rip_anomaly = (np.arctan2(R_SI_ringplane[1],R_SI_ringplane[0]) * dpr + 360.e0) % 360.e0 # compute model radius, radius residual, time residual e = ae_km/akm r_model_calc = akm * (1.e0 - e**2)/(1.e0 + e * np.cos(rip_anomaly * rpd)) ET_elements = spice.str2et(epoch_elements) dangle_rad = ((esec_ti_out + dtoff - ET_elements) * mode_OmegaP_degday/spice.spd() + mode_longitude_deg) * rpd arg = mode_wavenumber * (rip_londeg * rpd - dangle_rad) da_mode = -mode_da * np.cos(arg) r_model_calc += da_mode # print results print('Results of occgeom_example1.py') print('') print('PDART bundleID: ',bundleID) print('Occultation of star u14 by delta ring during 1982 April 22 occultation') print('as observed from OPMT at '+UTC_obs) print('') print('Comparison of results from this calulation and ring orbit fit:') print('') print('Ring plane intercept radius = ',spice.vnorm(R_SI_out)) print('compared to ring orbit fit result = ',R_obs) print('') print('sinB of ring plane = ',np.sin(BdegRP * rpd)) print('compared to ring orbit fit result = ',sinB) print('') print('Confirm that ring intercept vector has nearly zero z-element in ring plane:') print('R_SI_ringplane = ',R_SI_ringplane) print('') print('model radius = ',r_model_calc,' compared to ',R_model) print('ring longitude = ',rip_londeg,' compared to ',longitude_deg) print('ring anomaly = ',rip_anomaly,' compared to ',anomaly_deg) print('') print('Pole directions:') print('Uranus pole (J2000) = ',raU*dpr,decU*dpr) print('From ring orbit fit = ',RApole,DEpole) print('ring plane pole (J2000) = ',raRP*dpr,decRP*dpr) print('') print('Comparison of geocentric position calculated two ways:') print('Explicit calculation = ',Xobs,Yobs,Zobs) print('Using a site kernel = ',starg_obs[0],starg_obs[1],starg_obs[2])