bkusrunavgvisco.m

function [invqeff,veff,rhoup]=bkusrunavgvisco(z,vp,invqp,rho,win,om,omega0)
%        [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);

% Calculating the effective Q for each block
invqeff = zeros(1,length(z));
veff = zeros(1,length(z));
rhoup = zeros(1,length(z));
invqeff(1:floor(win/2)) = invqp(1:floor(win/2)); 
veff(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));
    layer = [vp(i:i+win-1) rho(i:i+win-1) dz(i:i+win-1)]; 
    invqeff(j+1) = qeffemt(layer, om./(2*pi), 1./invqp(i:i+win-1),omega0);
    veff(j+1) = veffemt(layer, om./(2*pi), 1./invqp(i:i+win-1), omega0 );
    rhoup(j+1) = sum(ratio'.*rho(i:i+win-1));
    j = j+1;
end

invqeff(length(z)-win+2: end) = invqp(length(z)-win+2:end); 
veff(length(z)-win+2: end) = vp(length(z)-win+2:end); 
rhoup(length(z)-win+2: end) = rho(length(z)-win+2:end); 

end