function [solarhour,solarminute]=solar_time(year, doy, lat, lon, hour, minute, Timezone, DSTcorr, DSTdata)
%Usage
%[solarhour,solarmin]=solar_time(doy, lat, lon, hour,minute,Timezone,DSTcorr, DSTdata)
%Input
%year       [m]   year for m records
%doy        [m]   julian day number of year
%lat        [m]   degree latitude, 30 for Northern Hemi, -30 for southern
%lon        [m]   degree longitude, 116 for East, -116 for West
%hour       [m]   time in 24h format
%minute     [m]   minutes
%Timezone         time zone in GMT (+8 for Perth) 
%DSTcorr    [1/0] Correct for daylight saving time (DST) 
%DSTdata    []    Optional vector with periods of DST. Default includes times for Perth: 
%                 [year juliandayend, juliandaystart;
%                  year juliandayend, juliandaystart];etc
%Output
%solarhour         solar time, hour
%solarmin          solar time, minute
%Tom Schut, 2008

if nargin<7 DSTcorr=0; end
if nargin<8
    DSTdata=[2006 daynumber(2006,12,31) daynumber(2006,12, 3);...
             2007 daynumber(2007, 3,25) daynumber(2007,10,28);... 
             2008 daynumber(2008, 3,30) daynumber(2008,10,26);
             2009 daynumber(2009, 3,29) daynumber(2009,10,01)];
end


nd=length(doy);
if DSTcorr
    for i=1:nd
        a=find(year(i)==DSTdata(:,1));
        if doy(i)<DSTdata(a,2) | doy(i)>DSTdata(a,3) hour(i)=hour(i)-1; end
        if year>=2009 & doy(i) > daynumber(2009,10,01) & nargin<8
            disp('check DST dates in code of solar_time.m')
        end
    end
end
hour=hour + minute ./ 60;

%----------------------------------------------
%from website: http://squ1.org/wiki/Solar_Position_Calculator
t = (279.134 + 0.985647 * doy) * (pi()/180.0);
fEquation = ( 5.0323- 100.976 * sin(t)+ 595.275 * sin(2*t)+ 3.6858 * sin(3*t)...
		      - 12.47 * sin(4*t) - 430.847 * cos(t) + 12.5024 * cos(2*t) + 18.25 * cos(3*t)  ) ./ 3600;
fLongitude=lon*pi()/180;
fTimeZone=8 * 15 *pi()/180;
fDifference = (((fLongitude - fTimeZone) * 180/pi() ) * 4) / 60.0;

hour=hour + fEquation + fDifference; %solar hour

solarhour=floor(hour);
solarminute = floor((hour - solarhour) * 60.0);
%----------------------------------------------