; v 2   2010-03-25   AJS  Modified to use PIXEL_CENTERS from the R-ALice IK

PRO ra_mike_add_geometry, header, test = test

COMMON MIKE_DATA

IF NOT keyword_set(test) THEN mike_log, 'Adding geometry keywords to FITS header.'

radeg = 180D0/!DPI
dtor = !DPI/180D0

badtime_flag = 0
; Get the total number of SPICE kernels that have been loaded and
; their names

CSPICE_KTOTAL, 'SPK CK PCK EK TEXT', count
loaded_kernels = ''
FOR i = 0, count-1 DO BEGIN
   CSPICE_KDATA, i, 'SPK CK PCK EK TEXT' , kernelfile, kerneltype, source, handle, found
   kernelfile = strmid(kernelfile, strpos(kernelfile, '/', /reverse_search) + 1)
   loaded_kernels = [loaded_kernels, kernelfile]
ENDFOR
loaded_kernels = loaded_kernels[1:*]

; Get Ephemeris Time (seconds since J2000 epoch) of midpoint of observations
starttime = sxpar(Header, 'STRTSCET')
stoptime = sxpar(Header, 'STOPSCET')
CSPICE_UTC2ET, starttime, startet
CSPICE_UTC2ET, stoptime, stopet
midpoint_et = (startet + stopet)/2d0
   
iref = 'J2000'
scname = 'ROSETTA'
scid = -226
abcorr = 'LT+S'
inst_frame = 'ROS_ALICE'

; Get the state vector with respect to the Sun
CSPICE_SPKEZR, 'SUN', startet, iref, abcorr, scname, state_sun, lt_sun
sunpos = state_sun[0:2]
sunvel = state_sun[3:5]

CSPICE_SUBPT, 'INTERCEPT', 'SUN', startet, abcorr, scname, subsolarpoint, alt
CSPICE_BODVAR, 10, 'RADII', sunradii
CSPICE_RECPGR, 'SUN', subsolarpoint, sunradii[0], $
               (sunradii[0]-sunradii[2])/sunradii[0], sunlon, sunlat, foo

sunlon *= radeg
IF sunlon LT 0 THEN sunlon += 360.
sunlat *= radeg

; Print out solar geometry keywords
sxaddhist, 'BEGIN GEOMETRY KEYWORD BLOCK', Header, /comment
sxaddhist, 'All values calulated for start time of observation', Header, /comment
sxaddhist, '', Header, /blank
sxaddpar, Header, 'SCSUNX', sunpos[0], ' J2000 EME S/C -> Sun X vector [km]'
sxaddpar, Header, 'SCSUNY', sunpos[1], ' J2000 EME S/C -> Sun Y vector [km]'
sxaddpar, Header, 'SCSUNZ', sunpos[2], ' J2000 EME S/C -> Sun Z vector [km]'
sxaddpar, Header, 'SCSUNVX', sunvel[0], ' X vel. of Sun wrt S/C in [km/s] (J2000 EME)'
sxaddpar, Header, 'SCSUNVY', sunvel[1], ' Y vel. of Sun wrt S/C in [km/s] (J2000 EME)'
sxaddpar, Header, 'SCSUNVZ', sunvel[2], ' Z vel. of Sun wrt S/C in [km/s] (J2000 EME)'
sxaddpar, Header, 'SCSUNLON', sunlon, ' Sub solar longitude [deg.] (planetographic)'
sxaddpar, Header, 'SCSUNLAT', sunlat, ' Sub solar latitude [deg.] (planetographic)'

; Get the state vector with respect to the target
ralice_target, '3', startet, phaseabbrev, phasename, obstype, target_id, target_name, $
               target_type, dataset_target_name, spice_target_name, target_desc 

IF strupcase(sxpar(header, 'APDOOR')) EQ 'CLOSED' THEN BEGIN
   target_desc = '"DARK"'
   target_type = '"CALIBRATION"'
   target_name = '"CALIBRATION"'
   spice_target_name = '"N/A"'
ENDIF

CSPICE_BODN2C, spice_target_name, spice_id, found_target

IF found_target THEN BEGIN

   CSPICE_SPKEZR, spice_target_name, startet, iref, abcorr, scname, state_target, lt_target
   targetpos = state_target[0:2]
   targetvel = state_target[3:5]
   cspice_recrad, targetpos, dist, targ_ra, targ_dec
   targ_ra *= radeg 
   targ_dec *= radeg

