; $Id: solarposR.pro,v 6.0 28/07/2004
;+
; NAME:
;   SOLARPOSR
;
; PURPOSE:
;   computes solar vector, azimuth and elevation for a given day, time,
;   longitude and latitude, implements refraction correction
;
; CATEGORY:
;   Astronomy
;   Solar energy
;
; CALLING SEQUENCE:
;   RESULT = solarposr(sundata,[LIMB=limb])
;
; INPUTS:
;    data structure SUNDATA:
;     sundata.year
;     sundata.JDay
;     sundata.time
;     sundata.latitude
;     sundata.longitude
;     sundata.timezone
;     sundata.altitude
;     sundata.temperature
; 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.
;
;; Solar position algorithm for solar radiation applications
;   Solar Energy, Volume 76, Issue 5, 2004, Pages 577-589
;   Ibrahim Reda and Afshin Andreas
;
;; 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.
;  jgc@geo.ed.ac.uk  http://www.geo.ed.ac.uk/~jgc/
;  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 solarposr, SUNDATA, LIMB=limb

On_error,2
nparams=N_PARAMS()
;** check that the argument exists
message=[ ['USAGE: Result = solarposR(sundata)'],$
   ['sundata={year:year,JDay:JD,time:time,timezone:tmz,latitude:latitude,longitude:longitude, altitude:altitude,temperature:temperature}'] ]
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 6) THEN BEGIN
    print,message
    print,'structure data imcomplete values:'+string(szdim)
        RETURN,!values.F_NaN
ENDIF

DRADEG = 180./!DPI
year=sundata.year
JDay=sundata.JDay+(sundata.time/24.)
sundata.temperature=sundata.temperature+273.15
; *****************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):
long = sundata.longitude
stndmeridian = sundata.timezone*15.
deltaLatTime=long-stndmeridian ;** angle to time zone
deltaLatTime = deltaLatTime * 24/360. ;** angle to hours

omega = !PI*( ( (sundata.time+deltaLatTime+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

;Refraction Correction According to Reda & Andreas, 2004
Psite = z2p(sundata.altitude)
deltaelev = (Psite/1010.)*(283./sundata.temperature)* $
         (1.02/( 60.*TAN((elevation+10.3/(elevation+5.11))/DRADEG)) )
elevation = elevation + deltaelev

SunVector(2) = SIN(elevation/DRADEG)
Sunvector = unitvector(SunVector)

;** 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+0.014544410433 ;(50 arcminutes)
sunrise = 12*(1-omeganul/!DPI)-deltaLatTime-EqTime ; **back to official time
sunset = 12*(1+omeganul/!DPI)-deltaLatTime-EqTime

declination = delta * DRADEG

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


RETURN, solpos

END