head_delay(hrad,src,mic,aim,sos);

function extratime = head_delay(hrad,src,mic,aim,sos);
%
%function head_delay(hrad,src,mic,aim,sos)

%  Computes the delay(seconds) for having to go around a spherical 
%  head of radius r. It really computes the extra distance the
%  sound has to travel for those directions in which the shape of the
%  head affects the path
%  INPUTS:
%  hrad      = the radius of the sphere(head) (meters)
%  src       = [x,y,z] of the mouth in absolute coordinates
%  mic       = [x,y,z] of the microphone in absolute coordinates
%              it MUST lie on the surface of the sphere.
%  aim       = [x,y,z] of the aiming point of the mouth in absolute coordinates
%  sos       = speed of sound
%  OUTPUT:   
%  extratime = the time delay(seconds) incurred by passing around 
%              rather than through the head

% First compute the coordinates of the center of the sphere

% src=src';
q = aim - src;
unit = q / sqrt(q * q');
center = src - unit * hrad;

direct = sqrt((src - mic) * (src - mic)');
f = 1.0 / (hrad*hrad);
nc = center - src;    % adjusted center
nmx = (mic - center) * nc' * f;
d2 = (mic - center)*(mic - center)' * f;
nmy = sqrt(d2 - nmx^2);
if nmx < -1
   dist = sqrt((nmx+1)*(nmx+1) + nmy*nmy) * hrad;
elseif d2 >= 1
   g = 1/d2; 
   f = sqrt(g);
   ux = nmx*f; 
   uy = nmy*f;
   g = sqrt(1-g);
   x = ux*f - uy*g;
   y = uy*f + ux*g;
   dist = sqrt((x-nmx)*(x-nmx)+(y-nmy)*(y-nmy)) * hrad;
   ang = atan2(x,y) + pi/2;
   dist = dist +  ang * hrad;
   
  
   
else
    dist = -1;
end   
extratime = (dist-direct) / sos;
return