This project is about doing tomography with Surface Waves, using a Finite Frequency forward theory to linearize the tomographic problem properly and implement the Subtractive Optimally Localized Averages (SOLA) inversion (Backus and Gilbery 1967, 1968, 1970; Pijpers and Thompson 1992, 1994; Zaroli 2016).

The treeview of this project is as follows :

                ---/observed_seismograms/
                |
                ---/stations/
	        |
   ---/inputs/-----/sources/
   | 
/-----/codes/------/MINEOS/
   |
   ---/outputs/------/measurement/
   |
   ---/figures/------/crust_measurements/

The directory ./MINEOS/ contains fortran routines to compute synthetic seismograms based on normal modes summation (developped by F. Gilbert, John Woodhouse and Guy Masters). The program consists of four successive routines. minos_bran.f computes normal modes given a radial model and eigen.f applies some post-processing to the results. The rountine green.f computes the green functions for a station given the source location and depth using the normal mode solution obtained from the two first routines. The program syndat computes the seismograms given the source solution using the green functions computed at the previous step. 

The directory ./FFSWK/ contains fortran routines to compute the Finite Frequency sensitivity Kernels. It uses some fortran routines (minos_bran.f and read_eigf.f) from Mineos to solve the eigenvalue-eigenfunction problem given a reference model of the Earth (see Mineos user manual Masters et al. 2014). Fortran routines (ffswk.f90 together with ffswk-sub.f90) have been developped (Zhou et al 2004, 2009) to compute the kernels given the eigenproblem solutions, a reference model, source solution and observer location. The kernels are for phase delay data based on parameters that are horizontally and vertically P and S wave velocities, intermediate direction wave speed and density.

normal_modes.py is a python script with functions that act as a wrapper around the C and fortran routines of MINEOS and FFSWK.

user_defined_io.py is python script that contains read/write functions to manage internal or external I/O with various formats.

utility_functions.py is a python script that contains functions designed to execute very specific tasks such as coordinate conversions, normalizations for example.

parameterizations.py is a python script that contains classes within various branches related to various kind of model parameterizations.

graphics.py is a python script written as object_oriented programing that contains simple functions to make various plots.

make_figures.py is a python script to build more complex figures using the functions in the graphics.py script.

colortomo.py is a python script that contains the definition of specific colorscales.

README is this ASCII file that describes the project.


------------------- Measurements -------------------------

Three python scripts in sequence measure phase delays.

get_data.py download station and observed seismogram data from IRIS.
	Arguments :
		The basename that link the data for the all process (ex. cpacific)
		The source name (ex. C201303011320A)
		The station name (ex. IU_KIP)

crust.py computes the path averaged radial model that contains the crust (from CRUST1.0)
         it also computes the synthetic seismograms for the reference mode and the 
	 reference model with the crust
	Arguments :
		The basename
		The source name
		The station name

The reference model to use is stored in /inputs/input_models/ and it is refered in the config file ffswk_sola_params.cfg.
The model with the crust is written in /inputs/input_models/
I/O for synthetic computations are specified in the part about MINEOS routines.
A plot for the reference model and reference with crust model is saved in /figures/crust_measurement/

measurement_MT.py make the multitapering measurement
	Arguments :
		1st possibility synthetic wo crust wrt synthetic w crust
			The basename
			The source name
			The station name
			1
			The crust basename {basename}crust{source}{network}{stationitself}
			Ex. 'cpacific' 'C201303011320A' 'IU_KIP' 1 'cpacificcrustC201303011320AIUKIP'
		2nd possibility synthetic w crust wrt observed
			The crust basename
			The source name
			The station name
			0

Measured data file is saved in /output/measurements/
Some plots (map, seismograms, dispersion curve) are saved in /figures/crust_measurement/


############# Crustal corrections ################

The crustal correction is based on the idea of Woodhouse that the local phase velocity laterally is the phase velocity of a radial Earth with the local properties. Therefore we discretise laterally the Earth and compute the phase velocity from normal mode theory for radial models made of the reference model with the local crustal properties on top. We then accumulate along a ray path the phases for each of these. The phase delay due to the crust is the phase accumulated along the ray from the reference model alone minus the former.

More formally:

r: radius in km
w: frequency in Hz
c: phase velocity in km/s
s: location along some line path
ds: length element along some line path [km]

m^r( r) : reference 1D radial velocity model.
m^rc_j( r) : reference 1D radial velocity model with 1D CRUST1.0 on top evaluated at lateral location j.

c^r( w) : phase velocity model
c^{rc}_j( w) : phase velocity model at lateral location j obtained from normal mode solution of reference model with crust on top.

Phase accumulated d [cycle] in phase velocity model c [km/s] is:
\int_{ s \in path} w/c ds

G_{ij}: length of path i in cell j [km]

d^r_i( w) : accumulated phase along source-receiver path i in phase velocity model obtained from reference radial model velocity model. In cycle.
d^r_i( w) = \sum_j[ G_{ij}] c^r( w)

d^rc_i( w) : accumulated phase along source-receiver path i in phase velocity model obtained from reference radial model with 3D crustal model on top. In cycle.
d^rc_i( w) = \sum_j[ G_{ij} c^rc_j( w)]

d^c_i(w) : phase delay due to crust. In cycle.
d^c_i( w) =  d^rc_i(w) - d^r_i(w).


## Radial reference with crust

crust_get_radial.py

This script takes a 1D radial reference model and CRUST1.0 (specifying crustal layers with varying lateral depths, velocities and density). For each lateral location, it extracts a 1D radial model specific for that location with the reference 1D radial model with the crust properties at this location. Writes the 1D radial models with crust on top for every location in files with MINEOS format. This script is sequential and rather fast (~3 minutes) thanks to a good optimisation in grid management (no use of whereami function).

Notes:
- A layer may appear when adding crust on top of reference. For example because reference decreases to a particularly small value upward but bottom of crust is high.

Reference  Crust       Reference + Crust

| |        |   |       |   |
| |        |   |       |   |
| |      + |   |    =  |   |  /_____ slow velocity layer!
|  \       |           |  \   \
|   \      |           |   \
|    \     |           |    \
|_______   |_______    |_______

- Honoring discontinuities. Mineos can honor discontinuities. Which one to consider as discontinuities in crust model? Solid-Fluid for sure (when an ocean is present).
- Attenuation not in CRUST1.0. How to feed MINEOS with it? Shall not use reference value because may need to be consistent with velocities (then taken from crust).

## Radial reference with crust normal modes

crust_paral.py

Computes the normal mode solution for every radial model (reference with crust on top) corresponding to all lateral locations, that have been obtained with crust_get_radial.py. This script takes a number of CPUs in argument to run in parallel. This is a reather expensive task (~6hrs on 20 CPUs). It runs minos_bran from MINEOS only to obtain frequency dependent phase velocities of modes, used later to compute predicted phase accumulated in each path length along a full path.

crust_compute_fwd.py

For a set of source-receiver locations, computes the lengths in each cell crossed by the ray path. A source-receiver pair is denoted i, a cell is denoted j, the script computes Gij where Gij is the length of path i in cell j. The phase velocity model obtained from reference model with crust on top being mj, the phase accumulated in the model with crust is di=Gijmj.
