## How to make the rotating globe

Many people have asked me how to make the rotating globe we have on our home page. This post will take you step by step how to make the rotating globe gif using “gnuplot”.

## Installing GNUPLOT

In Mac, you can install gnuplot using the homebrew:

``````ruby -e "\$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"
``````

This will install homebrew in your Mac in case you don’t have.
Then use the command:

``````brew install gnuplot --with-x11
``````

This should do the job and then you should be able to run the command
`gnuplot`

If this throws error then you can overwrite the symbolic links by
`brew link --overwrite gnuplot`

## Making the rotating globe gif

Now, you can download the codes for making the gif from the GitHub by clicking here (Data Source: gnuplotting). Run this bash script in your computer, it may take a few seconds and you should be able to get the “globe.gif” file.

## 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

## 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)

## 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