;+
; NAME:
;        HXA2HXAXY
; PURPOSE:
;        Tries to reconstruct HXA limb positions hidden by fiducial marks
;        and calculate the suncenter position (in HXA coordinates) from
;        the HXA info in the PNT files.
; CATEGORY:
; CALLING SEQUENCE:
;        sunc = hxa2hxaxy(hxa_i,hxa_d)
;        sunc = hxa2hxaxy(atr,atr.hxa)
; INPUTS:
;        hxa_i = pnt, atr, or ada_index structure.
;        hxa_d = array(4,*). hxa data, e.g., atr.hxa or gt_hxa(ada_data)
; KEYWORDS (INPUT):
; OUTPUTS:
;        sunc = array(2,*). suncenter coordinates in HXA coordinates.
;               sunc(0,*) = hxa_x, sunc(1,*) = hxa_y
; KEYWORDS (OUTPUT):
;        code = vector.  Indicates how each sun center position was obtained.
;               0 : error, no result, 2 : result deduced from 2 limbs,
;               3: result deduced from 3 limbs, 4 : all limbs used.
;        zon = vector, the 4 LSB indicate which limbs were hidden.
; COMMON BLOCKS:
;        None.
; SIDE EFFECTS:
; RESTRICTIONS:
; PROCEDURE:
;        Uses the solar radius to reconstruct hidden limbs.
; MODIFICATION HISTORY:
;	19-Jul-93 (JPW) - part of hxa_suncenter
;       16-Sep-96 (JPW) - modif. to allow for diff. solar r in x and y
;                         (to correct for diff. sensitiv. of HXA x and y)
;-

function hxa2hxaxy,hxa_i,hxa_d,code=code,zon=zon ;,hxa_out=hxa

;hdb_file = '$ys/ucon/idl/wuelser/hxa_dbase.genx'
hdb_file = concat_dir('DIR_GEN_STATUS','hxa_dbase.genx')

; various constants
hxa_parms, hpix=hpix, h_xx=h_xx, h_yy=h_yy
hpix = double(hpix)    ; hpix = HXA pixel size
h_xx = double(h_xx)    ; h_xx, h_yy = HXA detector crossing point coordinates
h_yy = double(h_yy)
maxerr = 55.0d    ; error margin for tests after limb reconstruction,
minerr = -5.0d    ; assumes that fiducial marks are narrower than 50 pixels
finerr = 3.0d    ; error margin for final test

; prepare variables to resolve for hidden limbs
n = n_elements(hxa_i)
sunc = dblarr(2,n)
code = intarr(n)
zon = intarr(n)   ; commented out for testing only, put back in!!
hxa = (hxa_d > 0) < 2047
sunr = get_rb0p(hxa_i,/radius)/hpix          ; solar radius in hxa pixels

;fx = findgen(n)+1.0   ; apparent sunr_x = fx * sunr and sunr_y = fy * sunr
;fy = findgen(n)+1.0   ; to be replaced with values from a database
restgen,hdb,file=hdb_file
hdbs = int2secarr(hdb,hdb(0))
hxas = int2secarr(hxa_i,hdb(0))
fx = interpol(hdb.hxa_fx*1e-4,hdbs,hxas)
fy = interpol(hdb.hxa_fy*1e-4,hdbs,hxas)

; create bit pattern of hidden limbs.
; commented out for testing only:  Put back in!!
forzones,xzon,yzon
whe = where(xzon(hxa(0,*)))
if whe(0) ge 0 then zon(whe) = zon(whe) or 1
whe = where(xzon(hxa(1,*)))
if whe(0) ge 0 then zon(whe) = zon(whe) or 2
whe = where(yzon(hxa(2,*)))
if whe(0) ge 0 then zon(whe) = zon(whe) or 4
whe = where(yzon(hxa(3,*)))
if whe(0) ge 0 then zon(whe) = zon(whe) or 8

hxa = double(hxa)
hxa(0:1,*) = hxa(0:1,*) - h_xx          ; 0 at crossing point
hxa(2:3,*) = hxa(2:3,*) - h_yy
sunc(0,*) = (hxa(0,*)+hxa(1,*))/2.0d
sunc(1,*) = (hxa(2,*)+hxa(3,*))/2.0d

