function [Elsun,Azsun,Elmoon,Azmoon,phase,DECLsun,DECLmoon,rs,rm]=sunmoon(year,month,day,hour,LATpos,LONpos);

% Calculate sun and moon position  based on code on www.stjarnhimlen.se
% Inculuding Atmosphere refraction correction, only for output variable.
% Original code doesn't include this correction.
% Written by kmark Oct 2004.
% input:
%    year
%    month 
%    day - day NOT day of year!!
%    hour - UTC time
%    LATpos - Lat coordinate
%    LONpos - Lon ---//----
% output:
%      Elsun- elevation of sun
%      Azsun- azimuth of sun
%      Elmoon - elevation of moon
%      Azmoon - azimuth of moon
%      phase - Moon phase (0-1) 1 is full moon
%      DECLsun - sun declination without refraction correction
%      DECLmoon - moon ----------//------------------------
%      rs - SUN-Earth distance in AU
%      rm - EARTH-MOON distance in Earth Radii UNIT!!!!

% The time scale in these formulae are counted in days. Hours, minutes, 
% seconds are expressed as fractions of a day. Day 0.0 occurs at 2000 Jan 0.0 UT 
% (or 1999 Dec 31, 0:00 UT). This "day number" d is computed as follows 
% (y=year, m=month, D=date, UT=UT in hours+decimals):
% Note that ALL divisions here should be INTEGER divisions. In Pascal, use "div" 
% instead of "/", in MS-Basic, use "\" instead of "/". In Fortran, C and C++ "/" 
% can be used if both y and m are integers. Finally, include the time of the day, by adding:

y=year;
m=month;
D=day;
d = 367*y - floor(7 * ( y + floor((m+9)/12 )) / 4) + floor(275*m/9) + D - 730530;
d=d+hour/24.0;       % (this is a floating-point division)

