; occgeom_example1.pro ; ; 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 Apr 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 ***************** program = 'occgeom_example1.pro' ; Specify all input information ; local directory for kernels 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.1493474d0 DEstar_J2000 = -20.80743248d0 pmRA_mas_yr = -1.16d0 pmDE_max_yr = 0.45d0 parallax_mas = -6.d0 star_epochJD = 'JD 2448349.0625' ; corrections to star direction, from ring orbit fit dRA_mas = 8.8145056956578607d0 dDE_mas = -21.8548569599863640d0 Obsname = 'OPMT' Obs = '586' latitude_deg = 42.d0 + 56.d0/60.d0 + 11.6d0/3600.d0 Elongitude_deg = 0.d0 + 8.d0/60.d0 + 32.3d0/3600.d0 altitude_km = 2.890d0 bundleID = 'uranus_ringocc_u14_opmt_200cm' UTC_obs = '1982 APR 22 01:37:00.9786' direction = 'ingress' ringname = 'delta' ; values from ring orbit model GM = 5.793951322d6 J2 = 3.510561d-03 RrefJ2 = 25559.0d0 RApole = 77.3111427895026821d0 DEpole = 15.1721876765447306d0 epoch_elements = 'UTC Jan 01, 1987 12:00:00' frame_heliocentric= 'J2000' dtoff = 3.6890647775855578d0 ; time offset correction in sec R_obs = 48295.9251448240538593d0 R_model = 48297.1535286743601318d0 longitude_deg = 36.1801969754847050d0 sinB = -0.9637606993140888d0 anomaly_deg = 325.0064839021191574d0 ; ring orbital elements and normal mode parameters, if present akm = 48300.4466287911563995d0 ae_km = 0.2878785923405290d0 pi0_deg = 125.8816916709482427d0 mode_wavenumber = 2 mode_da = 3.1621875079657320d0 mode_longitude_deg = 245.5225568967657921d0 mode_OmegaP_degday = 562.5159513898452133d0 ; ----------- END OF USER INPUT ---------------- ; load kernels cspice_kclear nkernels= n_elements(kernels) for nk=0,nkernels-1 do cspice_furnsh,kernels_dir + kernels[nk] ; get ID numbers for Earth and planet cspice_bodn2c,'EARTH',ID_Earth,found cspice_bodn2c,planet,ID_Planet,found ; get Earth radii, if computing geocentric observer coordinates cspice_bodvar,ID_earth,"RADII",radii a = radii[0] b = radii[2] ; compute ET of the observation cspice_str2ET,UTC_obs,ETsec ; Use site kernel to grab position of Earth observer (obscure NAIF syntax!) SPK = kernels_dir + site_kernel cspice_dafopr, SPK, handle cspice_dafbfs, handle cspice_daffna, found NI = 6 ND = 2 while (found) do begin cspice_dafgs, ND, NI, dc, ic siteID = ic[0] cspice_dafgn, site cspice_boddef,site,siteID cspice_daffna, found endwhile cspice_dafcls, handle this_obs = obs cspice_spkezr, this_obs, ETsec,frame_heliocentric,'NONE','EARTH',starg,ltime R_E0_t0 = starg[0:2] Rdot_E0_t0 = starg[3:5] ; For comparison, compute the geocentric body-fixed coordinates of the observatory ; two ways: e2 = (1.d0 - (b/a)^2) phi = latitude_deg * cspice_rpd() lambda = Elongitude_deg * cspice_rpd() Nphi = a/sqrt(1.d0 - e2*sin(phi)^2) h = altitude_km Xobs = (Nphi +h) * cos(phi) * cos(lambda) Yobs = (Nphi +h) * cos(phi) * sin(lambda) Zobs = ((b/a)^2 * Nphi + h) * sin(phi) ; now use the site kernel cspice_spkezr, this_obs, ETsec,'ITRF93','NONE','EARTH',starg_obs,ltime ; (Comparison is printed at the end) cspice_spkssb, ID_Earth,ETsec, frame_heliocentric, R_E_t0_state R_E_t0 = R_E_t0_state[0:2] Rdot_E_t0 = R_E_t0_state[3:5] ; heliocentric planet barycenter position at ti cspice_bodn2c,this_obs,ID_obs,found cspice_ltime, ETsec,ID_obs, "<-",ID_Planet,receive,ltime ETsec_ti0 = ETsec - ltime cspice_spkssb, ID_Planet,ETsec_ti0, frame_heliocentric, R_S_ti0_state R_S_ti0 = R_S_ti0_state[0:2] Rdot_S_ti0 = R_S_ti0_state[3:5] ; correct star position for proper motion and parallax km_per_AU = 149597870.691d0 spd = cspice_spd() dpy = 365.25d0 spy = dpy*spd dpr = cspice_dpr() rpd = cspice_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.d-3 * dDE_mas/3600.d0 RA_J2000 = RAstar_J2000 + 1.d-3 * dRA_mas/(3600.d0*cos(DEstar_J2000 * cspice_rpd())) ; and correct for proper motion pmRA_deg = pmRA/(cos(DEC_J2000*rpd)*3600.d0*1000.d0) pmDEC_deg = pmDEC/(3600.d0*1000.d0) cspice_str2ET,star_epochJD,ETsec_epc 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 cspice_radrec, 1.d0, RA_2[0]*rpd, DEC_2[0]*rpd, nhat_sun_star_ppm_corr R_scT = R_E_t0 + R_E0_t0 PLX = parallax_mas ; note that if PLX <0, no parallax correction is applied if PLX gt 0.d0 then begin D_Star = km_per_AU/sin(PLX*rpd/(3600.d0*1000.d0)) rect_sun_star = nhat_sun_star_ppm_corr*D_Star nhat_star_ppm_plx_corr = (rect_sun_star - R_scT)/cspice_vnorm(rect_sun_star - R_scT) endif else nhat_star_ppm_plx_corr = nhat_sun_star_ppm_corr ; compute Uranus pole direction in J2000 coordinates pole = [0.d0,0.d0,1.d0] cspice_pxform,'URANUS_EQUATORIAL',frame_heliocentric,ETsec_ti0,rotate cspice_mxv,rotate,pole,nhat_pole_eq ; confirm that equator pole direction agrees with ring orbit result cspice_recrad,nhat_pole_eq,radius,raU,decU ; compute pole of ring at observed time backdated by light travel time ringframe = 'URING_' + strupcase(ringname) cspice_pxform,ringframe,frame_heliocentric,ETsec_ti0,mm_RP2EEQ mm_EEQ2RP = transpose(mm_RP2EEQ) cspice_mxv,mm_RP2EEQ,pole,nhat_pole_delta_ring cspice_recrad,nhat_pole_delta_ring,radius,raRP,decRP ; compute B and P, needed for GR deflection calculation nhat_star = nhat_star_ppm_plx_corr cspice_recrad,nhat_star,len,ra,dec cspice_rotate, dec, 2, mm_temp1 cspice_rotate, -ra, 3, mm_temp2 cspice_mxm, mm_temp2,mm_temp1, mm_SKY2EEQ mm_EEQ2SKY = transpose(mm_SKY2EEQ) cspice_pxform,'URANUS_EQUATORIAL',frame_heliocentric,ETsec,mm_PL2EEQ mm_EEQ2PL = transpose(mm_PL2EEQ) cspice_mxv,mm_PL2EEQ,[0.d0,0.d0,1.d0],nhat_pole cspice_mxv,mm_EEQ2SKY,nhat_pole,nhat_pole_sky Pdeg = (atan(-nhat_pole_sky[1],-nhat_pole_sky[2]) * cspice_dpr() + 360.d0) mod 360.d0 Bdeg = asin(nhat_pole_sky[0]) * cspice_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.d0,0.d0,0.d0] 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 = cspice_clight() ; apply time offset dtoff = input_dtoff[0] ; ensure that it is a scalar case dtoff of 0: begin 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 end else: begin 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 end endcase n_iter = 3 cspice_vadd,R_E_t0,R_E0_t0,R_EE0_t0 cspice_vminus,R_EE0_t0,mR_EE0_t0 Delta = (esec_t0 - esec_ti0) * c ; initial value nhat_stard = nhat_star ; initial value D = 0.d0 rho = [0.d0,0.d0,0.d0] rho_mag = 0.d0 if GM[0] ne 0.d0 then GR_correc=1 else GR_correc=0 ; compute info needed for J2 term for GR correction if GR_correc then begin GRfactor = 4.d0*GM[0]/c^2 if J2 ne 0.d0 then begin rpd = cspice_rpd() J2fac = J2*(RrefJ2*cos(Bdeg*rpd))^2 cspice_recrad,nhat_star,radius,ra_star,dec_star cspice_eul2m,-ra_star,dec_star,Pdeg*rpd,3,2,1,mm_PSK2EEQ mm_EEQ2PSK = transpose(mm_PSK2EEQ) endif endif for n=1,n_iter do begin nhats_dot = cspice_vdot(nhat_stard,nhat_pole) esec_ti = esec_t0 - Delta/c cspice_vscl,(esec_ti - esec_ti0),Rdot_S_ti0,dR_S_ti0 cspice_vadd,R_S_ti0,dr_S_ti0,R_S_ti cspice_vadd,R_S_ti,dr_S,R_S_ti ; ephemeris offset cspice_vadd,R_S_ti,mR_EE0_t0,vecsum Delta = cspice_vdot(vecsum,nhat_pole)/nhats_dot cspice_vminus,R_S_ti,mR_S_ti cspice_vadd,R_EE0_t0,mR_S_ti,vecsum cspice_vscl,Delta,nhat_stard,Delta_nhat_stard cspice_vadd,vecsum,Delta_nhat_stard,R_SI factor_rho = 0.d0 factor_rho_sphere = 0.d0 if GR_correc and n eq 1 then begin cspice_vminus,vecsum,mvecsum D = cspice_vdot(mvecsum,nhat_stard) cspice_vscl,D,nhat_stard,D_nhat_stard cspice_vadd,D_nhat_stard,vecsum,rho rho_mag = cspice_vnorm(rho) factor0 =GRfactor/rho_mag^2 if J2 eq 0.d0 then begin ; spherical bending: cspice_vscl,factor0,rho,factor_rho endif else begin ; J2 bending: cspice_mxv,mm_EEQ2PSK,rho,vecxy x = vecxy[1] y = vecxy[2] xyterm1 = J2fac*((3.d0*y^2-x^2)/rho_mag^4) eps_x = factor0 * (1.d0 - xyterm1)*x xyterm2 = J2fac*((3.d0*x^2-y^2)/rho_mag^4) eps_y = factor0 * (1.d0 + xyterm2)*y cspice_vpack,0.d0,eps_x,eps_y,eps_vec cspice_mxv,mm_PSK2EEQ,eps_vec,factor_rho eps_norm= cspice_vnorm(eps_vec) endelse cspice_vadd,nhat_star,factor_rho,nhat_stard nhat_stard_norm= cspice_vnorm(nhat_stard) cspice_vscl,1.d0/nhat_stard_norm,nhat_stard,nhat_stard endif endfor ; end of iteration loop 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 cspice_mxv,mm_RP2EEQ,[0.d0,0.d0,1.d0],nhat_pole cspice_mxv,mm_EEQ2SKY,nhat_pole,nhat_pole_sky BdegRP = asin(nhat_pole_sky[0]) * cspice_dpr() cspice_mxv,mm_EEQ2RP,R_SI_out,R_SI_ringplane ; calculate ring longitude and anomaly cspice_mxv,mm_EEQ2PL,R_SI_out,R_SI_planet rip_londeg = (atan(R_SI_planet[1],R_SI_planet[0]) * cspice_dpr() + 360.d0) mod 360.d0 rip_anomaly = (atan(R_SI_ringplane[1],R_SI_ringplane[0]) * cspice_dpr() + 360.d0) mod 360.d0 ; compute model radius, radius residual, time residual e = ae_km/akm r_model_calc = akm * (1.d0 - e^2)/(1.d0 + e * cos(rip_anomaly * cspice_rpd())) cspice_str2et,epoch_elements,ET_elements dangle_rad = ((esec_ti_out + dtoff - ET_elements) * mode_OmegaP_degday/cspice_spd() + mode_longitude_deg) $ * cspice_rpd() arg = mode_wavenumber * (rip_londeg * cspice_rpd() - dangle_rad) da_mode = -mode_da * cos(arg) r_model_calc += da_mode ; print results print,'Results of occgeom_example1.pro' print print,'PDART bundleID: ',bundleID print,'Occultation of star u14 by the delta ring during the 22 Apr 1982 occultation' print,'as observed from OMPT at '+UTC_obs print print,'Comparison of results from this calculation and ring orbit fit:' print print,'Ring plane intercept radius = ',cspice_vnorm(R_SI_out),format='(A,2F20.14)' print,'compared to ring orbit fit result = ',R_obs,format='(A,2F20.14)' print print,'sinB of ring plane = ',sin(BdegRP * cspice_rpd()),format='(A,2F20.14)' print,'compared to ring orbit fit result = ',sinB,format='(A,2F20.14)' 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 = ',string(r_model_calc,format='(F20.14)'),' compared to ',string(R_model,format='(F20.14)') print,'ring longitude = ',string(rip_londeg,format='(F20.14)'),' compared to ',string(longitude_deg,format='(F20.14)') print,'ring anomaly = ',string(rip_anomaly,format='(F20.14)'),' compared to ',string(anomaly_deg,format='(F20.14)') print print,'Pole directions:' print,'Uranus pole (J2000) = ',raU*cspice_dpr(),decU*cspice_dpr(),format='(A,2F24.14)' print,'From orbit model = ',RApole,DEpole,format='(A,2F24.14)' print,'ring pole (J2000) = ',raRP*cspice_dpr(),decRP*cspice_dpr(),format='(A,2F24.14)' print print,'Comparison of geocentric position calculated two ways:' format = '(A,3F15.6)' print,'Explicit calculation = ',Xobs,Yobs,Zobs, format=format print,'Using a site kernel = ',starg_obs[0:2], format=format end