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):'
    read *, stln

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

    print *, 'Enter the position of the source (in m) (ex: 0.2) :'
    read *, xsrc
    
    print *, 'Enter the position of the receiver (in m) (ex: 0.7) :'
    read *, xrcvr

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

    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
        rcvrtm=sin(npiL*xrcvr)    !Receiver 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
load datafile
EOF

test1

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.

fwn_wnfwn_tauwn3-14

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.test1.png

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.

test1.png

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.

test1.png

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

screenshot-from-2016-12-14-15-42-21

,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:

screenshot-from-2016-12-14-15-42-16

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:

screenshot-from-2016-12-14-15-42-25

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?

Solutions can be downloaded here.

!**********************************************************************
!* *
!* 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

animate

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):'
        read *, ph1,del1,lb1
        !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).

Screenshot from 2016-12-05 14-34-26.png

output1

-Utpal Kumar (IES, Academia Sinica)

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
        xrad(i) = x(i)*pi/180 ! convert degrees to radians
        !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

travel_time

– 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”.

Screen Shot 2016-11-14 at 8.23.09 PM.png

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.

Screen Shot 2016-11-14 at 8.30.17 PM.png

Screen Shot 2016-11-14 at 8.32.52 PM.png

-Utpal Kumar, IES, Academia Sinica

Guessing a number (Fortran)

This program will randomly pick a value and the user is prompted to enter their guess. It will keep counting your number of attempts to reach the value stored by the system. The lower the number of guesses you have, the better guesser you are or your algorithm is better.

To run the program:

Screen Shot 2016-10-27 at 4.54.33 PM.png

You can play with the program and try to modify it to learn more. Let’s see how many attempts you need!

Program: guess.f (by Cedric Legendre IES, Academia Sinica)