; go through the various possible cases of hidden/nonhidden limbs

; 4 limbs visible : no fix required
whe = where(zon eq 0)
if whe(0) ge 0 then code(whe) = 4

; 3 limbs visible.  4 cases
whe = where(zon eq 1)
if whe(0) ge 0 then begin
   ; fix hxa0
;   aux = sunr(whe)*sunr(whe) - sunc(1,whe)*sunc(1,whe)
   aux = fx(whe)*fx(whe)*sunr(whe)*sunr(whe) - sunc(1,whe)*sunc(1,whe)
   ; test if root nonzero
   ww = where(aux ge 0.0d)
   if ww(0) ge 0 then begin
      whe = whe(ww)
      hxa0 = hxa(1,whe) - 2.0d*sqrt(aux(ww))
      ; test if not too far from hidden limb (fid. marks about 50 pixel wide)
      err1 = hxa(0,whe)-hxa0
      ww = where((err1 le maxerr) and (err1 ge minerr))
      if ww(0) ge 0 then begin
         whe = whe(ww)
         hxa(0,whe) = hxa0(ww)
         sunc(0,whe) = (hxa(0,whe)+hxa(1,whe))/2.0d
         code(whe) = 3
      endif
   endif
endif
whe = where(zon eq 2)
if whe(0) ge 0 then begin
   ; fix hxa1
;   aux = sunr(whe)*sunr(whe) - sunc(1,whe)*sunc(1,whe)
   aux = fx(whe)*fx(whe)*sunr(whe)*sunr(whe) - sunc(1,whe)*sunc(1,whe)
   ; test if root nonzero
   ww = where(aux ge 0.0d)
   if ww(0) ge 0 then begin
      whe = whe(ww)
      hxa1 = hxa(0,whe) + 2.0d*sqrt(aux(ww))
      ; test if not too far from hidden limb (fid. marks about 50 pixel wide)
      err1 = hxa1-hxa(1,whe)
      ww = where((err1 le maxerr) and (err1 ge minerr))
      if ww(0) ge 0 then begin
         whe = whe(ww)
         hxa(1,whe) = hxa1(ww)
         sunc(0,whe) = (hxa(0,whe)+hxa(1,whe))/2.0d
         code(whe) = 3
      endif
   endif
endif
whe = where(zon eq 4)
if whe(0) ge 0 then begin
   ; fix hxa2
;   aux = sunr(whe)*sunr(whe) - sunc(0,whe)*sunc(0,whe)
   aux = fy(whe)*fy(whe)*sunr(whe)*sunr(whe) - sunc(0,whe)*sunc(0,whe)
   ; test if root nonzero
   ww = where(aux ge 0.0d)
   if ww(0) ge 0 then begin
      whe = whe(ww)
      hxa2 = hxa(3,whe) - 2.0d*sqrt(aux(ww))
      ; test if not too far from hidden limb (fid. marks about 50 pixel wide)
      err1 = hxa(2,whe)-hxa2
      ww = where((err1 le maxerr) and (err1 ge minerr))
      if ww(0) ge 0 then begin
         whe = whe(ww)
         hxa(2,whe) = hxa2(ww)
         sunc(1,whe) = (hxa(2,whe)+hxa(3,whe))/2.0d
         code(whe) = 3
      endif
   endif
endif
whe = where(zon eq 8)
if whe(0) ge 0 then begin
   ; fix hxa3
;   aux = sunr(whe)*sunr(whe) - sunc(0,whe)*sunc(0,whe)
   aux = fy(whe)*fy(whe)*sunr(whe)*sunr(whe) - sunc(0,whe)*sunc(0,whe)
   ; test if root nonzero
   ww = where(aux ge 0.0d)
   if ww(0) ge 0 then begin
      whe = whe(ww)
      hxa3 = hxa(2,whe) + 2.0d*sqrt(aux(ww))
      ; test if not too far from hidden limb (fid. marks about 50 pixel wide)
      err1 = hxa3-hxa(3,whe)
      ww = where((err1 le maxerr) and (err1 ge minerr))
      if ww(0) ge 0 then begin
         whe = whe(ww)
         hxa(3,whe) = hxa3(ww)
         sunc(1,whe) = (hxa(2,whe)+hxa(3,whe))/2.0d
         code(whe) = 3
      endif
   endif
