### plotting the results of rayleigh.m
#################################################################################

clear;
load 'rayleigh_result.mat' model_* z vp vs nfreq freq kreal vphase vgroup;
load 'data.mat' phvel_data_*

nmodes = size(freq)(1);

figure(1);
clf;

hold on;

for i=1:nmodes

	if ( i <= 6 )
		plot(freq(i,:),vphase(i,:),sprintf('%d;mode %d;',i,i));
	elseif ( i <= 12 )
		plot(freq(i,:),vphase(i,:),sprintf('+%d;mode %d;',i-6,i));
	endif

endfor

plot(phvel_data_1(:,1),phvel_data_1(:,2),';obs 1;');	
xlabel('Frequency');
ylabel('Phase velocity');

axis([0 80 0 0.3]);