function [tau,gam,varr,earr]=fit_poschel(inputdata,modelfile,D,R,rhos,tauest,gamest)

% Finds the best fit viscous relaxation time and surface energy
% to data for normal restitution e versus relative impact velocity of two identical spheres
% using the model of N.V. Brilliantov, F. Spahn, J.-M. Hertzsch and T.
% P\"oschel:``Model for collisions in granular gases,'' PRE 53, 5382-5392 (1996).
% The model is implemented in Poschel_vdW_TN_dim-.nb, which generates
% model curves of e versus v/vc for different R/R0.
% The theory for adhesion is that of Johnson.
%
% MYL 11/11/06. Modified 8/3/08 to output the model curve.
%
%	input	-	inputdata	=	string with experimental data file name
%			-	modelfile	=	string of file name with model predictions of e
%                               compatible with ranges of R/R0 and v/vc generated by 
%                               [rr0,vvc]=meshgrid(3:(3/39):6,-6:.05:0);
%                               modify Poschel_vdW_TN_dim-.nb if other ranges are needed
%           -   D           =   (3/2) (1-nu^2)/E (1/Pa)
%           -   R           =   sphere radius (m)
%                               where nu =Poisson's ratio and E = Young's modulus (Pa)
%           -   rhos        =   material density (kg/m^3)
%           -   tauest      =   estimate of the best fit relaxation time (s)
%           -   gamest      =   ditto surface energy (J/m^2)
%
%	output	-	tau         =   best fit relaxation time (s)
%           -   gam         =   best fit surface energy (J/m^2)
%           -   varr        =   array of impact velocities (m/s)
%           -   earr        =   array of restitution coefficients (-)
%
%   global variables
%   
%           -   indat       =   experimental data in file inputdata
%                                   column 1: values of v (m.s)
%                                   column 2: values of e (-)
%                                       values of e for each element pair in rr0,vvc 
%                                       are entered in modelfile and called indat
%           -   model       =   model data in file modelfile
%                                   columns:    constant log(R/R0)
%                                   rows:       constant log(v/vc)
%           -   R0root      =   see definition below
%           -   vcroot      =   see definition below
%           -   rr0         =   array of log(R/R0) generated by meshgrid
%                                   along columns:	invariant   log(R/R0)
%                                   along rows:     varying     log(R/R0)
%           -   vvc         =   array of log(v/vc) generated by meshgrid;
%                                   along columns:	varying     log(v/vc)
%                                   along rows:     invariant   log(v/vc)
%           -   sizedat     =   number of data pts
%
%  Note: array dimensions of rr0, vvc and indat are identical.
%           
%   to access the proper directory, type
%   cd ~/Michel_files/Research/granular_flows_impact/project_Cherry/PDEs/optimization
%
% Example:
%
% clear all;[tau,gam]=fit_poschel('data','restit',7.27579E-12,0.0015,7920,1.4E-7,92);
%

global  indat model R0root vcroot rr0 vvc sizedat

%   internal variables
%
%           -   x           =   line vector of [tau,gam]
%           -   vc          =   reference velocity (m/s)
%

indat       =   load(inputdata)                 ;
model       =   load(modelfile)                 ;

sizedat     =   size(indat,1)                   ;

[rr0,vvc]   =   meshgrid(3:(3/39):6,-6:.05:0)   ;

if size(rr0) ~= size(model) | size(vvc) ~= size(model)
    disp(['warning: incompatible file sizes'])
    return
end


R0root  =   2^(15/14)*pi^(5/7)/3^(2/7)/5^(3/7)*...
                ( gamma(2/5)/gamma(9/10) )^(6/7)*...
                    D^(2/7)*rhos^(3/7)*R                ;   % R/R0 = R0root/tau^(6/7)/gam^(1/7)
                
vcroot  =   2^(51/14)*pi^(13/14)/3^(4/7)/5^(6/7)*...
                ( gamma(2/5)/gamma(9/10) )^(5/7)*...
                    D^(4/7)/rhos^(1/7)                  ;   % vc = vcroot*gam^(5/7)/tau^(5/7)
                
x       =	fminsearch(@minerror,[tauest,gamest])       ;

tau     =   x(1,1)  ;
gam     =   x(1,2)  ;

vc      =   vcroot*(gam/tau)^(5/7)          ;

varr    =   exp(-6:.05:0).*vc              	;
earr    =   max(0,interp2(rr0,vvc,model,...
                log(R0root/tau^(6/7)/gam^(1/7)),(-6:.05:0),'cubic'))    ;

plot((-6:.05:0),max(0,interp2(rr0,vvc,model,log(R0root/tau^(6/7)/gam^(1/7)),...
    (-6:.05:0),'cubic')),'r',log(indat(:,1)./vc),indat(:,2),'o')    %   plot the best fit and the data
xlabel('v/vc')
ylabel('e')
title([' gam = ',num2str(gam),'(J/m^2); tau = ',num2str(tau*1e6),' (microsec)'])
