[WD,Phi,k,stadium,mag]=physmoon(T)

function [WD,Phi,k,stadium,mag]=physmoon(T)

% function [WD,Phi,k,stadium,mag]=physmoon(T)
%
% Physical Ephemeris of the Moon
%
% Montenbruck/Pfleger, Acker/Jaschek
%
% 18.05.04 M.Penzkofer

% Constants
rho = pi/180;

% Moon
[lM,bM,rM] = moon(T);
% Sun
[lS,bS,rS] = sun2000(T);

% equatorial radius of the Moon (in km)
R_EQU = 1738.0;
% equatorial radius of the Earth (in km)
R_ERDE = 6378.14;
% 1 AU in km
AU     = 149597870.0;
% apparent equatorial radius of the Moon (in Degrees)
WD = 2.0*asin(R_EQU/(rM*R_ERDE)) / rho;
% phase angle/phase, stadium
delta = rM*R_ERDE;
R = rS*AU;
if (lM < lS)
	lM = lM+360;
end
El = lM-lS;
if (El < 180)
	stadium = +1;
elseif (El > 180)
	stadium = -1;
else
	stadium = 0;
end
rMS = sqrt(delta^2+R^2-2*delta*R*cos(El*rho));
Phi0 = acos((delta^2+rMS^2-R^2)/(2*rMS*delta)) / rho;
Phi = acos(cos(Phi0*rho)*cos(bM*rho)) / rho;
k = 0.5*(1+cos(Phi0*rho)*cos(bM*rho));
% magnitude
P = Phi/100.0;
MAG = 0.23 + ( 3.05 - ( 1.02 - 1.05*P ) * P ) * P;
mag = MAG + 5.0 * log( rS*delta/AU ) / log(10);