l_reflection_coefficients.m

地震、测井方面matlab代码,解释的比较详细

function wlog=l_reflection_coefficients(wlog,varargin)
% Function computes reflection acoustic reflection coefficients in depth and
% adds them as a new curve
% By default, the function assumes that acoustic impedance, density, interval 
% transit time and/or compressional velocity have the standard mnemonics. 
% This can be changed by temporarily redefining the standard mnemonics via 
% function arguments; e.g. {'DTp','DTCO'} specifies that the mnemonic of the 
% sonic log is 'DTCO' rather 'DTp'. Standard mnemonics are defined in 
% function "presets0" (global variable "CURVES")
% If a velocity curve does not exist it is created from the sonic (standard mnemonics) 
% Written by: E. R., February 16, 2001
% Last updated: February 17, 2001: Default mnemonic CURVES.arefl for new output curve
%
%         wlog=l_reflection_coefficients(wlog,varargin)
% INPUT
% wlog    log structure with density curve and sonic and/or velocity curve
% varargin one or more cell arrays; the first element of each cell array is 
%          a keyword string; the following arguments contains a parameter(s). 
%          Accepted keywords are:
%      'action' describes how to handle new curve. Possible values are:
%             'add'     gives an error message if a curve with this mnemonic already exists
%             'replace' gives an error message if a curve with this mnemonic does not exist
%             'add_ne'  (add --- no error) adds a new curve or replaces an existing one
%             Default: {'action','add_ne'}
%      'description'  string with curve description.
%             Default: {'description','Acoustic reflectivity'}
% OUTPUT
% wlog   log structure with reflectivity curve added (default mnemonic: 
%        string stored in "CURVES.arefl"); of course, this can be changed by 
%        temporarily redefining the standard mnemonics, 

%       Set defaults for input parameters
def.action='add_ne';
def.description='Acoustic reflectivity';

%       Decode and assign input arguments
[def,cm]=l_assign_input(def,varargin);

%       Check if necessary log curves exist
[ivp,ier]=curve_index1(wlog,cm.vp);
if ier
  [idtp,ier]=curve_index1(wlog,cm.dtp);
  if ier
    disp([' Neither sonic log "',cm.dtp,'" nor velocity log "',cm.vp,'" are present'])
    disp(' Available curves:')
    disp(cell2str(wlog.curve_info(:,1),', '))
    error(' Abnormal termination')
  end
  vp=1.0e6./wlog.curves(:,idtp);
  if strcmp(wlog.curve_info{idtp,2},'us/ft')
    vunits='ft/s';
  elseif strcmp(wlog.curve_info{idtp,2},'us/m')
    vunits='m/s';
  else
    error([' Unknown units of measurement for sonic log :',wlog.curve_info{idtp,2}])
  end
  wlog=l_curve(wlog,'add',cm.vp,vp,vunits,'Compressional velocity')
else
  vp=wlog.curves(:,ivp);
end

[iaimp,ier]=curve_index1(wlog,cm.aimp);
if ier
  [irho,ier]=curve_index1(wlog,cm.rho);
  if ier
    disp([' Neither density log "',cm.rho,'" nor acoustic impedance log "',cm.aimp,'" are present'])
    disp(' Available curves:')
    disp(cell2str(wlog.curve_info(:,1),', '))
    error(' Abnormal termination')
  end
  
  aimp=vp.*wlog.curves(:,irho);

else
  aimp=wlog.curves(:,iaimp);
end

%       Compute reflection coefficients
temp=diff(aimp)./(aimp(1:end-1)+aimp(2:end));
if strcmp(wlog.curve_info{1,2},'m') | strcmp(wlog.curve_info{1,2},'ft')
  refl=temp.*(vp(1:end-1)+vp(2:end))./diff(wlog.curves(:,1));
  if (strcmp(wlog.curve_info{1,2},'m') & strcmp(l_gu(wlog,cm.vp),'m/s')) | ...
     (strcmp(wlog.curve_info{1,2},'ft') & strcmp(l_gu(wlog,cm.vp),'ft/s'))
    refl=refl/4000;

  elseif strcmp(wlog.curve_info{1,2},'m') & strcmp(l_gu(wlog,cm.vp),'ft/s')
    refl=refl*0.0000762;

  elseif strcmp(wlog.curve_info{1,2},'ft') & strcmp(l_gu(wlog,cm.vp),'m/s')
    refl=refl/304.8; 

  else
    error([' Unknown depth (',wlog.curve_info{1,2},') or velocity (',l_gu(wlog,cm.vp),') units'])
  end
 
elseif strcmp(wlog.curve_info{1,2},'s') | strcmp(wlog.curve_info{1,2},'ms')
  if wlog.step > 0
      refl=temp/wlog.step;
  else
      refl=temp./diff(wlog.curves(:,1));
  end
  if strcmp(wlog.curve_info{1,2},'s')
      refl=refl/1000;
  end
 
else
  error(' Depth units must be either "m", "ft", "s", or "ms"')
end

%       Add reflectivity to (replace reflectivity in) log structure 
wlog=l_curve(wlog,def.action,cm.arefl,[refl;0],'n/a',def.description);