{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Course website**: http://www.leouieda.com/geofisica1\n",
    "\n",
    "**Note**: This notebook is part of the course \"Geofísica 1\" 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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Esse documento que você está usando é 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)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Prática 7 - Magnetometria - Transformada de Fourier e derivadas\n",
    "\n",
    "## Objetivos\n",
    "\n",
    "* Visualizar os efeitos e resultados da Transformada de Fourier em sinais simples.\n",
    "* Aplicar a Transformada de Fourier para calcular derivadas da anomalia magnética de campo total.\n",
    "* Observar os efeitos do erro aleatório nas derivadas calculadas."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Instruções\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": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from __future__ import division\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import ipywidgets as widgets\n",
    "from IPython.display import display\n",
    "import seaborn\n",
    "import fatiando\n",
    "from fatiando import mesher, utils\n",
    "from fatiando.gravmag import prism\n",
    "seaborn.set_context('talk')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print('Versão do Fatiando a Terra: {}'.format(fatiando.__version__))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tranformada de Fourier"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A primeira tarefa servirá para ilustrar o conceito de [Transformada de Fourier](http://en.wikipedia.org/wiki/Fourier_transform) que vimos na aula teórica.\n",
    "\n",
    "Como vimos em aula, uma função $h(t)$ pode ser escrita como:\n",
    "\n",
    "$$\n",
    "h(t) = \\int\\limits_{-\\infty}^{\\infty} H(f)\\ e^{i 2 \\pi f t}df\n",
    "$$\n",
    "\n",
    "em que $f$ é a frequência. $e^{i 2 \\pi f t}$ são \"coisas\" que oscilam com frequência $f$ pois\n",
    "\n",
    "$$\n",
    "e^{i 2 \\pi f t} = \\cos(2 \\pi f t) + i\\sin(2 \\pi f t)\n",
    "$$\n",
    "\n",
    "A equação da transformada nos diz que a função $h(t)$ pode ser escrita como uma \"soma\" infinita de oscilações de frequências diferentes.\n",
    "A oscilação de frequência $f$ possui uma amplitude $H(f)$. A função que descreve as amplitudes, $H(f)$, é a **Transformada de Fourier**. Podemos calcular a transformada $H(f)$ da função utilizando a fórmula abaixo:\n",
    "\n",
    "$$\n",
    "H(f) = \\int\\limits_{-\\infty}^{\\infty} h(t)\\ e^{-i 2 \\pi f t}dt\n",
    "$$\n",
    "\n",
    "$H(f)$ é um número complexo. É comum fazer um gráfico do módulo da transformada $|H(f)|$ pela frequência. Esse gráfico é chamado de **espectro de amplitudes**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Rode a célula abaixo para gerar uma figura interativa.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def dois_senos(f1, A1, f2, A2):\n",
    "    sample = 1/200\n",
    "    t = np.arange(0, 10, sample)\n",
    "    n = t.size\n",
    "    s = A1*np.sin(2*np.pi*f1*t) + A2*np.sin(2*np.pi*f2*t)\n",
    "    f = np.fft.fftfreq(n, sample)[:n//2]\n",
    "    S = np.fft.fft(s)[:n//2]/n\n",
    "    fig = plt.figure(figsize=(12, 6))\n",
    "    ax = plt.subplot(211)\n",
    "    plt.title(r'$h(t) = A_1 \\sin(2\\pi f_1 t) + A_2 \\sin(2\\pi f_2 t)$', fontsize=16)\n",
    "    plt.plot(t, s, '-k')\n",
    "    plt.xlim(0, 5)\n",
    "    plt.ylim(-200, 200)\n",
    "    ax.set_xlabel('tempo (s)', labelpad=-15)\n",
    "    plt.ylabel('Amplitude')\n",
    "    ax = plt.subplot(212)\n",
    "    plt.title(r'$|H(f)|$', fontsize=16)\n",
    "    plt.vlines(f, [0], np.abs(S), colors='k', linewidth=3)\n",
    "    plt.hlines(0, f.min(), f.max(), colors='k', linewidth=5)\n",
    "    plt.xlim(0, 50)\n",
    "    plt.ylim(0, 60)\n",
    "    ax.set_xlabel('frequencia (Hz)', labelpad=-10)\n",
    "    plt.ylabel('Amplitude')\n",
    "    plt.tight_layout()\n",
    "widgets.interactive(dois_senos, \n",
    "                    f1=widgets.IntSlider(min=1, max=40, step=1, value=1), \n",
    "                    A1=widgets.IntSlider(min=0, max=100, step=10, value=100), \n",
    "                    f2=widgets.IntSlider(min=1, max=40, step=1, value=10), \n",
    "                    A2=widgets.IntSlider(min=0, max=100, step=10, value=0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sobre a figura\n",
    "\n",
    "* A figura de **cima** mostra a função $h(t)$ que é composta por dois senos.\n",
    "* Os **botões** controlam a amplitude e frequência de cada seno.\n",
    "* A figura de **baixo** mostra o módulo Transformada de Fourier (espectro de amplitudes) da função $h(t)$. O espectro foi calculado a partir de dados da função $h(t)$ (pontos em uma tabela). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Derivadas\n",
    "\n",
    "Na gravimetria e magnetometria fazemos diversas operações que dependem das derivadas das  anomalias. Precisamos de um jeito de calcular essas derivadas sem sabermos a função que produz a anomalia. Felizmente, as derivadas parciais podem ser calculadas usando a Transformada de Fourier.\n",
    "\n",
    "Sabemos que a transformada de Fourier da derivada $W(f)$ é \n",
    "\n",
    "$$\n",
    "W(f) = H(f)\\ i 2 \\pi f\n",
    "$$\n",
    "\n",
    "Ou seja, podemos calcular a transformada da derivada através da transformada dos dados. Uma vez tendo a transformada da derivada, podemos calcular a derivada fazendo a transformada inversa:\n",
    "\n",
    "\n",
    "$$\n",
    "\\dfrac{\\partial h}{\\partial t} = \\int\\limits_{-\\infty}^{\\infty} W(f)\\ e^{i 2 \\pi f t}df\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Rode a célula abaixo para gerar uma figura interativa.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def dois_senos_derivada(f2, A2):\n",
    "    f1, A1 = 1, 100\n",
    "    sample = 1/200\n",
    "    t = np.arange(0, 10, sample)\n",
    "    n = t.size\n",
    "    h = A1*np.sin(2*np.pi*f1*t) + A2*np.sin(2*np.pi*f2*t)\n",
    "    f = np.fft.fftfreq(n, sample)\n",
    "    H = np.fft.fft(h)\n",
    "    # A linha abaixo calcula a transformada da derivada\n",
    "    W = H*1j*2*np.pi*f\n",
    "    # e essa calcula a derivada através da transformada inversa\n",
    "    w = np.fft.ifft(W)    \n",
    "    f = f[:n//2]\n",
    "    Habs = np.abs(H[:n//2]/n)\n",
    "    Wabs = np.abs(W[:n//2]/n)        \n",
    "    fig, axes = plt.subplots(2, 2, figsize=(14, 7))\n",
    "    # Plot the function and spectrum\n",
    "    ax = axes[0][0]\n",
    "    ax.set_title(r'$h(t) = A_1 \\sin(2\\pi f_1 t) + A_2 \\sin(2\\pi f_2 t)$', fontsize=16)\n",
    "    ax.plot(t, h, '-k')\n",
    "    ax.set_xlim(0, 2.5)\n",
    "    ax.set_ylim(-200, 200)\n",
    "    ax.set_xlabel('t (s)', labelpad=-15)\n",
    "    ax.set_ylabel('Amplitude')\n",
    "    ax = axes[1][0]\n",
    "    ax.set_title(r'$|H(f)|$', fontsize=16)\n",
    "    ax.vlines(f, [0], Habs, colors='k', linewidth=3)\n",
    "    ax.hlines(0, f.min(), f.max(), colors='k', linewidth=5)\n",
    "    ax.set_xlim(0, 50)\n",
    "    ax.set_ylim(0, 60)\n",
    "    ax.set_xlabel('f (Hz)', labelpad=-10)\n",
    "    ax.set_ylabel('Amplitude')\n",
    "    # Plot the derivative and spectrum\n",
    "    ax = axes[0][1]\n",
    "    ax.set_title(r'$\\partial h / \\partial t$', fontsize=16)\n",
    "    ax.plot(t, w, '-k')\n",
    "    ax.set_xlim(0, 2.5)\n",
    "    ax.set_xlabel('t (s)', labelpad=-15)    \n",
    "    ax = axes[1][1]\n",
    "    ax.set_title(r'$|W(f)|$', fontsize=16)\n",
    "    ax.vlines(f, [0], Wabs, colors='k', linewidth=3)\n",
    "    ax.hlines(0, f.min(), f.max(), colors='k', linewidth=5)\n",
    "    ax.set_xlim(0, 50)\n",
    "    ax.set_xlabel('f (Hz)', labelpad=-10)    \n",
    "    plt.tight_layout()\n",
    "widgets.interactive(dois_senos_derivada, \n",
    "                    f2=widgets.IntSlider(min=1, max=40, step=1, value=10), \n",
    "                    A2=widgets.IntSlider(min=0, max=50, step=10, value=0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sobre a figura\n",
    "\n",
    "* A coluna da esquerda na figura acima mostra os dois senos e seu espectro de amplitudes, como na figura anterior.\n",
    "* A coluna da direita mostra a derivada de $h(t)$ e seu espectro de amplitudes. A derivada foi calculada através da transformada de Fourier."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Derivada da anomalia magnética de campo total\n",
    "\n",
    "Vamos utilizar a tranformada de Fourier para calcular a derivada da anomalia de campo total ($\\Delta T$) de um **dique** magnetizado pelo campo da Terra. Nosso dique está **localizado em x = 0** e terá 100 m de espessura e 1 km de profundidade."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Rode a célula abaixo para gerar uma figura interativa.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "areacubo = [-50, 50, 50, 1050]\n",
    "cubo = mesher.Prism(areacubo[0], areacubo[1], -50000, 50000, areacubo[2], areacubo[3])\n",
    "x1, x2 = -2500, 2500\n",
    "xp = np.linspace(x1, x2, 200)\n",
    "yp = np.zeros_like(xp)\n",
    "zp = np.zeros_like(xp)\n",
    "n = xp.size\n",
    "sample = xp[1] - xp[0]\n",
    "f = np.fft.fftfreq(n, sample)\n",
    "\n",
    "def anomalia_derivada_cubo(inc, erro):\n",
    "    cubo.addprop('magnetization', utils.ang2vec(1, inc, 0))\n",
    "    tf = prism.tf(xp, yp, zp, [cubo], inc, 0)       \n",
    "    if erro > 0:\n",
    "        tf = utils.contaminate(tf, erro)\n",
    "    TF = np.fft.fft(tf)\n",
    "    # A linha abaixo calcula a transformada da derivada\n",
    "    W = TF*1j*2*np.pi*f\n",
    "    # e essa calcula a derivada através da transformada inversa\n",
    "    w = np.fft.ifft(W)\n",
    "    espectro_tf = np.abs(TF[:n//2]/n)\n",
    "    espectro_deriv = np.abs(W[:n//2]/n)        \n",
    "    \n",
    "    fig, axes = plt.subplots(2, 2, figsize=(14, 7))\n",
    "    # Plot the function and spectrum\n",
    "    ax = axes[0][0]\n",
    "    ax.set_title(r'$\\Delta T$', fontsize=16)\n",
    "    ax.plot(xp, tf, '-k')\n",
    "    ax.set_xlim(xp.min(), xp.max())\n",
    "    ax.set_xlabel('x (m)')\n",
    "    ax.set_ylabel(r'nT')\n",
    "    ax = axes[1][0]\n",
    "    ax.set_title(r'$|H(f)|$', fontsize=16)\n",
    "    ax.vlines(f[:n//2], [0], espectro_tf, colors='k', linewidth=3)\n",
    "    ax.hlines(0, 0, f.max(), colors='k', linewidth=5)\n",
    "    ax.set_xlim(0, f.max())\n",
    "    ax.set_xlabel('f (1/m)')\n",
    "    ax.set_ylabel('Amplitude')\n",
    "    # Plot the derivative and spectrum\n",
    "    ax = axes[0][1]\n",
    "    ax.set_title(r'$\\partial \\Delta T / \\partial x$', fontsize=16)\n",
    "    ax.plot(xp, w, '-k')\n",
    "    ax.set_xlim(xp.min(), xp.max())\n",
    "    ax.set_xlabel('x (m)')     \n",
    "    ax.set_ylabel(r'nT/m')\n",
    "    ax = axes[1][1]\n",
    "    ax.set_title(r'$|W(f)|$', fontsize=16)\n",
    "    ax.vlines(f[:n//2], [0], espectro_deriv, colors='k', linewidth=3)\n",
    "    ax.hlines(0, 0, f.max(), colors='k', linewidth=5)\n",
    "    ax.set_xlim(0, f.max())\n",
    "    ax.set_xlabel('f (1/m)')   \n",
    "    plt.tight_layout()  \n",
    "widgets.interactive(anomalia_derivada_cubo, \n",
    "                    inc=widgets.IntSlider(min=-90, max=90, step=15, value=-45),\n",
    "                    erro=widgets.IntSlider(min=0, max=20, step=1, value=0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[]()"
   ]
  }
 ],
 "metadata": {
  "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}