function [V alpha] = relative_vel(p,n,name,inc_wave)

% Returns specified bodies speed and direction relative to local flow in 
% body-fixed coordinate system.
%  Inputs:  p = rigid body specification.
%           n = timestep,  if n >= 0, uses current state vector, if n<0,
%           it specifies the timestep from the solution to use.
%           name = name of body that incident wave field is centered on.
%           inc_wave = incident wave object.


if (n >= 0)  %Use current state vector for position and velocities
  pos1 = p.sv(1:6);
  vel1 = p.sv(7:12); 
  sv2 = get(name,'sv');
  pos2 = sv2(1:6);
  vel2 = sv2(7:12);
  clear sv2;
else    %Use specified time-step for position and velocities
  n = -n;
  pos1 = p.x(n,:)';
  vel1 = p.xdot(n,:)';
  x2 = get(name,'x');
  v2 = get(name,'xdot');
  pos2 = x2(n,:)';
  vel2 = v2(n,:)';
  clear x2 v2;
end


dt = p.dt;

xoa = p.x0;
xob = get(name,'x0');

%Stuff for including incident wave effects.
vel_i = zeros(3,1);
if (isobject(inc_wave))     % Get incident wave velocity at body coordinate center (relative to
                            % surface buoy, specified as 'name')
  [vi1 vi2 vi3] = vel(inc_wave,(xoa(1)+pos1(1))-(xob(1)+pos2(1)),(xoa(3)+pos1(3))-(xob(3)+pos2(3)),n*dt);  
  vel_i(1) = vi1;
  vel_i(2) = vi2;
  vel_i(3) = vi3;
end

vel1(1:3) = vel1(1:3)-vel_i;

V2 = vel1(1).^2+vel1(3).^2;

theta = pos1(5);
beta = atan2(-vel1(3),vel1(1));

alpha = beta-theta;
V = sqrt(V2);

if(alpha > pi)
  alpha = alpha-2*pi;
end
if(alpha < -pi)
  alpha = alpha+2*pi;
end


%disp(['V = ' num2str(V)]);
%disp(['theta = ' num2str(theta*180/pi)]);
%disp(['beta = ' num2str(beta*180/pi)]);
%disp(['alpha = ' num2str(alpha*180/pi)]);



end