kennettQ2_tf.m

function [tf]=kennettQ2_tf(lyr,om,mopt,fs,Q,omega0)
% KENNET CAUSAL algorithm with Q for Synthetic seismograms for plane wave, normal incidence propagation
% in 1-D layered media. This algorithm uses a nearly constant Q model proposed by Kennett, with
% a velocity correction given by Aki and Richards.
% When Q becomes very large, the results should converge with the usual elastic Kennett algorithm.
%
%[TF]=kennettQ2_tf(lyr,om,mopt,fs,Q,omega0)
%
%	TF:    [frequency, reflectivity, transmissivity]
%	       Transfer function (complex) of the layered medium
%	LYR:   [velocity, density, thickness] of layers
%          velocities should correspond to frequency equal to infinite
%	       LYR is a matrix of 3 columns with
%             number of rows = number of layers
%   om:     Row vector of all omega (2*pi*frequency) values for which the
%   transfer function is calculated
%   MOPT: =0 for primaries only, =1 for 1st order multiples
%             anything else calculates all reverberations
%	FS:    option for free surface reverberations.
%             FS=0 for no free surface multiples, FS=1 for free surface effect
%  Q:     Quality factor (the inverse of the attenuation): it is a measure
%         of how dissipative the material is. Typical values are: 50,100, 200...
% OMEGA0: Reference omega (2*pi*frequency) for Q
%

% Adapted from Mukerji's Kennet, by Manuel Diaz to add finite Q
% May 17, 2001
%
% Increased vectorization for faster computation by Gary Mavko, Nov 2001
% Also increased the number of comments.
% Modified by Vishal Das to just output transfer functions
% June 2017

i=sqrt(-1);

[nlr,nc] = size(lyr);              % dimensions of input Earth model
ro = lyr(:,2);                     % column of layer densities
v  = lyr(:,1);                     % column of layer velocities
[v1,v2] = size(v);
d = lyr(:,3);                      % column of layer thicknesses
[om1,om2] = size(om);

v  = repmat(v,1,om2);              % repeat velocity column by the number of frequencies
om = repmat(om,v1,1);              % repeat frequency row by the number of layers
om(:,1)=0.000001;
ro = repmat(ro,1,om2);             % repeat density column by the number of frequencies

if length(Q)==1
    Q=Q.*ones(v1,1);               % make column of Q, in case only a scalar is input
end
Q=repmat(Q,1,om2);                 % repeat Q column by the number of frequencies

% Kennett nearly constant Q dispersion model
Qw=Q.*(ones(v1,om2)-((pi.*Q).^(-1)).*log(0.577.*om./omega0));
Qwliu=Qw;

Qw(isinf(Qw)) = mean(Qw(~isinf(Qw)));

% Modified velocity, from Aki and Richards equation 5.94
v=v.*(ones(v1,om2) +log(om./omega0)./(pi*Qw) -i./(2*Qw));

v(isinf(v)) = mean(v(~isinf(v)));

% % Calculate the Ray theory velocity
% f = d./sum(d);
% f = repmat(f,1,om2);
% vrt = (sum(f(1:end,:)./v(1:end,:))).^-1;
% [index,index] = min(abs(om(1,:)-omega0));
% vrt = vrt(index);

freq  = om(1,:)/(2*pi);

rdhat = zeros(size(om(1,:)));
tdhat = ones(size(om(1,:)));

deno = ones(size(v)); rd = ones(size(v)); td = ones(size(v));
deno(2:nlr,:) = ro(2:nlr,:).*(v(2:nlr,:)) + ro(1:nlr-1,:).*(v(1:nlr-1,:));
rd(2:nlr,:)   = (ro(1:nlr-1,:).*(v(1:nlr-1,:)) - ro(2:nlr,:).*(v(2:nlr,:)))./deno(2:nlr,:);
td(2:nlr,:)   = 2*sqrt(ro(2:nlr,:).*(v(2:nlr,:)).*ro(1:nlr-1,:).*(v(1:nlr-1,:)))./deno(2:nlr,:);

if fs==1,
    rd(1,:) = -1*rd(1,:);
    td(1,:) =  1*td(1,:);
else
    rd(1,:) = 0*rd(1,:);
    td(1,:) = 1*td(1,:);
end

ru=-rd; tu=td;

for j=nlr:-1:1
    ed=exp(i.*(d(j)./v(j,:)).*om(j,:));
    if mopt==0
        reverb=ones(size(om(j,:)));
    elseif mopt==1
        reverb=ones(size(om(j,:)))+ru(j,:).*ed.*rdhat.*ed;
    else
        reverb=1./(1-ru(j,:).*ed.*rdhat.*ed);
    end
    rdhat=rd(j,:)+tu(j,:).*ed.*rdhat.*ed.*reverb.*td(j,:);
    tdhat=tdhat.*ed.*reverb.*td(j,:);
    if mopt==1
        dx=find(real(rdhat)>1); rdhat(dx)=ones(size(dx))+i*imag(rdhat(dx));
        dx=find(imag(rdhat)>1); rdhat(dx)=real(rdhat(dx))+i*ones(size(dx));
        dx=find(real(tdhat)>1); tdhat(dx)=ones(size(dx))+i*imag(tdhat(dx));
        dx=find(imag(tdhat)>1); tdhat(dx)=real(tdhat(dx))+i*ones(size(dx));
    end
end;


tf=[freq(:), rdhat(:), tdhat(:)];


%dont bother plotting
return;