; position of the sun as seen from the target at the time a photon
; left the target towards the s/c
   CSPICE_SPKPOS, 'SUN', startet-lt_target, iref, abcorr, spice_target_name, sunpos_targ, $
                  lt_target_sun
   CSPICE_SPKPOS, 'ROSETTA', startet-lt_target, iref, 'XLT+S', spice_target_name, $
                  scpos_targ, LT_target_sc
   phase = CSPICE_VSEP(sunpos_targ, scpos_targ) * radeg

   CSPICE_EXPOOL, 'BODY' + strtrim(spice_id, 2) + '_RADII', found_radii
   IF found_radii THEN BEGIN
      CSPICE_BODVAR, spice_id, 'RADII', radii

      CSPICE_SUBPT, 'INTERCEPT', spice_target_name, startet, abcorr, scname, subpoint, scalt
      CSPICE_RECPGR, spice_target_name, subpoint, radii[0], (radii[0]-radii[2])/radii[0], $
                     sclon, sclat, foo
      sclon *= radeg
      IF sclon LT 0 THEN sclon += 360.
      sclat *= radeg
   ENDIF ELSE BEGIN
      scalt = '"N/A"'
      sclon = '"N/A"'
      sclat = '"N/A"'
   ENDELSE
ENDIF ELSE BEGIN
   targetpos = ['"N/A"', '"N/A"', '"N/A"']
   targetvel = ['"N/A"', '"N/A"', '"N/A"']
   scalt = '"N/A"'
   sclon = '"N/A"'
   sclat = '"N/A"'
   phase = '"N/A"'
   found_radii = 0
   targ_ra = '"N/A"'
   targ_dec = '"N/A"'
ENDELSE

sxaddhist, '', Header, /blank
sxaddhist, 'S/C -> Target Geometry', Header, /comment
sxaddpar, Header, 'OBS_TYPE', obstype, ' Epoch of observation, e.g., PC6 or Steins flyby' 
sxaddpar, Header, 'TARGNAME', target_name, ' Name of observation target'
sxaddpar, Header, 'TARGTYPE', target_type, ' Type of observation target'
sxaddpar, Header, 'TARGDESC', target_desc, ' Description of observation target'
sxaddpar, header, 'TARG_RA', targ_ra, ' Right ascension of target [deg]'
sxaddpar, header, 'TARG_DEC', targ_dec, ' Declination of target [deg]'
sxaddpar, Header, 'SCTARGX', targetpos[0], ' J2000 EME S/C -> target X vector [km]'
sxaddpar, Header, 'SCTARGY', targetpos[1], ' J2000 EME S/C -> target Y vector [km]'
sxaddpar, Header, 'SCTARGZ', targetpos[2], ' J2000 EME S/C -> target Z vector [km]'

IF size(targetpos[0], /type) EQ 7 THEN sxaddpar, Header, 'TARGDIST', '"N/A"', $
  ' Distance to target center [km]' ELSE sxaddpar, Header, 'TARGDIST', $
    sqrt(total(targetpos^2)), ' Distance to target center [km]'
sxaddpar, Header, 'SCTARGVX', targetvel[0], ' X vel. of target wrt S/C in [km/s] (J2000 EME)'
sxaddpar, Header, 'SCTARGVY', targetvel[1], ' Y vel. of target wrt S/C in [km/s] (J2000 EME)'
sxaddpar, Header, 'SCTARGVZ', targetvel[2], ' Z vel. of target wrt S/C in [km/s] (J2000 EME)'
sxaddpar, Header, 'SCLON', sclon, ' Sub-Spacecraft planetographic longitude [deg.]'
sxaddpar, Header, 'SCLAT', sclat, ' Sub-Spacecraft planetographic latitude [deg.]'
sxaddpar, Header, 'SCALT', scalt, ' Spacecraft altitude [km]'
sxaddpar, Header, 'TARGPHAS', phase, ' Sub-S/C Phase angle (Sun-Target-S/C) [deg]'

; Get the transformation vector from the instrument frame to the
; inertial reference frame at the mid-point of the observation. 
CSPICE_PXFORM, inst_frame, iref, startet, inst2inertial
boresight = [0d, 0d, 1d]

sxaddhist, '', Header, /blank
sxaddhist, 'Row-specific Geometry Information', Header, /comment
sxaddhist, 'The ALICE slit is imaged onto detector rows 5-23 (inclusive)', Header, /comment
sxaddhist, 'When the LOS intercepts the surface, LAT and LON are given for', Header, /comment
sxaddhist, 'the intercept point. If the LOS fails to intercept the target,', Header, /comment
sxaddhist, 'LAT and LON are given for the nearest point to the target surface', Header, /comment
sxaddhist, 'along the LOS.', Header, /comment

; Define the Target reference frame and check for its existence
CSPICE_CIDFRM, spice_id, frame_code, target_frame, found_frame

