## Normal Mode Summation for calculating the synthetic seismogram for a string

When we pluck a string fixed at both ends, then this will creating a standing wave. We can get some insights on the behavior of the propagating wave by considering normal modes or free oscillations of the string.

A wave propagating along the string follows the 1D scalar wave equation.

Since, the string has fixed ends, it satisfies the boundary conditions (zero displacement at the fixed ends) and this lead to the vibration of the string at only specific frequencies called eigenfrequencies. The eigenfrequencies are separated by pi*v/L, where v is the velocity of the wave and L is the length of the string. The velocity of the wave depends on the material of the string. Since, the eigenfrequencies are inversely proportional to the length of the string, the longer the string, closer the eigenfrequencies become.

Each eigenfrequencies correspond to a solution of space_term*time_term where the spatial term is known as the spatial eigenfunction.

A traveling wave is the weighted sum of the string’s normal modes. Hence, it is the sum of the eigenfunctions times its amplitude (function of the eigenfrequency). Thus, the real propagating wave is the interference of different modes of the string.

The amplitude for each eigenfunction depends on the position of the source that generated the waves and on the behaviour of the source as a function of time.

We can consider a very simple function to model the behavior of the source as a function of time

F(wn)=exp[-(wn*tau)^2/4]    (Stein and Wysession)

The shape of the source depends on the value of tau.

Let us make a program to calculate a synthetic seismogram using the normal mode summation (see Stein and Wysession (2003) pg. 466). We first keep the parameters same as the Stein and Wysession.

We consider a string of length 1m with a wavespeed 1m/s, a source at 0.2 m and the receiver at 0.7 m. Here, we add 200 modes to get the synthetic seismogram. The seismogram length is 1.25 sec and the time steps is 100.

```program synt_seis_strings
implicit none
real, parameter :: pi = 4*atan(1.0)
real :: stln,vel,xsrc,xrcvr,tseis,dt,tau,npiL,srctm,rcvrtm,wn,Fwn,spctm,coswt
integer :: i,n,j
integer, parameter :: nmode=200    !number of modes (nmode)
integer, parameter :: nt=100    !number of time steps(nt)
real(kind=16), dimension(1:nmode) :: Uxt
real(kind=16), dimension(1:nt) :: Umode
character*12 :: filename

print *, 'Enter the string length (ex: 1.0):'

print *, 'Enter the wave speed: (ex: 1.0):'

print *, 'Enter the position of the source (in m) (ex: 0.2) :'

print *, 'Enter the position of the receiver (in m) (ex: 0.7) :'

print *, 'Enter the duration of the seismogram (in sec) (ex: 1.25) :'

dt=tseis/nt    !time step in sec
tau=0.02    !source shape term

do 5 i=1,nt
Uxt(i)=0.0
5 continue
do 10 n=1,nmode
npiL=n*pi/stln
srctm=sin(npiL*xsrc)    !source term
wn=n*pi*vel/stln    !mode frequency
Fwn=exp(-((tau*wn)**2)/4)    !weighting factor of different frequencies

spctm=srctm*rcvrtm*Fwn        !space term

do 15 j=1,nt
coswt=cos(wn*dt*(j-1))    !time term
Uxt(j)=Uxt(j) + coswt*spctm    !displacement
15 continue

10 continue
open(unit=1,file='seism.dat')
write(1,'(f6.2,f12.6)') (dt*(j-1),Uxt(j), j=1,nt)    !writing the displacement
close(1)

stop
end program synt_seis_strings```

The synthetic seismogram data is saved in the file seism.dat. Now, we run the above program and plot the output.

```#!/bin/bash

prog="synt_seis_strings.f" #fortran program name
exect=`echo \$prog | cut -d"." -f1` #executable name

rm -f *.dat \$exect #remove files

gfortran -ffree-form \$prog -o \$exect #compiling the fortran program
./\$exect #running the executable

datascrpt='modesPlot.gp'
echo "plot [0:1.25][-8:8] 'seism.dat' using 1:2 w lines " >\$datascrpt

#######GNUPLOT##################
gnuplot << EOF
fignm='synth_seis.ps'
figtitle='Synthetic Seismogram for a string'
datafile=system("echo \$datascrpt")

xlb="Time"
ylb="Displacement"
set terminal postscript color solid #terminal type
set title figtitle #figure title
set xlabel xlb #xlabel
set ylabel ylb #ylabel
set output fignm #output fig name
unset key
EOF```

In the above figure, we can see three peaks. First peak at 0.5s is the direct arrival from the source to the receiver. The distance between the source and receiver is 0.5. So the time taken from the source to the receiver is 0.5/1=0.5 sec.

