function penalty = multivariable(x)

fmmob=x(1);
fmdry=x(2);
epdry=x(3);

krw0=0.79;
krg0=1.0;
S_wc=0.07;
S_gr=0;
nw=1.96;
ng=2.29;
visw=0.001*1e3; %cp
visg=0.00002*1e3; %cp

fg_range=[0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.95 0.99];
visfoam_range=[232.4612852 257.2143292 345.125721 380.9955291 420.4661673 333.9610169 263.0471435 182.4614989 122.7398356 76.12673182 34.15912117];
w=[1 1 1 500 500 500 1 1 1 1 1]; % weighting paramters

S_w=zeros(1,length(fg_range));
vis_foam=zeros(1,length(fg_range));

for i=1:length(fg_range)
S_w(i)=fzero(@(Sw) 1./(1+(krw0*((Sw-S_wc)/(1-S_gr-S_wc)).^nw).*visg./...
    ((krg0*(1-(Sw-S_wc)/(1-S_gr-S_wc)).^ng.*(1./(1+fmmob*(0.5+(atan(epdry*(Sw-fmdry)))/pi)))).*visw))-fg_range(i), [S_wc, 1-S_gr]);
end

Fwater=0.5+(atan(epdry*(S_w-fmdry)))/pi;
Fsurf=1;
krw=krw0*((S_w-S_wc)/(1-S_gr-S_wc)).^nw;
krg=krg0*(1-(S_w-S_wc)/(1-S_gr-S_wc)).^ng.*(1./(1+fmmob*Fwater*Fsurf));
vis_foam=1./(krw/visw+krg/visg);

leastsqr=0;
for k=1:length(fg_range)
    leastsqr=leastsqr+w(k).*((vis_foam(k)-visfoam_range(k))/visfoam_range(k))^2;
end

if fmdry>=S_wc
    phi=0;
else
    phi=1;
end

miu=0.1;
penalty=leastsqr+miu*((fmdry-S_wc)^2*phi);

% to run the multi-variable search
%A=fminsearch('multivariable',[1, 0.07, 1],optimset('MaxFunEvals',5000))
