;+
; Project     :	STEREO
;
; Name        :	WCS_SIP_TRANSLATE()
;
; Purpose     :	Add SIP convention to WCS structure
;
; Category    :	FITS, Coordinates, WCS
;
; Explanation :	This routine looks for the keywords relevant to the Simple
;               Image Polynomial (SIP) convention for image distortions, and
;               translates them to the distortion system used by the Calabretta
;               et al. distortion paper (WCS Paper IV), adding the resulting
;               DISTORTION structure to the WCS structure.
;
; Syntax      :	NewWCS = WCS_SIP_TRANSLATE(OrigWCS, HEADER)
;
; Examples    :	See FITSHEAD2WCS.
;
; Inputs      :	OrigWCS = The original WCS structure, without the DISTORTION
;                         structure.
;
;               HEADER  = The original FITS header from which OrigWCS was
;                         derived.
;
; Opt. Inputs :	None
;
; Outputs     :	The output of the function is the WCS structure with the
;               DISTORTION structure appended to it.
;
; Opt. Outputs:	None
;
; Keywords    :	SYSTEM  = Selects which alternate coordinate system should be
;                         used.  See wcs_find_system.pro for more information.
;                         Must be the same system that OrigWCS is based on.
;
; Calls       :	DATATYPE, FITSHEAD2STRUCT, IS_STRING, WCS_FIND_SYSTEM,
;               TAG_EXIST, ADD_TAG
;
; Common      :	None
;
; Restrictions:	None
;
; Side effects:	None
;
; Prev. Hist. :	None
;
; History     :	Version 1, 28-Apr-2025, William Thompson, GSFC
;
; Contact     :	WTHOMPSON
;-
;
function wcs_sip_translate, wcs, header, system=k_system
;
;  Distinguish between string headers and structures.
;
if datatype(header) eq 'STR' then hdr = fitshead2struct(header) else $
  hdr = header
;
;  Determine which coordinate system to use.
;
if not is_string(k_system) then system = '' else $
  system = wcs_find_system(hdr, k_system)
;
;  Step through the axes.
;
prefixes = ['A','B'] + '_'
for iaxis=0,1 do begin
    prefix = prefixes[iaxis]
;
;  Test if the ORDER keyword exists for this axis.
;
    test = tag_exist(hdr, prefix+'ORDER' + system, index=index)
    if test then begin
        n_order = hdr.(index)
        param = ['NAXES', 'AXIS.1', 'OFFSET.1', 'AXIS.2', 'OFFSET.2', 'NTERMS']
        value = [2, 1, hdr.crpix1, 2, hdr.crpix2, 0]
;
;  Step through the possible orders, and find the coefficients.
;
        nterms = 0
        for i=0,n_order do begin
            for j=0,n_order do begin
                keyword = prefix + ntrim(i) + '_' + ntrim(j) + system
                test = tag_exist(hdr, keyword, index=index)
                if test then begin
                    nterms = nterms + 1
                    coeff = hdr.(index)
                    term = 'TERM.' + ntrim(nterms) + '.'
                    param = [param, term + 'COEFF']
                    value = [value, coeff]
                    param = [param, term + 'VAR.1']
                    value = [value, i]
                    param = [param, term + 'VAR.2']
                    value = [value, j]
                endif
            endfor
        endfor
        value[5] = nterms
        dp = {param: param, value: value, cdis: 'Polynomial'}
        case iaxis of
            0: distortion = {DP1: dp}
            1: if n_elements(distortion) eq 0 then distortion = {DP2: dp} else $
              distortion = add_tag(distortion, dp, 'DP2')
        endcase
    endif
endfor
;
wcssip = wcs
if n_elements(distortion) ne 0 then $
  wcssip = add_tag(wcssip, distortion, 'DISTORTION')
return, wcssip
end