-
Notifications
You must be signed in to change notification settings - Fork 11
/
jcb_gss.m
50 lines (41 loc) · 919 Bytes
/
jcb_gss.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
function p = jcb_gss( t, s, tmin, tmax)
%JCB_GSS - Linear Least Square fitting of Jacob solution
%
% Syntax: p = jcb_gss(t,s,tmin,tmax)
%
% p(1) = a = slope of Jacob straight line
% p(2) = t0 = intercept of the Jacob straight line
%
% t = time
% s = drawdown
% tmin = optional argument: start of the fitting period
% tmax = optional argument: end of the fitting period
%
% Description:
% First guess for the parameters of Jacob solution
%
% See also: jcb_dim, jcb_dmo, jcb_rpt
%
if(nargin<=2) % Default value
tmin=t(1);
end
if(nargin<=3) % Default value
tmax=t(end);
end
i=find(isnan(s)); % The inversion will crash if there are NaN
s(i)=[];
t(i)=[];
i=find(t<tmin);
t(i)=[];
s(i)=[];
i=find(t>tmax);
t(i)=[];
s(i)=[];
%%%%%%%%%%%%%%%% This is the only calculation !
g=[log10(t),ones(size(t,1),1)];
p=inv(g'*g)*g'*s;
%%%%%%%%%%%%%%%%%
a=p(1);
c=p(2);
t0=10^(-c/a);
p(2)=t0;