function [az,el]=sunpos(year,day,hour,lat,long)
% calculate azymuth and sun elevation angle
% input:
%      year
%      day - julian day (DAY of Year)
%      hour - decimal time
%      lat - latitude in deg
%      lon - longitude in deg
% output:
%        az - sun azimuth angle
%        el - sun elevation angle  
 twopi=2.*pi;
 rad=pi/180.;

 delta=year-1949.;    % get the current Julian date
 leap=fix(delta/4.);
 jd=32916.5+delta*365.+leap+day+hour/24.;

 time=jd-51545.0;     % calculate ecliptic coordinates

 mnlong=280.460+0.9856474*time;      % force mean longitude between 0 and 360 degs
 mnlong=mod(mnlong,360.);
 ind=find(mnlong<0);
 if isempty(ind)==0
  mnlong(ind)=mnlong(ind)+360.;
 end

  mnanom=357.528+0.9856003*time;   % mean anomaly in radians between 0, 2*pi
  mnanom=mod(mnanom,360.);
  ind=find(mnanom<0);
  if isempty(ind)==0
   mnanom(ind)=mnanom(ind)+360.;
  end

  mnanom=mnanom*rad;
 % compute ecliptic longitude and obliquity of ecliptic
  eclong=mnlong+1.915*sin(mnanom)+0.020*sin(2.*mnanom);
  eclong=mod(eclong, 360.);
  ind=find(eclong<0);
  if isempty(ind)==0
    eclong(ind)=eclong(ind)+360.;
  end

  oblqec=23.429-0.0000004*time;
  eclong=eclong*rad;
  oblqec=oblqec*rad;

  num=cos(oblqec).*sin(eclong);   % calculate right ascention and declination
  den=cos(eclong);
  ra=atan(num./den);

  ind=find(den<0);
  if isempty(ind)==0       % force ra between 0 and 2*pi
    ra(ind)=ra(ind)+pi;
  end    
  ind=find(num<0);  
  if isempty(ind)==0
    ra(ind)=ra(ind)+twopi;
  end

dec=asin(sin(oblqec).*sin(eclong));     % dec in radians

gmst=6.697375+0.0657098242*time+hour;  % calculate Greenwich mean sidereal time in hours
% hour not changed to sidereal sine "time" includes the fractional day
  gmst=mod(gmst,24.);
  ind=find(gmst<0);
  if isempty(ind)==0
    gmst(ind)=gmst(ind)+24.;
  end

  lmst=gmst+long/15.;    % calculate local mean sidereal time in radians
  lmst=mod(lmst,24.);
  ind=find(lmst<0);
  if isempty(ind)==0
    lmst(ind)=lmst(ind)+24.;
  end

  lmst=lmst*15.*rad;
% calculate hour angle in radians between -pi, pi
  ha = lmst -ra;
  ind=find(ha<-pi);
  if isempty(ind)==0
     ha(ind)=ha(ind)+twopi;
  end
  ind=find(ha>pi);
  if isempty(ind)==0
     ha(ind)=ha(ind)-twopi;
  end
  lat=lat*rad;

% calculate azimuth and elevation
  el=asin(sin(dec).*sin(lat)+cos(dec).*cos(lat).*cos(ha));
  az=asin(-cos(dec).*sin(ha)./cos(el));

% dd J. W. Spencer code (next 5 lines)
 term=sin(dec) - sin(el).*sin(lat);
 ind=find(term>=0);
 if isempty(ind)==0
   index=find(sin(az(ind))<0);
   if isempty(index)==0
     az(ind(index))=az(ind(index))+twopi;
   end
 end
 ind=find(term<0);
 if isempty(ind)==0
   az(ind)=pi-az(ind);
 end
% end Spencer's corrections

%%%        elc=asin(sin(dec)/sin(lat))
%%%        if(el.ge.elc) az=pi-az
%%%        if(el.le.elc  .and. ha.gt.0.) az=twopi+az

% this puts azimuth between 0 and 2*pi radians
% calculate refraction correction for US stand. atm.
  el=el/rad;
   ind=find(el>=-0.56);
   refrac=zeros(size(el));
   if isempty(ind)==0
     refrac(ind)=3.51561*(0.1594+0.0196*el(ind)+0.00002*(el(ind)).^2)./(1.+0.505*el(ind)+0.0845*(el(ind)).^2);
   end     
   ind=find(el<-0.56);
   if isempty(ind)==0
     refrac(ind)=0.56;
   end
   el=el+refrac;
   az=az/rad;
   lat=lat/rad;
   return;

