pro norh_t2azel,time,az,el,delta=delta,hour_angle=hour_angle

;+
; NAME:
;       norh_t2azel
;
; PURPOSE:
;       Return the elevation and azimuth of the Sun at NoRH
;
; CALLING SEQUENCE:
;       norh_t2azel,time,az,el
;
; INPUTS:
;       time: UT time
; OUTPUTS:
;       el: elevation of the Sun at NoRH (in unit of arc-degree)
;       az: azimuth of the Sun at NoRH (in unit of arc-degree)
;
; OPTIONAL OUTPUTS:
;       delta: solar declination (in unit of arc-degree)
;       hour_angle: hour angle (in unit of arc-sec)
;
; HISTORY:
;       Modified 2000-09-10 TY
;       Modified 2000-12-20 TY added 180 to az
;
; Reference:
;       Christiansen & Hoegbom (1969) "Radiotelescopes" Appendix 4
;-

if not keyword_set(ephdir) then ephdir=getenv('DIR_NORH_EPH')

deg_per_sec=360./60./60./24.

nx=n_elements(time)
delta=fltarr(nx)
hour_angle=fltarr(nx)

timej=norh_ut2jst(time)

yymmdd=norh_t2ymd(timej,yy=yy,mm=mm,dd=dd)

if (n_elements(yy) ge 2) then yy0=yy(uniq(yy)) else yy0=yy

file=ephdir+'/track_table'+yy0+'.dat'
text=rd_tfile(file,/autocol,hskip=4)
mmdd=reform(text(0,*))

for n=0,nx-1 do begin

utime_time=anytim(time(n))
datej=(anytim(timej(n),/ints)).day
time_datejnoon= norh_jst2ut([3600000L*12,datej])
utime_noon=anytim(time_datejnoon)

whr=where(mmdd eq mm(n)+dd(n),cnt)
if (cnt eq 0) then begin
  print,'No data for ',time(n)
  goto,skip
endif

if (utime_time lt utime_noon) then begin ; before noon
    ; at 12:00 JST (3:00 UT)
    delta0=float(reform(text(1,whr(0)+[0,-1])))   ; in unit of (arc-degree)
    ha0=float(reform(text(2,whr(0)+[0,-1]))) ; in unit of (time-sec)
    pm=-1.
endif else begin  ; afternoon
    ; at 12:00 JST (3:00 UT)
    delta0=float(reform(text(1,whr(0)+[0,1])))   ; in unit of (arc-degree)   
    ha0=float(reform(text(2,whr(0)+[0,1]))) ; in unit of (time-sec)   
    pm=1.
endelse

delta(n)=delta0(0)+pm*(delta0(1)-delta0(0))*(utime_time-utime_noon)/86400.

ha_diff=ha0(0)+pm*(ha0(1)-ha0(0))*(utime_time-utime_noon)/86400.
hour_angle(n)=(ha_diff+(utime_time-utime_noon))*deg_per_sec

skip:
endfor

lonlat=norh_gt_lonlat()
lat=lonlat(1)

el=asin( $
      cos(delta/!radeg)*cos(hour_angle/!radeg)*cos(lat/!radeg) $
     +sin(delta/!radeg)*sin(lat/!radeg) $
    )*!radeg

zang=90-el

az=acos( $
         (sin(lat/!radeg)*cos(zang/!radeg)-sin(delta/!radeg)) $
         /(cos(lat/!radeg)*sin(zang/!radeg)) $
        )*!radeg

whr=where(hour_angle lt 0.,cnt)  ; before-noon
if (cnt ge 1) then az(whr)= -az(whr)

az=az+180.


return
end