;+
; $Id: cor1_point.pro,v 1.12 2023/08/15 22:14:31 secchia Exp $
;
; Project     :	STEREO - SECCHI
;
; Name        :	COR1_POINT
;
; Purpose     :	Correct pointing values in COR1 headers.
;
; Category    :	STEREO, SECCHI, COR1, Calibration, Coordinates
;
; Explanation :	This routine applies the latest pointing calibration to a COR1
;               image header.  The following keywords are affected:
;
;                       CRPIX1, CRPIX2, CRPIX1A, CRPIX2A,
;                       CRVAL1, CRVAL2, CRVAL1A, CRVAL2A,
;                       CDELT1, CDELT2, CDELT1A, CDELT2A,
;                       PC1_1,  PC1_2,  PC2_1,   PC2_2,
;                       PC1_1A, PC1_2A, PC2_1A,  PC2_2A,
;                       CROTA, XCEN, YCEN, ATT_FILE, HISTORY
;
; Syntax      :	COR1_POINT, HDR
;
; Examples    :	A = SCCREADFITS(Filename, HDR)
;               COR1_POINT, HDR
;
; Inputs      : HDR     = SECCHI header structure, or array of structures.  Can
;                         also be passed as a FITS text header, in which case
;                         it will be converted into a SECCHI header structure.
;
; Opt. Inputs :	None.
;
; Outputs     :	HDR     = The modified header is returned in-place.
;
; Opt. Outputs:	None.
;
; Keywords    :	Accepts the keywords GTDBFILE, DIR, NO_GRH, and ERROR for the
;               SCC_SUNVEC routine.  See that routine for more details.
;
; Calls       :	DATATYPE, SCC_FITSHDR2STRUCT, TAG_NAMES, GET_STEREO_HPC_POINT,
;               GET_STEREO_ROLL, CONVERT_STEREO_LONLAT, FITSHEAD2WCS, SCC_SUNVEC,
;               WCS_GET_COORD, SCC_UPDATE_HISTORY, GET_STEREO_ATT_FILE
;
; Common      :	None.
;
; Restrictions:	None.
;
; Side effects:	None.
;
; Prev. Hist. :	None.
;
; History     :	
; $Log: cor1_point.pro,v $
; Revision 1.12  2023/08/15 22:14:31  secchia
; roll08/12 correction
;
; Revision 1.11  2015/07/01 14:57:25  thompson
; Modified for post-solar-conjunction
;
; Revision 1.10  2012/12/03 18:21:47  nathan
; adjust syntax for array input so it works
;
; Revision 1.9  2008/01/17 15:47:49  thompson
; Corrected bug with how GT data is incorporated into the pointing.
; This mostly affects Ahead data before 4-Feb-2007, and Behind data
; before 22-Feb-2007.
;
; Revision 1.8  2007/12/17 21:00:19  thompson
; Changed to use R* rectified position keywords instead of P* equivalents
;
; Revision 1.7  2007/08/30 18:47:12  nathan
; added INS_[xyr]0 (Bug 91); added CVS Log
;
;               Version 4, 14-May-2007, William Thompson, GSFC
;                       Update Behind parameters based on LT+S correction
;   	    	Version 3, 4-May-2007, N.Rich, NRL
;   	    	    	Update ATT_FILE and HISTORY
;               Version 2, 2-May-2007, William Thompson, GSFC
;                       Add calls to SCC_SUNVEC for pointing history, and
;                       optional jitter.
;   	    	Version 1, 3-Apr-2007, William Thompson, GSFC
;
; Contact     :	WTHOMPSON
;-
;
pro cor1_point, hdr, gtdbfile=gtdbfile, dir=telemdir, no_grh=no_grh, $
                error=error, SILENT=silent
;
;  Make sure that the header is a SECCHI structure, and that it belongs to
;  COR1.
;
if datatype(hdr) ne 'STC' then hdr=scc_fitshdr2struct(hdr)
if tag_names(hdr, /structure_name) ne 'SECCHI_HDR_STRUCT' then message, $
  'Only SECCHI header structure allowed'
w = where(hdr.detector ne 'COR1', count)
if count gt 0 then message, 'Calibration for COR1 detector only'

