PRO sunangle,D, time, lat, lon, theta, azimuth


;***************************************************************;
;                               ;
;   PROCEDURE sunangle                  ;
;                               ;
;   Written by Angela Strub, October 1998           ;
;   IDL Version of radtran.for by W.W.Gregg         ;
;   Modified by Chuanmin Hu, October 1999
;                               ;
;---------------------------------------------------------------;
;                               ;
; Computes sun azimuth and zenith angles for a given        ;
; time, date, latitude and longitude.  This program     ;
; is from the NMFS ELAS computer code.  Modified for        ;
; standard coordinates (W long. negative), to correct       ;
; for dateline problem, and to correct coefficients (taken  ;
; from Iqbal, 1983, An Introduction to Solar Radiation).    ;
; Watson Gregg, Research and Data Systems, Corp.        ;
;                               ;
;---------------------------------------------------------------;
;                               ;
;  VARIABLES:                           ;
;                               
;   D   = julian day  1-365 as input 
;   time    = time of day in decimal GMT hour, e.g. 17.5129, as input
;   lat     = latitude in decimal degrees as input  ;
;   lon     = longitude in decimal degrees as input, negative for west  ;
;   theta   = solar zenith in degrees as output
;   azimuth = solar azimuth in degrees as output
;---------------------------------------------------------------
;  INTERNAL VARIABLES:
;   declin  = solar declination angle in radians        ;
;   theta_0 = theta zero orbital position in radians    ;
;   correct = time correction               ;
;   angle   = solar hour angle in degrees           ;
;                               ;
;***************************************************************;

;***************************************************************;
;                               ;
;   Compute solar declination angle             ;
;                               ;
;***************************************************************;

rad=double(180./!pi)    ;radiance to degree conversion factor

theta_0 = double(360.0 * double(D - 1) / 365.0) / rad

declin  = 0.396372 - 22.91327*COS(theta_0)          $
      + 4.02543*SIN(theta_0)                $
      - 0.387205*COS(2.0*theta_0)               $
      + 0.051967*SIN(2.0*theta_0)               $
      - 0.154527*COS(3.0*theta_0)               $
      + 0.084798*SIN(3.0*theta_0)

declin = declin / rad   ;convert to radians



;***************************************************************;
;                               ;
;   Time correction for solar hour angle            ;
;                               ;
;***************************************************************;

correct  = 0.004297 + 0.107029*COS(theta_0)             $
       - 1.837877*SIN(theta_0) - 0.837378*COS(2.0*theta_0)  $
       - 2.342824*SIN(2.0*theta_0)

angle = (time - 12.0) * 15.0 + lon + correct
IF angle GT  180.0 THEN angle = angle - 360.0
IF angle LT -180.0 THEN angle = angle + 360.0

lat_rad = lat / rad
lon_rad = lon / rad

angle_rad  = angle / rad



;***************************************************************;
;                               ;
;   Sun Zenith                      ;
;                               ;
;***************************************************************;

tmp1 = SIN(lat_rad)*SIN(declin) + COS(lat_rad)*COS(declin)*COS(angle_rad)
tmp2   = ABS(tmp1)

IF tmp2 GT 1.1 THEN BEGIN
    PRINT, 'Error in acos argument in sun zenith'
    PRINT, tmp1
ENDIF ELSE BEGIN                            $
    IF tmp2 GT 1.0 THEN BEGIN
    IF tmp1 GT 0.0 THEN tmp1=1.0
    IF tmp1 LT 0.0 THEN tmp1=-1.0
    ENDIF
ENDELSE

theta_rad = ACOS(tmp1)



;***************************************************************;
;                               ;
;   Sun Azimuth                     ;
;                               ;
;***************************************************************;

tmp1 = SIN(ABS(angle_rad)) *COS(declin) / SIN(theta_rad)
azimuth_rad  = ASIN(tmp1)

IF lat GT declin THEN azimuth_rad = 180.0 / rad - azimuth_rad
IF angle  GT 0.0  THEN azimuth_rad = 360.0 / rad - azimuth_rad



;***************************************************************;
;                               ;
;   Convert to degrees                  ;
;                               ;
;***************************************************************;

azimuth = azimuth_rad * rad
theta = theta_rad * rad


RETURN

END