{
 "metadata": {
  "name": "Faulted_model"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Faulted model"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Looking at ways to build faulted models. We'll start with the prelims."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import matplotlib.pyplot as plt\n",
      "import numpy as np"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Set up the model parameters."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "length = 100.# x range\n",
      "width = 100. # y range\n",
      "depth = 240. # z range"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Define the plotting function."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def plot_slices(data, interpolation='nearest', cmap='Greys', vmin=0., vmax=1.):\n",
      "    plt.figure(figsize=(16,12))\n",
      "    \n",
      "    for j in range(3):\n",
      "        plt.subplot(3,4,j+1)\n",
      "        plt.title('y = ' + str(int(j*(length/3))))\n",
      "        plt.imshow(data[:,:,int(j*(length/3))], interpolation=interpolation, cmap=cmap, vmax=vmax, vmin=vmin)\n",
      "    plt.subplot(3,4,4)\n",
      "    plt.title('y = ' + str(int(length-1)))\n",
      "    plt.imshow(data[:,:,-1], interpolation=interpolation, cmap=cmap, vmax=vmax, vmin=vmin)\n",
      "    \n",
      "    for j in range(3):\n",
      "        plt.subplot(3,4,5+j)\n",
      "        plt.title('x = ' + str(int(j*(width/3))))\n",
      "        plt.imshow(data[:,int(j*(width/3)),:], interpolation=interpolation, cmap=cmap, vmax=vmax, vmin=vmin)\n",
      "    plt.subplot(3,4,8)\n",
      "    plt.title('x = ' + str(int(width-1)))\n",
      "    plt.imshow(data[:,-1,:], interpolation=interpolation, cmap=cmap, vmax=vmax, vmin=vmin)\n",
      "    \n",
      "    for j in range(3):\n",
      "        plt.subplot(3,4,9+j)\n",
      "        plt.title('z = ' + str(1 + int(1 + depth/3 + j*depth/(3*3))))\n",
      "        plt.imshow(data[1 + int(depth/3 + j*depth/(3*3)),:,:], interpolation=interpolation, cmap=cmap, vmax=vmax, vmin=vmin)\n",
      "    plt.subplot(3,4,12)\n",
      "    plt.title('z = ' + str(int(2*depth/3 - 1)))\n",
      "    plt.imshow(data[int(2*depth/3 - 1),:,:], interpolation=interpolation, cmap=cmap, vmax=vmax, vmin=vmin)\n",
      "    \n",
      "    plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Build some layers"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's start with thicknesses. We'll just normalize to a sum of 100% and split the total thickness accordingly. When assigning rocks, I guess we just cycle through them, bottom to top. Let's go straight to 3D. "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# These are oldest > youngest\n",
      "layers = [4,1,1,3,1,2,4]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Start with 4 rocks\n",
      "rocks = [(2800, 2850), (2700, 2750), (2400, 2450), (2300,2350)]\n",
      "\n",
      "def make_ai(rocks):\n",
      "    rx = np.array(rocks)\n",
      "    return rx[:,0] * rx[:,1] / 10e6\n",
      "\n",
      "ai = make_ai(rocks)\n",
      "print ai"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 0.798   0.7425  0.588   0.5405]\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def get_depths(depth, layers):\n",
      "    l = [depth]\n",
      "    for layer in layers:\n",
      "        thickness = depth * float(layer) / sum(layers)\n",
      "        l.append(l[-1] - thickness)\n",
      "    return l[:-1]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "depths = get_depths(depth, layers)\n",
      "for i,d in enumerate(depths):\n",
      "    print ai[i%len(ai)]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.798\n",
        "0.7425\n",
        "0.588\n",
        "0.5405\n",
        "0.798\n",
        "0.7425\n",
        "0.588\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def make_earth(depth, length, width, layers, ai, throw=2):\n",
      "    '''\n",
      "    Given dimensions, a list of AIs, and a list of layers,\n",
      "    the function returns an array of AIs - a model earth.\n",
      "    '''\n",
      "    earth = np.empty((depth,length,width))\n",
      "    \n",
      "    # Make a new set of layers with real thicknesses\n",
      "    # Actually want these as list of bases, bottom first\n",
      "    depths = get_depths(depth, layers)\n",
      "    for i,d in enumerate(depths):\n",
      "        earth[:d] = ai[i%len(ai)]\n",
      "        \n",
      "    # Make the faulted layers (shifted version)\n",
      "    earth_f = np.empty_like(earth)\n",
      "    faulted_layers = layers[:]\n",
      "    faulted_layers[0] += throw\n",
      "    faulted_layers[-1] -= throw\n",
      "    depths = get_depths(depth, faulted_layers)\n",
      "    for i,d in enumerate(depths):\n",
      "        earth_f[:d] = ai[i%len(ai)]\n",
      "        \n",
      "    earth[:,length/2:,:] = earth_f[:,length/2:,:]\n",
      " \n",
      "    return earth"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Maybe there's a function to make a 'shifted' version of a vector?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit make_earth(depth, length, width, layers, ai)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10 loops, best of 3: 34.6 ms per loop\n"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "earth = make_earth(depth, length, width, layers, ai)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot_slices(earth, interpolation='nearest', cmap='jet')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}