function [ur,uy] = green(freq,vr,group,energy,z,r,offsets,Fx,Fy,Fz,Phi,s_depth,r_depth)
% This function calculates the Green's function of the horizontal and vertical
% Rayleigh wave displacement field using modal superposition.
% Copyright 1999 by Glenn J. Rix and Carlo G. Lai
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
% Initialize matrices for the horizontal and vertical displacements
ur = zeros(length(freq),length(offsets));
uy = zeros(length(freq),length(offsets));
% Define the vector of circular frequencies
om = 2*pi*freq;
% Loop through the frequencies
for j = 1:length(freq)
% Calculate the position of the source and receiver
index = find(z(:,j) >= s_depth);
s_index = index(1);
index = find(z(:,j) >= r_depth);
r_index = index(1);
% Determine the number of modes at each frequency
index = find(vr(j,:));
% Calculate the wavenumbers
k = om(j)./shiftdim(vr(j,index),1);
% Assign displacement vectors, velocities, and energy integrals to local variables
r1s = shiftdim(r(j,index,s_index,1),1);
r1r = shiftdim(r(j,index,r_index,1),1);
r2s = shiftdim(r(j,index,s_index,2),1);
r2r = shiftdim(r(j,index,r_index,2),1);
V = shiftdim(vr(j,index),1);
U = shiftdim(group(j,index),1);
I1 = shiftdim(energy(j,index),1);
% Loop through the offsets
for n = 1:length(offsets)
a = (Fz*r2s + i*(Fx*cos(Phi)+Fy*sin(Phi))*r1s)./(8*V.*U.*I1)./...
sqrt((pi*k*offsets(n))/2).*exp(-i*k*offsets(n)).*r1r*exp(-i*pi/4);
b = (Fz*r2s + i*(Fx*cos(Phi)+Fy*sin(Phi))*r1s)./(8*V.*U.*I1)./...
sqrt((pi*k*offsets(n))/2).*exp(-i*k*offsets(n)).*r2r*exp(i*pi/4);
ur(j,n) = sum(a);
uy(j,n) = sum(b);
end
end