% The primary orbital elements are here denoted as:
%    N = longitude of the ascending node
%    i = inclination to the ecliptic (plane of the Earth's orbit)
%    w = argument of perihelion
%    a = semi-major axis, or mean distance from Sun
%    e = eccentricity (0=circle, 0-1=ellipse, 1=parabola)
%    M = mean anomaly (0 at perihelion; increases uniformly with time)
% Related orbital elements are: 
%    w1 = N + w   = longitude of perihelion
%    L  = M + w1  = mean longitude
%    q  = a*(1-e) = perihelion distance
%    Q  = a*(1+e) = aphelion distance
%    P  = a ^ 1.5 = orbital period (years if a is in AU, astronomical units)
%    T  = Epoch_of_M - (M(deg)/360_deg) / P  = time of perihelion
%    v  = true anomaly (angle between position and perihelion)
%    E  = eccentric anomaly

% Orbital elements of the Moon:
Nm = 125.1228 - 0.0529538083 * d;
im = 5.1454;
wm = 318.0634 + 0.1643573223 * d;
am = 60.2666;                     % (Earth radii)
em = 0.054900;
Mm = 115.3654 + 13.0649929509 * d;
% conversion degree value to 0-360 range by rev functio
Nm=rev(Nm);
im=rev(im);
wm=rev(wm);
Mm=rev(Mm);

% Orbital elements of the Sun:
Ns = 0.0;
is = 0.0;
ws = 282.9404 + 4.70935E-5 * d;
as = 1.000000;                 % (AU)
es = 0.016709 - 1.151E-9 * d;
Ms = 356.0470 + 0.9856002585 * d;
% conversion degree value to 0-360 range by rev functio
Ns=rev(Ns);
is=rev(is);
ws=rev(ws);
Ms=rev(Ms);

oblecl=23.4393-3.563e-7*d;
oblecl=rev(oblecl);

% sun position calculations
% First, compute the eccentric anomaly E from the mean anomaly M and 
% from the eccentricity e (E and M in degrees):
E = Ms + es*(180/pi) * sin(deg2rad(Ms)) * ( 1.0 + es * cos(deg2rad(Ms)) );

% Note that the formulae for computing E are not exact; however they're 
% accurate enough here.
% Then compute the Sun's distance r and its true anomaly v from:
xv = cos(deg2rad(E)) - es;    % r * cos(v)
yv = sqrt(1.0 - es*es) * sin(deg2rad(E));    % r * sin(v) =
v = rad2deg(atan2( yv, xv ));
v=rev(v);
r = sqrt( xv*xv + yv*yv );
% Now, compute the Sun's true longitude:
lonsun = v + ws;
lonsun=rev(lonsun);
latsun=0;    % from definition;
% Convert lonsun,r to ecliptic rectangular geocentric coordinates xs,ys:
rs=r;
[DECLsun,RAsun]=eclip2DeclRa(latsun,lonsun,oblecl);
Ls=Ms+ws;
Ls=rev(Ls);
[Elsun,Azsun]=DeclRa2ElevAz(hour,Ls,LATpos,LONpos,DECLsun,RAsun);
Elsun=refrac(Elsun);
if nargout<=2
    % only sun position calulations
    return
end    

    
% Moon position calulation
% First, compute the eccentric anomaly, E, from M, the mean anomaly, and e, 
% the eccentricity. As a first approximation, do (E and M in degrees):
E = Mm + em*(180/pi) * sin(deg2rad(Mm)) * ( 1.0 + em * cos(deg2rad(Mm)) );
% If e, the eccentricity, is less than about 0.05-0.06, this approximation 
% is sufficiently accurate. If the eccentricity is larger, set E0=E 
% and then use this iteration formula (E and M in degrees):
% For each new iteration, replace E0 with E1. Iterate until E0 and E1 are 
% sufficiently close together (about 0.001 degrees). For comet orbits with 
% eccentricites close to one, a difference of less than 1E-4 or 1E-5 degrees 
% should be required.

if am>0.05
   E0=E;
  while 1  
    E1 = E0 - ( E0 - em*(180/pi) * sin(deg2rad(E0)) - Mm ) / ( 1 - em * cos(deg2rad(E0)) ); 
    if (abs(E0-E1)<1e-5)
      E=E1;
      break
    else
      E0=E1;
    end
  end 
end
E=rev(E);
% Now compute the planet's distance and true anomaly:
 x =  am * ( cos(deg2rad(E)) - em );       %=r * cos(v)
 y =  am * ( sqrt(1.0 - em*em) * sin(deg2rad(E)));  %  r * sin(v)
 r = sqrt( x*x + y*y);
 v = rad2deg(atan2( y, x ));
 v=rev(v);
 
% Compute the planet's position in 3-dimensional space:

xh = r * ( cos(deg2rad(Nm)) * cos(deg2rad(v+wm)) - sin(deg2rad(Nm)) * sin(deg2rad(v+wm)) * cos(deg2rad(im)) );
yh = r * ( sin(deg2rad(Nm)) * cos(deg2rad(v+wm)) + cos(deg2rad(Nm)) * sin(deg2rad(v+wm)) * cos(deg2rad(im)) );
zh = r * ( sin(deg2rad(v+wm)) * sin(deg2rad(im)) );

lonecl = rad2deg(atan2( yh, xh ));
lonecl=rev(lonecl);
latecl = rad2deg(atan2(zh, sqrt(xh*xh+yh*yh)));
r=sqrt(xh*xh+yh*yh+zh*zh); 

% Perturbations of the Moon
% If the position of the Moon is computed, and one wishes a better accuracy 
% than about 2 degrees, the most important perturbations has to be taken into 
% account. If one wishes 2 arc minute accuracy, all the following terms 
% should be accounted for. If less accuracy is needed, some of the smaller 
% terms can be omitted.

% First compute:

%    Ms, Mm             Mean Anomaly of the Sun and the Moon
%    Nm                 Longitude of the Moon's node
%    ws, wm             Argument of perihelion for the Sun and the Moon
%    Ls = Ms + ws       Mean Longitude of the Sun  (Ns=0)
%    Lm = Mm + wm + Nm  Mean longitude of the Moon
%    D = Lm - Ls        Mean elongation of the Moon
%    F = Lm - Nm        Argument of latitude for the Moon
%Orbital elements of the Sun:
%   N = 0.0
%    i = 0.0
%    w = 282.9404 + 4.70935E-5 * d
%    a = 1.000000  (AU)
%    e = 0.016709 - 1.151E-9 * d
%    M = 356.0470 + 0.9856002585 * d


Lm=Mm + wm + Nm;
Lm=rev(Lm);
F=Lm-Nm;
D=Lm-Ls;
% calculation corrections
dlon=-1.274 * sin(deg2rad(Mm - 2*D))+...          % (the Evection)
    +0.658 * sin(deg2rad(2*D))+...                % (the Variation)
    -0.186 * sin(deg2rad(Ms))+...                 % (the Yearly Equation)
    -0.059 * sin(deg2rad(2*Mm - 2*D))+...
    -0.057 * sin(deg2rad(Mm - 2*D + Ms))+...
    +0.053 * sin(deg2rad(Mm + 2*D))+...
    +0.046 * sin(deg2rad(2*D - Ms))+...
    +0.041 * sin(deg2rad(Mm - Ms))+...
    -0.035 * sin(deg2rad(D))+...                  %(the Parallactic Equation)
    -0.031 * sin(deg2rad(Mm + Ms))+...
    -0.015 * sin(deg2rad(2*F - 2*D))+...
    +0.011 * sin(deg2rad(Mm - 4*D));

% Add these terms to the Moon's latitude (degrees):


dlat=-0.173 * sin(deg2rad(F - 2*D))+...
    -0.055 * sin(deg2rad(Mm - F - 2*D))+...
    -0.046 * sin(deg2rad(Mm + F - 2*D))+...
    +0.033 * sin(deg2rad(F + 2*D))+...
    +0.017 * sin(deg2rad(2*Mm + F));

% Add these terms to the Moon's distance (Earth radii):

dr=-0.58 * cos(deg2rad(Mm - 2*D))+...
    -0.46 * cos(deg2rad(2*D));
long=lonecl+dlon;
lat=latecl+dlat;
dist=r+dr;
rm=dist;
[DECLmoon,RAmoon]=eclip2DeclRa(lat,long,oblecl);
[Elmoon,Azmoon]=DeclRa2ElevAz(hour,Ls,LATpos,LONpos,DECLmoon,RAmoon);
Elmoon=refrac(Elmoon);
% Moon phase calculation
% As usual we must use a different procedure for the Moon. Since the Moon 
% is so close to the Earth, the procedure above would introduce too big errors. 
% Instead we use the Moon's ecliptic longitude and latitude, mlon and mlat, 
% and the Sun's ecliptic longitude, mlon, to compute first the elongation, 
% then the phase angle, of the Moon:
slon=lonsun;
mlon=long;
mlat=lat;
elong=rad2deg(acos(cos(deg2rad(slon-mlon))*cos(deg2rad(mlat))));
FV=180-elong;
% When we know the phase angle, we can easily compute the phase:
phase=(1+cos(deg2rad(FV)))/2;
r=rs;
R=rm;
% Magnitude (brigthnes) of planet
mag=-21.62+5*log10(r*R)+0.026*FV+4.0E-9*FV^4;
% Light Intensity in Lux
L=10^((mag+14.05)/-2.5); % based on Tablice Astronomiczne

% Sun light intensity 127 kLux
% Full moon intensity 0.29 Lux


% ---------------------------------------------------
function [DECL,RA]=eclip2DeclRa(lat,lon,oblecl);
% converion ecliptic coordinates to Right Ascension and Declination
% (Equatiorial coordinate)

xeclip=cos(deg2rad(lon))*cos(deg2rad(lat));
yeclip=sin(deg2rad(lon))*cos(deg2rad(lat));
zeclip=sin(deg2rad(lat));

% rotate through oblecl=23.4 deg

xequat=xeclip;
yequat=yeclip*cos(deg2rad(oblecl))-zeclip*sin(deg2rad(oblecl));
zequat=yeclip*sin(deg2rad(oblecl))+zeclip*cos(deg2rad(oblecl));

RA=rad2deg(atan2(yequat,xequat));
RA=rev(RA);
DECL=rad2deg(atan2(zequat,sqrt(xequat^2+yequat^2)));

% ---------------------------------------------------
function [el,az]=DeclRa2ElevAz(hour,Ls,LATpos,LONpos,DECL,RA)
% conversion to Elevation and Evevation angle

UTC=hour;
GMSTO=(Ls+180)/15;
SIDTIME=GMSTO+UTC+LONpos/15;
HA=15*SIDTIME-RA;

x=cos(deg2rad(HA))*cos(deg2rad(DECL));
y=sin(deg2rad(HA))*cos(deg2rad(DECL));
z=sin(deg2rad(DECL));

xhor=x*cos(deg2rad(90-LATpos))-z*sin(deg2rad(90-LATpos));
yhor=y;
zhor=x*sin(deg2rad(90-LATpos))+z*cos(deg2rad(90-LATpos));

az=rad2deg(atan2(yhor,xhor))+180;
el=rad2deg(asin(zhor));
 % -----------------------------
 function res=rev(X)
 res=X-floor(X/360)*360;
 if res<0
     res=res+360;
 end    
 
 function res=refrac(el)
% calculate refraction correction for US stand. atm.
% el is in degree!!
 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
  res=el+refrac;




   







