{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# BurnMan - Compute seismic velocities - simple case"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Copyright (C) 2012 - 2015 by the BurnMan team, released under the GNU GPL v2 or later.\n",
    "Released under GPL v2 or later.\n",
    "\n",
    "This example script is intended for absolute beginners to BurnMan.\n",
    "We cover importing BurnMan modules, creating a composite material,\n",
    "and calculating its seismic properties at lower mantle pressures and\n",
    "temperatures.  Afterwards, we plot it against a 1D seismic model\n",
    "for visual comparison. Run the seperate cells of the code with shift+enter\n",
    "\n",
    "requires:\n",
    "\n",
    "- geotherms\n",
    "\n",
    "- seismic models\n",
    "\n",
    "- compute seismic velocities\n",
    "\n",
    "- composite mineral helpers\n",
    "\n",
    "teaches:\n",
    "\n",
    "- changing composition\n",
    "- changing geotherm\n",
    "- computing self-consistent depth to compare to seismology"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Import modules"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we import standard python modules that are required for \n",
    "usage of BurnMan.  In particular, numpy is used for handling\n",
    "numerical arrays and mathematical operations on them, and \n",
    "matplotlib is used for generating plots of results of calculations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline \n",
    "import os, sys, numpy as np, matplotlib.pyplot as plt\n",
    "sys.path.insert(1,os.path.abspath('../..'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we import the relevant modules from BurnMan.  The burnman\n",
    "module imports several of the most important functionalities of\n",
    "the library, including the ability to make composites, and compute\n",
    "thermoelastic properties of them.  The minerals module includes\n",
    "the mineral physical parameters for the predefined minerals in\n",
    "BurnMan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import burnman\n",
    "from burnman import minerals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Import seismic model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we create and load the PREM seismic velocity model, which will be\n",
    "used for comparison with the seismic velocities of the \"rock\" composite"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "seismic_model = burnman.seismic.PREM()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create an array of 20 depths at which we want to evaluate PREM, and then\n",
    "query the seismic model for the pressure, density, P wave speed, S wave\n",
    "speed, and bulk sound velocity at those depths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "depths = np.linspace(750e3, 2800e3, 20)\n",
    "pressure, seis_rho, seis_vp, seis_vs, seis_vphi = \\\n",
    "                        seismic_model.evaluate(['pressure','density','v_p','v_s','v_phi'],depths)\n",
    "\n",
    "plt.plot(pressure/1.e9,seis_vs/1.e3,'k',label='Vs')\n",
    "plt.plot(pressure/1.e9,seis_vp/1.e3,'b',label='Vp')\n",
    "plt.plot(pressure/1.e9,seis_rho/1.e3,'r',label='Vphi')\n",
    "plt.plot(pressure/1.e9,seis_vphi/1.e3,'g',label='rho')\n",
    "plt.xlabel('pressure (GPa)')\n",
    "plt.ylabel('velocity (km/s) density (kg/m^3)')\n",
    "plt.xlim(min(pressure)/1.e9,max(pressure)/1.e9)\n",
    "plt.title('PREM')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Input composition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define\n",
    "composite object and name it \"rock\".  A composite is made by\n",
    "giving burnman.composite a list of minerals and their molar fractions.\n",
    "Here \"rock\" has two constituent minerals: it is 80% Mg perovskite\n",
    "and 20% periclase.  More minerals may be added by simply extending\n",
    "the list given to burnman.composite"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rock = burnman.Composite([minerals.SLB_2011.mg_perovskite(), minerals.SLB_2011.periclase()], \\\n",
    "                         [0.8, 0.2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point we want to tell the rock which equation of state to use for\n",
    "its thermoelastic calculations. In general, we recommend the 'slb3'\n",
    "equation of state as the most self-consistent model.  The parameters from\n",
    "the SLB_2011 mineral library are fit using this model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rock.set_method('slb3')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Input temperature"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we get an array of temperatures at which will be used for computing\n",
    "the seismic properties of the rock.  Here we use the Brown+Shankland (1981)\n",
    "geotherm for mapping pressure to temperature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "temperature = burnman.geotherm.brown_shankland(depths)\n",
    "plt.plot(pressure/1.e9,temperature,'r')\n",
    "plt.xlim(min(pressure)/1.e9,max(pressure)/1.e9)\n",
    "plt.xlabel('pressure (GPa)')\n",
    "plt.ylabel('temperature (K)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Calculate velocities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    \n",
    "Here is the step which does the heavy lifting.  burnman.velocities_from_rock sets the state of the rock at each of the pressures and temperatures defined,then calculates the elastic moduli and density of each individual phase.  After that,it performs elastic averaging on the phases to get a single bulk and shear\n",
    "modulus for the rock.  This averaging scheme defaults to Voigt-Reuss-Hill,\n",
    "but see example_averaging.py for other options.  Finally, it calculates the seismic\n",
    "wave speeds for the whole rock.  It returns a tuple of density, p-wave velocity\n",
    "s-wave velocity, bulk sound speed, bulk modulus, and shear modulus."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "density, vp, vs, vphi, K, G = rock.evaluate(['density','v_p','v_s','v_phi','K_S','G'], pressure, temperature)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Plot results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All the work is done except the plotting!  Here we want to plot the seismic wave\n",
    "speeds and the density against PREM using the matplotlib plotting tools.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First, we plot the s-wave speed verses the PREM s-wave speed\n",
    "plt.figure(figsize=(15,5))\n",
    "plt.subplot(1,3,1)\n",
    "plt.plot(pressure/1.e9,vs/1.e3,color='b',linestyle='-',marker='o', markerfacecolor='b',markersize=4,label='computation')\n",
    "plt.plot(pressure/1.e9,seis_vs/1.e3,color='k',linestyle='-',marker='o', markerfacecolor='k',markersize=4,label='reference')\n",
    "plt.title(\"S wave speed (km/s)\")\n",
    "plt.xlim(min(pressure)/1.e9,max(pressure)/1.e9)\n",
    "plt.xlabel('pressure (GPa)')\n",
    "plt.legend(loc='lower right')\n",
    "plt.ylim(5,8.0)\n",
    "  \n",
    "\n",
    "# Next, we plot the p-wave speed verses the PREM p-wave speed\n",
    "plt.subplot(1,3,2)\n",
    "plt.plot(pressure/1.e9,vp/1.e3,color='b',linestyle='-',marker='o',markerfacecolor='b',markersize=4)\n",
    "plt.plot(pressure/1.e9,seis_vp/1.e3,color='k',linestyle='-',marker='o',markerfacecolor='k',markersize=4)\n",
    "plt.title(\"P wave speed (km/s)\")\n",
    "plt.xlabel('pressure (GPa)')\n",
    "plt.xlim(min(pressure)/1.e9,max(pressure)/1.e9)\n",
    "plt.ylim(10,16)\n",
    "    \n",
    "# Next, we plot the density verses the PREM density\n",
    "plt.subplot(1,3,3)\n",
    "plt.plot(pressure/1.e9,density/1.e3,color='b',linestyle='-',marker='o', markerfacecolor='b',markersize=4)\n",
    "plt.plot(pressure/1.e9,seis_rho/1.e3,color='k',linestyle='-',marker='o', markerfacecolor='k',markersize=4)\n",
    "plt.xlim(min(pressure)/1.e9,max(pressure)/1.e9)\n",
    "plt.xlabel(\"Pressure (GPa)\")\n",
    "\n",
    "plt.title(\"Density ($10^3$ kg/m$^3$)\");"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}