IF found_target AND found_frame THEN $ ;AND (target_frame NE 'STEINS_FIXED') THEN $
  CSPICE_SPKPOS, scname, startet-lt_target, target_frame, 'XLT+S', $
                 spice_target_name, pos_targ, lt_target2

; Get the row vectors in the instrument frame
CSPICE_GDPOOL, 'INS-226120_PIXEL_CENTERS', 0, 500, row_vectors, found
row_vectors = reform(row_vectors, 3, 32)

FOR i = 5, 23 DO BEGIN 
   LOS = row_vectors[*, i]
   LOS_inertial = reform(inst2inertial ## row_vectors[*, i])
   CSPICE_RECRAD, LOS_inertial, radius, ra, dec
   ra *= radeg
   dec *= radeg

   IF found_target AND found_radii AND found_frame THEN BEGIN
      CSPICE_PXFORM, inst_frame, target_frame, startet-lt_target, inst2targ
      LOS_targ = reform(inst2targ ## LOS)

      CSPICE_SURFPT, pos_targ, LOS_targ, radii[0], radii[1], radii[2], xpoint, intercept_found

      IF intercept_found THEN BEGIN
         CSPICE_ILLUM, spice_target_name, startet-lt_target, abcorr, scname, xpoint, $
                       phase, incidence, emission
         phase *= RADEG
         incidence *= RADEG
         emission *= RADEG

         CSPICE_RECPGR, spice_target_name, xpoint, radii[0], (radii[0]-radii[2])/radii[0], $
                        lon, lat, foo
         lon *= RADEG
         IF lon LT 0. THEN lon += 360.
         lat *= RADEG
      ENDIF ELSE BEGIN
; Don't calculate Lat and Lon unless the lOS intersects with the
; target body.
;         CSPICE_NPEDLN, radii[0], radii[1], radii[2], pos_targ, los_targ, nearpoint, rayheight
;         CSPICE_RECPGR, spice_target_name, nearpoint, radii[0], (radii[0]-radii[2])/radii[0], $
;                        lon, lat, foo
;         lon *= RADEG
;         IF lon LT 0. THEN lon += 360.
;         lat *= RADEG

         lat = '"N/A"'
         lon = '"N/A"'
         phase = '"N/A"'
         incidence = '"N/A"'
         emission = '"N/A"'
      ENDELSE
   ENDIF ELSE BEGIN
      lat = '"N/A"'
      lon = '"N/A"'
      phase = '"N/A"'
      incidence = '"N/A"'
      emission = '"N/A"'
   ENDELSE

; Write the values to the FITS header
   sxaddpar, Header, 'ROW' + strtrim(i, 2) + 'RA', ra, ' Right Ascension [deg.]'
   sxaddpar, Header, 'ROW' + strtrim(i, 2) + 'DEC', dec, ' Declination [deg.]'
   sxaddpar, Header, 'ROW' + strtrim(i, 2) + 'LON', lon, $
             ' Planetographic Lon. of target along LOS [deg.]'
   sxaddpar, Header, 'ROW' + strtrim(i, 2) + 'LAT', lat, $
             ' Planetographic Lat. of target along LOS [deg.]'
   sxaddpar, Header, 'ROW' + strtrim(i, 2) + 'PHA', phase, $
             ' Phase angle at LOS intercept [deg.]'
   sxaddpar, Header, 'ROW' + strtrim(i, 2) + 'INC', incidence, $
             ' Incidence angle at LOS intercept [deg.]'
   sxaddpar, Header, 'ROW' + strtrim(i, 2) + 'EMI', emission, $
             ' Emission angle at LOS intercept [deg.]'

; Add boresight RA and DEC keywords
   IF i EQ 15 THEN BEGIN
      sxaddpar, Header, 'RA', ra, ' Right Ascension of ALICE Boresight [deg.]', $
                before = 'SCSUNX'
      sxaddpar, Header, 'DEC', dec, ' Declination of ALICE Boresight [deg.]', after = 'RA'
      sxaddhist, ' ', Header, /blank, location = 'SCSUNX'
      sxaddhist, 'S/C -> Sun Geometry', Header, /comment, location = 'SCSUNX'
   ENDIF

ENDFOR

sxaddhist, '', Header, /blank
sxaddhist, 'SPICE files used in geometry calculations', Header, /comment
sxaddpar, Header, 'NUMBSPCK', n_elements(loaded_kernels), ' Number of loaded SPICE kernels'
FOR i = 0, n_elements(loaded_kernels)-1 DO BEGIN
   kernelnum = strtrim(i + 1, 2)
   WHILE strlen(kernelnum) LT 4 DO kernelnum = '0' + kernelnum
   sxaddpar, Header, 'SPCK' + kernelnum, loaded_kernels[i], ' SPICE kernel ' + kernelnum
ENDFOR


END