function v=vel0(x,y,a,b,vmax,vmin,theta0)
% compute an angle position of the point
al=atan2(y./b,x./a);
v = (-pi+theta0 <= al & al < -theta0).*vmax+ ...
    ((pi >= al & al >= -theta0) | al < -pi+theta0).*(vmin+...
    2.*abs((al>=-pi/2).*al + (al<-pi/2).*(al+2*pi)-pi/2).*(vmax-vmin)./pi);
v = v.*sqrt(x.^2./a.^2+y.^2./b.^2);