function [afit,bfit,mod_r,mod_u] = Model_fit(rvals,dvals,R0,nusc,num,verbose)
%Continuum_fit fits a displacement profile to the continuum skin model
% Calling sequence 
% [alp,beta,mod_r,mod_u]=Continuum_fit(rvals,d.rad_d,num)
% Input:
% rvals - n by 1 radial position list
% dvals - n by 1 radial displacment list
% R0 = SC radius
% nusc - Poisson's ratio of SC
% num - optional - sets the size of the parametric alpha and beta fit.
% Default is 100 for each parameter
% Output:
% afit - best fit value of alpha
% bfit - best fit value for beta
% mod_r - best fit model radial position
% mod_u - best fit model radial displacments
% verbose - set to 1 to plot fits as you go
% for speed, set to 0

% Code associated with journal article 'Measuring and modeling contractile
% drying in human stratum corneum' by X. Liu and G. K. German for the Journal
% of Visualized Experiments 
% August 15 2016

if nargin ==2
    num = 100;
end
dtop = find(dvals == max(dvals)); % Find maximum value of radial displacment
dvals2 =dvals(1:dtop);% region to fit displacement profile to 
rvals2 = rvals(1:dtop);
alpm = linspace(0.07,0.001,num); %Typical range of Alpha values for Human SC
betam = linspace(0.02,10,num); % Typical range of Beta values for Human SC
for imi = 1:length(alpm);
    for ini = 1:length(betam);
        ai = alpm(imi);bi=betam(ini);
        I1_r = besseli(1,bi*rvals2/R0);
        A = 1/(bi*besseli(0,bi)-(((1-(2*nusc))/(1-nusc))*besseli(1,bi))); 
        uvalfit = ai*R0*A*I1_r;
        if verbose==1 %plot the figure to determine the range
        plot(rvals,dvals,'ko');hold on
        plot(rvals2,uvalfit,'r');hold off
        axis([0 max(rvals) 0 1.5*max(dvals)]); 
        drawnow
        end
        lsqpoly(imi,ini) = lsfit(dvals2,uvalfit);% find the least square fit
    end
end
[rpi cpi] = find(lsqpoly==min((lsqpoly(:)))); 
afit=alpm(rpi);
bfit=betam(cpi);
I1_r = besseli(1,bfit*rvals2/R0);
A = 1/(bfit*besseli(0,bfit)-(((1-(2*nusc))/(1-nusc))*besseli(1,bfit))); 
mod_u = afit*R0*A*I1_r;
mod_r = rvals2;
end