The other two inverted peaks are the reflected arrivals. It is the time taken by the wave from the source to reach the end and reflected back to reach the receiver. The first reflected phase we observe is at 0.9s. It it due to the reflection from the left boundary. This can be calculated ((0.2-0)+0.7)/1=0.9 sec. Similarly, the other arrival (reflected phase from the right boundary) could also be calculated.

We have used the exponential source term (F(wn)=exp[-(wn*tau)^2/4] ) which depends on the value of eigenfrequencies and tau. We have modeled its dependency.

In the second figure, we have considered the value of wn=3.14.

Here, also we have considered the source term to be an exponential function. Instead, if we consider the source term to be abs(sin(x)) then we get the following synthetic seismogram.

Here, we can still see the peaks are 0.5,0.9 and 1.1 seconds but there are some wiggles are it. This is the effect of the source function. We have considered the same number of modes and other parameters.

When we consider the sinc function as the source function then the seismogram wiggles surrounding the peak decreases.

### Conclusion:

In this study, we have obtained the synthetic seismogram using the normal mode summation method. We now have two ways to think of the displacement of a traveling wave, propagating waves or normal modes. These two ways can give different interesting insights about the properties of the source and the medium. The difference between the two is that in the propagating wave formulation, the travel times are measured to infer the velocity of the medium and in the normal mode formulation, eigenfrequecies are measured to infer the velocity.

The normal mode solution gives us information about the medium in which the waves propagate and the source that generates them. The physical properties of the string control its velocity (or eigenfunctions and eigenfrequencies). The displacement obtained can be used to study the source that excited them.

In normal mode solution, we sum all the modes to obtain the synthetic seismogram. The individual modes, though, are not physically meaningful. Each mode actually start vibrating even before the wave reaches the end of the string. But when all the modes are summed, the resulting wave propagate at the correct velocity.

In this study, we have considered a uniform string to calculate the synthetic seismogram. We could generalize this idea to obtain the modes of the non-uniform string. We can consider this case in the future posts. We can even calculate the modes for an sphere.

• Utpal Kumar (IES, Academia Sinica)

## Modeling a wave on a string using finite differences

Based on problem 3.7 chapter 3 of Introduction to Seismology (Shearer)

(COMPUTER) In the case of plane-wave propagation in the x direction within a uniform medium, the homogeneous momentum equation (3.9) for shear waves can be expressed as

,where u is the displacement. Write a computer program that uses finite differences to solve this equation for a bar 100 km in length, assuming β = 4 km/s. Use dx = 1 km for the length spacing and dt = 0.1 s for the time spacing. Assume a source-time function  at u (50 km) of the form:

Apply a stress-free boundary condition at u(0 km) and a fixed boundary condition at u (100 km).  Approximate the second derivatives using the finite difference scheme:

Plot u(x) at 4 s intervals from 1 to 33 s. Verify that the pulses travel at velocities of 4 km/s. What happens to the reflected pulse at each endpoint? What happens when the pulses cross?

```!**********************************************************************
!* *
!* 1-D Plane wave propagation with fixed and free boundaries *
!* *
!**********************************************************************

!----initialzie t, dx, dt, tlen, beta, and u1, u2, u3 arrays----c

implicit none
integer:: i, counter, nx, nx1
integer:: it, ntmax
real(8):: t, dt, dt2, dx, dx2, beta, beta2, tlen, tmax, rhs
real*8, allocatable:: u1(:), u2(:), u3(:)
character(10):: data
character(4):: nf
!----parameters ----c
! dt must be less than (dx/beta) for numerical stability

nx = 100 ! number of grid
dt = 0.10d0 ! time interval
dx = 1.0d0 ! grid interval
beta = 4.0d0 ! wave velocity
tlen = 5.0d0 ! time length of wave source
ntmax = 1000 ! maximum time step (finish at it = ntmax, or)
tmax = 33.0d0 ! maximum calculation time (finish at t< tmax)

!---- allocate dimension variables ----c

nx1 = nx + 1
allocate (u1(nx1), u2(nx1), u3(nx1))

!---- initialize t, counter, u1, u2, u3 ----c

t = 0.0d0
counter = 1

do i = 1, nx1
u1(i) = 0.0d0
u2(i) = 0.0d0
u3(i) = 0.0d0
end do

!---- calculate c^2, dt^2, dx^2 ----c
beta2 = beta**2
dt2 = dt**2
dx2 = dx**2

! ============= Time marching loop ===========

do it = 1, ntmax

t = t + dt

!---- calculate u3(i) ----c

do i= 2, nx
rhs=beta2*2*(u2(i+1)-2.0d0*u2(i)+u2(i-1))/dx2
u3(i)=dt2*rhs+2.0d0*u2(i)-u1(i)
end do

!---- free boundary (du/dx=0) ----c
u3(1)=u3(2)
!---- fixed boundary (u=0) ----c
u3(nx1)=0.0d0
!---- source time function c----c

if (t.le.tlen) then
u3(51) = sin(3.1415927*t/tlen)**2
end if

!---- change new and old variables ----c

do i=1,nx1
u1(i) = u2(i)
u2(i) = u3(i)
end do
!---- make data file ----c
write(nf,'(i3.3)') counter
data = "data"//nf

open(7,file=data,status='replace')
do i=1, nx1
write(7,*) i, u2(i)
end do
!---- output u2 at desired intervals, stop when t is big enough

write(6,'(a5,1x,i3,1x,f6.3)') 'loop', counter, t

if (t.gt.tmax) exit

counter = counter + 1

end do

stop
end```

