a=5;
b=4;
v=3;
k=sqrt(v^2-1);
f=(0:100)/50*pi;
n=length(f);
xx=a*cos(f);yy=b*sin(f);
x1=reshape([nan(1,n);xx;xx-(a*sin(f)+k*b*cos(f))*k/v^2.*(a/b*sin(f).^2+b/a*cos(f).^2)],1,n*3);
y1=reshape([nan(1,n);yy;yy-(a*k*sin(f)-b*cos(f))*k/v^2.*(a/b*sin(f).^2+b/a*cos(f).^2)],1,n*3);
x2=reshape([nan(1,n);xx;xx+(a*sin(f)-k*b*cos(f))*k/v^2.*(a/b*sin(f).^2+b/a*cos(f).^2)],1,n*3);
y2=reshape([nan(1,n);yy;yy-(a*k*sin(f)+b*cos(f))*k/v^2.*(a/b*sin(f).^2+b/a*cos(f).^2)],1,n*3);
plot(xx,yy,'g',x1,y1,'b',x2,y2,'r');daspect([1 1 1])