endif

; 2 limbs visible.  4 cases
; assumes that hxa0 hxa2 always neg. and hxa1 and hxa3 always pos.
; since solution not unique
; correction for different sensitivity, and thus different apparent
; solar radii for hxa_x and hxa_y makes expressions somewhat complicated
whe = where(zon eq 5)
if whe(0) ge 0 then begin
   ; fix hxa0 and hxa2
   ww = where((hxa(1,whe) gt 1.0d) and (hxa(3,whe) gt 1.0d))
   if ww(0) ge 0 then begin
      whe = whe(ww)
      hxa1=hxa(1,whe)*hxa(1,whe)
      hxa3=hxa(3,whe)*hxa(3,whe)
      srad=sunr(whe)*sunr(whe)
;      hxa0 = hxa3*(4.0d*srad-hxa1-hxa3)/(hxa1+hxa3)
;      hxa2 = hxa1*(4.0d*srad-hxa1-hxa3)/(hxa1+hxa3)
      fxrad=fx(whe)*fx(whe)
      fyrad=fy(whe)*fy(whe)
      aa = hxa1+hxa3
      bb = (fxrad-fyrad)*srad
      cc = 2d0*(fxrad+fyrad)*srad
      b0 = 2d0*bb*hxa(1,whe)
      c0 = hxa3*(aa-cc)+bb*bb
      hxa0 = b0*b0-4d0*aa*c0
      b2 = -2d0*bb*hxa(3,whe)
      c2 = hxa1*(aa-cc)+bb*bb
      hxa2 = b2*b2-4d0*aa*c2
      ww = where((hxa0 ge 0.0d) and (hxa2 ge 0.0d))
      if ww(0) ge 0 then begin
         whe = whe(ww)
;         hxa0 = -sqrt(hxa0(ww))
;         hxa2 = -sqrt(hxa2(ww))
         hxa0 = (-b0(ww) - sqrt(hxa0(ww))) / (2d0*aa(ww))
         hxa2 = (-b2(ww) - sqrt(hxa2(ww))) / (2d0*aa(ww))
         ; test if not too far from hidden limb (fid. marks < 50 pixel wide)
         err1 = hxa(0,whe)-hxa0
         err2 = hxa(2,whe)-hxa2
         ww = where((err1 le maxerr) and (err1 ge minerr) and $
                    (err2 le maxerr) and (err2 ge minerr))
         if ww(0) ge 0 then begin
            whe = whe(ww)
            hxa(0,whe) = hxa0(ww)
            hxa(2,whe) = hxa2(ww)
            sunc(0,whe) = (hxa(0,whe)+hxa(1,whe))/2.0d
            sunc(1,whe) = (hxa(2,whe)+hxa(3,whe))/2.0d
            code(whe) = 2
         endif
      endif
   endif
endif
whe = where(zon eq 6)
if whe(0) ge 0 then begin
   ; fix hxa1 and hxa2
   ww = where((hxa(0,whe) lt -1.0d) and (hxa(3,whe) gt 1.0d))
   if ww(0) ge 0 then begin
      whe = whe(ww)
      hxa0=hxa(0,whe)*hxa(0,whe)
      hxa3=hxa(3,whe)*hxa(3,whe)
      srad=sunr(whe)*sunr(whe)
