## Introduction

Often in seismology, we require the information about the earthquake such as location, origin time, depth, magnitude, etc. To obtain this information, users generally, browse the web or if he/she is experienced then just go to the predetermined sites such as global CMT catalog. Even, after reaching the correct website, we need to do some dig up to get the desired information. This program is meant to scrape the global CMT website to obtain the earthquake information.

## How to use the program

The package you download from the given site consists of two python programs. These programs require pre-installed python-3.

You can obtain python from the following websites:
1. Python software

After installing the python 3, to run the program, type:
`python3 EQinfo_extractor.py`
Then enter the event recognition parameter such as
`1976,3,8,4,39`
And enter. It will prompt the event information.

## Using the program as a module

You can also use the package as a python module. For using it as a module, save the program in the current directory or in the computer-discoverable path.
`import EQinfo_extractor as eq_ext`
`eq_ext.eventinfo(1976,3,8,4,39)`

In seismology, we always need to download and check the event information about the events. This python script can download the event catalog from the website to the local computer for the given range of time.

Running this program is very simple:

The user just needs to input the time range for the earthquakes e.g., 2000/05-2009/08.

Requirements: Python 3

## Ray tracing through a 1-D velocity model

Refer to Chapter 4 of Shearer, Introduction to Seismology.

For a ray piercing through Earth, the ray parameter (or horizontal slowness) p is defined by several expressions:

where u = 1/v is the slowness, θ is the ray incidence angle, T is the travel time, X is the horizontal range and utp is the slowness at the ray turning point.

The vertical slowness is defined as:

and integral expressions for the surface-to-surface travel time are:

With these equations we can calculate the travel time (T) and horizontal distance (X) for a given ray with ray parameter p and velocity model v(z).

Apply these to a problem (Exercise 4.8 Shearer):

(COMPUTER) Consider MARMOD, a velocity-versus-depth model, which is typical of much of the oceanic crust (Table 4.1). Linear velocity gradients are assumed to exist at intermediate depths in the model; for example, the P velocity at 3.75 km is 6.9 km/s. Write a computer program to trace rays through this model and produce a P-wave T(X) curve, using 100 values of the ray parameter p equally spaced between 0.1236 and 0.2217 s/km. You will find it helpful to use subroutine LAYERXT (provided in Fortran in Appendix D and in the supplemental web material as a Matlab script), which gives dx and dt as a function of p for layers with linear velocity gradients. Your program will involve an outer loop over ray parameter and an inner loop over depth in the model. For each ray, set x and t to zero and then, starting with the surface layer and proceeding downward, sum the contributions, dx and dt, from LAYERXT for each layer until the ray turns. This will give x and t for the ray from the surface to the turning point. Multiply by two to obtain the total surface-to-surface values of X(p) and T(p). Now produce plots of: (a) T(X) plotted with a reduction velocity of 8 km/s, (b) X(p), and (c) τ(p). On each plot, label the prograde and retrograde branches. Where might one anticipate that the largest amplitudes will occur?

Using the LAYERXT subroutine (Appendix D of Shearer), the FORTRAN and MATLAB can be downloaded from here. In this example, we will use python to calculate and plot the result.

The main program is as follow:

```from ttcal import layerxt
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

p1 = 0.1236
p2 = 0.2217
n = 100

result = [['p','X','T']]
for p in np.linspace(p1, p2, n):
X = 0; T = 0;
for i in range(0, len(v_model)-1):
p = p
h = v_model.h[i+1] - v_model.h[i]
utop = 1/v_model.vp[i]
ubot = 1/v_model.vp[i+1]
dx, dt, irtr = layerxt(p,h,utop,ubot)
X += dx; T += dt
#If the ray reflected at layer the dx, dt calculation stops
if irtr == 2:
break
result.append([p,X,T])

result_df = pd.DataFrame(result, columns=result.pop(0))
result_df = result_df[::-1]

#Multiply by 2 to get total surface-to-surface value ot X(p) and T(p)
result_df['T'] = 2*result_df['T']
result_df['X'] = 2*result_df['X']

plt.figure(1)
result_df['RT'] = result_df['T']-result_df['X']/8
result_df.plot(x='X', y='RT', legend=False)
plt.title('Reduction time')
plt.ylabel('Reduce time T - X/8 (s)')
plt.xlabel('Distance (km)')
plt.savefig('ex4.8.1.png')

plt.figure(2)
result_df.plot(x='p', y='X', legend=False)
plt.title('X(p)')
plt.ylabel('Distance X (km)')
plt.xlabel('Ray parameter (km/s)')
plt.xlim([p1,p2])
plt.savefig('ex4.8.2.png')

plt.figure(3)
result_df['Taup'] = result_df['T'] - result_df['p']*result_df['X']
result_df.plot(x='p', y='Taup', legend=False)
plt.title('Tau(p)')
plt.ylabel('Tau(p) (s)')
plt.xlabel('Ray parameter (km/s)')
plt.xlim([p1,p2])
plt.savefig('ex4.8.3.png')

plt.figure(4)
v_model.h = v_model.h * (-1)
v_model.plot(x='vp', y='h', legend=False)
plt.ylabel('Depth (km)')
plt.xlabel('Velocity (km/s)')
plt.xlim([4,9])
plt.ylim([-10,0])
plt.savefig('ex4.8.4.png')

#plt.show()```

LAYERXT subroutine written in python:

