C
C     This program computes a few NEAR-Eros geometry parameters that
C     Tony Farnham used to verify the repackaged NEAR SPICE
C     archive by a check against Horizons.
C
C     BVS/NAIF, 02/06/13.
C
      PROGRAM CHECK_NEAR

      IMPLICIT NONE

C
C     Variables
C
      CHARACTER*(256)       LINE
      CHARACTER*(256)       MK
      CHARACTER*(256)       UTC
      CHARACTER*(256)       TRGFRM

      CHARACTER*(32)        TRGNAM
      INTEGER               TRGID
      CHARACTER*(32)        SCNAM
      CHARACTER*(32)        INFRM
      CHARACTER*(32)        ABCORR

      DOUBLE PRECISION      ET

      INTEGER               N
      DOUBLE PRECISION      RADII  ( 3 )
      DOUBLE PRECISION      RADSPH ( 3 )
      DOUBLE PRECISION      RE
      DOUBLE PRECISION      RP
      DOUBLE PRECISION      F

      DOUBLE PRECISION      ALT

      DOUBLE PRECISION      EPOS   ( 3 )
      DOUBLE PRECISION      LT

      DOUBLE PRECISION      RANGE
      DOUBLE PRECISION      RA
      DOUBLE PRECISION      DEC

      DOUBLE PRECISION      MAT    ( 3, 3 )
      DOUBLE PRECISION      Z      ( 3 )
      DOUBLE PRECISION      POLE   ( 3 )

      DOUBLE PRECISION      SPOINT ( 3 )
      DOUBLE PRECISION      TRGEPC ( 3 )
      DOUBLE PRECISION      SRFVEC ( 3 )
      DOUBLE PRECISION      R
      DOUBLE PRECISION      LON
      DOUBLE PRECISION      LAT
      
      DOUBLE PRECISION      NORMAL ( 3 )

      DOUBLE PRECISION      SUBSC  ( 3 )

      DOUBLE PRECISION      PHASE
      DOUBLE PRECISION      SOLAR
      DOUBLE PRECISION      EMISSN

C
C     SPICELIB functions.
C
      DOUBLE PRECISION      DPR
      DOUBLE PRECISION      VSEP
      DOUBLE PRECISION      TWOPI

      INTEGER               PLNSNS

C
C     Set mission specific items.
C
      TRGNAM = 'EROS'
      TRGID  = 2000433
      SCNAM  = 'NEAR'
      INFRM  = 'J2000'
      ABCORR = 'LT+S'
      
C
C     Process command line.
C
      CALL GETCML( LINE )
      IF ( LINE .EQ. ' ' ) THEN
         WRITE (*,*) 'Usage: check_near <mk> <utc> <eros_frame>'
         STOP
      END IF

      CALL NEXTWD ( LINE, MK,  LINE )
      CALL NEXTWD ( LINE, UTC, TRGFRM )

      CALL FURNSH ( MK )

      CALL STR2ET ( UTC, ET )

C
C     Get target RADII. Make a copy of RADII for spheroidal shape.
C     Compute flattening coefficient.
C
      CALL BODVRD ( TRGNAM, 'RADII', 3, N, RADII )
      CALL VPACK  ( RADII(1), RADII(1), RADII(3), RADSPH ) 

      RE  =  RADII(1)
      RP  =  RADII(3)
      F   =  ( RE - RP ) / RE

C
C     Compute ra/dec/range of target as seen from spacecraft in
C     inertial frame.
C
      CALL SPKPOS ( TRGNAM, ET, INFRM, ABCORR, SCNAM, EPOS, LT )

      CALL RECRAD ( EPOS, RANGE, RA, DEC )
      WRITE (*,*) 
      WRITE (*,'(A20,F10.2)') 'Target RA    ', RA*DPR()
      WRITE (*,'(A20,F10.2)') 'Target Dec   ', DEC*DPR()
      WRITE (*,'(A20,F10.2)') 'Target Range ', RANGE

C
C     Compute target pole RA and DEC in inertial frame.
C
      CALL PXFORM ( TRGFRM, INFRM, ET, MAT ) 
      CALL VPACK  ( 0.D0, 0.D0, 1.D0, Z )
      CALL MXV    ( MAT, Z, POLE )

      CALL RECRAD ( POLE, RANGE, RA, DEC )
      WRITE (*,*) 
      WRITE (*,'(A20,F10.2)') 'Pole RA      ', RA*DPR()
      WRITE (*,'(A20,F10.2)') 'Pole Dec     ', DEC*DPR()

