;+
; Name    : iris_isp2wcs
;
; Purpose : return header structure with IRIS WCS keywords added
;
; Use     : IDL> iris_isp2wcs,index
;
; Input/Output   :
;    index = FITS header index structure.  may be a vector
;
; Keywords input:
;    roll_angle = float.  Default: use hdstr.sat_rot
;    file = string. name of geny database file to be used.
;                   Default (for now): create database with iris_mk_pointdb.pro
;    version = long.  Use latest database version with db.version le version.
;                     Default: 99999L
;    /level1point5  If set modify WCS parameters appropriate for Level 1.5 data
;    temperatures = struct. Structure with spectrograph temperature tags.
;                           Default: use temps in index, if there, otherwise
;                                    use database of daily averages.
;
; Keywords output:
;
; Common  :
;
; Restrictions:
;    Assumes that all the IRIS ISP keywords exist in the input header
;    This routine doesn't fill in the WCS position keywords such as
;    hgln_obs,hglt_obs,dsun_obs.
;
; Side effects:
;
; Category   : Pipeline
;
; Prev. Hist.:
;
; Written: Jean-Pierre Wuelser, LMSAL, 31-May-2013
;          JPW 2013-08-19 support for summed data
;          JPW 2013-09-11 accommodate changed tags in pointdb
;          JPW 2013-09-21 fix bug with scroll in vector call
;          JPW 2013-10-04 add /level1point5
;          JPW 2013-10-14 added code for xcen,ycen for those who no speak WCS
;          JPW 2014-02-04 added wavelength corr. using temp/ophase/pointing 
;          JPW 2014-02-20 bug fix per Paul B. & (re-)adding roll in orb.phase
;          JPW 2014-02-26 handling of NaN in ophase and sat_rot
;          JPW 2014-02-28 temperature corr. of slit center in axis2
;          JPW 2014-03-03 fix for bimodal FW position during early mission
;                         new temp.format incl wave corr., temp.corr SJIx&y
;          JPW 2014-03-03 added separate fid_fu1 and fid_fu2
;          JPW 2017-04-20 added secular drift term to wav, fid and fi1
;-

pro iris_isp2wcs,hd,roll_angle=scroll,file=file,version=ver,level1point5=l1p5, $
                 temperatures=temps

common iris_temp_pointdb, tempdb    ; tempory implementation of wm_yoff fix

if n_elements(ver) ne 1 then ver=99999L

; get parameter database
if keyword_set(file) then begin
   restgenx,file=file,db
   w=where(db.version le ver,nw)
   if nw gt 0 then db=db[w]
endif else db=iris_mk_pointdb(ver=ver)

; set some constants
yr2s = 3.6525d2*2.4d1*3.6d3        ; # of seconds in a year
tref = anytim('2013-09-17T00:00')  ; ref.time = time of flight software fix

; Will need SG temperatures for wavelength correction.  3 sources:
; 1. optional temps keyword structure
; 2. image header
; 3. daily temperature database
if ~keyword_set(l1p5) then begin
 ; tag names of temperatures
 ttag = 'IT1'+['4','6','5','7','8','9']+'SP'+['PX','NX','PY','NY','PZ','NZ']
 ntag = n_elements(ttag)
 wtag = lonarr(ntag)-1L
 wtag2 = lonarr(ntag)-1L
 tempsrc = 0
 ; 1. temps keyword?
 sizt = size(temps)
 if sizt[sizt[0]+1] eq 8 and sizt[sizt[0]+2] eq n_elements(hd) then begin
   namtag = tag_names(temps)
   for j=0,ntag-1 do wtag[j]=where(namtag eq ttag[j])
   if min(wtag) ge 0 then tempsrc=1
 endif
 ; 2. image header?
 if tempsrc eq 0 then begin
   namtag = tag_names(hd)
   for j=0,ntag-1 do wtag[j]=where(namtag eq ttag[j])
   if min(wtag) ge 0 then tempsrc=2
 endif
 ; 3. daily temperature database?  might need it in any case of NaNs
 if n_elements(tempdb) eq 0 then $
   restgenx,file=concat_dir('$SSW_IRIS_DATA','iris_temp_pointdb.geny'),tempdb
 if n_elements(tempdb) gt 0 then begin
   namtag = tag_names(tempdb)
   for j=0,ntag-1 do wtag2[j]=where(namtag eq ttag[j])
   if min(wtag2) ge 0 and tempsrc eq 0 then tempsrc=3
 endif
