function V = ortho(x0,mue)
%
% input:
% 	x0(N,4)		pseudo orbit
% 	mue		parameter for mass ratio of the primaries
%
% output:
% 	V(N,16)		V(n,:) contains orthonormal basis for x0(n,:)
%

N=size(x0,1);
[v1,v2] = funcgrad(x0,mue);
for j=1:N,
  [Q,R]=qr([v1(j,:)', v2(j,:)']);
  J=diag([sign(R(1,1)),sign(R(2,2)),1,1]);
  Q=Q*J;
  R=J*R; 
  V(j,:)=reshape(Q,1,16);
end

return

%
% written   Sept/10/02   D.Stoffer
%








