; occgeom_example2.pro ; ; Example program to compute the geometry of an earth-based ; stellar occultation for a ring event observed on the KAO ; aircraft during the Uranus occultation of 1977 March 10. ; ; Revisions: ; 2020 Jul 21 - rfrench@wellesley.edu - development version ; 2020 Dec 05 - rfrench@wellesley.edu - new tf file, new orbit fit ; Requirements: ; NAIF toolkit implementation in IDL (ICY) ; kernels files listed below ; ********* MAIN ROUTINE ***************** program = 'occgeom_example2.pro' ; Specify all input information ; local directory for kernels kernels_dir = '../../../../kernels/' site_kernel = 'urkao_v1.bsp' ; contains KAO flight path 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 = 'u0' star_catalog = 'Hipparcos' starID = '71567' RAstar_J2000 = 219.5492129d0 DEstar_J2000 = -14.95473933d0 pmRA_mas_yr = -52.61d0 pmDE_max_yr = -12.21d0 parallax_mas = 1.77d0 star_epochJD = 'JD 2448349.0625' ; corrections to star direction, from ring orbit fit dRA_mas = 7.2893775575624078d0 dDE_mas = -5.948861066428534d0 Obsname = 'Kuiper Airborne Observatory' Obs = 'KAS' ; specific event to calculate bundleID = 'uranus_occ_u0_kao_91cm' UTC_obs = '1977 MAR 10 20:24:48.2297' direction = 'ingress' ringname = '6' ; 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 = 0.d0 ; time offset correction in sec R_obs = 41877.6303305763722165d0 R_model = 41877.1346172405610560d0 longitude_deg = 41.2277028983717173d0 sinB = -0.8045319672611818d0 anomaly_deg = 159.3809737855939090d0 ; ring orbital elements and normal mode parameters, if present akm = 41837.3190487974134157d0 ae_km = 42.5463454260044998d0 pi0_deg = 60.0156139308620098d0 mode_wavenumber = 2 mode_da = 0.d0 ; no mode, but keep code present mode_longitude_deg = 0.e0 mode_OmegaP_degday = 0.e0 ; ----------- 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 ; not used in example 2 ;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] 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_example2.pro' print print,'PDART bundleID: ',bundleID print,'Occultation of star u0 by ring 6 during 1977 March 10 occultation' print,'as observed from the Kuiper Airborne Observatory' print,'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,F20.14)' print,'compared to R_obs = ',R_obs,format='(A,F20.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 plane pole (J2000) = ',raRP*cspice_dpr(),decRP*cspice_dpr(),format='(A,2F24.14)' end