endif else tempsrc=0

; Can we re-compute xcen,ycen?  Check whether xcen and naxis1 exist
if tag_exist(hd,'xcen') and tag_exist(hd,'naxis1') then doxycen=1 $
                                                   else doxycen=0

; SC roll explicit?
if n_elements(scroll) eq 1 then hd.sat_rot=scroll[0]
if n_elements(scroll) eq n_elements(hd) then hd.sat_rot=scroll

; loop through header index (allows hd to be vectorized)
for j=0,n_elements(hd)-1 do begin

 if tag_exist(hd[j],'date_obs') then t=anytim(hd[j].date_obs) else begin
   if tag_exist(hd[j],'date_d$obs') then t=anytim(hd[j].date_d$obs) else t=tref
 endelse

 ; in case sat_rot is NaN
 if finite(hd[j].sat_rot,/nan) then hd[j].sat_rot=0.

 ; find appropriate entry in parameter database
 ndb = n_elements(db)
 if ndb gt 0 and t ne tref then begin      ; t eq tref: not a valid image
   w = where(t ge anytim(db.t_start) and t lt anytim(db.t_stop) and $
             ver ge db.version,nw)
   if nw gt 0 then begin
      ww=where(db[w].version eq max(db[w].version),nww)
      d = db[w[ww[nww-1]]]
   endif else d=db[ndb-1]
 endif else d=db[ndb-1]
 
 ; Get vector of SG temperatures for wavelength correction
 if tempsrc gt 0 then begin
   tempi = fltarr(ntag)
   if tempsrc eq 1 then for k=0,ntag-1 do tempi[k]=temps[j].(wtag[k])
   if tempsrc eq 2 then for k=0,ntag-1 do tempi[k]=hd[j].(wtag[k])
   if (max(finite(tempi,/nan)) gt 0 or min(tempi) lt -273.0 or tempsrc eq 3) $
     then if min(wtag2) ge 0 then for k=0,ntag-1 do $
          tempi[k]=interpol(tempdb.(wtag2[k]),tempdb.anytim,t)
 endif

 ; Use Ted's routine to get CRVAL
 crvalxy = iris_isp2solar(hd[j],roll_angle=hd[j].sat_rot,db=d)
 garad = (-d.slit_rot-hd[j].sat_rot) / !radeg
 
 ; summing factors and offsets used in crpixj and pci_j calculations
 sum1 = float(hd[j].sumsptrl)
 sum2 = float(hd[j].sumspat)
 off1 = (sum1-1.0)/(sum1*2.0)
 off2 = (sum2-1.0)/(sum2*2.0)

 ; set keywords: separate cases for FUV, NUV, SJI
 case hd[j].instrume of
  'FUV': begin
    hd[j].wcsaxes = 3
    if keyword_set(l1p5) then begin
       berad = 0.0
       alrad = 0.0
       cpx2  = d.cpx2_1p5
       cdlt2 = d.cdlt_1p5
    endif else begin
       berad = d.be_fu1 / !radeg
       alrad = d.al_fu1 / !radeg
       cpx2  = d.cpx2_fu1
       cdlt2 = d.cdlt2_f1
    endelse
    if tempsrc eq 0 or finite(hd[j].ophase,/nan) then dcpx1=0.0 $ ; incl. l1p5
    else begin
       ophase = 2.0*!pi*hd[j].ophase + hd[j].sat_rot/!radeg
       dcpx1 = hd[j].obs_vr*d.cvl1_fu1/(2.998e8*d.cdlt1_f1) $     ; rad. veloc.
             + d.wvl_foc1*sin(ophase) + d.wvl_foc2*cos(ophase) $  ; orb. corr.
             + d.wvl_fpc1*crvalxy[0]*cos(hd[j].sat_rot/!radeg) $  ; pnt. corr.
             + d.wvl_fsc1*exp(d.wvl_fsc2*(t-tref)/yr2s) $         ; sec. drift
             + d.wvl_ftc0 + total(d.wvl_ftc1*tempi)               ; temp.corr
    endelse
    if tempsrc eq 0 then dcpx2=0.0 $
       else dcpx2 = d.fid_fu1t0 + total(d.fid_fu1t1*tempi) $
                  + d.fid_fu1s1*exp(d.fid_fu1s2*(t-tref)/yr2s)
    hd[j].crpix1  = (d.cpx1_fu1 + dcpx1) / sum1 + off1
    hd[j].crpix2  = (cpx2 + dcpx2) / sum2 + off2
    hd[j].crval1  = d.cvl1_fu1
    hd[j].crval2  = crvalxy[1]
    hd[j].crval3  = crvalxy[0]
    hd[j].cdelt1  = d.cdlt1_f1 * sum1
    hd[j].cdelt2  = cdlt2 * sum2
    hd[j].cdelt3  = hd[j].cdelt2
    hd[j].ctype1  = 'WAVE    '
    hd[j].ctype2  = 'HPLT-TAN'
    hd[j].ctype3  = 'HPLN-TAN'
    hd[j].cunit1  = 'Angstrom'
    hd[j].cunit2  = 'arcsec  '
    hd[j].cunit3  = hd[j].cunit2
    hd[j].pc1_1   = ( cos(berad) / cos(alrad-berad))
    hd[j].pc1_2   = (-sin(berad) / cos(alrad-berad)) * sum2/sum1
    hd[j].pc2_1   = ( cos(garad) * sin(alrad) / cos(alrad-berad)) * sum1/sum2
    hd[j].pc2_2   = ( cos(garad) * cos(alrad) / cos(alrad-berad))
    hd[j].pc3_1   = (-sin(garad) * sin(alrad) / cos(alrad-berad)) * sum1/sum2
    hd[j].pc3_2   = (-sin(garad) * cos(alrad) / cos(alrad-berad))
    if keyword_set(l1p5) then begin
       berad = 0.0
       alrad = 0.0
       cpx2  = d.cpx2_1p5
       cdlt2 = d.cdlt_1p5
    endif else begin
       berad = d.be_fu2 / !radeg
       alrad = d.al_fu2 / !radeg
       cpx2  = d.cpx2_fu2
       cdlt2 = d.cdlt2_f2
    endelse
    if tempsrc eq 0 or finite(hd[j].ophase,/nan) then dcpx1=0.0 $ ; incl. l1p5
    else begin
       ophase = 2.0*!pi*hd[j].ophase + hd[j].sat_rot/!radeg
       dcpx1 = hd[j].obs_vr*d.cvl1_fu2/(2.998e8*d.cdlt1_f2) $     ; rad. veloc.
             + d.wvl_foc1*sin(ophase) + d.wvl_foc2*cos(ophase) $  ; orb. corr.
             + d.wvl_fpc1*crvalxy[0]*cos(hd[j].sat_rot/!radeg) $  ; pnt. corr.
             + d.wvl_fsc1*exp(d.wvl_fsc2*(t-tref)/yr2s) $         ; sec. drift
             + d.wvl_ftc0 + total(d.wvl_ftc1*tempi)               ; temp.corr
    endelse
    if tempsrc eq 0 then dcpx2=0.0 $
       else dcpx2 = d.fid_fu2t0 + total(d.fid_fu2t1*tempi) $
                  + d.fid_fu2s1*exp(d.fid_fu2s2*(t-tref)/yr2s)
    hd[j].crpix1a = (d.cpx1_fu2 + dcpx1) / sum1 + off1
    hd[j].crpix2a = (cpx2 + dcpx2) / sum2 + off2
    hd[j].crval1a = d.cvl1_fu2
    hd[j].crval2a = hd[j].crval2
    hd[j].crval3a = hd[j].crval3
    hd[j].cdelt1a = d.cdlt1_f2 * sum1
    hd[j].cdelt2a = cdlt2 * sum2
    hd[j].cdelt3a = hd[j].cdelt2a
    hd[j].ctype1a = hd[j].ctype1
    hd[j].ctype2a = hd[j].ctype2
    hd[j].ctype3a = hd[j].ctype3
    hd[j].cunit1a = hd[j].cunit1
    hd[j].cunit2a = hd[j].cunit2
    hd[j].cunit3a = hd[j].cunit3
    hd[j].pc1_1a  = ( cos(berad) / cos(alrad-berad))
    hd[j].pc1_2a  = (-sin(berad) / cos(alrad-berad)) * sum2/sum1
    hd[j].pc2_1a  = ( cos(garad) * sin(alrad) / cos(alrad-berad)) * sum1/sum2
    hd[j].pc2_2a  = ( cos(garad) * cos(alrad) / cos(alrad-berad))
    hd[j].pc3_1a  = (-sin(garad) * sin(alrad) / cos(alrad-berad)) * sum1/sum2
    hd[j].pc3_2a  = (-sin(garad) * cos(alrad) / cos(alrad-berad))
    hd[j].lonpole = 180.0
    hd[j].specsys = 'HELIOCEN'    ; or 'SOURCE  ' if corrected for solar rot
    ; update xcen,ycen (this is not an endorsement of those keywords)
    if doxycen then begin
     hd[j].xcen = hd[j].crval3 + hd[j].cdelt3 * $
                  (hd[j].pc3_1*((hd[j].naxis1+1)/2.-hd[j].crpix1) $
                  +hd[j].pc3_2*((hd[j].naxis2+1)/2.-hd[j].crpix2))
     hd[j].ycen = hd[j].crval2 + hd[j].cdelt2 *$
                  (hd[j].pc2_1*((hd[j].naxis1+1)/2.-hd[j].crpix1) $
                  +hd[j].pc2_2*((hd[j].naxis2+1)/2.-hd[j].crpix2))
    endif
  end
  'NUV': begin
    hd[j].wcsaxes = 3
    berad = d.be_nuv / !radeg
    alrad = d.al_nuv / !radeg
    if keyword_set(l1p5) then begin
       berad = 0.0
       alrad = 0.0
       cpx2  = d.cpx2_1p5
       cdlt2 = d.cdlt_1p5
    endif else begin
       berad = d.be_nuv / !radeg
       alrad = d.al_nuv / !radeg
       cpx2  = d.cpx2_nuv
       cdlt2 = d.cdlt2_nu
    endelse
    if tempsrc eq 0 or finite(hd[j].ophase,/nan) then dcpx1=0.0 $ ; incl. l1p5
    else begin
       ophase = 2.0*!pi*hd[j].ophase + hd[j].sat_rot/!radeg
       dcpx1 = hd[j].obs_vr*d.cvl1_nuv/(2.998e8*d.cdlt1_nu) $     ; rad. veloc.
             + d.wvl_noc1*sin(ophase) + d.wvl_noc2*cos(ophase) $  ; orb. corr.
             + d.wvl_npc1*crvalxy[0]*cos(hd[j].sat_rot/!radeg) $  ; pnt. corr.
             + d.wvl_nsc1*exp(d.wvl_nsc2*(t-tref)/yr2s) $         ; sec. drift
             + d.wvl_ntc0 + total(d.wvl_ntc1*tempi)               ; temp.corr
    endelse
    if tempsrc eq 0 then dcpx2=0.0 $
       else dcpx2 = d.fid_nuvt0 + total(d.fid_nuvt1*tempi) $
                  + d.fid_nuvs1*exp(d.fid_nuvs2*(t-tref)/yr2s)
    hd[j].crpix1  = (d.cpx1_nuv + dcpx1) / sum1 + off1
    hd[j].crpix2  = (cpx2 + dcpx2) / sum2 + off2
    hd[j].crval1  = d.cvl1_nuv
    hd[j].crval2  = crvalxy[1]
    hd[j].crval3  = crvalxy[0]
    hd[j].cdelt1  = d.cdlt1_nu * sum1
    hd[j].cdelt2  = cdlt2 * sum2
    hd[j].cdelt3  = hd[j].cdelt2
    hd[j].ctype1  = 'WAVE    '
    hd[j].ctype2  = 'HPLT-TAN'
    hd[j].ctype3  = 'HPLN-TAN'
    hd[j].cunit1  = 'Angstrom'
    hd[j].cunit2  = 'arcsec  '
    hd[j].cunit3  = hd[j].cunit2
    hd[j].pc1_1   = ( cos(berad) / cos(alrad-berad))
    hd[j].pc1_2   = (-sin(berad) / cos(alrad-berad)) * sum2/sum1
    hd[j].pc2_1   = ( cos(garad) * sin(alrad) / cos(alrad-berad)) * sum1/sum2
    hd[j].pc2_2   = ( cos(garad) * cos(alrad) / cos(alrad-berad))
    hd[j].pc3_1   = (-sin(garad) * sin(alrad) / cos(alrad-berad)) * sum1/sum2
    hd[j].pc3_2   = (-sin(garad) * cos(alrad) / cos(alrad-berad))
    hd[j].crpix1a = hd[j].crpix1
    hd[j].crpix2a = hd[j].crpix2
    hd[j].crval1a = hd[j].crval1
    hd[j].crval2a = hd[j].crval2
    hd[j].crval3a = hd[j].crval3
    hd[j].cdelt1a = hd[j].cdelt1
    hd[j].cdelt2a = hd[j].cdelt2
    hd[j].cdelt3a = hd[j].cdelt3
    hd[j].ctype1a = hd[j].ctype1
    hd[j].ctype2a = hd[j].ctype2
    hd[j].ctype3a = hd[j].ctype3
    hd[j].cunit1a = hd[j].cunit1
    hd[j].cunit2a = hd[j].cunit2
    hd[j].cunit3a = hd[j].cunit3
    hd[j].pc1_1a  = hd[j].pc1_1
    hd[j].pc1_2a  = hd[j].pc1_2
    hd[j].pc2_1a  = hd[j].pc2_1
    hd[j].pc2_2a  = hd[j].pc2_2
    hd[j].pc3_1a  = hd[j].pc3_1
    hd[j].pc3_2a  = hd[j].pc3_2
    hd[j].lonpole = 180.0
    hd[j].specsys = 'HELIOCEN'    ; or 'SOURCE  ' if corrected for solar rot
    ; update xcen,ycen (this is not an endorsement of those keywords)
    if doxycen then begin
     hd[j].xcen = hd[j].crval3 + hd[j].cdelt3 * $
                  (hd[j].pc3_1*((hd[j].naxis1+1)/2.-hd[j].crpix1) $
                  +hd[j].pc3_2*((hd[j].naxis2+1)/2.-hd[j].crpix2))
     hd[j].ycen = hd[j].crval2 + hd[j].cdelt2 *$
                  (hd[j].pc2_1*((hd[j].naxis1+1)/2.-hd[j].crpix1) $
                  +hd[j].pc2_2*((hd[j].naxis2+1)/2.-hd[j].crpix2))
    endif
  end
  'SJI': begin
    hd[j].wcsaxes = 2
    if keyword_set(l1p5) then begin
       berad = 0.0
       cdelt = d.cdlt_1p5
       hd[j].crpix1 = d.cpx1_1p5
       hd[j].crpix2 = d.cpx2_1p5       
    endif else begin
     case ((hd[j].ifwpos+10) mod 180)/30 of   ; or case hd[j].img_path of
      0 : begin                ; or 'SJI_????' : begin
        berad = d.be_fsi / !radeg
        cdelt = d.cdlt_fsi
        hd[j].crpix1 = d.cpx1_fsi
        hd[j].crpix2 = d.cpx2_fsi
        fid_t0 = d.fid_279t0
        fid_t1 = d.fid_279t1
        fid_s1 = d.fid_279s1
        fid_s2 = d.fid_279s2
        fi1_t0 = d.fi1_279t0
        fi1_t1 = d.fi1_279t1
        fi1_s1 = d.fi1_279s1
        fi1_s2 = d.fi1_279s2
        o_bifid = 6.0
        o_bifi1 = 0.0
      end
      1 : begin                ; or 'SJI_1330' : begin
        berad = d.be_133 / !radeg
        cdelt = d.cdlt_133
        hd[j].crpix1 = d.cpx1_133
        hd[j].crpix2 = d.cpx2_133
        fid_t0 = d.fid_133t0
        fid_t1 = d.fid_133t1
        fid_s1 = d.fid_133s1
        fid_s2 = d.fid_133s2
        fi1_t0 = d.fi1_133t0
        fi1_t1 = d.fi1_133t1
        fi1_s1 = d.fi1_133s1
        fi1_s2 = d.fi1_133s2
        o_bifid =  1.7
        o_bifi1 = -2.3
      end
      2 : begin                ; or 'SJI_2796' : begin
        berad = d.be_279 / !radeg
        cdelt = d.cdlt_279
        hd[j].crpix1 = d.cpx1_279
        hd[j].crpix2 = d.cpx2_279
        fid_t0 = d.fid_279t0
        fid_t1 = d.fid_279t1
        fid_s1 = d.fid_279s1
        fid_s2 = d.fid_279s2
        fi1_t0 = d.fi1_279t0
        fi1_t1 = d.fi1_279t1
        fi1_s1 = d.fi1_279s1
        fi1_s2 = d.fi1_279s2
        o_bifid = 6.0
        o_bifi1 = 0.0
      end
      3 : begin                ; or 'SJI_1400' : begin
        berad = d.be_140 / !radeg
        cdelt = d.cdlt_140
        hd[j].crpix1 = d.cpx1_140
        hd[j].crpix2 = d.cpx2_140
        fid_t0 = d.fid_140t0
        fid_t1 = d.fid_140t1
        fid_s1 = d.fid_140s1
        fid_s2 = d.fid_140s2
        fi1_t0 = d.fi1_140t0
        fi1_t1 = d.fi1_140t1
        fi1_s1 = d.fi1_140s1
        fi1_s2 = d.fi1_140s2
        o_bifid =  1.7
        o_bifi1 = -2.3
      end
      4 : begin                ; or 'SJI_2832' : begin
        berad = d.be_283 / !radeg
        cdelt = d.cdlt_283
        hd[j].crpix1 = d.cpx1_283
        hd[j].crpix2 = d.cpx2_283
        fid_t0 = d.fid_283t0
        fid_t1 = d.fid_283t1
        fid_s1 = d.fid_283s1
        fid_s2 = d.fid_283s2
        fi1_t0 = d.fi1_283t0
        fi1_t1 = d.fi1_283t1
        fi1_s1 = d.fi1_283s1
        fi1_s2 = d.fi1_283s2
        o_bifid = 6.0
        o_bifi1 = 0.0
      end
      5 : begin                ; or 'SJI_????' : begin
        berad = d.be_mir / !radeg
        cdelt = d.cdlt_mir
        hd[j].crpix1 = d.cpx1_mir
        hd[j].crpix2 = d.cpx2_mir
        fid_t0 = d.fid_133t0
        fid_t1 = d.fid_133t1
        fid_s1 = d.fid_133s1
        fid_s2 = d.fid_133s2
        fi1_t0 = d.fi1_133t0
        fi1_t1 = d.fi1_133t1
        fi1_s1 = d.fi1_133s1
        fi1_s2 = d.fi1_133s2
        o_bifid =  1.7
        o_bifi1 = -2.3
      end
     endcase
     ; fix bimodal FW positions prior to flight software fix
     if t lt tref and hd[j].ifwenc eq hd[j].ifwpos-1 $
        then begin
        hd[j].crpix2 = hd[j].crpix2 + o_bifid
        hd[j].crpix1 = hd[j].crpix1 + o_bifi1
     endif
    endelse
    if tempsrc eq 0 then begin
       dcpx2=0.0
       dcpx1=0.0
    endif else begin
       dcpx2 = fid_t0 + total(fid_t1*tempi) + fid_s1*exp(fid_s2*(t-tref)/yr2s)
       dcpx1 = fi1_t0 + total(fi1_t1*tempi) + fi1_s1*exp(fi1_s2*(t-tref)/yr2s)
    endelse
    hd[j].crpix1  = (hd[j].crpix1 + dcpx1)  / sum1 + off1
    hd[j].crpix2  = (hd[j].crpix2 + dcpx2)  / sum2 + off2
    hd[j].crval1  = crvalxy[0]
    hd[j].crval2  = crvalxy[1]
    hd[j].cdelt1  = cdelt * sum1
    hd[j].cdelt2  = cdelt * sum2
    hd[j].ctype1  = 'HPLN-TAN'
    hd[j].ctype2  = 'HPLT-TAN'
    hd[j].cunit1  = 'arcsec  '
    hd[j].cunit2  = 'arcsec  '
    hd[j].pc1_1   =  cos(berad+garad)
    hd[j].pc1_2   = -sin(berad+garad) * sum2/sum1
    hd[j].pc2_1   =  sin(berad+garad) * sum1/sum2
    hd[j].pc2_2   =  cos(berad+garad)
    hd[j].crpix1a = hd[j].crpix1
    hd[j].crpix2a = hd[j].crpix2
    hd[j].crval1a = hd[j].crval1
    hd[j].crval2a = hd[j].crval2
    hd[j].cdelt1a = hd[j].cdelt1
    hd[j].cdelt2a = hd[j].cdelt2
    hd[j].ctype1a = hd[j].ctype1
    hd[j].ctype2a = hd[j].ctype2
    hd[j].cunit1a = hd[j].cunit1
    hd[j].cunit2a = hd[j].cunit2
    hd[j].pc1_1a  = hd[j].pc1_1
    hd[j].pc1_2a  = hd[j].pc1_2
    hd[j].pc2_1a  = hd[j].pc2_1
    hd[j].pc2_2a  = hd[j].pc2_2
    hd[j].lonpole = 180.0
    ; update xcen,ycen (this is not an endorsement of those keywords)
    if doxycen then begin
     hd[j].xcen = hd[j].crval1 + hd[j].cdelt1 * $
                  (hd[j].pc1_1*((hd[j].naxis1+1)/2.-hd[j].crpix1) $
                  +hd[j].pc1_2*((hd[j].naxis2+1)/2.-hd[j].crpix2))
     hd[j].ycen = hd[j].crval2 + hd[j].cdelt2 * $
                  (hd[j].pc2_1*((hd[j].naxis1+1)/2.-hd[j].crpix1) $
                  +hd[j].pc2_2*((hd[j].naxis2+1)/2.-hd[j].crpix2))
    endif
  end
 endcase

endfor

end