```from math import sqrt, log

def layerxt(p,h,utop,ubot):
''' LAYERXT calculates dx and dt for ray in layer with linear velocity gradient

Inputs:   p     =  horizontal slowness
h     =  layer thickness
utop  =  slowness at top of layer
ubot  =  slowness at bottom of layer
Returns:  dx    =  range offset
dt    =  travel time
irtr  =  return code
= -1,  zero thickness layer
=  0,  ray turned above layer
=  1,  ray passed through layer
=  2,  ray turned within layer, 1 segment counted in dx,dt'''
# ray turned above layer
if p >= utop:
dx = 0
dt = 0
irtr = 0
return dx,dt,irtr

#Zero thickness layer
elif h == 0:
dx = 0
dt = 0
irtr = -1
return dx,dt,irtr

#Calculate some parameters
u1 = utop
u2 = ubot
v1 = 1.0/u1
v2 = 1.0/u2
b = (v2 - v1)/h
eta1 = sqrt(u1**2 - p**2)

#Constant velocity layer
if b == 0:
dx = h*p/eta1
dt = h*u1**2/eta1
irtr = 1
return dx,dt,irtr

x1 = eta1/(u1*b*p)
tau1=(log((u1+eta1)/p)-eta1/u1)/b

#Ray turned within layer, no contribution to integral from bottom point
if p >= ubot:
dx = x1
dtau = tau1
dt = dtau + p*dx
irtr = 2
return dx,dt,irtr

irtr=1
eta2=sqrt(u2**2-p**2)
x2=eta2/(u2*b*p)
tau2=(log((u2+eta2)/p)-eta2/u2)/b

dx=x1-x2
dtau=tau1-tau2

dt=dtau+p*dx
return dx,dt,irtr

def flat(zsph,vsph):
''' Calculate the flatten transformation from a spherical Earth.
zf, vf = flat(zsph,vsph)'''

a = 6371.0

rsph = a - zsph
zf = -a*log(rsph/float(a))
vf = (a/float(rsph))*vsph
return zf, vf```

Input velocity model:

h,vp,vs,den
0.0,4.5,2.4,2
1.5,6.8,3.75,2.8
6.0,7.0,3.5,2.9
6.5,8.0,4.6,3.1
10.0,8.1,4.7,3.1

The codes can be downloaded here: main program, python LAYERXT subroutine, velocity model.

You need to install numpy, pandas and matplotlib package to run this. Normally it can be installed in the terminal by (presume python has been installed):

`pip install numpy pandas matplotlib`

Nguyen Cong Nghia – IESAS

## Introduction to Python Part II

I. Type of objects in Python:

In Python, every object has its own class – or type of data. The in-depth tutorial can be found on the web, for example, https://jeffknupp.com/blog/2014/06/18/improve-your-python-python-classes-and-object-oriented-programming/. In this tutorial, I will introduce some basic type in Python.

To check type of a variable, data, you can use function type(variable)

+ Numbers: most frequently use is float and int type. The float type uses decimal while the int rounds number. Below is the example of using these type of number. To convert to int and float type, we use int(variable) and float(variable). A number can take numeric operations like +(add), – (subtract), *(multiply), /(divide), % (modulo), ** (exponential)

+ Strings: any character information, that in between ‘ ‘ or ” “.  Strings can be joined together by using + operation.

+ List: Python often uses compound data types, used to group together with other values. The most versatile is the list, which can be written as a list of comma-separated values (items) between square brackets. Lists might contain items of different types: numbers, string or list itself (a list of lists). This type is a basic type to analyze data in Python because it makes you able to access data with ordering.

II. Function and method

In Python,methods are associated with object instances or classes; functions aren’t. When Python dispatches (calls) a method, then it binds the first parameter of that call to the appropriate object reference. A function often has the form of function(argument) while an argument can be any kind of data type (number, string, list). A method must be associated with a type of object and often have the form of object.method(argument) with an object is the suitable type to do a method.

Let’s do some practice and take the list as an example.

The functions used in here are len (return the length of an object – how many objects in a list) and print (display an object on screen). The methods used in here are append(add an object to a list) and reverse (reverse the order of objects in a list).

## Introduction to Python (Part I)

Python is a high-level, general-purpose, interpreted programming language. It can be used to do some different purpose: from managing large data, simple to complex calculations to very specific objectives by using the pre-written packages.

There are 2 most common used versions of Python: Python 2 (current version is 2.7) and Python 3 (current version is 3.5). The 2 versions are different in some syntax, which Python 2 is the stable version while Python 3 is the developing version.

I. Basics syntax

A python script can be invoked by the first line ‘#!/usr/bin/python’ similar to using Perl.

Running a python program: we will start with how to run a python program Helloworld.py as below:

```#!/usr/bin/python

hello = 'Hello World!'

print hello```

To run a program in python we use:

`python Helloworld.py`

Addition flags can be used to calibrate the running process, details can be found here. A useful flag is -i (interactive mode) which allows you to control the variable and running process.

Command locals() is used to give a dictionary of local variables. All of the variables with __ are python pre-defined variables. This mode helps to check the variables during your program.

Using as a calculator: Basic mathematics operation can be used in Python, for example:

For more complex calculations, a package named numpy is used to complete these task. A tutorial in using Python package will be presented later.

II. Using external commands

Python allows running external commands using either module: os and subprocess.

```import os

os.system('ls -l')```

or

```from subprocess import call

call(['ls','-l'])```

Nguyen Cong Nghia – IESAS