;+
; NAME:
;        GET_ASPECTS
;
; PURPOSE:
;        Subroutine interpolates, or extrapolates, the aspect angles, right
;        ascension and declination (alpha, delta), for the Y- and Z-axes
;        given the initial and final value for these angles.  The Z-axis is
;        assumed to be the spin axis.
;
;        The interpolation (or extrapolation) is done by first converting
;        alpha,delta to celestial cartesian coordinates, then cartesian
;        coordinates to Euler angles.  The Euler angles are interpolated to
;        the desired number of intermediate points, then converted back to
;        cartesian coordinates, and finally to alphas and deltas.
;
;        Note that the aspects are interpolated, or extrapolated, to the
;        middle of each bin, and that the initial and final aspects
;        supplied in the input arguments bound the beginning of the first
;        bin and the end of the last bin, respectively.
;
;
; CATEGORY:
;        Interpolation
;
; CALLING SEQUENCE:
;
;        GET_ASPECTS, Mode, Bins, Ry,Dy, Rz,Dz, Y_asp, Z_asp
;
; INPUTS:
;         Mode:    The interpolation "mode", [integer].
;             If MODE=0, aspects are INTERPOLATED for values between
;             the initial and final aspects.
;             If MODE=1, aspects are EXTRAPOLATED beyond the final
;             aspects given.
;
;         Bins:    The number of intermediate or extrapolated bins
;             between the initial and final aspects, [integer].
;
;           Ry:    The initial and final values of the Y-axis
;             right ascension, [float(2)].
;
;           Dy:    The initial and final values of the Y-axis
;             declination, [float(2)].
;
;           Rz:    The initial and final values of the Z-axis
;             right ascension, [float(2)].
;
;           Dz:    The initial and final values of the Z-axis
;             declination, [float(2)].
;
; OPTIONAL INPUT KEYWORDS:
;       OFFSET:    The number of bins from the bin edge where the aspects
;             are calculated.  The default is 0.5 bin, namely, the center
;             of each bin, [float].
;
; OUTPUTS:
;        Y_asp:    The interpolated (or extrapolated) Y-axis right
;             ascension and declination for each of BINS points, [float(2,Bins)].
;
;        Z_asp:    The interpolated (or extrapolated) Z-axis right
;             ascension and declination for each of BINS points, [float(2,Bins)].
;
; MODIFICATION HISTORY:
;        Written by:    K.H. Fairfield,  02-DEC-1992.
;        09-FEB-93:     Changed the meaning of the BINS input argument to be
;                  the actual number of bins between the bounding aspects,
;                  rather than the number of intermediate points, so as to
;                  calculate the returned aspects at the center of each bin
;                  rather than the bin edge.
;
;                  N.B. Because of these changes, calling programs must be
;                       modified accordingly.
;        06-JUN-1994:   Adapted routine into an IDL procedure, H.C. Wen.
;        22-JUN-1994:   Vectorized routine to eliminate for loops.
;                       Added the OFFSET keyword.
;        15-AUG-1995    Comment bugfix: removed extraneous ;+ and ;-.
;-
pro GET_ASPECTS, Mode, Bins, Ry,Dy, Rz,Dz, Y_asp, Z_asp, OFFSET=Offset

         Y_asp = fltarr(2,Bins)
         Z_asp = fltarr(2,Bins)

;
;  Get the X-, Y- and Z-axis unit vectors for the initial aspect.
;
         slaCS2C, Ry(0),Dy(0), Yvec
         slaCS2C, Rz(0),Dz(0), Zvec
         slaVXV,  Yvec, Zvec, Xvec

;
;  Form the rotation matrix from the X-, Y-, and Z-axis direction cosines.
;
         rotmat = fltarr(3,3)
         rotmat(0,*) = xvec
         rotmat(1,*) = yvec
         rotmat(2,*) = zvec
;
;  Calculate the Euler angles for this rotation matrix.
;
         MTOEA, rotmat, theta_0, phi_0, psi_0

;
;  Get the X-, Y- and Z-axis unit vectors for the final aspect.
;
         slaCS2C, Ry(1),Dy(1), Yvec
         slaCS2C, Rz(1),Dz(1), Zvec
         slaVXV,  Yvec, Zvec, Xvec

;
;  Form the rotation matrix from the X-, Y-, and Z-axis direction cosines.
;
         rotmat(0,*) = xvec
         rotmat(1,*) = yvec
         rotmat(2,*) = zvec
;
;  Calculate the Euler angles for this rotation matrix.
;
         MTOEA, rotmat, theta_1, phi_1, psi_1
         IF (psi_1 LT 0.0) AND (psi_0 GT 0.0) then $
          psi_1 = psi_1 + 2.*!dpi

;
;  Calculate the delta-Euler angles.
;
         dth = (theta_1 - theta_0) / bins
         dph = (phi_1 - phi_0  ) / bins
         dps = (psi_1 - psi_0  ) / bins

;
;  Loop over then number of bins, interpolating/extroplating to the
;  "next" set of Euler angles, and for each set, convert back to rotation
;  matrix, then to aspect angles.  Initial theta, phi and psi are adjusted
;  to the middle of the preceding bin.
;
         if not keyword_set( OFFSET ) then Offset = 0.5
         IF (Mode EQ 0) then begin
           th1 = theta_0 - (Offset * dth)
           ph1 = phi_0   - (Offset * dph)
           ps1 = psi_0   - (Offset * dps)
         endif else begin
           th1 = theta_1 - (Offset * dth)
           ph1 = phi_1   - (Offset * dph)
           ps1 = psi_1   - (Offset * dps)
         endelse

         range = findgen(Bins)+1.0
         EATOM, th1+range*dth, ph1+range*dph, ps1+range*dps, rotmat

         yvec = REFORM( rotmat(1,*,*) )
         zvec = REFORM( rotmat(2,*,*) )
         slaCC2S, yvec, y_asp_RA, y_asp_DEC
         slaCC2S, zvec, z_asp_RA, z_asp_DEC

         Y_asp(0,*) = y_asp_RA
         Y_asp(1,*) = y_asp_DEC

         Z_asp(0,*) = z_asp_RA
         Z_asp(1,*) = z_asp_DEC

end