The following is used to plot the results:

```#!/usr/bin/bash
gnuplot << eof
set term gif animate
set output "animate.gif"
n=9 #n frames
set xrange [0:102]
set yrange [-2:2]
i=0
list = system('ls data*')
do for [file in list] {
plot file w lines lt 1 lw 1.5 t sprintf("t=%i sec",i/10)
i = i + 1
}
set output
eof```

Nguyen Cong Nghia – IESAS

## Calculating Auxiliary Fault Plane Solutions given the main fault (Fortran)

We have calculated the auxiliary fault plane solution using the input of main fault solution in perl. Here, we do the same in Fortran. And we will also plot to solution to visualize the results.

Fortran Code to get the auxiliary fault plane solutions:

```program auxiliary_fault_plane
! Program to calculate the strike, dip and rake of the auxiliary fault plane solutions
! given the strike, dip and rake of the fault plane.
! Authors: Utpal Kumar, Li Zhao
implicit none
real (kind=8):: ph1,del1,lb1,ph2,del2,lb2,sinlb2,coslb2,sinph1_ph2,cosph1_ph2,ph1_ph2,phI,delI,lbI
real, parameter :: pi = 2*asin(1.0)
integer :: n,lbsgn

print *, 'Enter strike, dip and rake respectively of the fault plane (in degrees) (e.g., 40,50,60):'
!ph1=40
!del1=80
!lb1=200
phI=ph1
delI=del1
lbI=lb1
!Converting the strike, dip and rake in radians
ph1=ph1*pi/180  !strike
del1=del1*pi/180        !dip
lb1=lb1*pi/180          !rake

!Adaptation to include the negative values of rake
lbsgn=1
if (lb1 < 0) then
lbsgn=-1
end if

! Calculation of Strike, sip and rake of the auxiliary fault plane
del2=acos(sin(lb1)*sin(del1))   !dip of auxialiary fault plane
sinlb2=cos(del1)/sin(del2)
coslb2=-(sin(del1)*cos(lb1)/sin(del2))
lb2=acos(coslb2)        !rake of auxilairy fault plane
sinph1_ph2=cos(lb1)/sin(del2)
cosph1_ph2=-1/(tan(del1)*tan(del2))
ph1_ph2=acos(cosph1_ph2)

!Checking for the quadrant of the strike angle
if (sinph1_ph2 >= 0 .and. cosph1_ph2 >=0) then
ph1_ph2=ph1_ph2
else if (sinph1_ph2 > 0 .and. cosph1_ph2 < 0) then
ph1_ph2=ph1_ph2
else if (sinph1_ph2 < 0 .and. cosph1_ph2 < 0) then
ph1_ph2=-ph1_ph2
else if (sinph1_ph2 < 0 .and. cosph1_ph2 > 0) then
ph1_ph2=-ph1_ph2
end if

ph2=ph1-ph1_ph2         !strike of auxialiary fault plane

! For dip > 90 degrees and less than 180 degrees
if (del2 > pi/2 .and. del2 < pi) then
ph2= pi + ph2
del2= pi - del2
lb2= 2*pi - lb2
end if

if (lbsgn < 0) then
lb2 = -(2*pi - lb2)
end if

!Adaptation to give the strike value in the range of 0 to 360 degrees
if (ph2 > 2*pi) then
ph2 = ph2 - 2*pi
end if

101 format("The strike of the auxiliary plane is: ",f9.4, " degrees")
102 format("The dip of the auxiliary plane is: ",f9.4, " degrees")
103 format("The rake of the auxiliary plane is: ",f9.4, " degrees")
print 101, ph2*180/pi
print 102,  del2*180/pi
print 103, lb2*180/pi
open(unit=10,file='plt.dat')
10 format("25  25   0 ",  f8.2, f8.2, f8.2, f8.2, f8.2, f8.2," MainFault: ",f7.2,"/",f7.2,"/",f7.2, &    !continuation of line
" AuxFault: ",f7.2,"/",f7.2,"/",f7.2)
write(10,10) phI,delI,lbI,ph2*180/pi, del2*180/pi, lb2*180/pi, phI,delI,lbI,ph2*180/pi, del2*180/pi, lb2*180/pi !writing in the file
end program auxiliary_fault_plane```

