25 lines
1.1 KiB
Matlab
25 lines
1.1 KiB
Matlab
function [vwc]=dobsen_inverse(f,T,er,bd,vsand,vclay,x0)
|
|
%--------------------
|
|
alpha=0.65;
|
|
sd=2.65;
|
|
% dcs=4.7;
|
|
dcs=(1.01+0.44*sd)^2-0.062;
|
|
dc0=0.008854;
|
|
dcw0 = 88.045-0.4147.*T+6.295e-4.*T.^2+1.075e-5.*T.^3; % Refer to [3]
|
|
tpt=0.11109-3.824e-3.*T+6.938e-5.*T.^2-5.096e-7.*T.^3; %2*pi*tao, Refer to [4]
|
|
dcwinf=4.9;
|
|
if f>=1.4
|
|
sigma = -1.645+1.939*bd-0.0225622*vsand +0.01594*vclay; % Refer to [1]
|
|
else
|
|
sigma = 0.0467+0.2204*bd-0.004111*vsand +0.006614*vclay;% Refer to [2]
|
|
end
|
|
dcwr=dcwinf+((dcw0-dcwinf)./(1+(tpt.*f).^2));
|
|
% dcwi=(tpt.*f.*(dcw0-dcwinf))./(1+(tpt.*f).^2)...
|
|
% +sigma*(1.0-(bd/sd))./(8.0*atan(1.0)*dc0.*f.*vwc);% Refer to [1]
|
|
betar=1.2748-0.00519*vsand-0.00152*vclay;
|
|
betai=1.33797-0.00603*vsand-0.00166*vclay;
|
|
% dcsr=(1.0+(bd/sd)*((dcs^alpha)-1.0)+(vwc.^betar)*(dcwr^alpha)-vwc).^(1/alpha);% Result real part
|
|
% dcsi=(vwc.^(betai/alpha)).*dcwi;%Result imaginary part
|
|
vwc=fminsearch(@(vwc)abs(sqrt(((1.0+(bd/sd)*((dcs^alpha)-1.0)+(vwc.^betar)*(dcwr^alpha)-vwc).^(1/alpha))^2+((vwc.^(betai/alpha)).*((tpt.*f.*(dcw0-dcwinf))./(1+(tpt.*f).^2)...
|
|
+sigma*(1.0-(bd/sd))./(8.0*atan(1.0)*dc0.*f.*vwc)))^2)-er),x0);
|
|
end |