clear 
close all

% this program reads in data from the Wetzell observatory of the 26-29
% december 2004 with a sampling interval of 1s
% the goal is to see whether eigenmodes of the Earth can be identified
% the data represent transverse velocity

% load day traces
load wet26.dat
load wet27.dat
load wet28.dat
load wet29.dat

% concatenate to one long vector
s=[wet26; wet27; wet28; wet29];
t=1:length(s);
t=t/3600;  % t is now in hrs


% Plot seismogram
subplot(211),plot(t,s,'k-')
xlabel(' Time (s) ')
ylabel(' Amplitude ')


% calculate with spectrum with program myspectrum.m
[fs,f]=myspectrum(s,1);

% plot spectrum
subplot(212)
ff=abs(fs(1:length(f)));
plot(f*1000,ff)
%axis([0.1 6.5 0 Inf]),hold on
xlabel(' Frequency (mHz) ')
ylabel(' Spectral amplitude ')