SsPmP/plot_sac.m at master · jinwar/SsPmP · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
% script to plot sspmp
% written by Ge Jin

clear;

beforetime = 10;
aftertime = 20;
filt = [0.02 0.5];
Vp = 6.5;

event = '20102050535'

sacfiles = dir('data/',event,'/*.sac.z');
for i = 1:length(sacfiles)
	sac(i) = readsac(sacfiles(i).name);
end

stlas = [sac.STLA];
stlos = [sac.STLO];
avgdist = distance(mean(stlas),mean(stlos),sac(1).EVLA,sac(1).EVLO);
taup_com = ['taup_time -ph S -rayp -deg ',num2str(avgdist),' -h ',num2str(sac(1).EVDP),'> taup_temp'];
system(taup_com);
rayp = load('taup_temp');
rayp = rayp/deg2km(1);

sort_key = [sac.STLA];
ids = 1:length(sac);

mat = [ids(:),sort_key(:)];

mat = sortrows(mat,2);

figure(43)
clf
hold on
offset = 0;
for i=1:size(mat,1)
	id = mat(i,1);
	if ~isnan(sac(id).T2)
		taxis = sac(id).B:sac(id).DELTA:sac(id).B+sac(id).DELTA*(sac(id).NPTS-1);
		ind = find(taxis > sac(id).T2 - beforetime & taxis < sac(id).T2 + aftertime);
		data = sac(id).DATA1(ind);
		fN = 1/sac(id).DELTA/2;
		[b,a] = butter(2,[filt(1)/fN filt(2)/fN]);
		data = filtfilt(b,a,data);
		data = data./max(abs(data));
%		data = data.*(max(mat(:,2))-min(mat(:,2)))/20;
		t = -beforetime:sac(id).DELTA:aftertime;
		syndt = sac(id).T1 - sac(id).T2;
		if length(t) > length(data)
			t = t(1:length(data));
		end
		h = t./2./(Vp^(-2) - rayp^2);
%		offset = mat(i,2);
		offset = offset+1;
		plot(t,data + offset);
		plot(syndt,offset,'rx','markersize',15);
		data(find(data<0)) = 0;
		area(t,data + offset,offset);
		text(t(1)-3,offset,sac(id).KSTNM);
	end
end

figure(44)
clf
hold on
offset = 0;
for i=1:size(mat,1)
	id = mat(i,1);
	if ~isnan(sac(id).T2)
		taxis = sac(id).B:sac(id).DELTA:sac(id).B+sac(id).DELTA*(sac(id).NPTS-1);
		ind = find(taxis > sac(id).T2 - beforetime & taxis < sac(id).T2 + aftertime);
		data = sac(id).DATA1(ind);
		fN = 1/sac(id).DELTA/2;
		[b,a] = butter(2,[filt(1)/fN filt(2)/fN]);
		data = filtfilt(b,a,data);
		data = data./max(abs(data));
%		data = data.*(max(mat(:,2))-min(mat(:,2)))/20;
		t = -beforetime:sac(id).DELTA:aftertime;
		syndt = sac(id).T1 - sac(id).T2;
		if length(t) > length(data)
			t = t(1:length(data));
		end
		h = t./2./((Vp^(-2) - rayp^2).^.5);
%		offset = mat(i,2);
		offset = offset+1;
		plot(h,data + offset);
		data(find(data<0)) = 0;
		area(h,data + offset,offset);
		text(-25,offset,sac(id).KSTNM);
	end
end
xlim([-10 100])