;+
; $Id: scc_roll_image.pro,v 1.7 2016/02/29 16:53:07 thompson Exp $
; Project     :	STEREO - SECCHI
;
; Name        :	SCC_ROLL_IMAGE
;
; Purpose     :	Roll to align image to a specified coordinate system
;
; Category    :	SECCHI, coordinates
;
; Explanation : Rotate the image by the required amount to align it with one of
;               the specified coordinate systems.  The image is rotated about
;               the reference pixel to retain the ability to use the coordinate
;               information in the header.
;
; Syntax      :	SCC_ROLL_IMAGE, HEADER, IMAGE  [, SYSTEM=SYSTEM ]
;
; Examples    :	SECCHI_PREP, FILENAME, HEADER, IMAGE
;               SCC_ROLL_IMAGE, HEADER, IMAGE           ;Put solar north up
;
;               Use the following to align in R.A. & Dec.
;
;               SCC_ROLL_IMAGE, HEADER, IMAGE, SYSTEM='GEI'
;
;               Or this for ecliptic coordinates.
;
;               SCC_ROLL_IMAGE, HEADER, IMAGE, SYSTEM='HAE'
;
; Inputs      :	HEADER = SECCHI header structure
;               IMAGE  = Image array
;
; Opt. Inputs :	None.
;
; Outputs     :	HEADER and IMAGE are modified.
;
; Opt. Outputs:	None.
;
; Keywords    :	SYSTEM = Character string, giving one of the following
;                        standard coordinate systems:
;
;                               GEI     Geocentric Equatorial Inertial
;                               GEO     Geographic
;                               GSE     Geocentric Solar Ecliptic
;                               MAG     Geomagnetic
;                               GSM     Geocentric Solar Magnetospheric
;                               SM      Solar Magnetic
;                               HCI     Heliocentric Inertial
;                               HAE     Heliocentric Aries Ecliptic
;                               HEE     Heliocentric Earth Ecliptic
;                               HEEQ    Heliocentric Earth Equatorial (or HEQ)
;                               Carrington (can be abbreviated)
;                               HGRTN   Heliocentric Radial-Tangential-Normal
;                               RTN     Radial-Tangential-Normal (default)
;                               HPC     Helioprojective-Cartesian
;                               HERTN   Heliocentric Ecliptic RTN
;                               STPLN   STEREO Mission Plane
;
;                        For the purposes of this routine, HPC is treated as a
;                        synonym of RTN.
;
;               NAN_OFF = Normally, pixels outside the image are marked as
;                         Not-A-Number after rolling.  Use /NAN_OFF to mark
;                         missing pixels with HEADER.BLANK (usually zero).
;
;               INTERP  = Set to use bilinear interpolation.  The default is
;                         nearest neighbor sampling.
;
;               OFFSET  = Two-element pixel offset to apply to the image as
;                         part of the interpolation.  This option supports the
;                         /DEJITTER_ON feature in EUVI_PREP.
;
;               CUBIC   = Set -1 <= CUBIC < 0 to use cubic interpolation.  See
;                         the ROT documentation for more details.
;
;               ERRMSG = If defined and passed, then any error messages will be
;                        returned to the user in this parameter rather than
;                        depending on the MESSAGE routine in IDL.  If no errors
;                        are encountered, then a null string is returned.  In
;                        order to use this feature, ERRMSG must be defined
;                        first, e.g.
;
;                               ERRMSG = ''
;                               SCC_ROLL_IMAGE, ERRMSG=ERRMSG, ...
;                               IF ERRMSG NE '' THEN ...
;
; Calls       :	FITSHEAD2WCS, GET_STEREO_ROLL, ROT
;
; Common      :	None.
;
; Restrictions:	Currently only works for the SCIP telescopes.
;
; Side effects:	FITS string header arrays are converted into SECCHI header
;               structures.
;
; Prev. Hist. :	None.
;
; $Log: scc_roll_image.pro,v $
; Revision 1.7  2016/02/29 16:53:07  thompson
; Added keyword OFFSET to support EUVI_PREP,/DEJITTER_ON
;
; Revision 1.6  2012/08/20 20:20:13  thompson
; For complete consistency, use WCS software to calculate XCEN, YCEN
;
; Revision 1.4  2011/08/03 17:59:11  nathan
; make sure crota LT 180
;
; Revision 1.3  2008/09/24 17:28:46  nathan
; update xcen/ycen; added Id and Log cvs tags
;
; History     :	Version 1, 30-Jul-2007, William Thompson, GSFC
;
; Contact     :	WTHOMPSON
;-
;
pro scc_roll_image, header, image, system=k_system, nan_off=nan_off, $
                    errmsg=errmsg, interp=interp, offset=offset, _extra=_extra