;      hxa1 = hxa3*(4.0d*srad-hxa0-hxa3)/(hxa0+hxa3)
;      hxa2 = hxa0*(4.0d*srad-hxa0-hxa3)/(hxa0+hxa3)
      fxrad=fx(whe)*fx(whe)
      fyrad=fy(whe)*fy(whe)
      aa = hxa0+hxa3
      bb = (fxrad-fyrad)*srad
      cc = 2d0*(fxrad+fyrad)*srad
      b1 = 2d0*bb*hxa(0,whe)
      c1 = hxa3*(aa-cc)+bb*bb
      hxa1 = b1*b1-4d0*aa*c1
      b2 = -2d0*bb*hxa(3,whe)
      c2 = hxa0*(aa-cc)+bb*bb
      hxa2 = b2*b2-4d0*aa*c2
      ww = where((hxa1 ge 0.0d) and (hxa2 ge 0.0d))
      if ww(0) ge 0 then begin
         whe = whe(ww)
;         hxa1 = +sqrt(hxa1(ww))
;         hxa2 = -sqrt(hxa2(ww))
         hxa1 = (-b1(ww) + sqrt(hxa1(ww))) / (2d0*aa(ww))
         hxa2 = (-b2(ww) - sqrt(hxa2(ww))) / (2d0*aa(ww))
         ; test if not too far from hidden limb (fid. marks < 50 pixel wide)
         err1 = hxa1-hxa(1,whe)
         err2 = hxa(2,whe)-hxa2
         ww = where((err1 le maxerr) and (err1 ge minerr) and $
                    (err2 le maxerr) and (err2 ge minerr))
         if ww(0) ge 0 then begin
            whe = whe(ww)
            hxa(1,whe) = hxa1(ww)
            hxa(2,whe) = hxa2(ww)
            sunc(0,whe) = (hxa(0,whe)+hxa(1,whe))/2.0d
            sunc(1,whe) = (hxa(2,whe)+hxa(3,whe))/2.0d
            code(whe) = 2
         endif
      endif
   endif
endif
whe = where(zon eq 9)
if whe(0) ge 0 then begin
   ; fix hxa0 and hxa3
   ww = where((hxa(1,whe) gt 1.0d) and (hxa(2,whe) lt -1.0d))
   if ww(0) ge 0 then begin
      whe = whe(ww)
      hxa1=hxa(1,whe)*hxa(1,whe)
      hxa2=hxa(2,whe)*hxa(2,whe)
      srad=sunr(whe)*sunr(whe)
;      hxa0 = hxa2*(4.0d*srad-hxa1-hxa2)/(hxa1+hxa2)
;      hxa3 = hxa1*(4.0d*srad-hxa1-hxa2)/(hxa1+hxa2)
      fxrad=fx(whe)*fx(whe)
      fyrad=fy(whe)*fy(whe)
      aa = hxa1+hxa2
      bb = (fxrad-fyrad)*srad
      cc = 2d0*(fxrad+fyrad)*srad
      b0 = 2d0*bb*hxa(1,whe)
      c0 = hxa2*(aa-cc)+bb*bb
      hxa0 = b0*b0-4d0*aa*c0
      b3 = -2d0*bb*hxa(2,whe)
      c3 = hxa1*(aa-cc)+bb*bb
      hxa3 = b3*b3-4d0*aa*c3
      ww = where((hxa0 ge 0.0d) and (hxa3 ge 0.0d))
      if ww(0) ge 0 then begin
         whe = whe(ww)
;         hxa0 = -sqrt(hxa0(ww))
;         hxa3 = +sqrt(hxa3(ww))
         hxa0 = (-b0(ww) - sqrt(hxa0(ww))) / (2d0*aa(ww))
         hxa3 = (-b3(ww) + sqrt(hxa3(ww))) / (2d0*aa(ww))
         ; test if not too far from hidden limb (fid. marks < 50 pixel wide)
         err1 = hxa(0,whe)-hxa0
         err2 = hxa3-hxa(3,whe)
         ww = where((err1 le maxerr) and (err1 ge minerr) and $
                    (err2 le maxerr) and (err2 ge minerr))
         if ww(0) ge 0 then begin
            whe = whe(ww)
            hxa(0,whe) = hxa0(ww)
            hxa(3,whe) = hxa3(ww)
            sunc(0,whe) = (hxa(0,whe)+hxa(1,whe))/2.0d
            sunc(1,whe) = (hxa(2,whe)+hxa(3,whe))/2.0d
            code(whe) = 2
         endif
      endif
   endif
