## loading picked phase velocity values 
#########################################

clear;

# phase-velocity curves for different locations on the line:

phvel_data_1 = [
4.05806,0.12917
4.66232,0.123003
5.26659,0.117836
6.87796,0.112168
9.49645,0.107
11.7121,0.102666
15.5391,0.0988306
21.5818,0.097161
26.2145,0.0956587
31.25,0.0948229
36.8898,0.0944868
42.1268,0.0964843
48.7737,0.098481
55.4206,0.100144
61.0604,0.101475
];

phvel_data_2 = [
8.18935,0.110719
10.2012,0.107716
11.4083,0.104714
13.5207,0.10288
15.8343,0.100128
17.7456,0.0977932
20.0592,0.095542
22.8757,0.0940418
26.3964,0.092959
30.4201,0.0925439
34.8462,0.0926295
40.1775,0.0924653
44.6036,0.0937188
49.7337,0.0942218
53.858,0.0947243
57.5799,0.0949764
61.9053,0.0943946
65.6272,0.0940628
68.8462,0.094398
71.4615,0.0938988
74.8817,0.0928994
115.923,0.0779039
];

phvel_data_3 = [
2.56281,0.123674
3.76884,0.119673
5.9799,0.115504
7.18593,0.111503
8.59296,0.108835
10.603,0.103833
12.4121,0.101164
15.0251,0.0978281
19.0452,0.0954906
23.8693,0.0946522
28.6935,0.0938139
34.3216,0.092808
39.7487,0.0926357
45.1759,0.09213
];

# add columns for wavenumbers:

phvel_data_1(:,3) = phvel_data_1(:,1) * 2*pi ./ phvel_data_1(:,2);
phvel_data_2(:,3) = phvel_data_2(:,1) * 2*pi ./ phvel_data_2(:,2);
phvel_data_3(:,3) = phvel_data_3(:,1) * 2*pi ./ phvel_data_3(:,2);

save 'data.mat' phvel_data_*

figure(1);
clf;

hold on;

plot(phvel_data_1(:,1),phvel_data_1(:,2),'b;1;');
plot(phvel_data_2(:,1),phvel_data_2(:,2),'r;2;');
plot(phvel_data_3(:,1),phvel_data_3(:,2),'g;3;');

xlabel('Frequency');
ylabel('Vphase');

figure(2);
clf;

hold on;

plot(pi*phvel_data_1(:,3) .^-1,phvel_data_1(:,2),'b;1;');	# argument = 0.5*wavelength (apparent depth)
plot(pi*phvel_data_2(:,3) .^-1,phvel_data_2(:,2),'r;2;');
plot(pi*phvel_data_3(:,3) .^-1,phvel_data_3(:,2),'g;3;');

xlabel('Apparent depth');
ylabel('Vphase');