;P090421_pluto_occpred_v4.pro ; ; Simple routine to predict occultation for a given observatory ; ; Revisions: ; v1: ; 2007 March 28 - rfrench - original version for this event ; v2: ; 2007 July 17 - rfrench - updated star positions, observatories ; v3: ; 2007 July 18 - rfrench - updated star positions. ; v4: ; 2007 July 19 - rfrench - add Lowell Pluto offset ; P598.2_v4: ; 2008 Aug 23 - rfrench - Aug 25 occultation ; P090421: ; 2009 Apr 13 - rfrench - predictions for P090421 ; 2009 Apr 14 - rfrench - add Lowell2 star position function hms2rad,h,m,s return,cspice_rpd() * 15.d0 * (double(h) + double(m)/60.d0 + double(s)/3600.d0) end function dms2rad,h,m,s ; if < 0, ALL args must be <0 return,cspice_rpd() * (double(h) + double(m)/60.d0 + double(s)/3600.d0) end ; ******** MAIN ROUTINE ************* @set_tek_color kernels_file='P090421event_v4.ker' cspice_kclear cspice_furnsh,kernels_file UTC_ref = '2009-Apr 21 22:30' dtsec = 1.d0 n_times = 7201L target = 'PLUTO' ref_frame='J2000' ab_corr = 'CN' target_radius0=1400. ; use the larger of this value and IAU target radius ; empirically, Jim uses ~1400 km for time of occ ;;;target_radius0=1200. ; use the larger of this value and IAU target radius ; therefore, to use default radius, set this to zero ; ; Jim Elliot uses 1200 km, IAU value is 1195 rlimit = 5. ; in units of Pluto radius rlimit = 2.5 ; in units of Pluto radius xrange0 = [1,-1]*rlimit yrange0 = -xrange0 athick = 3 thick = 2 xtitle= 'E (km)' ytitle= 'N (km)' charsize= 1.5 charsize_label= 0.75 ; star positions (J2000, at epoch of event) ra_stars0 = [hms2rad(18,12,53.5913) $ ; 2MASS ,hms2rad(18,12,53.5913) $ ; USNO B1.0 ,hms2rad(18,12,53.5862) $ ; Brazil ,hms2rad(18,12,53.6002) $ ; Lowell ,hms2rad(18,12,53.5917) $ ; Lowell2 ,hms2rad(18,12,53.5903) $ ; Mean ] ; ra_star_offsets_mas = [ $ 0.d0 $ ,0.d0 $ ,0.d0 $ ,0.d0 $ ,0.d0 $ ,0.d0 $ ] dec_stars0= [ dms2rad(-17,-37,-56.923d0) $ ; 2Mass ,dms2rad(-17,-37,-56.958d0) $ ; USNO B1.0 ,dms2rad(-17,-37,-56.921d0) $ ; Brazil ,dms2rad(-17,-37,-56.858d0) $ ; Lowell ,dms2rad(-17,-37,-56.939d0) $ ; Lowell2 ,dms2rad(-17,-37,-56.905d0) $ ; Mean ] dec_star_offsets_mas = [ $ 0.d0 $ ,0.d0 $ ,0.d0 $ ,0.d0 $ ,0.d0 $ ,0.d0 $ ] ra_star_offsets_rad =( ra_star_offsets_mas/3600000.d0)*cspice_rpd() dec_star_offsets_rad=(dec_star_offsets_mas/3600000.d0)*cspice_rpd() dec_stars = dec_stars0 +dec_star_offsets_rad ra_stars = ra_stars0 + ra_star_offsets_rad/cos(dec_stars) star_sources=['2MASS','USNO B1.0','Brazil','Lowell','Lowell2','Mean'] which_star_sources = [0,1,2,3,4,5] ; which star sources to use this time lss = which_star_sources ; shorthand star_ref = 0 ; index of star position to use ; Pluto offsets: ra_target_offsets_mas =[ $ 0.d0 $ ; DE413 + PLU-017 , -69.0 $ ; DE418 , -22.0 $ ; PHOT07 , -11.0 $ ; PHOT08 , -94.0 $ ; Sicardy ] dec_target_offsets_mas =[ $ 0.d0 $ ; DE413 + PLU-017 , 30.0 $ ; DE418 , 38.0 $ ; PHOT07 , 105.0 $ ; PHOT07 , 140.0 $ ; Sicardy ] offset_sources=' '+['DE-413+PLU-017','DE418','PHOT07','PHOT08','Sicardy'] which_offsets=[0,1,2,3,4] loff = which_offsets ; shorthand OBSlist =['SAAOLY' $ ,'BOYDEN' $ ,'PRETORIA' $ ,'HESSATOM' $ ,'ETOSHA' $ ,'REUNION' $ ] Obsnames =['SAAO' $ ,'Boyden' $ ,'Pretoria' $ ,'HESS/ATOM' $ ,'Etosha' $ ,'Reunion' $ ] OBScodes=[[ $ 989 $ ,988 $ ,987 $ ,986 $ ,985 $ ,984 $ ] $ + 399000L] which_obs = [0,1,2,3,4,5] program = 'P090421_occpred_v4.pro' psfile = 'P090421_occpred_v4.ps' outfile = 'P090421_occpred_v4.txt' ; ********* END OF USER INPUT ********** n_obs = n_elements(which_obs) n_offsets = n_elements(which_offsets) n_star_sources = n_elements(which_star_sources) ; make storage for times so that we can format a nice table at the end et_immersion = dblarr(n_obs,n_offsets,n_star_sources) et_emersion = dblarr(n_obs,n_offsets,n_star_sources) et_midtime = dblarr(n_obs,n_offsets,n_star_sources) distance_km = dblarr(n_obs,n_offsets,n_star_sources) hostname = hostname() if n_elements(whoami) eq 0 then whoami = mywhoami() if isPS() then begin setps device,/color,file=psfile endif close,1 openw,1,outfile printf,1,'Predictions for 2009 Apr 21 Occultation of P090421 by Pluto' printf,1,'by Richard G. French, Wellesley College ' cd,CURR=cwd printf,1,hostname+':'+cwd+'/'+program printf,1,systime() printf,1 printf,1,'Predictions for the following observatories:' for mob=0,n_obs-1 do begin nob = which_obs[mob] printf,1,' '+OBSnames[nob] endfor printf,1 printf,1,'Key:' printf,1,' Pluto Ephemeris : JPL barycentric Pluto/Charon ephemeris and Pluto ephemeris' printf,1,' dRA, dDec mas : RA, Dec offsets added to the Pluto Ephemeris in mas' printf,1,' Source : Origin of the ephemeris offsets printf,1,' StarPosition : Origin of the adopted star position for a given prediction' ;printf,1,' ; UCAC2 17 53 27.1056 -17 31 31.655 J2000 ;printf,1,' Lowell 31-in 17 53 27.0981 -17 31 31.664 J2000 ;printf,1,' Brazil (Bruno) 17 53 27.1032 -17 31 31.647 J2000 ;printf,1,' MIT 17 53 27.0975 -17 31 31.613 J2000 printf,1,' Distance (km) : Closest approach distance in the sky plane to center of Pluto''s shadow' printf,1,' Immersion (UTC) : Immersion time at R=1400 km radius' printf,1,' Midtime (UTC) : Time at closest approach distance' printf,1,' Emersion (UTC) : Emersion time at R=1400 km radius' printf,1 printf,1,'For all predictions for a given observatory: printf,1,' Closest : Closest Distance, at Midtime, to center of Pluto''s shadow, in sky plane printf,1,' Furthest : Furthest " " " " " " " " " " " ' printf,1,' Earliest : Earliest Immersion, Midtime, Emersion printf,1,' Latest : Latest " " " ' printf,1 printf,1,$ "--------------------------------------------------------------------------------------------------" cspice_bodn2c,target,target_code,found cspice_bodvar, target_code,'RADII',radii target_radius = radii[0]>target_radius0 ; use larger of the two xrange = xrange0*target_radius yrange = yrange0*target_radius theta = findgen(361)*cspice_rpd() xcircle = cos(theta) ycircle = sin(theta) cspice_utc2et,UTC_ref, et_ref ra_vals_spk = dblarr(n_times) dec_vals_spk = dblarr(n_times) ranges = dblarr(n_times) ; for each observatory.... for mob = 0,n_obs-1 do begin nob = which_obs[mob] obsentry= OBSlist[nob] obscode = OBScodes[nob] obs = OBSnames[nob] cspice_boddef,obsentry,OBScode ; compute topocentric Pluto positions at requested times, no corrections yet for n=0,n_times-1 do begin et = et_ref + n*dtsec cspice_spkpos, target, et,ref_frame,ab_corr,obsentry,ptarg,ltime cspice_recrad,ptarg,dist,ra,dec ra_vals_spk[n] = ra dec_vals_spk[n] = dec ranges[n] = ltime * cspice_clight() ; ra_deg = ra * cspice_dpr() ; dec_deg = dec * cspice_dpr() ; cspice_et2utc, et, 'C',0,utcstr ; printf,1,utcstr,ra_deg,dec_deg,format='(A,2X,2(F15.7))' endfor ; end of loop over n times printf,1 printf,1,'P040921 predictions for ',obs printf,1 printf,1,$ "Pluto Ephemeris dRA dDec mas Source StarPosition Distance Immersion Midtime Emersion" printf,1,$ "----------------------------------------------------------------------------------------------------------" ; for each offset... for m_offset=0,n_offsets-1 do begin n_offset = which_offsets[m_offset] ra_target_offset_mas = ra_target_offsets_mas[n_offset] dec_target_offset_mas =dec_target_offsets_mas[n_offset] eph_string = 'DE-413 Plu017 + ' + string( ra_target_offset_mas,$ dec_target_offset_mas,format='("(",F6.1,", ",F6.1,") mas") ')+$ offset_sources[n_offset] ra_target_offset_rad =( ra_target_offset_mas/3600000.d0)*cspice_rpd() dec_target_offset_rad=(dec_target_offset_mas/3600000.d0)*cspice_rpd() ; ra_vals = ra_vals_spk + ra_target_offset_rad ; dec_vals = dec_vals_spk + dec_target_offset_rad dec_vals = dec_vals_spk + dec_target_offset_rad ; Corrected the RA for longitude convergence 2007 July 17 - rfrench ra_vals = ra_vals_spk + ra_target_offset_rad /cos(dec_vals) ; for each star position.... for mstar=0,n_star_sources-1 do begin nstar = which_star_sources[mstar] ra_star = ra_stars[nstar] dec_star =dec_stars[nstar] star_source = star_sources[nstar] rd2xieta, ra_vals, dec_vals, ra_star,dec_star, xi_rad, eta_rad xi_km = xi_rad * ranges eta_km = eta_rad * ranges skysep_km = sqrt(xi_km^2 + eta_km^2) lskysep = where(skysep_km lt rlimit*target_radius,nlskysep) min_skysep_km = min(skysep_km) lclose = where(skysep_km eq min_skysep_km) et_close = et_ref + lclose*dtsec et_midtime[mob,m_offset,mstar] = et_close distance_km[mob,m_offset,mstar] = min_skysep_km cspice_et2utc, et_close, 'C', 0, utcstr str_close = strmid(utcstr,12,8) if mstar eq 0 then begin plot,xcircle*target_radius,ycircle*target_radius,$ xrange=xrange,$ yrange=yrange,position=aspect(),$ xtitle=xtitle,ytitle=ytitle,$ title=obs,/nodata,$ xthick=athick,ythick=athick,charsize=charsize oplot,xcircle*target_radius,ycircle*target_radius,$ thick=thick,color=tek_red plots,0,0,psym=1,symsize=3,color=tek_red yloc = !Y.CRANGE[0]+ 0.90*(!Y.CRANGE[1]-!Y.CRANGE[0]) xyouts,0,yloc,eph_string,align=0.5,$ color=tek_red ; make a simple plot of star and planet offsets km_per_mas = mean(ranges)/(cspice_dpr()*60.d0*60.d0*1000.d0) x0=-0.35 * (!X.CRANGE[1] - !X.CRANGE[0]) + !X.CRANGE[0] y0=0.85 * (!Y.CRANGE[1] - !Y.CRANGE[0]) + !Y.CRANGE[0] x1=1.30 * (!X.CRANGE[1] - !X.CRANGE[0]) + !X.CRANGE[0] y1=0.85 * (!Y.CRANGE[1] - !Y.CRANGE[0]) + !Y.CRANGE[0] dx_mas = 100. dy_mas = 100. dx_km = dx_mas * km_per_mas dy_km = dy_mas * km_per_mas plots,x0+[-1,1]*dx_km,y0+[0,0],color=tek_red,thick=thick plots,x0+[0,0],y0+[-1,1]*dy_km,color=tek_red,thick=thick plots,x1+[-1,1]*dx_km,y1+[0,0],color=tek_blue,thick=thick plots,x1+[0,0],y1+[-1,1]*dy_km,color=tek_blue,thick=thick plots,x0+ [-4,-3,-2,-1,1,2,3,4]*25.*km_per_mas,$ y0+0*[-4,-3,-2,-1,1,2,3,4]*25.*km_per_mas,$ color=tek_red,psym=1 plots,x0+0*[-4,-3,-2,-1,1,2,3,4]*25.*km_per_mas,$ y0+ [-4,-3,-2,-1,1,2,3,4]*25.*km_per_mas,$ color=tek_red,psym=1 plots,x1+ [-4,-3,-2,-1,1,2,3,4]*25.*km_per_mas,$ y1+0*[-4,-3,-2,-1,1,2,3,4]*25.*km_per_mas,$ color=tek_blue,psym=1 plots,x1+0*[-4,-3,-2,-1,1,2,3,4]*25.*km_per_mas,$ y1+ [-4,-3,-2,-1,1,2,3,4]*25.*km_per_mas,$ color=tek_blue,psym=1 xyouts,x0,y0+dy_km,' 100 mas',color=tek_red xyouts,x0,y0-1.2*dy_km,'Pluto',color=tek_red,align=.5 xyouts,x1,y1+dy_km,' 100 mas',color=tek_blue xyouts,x1,y1-1.2*dy_km,'Star',color=tek_blue,align=.5 ra_target_offsets_km = ra_target_offsets_mas * km_per_mas dec_target_offsets_km =dec_target_offsets_mas * km_per_mas plots,x0+ra_target_offsets_km[loff],y0+dec_target_offsets_km[loff],$ psym=2,color=tek_red xyouts,x0+ra_target_offsets_km[loff],y0+dec_target_offsets_km[loff],$ offset_sources[loff],color=tek_red,charsize=charsize_label dra_stars = ra_stars - ra_stars[star_ref] ddec_stars = dec_stars - dec_stars[star_ref] dra_stars_mas = dra_stars*cspice_dpr()*60.d0*60.d0*1000.d0 $ * cos(dec_stars) ddec_stars_mas = ddec_stars*cspice_dpr()*60.d0*60.d0*1000.d0 dra_stars_km = dra_stars_mas * km_per_mas ddec_stars_km = ddec_stars_mas * km_per_mas plots,x1+dra_stars_km[lss],y1+ddec_stars_km[lss],$ psym=2,color=tek_blue xyouts,x1+dra_stars_km[lss],y1+ddec_stars_km[lss],$ star_sources[lss],color=tek_blue,charsize=charsize_label if isPS() then mytimestamp,program,psfile,/LANDSCAPE,$ HOSTNAME=hostname,WHOAMI=whoami endif ; end of mstar=0 case str_immersion='--------' ; reset each time str_emersion ='--------' if nlskysep gt 0 then begin lmax = max(lskysep) oplot,xi_km[lskysep],eta_km[lskysep],color=tek_blue,thick=thick plots,xi_km[lskysep[0]],eta_km[lskysep[0]],psym=4,color=tek_blue plots,xi_km[lmax],eta_km[lmax],psym=4,color=tek_blue et_start = et_ref + lskysep[0]*dtsec cspice_et2utc, et_start, 'C', 0, utcstr xyouts,xi_km[lskysep[0]],eta_km[lskysep[0]]+200.,$ strmid(utcstr,12,8),align=1,charsize=charsize_label et_end = et_ref + lmax*dtsec cspice_et2utc, et_end, 'C', 0, utcstr xyouts,xi_km[lmax],eta_km[lmax]+200.,$ strmid(utcstr,12,8),align=0,charsize=charsize_label xyouts,xi_km[lskysep[0]],eta_km[lskysep[0]]-200.,' '+star_source,$ align=0,color=tek_blue,charsize=charsize_label,$ orient=-15. plots,xi_km[lclose],eta_km[lclose],psym=2,color=tek_blue xyouts,xi_km[lclose],eta_km[lclose]+200.,str_close,$ align=0.5,charsize=charsize_label ; compute immersion and emersion times if min_skysep_km lt target_radius then begin locc=where(skysep_km lt target_radius,nlocc) lim=locc[0] lem=max(locc) et_imm = et_ref + lim*dtsec cspice_et2utc, et_imm, 'C', 0, utcstr str_immersion = strmid(utcstr,12,8) et_em = et_ref + lem*dtsec cspice_et2utc, et_em, 'C', 0, utcstr str_emersion = strmid(utcstr,12,8) plots,xi_km[[lim,lem]],eta_km[[lim,lem]],$ psym=1,color=tek_blue xyouts,xi_km[lim],eta_km[lim],$ align=1,str_immersion,charsize=charsize_label xyouts,xi_km[lem],eta_km[lem],$ align=0,str_emersion,charsize=charsize_label et_immersion[mob,m_offset,mstar] = et_imm et_emersion[mob,m_offset,mstar] = et_em endif else begin str_immersion='--------' str_emersion ='--------' endelse endif ; end of nlskysep gt 0 str_miss = string(round(min(skysep_km)),format='(I6," km")') if mstar eq 0 then this_eph_string=eph_string else this_eph_string=' ' if mstar eq 0 and m_offset ne 0 then printf,1 printf,1,this_eph_string,star_source,str_miss,$ str_immersion,str_close,str_emersion,$ format='(A-52,1x,A-15,1X,A10,3(1X,A8))' endfor ; end of loop over nstar ; if ~isPS() then STOP endfor ; end of loop over target offsets ; compute earliest, latest times these_midtimes = et_midtime[mob,*,*] midtime_earliest = min(these_midtimes) cspice_et2utc, midtime_earliest, 'C', 0, utcstr_midtime_earliest str_midtime_earliest = strmid(utcstr_midtime_earliest,12,8) midtime_latest = max(these_midtimes) cspice_et2utc, midtime_latest, 'C', 0, utcstr_midtime_latest str_midtime_latest = strmid(utcstr_midtime_latest,12,8) these_distances = distance_km[mob,*,*] these_immersions = et_immersion[mob,*,*] these_emersions = et_emersion[mob,*,*] l =where(these_immersions ne 0,nl) if nl eq 0 then begin str_immersion_earliest='-------- ' str_immersion_latest=' --------' str_emersion_earliest='-------- ' str_emersion_latest=' --------' endif else begin immersion_earliest = min(these_immersions[l]) cspice_et2utc, immersion_earliest, 'C', 0, utcstr_immersion_earliest str_immersion_earliest=strmid(utcstr_immersion_earliest,12,8)+' ' immersion_latest = max(these_immersions[l]) cspice_et2utc, immersion_latest, 'C', 0, utcstr_immersion_latest str_immersion_latest=strmid(utcstr_immersion_latest,12,8)+' ' emersion_earliest = min(these_emersions[l]) cspice_et2utc, emersion_earliest, 'C', 0, utcstr_emersion_earliest str_emersion_earliest=' '+strmid(utcstr_emersion_earliest,12,8) emersion_latest = max(these_emersions[l]) cspice_et2utc, emersion_latest, 'C', 0, utcstr_emersion_latest str_emersion_latest=' '+strmid(utcstr_emersion_latest,12,8) endelse printf,1 printf,1,$ " Distance" printf,1,string(string(round(min(these_distances)),format='(I6," km")'),format='(61x,"Closest ",A9)') printf,1,string(string(round(max(these_distances)),format='(I6," km")'),format='(61x,"Furthest ",A9)') printf,1 printf,1,$ " Immersion Midtime Emersion" printf,1,$ " Earliest "+ $ str_immersion_earliest+str_midtime_earliest+str_emersion_earliest printf,1,$ " Latest "+ $ str_immersion_latest+str_midtime_latest+str_emersion_latest printf,1,$ "----------------------------------------------------------------------------------------------------------" printf,1 endfor ; end of loop over observatories close,1 print,'Saved predictions in '+outfile spawn,'cat '+outfile if isPS() then begin device,/close print,'Saved plot as '+psfile endif cspice_unload,kernels_file end