[x,y,z]=position(PNr,T)

function [x,y,z]=position(PNr,T)

% function [x,y,z]=position(PNr,T)
%
% position of the planets based on Kepler orbits
% mean orbital elements of epoche J2000 for Mercury to Mars,
% osculating orbital elements for epoche 10.10.1995 (JD 2450000.5) for Jupiter to Pluto
% relative accuracy ca. 0.001 for years 1990 to 2000
% algorithms by Montenbruck/Pfleger: "Astronomy with the Personal Computer"
%
% 05.05.04 M.Penzkofer

% check input values
if (PNr == 0)
   disp('Moon?')
   x = NaN;
   y = NaN;
   z = NaN;
   return
end
if (PNr < 1 | PNr > 9)
   disp('Enter number between 1 and 9!');
   x = NaN;
   y = NaN;
   z = NaN;
   return
end

% heliocentric ecliptical orbital elements in respect of equinox J2000

switch PNr

	case 1 % Mercury
                  A= 0.387099; E=0.205634; M=174.7947; N=149472.6738;
                  O= 48.331;   I= 7.0048;  W= 77.4552;

	case 2 % Venus
                  A= 0.723332; E=0.006773; M= 50.4071; N= 58517.8149;
                  O= 76.680;   I= 3.3946;  W=131.5718;

	case 3 % Earth
                  A= 1.000000; E=0.016709; M=357.5256; N= 35999.3720;
                  O=174.876;   I= 0.0000;  W=102.9400;

	case 4 % Mars
                  A= 1.523692; E=0.093405; M= 19.3879; N= 19140.3023;
                  O=49.557;    I= 1.8496;  W=336.0590;

	case 5 % Jove
                  A= 5.202437; E=0.048402; M=250.3274; N=  3035.2275;
                  O=100.4683;  I=1.3047;   W= 15.7192;

	case 6 % Saturn
                  A= 9.551712; E=0.052340; M=267.2465; N=  1219.6465;
                  O=113.6439;  I=2.4855;   W= 90.9682;

	case 7 % Uranus
                  A=19.293108; E=0.044846; M=118.4320; N=   424.8150;
                  O= 74.0903;  I=0.7733;   W=176.6152;

	case 8 % Neptune
                  A=30.257162; E=0.007985; M=292.4716; N=   216.3047;
                  O=131.7750;  I=1.7700;   W=  3.0962;

	case 9 % Pluto
                  A=39.783607; E=0.254351; M=  8.2304; N=   143.4629;
                  O=110.3865;  I=17.1201;  W=224.7424;

end

if (PNr >= 1 & PNr <= 4)
   T0= 0.0;
elseif (PNr >=5 & PNr <= 9)
   T0=-0.042286105;
end

M = M + N*(T-T0);

% cartesian coordinates in respect of mean equinox J2000
PQR = gaussvec(O,I,W-O);
[XX,YY,VVX,VVY] = ellip(M,A,E);
[x,y,z] = orbecl(XX,YY,PQR);