:===============================================================================:
: FD2D = Finite Difference (elastic) seismic wave propagation in 2D. both for SH
: waves (computationally equivalent to acoustic wave propagation) and P-SV waves
: adjoint = it is also possible to run tomographic inversions using this code.
:===============================================================================:

copyright:
Nienke Blom
Christian Boehm
Andreas Fichtner

You're free to use this code, but it would be really kind to acknowledge our work.
The publication relating to this code is Blom, Boehm, Fichtner (GJI 2017) - 
"Synthetic inversions for density using seismic and gravity data" -- doi link:
https://doi.org/10.1093/gji/ggx076

The following is what you can do.
- run wave propagation in 2D for a single or multiple simultaneous sources.
- make a movie of said wave propagation
- calculate gravity data for your model.
- run a tomographic inversion of 2D P-SV and/or SH wave propagation using 
  waveform inversion (adjoint method). This can be done with multiple events.
- include gravity measurements into the inversion.
- look at the results of various steps
- tune the inversion to your specific needs using input files.

INVERSION SCHEME
The inversion scheme implemented is multi-scale (Bunks 1995) and makes use of a 
L-BFGS update algorithm (see e.g. Nocedal & Wright 1999). The multi-scale aspect 
of the observed data is implemented in a rather simplistic way: instead of 
generating high-frequency data once and then filtering them as needed, the wave
propagation is run separately with differently filtered source time functions. This
is computationally not extremely efficient, but otherwise fine.

The most important files for running an inversion will be 
./input/input_parameters.m                        --> change all the input here.
./InvTlbx_callback_fn/inversion_routine_InvTlbx.m --> the wrapper script that does 
                                                      it all.

OUTPUT DIRECTORY:
The inversion will generate output (or load it, sometimes) which will be put in a
directory ./output/[project-name]. This project name is defined in the input params.

RUNTHROUGH:
Inside the inversion_routine file:

Some preparatory steps:
- either load or generate "observed" data for the target model. This target model
  is defined in input_parameters. Obs data may also be present in a file (I ran
  many inversions with the same obs data, no use in running the wave propagation 
  every time again) which is expected to be in [projectfolder]/obs.all-vars.mat
- preparing starting model (may be loaded from ./models or generated) and 
  calculating the initial misfit. In a multi-scale approach, misfit is calculated
  for each frequency band separately.

The inversion consists of the following steps, which are executed until the max
number of steps is reached.
- running a forward propagation for SH and/or P-SV wave propagation
- (calculating gravity data, if required)
- generating adjoint sources using a misfit functional based on the forward 
  simulation results and the observed data.
- running an adjoint simulation for SH and/or P-SV wave propagation in order to
  obtain gradients for the seismic data
- (calculating the gravity gradient, if required, and combining the gradients,
  if required)
- calculate a model update using the L-BFGS algorithm. For this, it may be 
  necessary to run calculate the misfit for additional models, or calculate their
  gradient. This means that additional forward and adjoint wave propagations may
  be executed.

At each step, the process can be visualised:
- the model (various parametrisations)
- the wave propagation (real-time)
- the seismograms (and difference with observed seismograms, if applicable)
- the adjoint kernels (in various parametrisations)
- a summary figure describing the misfit development, differences with the target
  model and information on the kernels.


---------------------------------------------------------------------------------

How to do those things?

Tuning parameters (in ./input/input_parameters)
This script determines all the properties of both domain and inversion. Walk through it carefully to see what all the parameters do.

Run a forward simulation.
(in ./code/run_forward)
In order to do that, you edit the file input/input_parameters.m to your liking. You can choose any of a number models, either or both of wave propagation SH / PSV, source and receiver configurations, and a lot of other stuff too. Just take a look at it. Beware with grid size & stability issues though :) Just take a look at the header of the run_forward file to see how to run it.
Oh you can also make a movie! 


Run an adjoint simulation.
(in ./code/run_adjoint)
Backpropagation of receiver residuals and interaction with the forward field! Once you've run the forward simulation, you'll see that in the output there are two gigantic variable called v_forward and u_forward. (the name is a bit tricky since the forward fields are stored _backwards_ in time). You also have those sources x y and z from the make adjoint source thingy, and now you can backpropagate it all, and compute kernels (with the help of that v_forward and u_forward) on the fly. Cool.
A movie can be made of this process, too, and the resulting kernels can be plotted in all sorts of parametrisations. (see Andreas' book for more detail on parametrisations)

-------------------------------------------------------------------------------------

Pour le reste:
I would like to stress that it's ugly and barely functional but it's kind of cool to be able to play around with wave propagation AND SEE IT. Note the absorbing boundaries which are a Gaussian taper: just multiply the wavefields by smaller and smaller number as you approach the boundaries of the field. Also note that you can do 2nd order or 4th order accuracy of the wave eq. I do everything with 4th order - Andreas says 2nd is not really usable in any way.

note to self:
Like Gerhard explains in his course, you COULD use a combination of normal and 45 degrees rotated grid to do the forward calculations, which'll hugely reduce the nr of grid points needed. Might be worth looking into, but it'll take quite a bit of coding effort to implement so probably I'll leave it at this. That means that there's still considerable numerical dispersion -- as you'll see when you look at the seismograms.