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

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)

GNUPLOT (making animation)

GNUPLOT is very powerful set of programs for making graphs. Here, we see how can we use it for making animation and saving it in gif format.

#!/bin/bash
rm -f gnuplot.rot
###Making input file for the gnuplot
for ((z=0; z<=360; z=z+5)) #for varying the longitude
do
cat >> gnuplot.rot <<!
zrot=$z
zview(zrot)=zrot
set view xview(xrot), zview(zrot), 2, 1
set size square
set view xview(xrot), zview(zrot), 2, 1
set size square
splot cos(u)*cos(v),cos(u)*sin(v),sin(u) notitle with lines lt 5, \
      'world10.dat' notitle with lines lt 3 lw 1     
!
done

echo "Now Plotting"


gnuplot<<!
set term gif transparent nocrop enhanced animate delay 20 loop 0 nooptimize #setting the terminal size 600,600 background rgb 'white' font "verdana,12"
set output "globe.gif" #setting the output name
unset title #notitle
unset key #no key
unset xtics #no x axis ticks
unset ytics #no y axis ticks
set border 0 #no border
set hidden3d nooffset #hidden line removal for surface plotting
set parametric
set angles degrees
set samples 128,128 #no of points
set isosamples 13,13
set mapping spherical
set dummy u,v
set urange [ -90.0000 : 90.0000 ] noreverse nowriteback
set vrange [ 0.00000 : 360.000 ] noreverse nowriteback
set style data line
xrot=90
xview(xrot)=xrot
load "gnuplot.rot"
!

globe.gif

Data File: world10.dat