C
C     Compute planetocentric LON/LAT of spacecraft sub-spacecraft point
C     on target using SUBPNT.
C
      CALL SUBPNT ( 'INTERCEPT: ELLIPSOID',
     .              TRGNAM, ET, TRGFRM,  ABCORR, SCNAM, 
     .              SPOINT, TRGEPC, SRFVEC )

      CALL RECLAT ( SPOINT, R, LON, LAT )
      IF ( LON .LT. 0.D0 ) THEN
         LON = LON + TWOPI()
      END IF
      WRITE (*,*) 
      WRITE (*,*)             'Planetocentric:'
      WRITE (*,'(A20,F10.2)') 'Sub-obs lat  ', LAT*DPR()
      WRITE (*,'(A20,F10.2)') 'Sub-obs lon  ', LON*DPR()

      CALL VEQU ( SPOINT, SUBSC )

C
C     Compute platenographic sub-s/c LON/LAT that match Horizons.
C     Because Horizons uses spheroidal model, we need to replace
C     ellipsoidal RADII currently in the POOL with spheroidal RADII,
C     compute intercept, convert it to platenographic LON/LAT, and then
C     restore ellipsoidal RADII in the POOL.
C
      CALL PDPOOL ( 'BODY2000433_RADII', 3, RADSPH )

      CALL SUBPNT ( 'INTERCEPT: ELLIPSOID',
     .              TRGNAM, ET, TRGFRM,  ABCORR, SCNAM, 
     .              SPOINT, TRGEPC, SRFVEC )

      CALL RECPGR ( 'EROS', SPOINT, RE, F, LON, LAT, ALT )

      WRITE (*,*) 
      WRITE (*,*)             'Horizons:'
      WRITE (*,'(A20,F10.2)') 'Sub-obs lat  ', LAT*DPR()
      WRITE (*,'(A20,F10.2)') 'Sub-obs lon  ', LON*DPR()

      CALL PDPOOL ( 'BODY2000433_RADII', 3, RADII )

C
C     Compute planetocentric LON/LAT of sub-solar point on target using
C     SUBSLR.
C
      CALL SUBSLR ( 'INTERCEPT: ELLIPSOID',
     .              TRGNAM, ET, TRGFRM,  ABCORR, SCNAM, 
     .              SPOINT, TRGEPC, SRFVEC )

      CALL RECLAT ( SPOINT, R, LON, LAT )
      IF ( LON .LT. 0.D0 ) THEN
         LON = LON + TWOPI()
      END IF
      WRITE (*,*) 
      WRITE (*,*)             'Planetocentric:'
      WRITE (*,'(A20,F10.2)') 'Sub-sol lat  ', LAT*DPR()
      WRITE (*,'(A20,F10.2)') 'Sub-sol lon  ', LON*DPR()

C
C     Compute platenographic sub-solar LON/LAT that match Horizons.
C     Because Horizons uses spheroidal model, we need to replace
C     ellipsoidal RADII currently in the POOL with spheroidal RADII,
C     compute intercept, convert it to platenographic LON/LAT, and then
C     restore ellipsoidal RADII in the POOL.
C
      CALL PDPOOL ( 'BODY2000433_RADII', 3, RADSPH )

      CALL SUBSLR ( 'INTERCEPT: ELLIPSOID',
     .              TRGNAM, ET, TRGFRM,  ABCORR, SCNAM, 
     .              SPOINT, TRGEPC, SRFVEC )

      CALL RECPGR ( 'EROS', SPOINT, RE, F, LON, LAT, ALT )

      WRITE (*,*) 
      WRITE (*,*)             'Horizons:'
      WRITE (*,'(A20,F10.2)') 'Sub-sol lat  ', LAT*DPR()
      WRITE (*,'(A20,F10.2)') 'Sub-sol lon  ', LON*DPR()

      CALL PDPOOL ( 'BODY2000433_RADII', 3, RADII )

C
C     Compute illumination angles at sub-s/c point. Since this point is
C     on the target-spacecraft line and target is small and far away
C     from the Sun this point's phase angle is nearly the same as the
C     phase angle for the target center.
C
      CALL ILUMIN ( 'Ellipsoid', TRGNAM, ET, TRGFRM,  
     .              ABCORR, SCNAM, SUBSC, TRGEPC,  
     .              SRFVEC, PHASE,  SOLAR,  EMISSN  )

      WRITE (*,*) 
      WRITE (*,'(A20,F10.2)') 'Solar phase  ', PHASE*DPR()

      WRITE (*,*) 
      END