%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [nx,ny] = evaluate_normal(coefs,s)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This function finds the unit normal vector at each point in s on
% the B-splines in coefs.  
%
% Input:
% coefs: A 2 x (L+3) x N matrix, defining N B-splines.  
% s: A ns x N matrix, where normals at each s(j) are to be evaluated
% on each spline, defined by coefs(:,:,i).  
%
% Global variables required:
% BS_sigma and G_sigma: precomputed matrices needed for *any* cubic
% B-spline with L curve pieces.  These are precomputed to speed up
% computation.
%
% Output:
% [nx,ny]: unit normal vectors of the splines at the desired
% points.  Each is a ns x N matrix.  
%
% No global variables set.
%
% Most of the work is done by evaluate spline_prime, which computes
% the tangents to the splines at the desired points.  The unit
% vectors perpendicular to the tangents are computed from the
% tangent vectors returned by evaluate_spline_prime.  
%
% Functions called by this function:
% evaluate_spline_prime
%
function [nx,ny] = evaluate_normal(coefs,s)

L = size(coefs,2) - 3;
N = size(coefs,3);

% reshape coefficients to be (2*nspline) x (L+3)
coefsallx = permute(coefs(1,:,:),[3,2,1]);
coefsally = permute(coefs(2,:,:),[3,2,1]);

% Find the tangent 
% xprime and yprime are each ns x nsplines
xprime = evaluate_spline_prime(coefsallx,s);
yprime = evaluate_spline_prime(coefsally,s);

% Normal is perpendicular to tangent
% nx and ny are each ns x nsplines
nx = yprime;
ny = -xprime;

% Normalize to unit length
l = sqrt(nx.^2 + ny.^2);
nx = nx./l;
ny = ny./l;