% Fit capillary pressure data using Thomeer model
% At end of run, press ENTER to view plot window
% Input file thomd.dat has saturations as wetting saturations
% Saturations in code are Hg saturations, except for thomd.dat
% 0.01 < S/S(inf) < 0.97
% x = S/S(inf)
% y = Pc/Pd
clear all
clf reset
ninc=100;
inc=(log10(2)-log10(-log10(.97)))/ninc;
x=zeros(1,ninc+1);
ex=log10(-log10(.97));
x(1)=10^(-10^ex);
for i=2:ninc+1
  ex=ex+inc;
  x(i)=10^(-10^ex);
end
pd=input('pd = ?');
sinf=input('SHg(inf) = ?');
g=input('G =  ?');
y=zeros(size(x));
s=zeros(size(x));
pc=zeros(size(x));
for i=1:ninc+1
  y(i)=exp(2.303*g/(-log(x(i))));
  s(i)=sinf*x(i);
  pc(i)=pd*y(i);
  if pc(i) > 2000.
     s(i)=[];
     pc(i)=[];
  end
end
%logy(s,pc)
semilogy(s,pc,'k-')
axis([1.e-2 1 .1 1.e3])
set(gca,'xdir','reverse');
xlabel('Mercury Saturation')
ylabel('Capillary Pressure')
title(['Pd=',num2str(pd),' Sinf=',num2str(sinf),' G=',num2str(g)])
[fid, message]=fopen('thomd.dat','r');
load thomd.dat
status=fclose(fid);
xdat=1-thomd(:,1);
ydat=thomd(:,2);
hold
semilogy(xdat,ydat,'kx')
%  print fig1 -dbitmap
