function [Vinf1,Vinf2,DV,index,E,F,G,H]=trajcalc(X,B,Nrev)
%[Vin,Vout,DV,E,F,G,H]=trajcalc([Te Rm,Tm],[CB Bi],Nrev)
%patched conic method
%OUTPUTS
%Vin,Vout = incoming ant outgoing Vinf vectors
%           [Vin(1) = incoming at 2nd body, Vout(1) = out at 1st body]
%DV = Velocity after & before impulse [DV(1:3)=V(after), DV(4:6)=V(before)]
%index = [index(bodies) index(maneuvers)], index for path event order
%E,F,G,H = matrices for partial derivatives
%dV0/dRF = E, dV0/dTF = -E*VF 
%dV0/dR0 = F, dV0/dT0 = -F*V0 + accel(0)
%dVF/dRF = G, dVF/dTF = -G*VF + accel(F)
%dVF/dR0 = H, dV0/dRF = -H*V0, 

%INPUTS
%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

DV=[];%initialize
CB=B(1);B=B(2:end);%distinguish bodies
muC=planetdata(CB,0,1);%gravitational parameter
ne=length(B);nm=(length(X)-ne)/4;%# of events & maneuvers
if nargin <3;Nrev=zeros(1,ne+nm-1);end%default Nrev

if nm>0%maneuvers present
    Te=X(1:length(B));%Time of body encounter
    Xm=reshape(X(ne+1:end),4,nm);
    Rm=Xm(1:3,:);%maneuver position
    Tm=Xm(end,:);%maneuver time
    T=union(Te,Tm);%all times
    ie=find(ismember(T,Te));im=find(ismember(T,Tm));%distinguish between bodies and maneuvers
    [Re Ve]=ephem(CB,B,Te);%body position and velocity
    R(:,ie)=Re;R(:,im)=Rm;%all positions
    TOF=diff(T)*86164;%Time of flight on arcs
	[V1 V2]=lambert(R(:,1:end-1),R(:,2:end),TOF,muC,Nrev);%solve Lambert's problem
	Vinf1=V2(:,ie(2:end)-1)-Ve(:,2:end);Vinf2=V1(:,ie(1:end-1))-Ve(:,1:end-1);%Vinf vectors
    DV=[V1(:,im);V2(:,im-1)];%Velocity ad maneuver
    index={ie,im};%ouput path event order
else
    T=X;%body event times
    [R Ve]=ephem(CB,B,X);%body position and velocity
	TOF=diff(T)*86164;%Time of Flight on arcs
	[V1 V2]=lambert(R(:,1:end-1),R(:,2:end),TOF,muC,Nrev);%solve Lambert's problem
	Vinf1=V2-Ve(:,2:end);Vinf2=V1-Ve(:,1:end-1);%Vinf vectors
	DV=zeros(6,0);%no DV
    index={[1:ne],[]};%ouput path event order (no maneuvers)
end
    
if nargout>4
%Calculate partials    
for i=1:length(TOF)
%state transition matrix
[phi]=stm([R(:,i);V1(:,i)],TOF(i),muC,[R(:,i+1);V2(:,i)]);
%define outputs
E(:,:,i)=inv(phi(1:3,4:6));
F(:,:,i)=-E(:,:,i)*phi(1:3,1:3);
G(:,:,i)=phi(4:6,4:6)*E(:,:,i);
H(:,:,i)=phi(4:6,1:3)-G(:,:,i)*phi(1:3,1:3);
end
end