Then we save this in the file called “auxiliary_fault_plane.f”. This script basically aim for calculation of the fault plane solution. For visualizing the results, we can execute the program and plot the results using the psmeca command of GMT.

The bash script for compiling and executing the above program and then plotting it is following:

```#!/bin/bash
gfortran -ffree-form auxiliary_fault_plane.f -o auxiliary_fault_plane
./auxiliary_fault_plane
psbasemap -Bwsne -R0/50/0/50 -Jm1.0 -Xc -Yc -P -K >output.ps
echo "24.6 30 20 0 2 MC Focal Mechanism Plot"| pstext -Jm -R -K -O -P -N >>output.ps
psmeca plt.dat -Jm -R -Sa0.40/14/-0.2i -Gbrown -Fa0.2c/cd -Egray -P -O -N -L -V >>output.ps

gs output.ps
ps2pdf output.ps```

The P axis is plotted using circle and T axis using the diamond symbol (-Fa0.2c/cd).

## Travel time curve calculation for spherical, isotropic, homogeneous Earth model

A travel time curve is a graph of the time it takes a seismic ray to reach the station from an earthquake versus the distance of the earthquake to the station (epicentral distance).

Here, we calculate the travel time curve for a spherically homogeneous Earth. This means the seismic velocity is constant for the whole earth. We have considered the three cases,

1. When the seismic ray travel through the body from the surface event to the station (body waves).
2. When the seismic ray for a deeper event (300 km depth) reaches the station.

Fortran code for calculation of the Model:

```program travel_time_calc
!Program to calculate seismic travel times in a spherical, homogeneous Earth.
implicit none
real, parameter :: pi = 2*asin(1.0), R = 6371    ! radius of Earth (km
real :: V, xx, h1, h, dx
integer, parameter :: n = 89
integer :: i
real(kind=16), allocatable, dimension(:) :: x, xrad, Di, Ti, Tim, Dhi, Thi, Thim

V = 12;        ! assumed constant seismic wave velocity (km/s)
dx=2 !increment
allocate(x(0:n),xrad(0:n),Di(0:n),Ti(0:n),Tim(0:n),Dhi(0:n),Thi(0:n),Thim(0:n)) !allocating the size of the arrays
write(*,*) "Distance (in degrees)", "  ","T-time surface (in mins)", "  ", "T-time interior (in mins)"
open(1,file="Tcalc.dat")
xx=0
do i=0,n !do loop
x(i) = xx + dx    !distance
!Path along the body
Di(i) = 2*R*sin(0.5*xrad(i))    ! distance of the surface event along the body in km
Ti(i) = Di(i)/V               ! travel time in seconds
Tim(i) = Ti(i)/60             ! travel time in minutes

!For deeper events (300 km)
h1=300
h=R-h1
Dhi(i) = sqrt(R**2 + h**2 - 2*R*h*cos(xrad(i)))    !distance of the interior events in km
Thi(i) = Dhi(i)/V                ! travel time in seconds
Thim(i)= Thi(i)/60                ! travel time in minutes

write(1,10) x(i), Tim(i), Thim(i) !printing output on the terminal
write(*,10) x(i), Tim(i), Thim(i) !printing output on the terminal
xx = x(i)
end do
10 format(f7.2, "  ",f8.4, "  ",f8.4)
deallocate(x, xrad, Di, Ti, Tim, Dhi, Thi, Thim) ! deallocates the array
end program travel_time_calc```

For the above program, we have considered the Earth to be isotropic and homogeneous with constant velocity of 12 km/s. Then we parameterize for 89 events each with 2 degrees interval from 2 to 180 degrees. The radius of the earth is assumed to be 6371 km.

For compiling and plotting (bash script):

