Seismogram Data Generation Program

Seismogram Data Generation Program

Updated 08/12/04

The following program was written to generate the seismogram data files for each site. At this point, the sites are entered one at a time, files are generated and printed. These files are renamed by the program operator before data for the next site is generated. In a future program, the functions of the Distance Calculation program and this will be combined into one program and all files will be printed out with separate names.

This program uses two sine waves added together to generate the earhquake waves. In addition, an exponential function is attached to the first wave so that there is a decay of the signal with time. The p wave is generated from time t=1 up to t=1+the sp time difference. The s wave is generated from that time forward. A one second random signal is added to the front of the signal.

Advanced students could be invited to work with the various frequencies in each sine function, as well as the decay constant, to see the effect on the resulting wave. Optimally, a student would discover a "better" combination of factors than those listed here to give a "truer" earthquake wave. 


% Program to calculate distance between epicenter of quake and
% each participating MVHS site.  Create seismic data and print
% to file

clear

load info.dat
z = input('Number of the site to calculate ');
f = input('Frequency for signal recording ');

amp = info(z,3);
nsamp = info(z,4);
ewamp = info(z,5);
tdiff = info(z,6);
ttot = info(z,7);
p=1/f;

% Time for p wave
tp=1+p:p:tdiff+1+p;

% Time for s wave
ts=tdiff+1+2*p:p:ttot;

% Time for total signal
txp=0:p:(ttot-tdiff-(1+2*p));

% Calculate p wave for z
% x = main signal, r = random signal
% 0.8 and 0.2 to determine amplitude of these two
% exp function for decay so signal dies off
% sine is main function - values between pi and tp are frequencies

xpz=0.8*exp(-.05*tp) .*sin(2*pi*2.0*tp);
rpz=0.2*sin(2*pi*2.87*tp);
ypz=(amp/4)*(xpz+rpz)';

% Calculate p wave for ns
xpns=0.8*exp(-.05*tp) .*sin(2*pi*2.0*tp);
rpns=0.2*sin(2*pi*2.85*tp);
ypns=(amp/(5-nsamp))*(xpns+rpns)';

% Calculate p wave for ew
xpew=0.8*exp(-.05*tp) .*sin(2*pi*2.0*tp);
rpew=0.2*sin(2*pi*2.89*tp);
ypew=(amp/(5-ewamp))*(xpew+rpew)';


% Calculate s wave for z
xsz=0.7*exp(-.05*txp) .*sin(2*pi*1*ts);
rsz=0.3*sin(2*pi*1.54*ts);
ysz=(amp)*(xsz+rsz)';

% Calculate s wave for ns
xsns=0.7*exp(-.05*txp) .*sin(2*pi*1*ts);
rsns=0.3*sin(2*pi*1.52*ts);
ysns=(amp)*(xsns+rsns)';

% Calculate s wave for ew
xsew=0.7*exp(-.05*txp) .*sin(2*pi*1*ts);
rsew=0.3*sin(2*pi*1.56*ts);
ysew=(amp)*(xsew+rsew)';

tn=0.0:p:1.0;
yn=(amp/50)*randn(size(tn));
t=0:p:ttot;

yew(1:(f+1))=yn;
yew((f+2):(length(tp)+(f+1)))=ypew;
yew((length(tp)+(f+2)):length(t))=ysew;

yns(1:(f+1))=yn;
yns((f+2):(length(tp)+(f+1)))=ypns;
yns((length(tp)+(f+2)):length(t))=ysns;

yz(1:(f+1))=yn;
yz((f+2):(length(tp)+(f+1)))=ypz;
yz((length(tp)+(f+2)):length(t))=ysz;


fid=fopen('SeisV.dat','w');
fprintf(fid,'%3.4f\n',yz);
fclose(fid);

fid=fopen('SeisNS.dat','w');
fprintf(fid,'%3.4f\n',yns);
fclose(fid);

fid=fopen('SeisEW.dat','w');
fprintf(fid,'%3.4f\n',yew);
fclose(fid);