;+
;
; Project: STEREO-SECCHI
;
; Name: hi_hae_point
;
; Purpose: Calculate the HI yaw pitch and roll in the HAE coordinate system
;
; Category: STEREO, SECCHI, HI, Coordinates
;
; Explanation: 
;
; Syntax: hi_hae_point,index
;
; Example: roll = hi_hae_point(index,yaw=yaw,pitch=pitch)
;
; Inputs: index - one HI header structure
;
; Opt. Inputs: None
;
; Outputs: index - one HI header structure containing calibrated
;                  pointing information
;
; Opt. Outputs: yaw
;               pitch
;
; Keywords: None
;
; Calls: get_hi_params
;
; Common: None
;
; Restrictions: None 
;
; Side effects: None
;
; Prev. Hist.: None
;
; History: 20-Dec-2007, Danielle Bewsher
;
; Contact: Danielle Bewsher (d.bewsher@rl.ac.uk)
;
;-
FUNCTION hi_hae_point,index,yaw=yaw,pitch=pitch

;determine the spacecraft HAE roll angle
hae_roll = get_stereo_roll(index.date_obs,index.obsrvtry,hae_yaw,hae_pitch,system='HAE')

;This part is all from hi_calib_roll
;----------------------------------

xv = [0.5*float(index.naxis1),float(index.naxis1)]
yv = [0.5*float(index.naxis2),0.5*float(index.naxis2)]

ccdosx=0.
ccdosy=0.
pmult=float(index.naxis1)/2.0

get_hi_params,index,pitch_hi,offset_hi,roll_hi,mu,d,_extra=extra

ang = (90.-0.5*d)*!dpi/180.
rng = (1.0+mu)*cos(ang)/(sin(ang)+mu)

vfov = transpose([[reform(xv)],[reform(yv)]])

nst = n_elements(vfov(0,*))

vv4 = dblarr(4,nst)
vv4(0,*) = ((vfov(0,*)-ccdosx)/pmult-1.0)*rng
vv4(1,*) = ((vfov(1,*)-ccdosy)/pmult - 1.0)*rng
vv4(2,*) = 1.0
vv4(3,*) = 1.0

vv3=azp2cart(vv4,mu)
vv3(2,*) = 1.0

vv2=hi2sc(vv3,roll_hi,pitch_hi,offset_hi)

xy=sc2cart(vv2,hae_roll,hae_pitch,hae_yaw)

;generate vector of transformed points
z=xy(2,1)-xy(2,0)
x=-(xy(1,1)-xy(1,0))
y=(xy(0,1)-xy(0,0))

; generate comparative tangent of circle at fixed height
tx=xy(0,0)
ty=xy(1,0)
tz=0.0

; use inverted dot product to calculate angle
a = sqrt(tx^2 + ty^2)
b = sqrt(x^2 + y^2 + z^2)
ab = x*tx + y*ty

; noting that ambiguity is resolved with up/down direction of vector
val=min([max([ab/(a*b),-1.0]),1.0])
if (z ge 0.0) then begin
   oroll=-acos(val)*180./!dpi
endif else begin
   oroll=acos(val)*180./!dpi
endelse

;------------------------------------

; sort out pointing in RA/DEC
xv=[0.5*float(index.naxis1),float(index.naxis1)]
yv=[0.5*float(index.naxis2),float(index.naxis2)]

;this part is all from hi_calib_point
;-----------------------------------
vfov=transpose([[reform(xv)],[reform(yv)]])

nst=n_elements(vfov(0,*))

vv4=fltarr(4,nst)
vv4(0,*) = ((vfov(0,*)-ccdosx)/pmult - 1.0)*rng
vv4(1,*) = ((vfov(1,*)-ccdosy)/pmult - 1.0)*rng
vv4(3,*) = 1.0

vv3=azp2cart(vv4,mu)
vv2=hi2sc(vv3,roll_hi,pitch_hi,offset_hi)
vv=sc2cart(vv2,hae_roll,hae_pitch,hae_yaw)

radec = fltarr(2,nst)

for i=0,nst-1 do begin
    th=acos(vv(2,i))
    radec(1,i) = 90. - th*!radeg

    cphi=vv(0,i)/sin(th)
    sphi=vv(1,i)/sin(th)

    th1=acos(cphi)*!radeg
    th2=asin(sphi)*!radeg
    if (th2 gt 0) then begin
        ra = th1
    endif else begin
        ra = -th1
    endelse
    radec(0,i) = ra
endfor
;--------------------------------

yaw = radec(0,0)
pitch = radec(1,0)

return,oroll

END