[R,D,ELONG,PHI,K]=illum(x,y,z,xE,yE,zE)

function [R,D,ELONG,PHI,K]=illum(x,y,z,xE,yE,zE)

% function [R,D,ELONG,PHI,K]=illum(x,y,z,xE,yE,zE)
%
% illumination parameters of a planet
% algorithms by Montenbruck/Pfleger: "Astronomy with the Personal Computer"
%
% 05.05.04 M.Penzkofer

% constants
rho = pi/180;

% geocentric position of the planet
xP=x-xE;
yP=y-yE;
zP=z-zE;

% distances in the triangle Sun-Earth-Planet
R  = sqrt(  x*x  +  y*y  +  z*z  );   % distance Sun-Planet
RE = sqrt( xE*xE + yE*yE + zE*zE );   % distance Sun-Earth
D  = sqrt( xP*xP + yP*yP + zP*zP );   % distance Earth-Planet

% elongation, phase angle and phase
ELONG = acos( ( D*D + RE*RE - R*R ) / ( 2.0*D*RE ) ) / rho;

C_PHI = ( D*D + R*R - RE*RE ) / ( 2.0*D*R ) ;
PHI   = acos( C_PHI ) / rho;
K     = 0.5*(1.0+C_PHI);