{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5 - Sísmica de reflexão: Reflexão e seções Common Mid Point (CMP)\n",
    "\n",
    "Nessa aula, vamos começar a estudar a sísmica de reflexão. Veremos como identificar a reflexão na simulação de ondas P e no gráfico de tempo x distância. Vamos simular também uma seção Common Mid Point, na qual movemos a fonte e o receptor mantendo o ponto médio fixo.\n",
    "\n",
    "\n",
    "Utilizaremos as simulações de ondas da biblioteca [Fatiando a Terra](http://www.fatiando.org). Essas simulações utilizam o [método de diferenças finitas](http://en.wikipedia.org/wiki/Finite_difference_method) para calcular soluções da equação da onda."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Objetivos\n",
    "\n",
    "* Identificar a reflexão no grafico tempo x distância.\n",
    "* Explicar por que a reflexão aparece como uma hiperbole nesse gráfico.\n",
    "* Identificar a reflexão na seção CMP."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Questão para entregar\n",
    "\n",
    "<div class=\"alert alert-info\" style=\"font-size:12pt; margin-top:20px\">\n",
    "<b>Explique (deduza) por que a reflexão aparece como uma hipérbole na seção CMP e por que utilizamos essa seção ao invés de um único tiro com offset zero.</b>\n",
    "</div>\n",
    "\n",
    "### Regras para a resposta\n",
    "\n",
    "* Coloque **nome, data e o número da prática** em sua resposta. \n",
    "* A resposta pode ter no **máximo 1 página** (não uma folha).\n",
    "* **Execute o notebook** antes de responder. As simulações abaixo foram feitas para te ajudar.\n",
    "* **Pense e organize** sua resposta andtes de começar a escrever."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Instruções\n",
    "\n",
    "Esse documento é um [Jupyter notebook](http://jupyter.org/), um documento interativo que mistura texto (como esse), código (como abaixo), e o resultado de executar o código (números, texto, figuras, videos, etc).\n",
    "\n",
    "O notebook te fornecerá exemplos interativos que trabalham os temas abordados no questionário. Utilize esses exemplos para responder as perguntas.\n",
    "\n",
    "As células com números ao lado, como `In [1]:`, são código [Python](http://python.org/). Algumas dessas células não produzem resultado e servem de preparação para os exemplos interativos. Outras, produzem gráficos interativos. **Você deve executar todas as células, uma de cada vez**, mesmo as que não produzem gráficos.\n",
    "\n",
    "**Para executar uma célula**, clique em cima dela e aperte `Shift + Enter`. O foco (contorno verde ou cinza em torno da célula) deverá passar para a célula abaixo. Para rodá-la, aperte `Shift + Enter` novamente e assim por diante. Você pode executar células de texto que não acontecerá nada."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup\n",
    "\n",
    "Rode as células abaixo para carregar os módulos necessários para essa prática."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from __future__ import division, print_function\n",
    "from multiprocessing import Pool, cpu_count\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.animation import FuncAnimation\n",
    "import ipywidgets as ipw\n",
    "from fatiando.seismic import RickerWavelet, FDAcoustic2D\n",
    "from fatiando.vis import anim_to_html\n",
    "from fatiando.vis.mpl import seismic_image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reflexão no gráfico tempo x distância\n",
    "\n",
    "Vamos utilizar as simulações de onda P para gerar dados sintéticos. Dessa vez também vamos ignorar as ondas S para não complicar a vida mais do que o necessário. \n",
    "\n",
    "O nosso modelo será uma interface plana separando dois meios."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "shape = (300, 1000)\n",
    "spacing = 10\n",
    "extent = [0, shape[1]*spacing, shape[0]*spacing, 0]\n",
    "interface = 100\n",
    "density = np.zeros(shape, dtype='float32') + 2000\n",
    "density[interface:,:] = 2500\n",
    "velocity = np.zeros(shape, dtype='float32') + 4000\n",
    "velocity[interface:,:] = 5000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "shot = FDAcoustic2D(velocity, density, spacing=spacing, taper=0.004, padding=60)\n",
    "fonte = (0, 100)\n",
    "shot.add_point_source(fonte, RickerWavelet(1, 60))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "shot.run(3000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rode a célula abaixo para gerar uma animação da aquisição."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "every = 20\n",
    "frames = shot.simsize//every\n",
    "receptores = np.arange(fonte[1] + 50, shape[1], 15, dtype='int')\n",
    "x = receptores*spacing\n",
    "dados = shot[:, 0, receptores]\n",
    "times = np.linspace(0, shot.dt*shot.simsize, shot.simsize)\n",
    "dados_dummy = np.zeros_like(dados)\n",
    "\n",
    "fig, axes = plt.subplots(2, 1, sharex=True, sharey=False, figsize=(12, 8), facecolor='white')\n",
    "ax1 = axes[0]\n",
    "ax1.set_title('iteration: 0')\n",
    "ax1.set_ylabel('time (s)')\n",
    "cutoff = 0.4\n",
    "section = ax1.imshow(dados_dummy, extent=[x.min(), x.max(), times.max(), times.min()], \n",
    "                     aspect=1400, cmap='Greys', vmin=-cutoff, vmax=cutoff,\n",
    "                     interpolation='nearest')\n",
    "ax1.set_ylim(0, times.max())\n",
    "ax2 = axes[1]\n",
    "ax2.set_xlabel('x (m)')\n",
    "ax2.set_ylabel('depth (m)')\n",
    "cutoff = 1\n",
    "wavefield = ax2.imshow(shot[0, :, :], extent=extent, vmin=-cutoff, vmax=cutoff, cmap='Greys')\n",
    "ax2.plot(fonte[1]*spacing, 0, '*y', markersize=15)\n",
    "ax2.plot(x, np.zeros_like(receptores), 'vb', markersize=10)\n",
    "ax2.hlines(interface*spacing, 0, shape[1]*spacing)\n",
    "fig.tight_layout(h_pad=0)\n",
    "\n",
    "def anim_shot(frame):\n",
    "    ax1.set_title('iteration: {:d}'.format(frame*every))\n",
    "    u = shot[frame*every]\n",
    "    wavefield.set_array(u)\n",
    "    dados_dummy[:frame*every, :] = dados[:frame*every, :]\n",
    "    section.set_array(dados_dummy)\n",
    "    return wavefield\n",
    "\n",
    "anim = FuncAnimation(fig, anim_shot, frames=frames)\n",
    "anim_to_html(anim, fps=6, dpi=60)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Na animação acima:\n",
    "\n",
    "* Os triângulo azuis são os receptores.\n",
    "* A estrela amarela é a fonte.\n",
    "* O painel de baixo mostra a simulação da propagação da onda.\n",
    "* O painel de cima mosta os dados registrados em cada receptor ao longo do tempo. A cor branca significa movimentação para baixo e a cor preta movimentação para cima."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Para pensar\n",
    "\n",
    "* Tente identificar a onda direta, a refletida e a refratada no ângulo crítico na animação.\n",
    "* Repare que a frente de onda refletida se aproxima da direta (no painel de baixo).\n",
    "* Isso é condizente com a fórmula para o tempo de chegada da refletida (o que acontece com a fórmula quando x fica grande)?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulando um CMP\n",
    "\n",
    "Agora vamos simular a aquisição de uma seção Common Mid Point. A diferença é que dessa vez utilizaremos vários pares fonte-receptor a distâncias diferentes. Na simulação anterior, utilizamos apenas uma fonte e vários receptores. Cada fonte será uma simulação da qual vamos extrair as medições de um único receptor. A animação abaixo deixará isso mais claro."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "shape = (100, 200)\n",
    "spacing = 10\n",
    "extent = [0, shape[1]*spacing, shape[0]*spacing, 0]\n",
    "interface = shape[0]//2\n",
    "density = np.zeros(shape, dtype='float32') + 2000\n",
    "density[interface:,:] = 2500\n",
    "velocity = np.zeros(shape, dtype='float32') + 4000\n",
    "velocity[interface:,:] = 5000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "step = 3\n",
    "fontes = np.array(list(reversed(range(55, shape[1]//2 - step//2, step))))\n",
    "recep = np.array([shape[1] - s for s in fontes])\n",
    "offsets = (recep - fontes)*spacing\n",
    "print(\"Utilizando {} fontes e {} receptores.\".format(len(fontes), len(recep)))\n",
    "print('Fontes (m):     {}'.format(fontes))\n",
    "print('Receptores (m): {}'.format(recep*spacing))\n",
    "print('Offsets (m):    {}'.format(offsets))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rode as células abaixo para rodar uma simulação para cada fonte acima. Não aparecerá a barrinha de progresso dessa vez pois vamos rodar as simulações em paralelo para agilizar o processo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def simul_shot(fonte, its=800, verbose=False):\n",
    "    shot = FDAcoustic2D(velocity, density, spacing=spacing, taper=0.005, padding=50, verbose=verbose)\n",
    "    shot.add_point_source((0, fonte), RickerWavelet(1, 100))\n",
    "    shot.run(its)\n",
    "    return shot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "print('Simulando...')\n",
    "pool = Pool(processes=cpu_count())\n",
    "shots = pool.map(simul_shot, fontes)\n",
    "pool.close()\n",
    "print('Terminado.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rode a célula abaixo para gerar uma animação da aquisição do CMP."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "every = 30\n",
    "frames = shots[0].simsize//every\n",
    "dt = shots[0].dt\n",
    "times = np.linspace(0, dt*shots[0].simsize, shots[0].simsize)\n",
    "CMP = np.empty((shots[0].simsize, len(recep)))\n",
    "for i, s in enumerate(shots):\n",
    "    CMP[:, i] = s[:, 0, recep[i]]\n",
    "CMP_dummy = np.zeros_like(CMP)\n",
    "\n",
    "fig, axes = plt.subplots(2, 1, sharex=False, sharey=False, figsize=(7, 8), facecolor='white')\n",
    "ax1 = axes[1]\n",
    "ax1.set_title('shot: 1')\n",
    "ax1.set_xlabel('x (m)')\n",
    "ax1.set_ylabel('depth (m)')\n",
    "cutoff = 1\n",
    "wavefield = ax1.imshow(shots[0][500, :, :], extent=extent, vmin=-cutoff, vmax=cutoff, cmap='Greys')\n",
    "src, = ax1.plot(fontes[0]*spacing, 0, '*y', markersize=15)\n",
    "rec, = ax1.plot(recep[0]*spacing, 0, 'vb', markersize=10)\n",
    "ax1.hlines(interface*spacing, 0, shape[1]*spacing)\n",
    "\n",
    "ax2 = axes[0]\n",
    "ax2.set_title('CMP: 1')\n",
    "ax2.set_xlabel('offset (m)')\n",
    "ax2.set_ylabel('times (s)')\n",
    "cutoff = 1\n",
    "section = ax2.imshow(CMP_dummy, extent=[offsets.min(), offsets.max(), times.min(), times.max()], \n",
    "                     aspect=690, cmap='Greys', vmin=-cutoff, vmax=cutoff, interpolation='nearest')\n",
    "\n",
    "fig.tight_layout()\n",
    "\n",
    "def anim_shot(i):\n",
    "    shot = i//frames\n",
    "    frame = i%frames\n",
    "    ax1.set_title('shot: {:d}'.format(shot + 1))\n",
    "    u = shots[shot][frame*every]\n",
    "    wavefield.set_array(u)\n",
    "    src.set_xdata(fontes[shot]*spacing)\n",
    "    rec.set_xdata(recep[shot]*spacing)\n",
    "    ax2.set_title('CMP: {:d}'.format(shot + 1))\n",
    "    CMP_dummy[:, :shot] = CMP[:, :shot]\n",
    "    CMP_dummy[:frame*every, shot] = CMP[:frame*every, shot]\n",
    "    section.set_array(CMP_dummy)\n",
    "    return wavefield, CMP_dummy\n",
    "\n",
    "anim = FuncAnimation(fig, anim_shot, frames=frames*len(shots))\n",
    "anim_to_html(anim, fps=6, dpi=60)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Na animação acima:\n",
    "\n",
    "* O painel de baixo mostra a simulação de cada tiro.\n",
    "* A estrela representa a fonte e o triângulo representa o receptor.\n",
    "* O painel de cima mostra a seção CMP. Cada coluna nesse gráfico representa os dados medidos por um receptor."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Para pensar\n",
    "\n",
    "* Identifique a onda direta e a reflexão no CMP.\n",
    "* De quais parâmetros depende a curva do tempo de chegada da onda refletida?\n",
    "* O que conseguimos obter com um CMP que não podemos obter com um único par fonte-receptor?\n",
    "* Por que utilizar um CMP ao invés de um único tiro com vários receptores, como na primeira simulação?\n",
    "* O CMP nos dá informação sobre toda a subsuperfície? Ou seja, com 1 único CMP podemos obter seção sísmica que vocês viram na estratigrafia?\n",
    "* É prático fazer a aquisição de CMPs da forma como simulamos (dê um tiro, mova fonte e receptor, dê um tiro, ...)?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## License and information\n",
    "\n",
    "**Course website**: https://github.com/leouieda/geofisica2\n",
    "\n",
    "**Note**: This notebook is part of the course \"Geofísica 2\" of Geology program of the \n",
    "[Universidade do Estado do Rio de Janeiro](http://www.uerj.br/). \n",
    "All content can be freely used and adapted under the terms of the \n",
    "[Creative Commons Attribution 4.0 International License](http://creativecommons.org/licenses/by/4.0/).\n",
    "\n",
    "![Creative Commons License](https://i.creativecommons.org/l/by/4.0/88x31.png)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}