function [residual] = wave_fit_error_LH82(params,Xr,Zr, mu_range)
%discretize the parameter space for mu, and compute the minimum distance
%for each data point to the discretized LH82 curve. Then sum these distances to
%compute the total residual

%Inputs
%params=[mag, thetaD, X0, Z0]
%Xr = vector of X data points;
%Zr = vector of Z data points

%define discretized parameter space


%number of data points
nData=length(Xr);
%run LH82
[X,Z]=LH82functionFF(params(1),params(2),params(3),params(4),params(5),mu_range);
%compute residuals

nmu = length(X);

npts = 1;

dist=pdist2([X Z],[Xr Zr],'euclidean','smallest',npts);

distsq_avg = sum(dist.^2,1)/npts;
%distsq_avg = sum(dist,1)/nData;

residual =  sqrt(sum( distsq_avg)./nmu );


end