;
;  Check the input parameters.
;
if n_params() ne 2 then begin
    message = 'Syntax: SCC_ROLL_IMAGE, HEADER, IMAGE  [, SYSTEM=SYSTEM]'
    goto, handle_error
endif
;
;  Get the name of the system to align to.
;
if n_elements(k_system) eq 1 then system=strupcase(k_system) else system='RTN'
if system eq 'HPC' then system = 'RTN'
;
;  If necessary, convert a string header into a SECCHI header structure.
;
if (datatype(header) ne 'STC') then header = scc_fitshdr2struct(header)
;
;  Get the roll angle to apply, taking into account the current image
;  orientation.
;
wcs = fitshead2wcs(header)
roll = -wcs.roll_angle
if system ne 'RTN' then begin
    date = header.date_obs
    spacecraft = header.obsrvtry
    message = ''
    roll = roll + get_stereo_roll(date, spacecraft, system=system, $
                                  errmsg=message)
    if message ne '' then goto, handle_error
    roll = roll - get_stereo_roll(date, spacecraft, system='RTN')
endif
;
;  Rotate the image, and apply any offset.  This section includes software
;  taken from the standard IDL routine rot.pro.
;
case datatype(image,2) of
    4: missing = !values.f_nan
    5: missing = !values.d_nan
    else: missing = header.blank
endcase
if keyword_set(nan_off) then missing = header.blank
;
theta = -roll * !dpi / 180.d0
cosr = cos(theta)
sinr = sin(theta)
;
x0 = wcs.crpix[0] - 1
y0 = wcs.crpix[1] - 1
if n_elements(offset) eq 2 then begin
    xc = x0 + offset[0]
    yc = y0 + offset[1]
end else begin
    xc = x0
    yc = y0
endelse
;
kx = -xc + cosr*x0 - sinr*y0
ky = -yc + sinr*x0 + cosr*y0
kk = 1.d0 / (1.d0 + sinr^2/cosr^2)
;
cx = kk * [ ky*sinr/cosr^2 + kx/cosr,  sinr/cosr^2, 1/cosr, 0.d0]
cy = kk * [-kx*sinr/cosr^2 + ky/cosr, 1/cosr, -sinr/cosr^2, 0.d0]
;
image = poly_2d(image, cx, cy, keyword_set(interp), missing=missing, $
                _extra=_extra)
;
;  Update the primary (HPC) header info with the new roll angle.
;
header.crota = header.crota + roll
IF header.crota GT 180. THEN header.crota = header.crota-360.
lambda = header.cdelt2 / header.cdelt1
cos_a = cos(header.crota * !dpi / 180.d0)
sin_a = sin(header.crota * !dpi / 180.d0)
header.pc1_1 = cos_a
header.pc1_2 = -lambda * sin_a
header.pc2_1 = sin_a / lambda
header.pc2_2 = cos_a
;
;  And the offset, if any.
;
if n_elements(offset) eq 2 then begin
    header.crpix1  = header.crpix1  + offset[0]
    header.crpix2  = header.crpix2  + offset[1]
    header.crpix1a = header.crpix1a + offset[0]
    header.crpix2a = header.crpix2a + offset[1]
endif
;
;  Update the XCEN, YCEN keywords, if found in the header.
;
if tag_exist(header,'xcen') then begin
    wcs1 = fitshead2wcs(header)
    center = (wcs1.naxis - 1) / 2.
    center = wcs_get_coord(wcs1, center)
    header.xcen = center[0]
    header.ycen = center[1]
endif
;
;  Update the secondary (RA/Dec) header info.
;
wcsa = fitshead2wcs(header, system='a')
crota = wcsa.roll_angle - roll
lambda = header.cdelt2a / header.cdelt1a
cos_a = cos(crota * !dpi / 180.d0)
sin_a = sin(crota * !dpi / 180.d0)
header.pc1_1a = cos_a
header.pc1_2a = -lambda * sin_a
header.pc2_1a = sin_a / lambda
header.pc2_2a = cos_a
;
return
;
;  Error handling point.
;
handle_error:
if n_elements(errmsg) eq 0 then message, message else $
  errmsg = 'SCC_ROLL_IMAGE: ' + message
;
end