clear all
clc
mmm=150;
fmmob_range=linspace(28000,38000,mmm);
fmdry_range=linspace(0.07,0.08,mmm);

global S_wc S_gr krw0 krg0 nw ng visw visg epdry     

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

for I=1:mmm     
    fmmob=fmmob_range(I); %fmmob as x1 (x coordinate)
    for J=1:mmm
        fmdry=fmdry_range(J); %fmdry as x2 (y coordinate)
        global fmmob fmdry    

S_w=fzero(@calSw,[S_wc 1]); % calculate water saturationat a given fg

Fwater=0.5+(atan(epdry*(S_w-fmdry)))/pi;

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));

vis_foam=1./(krw/visw+krg/visg);
VISFOAM(I,J)=vis_foam;
    end
end

figure('Position',[1 31 1280 696]);
h = gca; 
set(h,'FontSize',30);
[AA,BB]=contour(fmdry_range,fmmob_range,VISFOAM,'k','LineWidth',2)
clabel(AA,BB,'FontSize',20,'LabelSpacing',1000)
xlabel('fmdry')
ylabel('fmmob')
%title('Contour of foam apparent viscosity (cp)')