function fig=trajplot(X,B,Nrev)
%trajplot([Te Rm,Tm],[CB Bi],Nrev)
%Te = Time of encounter at body
%[Rm Tm] = Position and time of maneuver
%CB = designation for central body (0=sun,3=Earth,...)
%Bi = designation for bodies to encounter
%Nrev = # of revolutions on arce connecting events

%specify figure # & set up axes
fig=101;
figure(fig);
clf
hold on
axis equal
axis off

%gravitairtonal parameter
mu=planetdata(B(1),0,1);

%target trajectory
[Vi1 Vi2 DV ind]=trajcalc(X,B,Nrev);
%distinguish between bodies and maneuvers
ie=ind{1};
im=ind{2};
Ne=length(ie);
Nm=length(im);
%Dates
JD([ie im])=[X(1:Ne) X(Ne+4*[1:Nm])];
dJD=diff(JD);
%calulate initial state for each arc
for i=1:length(dJD)
    if any(ie==i)%planet
       [R0 V0]=ephem(B(1),B(find(ie==i)+1),JD(i));
       V0=V0+Vi2(:,find(ie==i));
    else%maneuver
       R0=X(Ne-4+4*find(im==i)+[1:3])';
       V0=DV(1:3,find(im==i));
    end
    %propigate as two-body system
    Xt=twobody([R0;V0],[0:1:dJD(i)]*86164,mu);
    plot3(Xt(1,:),Xt(2,:),Xt(3,:),'k','linewidth',1.3)%plot path
    plot3(Xt(1,1),Xt(2,1),Xt(3,1),'k.','markersize',20)%plot event
    text(Xt(1,1),Xt(2,1),Xt(3,1),datestr(JD(i),26),'fontsize',11,'horizontalalignment','left',...
        'verticalalignment','bottom')%add date
end
plot3(Xt(1,end),Xt(2,end),Xt(3,end),'k.','markersize',20)%plot last event
text(Xt(1,end),Xt(2,end),Xt(3,end),datestr(JD(end),26),'fontsize',11,'horizontalalignment','left',...
    'verticalalignment','bottom');%last date
axlim=axis;%get size of axis for trajectory only
%add planets
for b=B(2:end)
    %position and velocity
    [R0 V0 a]=ephem(B(1),b,datenum('1/1/2000'));
    period=2*pi*sqrt(a(1)^3/mu);
    %propigate (# of points inv. prop. to a
    Xt=twobody([R0;V0],linspace(0,period,round(period*1.5e7/a(1)/86164)),mu);
    plot3(Xt(1,:),Xt(2,:),Xt(3,:),'.')
end
% show vernal equinox
text(axlim(2),0,0,'\Upsilon\rightarrow','fontsize',13)
%hold off
%axis(axlim*1.2)%close in on trajectory if desired