; $Id: solarpos.pro
;+
; NAME:
;   SOLARPOS
;
; PURPOSE:
;   computes solar vector, azimuth and elevation for a given day, time,
;   longitude and latitude
;
; CATEGORY:
;   Astronomy
;   Solar energy
;
; CALLING SEQUENCE:
;   RESULT = SolarPos(sundata,[LIMB=limb])
;
; INPUTS:
;    data structure SUNDATA:
;     sundata.year
;     sundata.JDay
;     sundata.time
;     sundata.latitude
;     sundata.longitude
;     sundata.timezone
; KEYWORDS
;   LIMB: set this keyword to add 50 arcminutes to compute the occultation
;     and rise of the sun's upper limb: 16' sun radius + 34' refraction
; OUTPUTS:
;
;   data structure SOLPOS:
;     solpos.vector
;     solpos.azimuth
;     solpos.elevation
;     solpos.declination
;     solpos.sunrise
;     solpos.sunset
;     solpos.EqTime
;     solpos.omega
;
;
; REFERENCES
;; Bourges, B.: 1985, Improvement in solar declination computation,  Solar
;  Energy 35(4),367--369.
;
;; Corripio, J.G. 2002: "Vectorial algebra algorithms for calculating terrain
;  parameters from DEMs and the position of the sun for solar radiation
;  modelling in mountainous terrain" International Journal of Geographical
;  Information Science. In press.
;
;; Spencer, J.W.: 1971, Fourier series representation of the position of the sun,
;  2,172.
;
; MODIFICATION HISTORY:
;    Javier Corripio,  jgc@geo.ed.ac.uk
;     Created   August, 2001
;     modified:    September 2002
;
;-
;  Copyright (c) 2002 Javier G Corripio.
;  corripio@members.ises.org corripio@ifu.baug.ethz.ch
;
;  This software may be used, copied, or redistributed as long as it is not
;  sold and this copyright notice is reproduced on each copy made.  This
;  routine is provided as is without any express or implied warranties
;  whatsoever.


FUNCTION solarpos, SUNDATA, LIMB=limb, TWILIGHT=twilight

On_error,2
nparams=N_PARAMS()
;** check that the argument exists
message=[ ['USAGE: Result = SolarPos(sundata [,/Limb] [,TWILIGHT=Twilight])'],$
   ['sundata={year:year,JDAY:jday,TIME:time,LATITUDE:latitude,LONGITUDE:longitude,TIMEZONE:timezone}'] ]
IF (nparams EQ 0) THEN BEGIN
    print,message
        RETURN,!values.F_NaN
ENDIF

;** check that the argument has the right dimensions
sz=size(sundata,/structure)
szdim=size(sz.dimensions,/dimensions)
szdim=fix(szdim[0])
IF (szdim LT 5) THEN BEGIN
    print,message
    print,'structure data has wrong dimensions:'+string(szdim)
        RETURN,!values.F_NaN
ENDIF

DRADEG = 180D/!DPI
year=sundata.year
JDay=sundata.JDay+(sundata.time/24.)

; *****************computes declination

 w=(360.0/365.25)
 J0=78.801+0.2422*(Year-1969)-FIX(0.25*(Year-1969))
 t=Jday-0.5-J0
 dayang=w*t

dayang = dayang/DRADEG

;*** Declination: Bourgues, 1985

 delta = 0.3723 $
       + 23.2567*SIN(dayang)  - 0.758*COS(dayang) $
       + 0.1149*SIN(2*dayang) + 0.3656*COS(2*dayang) $
       - 0.1712*SIN(3*dayang) + 0.0201*COS(3*dayang)

delta = delta/DRADEG

;** Computes equation of time according to Spencer(1971).

daynumberS = 2*!DPI*(JDay-1)/365.
EqTime = 0.000075 + $
       0.001868*cos(daynumberS) - 0.032077*sin(daynumberS) - $
       0.014615*cos(2*daynumberS) - 0.040849*sin(2*daynumberS)

EqTime = EqTime * 12/!PI    ;** EqTime in hours
;** hour_angle (omega):
longitude = sundata.longitude
stndmeridian = sundata.timezone*15.
deltaLongTime=longitude-stndmeridian ;** angle to time zone
deltaLongTime = deltaLongTime * 24/360. ;** angle to hours

omega = !PI*( ( (sundata.time+deltaLongTime+EqTime)/12 ) - 1)

lat=sundata.latitude / DRADEG

SunVector = findgen(3)
SunVector(0)= -sin(omega)*cos(delta)
SunVector(1)= sin(lat)*cos(omega)*cos(delta)-cos(lat)*sin(delta)
SunVector(2)= cos(lat)*cos(omega)*cos(delta)+sin(lat)*sin(delta)

;** Solar elevation:
elevation = ASIN(SunVector(2)) * DRADEG

;** Solar azimuth

IF (SunVector(0) EQ 0) AND (SunVector(1) EQ 0) THEN BEGIN
    azimuth = 0.0
    ENDIF ELSE BEGIN
      azimuth = ( !PI - ATAN(SunVector(0),SunVector(1)) )  * DRADEG
    ENDELSE

;** Sunrise and sunset
omeganul = ACOS(-TAN(lat)*TAN(delta)) ; **Local Solar Time

;** next add/substract 50 arcminutes to compute the occultation
    ;  and rise of the sun's limb 16' diameter + 34' refraction
IF keyword_set(limb) THEN omeganul=omeganul+((.5*100.0/60.0)/DRADEG)
sunrise = 12*(1-omeganul/!DPI)-deltaLongTime-EqTime ; **back to official time
sunset = 12*(1+omeganul/!DPI)-deltaLongTime-EqTime

IF ARG_PRESENT(twilight) THEN BEGIN
sixbelowh = cos(96./dradeg)
omega_6 = ACOS((sixbelowh - sin(lat)*sin(delta))/(cos(lat)*cos(delta)))
srtwt = 12*(1-omega_6/!DPI)-deltaLongTime-EqTime ; **back to official time
sstwt = 12*(1+omega_6/!DPI)-deltaLongTime-EqTime
twilight = [srtwt,sstwt]
ENDIF


declination = delta * DRADEG

solpos=CREATE_STRUCT('vector',sunVector,'azimuth',azimuth,'elevation',elevation,$
'declination',declination,'sunrise',sunrise,'sunset',sunset,'EqTime',EqTime,$
    'omega',omega)


RETURN, solpos

END