info="$Id: cor1_point.pro,v 1.12 2023/08/15 22:14:31 secchia Exp $"
IF ~keyword_set(SILENT) THEN print,info
;
;  Define the default parameters based on the observatory.
;
nhdr   = n_elements(hdr)
crpix  = fltarr(2,nhdr)         ;Reference pixel
crval  = fltarr(2,nhdr)         ;Reference value
cdelt  = fltarr(nhdr)           ;Plate scale
spoint = fltarr(3,nhdr)         ;Solar pointing (yaw, pitch, roll)
drot   = fltarr(nhdr)           ;Rotation correction
croll  = fltarr(nhdr)           ;Celestial roll
cpoint = fltarr(2,nhdr)         ;Celestial yaw and pitch
wa = where(hdr.obsrvtry eq 'STEREO_A', na, complement=wb, ncomplement=nb)
pcdate0 = '2015-05-19'          ;Start of post-conjunction
pcdate1 = '2023-08-12'          ;End of post-conjunction
;
;  Define the defaults for Ahead.  The pitch and yaw parts of GET_STEREO_ROLL
;  and GET_STEREO_HPC_POINT are ignored in favor of the Guide Telescope data.
;
if na gt 0 then begin
    crpix[0,wa] = 1070.39 + 1 ;occulter center
    crpix[1,wa] = 1025.58 + 3 ;occulter center
    crval[0,wa] = -32.899852
    crval[1,wa] =  97.096261
    cdelt[wa] = 3.75215
    drot[wa]   = 0.203826
    spoint[2,wa] = (get_stereo_hpc_point(hdr[wa].date_obs, 'A', /degrees, $
                                        /post_conjunction))[2,*]
    date_obs = hdr[wa].date_obs
    croll[wa] = get_stereo_roll(date_obs, 'A', system='GEI', /degrees, $
                               /post_conjunction)
    att_file = get_stereo_att_file(date_obs,'A')  ; much faster
;
;  Fold in the Guide Telescope data.
;
    date0 = '2007-03-01T17:10:35.374'
    dsun0 = (get_stereo_lonlat(date0, 'a', /au))[0]
    d0 = scc_sunvec(date0, dsun0, obs='a', /quiet, gtdbfile=gtdbfile)
    dsun = (get_stereo_lonlat(hdr[wa].date_obs, 'a', /au))[0,*]
    for ia=0,na-1 do begin
        ja = wa[ia]
        d =scc_sunvec(hdr[ja].date_obs, dsun[ia], obs='a', gtdbfile=gtdbfile, $
                      dir=telemdir, no_grh=no_grh, error=error, quiet=silent)
        if (hdr[ja].date_obs ge pcdate0) and (hdr[ja].date_obs lt pcdate1) $
                then sgn = 1 else sgn = -1
        spoint[0,ja] = sgn * (d[1]-d0[1]) / 3600.d0
        spoint[1,ja] = sgn * (d[0]-d0[0]) / 3600.d0
    endfor
;
    point = spoint[0:1,wa]
    convert_stereo_lonlat, date_obs, point, 'HPC', 'GEI', spacecraft='A', $
      /degrees
    cpoint[*,wa] = point
endif
;
;  Define the defaults for Behind.
;
if nb gt 0 then begin
    crpix[0,wb] = 1069.52 + 1 ;occulter center
    crpix[1,wb] = 1017.72 + 3 ;occulter center
    crval[0,wb] = 0.24846208
    crval[1,wb] = -135.95599
    cdelt[wb] = 3.75088
    drot[wb]  = 0.0983413
    spoint[2,wb] = (get_stereo_hpc_point(hdr[wb].date_obs, 'B', /degrees))[2,*]
    date_obs = hdr[wb].date_obs
    croll[wb] = get_stereo_roll(date_obs, 'B', system='GEI', /degrees)
    att_file = get_stereo_att_file(date_obs,'B')  ; much faster
