;+
; NAME:
;        sb_p4display
; PURPOSE:
;        Calculate position angle of solar N relative to
;        direction to sbrowser display UP.
;        sbrowser display shall have its axis1 parallel to the
;        sun-stereoA-stereoB plane.  Initial version assumes STEREO
;        to be in the ecliptic plane.  May need to be revisited to
;        include actual inclinations of both S/C orbits
; CATEGORY:
;        image processing
; CALLING SEQUENCE:
;        p_angle = sb_p4display(t=t,h=h)
;        p_angle = sb_p4display(owcs=owcs)
;        p_angle = sb_p4display(owcs=owcs,swcs=swcs,/sc)
; INPUTS:
; KEYWORDS (INPUT):
;        time = time.  "anytim" compatible format
;        hgln_obs = heliographic longitude of observer
;        /sc : return position angle of solar N from normal to SC plane.
;              Default is ecliptic plane.  Requires o_wcs, s_wcs
;        owcs = WCS of observer
;        swcs = WCS of other SC defining observer plane
;            Note: WCSs must have position.hgln_obs tag
;        /helio : return 0.
; OUTPUTS:
;        p_angle = position angle of Solar N ccw from display axis2
; KEYWORDS (OUTPUT):
; COMMON BLOCKS:
; SIDE EFFECTS:
; RESTRICTIONS:
;        Accuracy limited to about 1 arcmin
;        Assumes STEREO in ecliptic plane unless /sc set
; PROCEDURE:
;        Simplified algorithms from sun_pos.pro and pb0r.pro
; MODIFICATION HISTORY:
;        JPW, 4/11/2006
;        JPW, 2/27/2007 added /sc and /helio options
;-

FUNCTION sb_p4display,time=tin,hgln_obs=hgln,helio=helio,sc=sc, $
                      owcs=owcs,swcs=swcs

if keyword_set(helio) then return,0d0

; sc plane:
if keyword_set(sc) and n_elements(owcs) gt 0 and n_elements(swcs) gt 0 $
 then begin
  ; vector to other SC (z-axis pointed to observer, y-axis towards solar N)
  hg = [swcs.position.hgln_obs,swcs.position.solar_b0]
  v2 = sb_hg2hc(hg,pos=owcs.position,/earth)  ; /earth since hgln rel. earth
  ; normal to SC plane v1 x v2, assuming v1 = [0,0,1] = vector to observer
  vn = [-v2[1],v2[0],0d0]
  if vn[1] lt 0.0 then vn=-vn              ; limit pos.angle to +/-90 degrees
  if sqrt(total(vn*vn)) gt 1e-3 then begin ; min. separation 1 mrad)
     p = atan(vn[0],vn[1]) * !radeg        ; pos.angle solar-N ccw from normal
     return,p
  endif
endif

; ecliptic plane
if n_elements(tin) eq 0 then tin=wcs_get_time(owcs)
if n_elements(hgln) eq 0 then hgln=owcs.position.hgln_obs

;---------------------------------------------------------------------------
;  number of Julian days since 2415020.0
;  ignoring difference between UT and ET
;  JPW approach
;---------------------------------------------------------------------------

mjd = anytim(tin,/mjd)
de  = double(mjd.mjd)+double(mjd.time)/(24d0*3.6d6)-15019.5d0

;---------------------------------------------------------------------------
;  get the longitude of the sun
;---------------------------------------------------------------------------
; from sun_pos.pro

;  form time in Julian centuries from 1900.0
   t = de/36525.0d0

;  form sun's mean longitude
   l = (279.696678d0+((36000.768925d0*t) mod 360.0d0))*3600.0d0

;  allow for ellipticity of the orbit (equation of centre)
;  using the Earth's mean anomoly ME
   me = 358.475844d0 + ((35999.049750D0*t) mod 360.0d0)
   ellcor  = (6910.1d0 - 17.2D0*t)*sin(me*!dtor) + 72.3D0*sin(2.0D0*me*!dtor)
   l = l + ellcor

; JPW: removed disturbances since << 1 arcmin (error << 0.2 arcsec on sun)

   longmed = l/3600.0d0                                         

;---------------------------------------------------------------------------
; get the p-angle relative to ecliptic N
;---------------------------------------------------------------------------
; from pb0r.pro , but without inclination of earth's axis

;  form aberrated longitude
   lambda = longmed - (20.5D0/3600.0D0)

;  form longitude of ascending node of sun's equator on ecliptic
   node = 73.666666D0 + (50.25D0/3600.0D0)*( (de/365.25d0) + 50.0d0 )
   arg = lambda - node

; JPW: add heliographic longitude of observer
   arg = arg + hgln

;  calculate P, the position angle of the pole (jpw:ccw from eclipt. pole)
   p = (ATAN( -0.12722D0 * COS(arg*!dtor))) * !radeg

;  ... and B0 the tilt of the axis
;   b = ASIN( 0.12620D0 * SIN(arg*!dtor) ) * !radeg

RETURN, p

END