diagen.m

Stanford的SRB实验室Quantitative Seismic Interpretation的免费MATLAB程序

%DIAGEN Diagenetic trend model for moduli vs. porosity
%[k, g, rho, phi]=diagen(phic, kmin, gmin, rhomin, kfl, rhofl, shearfact, cemflag)
%
%
% Computes the bulk and shear modulus vs. porosity for a diagnetic trend.
% The trend is a composite, beginning with Dvorkin's cement line, and finishing
% with the modified upper Hashin-Shtrikman.
%
% inputs:
%      phic       = critical porosity
%      kmin, gmin = moduli of grains
%      rhomin     = density of grains
%      kfl        = modulus of fluid
%      rhofl      = density of fluid
%      shearfact  = empirical correction to shear modulus from cement 
%                                               %model (default = 0.5)
%      cemflag    = 1 (default) to use Dvorkin's contact cement; =2 
%                                               for Hashin-Shtrikman only
% outputs:
%      k          = bulk modulus along diagenetic trend
%      g          = shear modulus along diagnetic trend
%      rho        = bulk density along diagnetic trend
%      phi        = porosity array where k, g, rho are computed

if nargin < 7, shearfact = 0.5; end;
if nargin < 8, cemflag=1; end;
if nargin < 1,
     prompt = {'critical porosity',...
              'mineral bulk modulus', ...
              'mineral shear modulus', ...
              'mineral density', ...
              'Shear Reduction Factor', ...
              'Fluid Bulk Modulus',...
              'Fluid Density', ...
              'Contact Cement? (yes/no)'};
     defans = {'.4', ...
               '36e9', ...
              '40e9',...
              '2650',...
              '.5', ...
              '2.5e9',...
              '1020' ,...
              'yes'};
     answer    = inputdlg(prompt,'Diagenesis Model',1,defans);
     phic      = str2num(answer{1});
     kmin      = str2num(answer{2});
     gmin      = str2num(answer{3});
     rhomin    = str2num(answer{4});
     shearfact = str2num(answer{5});
     kfl       = str2num(answer{6});
     rhofl     = str2num(answer{7});
     choice    =         answer{8};
     if strcmp(choice,'yes'),
         cemflag=1;
     elseif strcmp(choice,'no'),
         cemflag=0;
     else
         warndlg('must answer yes or no to contact cement prompt');
         return;
     end;
end;

%
% set up porosities and densities

phi = [0:.0001:phic]; phi(1) = 1e-7; phi(end) = phic;     % array of porosities
rho = (1-phi)*rhomin + phi*rhofl;                        % array of rock bulk density

if cemflag==1,
     % compute and plot Dvorkin's cement line
     %indcement = find(phi >= max(phic-.03, 0));     %indices for range of porosities for cement
     indcement = find(phi >= .9*phic);     %indices for range of porosities for cement
     phicement = phi(indcement);
     rhocement = rho(indcement);
     [kcement, gcement]= cement(phicement,phic ,8 ,gmin ,kmin ,gmin, kmin, kfl, 2); % Dvorkin model
     gcement = shearfact*gcement;

% connect low porosity end of Dvorkin cement line with modified Hashin-Shtrikman upper bound
     indcem = find(phi< phicement(1)); 
		% indices for porosities less than Dvorkin's cement calculation
     phicem = phi(indcem); 
		% porosities less than  Dvorkin's cement calculation
     rhocem = rho(indcem); 
		% densities less than  Dvorkin's cement calculation
 
[kuc,klc,guc,glc]=hsbound(kmin,gmin,kcement(1),gcement(1),phicem/phicement(1));
% combine Dvorkin and Hashin-Shtrikman arrays into one
     kk = [kuc kcement];
     gg = [guc gcement];

elseif cemflag==0,
     kc = 1./(phic/kfl + (1-phic)/kmin);
     gc = 0;
     [kuc,klc,guc,glc]=hsbound(kmin,gmin,kc,gc,phi/phic);
     kk = kuc;
     gg = guc;
end;


k = kk;
g = gg;
% smooth a little
%k = runavg(kk, 9); k(end) = kk(end);
%g = runavg(gg, 9); g(end) = gg(end);