;
;  Fold in the Guide Telescope data.
;
    date0 = '2007-02-25T03:33:52.383'
    dsun0 = (get_stereo_lonlat(date0, 'b', /au))[0]
    d0 = scc_sunvec(date0, dsun0, obs='b', /quiet, gtdbfile=gtdbfile)
    dsun = (get_stereo_lonlat(hdr[wb].date_obs, 'b', /au))[0,*]
    for ib=0,nb-1 do begin
        jb = wb[ib]
        d =scc_sunvec(hdr[jb].date_obs, dsun[ib], obs='b', gtdbfile=gtdbfile, $
                      dir=telemdir, no_grh=no_grh, error=error, quiet=silent)
        if (hdr[jb].date_obs ge pcdate0) and (hdr[jb].date_obs le pcdate1) $
                then sgn = -1 else sgn = 1
        spoint[0,jb] = sgn * (d[1]-d0[1]) / 3600.d0
        spoint[1,jb] = sgn * (d[0]-d0[0]) / 3600.d0
    endfor
;
    point = spoint[0:1,wb]
    convert_stereo_lonlat, date_obs, point, 'HPC', 'GEI', spacecraft='B', $
      /degrees
    cpoint[*,wb] = point
endif
;
; Adjust the default values for post-conjunction.
;
wdate = where((hdr.date_obs ge pcdate0) and (hdr.date_obs lt pcdate1), ndate)
if ndate gt 0 then crval[*,wdate] = -crval[*,wdate]
;
; Set instrument offset keywords
;
hdr.ins_x0 = reform(crval[0,*]) ; arcsec
hdr.ins_y0 = reform(crval[1,*]) ; arcsec
hdr.ins_r0 = drot   	; deg
;
;  Set SC_ values with GT values before transform
;
hdr.sc_yaw = reform(spoint[0,*])
hdr.sc_pitch=reform(spoint[1,*])
hdr.sc_yawa = reform(cpoint[0,*])
hdr.sc_pita = reform(cpoint[1,*])
;
;  Take into account the image rotation, if any.
;
wr = where(hdr.rectify eq 'T', countr)
if countr gt 0 then begin
    crpix0 = crpix
    wwa = where(hdr[wr].obsrvtry eq 'STEREO_A', nna, complement=wwb, $
                ncomplement=nnb)
    if nna gt 0 then begin
        wwa = wr[wwa]
        crpix[0,wwa] =        crpix0[1,wwa]
        crpix[1,wwa] = 2179 - crpix0[0,wwa]
    endif
    if nnb gt 0 then begin
        wwb = wr[wwb]
        crpix[0,wwb] = 2179 - crpix0[1,wwb]
        crpix[1,wwb] =        crpix0[0,wwb]
    endif
;
;  Adjust any post-conjunction data.
;
    wdate = where((hdr[wr].date_obs ge pcdate0) and $
                  (hdr[wr].date_obs lt pcdate1), ndate)
    if ndate gt 0 then begin
        wdate = wr[wdate]
        crpix[*,wdate] = 2179 - crpix[*,wdate]
    endif
endif
;
;  Take into account the area of the CCD used.
;  (If unrectified, then RCOL=PCOL)
;
crpix[0,*] = crpix[0,*] - hdr.r1col
crpix[1,*] = crpix[1,*] - hdr.r1row
;
;  Correct the plate scale and reference pixels for pixel summing.
;
summed = 2.^(hdr.summed - 1)    ;Amount of image summing
IF ~keyword_set(SILENT) THEN help,summed
w = where(summed gt 1, count)
if count gt 0 then begin
    cdelt[w] = cdelt[w] * summed[w]
    crpix[0,w] = (crpix[0,w] - 0.5) / summed[w] + 0.5
    crpix[1,w] = (crpix[1,w] - 0.5) / summed[w] + 0.5
endif
;
;  Step through the dates and transform the coordinates based on the spacecraft
;  pointing information.
;
dr = !dpi / 180.d0              ;Degrees to radians
ddr = dr / 3600.d0              ;Arcseconds to radians
;
for ihdr = 0,nhdr-1 do begin
;
;  Form the transformation matrix from the sines and cosines of the spacecraft
;  sun-pointing parameters.
;
    sy = sin(spoint[0,ihdr]*dr)  &  cy = cos(spoint[0,ihdr]*dr)
    sp = sin(spoint[1,ihdr]*dr)  &  cp = cos(spoint[1,ihdr]*dr)
    sr = sin(spoint[2,ihdr]*dr)  &  cr = cos(spoint[2,ihdr]*dr)
