bkusrunavgnew.m

function [vpeff,vseff,rhoup]=bkusrunavgnew(z,vp,vs,rho,win)
%        [invqeff,veff,rhoup]=bkusrunavgvisco(z,vp,invqp,rho,win,om,omega0)
%
% Running Backus average of log curves to yield an array of effective anisotropic
% stiffness tensors, c(6,6,N), where N is the number of depth points. At a
% particular depth point, j, the stiffness tensor is c(6,6,j).  
% Note that the output stiffness array is 3-dimensional.  If an application
% requires a 2-dimensional array at a single depth point, then use Matlab
% function SQUEEZE, e.g.  
%          squeeze(c(6,6,j))
% yields a 6x6 array.
% The Backus average is computed in a running window, and it is assumed that the
% medium consists of thin isotropic layers.
%
% input:   z       - Vector of depths, as in a well log. Can be in any units.  
%                    Depth points should be equally spaced.
%          vp      - Vector of P velocity, as in a well log -- same length
%                    as vector of depths.
%          invqp   - Vector of qp values -- same length
%                    as vector of depths.
%          rho     - Vector of densities, as in a well log -- same length
%                    as vector of depths.
%          win     - length of window for averaging - in the same length units as z
%          om      - Frequency of the upscaled log (seismic frequency)
%          om0     - Frequency of the input log (well log frequency)
% output:  c       - 6x6xN stiffness array expressed in Voigt notation -- as 
%                    defined in The Rock Physics Handbook
%          rho     - Average density at each depth point, averaged over the
%                    same window as the elastic constants
%
% see also BKUS, BKUSLOG, BKUSC, BKUSRUNAVG

% Modified from bkusrunavg by Vishal Das, October 2016
%

nwin = floor(win/mean(diff(z)));

% Create layer thicknesses from log depths
dz = diff(z); dz(end+1) = dz(end);

M = rho.*(vp.^2);
G = rho.*(vs.^2);
Meff = zeros(1,length(z));
Geff = zeros(1,length(z));
vpeff = zeros(1,length(z));
vseff = zeros(1,length(z));
rhoup = zeros(1,length(z));
Meff(1:floor(win/2)) = M(1:floor(win/2)); 
vpeff(1:floor(win/2)) = vp(1:floor(win/2)); 
rhoup(1:floor(win/2)) = rho(1:floor(win/2)); 

j = floor(win/2);
for i = 1:length(z)-win+1
    ratio = (1./win).*(ones(1,win));
    Meff(j+1) = (sum(ratio'./M(i:i+win-1))).^-1;
    Geff(j+1) = (sum(ratio'./G(i:i+win-1))).^-1;
    rhoup(j+1) = sum(ratio'.*rho(i:i+win-1));
    vpeff(j+1) = sqrt(Meff(j+1)./rhoup(j+1));
    vseff(j+1) = sqrt(Geff(j+1)./rhoup(j+1));
    j = j+1;
end

Meff(length(z)-win+2: end) = M(length(z)-win+2:end); 
vpeff(length(z)-win+2: end) = vp(length(z)-win+2:end); 
vseff(length(z)-win+2: end) = vs(length(z)-win+2:end); 
rhoup(length(z)-win+2: end) = rho(length(z)-win+2:end); 

end