```#!/bin/bash

gfortran -ffree-form travel_time_calc.f -o travel_time_calc
./travel_time_calc

gnuplot << EOF
set terminal postscript color solid #terminal type
set title 'Travel Time Curve' #figure title
set xlabel 'Delta (in degrees)' #xlabel
set ylabel 'Travel Time (in mins)' #ylabel
set output 'travel_time.ps' #output fig name
set grid #grid
plot [:] [:] "Tcalc.dat" using 1:2 title "Calculated body waves travel time (surface event)" with lines lc rgb "blue", \
"Tcalc.dat" using 1:3 title "Calculated body wave travel time (300 km deep event)" with lines lc rgb "red"
EOF
```

– Utpal Kumar (IES, Academia Sinica)

## Simple Wave Plot (Fortran and Gnuplot demo)

We can use the power of Fortran to do computation of large data set and then we can use Gnuplot to visualize the results. Here, we take a simple example to see how we can do that.  Let us take the first example, which we did using the MATLAB. We can do the similar calculation using the Fortran and Gnuplot and Bash.

```! Program to calculate the harmonic wave function
program harmonic_wave
implicit none ! do not assume the data type itself
integer :: i, n !defining integer data type
real :: dx, tm !defining real data type
real, parameter :: pi = 4*atan(1.0) !defining a parameter with fixed value
real(kind=16), allocatable, dimension(:) :: x, y !defining real array data type with 16 bit size, and allocatable dimension
real(kind=16), dimension(0:1) :: f, w, TP, lb, k, a !defining real array of size 2

n=1000 !number of data points
allocate(x(0:n), y(0:n)) !allocating the size of the arrays
a(0)=2 !amplitude1
a(1)=3 !amplitude1
f(0)=5 !frequency1
f(1)=10 !frequency2
TP(0:1)=1/f(0:1) !time period
w(0:1)=2*pi*f(0:1) !angular frequency
lb(0:1)=2*TP(0:1) !wavelength
k(0:1)=1/lb(0:1) !wavenumber
tm=2 !time

dx=pi/200 !increment
x(0:n) =[(i*dx, i=0,n)] !x values
y(0:n)=a(0)*sin(k(0)*x(0:n)-w(0)*tm) + a(1)*cos(k(1)*x(0:n)-w(1)*tm); !waveform
open(unit=1, file='wavedata.dat') !opening a file for writing data
do i=0,n !do loop
!print 11, i, x(i), y(i) !printing output on the terminal
write(1,10) x(i), y(i) !writing in the file
10 format(f8.4, f8.4) !format for writing in the file
11 format(i4, " x: ", f8.4, " y: ", f8.4) !format for output on the terminal
end do
close(1) !close the file
deallocate(x, y) ! deallocates the array
end program harmonic_wave```

Let us save this in the file called “simple_wave_model.f”.

For compiling, we can use gfortran and use free-form source. If the file extension is .f90, then gfortran by default consider it to be free form source, otherwise it consider it to be a fixed form. We can overwrite this by giving extra argument of “-ffree-form”.

In the command line, we typed first command to compile the program and the next one to run it. We can write these sets of program in bash script and make it more general. We can also add the commands for plotting the results using gnuplot in the bash script.

```#!/bin/bash
prog="simple_wave_model.f" #fortran program name
exect=`echo \$prog | cut -d"." -f1` #executable name

rm -f *.ps *.png *.eps \$exect #remove files

gfortran -ffree-form \$prog -o \$exect #compiling the fortran program
./\$exect #running the executable

###GNUPLOT parameters
fignm="waveplot.ps" #figure name for the plot
figtitle="Wave Plot" #figure title
datafile="wavedata.dat" #data file name

##Deciding output figure format (ps or png) NOTE: there are many more formats available
figtype=`echo \$fignm | cut -d"." -f2`
if [ "\$figtype" == "ps" ]; then
figformat=postscript
elif [ "\$figtype" == "png" ]; then
figformat=png
fi

#######GNUPLOT
gnuplot << EOF
fignm=system("echo \$fignm")
figtitle=system("echo \$figtitle")
datafile=system("echo \$datafile")
figformat=system("echo \$figformat")

subt1='Harmonic wave'
xlb="x -values"
ylb="Amplitude"
set terminal postscript color solid #terminal type
set title figtitle #figure title
set xlabel xlb #xlabel
set ylabel ylb #ylabel
set output fignm #output fig name
set grid #grid
plot [:] [:] datafile using 1:2 title subt1 with lines lc rgb 'blue' #plotting the figure
EOF

open \$fignm #opening the figure for Mac user
#gs \$fignm #opening figure using ghostviewer

```

Save the above file as “wavePlot.sh” and change the permission using the “chmod” command.