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.);
 if (mnlong < 0)
  mnlong=mnlong+360.;
 end

  mnanom=357.528+0.9856003*time;   % mean anomaly in radians between 0, 2*pi
  mnanom=mod(mnanom,360.);
  if(mnanom < 0.)
   mnanom=mnanom+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.);
  if(eclong < 0)
    eclong=eclong+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);

  if(den < 0.)       % force ra between 0 and 2*pi
    ra=ra+pi;
  elseif(num < 0.)
    ra=ra+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.);
  if(gmst < 0.)
    gmst=gmst+24.;
  end

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

  lmst=lmst*15.*rad;
% calculate hour angle in radians between -pi, pi
  ha = lmst -ra;
  if(ha < -pi)
     ha=ha+twopi;
  end
  if(ha > pi)
     ha=ha-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)
 if(sin(dec) - sin(el)*sin(lat) >= 0. )
   if(sin(az) <  0.)
     az=az+twopi;
   end
  else
   az=pi-az;
 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;
   if(el >= -0.56)
     refrac=3.51561*(0.1594+0.0196*el+0.00002*el^2)/(1.+0.505*el+0.0845*el^2);
   else
     refrac=0.56;
   end
   el=el+refrac;
   az=az/rad;
   lat=lat/rad;
   return;

