# Function returning the phase and group velocities for a Rayleigh wave 
# in the globally-defined velocity/density model.
#	'kreal' is the wavenumber for which the eigenfrequency and velocity are
#	determined
######################################################################

function [ vphase vgroup ] = velrayleigh( kreal )

global n2 rho lambda myu l2myu

I1 = zeros(n2,n2);	# This will contain the kinetic-energy matrix
I2 = zeros(n2,n2);	# Bulk-energy matrix
I3 = zeros(n2,n2);	# Shear-energy matrix

#z
#kreal

# Evaluate the energy matrices for pairs of basis functions
# Basis functions are constructed so that they satisfy the boundary conditions
# on the surface, internal interfaces, and at teh bottom of the model.

for i=1:n2

  [ ir1 r1i ir2 r2i il1 il2 ] = rayleigh_basis(i,kreal);

#ir1
#r1i

  jstart = i-8;

  if ( jstart < 1 )
	jstart = 1;
  endif

  jend = i+8;

  if ( jend > n2 )
	jend = n2;
  endif

  for j=jstart:jend

  	[ jr1 r1j jr2 r2j jl1 jl2 ] = rayleigh_basis(j,kreal);

	if ( jl1 < il1 )
		jl1 = il1;
	endif

	if ( jl2 > il2 )
		jl2 = il2;
	endif

#jl1
#jl2

  	for layer = jl1:jl2

  		r1_1 = r1i*bfunpoly(layer,ir1,0);
  		r1d_1 = r1i*bfunpoly(layer,ir1,1);
  		r2_1 = r2i*bfunpoly(layer,ir2,0);
  		r2d_1 = r2i*bfunpoly(layer,ir2,1);

  		r1_2 = r1j*bfunpoly(layer,jr1,0);
  		r1d_2 = r1j*bfunpoly(layer,jr1,1);
  		r2_2 = r2j*bfunpoly(layer,jr2,0);
  		r2d_2 = r2j*bfunpoly(layer,jr2,1);

		Ekin = 0.5 * rho(layer) * \
				( 	integral(layer,conv(r1_1,r1_2)) +
					integral(layer,conv(r2_1,r2_2)) 
				);

		p1 = polyadd(kreal*r1_1,r2d_1);
		p2 = polyadd(kreal*r1_2,r2d_2);

		Elambda = 0.5*lambda(layer) * integral(layer,conv(p1,p2));

		p1 = polyadd(r1d_1,-kreal*r2_1);
		p2 = polyadd(r1d_2,-kreal*r2_2);

		Emyu = 0.5*myu(layer) * \
			(
				integral(layer,conv(p1,p2)) + 
				2.0*kreal*kreal*integral(layer,conv(r1_1,r1_2)) +
				2.0*integral(layer,conv(r2d_1,r2d_2))
			);

#Ekin 
#Elambda 
#Emyu

		I1(i,j) += Ekin;
		I2(i,j) += Elambda;
		I3(i,j) += Emyu;
	endfor

  endfor

endfor

#I1

A = inv(I1)*(I2+I3);

	# sort eigenvectors in descending order

[evect1,evalarr] = eig(A);
eval1 = diag(evalarr);
[eval,indeval] = sort(eval1,'ascend');
evect=evect1(:,indeval);

# test :
#eval
#evect
#wl = 2*pi/kreal		# wavelength
# l1 = evect(:,2);
# A*l1 ./ l1 / eval(2)	 	# test: these values should all equal 1
	
	# 'evel' now contains squared eigenfrequencies.
	# Determine the phase velocities from them:

vphase = sqrt(eval)/kreal;

	# determine group velocities:

for i=1:n2
	l1 = evect(:,i);
	vgroup(i) = l1' * (I2 + I3/2.0/kreal) * l1 / (l1' * I1 * l1) / vphase(i);
endfor

endfunction