function [r,v,ele]=ephem(CB,Bi,dateenc)
%[R,V,elements]=ephem(CentralBody,Body,ML Date#)
%R(km) & V(km/s) Heliocentric ecliptic of 2000
%R & V column vectors
%elements[a e i W w	q*]

%Data from JPL website (ssd.jpl.nasa.gov) 
if diff(size(Bi))>0;CB=CB';Bi=Bi';dateenc=dateenc';end
%(gravitational param of sun and planets (possible central bodies)
mu=[1.3271244e11 2.20320808e4 3.2485877e5 3.98600485043e5 4.282828685e4 1.267120008e8 3.7934000093e7 5.803200946e6 6.8713006597e6 9.47209e2];
AU=1.49597870691e8;
nt=length(dateenc);

switch CB
    
case 0
	   epoch=datenum(2000,01,01,12,0,0);
       time=dateenc-epoch;%days
       %ele( 	a(au) 		e 			i(deg) 	W(deg) 		w(deg) 		M(deg)
       ele=[    .38709893	.20563069	7.00487	 48.33167	 29.12478	174.79439
                .72333199	.00677323	3.39471	 76.68069	 54.85229	 50.44675 
               1.00000011	.01671022	0.00005	-11.26064	114.20783	357.51716 
               1.52366231	.09341233	1.85031	 49.57854	286.46230	 19.41248 
               5.20336301	.04839266	1.30530	100.55615	274.14035	 19.65053
               9.53707032	.05415060	2.48446	113.71504	338.71690	317.51238 
              19.19126393	.04716771	0.76986	 74.22988	 96.73436	142.26794 
              30.06896348	.00858587	1.76917	131.72169	273.24966	259.90868 
              39.48168677	.24880766	17.14175	110.30347	113.76329	 14.86205];
       ele=(ones(size(ele,1),1)*[AU 1 pi/180 pi/180 pi/180 pi/180]).*ele;%convert to km & rad
       
       %/Julian Century, angles in arcseconds
       eledot=[ .00000066	0.00002527	-23.51	 -446.30		1019.87	538101054.72
                .00000092	-.00004938	 -2.86	 -996.89		 888.09	210664244.86
              	-.00000005	-.00003804	-46.94	-18228.25      19426.53	129596542.35
              	-.00007221	0.00011902	-25.47	-1020.19		2580.97	 68903543.00
                .00060737	-.00012880	 -4.15	 1217.17		-377.24	 10924238.42
              	-.00301530	-.00036762	  6.11	-1591.05		-357.84	  4403001.84
              	 .00152025	-.00019150	 -2.09	-1681.40		2996.96	  1541235.23
              	-.00125196	0.00002510	 -3.64	 -151.25	   -693.18		787293.64
               -.00076912	0.00006465	 11.07	  -37.33		-94.92	   522880.15];
		 eledot=(ones(size(eledot,1),1)*[AU 1 pi/180/3600 pi/180/3600 pi/180/3600 pi/180/3600]).*eledot/36525;%convert km, rad, day
case 3
        epoch=datenum(1976,08,10,0,0,0);
        time=dateenc-epoch;%days
        %ele( 	a(km) 		e 			i(deg) 	W(deg) 		w(deg) 		M(deg)
        ele=[   384400      .0554       5.16    125.08     318.15      135.27];
        ele=(ones(size(ele,1),1)*[1 1 pi/180 pi/180 pi/180 pi/180]).*ele;%convert to rad
        
        %/day, angles in deg
        eledot=[ 0          0           0       .05299      .1644       13.176358];
  		eledot=(ones(size(eledot,1),1)*[1 1 pi/180 pi/180 pi/180 pi/180]).*eledot;%convert km, rad, day
    
case 5
        epoch=datenum(1976,08,10,0,0,0);
        time=dateenc-epoch;%days
        %ele( 	a(km) 		e 			i(deg) 	W(deg) 		w(deg) 		M(deg)
        ele=[   421800      .0041       .037    312.981     144.754     8.989
                671100      .0093       .464    101.087     54.425      20.865
               1070400      .0016       .186    119.841     68.990      292.375
               1882800      .0074       .253    323.265     12.668      109.158];
        ele=(ones(size(ele,1),1)*[1 1 pi/180 pi/180 pi/180 pi/180]).*ele;%convert to rad
        
        %/day, angles in deg
        eledot=[ 0          0           0       .1328       .2942       203.4889558
                 0          0           0       .03263      .07989      101.3747247
                 0          0           0       .007177     .01430      50.3176092
                 0          0           0       .001759     .003599     21.5710712];
  		eledot=(ones(size(eledot,1),1)*[1 1 pi/180 pi/180 pi/180 pi/180]).*eledot;%convert km, rad, day
         
end
     
     
if length(Bi)==1
      ele=ones(nt,1)*ele(Bi,:)+time*eledot(Bi,:);
else
      ele=ele(Bi,:)+time*ones(1,6).*eledot(Bi,:);
end
       

E=kepler(ele(:,6),ele(:,2));
q=2*atan(sqrt((1+ele(:,2))./(1-ele(:,2))).*tan(E/2))+ele(:,5);
R=[ele(:,1).*(1-ele(:,2).*cos(E)) zeros(nt,1) zeros(nt,1)];
V=[mu(CB+1)*ele(:,2)./sqrt(mu(CB+1)*ele(:,1).*(1-ele(:,2).^2)).*sin(q-ele(:,5)) sqrt(mu(CB+1)*ele(:,1).*(1-ele(:,2).^2))./R(:,1) zeros(nt,1)];%(mu*e/h)sin(q*) h/R
o=ele(:,4);I=ele(:,3);
cx=[cos(o).*cos(q)-sin(o).*cos(I).*sin(q) -cos(o).*sin(q)-sin(o).*cos(I).*cos(q) sin(o).*sin(I) ];
cy=[sin(o).*cos(q)+cos(o).*cos(I).*sin(q) -sin(o).*sin(q)+cos(o).*cos(I).*cos(q) -cos(o).*sin(I)];
cz=[sin(I).*sin(q) 								sin(I).*cos(q) 							  cos(I)];
r=[dot(R',cx')' dot(R',cy')' dot(R',cz')'];
v=[dot(V',cx')' dot(V',cy')' dot(V',cz')'];

ele(:,6)=mod(q-ele(:,5),2*pi);

r=r';v=v';
return

function E=kepler(M,e)
E=M;
err=max(abs((E-e.*sin(E)-M)./(1-e.*cos(E))));
%newton's method to solve kepler's equation
while err>1e-10
       E=E - (E-e.*sin(E)-M)./(1-e.*cos(E));
       err=max(abs((E-e.*sin(E)-M)./(1-e.*cos(E))));
end
return
