[xP,yP,zP,xS,yS,zS,x,y,z,delta0]=geocen(T,lP,bP,rP,lS,bS,rS,PNr,Modus)

function [xP,yP,zP,xS,yS,zS,x,y,z,delta0]=geocen(T,lP,bP,rP,lS,bS,rS,PNr,Modus)

% function [xP,yP,zP,xS,yS,zS,x,y,z,delta0]=geocen(T,lP,bP,rP,lS,bS,rS,PNr,Modus)
%
% geocentric coordinates from heliocentric coordinates
% algorithms by Montenbruck/Pfleger: "Astronomy with the Personal Computer"
%
% 05.05.04 M.Penzkofer

% constants
P2 = 2*pi;

% Init
dl = 0;
db = 0;
dr = 0;
dlS = 0;
dbS = 0;
drS = 0;

if (Modus > 0)
   
   	% Sun
      M = P2*frac(0.9931266+ 99.9973604*T);
      dlS = 172.00+5.75*sin(M);
      drS = 2.87*cos(M);
      dbS = 0.0;

      % dl,db in 1e-4 rad/d, dr in 1e-4 AE/d
      switch PNr
      	case 0
         	% Sun
            dl=0.0;
            dr=0.0;
            db=0.0;
         case 1
            % Mercury
            M = P2*frac(0.4855407+415.2014314*T);
            dl = 714.00+292.66*cos(M)+71.96*cos(2*M)+18.16*cos(3*M)+4.61*cos(4*M)+3.81*sin(2*M)+2.43*sin(3*M)+1.08*sin(4*M);
            dr = 55.94*sin(M)+11.36*sin(2*M)+2.60*sin(3*M);
            db = 73.40*cos(M)+29.82*cos(2*M)+10.22*cos(3*M)+3.28*cos(4*M)-40.44*sin(M)-16.55*sin(2*M)-5.56*sin(3*M)-1.72*sin(4*M);
         case 2
            % Venus
            M = P2*frac(0.1400197+162.5494552*T);
            dl = 280.00+3.79*cos(M);
            dr = 1.37*sin(M);
            db = 9.54*cos(M)-13.57*sin(M);
         case 3
            % Earth
            dl=dlS;
            dr=drS;
            db=-dbS;
         case 4
            % Mars
            M = P2*frac(0.0538553+53.1662736*T);
            dl = 91.50+17.07*cos(M)+2.03*cos(2*M);
            dr = 12.98*sin(M)+1.21*cos(2*M);
            db =  0.83*cos(M)+2.80*sin(M);
         case 5
            % Jove
            M = P2*frac(0.0565314+8.4302963*T);
            dl = 14.50+1.41*cos(M);
            dr = 3.66*sin(M);
            db = 0.33*sin(M);
         case 6
            % Saturn
            M = P2*frac(0.8829867+3.3947688*T);
            dl = 5.84+0.65*cos(M);
            dr = 3.09*sin(M);
            db = 0.24*cos(M);
         case 7
            % Uranus
            M = P2*frac(0.3967117+1.1902849*T);
            dl = 2.05+0.19*cos(M);
            dr = 1.86*sin(M);
            db =-0.03*sin(M);
         case 8
            % Neptune
            M = P2*frac(0.7214906+0.6068526*T);
            dl = 1.04+0.02*cos(M);
            dr = 0.27*sin(M);
            db = 0.03*sin(M);
         case 9
            % Pluto
            M = P2*frac(0.0385795+0.4026667*T); 
            dl = 0.69+0.34*cos(M)+0.12*cos(2*M)+0.05*cos(3*M);
            dr = 6.66*sin(M)+1.64*sin(2*M);
            db = -0.08*cos(M)-0.17*sin(M)-0.09*sin(2*M);
      end

end

[xS,yS,zS,vxS,vyS,vzS]=posvel(lS,bS,rS,dlS,dbS,drS);
[xP,yP,zP,vx, vy, vz ]=posvel(lP,bP,rP,dl ,db ,dr );
x = xP+xS;
y = yP+yS;
z = zP+zS;
delta0 = sqrt(x*x+y*y+z*z);

if (PNr == 3)
   x = 0;
   y = 0;
   z = 0;
   delta0 = 0;
end

FAC = 0.00578*delta0*0.0001;
if (Modus == 1)
   x = x-FAC*vx;
   y = y-FAC*vy;
   z = z-FAC*vz;
elseif (Modus == 2)
   x = x-FAC*(vx+vxS);
   y = y-FAC*(vy+vyS);
   z = z-FAC*(vz+vzS);
end