hvs_rpt.m

function hvs_rpt( p, t, s, d, ttle )
%HSV_RPT - Produces the final figure and results for the Hvorslev (1951) model
%                
% Syntax:  hvs_rpt( p, t, s, d, ttle )
%
%   p = parameters of the Hvorslev model
%   t = measured time
%   s = measured drawdown
%   d(1) = rw = Radius of well screen
%   d(2) = rc = Radius of the casing
%   ttle = Title of the figure 
%
% Description:
%   Produces the final figure and results for the Hvorslev (1951) model.
%
% See also: hvs_dmo, hvs_dim, hvs_gss
%


[t,i]=sort(t);
s=s(i);
clf
if(nargin==4) % Default value for d if not given by user
  ttle='Slug test';
end

% Fit the model and calculate interval of confidence
rw=d(1); rc=d(2);
p0=hvs_gss(t,s);
[p,residual,jacobian]=nlinfit(t,s,'hvs_dim',p0);
ci=nlparci(p,residual,jacobian);
% 
T=0.5*rc^2/p;
Tmax=0.5*rc^2/ci(1,2);
accu=abs((Tmax-T)/T*100);

% Calculate the drawdown and derivative
[td,sd]=ldiffs(t,s,20);
sd=-sd;
tplot=logspace(log10(t(1)),log10(t(end)));
[sc,delta] = nlpredci('hvs_dim',t,p,residual,jacobian);
[tdc,dsc]=ldiff(t,sc');
dsc=-dsc;

% Select the area for the plot
xo=0.1; yo=0.5;
dy=0.4; dx=0.8;

% Plot the data and the calculated drawdown
clf
hax=axes('position',[xo,yo,dx,dy]);
plot(log10(t),s,'o',log10(td),sd,'+',log10(t),sc,'-',log10(tdc),dsc,'-.')
hold on
plot(log10(t),sc+delta,'--b',log10(t),sc-delta,'--b')
hold off
set(hax,'Visible','off')
axis([log10(t(1)),log10(t(end)),0,1])
xl=get(hax,'XLim');
yl=get(hax,'YLim');
legend('Drawdown','Derivative','Model','Model derivative')

% Coordinate transformation to position the legend
ax=(xl(2)-xl(1))./dx;
ay=(yl(2)-yl(1))./dy;


% Text of the legend
legend1=sprintf('Fitting parameters:\n t_0: %0.2g s\n\nHydraulic parameters:\n Transmissivity T:\n %3.1e m^2/s +/- %.2f per cent',p,T,accu);
ht=text(ax*(0.1-xo)+xl(1),ay*(yo-0.1-yo)+yl(1),legend1);
set(ht,'VerticalAlignment','top')
set(ht,'HorizontalAlignment','left')

ht=text(ax*(0.5-xo)+xl(1),ay*(yo-0.12-yo)+yl(1),'Spherical Model');
set(ht,'VerticalAlignment','top')
set(ht,'HorizontalAlignment','left')

legend2=sprintf('Test data:\n Well radius: %0.2g m\n Casing radius: %0.2g m',rw,rc);
ht=text(ax*(0.5-xo)+xl(1),ay*(yo-0.2-yo)+yl(1),legend2);
set(ht,'VerticalAlignment','top')
set(ht,'HorizontalAlignment','left')

copyright=hycoop;
hc=text(ax*(xo+dx-xo)+xl(1),ay*(0.1-yo)+yl(1),copyright);
set(hc,'HorizontalAlignment','right')
set(hc,'FontSize',8)

% Draw the axes in log scale
halog=axes('position',[xo,yo,dx,dy]);
set(halog,'Color','none')
set(halog,'XScale','log')
set(halog,'YLim',[yl(1),yl(2)])
set(halog,'XLim',[10^xl(1),10^xl(2)])
set(halog,'Box','on')

title(ttle)
xlabel('Time in seconds')
ylabel('Relative Drawdown')