;;File Version: 3.5 ;;Last Modified: April 22, 2009 ;;Modified by: JEC PRO CG_GEOMETER_ENGINE_6_5, time, instrument, target, phase, emission, incidence, $ ra, dec, lat, lon, ringplanerad, scloc, scvel, ra_targ, dec_targ, $ lat_subpt, lon_subpt, lat_sunpt, lon_sunpt, alt, phase_targ,emission_targ,$ incidence_targ, sun_elevation, earth_elevation, sunearth, ckerror, $ rayheight, latocc, lonocc, losdist, ra_full, dec_full, lat_full, $ lon_full, rrad_full, sc_wrt_plan, state_sun, state_planet, inc_full, $ sc_sun, occpoint_i, ringoccphi, f_peri, rayradius, ramlat, ramlon, $ ringsolar_full, PIXEL = PIXEL, ABCORR = abcorr, $ CORNERS = CORNERS, SLIT = SLIT, MOON = MOON, $ SOLAR = SOLAR , geom_ref; , SUBTYPE = SUBTYPE, DETIC = DETIC ;;+ ;; ;; This program computes various geometric properties at an input ;; time for a specified instrument relative to a specified target. ;; The resulting values are passed back in the parameter list. ;; ;; << Inputs >> ;; ;; Time: the time at which computations will be made, in NAIF format. ;; Instrument: a string identifying the instrument whose boresight ;; will be used for pointing-dependent calculations. Must be identified ;; in NAIF format. A value of 'NONE' indicates no instrument data ;; will be computed. ;; Target: a string identifying the intended target of the observation. ;; ra, dec: input with a physical value to indicate inertial pointing ;; rather than getting pointing from a C kernel. Otherwise, input with ;; a value of -1000D0. ;; ;; << Outputs >> ;; ;; phase: phase angle on the specified target at the specified time at ;; the point intersected by the instrument boresight or the nearest ;; point on the line of sight to the target body (deg). ;; emission: emission angle on the specified target at the specified ;; time at the point intersected by the instrument boresight or the nearest ;; point on the line of sight to the target body (deg). ;; incidence: solar incidence angle, measured from the local target surface ;; normal, at the point on the target that is intersected by the ;; instrument boresight at the specified time or the nearest ;; point on the line of sight to the target body (deg). ;; ringplanerad: Radius of occultation path in Saturn's ring plane (km). ;; scloc: 3-element vector with postion of S/C wrt to target (km) in ;; inertial frame. ;; scvel: 3-element vector with velocity of S/C wrt to target (km/sec) in ;; inertial frame. ;; ra_targ: Right ascension of center of target body (deg). ;; dec_targ: Declination of center of target body (deg). ;; lat_subpt: latitude of sub-spacecraft point (using nearpoint method) (deg). ;; lon_subpt: longitude of sub-spacecraft point (using nearpoint method) (deg). ;; lat_sunpt: latitude of sub-solar point (using nearpoint method) (deg). ;; lon_sunpt: longitude of sub-solar point (using nearpoint method) (deg). ;; alt: distance above target surface (km). ;; phase_targ: phase angle at the sub-S/C point (deg). ;; emission_targ: emission angle at the sub-S/C point (deg). ;; incidence_targ: solar incidence angle at the sub-S/C point (deg). ;; ckerror: set to 1 if C kernel data not found, otherwise 0. ;; rayheight: minimum separation of the look vector and the NAIF-defined ;; surface of the body (km). ;; rayradius: minimum distance between the look vector and the center of the ;; target body (km). ;; latocc and lonocc: latitude and longitude of the point on the NAIF-defined ;; surface of the target body that is closest to the look vector. ;; ra, dec: right ascension and declination of the instrument axis. ;; sc_wrt_plan: 6-element state of spacecraft with respect to planet ;; (not the target ;; if the target is a satellite, but with respect to the planet that satellite ;; orbits). ;; state_sun: 6-element position/velocity of Sun with respect to the target ;; body. ;; state_planet: 6-element position/velocity of planet with respect to the ;; solar system barycenter. ;; ramlat, ramlon: the latitude and longitude of the point on the body ;; that is in the ram direction as it orbits the central body (planet ;; or Sun). (deg). ;; ;; << Restrictions >> ;; ;; Uses the CSPICE DLM for IDL. ;; Requires existence of 'geometerkernels.txt' in the same directory as ;; this procedure with a list of all SPICE kernels for instruments and ;; the time period in question. ;; ra must be passed in with a value of -1000D0 if using C kernel pointing, ;; otherwise ra and dec values passed in are assumed to be the fixed pointing. ;; ;; << History >> ;; ;; Joshua Colwell ;; Begun: February 27, 2003. ;; February 28, 2003: Getting RA/Dec of boresight correctly, tested against ;; CASPER output. ;; March 5, 2003: Add NO_KERNEL_LOAD keyword. Phase angle is working, solar and ;; emission angles do not check with my expectations or with SOA output. RA ;; and Dec of boresight checked with SOA and is okay. ;; June 3, 2003: changing calls to PXFORM to CKGP so that we can handle ;; dropouts in C kernel data. June 4: checked against an Alpha PsA observation: ;; get correct RA/Dec values for boresight. ;; June 18, 2003: Add CKERNEL keyword. ;; June 30, 2003: I have lat/lon working, but phase, emission, incidence ;; still need to be checked. ;; Checked lat/lon visually against SOA display for one C kernel. ;; September 27, 2003: Changed NO_KERNEL_LOAD keyword to KERNEL_LOAD. Compute ;; ring plane intercept radius (ringplanerad) ;; September 30, 2003: retrieve basic state vector information on spacecraft. ;; Add computation of non-pointing related geometery. Checked against SOA using ;; both trajectory example and a C kernel. Emission angle disagrees with SOA. ;; Small discrepancy with S/C altitude. Everything else checks out. ;; October 1, 2003: rename GEOMETER_ENGINE (formerly GEOMETER) ;; October 2, 2003: add phase_targ to output. ;; October 3, 2003: simplify angle calcs by using ILLUM. ;; Calculate angles for ring intercept. ;; All Ring parameters have NOT been checked yet. ;; October 27, 2003: error handling for c kernel. ;; January 14, 2004: Replace RING keyword with "Rings" value for target input. ;; January 15, 2004: add rayheight, latocc, and lonocc output. ;; January 16, 2004: Pixel values, rayheight, latocc, and lonocc all need ;; verification. ;; January 20, 2004: some bug fixing. ;; January 22, 2004: Test case comparison with Jupiter observation, C kernel ;; display from CASPER for 2000 Dec 29 2:46. Have row ordering consistent ;; with CASPER. But, more pixels are on Jupiter than CASPER shows can fit ;; (i.e., Jupiter too small for that many pixels). ?? Plus, getting same ;; results for both EUV and FUV. ;; January 23, 2004: Getting different results for different channels ;; (new Frames kernel). ;; Test comparison with CASPER VERY close to perfect agreement. Getting about ;; one more pixel intersecting Jupiter than is shown in CASPER plots, ;; but otherwise excellent agreement. Possible visual display error in CASPER? ;; Test case with CASPER numbers from A. Hendrix: 2001 Jan 1 11:10:15, ;; FUV row 38: ARH: Lat = 15N, Lon = 250W, GEOMTER: Lat =18N, Lon = 244W. ;; January 28, 2004: Use LT+S LT correction, and subtract LT time from ET for ;; body conversion. Discovered cause of extra pixel in earlier Jupiter tests: ;; faulty Jupiter PCK values used. Now have perfect agreement with CASPER ;; on Jupiter test of 2000 Dec 29 2:46. ;; February 2, 2004: Using 'LT' LT correction and fixed some errors got ;; excellent agreement between ingress and egress for Jupiter Sigma Leo ;; stellar occultation. Need to check if sometimes 'LT+S' should be used ;; instead of 'LT'. ;; February 3, 2004: Add ABERR keyword. ;; February 14, 2004: initialize latocc and lonocc with dummy values. ;; February 25, 2004: Add losdist. Fix lon computation when miss target. ;; March 16, 2004: Add CORNERS keyword and SLIT keyword: have agreement ;; between two corners of adjacent pixels. Full mapping test of slit looks ;; correct. ;; March 24, 2004: Allow for solar port. ;; April 1, 2004: add sc_wrt_plan, state_planet, and state_sun to output. ;; April 23, 2004: Prevent corner calculations for inertial pointing case. ;; April 26, 2004: Fix and simplify handling of off-body longitude values. ;; June 14-16, 2004: Fix bug in computation of losdist. Change corner values ;; of lat and lon to -1000 if off-body, instead of using occultation values. ;; Fix handling of Sun so that planet = 'Sun' if target = 'Sun' ;; June 21, 2004: was using light travel time to planet, not target.fixed that ;; June 28, 2004: Fixed bug in line of sight distance for ring target. ;; August 4, 2004: Set ckerror=0 when using star pointing. ;; September 16, 2004: Fixed bug in sub-solar longitude: was East longitude. ;; September 27, 2004: calculate longitude of ring intercept point in body ;; coordinate system (not ring coordinate system). ;; October 15, 2004: Calculate corners for HSP. ;; November 23, 2004: Change "boresight" pixel corners to the corners of ;; the entire field of view. ;; January 3, 2005: Add the MOON keyword to get ringplane intercept ;; longitude relative to a moon whose name is passed in through the MOON ;; keyword. ;; January 6, 2005: Change pixel calculation to work with new cas_v39.tf ;; frames kernel. Make default ring longitude relative to the intercept ;; of the Saturn equatorial plane on the J2000 X-Y plane. ;; April 12, 2005: Make zerolongitude adjustment for ring plane longitudes ;; when doing C kernel pointing. ;; April 27, 2005: Miodrag Sremcevic caught bug in ring longitude calculation ;; (when not referenced to Moon longitude). Include phase, emission, incidence ;; angles for lines of sight that miss the body, using the nearpoint as the ;; intercept. Algorithm from Kris Larsen. ;; April 29, 2005: Minor bug fix in checking ABERR keyword. ;; May 10, 2005: Minor cosmetic changes to programming. ;; May 23, 2005: replace !RADEG and !DTOR with double precision values. ;; Calculate ringoccphi. Remove TRANSPOSEs and replace with ## for SPICE ;; rotation matrices. Compute stellar aberration here instead of in ;; calling routine. Change ABERR keyword to ABCORR. Ring plane intercepts ;; verified by comparison with Dick French program. ;; June 15, 2005: add solar keyword for solar occs. Need to check light ;; travel time corrections throughout program. ;; October 24, 2005: Fixed losdist calculation for rays that miss target, ;; per Feng Tian bug catch. Fix value of rayheight inside body (case where ;; ray intercepts body). ;; November 10, 2005: Allow for "F RING" target. RINGS flag = 1 for main rings, ;; RINGS flag = 3 for F ring. Added f_peri output which is the longitude of ;; the pericenter of the F ring. True anomaly of the point on the F ring is ;; longitude-f_peri. ;; April 19, 2006: Changed missing target error handling per Dave Judd. ;; July 14, 2006: Include General Relativity correction (spherical target ;; body assumed; error is < 1 m for observations made by Cassini). ;; August 14, 2006: Fix bug in GR correction, and make it do one iteration. ;; September 17, 2006: No iteration in GR correction. Verified with RGF. ;; Use 3 degree square for HDAC (no recognized FOV selected). ;; December 5, 2006: Incorporate additional output requested by Don ;; Shemansky and implemented by Jean Yoshii at Space Environment. ;; January 30, 2007: Include CIRS footprints. ;; February 5, 2007: Fix sign error in F ring frame. ;; February 27, 2007: Include ISS footprints. ;; April 3, 2007: Fix bug for rayheights when looking away from minimum ;; line of sight distance. ;; June 11, 2007: Minor mod to be able to handle G Ring Arc frames kernel. ;; September 24, 2007: Add ramlat, ramlon calculation. ;; February 12, 2009: Add ringsolar output parameter. ;; March 25, 2009: Correct ringsolar bug. ;; April 22, 2009: add the SUBTYPE keyword. Add the DETIC keyword. If ;; the DETIC keyword is set, then latitudes and longitudes are ;; planetodetic, while otherwise they are planetocentric. ;;- COMMON basics, target_copy, et, lt_corr_targ, abcorr_copy, scname, radeg, dtor radeg = 180D0/!DPI dtor = !DPI/180D0 IF NOT KEYWORD_SET(DETIC) THEN detic = 0 ELSE detic = 1 ;;Make sure we have the right frames kernel loaded: CSPICE_GCPOOL, 'TEXT_KERNEL_ID', 0, 50, 120, vals, found ;CSPICE_GCPOOL,'TEXT_KERNEL_ID',0, vals, found properframeloaded = 0B FOR index = 0, N_ELEMENTS(vals)-1 DO BEGIN IF STRMID(vals[index], 0, 16) EQ 'CASSINI_FRAMES V' THEN $ IF FLOAT(STRMID(vals[index], 16, 3)) GE 3.9 THEN $ properframeloaded = 1B ENDFOR IF NOT properframeloaded THEN $ STOP,'You must load cas_v39.tf frames kernel or later to continue.' ;;Identify target body: upcasetarget = STRUPCASE(target) IF upcasetarget EQ "RINGS" OR upcasetarget EQ "F RING" THEN BEGIN target = 'SATURN' IF upcasetarget EQ "RINGS" THEN rings = 1B ELSE rings = 3B ENDIF ELSE rings = 0B target_copy = target CSPICE_BODN2C, target, targetid, found IF NOT found THEN BEGIN stopstring = 'Target '+target+$ ' not found. Stopping. Use .c to continue if you fix the problem.' MESSAGE, stopstring GOTO, skipc ENDIF IF STRUPCASE(target) EQ 'SUN' THEN planetid = 10L ELSE $ planetid = FIX(targetid/100)*100L+99L CSPICE_BODC2N, planetid, planet, found ;stop CSPICE_BODVAR, planetid, 'GM', gmplanet gmplanet = gmplanet[0] iref = 'J2000' scname = 'CASSINI' ;CSPICE_BODN2C, scname, sc_id, found sc_id = -82L IF NOT KEYWORD_SET(ABCORR) THEN abcorr = 'LT' abcorr_copy = abcorr IF NOT KEYWORD_SET(SUBTYPE) THEN subtype = 'nearpoint' framestring = 'IAU_'+target CSPICE_BODVAR, targetid, 'RADII', radii flat = (radii[0]-radii[2])/radii[0] origin = [0D, 0D, 0D] ringplanerad = -1000D0 ringoccphi = -1000D0 rayheight = -1000D0 rayradius = -1000D0 losdist = -1000D0 phase = -1000D0 emission = -1000D0 incidence = -1000D0 latocc = -1000D0 lonocc = -1000D0 lat = -1000D0 lon = -1000D0 f_peri = -1000D0 ramlat = -1000D0 ramlon = -1000D0 ringsolar = -1000D0 ringsolar1 = -1000D0; ETB 3/9/2010 ringsolar2 = -1000D0; ETB 3/9/2010 ringsolar3 = -1000D0; ETB 3/9/2010 ringsolar4 = -1000D0; ETB 3/9/2010 ram_sun_angle = -1000D0; ETB 6/14/2012 ;; values are non-physical in case no pointing information is supplied. lat1 = -1000D0 ;; these place-holder values for the corners will retain these values lat2 = -1000D0 ;; if the corners keyword is not set. lat3 = -1000D0 lat4 = -1000D0 lon1 = -1000D0 lon2 = -1000D0 lon3 = -1000D0 lon4 = -1000D0 ra1 = -1000D0 ra2 = -1000D0 ra3 = -1000D0 ra4 = -1000D0 dec1 = -1000D0 dec2 = -1000D0 dec3 = -1000D0 dec4 = -1000D0 ringplanerad1 = -1000D0 ringplanerad2 = -1000D0 ringplanerad3 = -1000D0 ringplanerad4 = -1000D0 inc1 = -1000D0 inc2 = -1000D0 inc3 = -1000D0 inc4 = -1000D0 ;;Get ephemeris time of observation: CSPICE_UTC2ET, time, et CSPICE_SCE2T, sc_id, et, ticks ;;Get the position of the spacecraft with respect to the target body: CSPICE_SPKEZR, target, et, iref, abcorr, scname, state, lt_corr_targ scloc = state[0:2] ;;Get state relative to planet (in case target is a satellite): CSPICE_SPKEZR, planet, et, iref, abcorr, scname, sc_wrt_plan, lt_corr_plan ;;Get planet position wrt solar system barycenter: CSPICE_SPKSSB, planetid, et, iref, state_planet IF target EQ 'SATURN' THEN BEGIN CSPICE_SPKEZR, 'SATURN', et, iref, abcorr, 'SUN', state_sat_sun, $ lt_corr_sat_sun CSPICE_PXFORM, iref, 'IAU_SATURN', et, inert2sat ramvec = REFORM(inert2sat##state_sat_sun[3:5]) IF detic THEN $ CSPICE_RECGEO, xpoint, radii[0], flat, lon, lat, deticoffset ELSE $ CSPICE_RECLAT, ramvec, junk, ramlon, ramlat ramlat *= !RADEG ramlon *= !RADEG ENDIF ELSE IF planetid EQ 699L THEN BEGIN targframe = 'IAU_'+target CSPICE_SPKEZR, target, et, iref, abcorr, 'SATURN', $ targstate_saturn, lt_corr_targsat CSPICE_PXFORM, iref, targframe, et, inert2targ ramvec = REFORM(inert2targ##targstate_saturn[3:5]) IF detic THEN $ CSPICE_RECGEO, xpoint, radii[0], flat, lon, lat, deticoffset ELSE $ CSPICE_RECLAT, ramvec, junk, ramlon, ramlat ramlat *= !RADEG ramlon *= !RADEG ENDIF ;CSPICE_PXFORM, iref, framestring, et-lt_corr_targ, inert2p CSPICE_PXFORM, iref, framestring, et, inert2p scloc_p = REFORM(-inert2p##scloc) ;;(see documentation of surfpt in naif toolkit.) ;;Get position of the Sun with respect to target in the target frame to get ;;the ring plane Sun angle: CSPICE_SPKEZR, 'SUN', et, framestring, abcorr, target, state_sun, lt_corr pos_sun =state_sun[0:2] sun_normal = pos_sun/CSPICE_VNORM(pos_sun) sun_elevation = ASIN(sun_normal[2])*RADEG CSPICE_SPKEZR, 'EARTH', et, framestring, abcorr, target, state_earth, lt_corr state_earth = state_earth[0:2] earth_normal = state_earth/CSPICE_VNORM(state_earth) earth_elevation = ASIN(earth_normal[2])*RADEG sunearth = CSPICE_VSEP(sun_normal, earth_normal) ;;Get S/C-Sun vector: CSPICE_SPKEZR, 'SUN', et, iref, abcorr, scname, sc_sun, lt_dist ;;Get boresight vector in inertial reference frame. IF KEYWORD_SET(PIXEL) THEN BEGIN IF pixel EQ -1 THEN px = 0 ELSE px = pixel CSPICE_AXISAR, [1L,0L,0L], (DOUBLE(px)-31.5)/1000D0, i2px Endif ELSE BEGIN i2px = [[1D0,0D0,0D0],[0D0,1D0,0D0],[0D0,0D0,1D0]] pixel=64 ENDELSE IF instrument NE 'NONE' THEN BEGIN IF KEYWORD_SET(CORNERS) AND KEYWORD_SET(SLIT) THEN BEGIN channel = STRMID(instrument, STRLEN(instrument)-3, STRLEN(instrument)-1) ;; This is confusing, and some day I will change it, but the units ;; in the CASE block below are mrad for width, and radians for hangle. CASE channel OF 'FP1': BEGIN ;; 0.22 degree circular FOV. hangle = 1.92E-3 ;;radians width = 3.84 ;;mrad END 'FP3': BEGIN hangle = 1.48E-3 ;;radians ;;"height" from center of pixel to top or bottom edge of pixel. width = .35 ;;mrad END 'FP4': BEGIN hangle = 1.48 ;;radians ;;"height" from center of pixel to top or bottom edge of pixel. width = .35 ;;mrad END 'NAC': BEGIN hangle = 0.003054 ;; radians of half the height of the NAC FOV width = 6.1 ;; mrad full width of NAC FOV END 'WAC': BEGIN hangle = 0.02967 ;; radians of half the height of the WAC FOV width = 60.74 ;; mrad full width of WAC FOV END 'FUV': BEGIN ;;If the pixel is the boresight pixel, then use the full ;;slit length rather than a single pixel height. hangle = 5D-4*(64.*((pixel EQ 64) ? 1. : 1./64.)) ;;"height" from center of pixel to top or bottom edge of pixel. CASE slit OF 'LO':width = 1.5 ;;mrad 'HI':width = 0.75 'OCC':width = 8. ENDCASE END 'EUV': BEGIN ;;If the pixel is the boresight pixel, then use the full ;;slit length rather than a single pixel height. hangle = 5D-4*(64.*((pixel EQ 64) ? 1. : 1./64.)) ;;"height" from center of pixel to top or bottom edge of pixel. CASE slit OF 'LO':width = 2. ;;mrad 'HI':width = 1. 'OCC':width = 8. ENDCASE END 'LAR': BEGIN width = 8. ;;last three letters of the 'SOLAR' instrument. ;;If the pixel is the boresight pixel, then use the full ;;slit length rather than a single pixel height. hangle = 5D-4*(64.*((pixel EQ 64) ? 1. : 1./64.)) END 'HSP': BEGIN width = 6. ;;mrad hangle = 3D0/1D3 ;;radians END ELSE: BEGIN width = 52. hangle = 0.052 END ENDCASE ;;Get the rotation matrices that convert the vector along the pixel ;;center to the corners of the pixel: wangle = 0.5*width/1D3 ;;half of the width of the pixel. CSPICE_AXISAR, [0L, 1L, 0L], wangle, w_pos_mat CSPICE_AXISAR, [0L, 1L, 0L], -wangle, w_neg_mat CSPICE_AXISAR, [1L, 0L, 0L], hangle, h_pos_mat CSPICE_AXISAR, [1L, 0L, 0L], -hangle, h_neg_mat corner1mat = w_pos_mat#h_pos_mat corner2mat = w_neg_mat#h_pos_mat corner3mat = w_neg_mat#h_neg_mat corner4mat = w_pos_mat#h_neg_mat ENDIF ;;Corners and Slit set. IF ra EQ -1000D0 THEN BEGIN ;;This is the case where we get pointing from a C kernel. inertial_pointing = 0B ;stop ckerror = 0B CATCH, /CANCEL CATCH, pxform_error_status IF STRMID(!ERROR_STATE.MSG,0,13) EQ 'CSPICE_PXFORM' THEN BEGIN ;;Don't have data for pointing at that time. Set a flag and return. ckerror = 1B CATCH, /CANCEL MESSAGE, /RESET_ERROR_STATE GOTO, skipc ENDIF CSPICE_PXFORM, instrument, iref, et, i2inert boresight=i2px#[0D,0D,1D] ;;boresight in instrument frame boresight_inert = REFORM(i2inert##boresight) CSPICE_RECRAD, boresight_inert, radius, ra, dec Ra_save = ra*RADEG dec_save = dec*RADEG ;;Need the rotation matrix from instrument to target frame. CSPICE_PXFORM, instrument, framestring, et-lt_corr_targ, i2p boresight_p = REFORM(i2p ## boresight) IF KEYWORD_SET(CORNERS) THEN BEGIN CORNER_POINTING_temp, corner1mat, i2inert, inert2p, scloc_p, losdist, boresight, radii, ra1, dec1, lat1, lon1, inc1, detic CORNER_POINTING_temp, corner2mat, i2inert, inert2p, scloc_p, losdist, $ boresight, radii, ra2, dec2, lat2, lon2, inc2, detic CORNER_POINTING_temp, corner3mat, i2inert, inert2p, scloc_p, losdist, $ boresight, radii, ra3, dec3, lat3, lon3, inc3, detic CORNER_POINTING_temp, corner4mat, i2inert, inert2p, scloc_p, losdist, $ boresight, radii, ra4, dec4, lat4, lon4, inc4, detic ENDIF ENDIF ELSE BEGIN ;;This is for the case where the pointing is specified in inertial space, ;;and is not obtained from a C kernel, so computing corners for the FOV ;;makes no sense. inertial_pointing = 1B ;; Save values passed in by calling routine so they will not be altered. ra_save = ra dec_save = dec ra *= DTOR dec *= DTOR IF KEYWORD_SET(SOLAR) THEN BEGIN ;;We're looking at the Sun. Get the RA/Dec of the Sun at the ;;time solar photons reached the target body (check this!) CSPICE_SPKPOS, 'SUN', et-lt_corr_targ, iref, abcorr, scname, $ sunrapos, lt_corr CSPICE_RECRAD, sunrapos, unit, ra, dec ENDIF ;; Correction for stellar aberration due to motion of spacecraft: CSPICE_SPKEZR, 'SATURN', et, iref, 'NONE', 'SSB', satstatessb, lt_corr_temp CSPICE_RADREC, 1D0, ra, dec, boresight_inert clight = CSPICE_CLIGHT() boresight_inert += satstatessb[3:5]/clight boresight_inert /= CSPICE_VNORM(boresight_inert) ;;Correct for General Relativity change in apparent star position: CSPICE_SPKPOS, planet, et, iref, abcorr, scname, scloc_i, lt_corr CSPICE_NPEDLN, radii[0], radii[1], radii[2], scloc_i, $ boresight_inert, nearpoint, rayheight nearpointmag = CSPICE_VNORM(nearpoint) neardist = nearpointmag+rayheight impact_param = -nearpoint/nearpointmag*neardist grfactor = 4D0*gmplanet/(clight*clight) ;; Compute deflection assuming a spherical planet. Although this ;; is generally not a good assumption, the magnitude of the change ;; for spacecraft based occultations is much smaller than the measurement ;; resolution. Joe Spitale calculates 17 cm change due to J2 for the ;; 2004 Xi(2) Ceti occultation, which is the occultation at the greatest ;; distance from Saturn. boresight_inert += grfactor/(neardist*neardist)*impact_param boresight_inert /= CSPICE_VNORM(boresight_inert) ;; Get new aberrated/GR-corrected RA and Dec for the rest of the program. CSPICE_RECRAD, boresight_inert, unit, ra, dec Boresight_p = REFORM(inert2p ## boresight_inert) ckerror = 0B ENDELSE CSPICE_SURFPT, scloc_p, boresight_p, radii[0], radii[1], radii[2], $ xpoint, found IF found THEN BEGIN CSPICE_ILLUM, target, et-lt_corr_targ, abcorr, scname, xpoint, $ phase, incidence, emission phase *= RADEG incidence *= RADEG emission *= RADEG ;stop ;;Get latitude of xpoint: IF detic THEN $ CSPICE_RECGEO, xpoint, radii[0], flat, lon, lat, deticoffset ELSE $ CSPICE_RECLAT, xpoint, bodintercept, lon, lat losdist = CSPICE_VNORM(scloc_p-xpoint) CSPICE_NPLNPT, scloc_p, boresight_p, origin, nearpoint, rayradius ;;the above call to nplnpt finds the distance, rayradius, ;;between the planet origin and the line containing ;;the spacecraft position and the vector along the line of sight. ;;Next, get where the vector from planet center in the ;;direction of the minimum ray height pierces the body: CSPICE_SURFPT, origin, 1000D0*nearpoint, radii[0], radii[1], $ radii[2], xpoint, found IF NOT found THEN STOP,'Vector did not reach body surface.' ;;Now get impact parameter: rayheight = rayradius - CSPICE_VNORM(xpoint) ;; Handle the case where we're looking away from the nearpoint. In ;; this case, rayheight should be the altitude and rayradius should ;; be the distance from body center: IF CSPICE_VSEP(boresight_p, scloc_p) LT !DPI/2D0 THEN BEGIN rayradius = CSPICE_VNORM(scloc_p) CSPICE_SUBPT, subtype, target, et, abcorr, scname, point, $ rayheight ENDIF ENDIF ELSE BEGIN ;;handle case where it misses CSPICE_NPEDLN, radii[0], radii[1], radii[2], scloc_p, $ boresight_p, nearpoint, rayheight rayradius = CSPICE_VNORM(nearpoint)+rayheight ;; Get where the vector from planet center to the minimum ;; ray height pierces the body: CSPICE_SURFPT, origin, nearpoint, radii[0], radii[1], $ radii[2], xpoint, found losdist = SQRT(CSPICE_VNORM(scloc_p)^2-rayradius^2) ;; Handle the case where we're looking away from the nearpoint. In ;; this case, rayheight should be the altitude and rayradius should ;; be the distance from body center, and losdist should be the ;; altitude: IF CSPICE_VSEP(boresight_p, scloc_p) LT !DPI/2D0 THEN BEGIN rayradius = CSPICE_VNORM(scloc_p) CSPICE_SUBPT, subtype, target, et, abcorr, scname, point, $ rayheight losdist = rayheight ENDIF ENDELSE ;; Calculate illumination angles for lines of sight that do not ;; intersect the body, using the nearest point on the line of sight ;; to the body. IF phase EQ -1000 THEN BEGIN CSPICE_SURFNM, radii[0], radii[1], radii[2], nearpoint, out_norm out_norm *= rayheight np1 = nearpoint + out_norm CSPICE_ILLUM, target, et-lt_corr_targ, abcorr, scname, np1, $ phase, incidence, emission phase *= RADEG incidence *= RADEG emission *= RADEG ;stop ENDIF ;; Get latitude of nearpoint: IF detic THEN $ CSPICE_RECGEO, xpoint, radii[0], flat, lon, lat, deticoffset ELSE $ CSPICE_RECLAT, nearpoint, bodyintercept, lonocc, latocc ;;This is the latitude and longitude of the point on the ;;body below the minimum ray height. These are not the same ;;as lat and lon. ra *= RADEG dec *= RADEG lat = lat*RADEG > (-1000.0) IF lon NE -1000. THEN lon = (2*!DPI-lon)*RADEG MOD 360D0 latocc *= RADEG lonocc = (2*!DPI-lonocc)*RADEG MOD 360D0 IF rings THEN BEGIN ;; Get the Saturn pole vector. sat_pole_body = [0D0, 0D0, 1D0] CSPICE_MTXV, inert2p, sat_pole_body, pole CSPICE_RECRAD, pole, unitvec, pole_ra, pole_dec ;; If it's the F ring, get the pole of the plane of the F ring ;; instead: IF rings EQ 3B THEN BEGIN ;; Bosh et al. (2002) F Ring model uses the following Saturn ;; pole. Need to use it to have the appropriate F Ring pole. ; nyears = et/3.15576D7 ; sat_pole_ra = 40.59287D0*DTOR - 6.1772D-4*DTOR*nyears ; sat_pole_dec = 83.53833D0*DTOR -6.42D-5*DTOR*nyears ;; Get a three-vector for the Saturn pole: ; CSPICE_RADREC, 1D0, sat_pole_ra, sat_pole_dec, sat_pole ;; Or use NAIF values for consistency with other ring solutions: sat_pole_ra = pole_ra sat_pole_dec = pole_dec sat_pole = pole ;; Bosh et al. (2002) Table III, Fit #3: fring_inc = 6.5D-3*DTOR fring_node = (16.1*DTOR - 2.6876*DTOR*et/86400D0) MOD (2*!DPI) IF fring_node LT 0D0 THEN fring_node += (2*!DPI) f_node = fring_node*RADEG ;; fring_node is the angle from the Saturn ascending node to ;; the F Ring ascending node. fring_pericenter = (24.1*DTOR + 2.7001*DTOR*et/86400D0) MOD $ (2*!DPI) f_peri = fring_pericenter*RADEG ;; fring_pericenter is the angle from the Saturn ascending node ;; to the pericenter of the F ring. ;; Get a three-vector that points to the ascending node ;; of Saturn: sat_node_vec = [-sat_pole[1], sat_pole[0], 0D0] CSPICE_TWOVEC, sat_pole, 3, sat_node_vec, 1, inert2sat_node ;; inert2sat_node rotates from J2000 to a Saturn frame where ;; the X axis is along the ascending node of Saturn. ;; Now get the matrix to go from this Saturn frame to the F ;; ring frame. The x axis of the F ring frame points to the ;; ascending node of the F ring. CSPICE_EUL2M, -fring_node, -fring_inc, 0D0, 3L, 1L, 3L, sat2f ;; minus signs added in above line per bug catch by ;; Nicole Albers, January 2007. CSPICE_MTXM, sat2f, inert2sat_node, inert2f CSPICE_MTXV, inert2f, sat_pole_body, pole CSPICE_RECRAD, pole, unitvec, pole_ra, pole_dec ;; pole is the pole of the F ring in the inertial frame, and ;; pole_ra and pole_dec are the RA and Dec of the F ring pole. ENDIF ;; Convert the pole vector to a CSPICE plane using nvc2pl. ;; Define origin of our system at center of Saturn. CSPICE_NVC2PL, pole, 0D0, ringplane ;;Get ring plane intercept radius. ;;Get a vector from the planet center to the spacecraft: CSPICE_SPKPOS, target, et, iref, abcorr, $ scname, scloc_i, lt_corr scloc_i = -scloc_i ;;This returns occpoint_i, a vector from Saturn center to the ;;ring plane intercept point in the inertial frame. CSPICE_INRYPL, scloc_i, boresight_inert, ringplane, nxpts, occpoint_i ;; Get the solar ring hour angle of the intercept point: ringsolar = CSPICE_VSEP(occpoint_i, state_sat_sun[0:2])*!RADEG vtemp = CROSSP(occpoint_i, state_sat_sun[0:2]) IF vtemp[2] LT 0. THEN ringsolar = 360.-ringsolar; ETB 3/9/2010 (changed GT to LT) ;; Get the intersection of the Saturn equatorial plane ;; with the J2000 X-Y plane to get an inertial reference ;; direction for zero longitude for the rings: xzerolon = CROSSP([0D0,0D0,1D0], pole) xzerolon = xzerolon/CSPICE_VNORM(xzerolon) IF xzerolon[1] GT 0. THEN $ zerolon = ACOS(xzerolon[0])*RADEG ELSE $ zerolon = 360D0-ACOS(xzerolon[0])*RADEG IF NXPTS EQ 1 THEN BEGIN ringplanerad = CSPICE_VNORM(occpoint_i) ;print,pole_ra*radeg,format = '(f20.14)' ;print,pole_dec*radeg,format = '(f20.14)' ;;French's algorithm: ;cspice_radrec,1d0,ra_save*dtor,dec_save*dtor, nhat_star ;cspice_spkssb,-82l,et,iref,r_e_t0_state ;r_e_t0 = r_e_t0_state[0:2] ;cspice_spkssb,targetid, et, iref, r_s_state ;r_s = r_s_state[0:2] ;rdot_s = r_s_state[3:5] ;nhat_starp = nhat_star+rdot_s/cspice_clight() ;nhat_starp /=cspice_vnorm(nhat_starp) ;r_sc = r_e_t0-r_s ;delta2 = cspice_vdot(-r_sc,pole)/cspice_vdot(nhat_starp,pole) ;r_si_pc = r_sc+delta2*nhat_starp ;print,'RGF minus JEC: ',cspice_vnorm(r_si_pc)-ringplanerad CSPICE_TIPBOD, iref, targetid, et, i2targ occpoint_p = i2targ##occpoint_i ;; Calculate clock angle phi of occultation (angle between ;; star look vector projected into ring plane and occpoint ;; vector): occp_3vec = [occpoint_p[0], occpoint_p[1], 0D0] occp_3vec /= CSPICE_VNORM(occp_3vec) b_p_proj = [boresight_p[0], boresight_p[1], 0D0] ; b_p_proj = [boresight_p[0], boresight_p[1], boresight_p[2]] b_p_proj /= CSPICE_VNORM(b_p_proj) vcross = CROSSP(occp_3vec, b_p_proj) ringoccphi = CSPICE_VSEP(occp_3vec, b_p_proj)*RADEG IF vcross[2] LT 0. THEN ringoccphi = 360D0-ringoccphi IF KEYWORD_SET(MOON) THEN BEGIN CSPICE_SPKPOS, STRING(moon), et, iref, abcorr, target, $ moon_i, lt_corr occ_norm = REFORM(occpoint_p/ringplanerad) moon_p = REFORM(i2targ##moon_i) moon_norm = moon_p/CSPICE_VNORM(moon_p) lon = CSPICE_VSEP(occ_norm,moon_norm)*RADEG IF TOTAL(pole*CROSSP(moon_norm,occ_norm)) LT 0D0 THEN $ lon = 360D0-lon ; Previous two lines from Miodrag Sremcevic are a more ; elegant solution to the longitude than what I originally ; had which are the commented lines below, which give the ; same result. ; vcross = CROSSP(moon_norm, occ_norm) ; sinlon = CSPICE_VNORM(vcross) ; signsinlon = vcross[2]/ABS(vcross[2]) ; coslon = TRANSPOSE(moon_norm) # occ_norm ; IF signsinlon GT 0. THEN $ ; lon = ACOS(coslon)*!RADEG ELSE $ ; lon = 360D0-ACOS(coslon)*!RADEG ENDIF ELSE BEGIN lon = CSPICE_VSEP(occpoint_i,xzerolon)*RADEG IF TOTAL(pole*CROSSP(xzerolon,occpoint_i)) LT 0D THEN $ lon = 360D0-lon ENDELSE ;;Get lighting angles for boresight in the ring plane. Note that ;;these overwrite the body lighting angles calculated above if ;;it hit the body because in this loop the RING keyword has ;;been set. That overwrites body values. scfromring = scloc_i-occpoint_i losdist = CSPICE_VNORM(scfromring) CSPICE_SPKPOS, 'SUN', et, iref, abcorr, scname, sctosun, lt_corr sunfromring = scfromring+sctosun phase = CSPICE_VSEP(scfromring, sunfromring)*RADEG emission = CSPICE_VSEP(scfromring, pole)*RADEG incidence = CSPICE_VSEP(sunfromring, pole)*RADEG ;; Return original RA and Dec values to calling routine. ra = ra_save dec = dec_save ENDIF IF KEYWORD_SET(CORNERS) AND NOT inertial_pointing THEN BEGIN CSPICE_RADREC, 1D0, ra1*DTOR, dec1*DTOR, corner1_inert CSPICE_RADREC, 1D0, ra2*DTOR, dec2*DTOR, corner2_inert CSPICE_RADREC, 1D0, ra3*DTOR, dec3*DTOR, corner3_inert CSPICE_RADREC, 1D0, ra4*DTOR, dec4*DTOR, corner4_inert CSPICE_INRYPL, scloc_i, corner1_inert, ringplane, nxpts1, $ occpoint1_i ;;JOSH HERE! What about relative to moon longitudes for the corners? IF nxpts1 THEN BEGIN ringplanerad1 = CSPICE_VNORM(occpoint1_i) coslon = occpoint1_i[0]/ringplanerad1 sinlon = occpoint1_i[1]/ringplanerad1 IF sinlon GT 0. THEN $ lon1 = ACOS(coslon)*RADEG ELSE $ lon1 = (!DPI*2D0-ACOS(coslon))*RADEG lon1 -= zerolon ringsolar1 = CSPICE_VSEP(occpoint1_i, state_sat_sun[0:2])*!RADEG; ETB 3/9/2010 vtemp = CROSSP(occpoint1_i, state_sat_sun[0:2]) IF vtemp[2] LT 0. THEN ringsolar1 = 360.-ringsolar1 ENDIF CSPICE_INRYPL, scloc_i, corner2_inert, ringplane, nxpts2, $ occpoint2_i IF nxpts2 THEN BEGIN ringplanerad2 = CSPICE_VNORM(occpoint2_i) coslon = occpoint2_i[0]/ringplanerad2 sinlon = occpoint2_i[1]/ringplanerad2 IF sinlon GT 0. THEN $ lon2 = ACOS(coslon)*RADEG ELSE $ lon2 = (!DPI*2D0-ACOS(coslon))*RADEG lon2 -= zerolon ringsolar2 = CSPICE_VSEP(occpoint2_i, state_sat_sun[0:2])*!RADEG; ETB 3/9/2010 vtemp = CROSSP(occpoint2_i, state_sat_sun[0:2]) IF vtemp[2] LT 0. THEN ringsolar2 = 360.-ringsolar2 ENDIF CSPICE_INRYPL, scloc_i, corner3_inert, ringplane, nxpts3, $ occpoint3_i IF nxpts3 THEN BEGIN ringplanerad3 = CSPICE_VNORM(occpoint3_i) coslon = occpoint3_i[0]/ringplanerad3 sinlon = occpoint3_i[1]/ringplanerad3 IF sinlon GT 0. THEN $ lon3 = ACOS(coslon)*RADEG ELSE $ lon3 = (!DPI*2D0-ACOS(coslon))*RADEG lon3 -= zerolon ringsolar3 = CSPICE_VSEP(occpoint3_i, state_sat_sun[0:2])*!RADEG; ETB 3/9/2009 vtemp = CROSSP(occpoint3_i, state_sat_sun[0:2]) IF vtemp[2] LT 0. THEN ringsolar3 = 360.-ringsolar3 ENDIF CSPICE_INRYPL, scloc_i, corner4_inert, ringplane, nxpts4, $ occpoint4_i IF nxpts4 THEN BEGIN ringplanerad4 = CSPICE_VNORM(occpoint4_i) coslon = occpoint4_i[0]/ringplanerad4 sinlon = occpoint4_i[1]/ringplanerad4 IF sinlon GT 0. THEN $ lon4 = ACOS(coslon)*RADEG ELSE $ lon4 = (!DPI*2D0-ACOS(coslon))*RADEG lon4 -= zerolon ringsolar4 = CSPICE_VSEP(occpoint4_i, state_sat_sun[0:2])*!RADEG; ETB 3/9/2010 vtemp = CROSSP(occpoint4_i, state_sat_sun[0:2]) IF vtemp[2] LT 0. THEN ringsolar4 = 360.-ringsolar4 ENDIF ENDIF ENDIF ENDIF ;;No instrument. skipc: ;;label. jump here if we had a problem with the C kernel. ;;Now get the basic geometery data based just on position of spacecraft ;;and target: ;;Already have position and velocity (scloc and scvel) above in planet frame. ;;Need them in inertial frame: CSPICE_SPKEZR, target, et, iref, abcorr, scname, state_i, lt_corr scloc = state_i[0:2] scvel = state_i[3:5] CSPICE_RECRAD, scloc, dist_targ, ra_targ, dec_targ ra_targ *= RADEG dec_targ *= RADEG CSPICE_SUBPT, subtype, target, et, abcorr, scname, point, alt IF detic THEN $ CSPICE_RECGEO, xpoint, radii[0], flat, lon, lat, deticoffset ELSE $ CSPICE_RECLAT, point, radius, lon_subpt, lat_subpt lon_subpt = (2*!DPI-lon_subpt)*RADEG MOD 360D0 lat_subpt *= RADEG CSPICE_SUBSOL, subtype, target, et, abcorr, scname, sunpoint IF detic THEN $ CSPICE_RECGEO, xpoint, radii[0], flat, lon, lat, deticoffset ELSE $ CSPICE_RECLAT, sunpoint, radius, lon_sunpt, lat_sunpt lon_sunpt = (2*!DPI-lon_sunpt)*RADEG MOD 360D0 lat_sunpt *= RADEG ram_sun_angle = cspice_vsep(sunpoint,ramvec)*180./!pi; ETB 6/14/2012 CSPICE_ILLUM, target, et-lt_corr_targ, abcorr, scname, point, $ phase_targ, incidence_targ, emission_targ phase_targ *= RADEG incidence_targ *= RADEG emission_targ *= RADEG ra_full = [ra, ra1, ra2, ra3, ra4] dec_full = [dec, dec1, dec2, dec3, dec4] lat_full = [lat, lat1, lat2, lat3, lat4] lon_full = [lon, lon1, lon2, lon3, lon4] rrad_full = [ringplanerad, ringplanerad1, ringplanerad2, ringplanerad3, $ ringplanerad4] inc_full = [incidence, inc1, inc2, inc3, inc4] ringsolar_full=[ringsolar,ringsolar1,ringsolar2,ringsolar3,ringsolar4]; ETB 3/9/2010 CASE rings OF 0B: 1B: target = 'RINGS' 3B: target = 'F RING' ENDCASE ;;Gotta reset this parameter for the calling routine. ;STOP,'Stopping normally at end of geometer_engine.' END PRO CORNER_POINTING_temp, cornermatrix, i2inert, inert2p, scloc_p, losdist, boresight, radii, ra, dec, lat, lon, incidence, detic COMMON basics, target_copy, et, lt_corr_targ, abcorr_copy, scname, radeg, dtor cornervec = cornermatrix#boresight corner_inert = REFORM(i2inert##cornervec) CSPICE_RECRAD, corner_inert, radius, ra, dec ra *= RADEG dec *= RADEG corner_p = REFORM(inert2p##corner_inert) CSPICE_SURFPT, scloc_p, corner_p, radii[0], radii[1], radii[2], xpoint, found IF found THEN BEGIN IF detic THEN $ CSPICE_RECGEO, xpoint, radii[0], flat, lon, lat, deticoffset ELSE $ CSPICE_RECLAT, xpoint, bodyintercept, lon, lat losdist = CSPICE_VNORM(scloc_p-xpoint) CSPICE_ILLUM, target_copy, et-lt_corr_targ, abcorr_copy, scname, $ xpoint, phase, incidence, emission incidence *= RADEG lat *= RADEG lon = (2*!DPI-lon)*RADEG MOD 360D0 ENDIF ELSE BEGIN CSPICE_NPEDLN, radii[0], radii[1], radii[2], scloc_p, corner_p, $ nearpoint, rayheight lat = -1000.0 lon = -1000.0 incidence = -1000.0 losdist = SQRT(CSPICE_VNORM(scloc_p)^2-(CSPICE_VNORM(nearpoint)+$ rayheight)^2) ENDELSE END