:===============================================================================:
: FD2D = Finite Difference (elastic) seismic wave propagation in 2D. both for SH
: waves (computationally equivalent to acoustic wave propagation) and P-SV waves
:===============================================================================:

The following is what you can do.
- run a single wave propagation, for example to generate 'observed data'
- make a movie of said wave propagation
- run a tomographic inversion of 2D P-SV and/or SH wave propagation, either with
  --> cross-correlation travel time residuals, or
  --> waveform inversion (adjoint method)
- look at the results of various steps
- tune the inversion to your specific needs using input files.

The inversion consists of the following steps, which are executed until the max
number of steps is reached (future work: or until some misfit criterion is met)
- running a forward propagation for SH and/or P-SV wave propagation
- 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
- calculating the descent step length with a line search, based on the previous 
  model and adjoint kernels
- applying hard constraints of total mass and moment of inertia, if requested

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)
- the density update as a result of the application of hard constraints

Future work:
- adding gravity to the inversion. Currently, some gravity information is used
  through the mass and moment of inertia hard constraints, but I want to add
  gravity receivers and add the difference between observed & model gravity as
  separate data to the misfit functional.

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

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.

Running an inversion (in ./code/inversion_routine)
This is the most important thing, for my aims at least. The inversion can be adjusted using the script inversion_routine itself, but most importantly in the file input_parameters.
The first step here is to generate observables with the function prepare_obs(modelnr). The input modelnr determines which model (found in define_material_parameters) is the 'real' model, and the rest of the model properties are taken from input_parameters.
Next, one can simply execute the inversion script. Depending on how much you visualise, how big your domain is, how many timesteps you have and how many iterations you do, this can take a while. An inversion with no wave propagation visualisation (but plotting & saving all other output), a model size (450x390), 3500 time steps and 7 iterations takes about 40 minutes on my machine (32GB RAM)


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! That is basically what I'm so disproportionately excited about. If you just fiddle around with the stf you suddenly see what happens to the X and Z components when it points in a different direction. SO COOL.


Look at the seismograms.
(in ./tools/plot_recordings)
What it says. You can look at the seismograms. Check the code to see what input is required.


Make adjoint sources.
(in ./tools/make_adjoint_sources)
This function is the core of the inversion, as it determines the type of misfit and therefore the type of adjoint source. You can choose the misfit and what data to compare it to. From the forward field you got all the receiver information in X Y and Z directions, observables (if present) are compared to that and now we make sources from that by picking the wavefields.
Currently, the following misfit types can be used:
- a cross-correlation travel time difference
- a waveform misfit using an L2 norm of the waveform difference
Future work: adding other misfits that are more targeted at density, and adding gravity to at least the waveform misfit.



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 just multiply the wavefields by smaller and smaller number as you approach the boundaries of the field. Apparently only works in FD. 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.
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.