;
    mat = [[cp*cy, cp*sy, -sp], $
           [-cr*sy+sr*sp*cy,  cr*cy+sr*sp*sy, sr*cp], $
           [ sr*sy+cr*sp*cy, -sr*cy+cr*sp*sy, cr*cp]]
;
;  Form the vector representing the pointing of the instrument.
;
    yaw   = -crval[0,ihdr] * ddr
    pitch =  crval[1,ihdr] * ddr
    vec = [cos(yaw)*cos(pitch), sin(yaw)*cos(pitch), sin(pitch)]
;
;  Transform the pointing and recalculate the pointing.
;
    newvec = mat ## vec
    crval[0,ihdr] = -atan(newvec[1], newvec[0]) / ddr
    crval[1,ihdr] =  asin(newvec[2])            / ddr
;
;  Convert the pointing into celestial coordinates.
;
    point = crval[0:1,ihdr] / 3600.d0
    convert_stereo_lonlat, hdr[ihdr].date_obs, point, 'HPC', 'GEI', /degrees, $
      spacecraft=hdr[ihdr].obsrvtry
    cpoint[0:1,ihdr] = point
endfor
;
;  Update the information in the header.
;
hdr.crpix1  = reform(crpix[0,*])
hdr.crpix1a = hdr.crpix1
hdr.crpix2  = reform(crpix[1,*])
hdr.crpix2a = hdr.crpix2
;
hdr.cdelt1  = cdelt
hdr.cdelt2  = cdelt
ccdelt = cdelt / 3600.0         ;Convert arcseconds to degrees
hdr.cdelt1a = -ccdelt           ;Negative for celestial coordinates
hdr.cdelt2a =  ccdelt
;
hdr.crval1  = reform(crval[0,*])
hdr.crval2  = reform(crval[1,*])
hdr.crval1a = reform(cpoint[0,*])
hdr.crval2a = reform(cpoint[1,*])
;
;  Add either 90 or 270 degrees to the rotation angle for unrectified images.
;
gamma = reform(spoint[2,*]) + drot
wa = where((hdr.rectify ne 'T') and (hdr.obsrvtry eq 'STEREO_A'), counta)
if counta gt 0 then gamma[wa] = gamma[wa] + 270
wb = where((hdr.rectify ne 'T') and (hdr.obsrvtry eq 'STEREO_B'), countb)
if countb gt 0 then gamma[wb] = gamma[wb] + 90
;
;  Put the correct rotation angle and matrices in the headers.
;
hdr.crota = gamma
gamma = gamma * dr              ;Convert degrees to radians
cgamma = cos(gamma)
sgamma = sin(gamma)
hdr.pc1_1 =  cgamma
hdr.pc1_2 = -sgamma
hdr.pc2_1 =  sgamma
hdr.pc2_2 =  cgamma
;
;  Calculate XCEN, YCEN
;
di = (hdr.naxis1 + 1.)/2. - reform(crpix[0,*])
dj = (hdr.naxis2 + 1.)/2. - reform(crpix[1,*])
hdr.xcen = reform(crval[0,*]) + cdelt*(di*cgamma - dj*sgamma)
hdr.ycen = reform(crval[1,*]) + cdelt*(di*sgamma + dj*cgamma)
;
;  Do the same for celestial coordinates.
;
gamma = (croll - drot) * dr
cgamma = cos(gamma)
sgamma = sin(gamma)
hdr.pc1_1a =  cgamma
hdr.pc1_2a =  sgamma
hdr.pc2_1a = -sgamma
hdr.pc2_2a =  cgamma

IF error GT 0 THEN errstr=trim(fix(error)) ELSE errstr=''
att_file = att_file+'+'+errstr+'GT'
; update att_file
if tag_exist(hdr,'att_file') then hdr.att_file = att_file
if ~keyword_set(silent) then help,error,att_file

len=strlen(info)
histinfo=strmid(info,1,len-2)
if tag_exist(hdr,'history') then hdr=scc_update_history(hdr, histinfo)

return
end