; ----------------------------------------------------------------------
; Check header quantities such as the parallactic angle.
;
; Compute Position Angle (paralactic angle)..
; taken from SKY() subroutine in astro_lib.f in makee package.
;
;  HA : Hour Angle  (0h at meridian, negative in east)
;  AZI: Azimuth  (0 degrees due north, 90 deg. due east, ...)
;  ZT : Zenith angle (True)  (0 degrees at zenith)
;  PA : Position Angle  (0 deg. red end toward due south, 90 deg due west,...)
;(In the triangle zenith—object—celestial pole, the parallactic angle will be 
; the position angle of the zenith at the celestial object.  So the PA points
; the blue end at the zenith. On a N(up) and E(right) plot, PA=0 degrees would
; be pointing up (or north) and PA=90 would be pointing right (or east).)
;
; NOTE: NIRES is mounted on the Keck 2 telescope.
;


pro nsx_header_check, IMG
  common nsx 

; Echo and set.
    print, " _ _ _  _ _ _  _ _ _  Header check _ _ _  _ _ _  _ _ _ "
    Dec   = IMG.dec
    HAdeg = IMG.ha

    if (Dec lt -90.) then begin
       print, "No RA,Dec info.."
       return
    endif

; Calculate AZImuth and Zenith angle (True)
    deg = !dpi/180.d
    zp = sin(deg*(Dec))                                                
    xp = cos(deg*(Dec)) * sin(deg*(HAdeg))
    yp = cos(deg*(Dec)) * cos(deg*(HAdeg))
    x  = xp
    z  = ( zp * sin(deg*(KECK2LAT)) ) + ( yp * cos(deg*(KECK2LAT)) )
    y  = ( yp * sin(deg*(KECK2LAT)) ) - ( zp * cos(deg*(KECK2LAT)) )
    if (z ne 0.) then zt = ( atan( sqrt((x*x)+(y*y)) / z ) )/deg else zt = 90.0
    if (y ne 0.) then begin
       azi = ( atan(x/y) )/deg
    endif else begin
       if (x gt 0.) then azi = 270. else azi = 90.
    endelse
  
    if (y gt 0.) then azi = 180. + azi
    if (azi lt 0.) then azi = 360. + azi

;
;  Now calculate position angle... (0 deg.=red end toward south, 90deg.=west)
;  Need to use x,y,z alt-azimuth euclidean coordinates. and the spherical
;  triangle rule: cos(a)=cos(b)cos(c)+sin(b)sin(c)cos(ANGLE) where ANGLE is
;  opposite "side" a.    ..a,b,c are in radians..
;
    a = deg*(90. - KECK2LAT)
    b = acos( (y*(-1.) * cos(deg*(KECK2LAT))) + (z * sin(deg*(KECK2LAT))) )
    c = acos(z)
    if ( (sin(b) eq 0.)||(sin(c) eq 0.) ) then begin
       pa = 999.0
    endif else begin
       r = ( cos(a) - (cos(b) * cos(c)) ) / ( sin(b) * sin(c) )
; adjust for possible numerical(machine) problems */
       if (ABS((r)) gt 1.0) then begin
          pa = 999.0
          if ((r ge  1.0)&&(r lt  1.001)) then pa=0.
          if ((r le -1.0)&&(r gt -1.001)) then pa=180.0
       endif else begin
          pa = (( acos(r) ))/deg
       endelse
       if (HAdeg lt 0.) then pa = -1. * pa
    endelse
    za = gal_apparent_zenith_angle( zt )

    if (za lt 0.) then za=zt                                                               
    air = gal_airmass( za )     
    print, IMG.ra,IMG.ha,IMG.dec, format='("  RA=",F9.5,"  HA=",F9.5,"  Dec=",F9.5)'
    print, format='("  zt   = ",F9.5,"  [zenith]=",F9.5,"  dif=",F9.5)',zt,90.-IMG.el,zt-(90.-IMG.el)
    print,format='("  za   = ",F9.5,"  [zenith]=",F9.5,"  dif=",F9.5)',za,90.-IMG.el,za-(90.-IMG.el)
    print,format='("  azi  = ",F9.5,"  .az     =",F9.5,"  dif=",F9.5)',azi,IMG.az,azi-IMG.az
    print,format='("  pa   = ",F9.5,"  .parang =",F9.5,"  dif=",F9.5)',pa,IMG.parang,pa-IMG.parang

; 
;   Now calculate differential refraction relative to 12000 angstroms in arcsecs.
;   (pressure=760 mm Hg, T=15 deg Celsius, dry air) from CRC and other sources
;
    wv = [5000.,8000.,9000.,10000.,11000.,12000.,13000.,14000.,15000.,16000.,17000.,18000.,19000.,20000.,21000.,22000.,23000.,24000.]
    for ii=0, 17 do begin
       dar = nsx_difatmref( wv[ii], zt ) 
       print, format='(" ",F8.1," ang. rel. to. ",F8.1," ang (ref) : dar=",F8.4," (",F8.4," px)")',wv[ii],DAR_refwave,dar,dar/ARCSEC_PER_PIXEL
    endfor

    print, " _ _ _  _ _ _  _ _ _  _ _ _ _ _ _  _ _ _  _ _ _  _ _ _ "


end