endif
whe = where(zon eq 10)
if whe(0) ge 0 then begin
   ; fix hxa1 and hxa3
   ww = where((hxa(0,whe) lt -1.0d) and (hxa(2,whe) lt -1.0d))
   if ww(0) ge 0 then begin
      whe = whe(ww)
      hxa0=hxa(0,whe)*hxa(0,whe)
      hxa2=hxa(2,whe)*hxa(2,whe)
      srad=sunr(whe)*sunr(whe)
;      hxa1 = hxa2*(4.0d*srad-hxa0-hxa2)/(hxa0+hxa2)
;      hxa3 = hxa0*(4.0d*srad-hxa0-hxa2)/(hxa0+hxa2)
      fxrad=fx(whe)*fx(whe)
      fyrad=fy(whe)*fy(whe)
      aa = hxa0+hxa2
      bb = (fxrad-fyrad)*srad
      cc = 2d0*(fxrad+fyrad)*srad
      b1 = 2d0*bb*hxa(0,whe)
      c1 = hxa2*(aa-cc)+bb*bb
      hxa1 = b1*b1-4d0*aa*c1
      b3 = -2d0*bb*hxa(2,whe)
      c3 = hxa0*(aa-cc)+bb*bb
      hxa3 = b3*b3-4d0*aa*c3
      ww = where((hxa1 ge 0.0d) and (hxa3 ge 0.0d))
      if ww(0) ge 0 then begin
         whe = whe(ww)
;         hxa1 = +sqrt(hxa1(ww))
;         hxa3 = +sqrt(hxa3(ww))
         hxa1 = (-b1(ww) + sqrt(hxa1(ww))) / (2d0*aa(ww))
         hxa3 = (-b3(ww) + sqrt(hxa3(ww))) / (2d0*aa(ww))
         ; test if not too far from hidden limb (fid. marks < 50 pixel wide)
         err1 = hxa1-hxa(1,whe)
         err2 = hxa3-hxa(3,whe)
         ww = where((err1 le maxerr) and (err1 ge minerr) and $
                    (err2 le maxerr) and (err2 ge minerr))
         if ww(0) ge 0 then begin
            whe = whe(ww)
            hxa(1,whe) = hxa1(ww)
            hxa(3,whe) = hxa3(ww)
            sunc(0,whe) = (hxa(0,whe)+hxa(1,whe))/2.0d
            sunc(1,whe) = (hxa(2,whe)+hxa(3,whe))/2.0d
            code(whe) = 2
         endif
      endif
   endif
endif

; do final testing of the reconstructed values.  Helps avoid garbage results
; in case of bad detector values or ambiguous reconstructions.
whe = where((zon/4 ne 3) and (zon mod 4 ne 3))
if whe(0) ge 0 then begin
   xoff = sunc(0,whe)
   yoff = sunc(1,whe)
   xrad = (hxa(1,whe)-hxa(0,whe))/2.0d
   yrad = (hxa(3,whe)-hxa(2,whe))/2.0d
;   err1 = abs(sqrt(xrad*xrad + yoff*yoff) - sunr(whe))
;   err2 = abs(sqrt(yrad*yrad + xoff*xoff) - sunr(whe))
   err1 = abs(sqrt(xrad*xrad + yoff*yoff) - fx(whe)*sunr(whe))
   err2 = abs(sqrt(yrad*yrad + xoff*xoff) - fy(whe)*sunr(whe))
   werr = where((err1 gt finerr) or (err2 gt finerr))
   if werr(0) ge 0 then code(whe(werr)) = 0    ; set error code
endif

whe = where(code eq 0)
if whe(0) ge 0 then sunc(*,whe) = 0.0d

sunc(0,*) = sunc(0,*) + h_xx    ; re-add offsets
sunc(1,*) = sunc(1,*) + h_yy
sunc = float(sunc)

; for testing:
;hxa(0:1,*) = hxa(0:1,*) + h_xx
;hxa(2:3,*) = hxa(2:3,*) + h_yy

return,sunc
end