########## Modeling Rayleigh-wave velocities in a layered model
##################################################################################

clear;

# several (optional) versions of the model shown by different colors in the plots.
# Columns of these matrices are:
#	1) depth to the bottom of layer
#	2) density
#	3) P-wave velocity
#	4) S-wave velocity
#	4) Qs

global model_blue = [
	0.0005	15/9.8	0.19	0.095	1000		
	0.002	15/9.8	0.22	0.11	1000			
	0.005	11/9.8	0.23	0.115	1000	
	0.008	11/9.8	0.24	0.12	1000	
	0.015	11/9.8	0.3	0.15	1000	
	0.035	15/9.8	0.3	0.15	1000
];

global model_red = [
	0.0005	15/9.8	0.18	0.09	1000		
	0.002	15/9.8	0.21	0.105	1000			
	0.005	11/9.8	0.24	0.12	1000	
	0.008	11/9.8	0.27	0.135	1000	
	0.015	11/9.8	0.3	0.15	1000	
	0.035	15/9.8	0.3	0.15	1000
];

global model_yellow = [
	0.0005	15/9.8	0.188	0.094	1000		
	0.002	15/9.8	0.2	0.1	1000			
	0.005	11/9.8	0.23	0.115	1000	
	0.008	11/9.8	0.238	0.119	1000	
	0.015	11/9.8	0.3	0.15	1000	
	0.035	15/9.8	0.3	0.15	1000
];

global model_green = [
	0.0005	15/9.8	0.192	0.096	1000		
	0.002	15/9.8	0.222	0.111	1000			
	0.005	11/9.8	0.232	0.116	1000	
	0.008	11/9.8	0.242	0.121	1000	
	0.015	11/9.8	0.3	0.15	1000	
	0.035	15/9.8	0.3	0.15	1000
];

# the model actually used in forward modeling of Rayleigh waves:

global model = model_green;

global nlayers = size(model)(1);
global nbfun = 2+(nlayers-1)*2;   # two basis functions for r1 or r2 on each boundary except bottom, 
				# plus two basis functions on the free surface 
global n2 = 2*nbfun;		# number of parameters in r1 and r2 combined

global rho = model(:,2);	# density
global vp = model(:,3);		# P-wave velocity
global vs = model(:,4);		# S-wave velocity

global myu = vs.^2 .* rho;
global lambda = vp.^2 .* rho - 2.0 * myu;
global bulk = lambda - 2.0/3.0 * myu;
global l2myu = lambda + 2.0*myu;

lambda2 = [ 0.0005 : 0.00025 : 0.02 ];	# desired half-wavelengths for output
kreal = pi*lambda2.^-1;			# real part of the wavenumber for each lambda2(i)

nfreq = length(kreal);

freq = zeros(n2,nfreq);			# output frequencies, first index is mode number
vphase = zeros(n2,nfreq);		# phase velocities
vgroup = zeros(n2,nfreq);		# group velocities

# scan the range of frequencies and solve the eigenvalue probem for each:

for i=1:nfreq
	[ vphase(:,i) vgroup(:,i) ] = velrayleigh(kreal(i) );
	freq(:,i) = kreal(i)*vphase(:,i)/2/pi;
endfor

## save all results in a .mat file 

save 'rayleigh_result.mat' model_* vp vs nfreq freq kreal vphase vgroup