{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# BurnMan - Compute seismic velocities - advanced 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 an extension the simple case. It gives suggestions on how to change the geotherm and composition, and how to compare to the seismic model with depth instead of pressures. One can choose to run  the different cells for the various options (press shift+enter on the cell of choice). \n",
    "\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",
    "\n",
    "- changing geotherm\n",
    "\n",
    "- changing averaging scheme\n",
    "\n",
    "- computing self-consistent depth to compare to seismology"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Import modules"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import of relevant modules."
   ]
  },
  {
   "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('../..'))\n",
    "import burnman\n",
    "from burnman import minerals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Import seismic model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we load the PREM seismic velocity model for comparison."
   ]
  },
  {
   "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",
    "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()\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Input composition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Beneath are three different methods to input a composition. Choose and run one of the cells. \n",
    "Here we show how to print the minerals in a library (needs improvement) if one wants to choose a different one.\n",
    "\n",
    "Available mineral libraries:\n",
    " \n",
    "- SLB_2005\n",
    "- SLB_2011\n",
    "- Matas_etal_2007\n",
    "- Murakami_etal_2012"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "minlib=minerals.SLB_2011\n",
    "print(dir(minlib))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Composite I.  Fixed minerals"
   ]
  },
  {
   "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 than two minerals can simply be added to the composite."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "amount_of_perovskite=0.8\n",
    "rock = burnman.Composite([minerals.SLB_2011.mg_perovskite(), minerals.SLB_2011.periclase()],\\\n",
    "                         [amount_of_perovskite, 1.0-amount_of_perovskite]  )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Composite II. Including iron"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "amount_perovskite = 0.95\n",
    "fe_num_perovskite = 0.08\n",
    "fe_num_periclase  = 0.2\n",
    "perovskite=minerals.SLB_2011.mg_fe_perovskite()\n",
    "ferropericlase=minerals.SLB_2011.ferropericlase()\n",
    "perovskite.set_composition([1.0 - fe_num_perovskite, fe_num_perovskite, 0.0])\n",
    "ferropericlase.set_composition([1.0 - fe_num_periclase, fe_num_periclase])\n",
    "rock = burnman.Composite([perovskite, ferropericlase], \\\n",
    "                         [amount_perovskite, 1.0-amount_perovskite] )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Set method of composite"
   ]
  },
  {
   "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. Methods available are:\n",
    "\n",
    "- bm2 -- third order method with second order shear modulus (Stixrude & Lithgow-Berteloni 2005)- no temperature effects\n",
    "- bm3 -- third order method (Stixrude & Lithgow-Berteloni 2005) - no temperature effects\n",
    "- slb2 -- third order mehtod with second order shear modulus (Stixrude & Lithgow-Berteloni 2005)\n",
    "- slb3 -- third order method (Stixrude & Lithgow-Berteloni 2005)\n",
    "- mgd2 -- second order method (Matas et al. 2007) non-benchmarked\n",
    "- mgd3 -- third order method (Matas et al. 2007) non-benchmarked\n",
    "\n",
    "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 and we suggest not changing this method when using this library. \n"
   ]
  },
  {
   "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 one can choose between the Brown and Shankland (1981) model or an adiabat."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Geotherm I. Brown & Shankland"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "temperature = burnman.geotherm.brown_shankland(depths)\n",
    "title='Brown&Shankland'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Geotherm II. Adiabat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T0=1900 #K anchor temperature\n",
    "temperature = burnman.geotherm.adiabatic(pressure, T0, rock)\n",
    "title='adiabat'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot temperature\n",
    "plt.plot(pressure/1.e9,temperature,'r')\n",
    "plt.xlim(min(pressure)/1.e9,max(pressure)/1.e9)\n",
    "plt.title(title)\n",
    "plt.xlabel('pressure (GPa)')\n",
    "plt.ylabel('temperature (K)')\n",
    "plt.show();"
   ]
  },
  {
   "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. \n",
    "First we plot the results with pressure and then convert to depths to plot with depth."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Plots I. As a function of pressure"
   ]
  },
  {
   "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=(16,4.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
}