FQ_walkthrough.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"border:2px solid gray\"> </hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center>FakeQuakes Tutorial</center>\n",
    "## <center>Generates synthetic rupture models and waveforms </center>\n",
    "### <center>Tara Nye</center>\n",
    "### <center>PhD Candidate | University of Oregon</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"border:2px solid gray\"> </hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "FakeQuakes is a set of semistochastic foward modeling codes from the github repository [MudPy](https://github.com/dmelgarm/MudPy).  FakeQuakes was created by __Dr. Diego Melgar__ (professor at the Unviersity of Oregon), and the components are outlined in the paper [__Kinematic rupture scenarios and synthetic displacement data: An example application to the Cascadia subduction zone__ (Melgar et al., 2016)](http://doi.wiley.com/10.1002/2016JB013314).  \n",
    "\n",
    "This set of codes generates synthetic rupture models and can generate low frequency displacement waveforms and broadband acceleration waveforms.\n",
    "\n",
    "For more information on how to use FakeQuakes, you can refer to this [wiki page](https://github.com/taranye96/tsuquakes/wiki/FakeQuakes-Script) on my github account.  This page walks you through how to set up the parameter file (__.fq.py__) listed in the [tsuquakes](https://github.com/taranye96/tsuquakes) respository and how to run the different steps in FakeQuakes.  If you have any other questions, feel free to contact me or Diego Melgar.\n",
    "\n",
    "Tara Nye: tnye@uoregon.edu\n",
    "<br>Diego Melgar: dmelgarm@uoregon.edu"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"./figures/rupture_models.png\"/>\n",
    "\n",
    "<br>__Figure 1.__ Mean slip model for the 2010 M7.8 Mentawai event on the left (modified from Yue et al., 2014), and two example synthetic rupture models patterned after this event in the center and on the right."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"border:1px solid gray\"> </hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center>Table of Contents</center>\n",
    "## 1. <a href=#installation>FakeQuakes Installation</a>\n",
    "### &emsp;&emsp;1.1 <a href=#python>Setup Python Environment</a>\n",
    "### &emsp;&emsp;1.1 <a href=#mudpy>Install MudPy</a>\n",
    "## 2. <a href=#setup>Project Setup</a>\n",
    "### &emsp;&emsp;2.1 <a href=#directories>Setup Directories</a>\n",
    "### &emsp;&emsp;2.2 <a href=#parameters>Define Parameters</a>\n",
    "## 3. <a href=#generate_ruptures>Generate Rupture Models</a>\n",
    "## 4. <a href=#LF>Generate Low Frequency Displacement Waveforms</a>\n",
    "###  &emsp;&emsp;4.1 <a href=#GFs>Compute Green's Functions</a>\n",
    "### &emsp;&emsp;4.2  <a href=#lf_waveforms>Make Low Frequency Waveforms </a>\n",
    "## 5. <a href=#BB>Generate Broadband Acceleration Waveforms</a>\n",
    "### &emsp;&emsp;5.1 <a href=#HF>Make High Frequency Waveforms</a>\n",
    "### &emsp;&emsp;5.2 <a href=#match_filt>Combine Waveforms</a>\n",
    "## 5. <a href=#notes>Final Notes</a>\n",
    "## 6. <a href=#references>References</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='installation'></a>\n",
    "<hr style=\"border:1px solid gray\"> </hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center>FakeQuakes Installation</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='python'></a>\n",
    "***\n",
    "## <center>Setup Python Environment</center>\n",
    "### MudPy requires the following dependencies:\n",
    "- __ObsPy__ **`$ conda install -c conda-forge obspy`**\n",
    "- __Pyproj__ **`$ conda install -c conda-forge pyproj`**\n",
    "- __UTM__ **`$ conda install -c conda-forge utm`**\n",
    "- __mpi4py__ **`$ conda install -c anaconda mpi4py`**\n",
    "- __GCC and GFortgran__\n",
    "    - If you are using a Mac, you can just install __Xcode Developer Tools__ (Apple app store)\n",
    "    - If you are on a Linux machine:\n",
    "        **`$ sudo apt install build-essential`**\n",
    "        **`$ sudo apt-get install gfortran`**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='mudpy'></a>\n",
    "***\n",
    "## <center>Install MudPy</center>\n",
    "### 1. Clone MudPy repository \n",
    "**`$ git clone https://github.com/dmelgarm/MudPy.git`**\n",
    "### 2. Build the fk Green's function code\n",
    "#### - Inside *MudPy/src/fk/* run:\n",
    "**`$ make clean`**\n",
    "\n",
    "**`$ make all`**\n",
    "### 3. Set Paths (.bash_profile or .bashrc script)\n",
    "### - Add the Mudpy src/fk folder to your PATH variable \n",
    "#### i.e. in my .bash_profile script I have:\n",
    "**`export PATH=/Users/tnye/code/MudPy/src/fk:$PATH`**\n",
    "### - Add the Mudpy src/python folder to your PYTHONPATH \n",
    "#### i.e. in my .bash_profile script I have:\n",
    "**`export PYTHONPATH=/Users/tnye/code/MudPy/src/python:$PYTHONPATH`**\n",
    "### - Define the MUD environment variable\n",
    "#### i.e. in my .bash_profile script I have:\n",
    "**`export MUD=/Users/tnye/code/MudPy`**\n",
    "\n",
    "### 4. Make sure everything worked\n",
    "#### If you type these into the terminal, a help screen should appear:\n",
    "\n",
    "**`$ fk.pl`**\n",
    "\n",
    "**`$ syn`**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='setup'></a>\n",
    "<hr style=\"border:1px solid gray\"> </hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center>Project Setup</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='directories'></a>\n",
    "***\n",
    "## <center>Setup Directories</center>\n",
    "### In the FakeQuakes_tutorial folder, there is a files folder with these types of files:\n",
    "- __.gflist__ (list of the stations and their coordinates)\n",
    "- __.mod__ (velocity profile model for the region where the synthetic scenarios are being generated)\n",
    "- __.fault__ (fault model with subfault dimensions that will be used to generate the rupture scenarios)\n",
    "\n",
    "__Note: These are example files, and you will need to create your own for your specific region of study and dataset of stations.__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### __.gflist__\n",
    "\n",
    "The first three columns are the station name, station longitude, and station latitude. The next five columns are flags for what functions you need depending on the application. If you are just generating low-frequency or high-frequency waveforms, you will only want the \"disp\" column to have 1's, and the other columns will just be 0's. The final column is an optional column for kappa if you happen to have station-specific kappa values. If not, FakeQuakes will just use a default value of 0.04 s, or you can specify a singular value in the parameters discussed later on. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!head ./files/sm.gflist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### __.mod__ \n",
    "(Obtained from Yue et al., 2014)\n",
    "\n",
    "The columns are: Layer thickness(km), Vs(km/s), Vp(km/s), Density(g/cm^3), Qs, and Qp.\n",
    "\n",
    "Note: The last row has to be a thickness of 0 (it's assumed to be infinite below this)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!head ./files/mentawai.mod"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### __.fault__ \n",
    "(Modified from Yue et al., 2014)\n",
    "\n",
    "The columns are: Subfault #, longitude, latitude, depth(km), strike, dip, type (*doesn't do anything*), risetime (*doesn't do anything*), subfault length(m), and subfault width(m).\n",
    "\n",
    "Note: The columns that don't do anything are left over from slip inversion file setup, but need to be included in the structure of the file.\n",
    "\n",
    "Note: If you are using a 2D-fault, these subfaults are assumed to be squares, and if you are using a 3D-fault, these subfaults are assumed to be triangles. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!head ./files/mentawai.fault"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### We need to initialize the project folder and put these files into their respective folders"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Standard Library Imports\n",
    "import os\n",
    "import numpy as np\n",
    "from obspy.core import UTCDateTime\n",
    "from obspy import read\n",
    "import time\n",
    "\n",
    "# MudPy Imports\n",
    "from mudpy import fakequakes,runslip,forward,view"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define variables for project location \n",
    "home = os.getcwd() + '/' # This sets the current working directory as the home path\n",
    "project_name='test' # Name of project folder that will be set up in home directory "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you don't have the dependencies installed, or you want to just use the data I created, set __project_name=demo__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initalize project folders\n",
    "if project_name != 'demo':\n",
    "    fakequakes.init(home,project_name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Copy the files into these folders \n",
    "- __.mod__ should be placed in __*/home/project_name/structure/*__ \n",
    "- __.fault__ should be placed in __*/home/project_name/data/model_info/*__\n",
    "- __.gflist__ file should be placed in __*/home/project_folder/data/station_info/*__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cp ./files/mentawai.mod $home/$project_name/structure\n",
    "!cp ./files/mentawai_coarse.fault $home/$project_name/data/model_info\n",
    "!cp ./files/sm.gflist $home/$project_name/data/station_info"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='parameters'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='parameters'></a>\n",
    "***\n",
    "## <center>Define Parameters</center>\n",
    "### Rupture parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Runtime parameters \n",
    "ncpus=2                                        # how many CPUS you want to use for parallelization (needs ot be at least 2)\n",
    "Nrealizations=2                                # Number of fake ruptures to generate per magnitude bin\n",
    "hot_start=0                                    # If code quits in the middle of running, it will pick back up at this index\n",
    "\n",
    "# File parameters\n",
    "model_name='mentawai.mod'                      # Velocity model file name\n",
    "fault_name='mentawai_coarse.fault'             # Fault model name\n",
    "mean_slip_name=None                            # Set to path of .rupt file if patterning synthetic runs after a mean rupture model\n",
    "run_name='mentawai'                            # Base name of each synthetic run (i.e. mentawai.000000, mentawai.000001, etc...)\n",
    "rupture_list='ruptures.list'                   # Name of list of ruptures that are used to generate waveforms.  'ruptures.list' uses the full list of ruptures FakeQuakes creates. If you create file with a sublist of ruptures, use that file name.\n",
    "distances_name='mentawai'                      # Name of matrix with estimated distances between subfaults i and j for every subfault pair\n",
    "load_distances=0                               # This should be zero the first time you run FakeQuakes with your fault model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Source parameters\n",
    "UTM_zone='47M'                                 # UTM_zone for rupture region \n",
    "time_epi=UTCDateTime('2010-10-25T14:42:12Z')   # Origin time of event (can set to any time, as long as it's not in the future)\n",
    "target_Mw=np.array([7.8])                      # Desired magnitude(s), can either be one value or an array\n",
    "hypocenter=None                                # Coordinates of subfault closest to desired hypocenter, or set to None for random\n",
    "force_hypocenter=False                         # Set to True if hypocenter specified\n",
    "rake=90                                        # Average rake for subfaults\n",
    "scaling_law='T'                                # Type of rupture: T for thrust, S for strike-slip, N for normal\n",
    "force_magnitude=True                           # Set to True if you want the rupture magnitude to equal the exact target magnitude\n",
    "force_area=True                                # Set to True if you want the ruptures to fill the whole fault model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Correlation function parameters\n",
    "# Don't change this\n",
    "hurst=0.4                                      # Hurst exponent form Melgar and Hayes 2019\n",
    "Ldip='auto'                                    # Correlation length scaling: 'auto' uses Melgar and Hayes 2019, 'MB2002' uses Mai and Beroza 2002\n",
    "Lstrike='auto'                                 # Same as above\n",
    "slip_standard_deviation=0.9                    # Standard deviation for slip statistics: Keep this at 0.9\n",
    "lognormal=True                                 # Keep this as True to solve the problem of some negative slip subfaults that are produced"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rupture propagation parameters\n",
    "rise_time = 'MH2017'                           # Rise time scaling to use. 'GP2010' uses Graves and Pitarka (2010), 'GP2015' uses Graves and Pitarka (2015), 'S1999' uses Sommerville (1999), and 'MH2017' uses Melgar and Hayes (2017).  \n",
    "rise_time_depths=[10,15]                       # Transition depths for rise time scaling (if slip shallower than first index, rise times are twice as long as calculated)\n",
    "max_slip=40                                    # Maximum sip (m) allowed in the model\n",
    "max_slip_rule=False                            # If true, uses a magntidude-depence for max slip\n",
    "shear_wave_fraction_shallow=0.49               # Shear wave fraction for depths shallower than rise_time_depths[0]\n",
    "shear_wave_fraction_deep=0.8                   # Shear wave fraction for depths depper than rise_time_depths [1] (0.8 is a standard value (Mai and Beroza 2002))\n",
    "source_time_function='dreger'                  # options are 'triangle' or 'cosine' or 'dreger'\n",
    "stf_falloff_rate=4                             # Only affects Dreger STF, 4-8 are reasonable values\n",
    "num_modes=72                                   # Number of modes in K-L expansion\n",
    "slab_name=None                                 # Slab 2.0 Ascii file for 3D geometry, set to None for simple 2D geometry\n",
    "mesh_name=None                                 # GMSH output file for 3D geometry, set to None for simple 2D geometry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Waveform Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Green's Functions parameters\n",
    "GF_list='sm.gflist'                            # Stations file name\n",
    "G_name='sm'                                    # Basename you want for the Green's functions matrices\n",
    "make_GFs=1                                     # This should be 1 to run Green's functions\n",
    "make_synthetics=1                              # This should be 1 to make the synthetics\n",
    "G_from_file=0                                  # This should be zero the first time you run FakeQuakes with your fault model and stations.\n",
    "\n",
    "# fk parameters\n",
    "# used to solve wave equation in frequency domain \n",
    "dk=0.1 ; pmin=0 ; pmax=1 ; kmax=20             # Should be set to 0.1, 0, 1, 20\n",
    "custom_stf=None                                # Assumes specified source time function above if set to None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Low frequency waveform parameters\n",
    "dt=0.5                                         # Sampling interval of LF data \n",
    "NFFT=512                                       # Number of samples in LF waveforms (should be in powers of 2)\n",
    "# dt*NFFT = length of low-frequency dispalcement record\n",
    "# want this value to be close to duration (length of high-frequency record)\n",
    "\n",
    "# High frequency waveform parameters\n",
    "stress_parameter=50                            # Stress drop measured in bars (standard value is 50)\n",
    "moho_depth_in_km=30.0                          # Average depth to Moho in this region \n",
    "Pwave=True                                     # Calculates P-waves as well as S-waves if set to True, else just S-Waves\n",
    "kappa=None                                     # Station kappa values: Options are GF_list for station-specific kappa, a singular value for all stations, or the default 0.04s for every station if set to None\n",
    "hf_dt=0.01                                     # Sampling interval of HF data\n",
    "duration=250                                   # Duration (in seconds) of HF record\n",
    "\n",
    "# Match filter parameters\n",
    "zero_phase=True                                # If True, filters waveforms twice to remove phase, else filters once\n",
    "order=4                                        # Number of poles for filters\n",
    "fcorner_low=0.998                              # Corner frequency at which to filter waveforms (needs to be between 0 and the Nyquist frequency)\n",
    "fcorner_high=0.01                              # Corner frequency at which to filter waveforms (needs to be between 0 and the Nyquist frequency)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Whew...we're finally done setting up the parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='generate_ruptures'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='generate_ruptures'></a>\n",
    "<hr style=\"border:1px solid gray\"> </hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center>Generate Rupture Models</center>\n",
    "This step uses the __fakequakes.generate_ruptures__ function to generate the different rupture models.  The outputs are .rupt files which are the rupture models with information on each subfault, and .log files which give details on the rupture, such as magntidue, hypocenter, origin time, etc...  Most of this is taken from Mai and Beroza (2002).  \n",
    "<br>There are two main types of rupture models you can make:\n",
    "- a rupture model patterned after a real or average rupture model.\n",
    "- a rupture model with completely random slip patterns\n",
    "\n",
    "<br> For this tutorial, we will be making random rupture models for the Mentawai region. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if project_name != 'demo':\n",
    "    #Generate rupture models\n",
    "    fakequakes.generate_ruptures(home,project_name,run_name,fault_name,slab_name,mesh_name,load_distances,\n",
    "        distances_name,UTM_zone,target_Mw,model_name,hurst,Ldip,Lstrike,num_modes,Nrealizations,rake,rise_time,\n",
    "        rise_time_depths,time_epi,max_slip,source_time_function,lognormal,slip_standard_deviation,scaling_law,\n",
    "        ncpus,mean_slip_name=mean_slip_name,force_magnitude=force_magnitude,force_area=force_area,\n",
    "        hypocenter=hypocenter,force_hypocenter=force_hypocenter,shear_wave_fraction_shallow=shear_wave_fraction_shallow,\n",
    "        shear_wave_fraction_deep=shear_wave_fraction_deep,max_slip_rule=max_slip_rule)   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rupture Dimensions\n",
    "<img src=\"./figures/fault.png\" style=\"width:400px;height:182px\"/>\n",
    "The dimensions of the fault rupture are determined by the target magnitude and rupture dimension scaling laws $\\eqref{eq:scaling law length}$$\\eqref{eq:scaling law width}$ (Blaser et al., 2010).\n",
    "\n",
    "\\begin{equation*}\n",
    "log_{10}L = -2.37 + 0.57M_{w}\n",
    "\\label{eq:scaling law length} \\tag{1}\n",
    "\\end{equation*}\n",
    "\n",
    "\\begin{equation*}\n",
    "log_{10}W = -1.86 + 0.46M_{w}\n",
    "\\label{eq:scaling law width} \\tag{2}\n",
    "\\end{equation*}\n",
    "\n",
    "To add some variability in the rupture dimensions of events of the same magnitude, FakeQuakes uses a stochastic approach to the scaling laws $\\eqref{eq:stochastic scaling law}$.\n",
    "\n",
    "\\begin{equation*}\n",
    "log_{10}L \\sim N(a_{1} + b_{1}M_{w}, s^{2}_{xy})\n",
    "\\label{eq:stochastic scaling law} \\tag{3}\n",
    "\\end{equation*}\n",
    "\n",
    "### Slip Statistics\n",
    "<img src=\"./figures/slip.png\" style=\"width:400px;height:182px\"/>\n",
    "\n",
    "Unless a hypocenter is specified, a random subfault is chosen at the hypocenter.  \n",
    "<br>The slip on each subfault is normally distributed, and the vector __s__ contains all of the subfault slips $\\eqref{eq:slip vector}$.\n",
    "\n",
    "\\begin{equation*}\n",
    "s \\sim N(\\mu, \\hat{C})\n",
    "\\label{eq:slip vector} \\tag{4}\n",
    "\\end{equation*}\n",
    "\n",
    "<img src=\"./figures/slip_pdf.png\" style=\"width:440px;height:350px\"/>\n",
    "\n",
    "To produce a synthetic event with the desired magnitude, the mean vector, $\\mu$, has uniform slip and enough scalar moment to match the target magnitude we specified $\\eqref{eq:moment}$.\n",
    "\n",
    "\\begin{equation*}\n",
    "M_{0} = \\rho LW\\mu\n",
    "\\label{eq:moment} \\tag{5}\n",
    "\\end{equation*}\n",
    "\n",
    "$\\hat{C}$ is the covariance matrix of the slip distribution and is a function of the standard deviation of slip at each subfault and the correlation between each ith and jth subfault pair.  FakeQuakes uses the __Von Karman correlation function__ $\\eqref{eq:correlation function}$ because it represents fault slip patterns well. \n",
    "\n",
    "\\begin{equation*}\n",
    "C_{ij}(r_{ij}) = \\frac{G_H(r_{ij})}{G_0(r_{ij})}\n",
    "\\label{eq:correlation function} \\tag{6}\n",
    "\\end{equation*}\n",
    "\n",
    "\\begin{equation*}\n",
    "G_{H} = r^{H}_{ij}K_{H}(r_{ij})\n",
    "\\label{eq:GH} \\tag{7}\n",
    "\\end{equation*}\n",
    "\n",
    "- r$_{ij}$ is a length measure between subfaults i and j<br>(__Note__: This is not the actual distance.  It is a value dependent on the along strike and along dip distances between subfaults)\n",
    "- __H__ is the hurst exponent that we defined in the parameters\n",
    "- K$_{H}$ is the modified Bessel function of the second kind\n",
    "    - solution to Bessel differential equation which is used when solving the 3D wave equation\n",
    "\n",
    "If $\\lambda$$_{k}$ and v$_{k}$ are the eigenvalues and eigenvectors of the correlation function C$_{ij}$, then the slip vector, s, can be expressed by the __Karhunen Loeve (K-L) expansion__ $\\eqref{eq:K-L expansion}$.  \n",
    "\n",
    "\\begin{equation*}\n",
    "s = \\mu + \\sum\\limits_{k=1}^N z_{k}\\sqrt{\\lambda_{k}v_{k}}\n",
    "\\label{eq:K-L expansion} \\tag{8}\n",
    "\\end{equation*}\n",
    "\n",
    "- N is the number of modes specified, and it can be any number up to the number of subaults in your mode.  The more modes you have, the better you can model the high frequency stuff, but it also takes longer to compute. \n",
    "- z$_{k}$ are normally distrubuted random numbers between 0-1 that are used as weights for each eigenmode. This is what gives some randomness to the slip patterns.\n",
    "\n",
    "\n",
    "### Kinematics \n",
    "<img src=\"./figures/rupture.gif\" style=\"width:400px;height:182px\"/>\n",
    "Now that we have the static slip distribution, we need to determine the kinematics of the rupture so that we can generate the waveforms.  \n",
    "\n",
    "The rupture onset time of each subfault is determined __(2-stage)__\n",
    "1. Background rupture distribution $\\eqref{eq:Vrupt}$ is determined using the shear-wave speed from the velocity model (V$_{s}$) and the shear-wave fraction parameter we set (SF)\n",
    "\\begin{equation*}\n",
    "V_{r} = (SF)(V_{s})\n",
    "\\label{eq:Vrupt} \\tag{9}\n",
    "\\end{equation*}\n",
    "    - There is a reduction in rupture speed at shallow depths, as represented by the small shear_wave_fraction_shallow compared to shear_wave_fraction_deep.  This is to account for the low shear-wave velocity at shallow depths in the accretionary wedge.\n",
    "    - Subfault onset times from this background distribution are perturbed so that there's faster propagation where slip is large (Graves and Pitarka, 2010).\n",
    "    - Risetime of each subfault (duration of slip) scaled by __√__(total slip on subfault) so that larger slips have longer risetimes\n",
    "2. A slip rate function is applied to determine how the slip evolved over time. \n",
    "\n",
    "\n",
    "### Output \n",
    "This step results in two output file types: .log and .rupt files. \n",
    "- .log file: A file for one run containing the overall parameters of the rupture (e.g. Project name, earthquake magnitude and coordinates, etc..)\n",
    "- .rupt file: A file for one run containing the kinematics of the rupture (e.g. slip on each subfault, rupture onset time, slip duration, etc...)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### View .log file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logpath=os.getcwd() +'/' + project_name + '/output/ruptures/mentawai.000000.log'\n",
    "var=!cat $logpath\n",
    "var"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### View .rupt file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ruptpath=f'{home}/{project_name}/output/ruptures/mentawai.000000.rupt'\n",
    "!head $ruptpath"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quick view of the rupture models: Coarse Fault Model\n",
    "\n",
    "We are using a coarse (large subfaults) fault model for this tutorial to save computation time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "#%matplotlib notebook\n",
    "view.quick_model(f'{home}/{project_name}/output/ruptures/mentawai.000000.rupt')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quick view of example rupture models: Fine Fault Model\n",
    "\n",
    "This is an example of what are more fine-scale resolution (smaller subfaults) fault model would look like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "#%matplotlib notebook\n",
    "view.quick_model(f'{os.getcwd()}/figures/mentawai_fine.000000.rupt')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='LF'></a>\n",
    "<hr style=\"border:1px solid gray\"> </hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center>Generate Low Frequency Displacement Waveforms</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='GFs'></a>\n",
    "***\n",
    "## <center>Compute Green's Functions</center>\n",
    "This step uses the __runslip.inversion__ function to generate the Green's Functions, which are seismograms for each subfault-station pair...one step closer to getting the final waveforms!  This uses the frequency-waveform (fk) integration code from Zhu and Rivera (2002)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if project_name != 'demo':\n",
    "    # I'm increasing the ncpus here to speed up computation time. If I had set ncpus to 4 initially, it would have tried to make 4 rupture models.\n",
    "    ncpus=4  \n",
    "    runslip.inversionGFs(home,project_name,GF_list,None,fault_name,model_name,\n",
    "        dt,None,NFFT,None,make_GFs,make_synthetics,dk,pmin,\n",
    "        pmax,kmax,0,time_epi,hot_start,ncpus,custom_stf,impulse=True) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Green's Functions (GFs)\n",
    "Green's functions are what the seismogram would look like for 1 meter of displacement, assuming a delta source.   \n",
    "<br> To calculate the GFs, we assume a layered half-space model and use a propagator matrix, which is what we use to go from stress changes along the fault to displacement at the surface.\n",
    "<img src=\"./figures/half_space.png\"/>\n",
    "<br>__Figure 2.__ Half-space model from Zhu and Rivera, 2002.\n",
    "\n",
    "### Synthetics\n",
    "The GFs provide the static displacement, but we need to know how displcement evolved over time.  To obtain this, we convolve the GFs (*G*(*t*)) with a STF (*X*(*t*)) $\\eqref{eq:synthetic}$.  FakeQuakes refers to these products as synthetics (*S*(*t*)).\n",
    " \n",
    "\\begin{equation*}\n",
    "S(t) = X(t)\\circledast G(t)\n",
    "\\label{eq:synthetic} \\tag{9}\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"./figures/subfault_seismograms.png\" style=\"width:400px;height:250px\"/>\n",
    "\n",
    "<br>__Figure 3.__ Example subfault seismograms for one station (blue triangle)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='lf_waveforms'></a>\n",
    "***\n",
    "## <center>Make Low Frequency Waveforms</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the __waveforms_fakequakes__ function to get the displacement waveforms at the different stations in our stations file using the synthetics we obtained in the last step.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if project_name != 'demo':\n",
    "    #Make low frequency waveforms\n",
    "    forward.waveforms_fakequakes(home,project_name,fault_name,rupture_list,GF_list,\n",
    "                model_name,run_name,dt,NFFT,G_from_file,G_name,source_time_function,\n",
    "                stf_falloff_rate)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The generation of LF waveforms is deterministic, which means it is based off of our understanding of how seismic waves propagate through a media.  This is done using a 1D velocity model and GFs.\n",
    "\n",
    "The GFs assume slip of 1 m, so we have to scale them to match the actual slip on the fault by multiplying each synthetic by the slip on that subfault.  The final waveforms for each station are the sum of the subfault synthetics associated with the station."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## Results...Finally!!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read in two of the displacement waveforms\n",
    "\n",
    "The 'LY' channel code refers to low-frequency displacement waveforms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "st1 = read(f'{home}{project_name}/output/waveforms/mentawai.000000/MNSI.LYE.sac')\n",
    "tr1 = st1[0]\n",
    "st2 = read(f'{home}{project_name}/output/waveforms/mentawai.000000/PPSI.LYE.sac')\n",
    "tr2 = st2[0]\n",
    "\n",
    "tr1.plot()\n",
    "tr2.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PPSI has a static offset (final displacement does not equal initial) because it is closer to the source."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='BB'></a>\n",
    "<hr style=\"border:1px solid gray\"> </hr>\n",
    "\n",
    "# <center>Generate Broadband Acceleration Waveforms</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='HF'></a>\n",
    "***\n",
    "## <center>Make High Frequency Waveforms</center>\n",
    "We can use the __hf_waveforms__ function to get the high frequency acceleration waveforms at the different stations in our stations file.  Most of this is taken from Graves and Pitarka (2010, 2015)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if project_name != 'demo':    \n",
    "    # Make high-frequency waveforms\n",
    "    forward.hf_waveforms(home,project_name,fault_name,rupture_list,GF_list,\n",
    "                    model_name,run_name,dt,NFFT,G_from_file,G_name,rise_time_depths,\n",
    "                    moho_depth_in_km,ncpus,source_time_function=source_time_function,\n",
    "                    duration=duration,stf_falloff_rate=stf_falloff_rate,hf_dt=hf_dt,\n",
    "                    Pwave=Pwave,hot_start=hot_start,stress_parameter=stress_parameter)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The generation of HF waveforms is done stochastically.  This means that there is some randomness to the process, but it is based off of probability distributions.  It would be computationally-expensive to run the HF part determinisitically and would require a much more detailed model of the subsurface.  \n",
    "<br>This step works in the frequency domain to generate the waveforms.  Each subfault contributes an acceleration amplitude spectrum for each station $\\eqref{eq:hf_spectra}$.\n",
    "\n",
    "\\begin{equation*}\n",
    "A_{i}(f) = \\sum\\limits_{j=1,M} C_{ij}S_{i}(f)G_{ij}(f)P(f)\n",
    "\\label{eq:hf_spectra} \\tag{10}\n",
    "\\end{equation*}\n",
    "\n",
    "- C$_{ij}$ is the radiation scale factor\n",
    "- S$_{i}$(f) is the source radiation spectrum for subfault *i*\n",
    "- G$_{ij}$(f) is the path term, which accounts for intrinsic attenuation along the path\n",
    "- P(f) is the high-frequency spectral decay from site effects\n",
    "\n",
    "To add some \"randomness\", the spectra are combined with the random phase of a Gaussian white noise time series.  <br>Then the spectra are Fourier-transformed back into the time domain, giving an acceleration time series for each subfault-station pair. \n",
    "<br>As with the low-frequencies, all of the subfault seismograms for each station are then summed up to produce the final waveforms for each station. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read in two of the acceleration waveforms\n",
    "\n",
    "The 'HN' channel code refers to high-frequency acceleration waveforms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "st1 = read(f'{home}/{project_name}/output/waveforms/mentawai.000000/PPSI.HNE.mpi.sac')\n",
    "tr1 = st1[0]\n",
    "st2 = read(f'{home}/{project_name}/output/waveforms/mentawai.000000/MNSI.HNE.mpi.sac')\n",
    "tr2 = st2[0]\n",
    "\n",
    "tr1.plot()\n",
    "tr2.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='match_filt'></a>\n",
    "***\n",
    "## <center>Combine Waveforms</center>\n",
    "We can use the __match_filter__ function to combine the low frequency and high frequency waveforms to get final broadband waveforms for the stations in our station file. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if project_name != 'demo':\n",
    "    # Combine LF and HF waveforms with match filter                              \n",
    "    forward.match_filter(home,project_name,fault_name,rupture_list,GF_list,zero_phase,order,\n",
    "                         fcorner_high=fcorner_high,fcorner_low=fcorner_low)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* The LF displacement waveforms are double-differentiated into acceleration waveforms.\n",
    "* The LF waveforms are resampled to have same sampling rate as the HF waveforms.\n",
    "* The LF and HF waveforms are trimmed to have the same length.\n",
    "* The LF waveforms are lowpass filtered and the HF waveforms are highpass filtered.\n",
    "* The LF and HF waveforms are added together to get final broadband waveforms. \n",
    "\n",
    "<img src=\"./figures/match_filter_time.png\"/>\n",
    "\n",
    "<br>__Figure 3.__ Example low frequency and high frequency acceleration waveforms before filtering (left) and after filtering (right).  The final broadband waveforms is a sum of the two filtered waveforms.\n",
    "<img src=\"./figures/match_filter.png\"/>\n",
    "\n",
    "<br>__Figure 4.__ Example low frequency and high frequency acceleration Fourier amplitude spectra before filtering (left) and after filtering (right).  The final broadband waveforms is a sum of the two filtered waveforms."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read in two of the broadband acceleration waveforms\n",
    "\n",
    "The 'bb.HN' channel code refers to broadband acceleration waveforms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "st1 = read(f'{home}/{project_name}/output/waveforms/mentawai.000000/PPSI.bb.HNE.sac')\n",
    "tr1 = st1[0]\n",
    "st2 = read(f'{home}/{project_name}/output/waveforms/mentawai.000000/MNSI.bb.HNE.sac')\n",
    "tr2 = st2[0]\n",
    "\n",
    "tr1.plot()\n",
    "tr2.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='notes'></a>\n",
    "<hr style=\"border:1px solid gray\"> </hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center> Final Notes </center>\n",
    "* All rupture files are stored in __/output/ruptures/__\n",
    "* All waveforms are stored in __/output/waveforms/__\n",
    "    * __.LY*.sac__ files are the low-frequency displacement files\n",
    "    * __.HN*mpi.sac__ files are the high-frequency acceleration files\n",
    "    * __.bb.HN*.sac__ files are the final broadband acceleration waveforms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='references'></a>\n",
    "<hr style=\"border:1px solid gray\"> </hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References\n",
    "Blaser, L., F. Krüger, M. Ohrnberger, and F. Scherbaum (2010), Scaling relations of earthquake source parameter estimates with special focus on subduction environment, Bull. Seism. Soc. Am., 100, 2914–2926.\n",
    "\n",
    "Graves, R., and A. Pitarka (2015), Refinements to the Graves and Pitarka (2010) broadband ground-motion simulation method, Seismol. Res. Lett., 86(1), 75–80.\n",
    "\n",
    "Graves, R. W., & Pitarka, A. (2010), Broadband Ground-Motion Simulation Using a Hybrid Approach. Bulletin of the Seismological Society of America, 100(5A), 2095-2123. doi:10.1785/0120100057\n",
    "\n",
    "Mai, P. M., and Beroza, G. C. (2002), A spatial random field model to characterize complexity in earthquake slip, J. Geophys. Res., 107(B11), 2308, doi:10.1029/2001JB000588\n",
    "\n",
    "Melgar, D., Leveque, R. J., Dreger, D. S., & Allen, R. M. (2016), Kinematic rupture scenarios and synthetic displacement data: An example application to the Cascadia subduction zone. Journal of Geophysical Research: Solid Earth, 121(9), 6658-6674, doi:10.1002/2016jb013314\n",
    "\n",
    "Yue, H., Lay, T., Rivera, L., Bai, Y., Yamazaki, Y., Cheung, K. F., … Muhari, A. (2014), Rupture process of the 2010 Mw 7.8 Mentawai tsunami earthquake from joint inversion of near-field hr-GPS and teleseismic body wave recordings constrained by tsunami observations. AGU: Journal of Geophysical Research, Solid Earth, 119, 5574–5593, doi:10.1002/2014JB011082\n",
    "\n",
    "Zhu, L., & Rivera, L. A. (2002), A note on the dynamic and static displacements from a point source in multilayered media. Geophysical Journal International, 148(3), 619-627, doi:10.1046/j.1365-246x.2002.01610.x"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}