notebooks/temp/pratica3.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:1a3420b9cbcccd54f98b6c03530ed8f6b5eaabb616bc8451cce51a6e1befa6ea"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Course website**: http://www.leouieda.com/geofisica1\n",
      "\n",
      "**Note**: This notebook is part of the course \"Geof\u00edsica 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": [
      "# Pr\u00e1tica 3 - Magnetiza\u00e7\u00e3o remanente\n",
      "\n",
      "Objetivos:\n",
      "\n",
      "* Visualizar os efeitos da magnetiza\u00e7\u00e3o remanente na anomalia magn\u00e9tica de campo total.\n",
      "* Analisar o comportamento da derivada total e da redu\u00e7\u00e3o ao polo na presen\u00e7a de magnetiza\u00e7\u00e3o remanente.\n",
      "* Utilizar a redu\u00e7\u00e3o ao polo como meio de verificar se h\u00e1 magnetiza\u00e7\u00e3o remanente (e possivelmente estimar sua dire\u00e7\u00e3o)."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Prepara\u00e7\u00e3o\n",
      "\n",
      "Esse documento que voc\u00ea est\u00e1 usando \u00e9 um [IPython notebook](http://ipython.org/notebook.html). \u00c9 um documento interativo que mistura texto (como esse), c\u00f3digo (como abaixo), e o resultado de executar o c\u00f3digo (que pode ser n\u00fameros, texto, figuras, videos, etc). Esta pr\u00e1tica usar\u00e1 a biblioteca [Fatiando a Terra](http://fatiando.org) de modelagem geof\u00edsica e tamb\u00e9m o m\u00f3dulo [numpy](http://www.numpy.org/).\n",
      "\n",
      "O notebook \u00e9 divido em c\u00e9lulas (como esta). Para editar o conte\u00fado de uma c\u00e9lula, clique nela (clique nesta para editar esse texto). Para executar uma c\u00e9lula, aperte `Shift + Enter`. Execute as duas c\u00e9lulas abaixo."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "from __future__ import division\n",
      "from IPython.html import widgets\n",
      "import numpy as np\n",
      "from fatiando.gravmag import prism, fourier\n",
      "from fatiando import utils, gridder, mesher\n",
      "from fatiando.vis import mpl\n",
      "import fatiando\n",
      "mpl.rc('lines', linewidth=2)\n",
      "mpl.rc('font', size=12)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print('Vers\u00e3o do Fatiando a Terra: {}'.format(fatiando.__version__))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Tarefa 1: Magnetiza\u00e7\u00e3o induzida e remanente"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "areacubo = [-400, 400, 400, 1200]\n",
      "cubo = mesher.Prism(areacubo[0], areacubo[1], -50000, 50000, areacubo[2], areacubo[3])\n",
      "area = (-1200, 1200, -400, 2000)\n",
      "x, y, z = gridder.regular(area, (13, 13), z=0)\n",
      "def total_cubo(remanente, inc):\n",
      "    def mask(area, v):\n",
      "        x1, x2, y1, y2 = area\n",
      "        v[(x >= x1) & (x <= x2) & (y >= y1) & (y <= y2)] = 0\n",
      "        return v\n",
      "    terra = utils.ang2vec(5000, 45, 0)\n",
      "    remanente = utils.ang2vec(remanente, inc, 0)\n",
      "    induzido = utils.ang2vec(20, 45, 0)\n",
      "    total = induzido + remanente\n",
      "    cubo.addprop('magnetization', total)\n",
      "    bxi = mask(areacubo, prism.bx(x, z, y, [cubo], pmag=induzido))\n",
      "    bzi = mask(areacubo, prism.bz(x, z, y, [cubo], pmag=induzido))\n",
      "    bxr = mask(areacubo, prism.bx(x, z, y, [cubo], pmag=remanente))\n",
      "    bzr = mask(areacubo, prism.bz(x, z, y, [cubo], pmag=remanente))\n",
      "    xp = np.linspace(area[0], area[1], 100)\n",
      "    zp = np.zeros_like(xp)\n",
      "    yp = zp\n",
      "    tf = prism.tf(xp, yp, zp, [cubo], 45, 0)\n",
      "    tx =  bxi + bxr + terra[0]\n",
      "    tz =  bzi + bzr + terra[2]\n",
      "    fig, axes = mpl.subplots(2, 2, sharex='col', figsize=(10, 8))\n",
      "    ax1, ax2, ax3, ax4 = axes.ravel()\n",
      "    for ax in [ax1, ax2, ax3]:\n",
      "        ax.set_aspect('equal')\n",
      "        mpl.sca(ax)\n",
      "        mpl.set_area(area)\n",
      "        ax.invert_yaxis()        \n",
      "    fig.text(0.05, 0.96, '(a)', fontsize=16)        \n",
      "    fig.text(0.52, 0.96, '(b)', fontsize=16)        \n",
      "    fig.text(0.05, 0.51, '(c)', fontsize=16)        \n",
      "    fig.text(0.52, 0.51, '(d)', fontsize=16)\n",
      "    mpl.sca(ax1)\n",
      "    mpl.square(areacubo, fill='grey')\n",
      "    mpl.quiver(x, y, bxi, bzi, linewidth=0.5,\n",
      "               scale=50, pivot='middle', angles='xy', scale_units='xy')\n",
      "    mpl.hlines(0, area[0], area[1], linewidth=3)\n",
      "    mpl.arrow(0, 800, 5*induzido[0], 5*induzido[2], width=2, color='k')\n",
      "    mpl.ylabel('z (m)')\n",
      "    mpl.title('Campo induzido')\n",
      "    mpl.sca(ax2)\n",
      "    mpl.square(areacubo, fill='grey')\n",
      "    mpl.quiver(x, y, bxr, bzr, linewidth=0.5,\n",
      "               scale=50, pivot='middle', angles='xy', scale_units='xy',\n",
      "               color='r')\n",
      "    mpl.hlines(0, area[0], area[1], linewidth=3)\n",
      "    mpl.arrow(0, 800, 5*remanente[0], 5*remanente[2], width=2, color='r')\n",
      "    mpl.title('Campo remanente')\n",
      "    mpl.sca(ax3)\n",
      "    mpl.square(areacubo, fill='grey')\n",
      "    mpl.quiver(x, y, tx, tz, linewidth=0.5, color='b',\n",
      "               scale=50, pivot='middle', angles='xy', scale_units='xy')\n",
      "    mpl.hlines(0, area[0], area[1], linewidth=3)\n",
      "    mpl.arrow(0, 800, 5*total[0], 5*total[2], width=2, color='b')\n",
      "    mpl.xlabel('x (m)')\n",
      "    mpl.ylabel('z (m)')\n",
      "    mpl.title('Campo total = corpo + Terra')\n",
      "    mpl.sca(ax4)\n",
      "    mpl.plot(xp, tf, '-k')\n",
      "    mpl.title('Anomalia de campo total')\n",
      "    mpl.grid(True)\n",
      "    mpl.hlines(0, area[0], area[1], linewidth=1, linestyles='dashed')\n",
      "    mpl.xlabel('x (m)')\n",
      "    mpl.tight_layout()\n",
      "widgets.interactive(total_cubo, \n",
      "                    remanente=widgets.IntSliderWidget(min=0, max=50, step=5, value=0),                    \n",
      "                    inc=widgets.IntSliderWidget(min=-180, max=180, step=5, value=45))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Sobre a figura\n",
      "\n",
      "* Os paineis mostram:\n",
      "    * (a) Campo magn\u00e9tico causado pela magnetiza\u00e7\u00e3o **induzida** no quadrado pelo campo da Terra.\n",
      "    * (b) Campo magn\u00e9tico causado pela magnetiza\u00e7\u00e3o **remanente**.\n",
      "    * (c) Campo magn\u00e9tico total = induzido + remanente + Terra.\n",
      "    * (d) Anomalia magn\u00e9tica de campo total que seria medida na superf\u00edcie (z = 0).\n",
      "* Vetores coloridos no centro dos quadrados mostram a respectiva magnetiza\u00e7\u00e3o: (a) induzida, (b) remanente, (c) total.\n",
      "* **A inclina\u00e7\u00e3o do campo da Terra \u00e9 fixa em** $45^\\circ$.\n",
      "* O bot\u00e3o `remanente` controla a intensidade da magnetiza\u00e7\u00e3o remanente.\n",
      "* O bot\u00e3o `inc` controla a inclina\u00e7\u00e3o da magnetiza\u00e7\u00e3o remanente. Note que ela pode variar de -180 a 180."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Tarefas e perguntas\n",
      "\n",
      "1. Aumente a intensidade da magnetiza\u00e7\u00e3o remanente e mantenha a inclina\u00e7\u00e3o fixa em $45^\\circ$.\n",
      "    * O que acontece com o campo total? A dire\u00e7\u00e3o do campo magn\u00e9tico muda?\n",
      "    * O que acontece com a anomalia de campo total?\n",
      "2. Mude a inclina\u00e7\u00e3o da magnetiza\u00e7\u00e3o remanente dentro do intervalo [-90, 90].\n",
      "    * Como muda a anomalia de campo total? Dica: perceba a propor\u00e7\u00e3o entre a parte positiva e negativa da anomalia.\n",
      "3. Coloque uma magnetiza\u00e7\u00e3o remanente com inc > 90 ou inc < -90.\n",
      "    * Como muda a anomalia de campo total?\n",
      "    * Em que sita\u00e7\u00e3o/situa\u00e7\u00f5es poder\u00edamos encontrar esse tipo de magnetiza\u00e7\u00e3o remanente no mundo real?"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Tarefa 2: Derivada total na presen\u00e7a de magnetiza\u00e7\u00e3o remanente"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "shape = (100, 100)\n",
      "x, y, z = gridder.regular((-5000, 5000, -5000, 5000), shape, z=0)\n",
      "dx, dy = 500, 5000\n",
      "cubo = mesher.Prism(-dy/2, dy/2, -dx/2, dx/2, 400, 4000)\n",
      "cubo_area = cubo.get_bounds()[:4][::-1]\n",
      "def mapa_cubo(remanente, inc, dec):\n",
      "    induzida = utils.ang2vec(10, 45, 0)\n",
      "    remanente = utils.ang2vec(remanente, inc, dec)\n",
      "    tf = prism.tf(x, y, z, [cubo], 45, 0, pmag=induzida + remanente)\n",
      "    dt = fourier.ansig(x, y, tf, shape)\n",
      "    fig, axes = mpl.subplots(1, 2, figsize=(14, 6))\n",
      "    titles = ['Induzida: inc={} dec={} \\n Remanente: int={} inc={} dec={}'.format(45, 0, np.linalg.norm(remanente), inc, dec),\n",
      "              'Derivada total']\n",
      "    for data, ax, title in zip([tf, dt], axes, titles):\n",
      "        ax.set_aspect('equal')\n",
      "        mpl.sca(ax)\n",
      "        mpl.title(title)\n",
      "        #mpl.square(cubo_area)\n",
      "        scale = np.abs([data.min(), data.max()]).max()\n",
      "        mpl.contourf(y, x, data, shape, 30, cmap=mpl.cm.RdBu_r, vmin=-scale, vmax=scale)\n",
      "        mpl.colorbar(pad=0).set_label('nT')\n",
      "        mpl.m2km()\n",
      "        mpl.xlabel('y (km)')\n",
      "        mpl.ylabel('x (km)')\n",
      "    mpl.tight_layout(h_pad=0, w_pad=0)    \n",
      "widgets.interactive(mapa_cubo,\n",
      "                    remanente=widgets.FloatSliderWidget(min=0, max=50, step=1, value=0),\n",
      "                    inc=widgets.IntSliderWidget(min=-180, max=180, step=5, value=45),\n",
      "                    dec=widgets.IntSliderWidget(min=-90, max=90, step=5, value=0))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Sobre a figura\n",
      "\n",
      "* A figura da esquerda mostra a anomalia de campo total produzida por um corpo alongado em subsuperf\u00edcie.\n",
      "* A figura da direita mostra a derivada total da anomalia a esquerda.\n",
      "* Ambos mapas representam a anomalia vista de cima.\n",
      "* A inclina\u00e7\u00e3o e declina\u00e7\u00e3o do campo da Terra s\u00e3o fixos em inc=45 e dec=0.\n",
      "* Sempre h\u00e1 magnetiza\u00e7\u00e3o induzida.\n",
      "* O bot\u00e3o `remanente` controla a intensidade da magnetiza\u00e7\u00e3o remanente.\n",
      "* Os bot\u00f5es `inc` e `dec` controlam, respectivamente, a inclina\u00e7\u00e3o e declina\u00e7\u00e3o da magnetiza\u00e7\u00e3o remanente."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Tarefas e perguntas\n",
      "\n",
      "1. Aumente a intensidade da magnetiza\u00e7\u00e3o remanente e mantenha a inc=45 e dec=0.\n",
      "    * A anomalia de campo total se altera? Por que?\n",
      "    * A derivada total se altera? Por que?\n",
      "2. Varie a inclina\u00e7\u00e3o e declina\u00e7\u00e3o.\n",
      "    * Como a derivada total poderia auxiliar na interpreta\u00e7\u00e3o de anomalias magn\u00e9ticas complexas?\n",
      "    * Em quais situa\u00e7\u00f5es a derivada total tem um pior desempenho em delinear o corpo causador da anomalia?"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Tarefa 3: Redu\u00e7\u00e3o ao polo na presen\u00e7a de magnetiza\u00e7\u00e3o remanente"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "shape = (100, 100)\n",
      "x, y, z = gridder.regular((-5000, 5000, -5000, 5000), shape, z=0)\n",
      "dx, dy = 500, 5000\n",
      "cubo = mesher.Prism(-dy/2, dy/2, -dx/2, dx/2, 400, 4000)\n",
      "cubo_area = cubo.get_bounds()[:4][::-1]\n",
      "def reducao_polo(inc, dec, remanente):\n",
      "    induzida = utils.ang2vec(5, 45, 0)\n",
      "    remanente = utils.ang2vec(remanente, -45, 35)\n",
      "    total = induzida + remanente\n",
      "    tf = prism.tf(x, y, z, [cubo], 45, 0, pmag=total)\n",
      "    rtp = fourier.reduce_to_pole(x, y, tf, shape, 45, 0, inc, dec)\n",
      "    positivo = 100*np.sum(rtp >= 0)/rtp.size\n",
      "    titles = ['Induzida: inc={} dec={} \\n Remanente: int={}'.format(45, 0, np.linalg.norm(remanente)),\n",
      "              u'Reduzida ao polo com inc={} dec={} \\n Propor\u00e7\u00e3o positiva/negativa = {}%'.format(inc, dec, positivo)]\n",
      "    fig, axes = mpl.subplots(1, 2, figsize=(14, 6))\n",
      "    for ax, data, title in zip(axes.ravel(), [tf, rtp], titles):\n",
      "        ax.set_aspect('equal')\n",
      "        mpl.sca(ax)\n",
      "        mpl.title(title)\n",
      "        #mpl.square(cubo_area)\n",
      "        scale = np.abs([data.min(), data.max()]).max()\n",
      "        mpl.contourf(y, x, data, shape, 30, cmap=mpl.cm.RdBu_r, vmin=-scale, vmax=scale)\n",
      "        mpl.colorbar(pad=0).set_label('nT')\n",
      "        mpl.m2km()\n",
      "    mpl.tight_layout(h_pad=0, w_pad=0)\n",
      "    return utils.vec2ang(total)\n",
      "w = widgets.interactive(\n",
      "    reducao_polo, \n",
      "    remanente=widgets.FloatSliderWidget(min=0, max=50, step=5, value=0),\n",
      "    inc=widgets.IntSliderWidget(min=-90, max=90, step=5, value=20),\n",
      "    dec=widgets.IntSliderWidget(min=-90, max=90, step=5, value=-30))\n",
      "w"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Sobre a figura\n",
      "\n",
      "* A figura da esquerda mostra a anomalia de campo total produzida por um corpo alongado em subsuperf\u00edcie.\n",
      "* A figura da direita mostra a anomalia da esquerda reduzida ao polo.\n",
      "* A redu\u00e7\u00e3o foi feita utilizando os valores de inclina\u00e7\u00e3o e declina\u00e7\u00e3o da magnetiza\u00e7\u00e3o mostrados no t\u00edtulo.\n",
      "* O t\u00edtulo da figura da direita  tamb\u00e9m mostra a  propor\u00e7\u00e3o entre a \u00e1rea do mapa com valores positivos e a \u00e1rea com valores negativos.\n",
      "* Ambos mapas representam a anomalia vista de cima.\n",
      "* A inclina\u00e7\u00e3o e declina\u00e7\u00e3o do campo da Terra s\u00e3o fixos em inc=45 e dec=0.\n",
      "* Sempre h\u00e1 magnetiza\u00e7\u00e3o induzida.\n",
      "* O bot\u00e3o `remanente` controla a intensidade da magnetiza\u00e7\u00e3o remanente.\n",
      "* A anomalia da esquerda \u00e9 causada pela magnetiza\u00e7\u00e3o total = induzida + remanente.\n",
      "* A dire\u00e7\u00e3o (inclina\u00e7\u00e3o e declina\u00e7\u00e3o) da magnetiza\u00e7\u00e3o remanente \u00e9 fixa.\n",
      "* Os bot\u00f5es `inc` e `dec` controlam, respectivamente, a inclina\u00e7\u00e3o e declina\u00e7\u00e3o utilizadas na redu\u00e7\u00e3o ao polo.\n",
      "* **Voc\u00ea deve usar a dire\u00e7\u00e3o (inc e dec) da magnetiza\u00e7\u00e3o total para fazer a redu\u00e7\u00e3o ao polo.**"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Tarefas e perguntas\n",
      "\n",
      "1. Mantenha a magnetiza\u00e7\u00e3o remanente = 0. Nesse caso a anomalia magn\u00e9tica \u00e9 causada somente por magnetiza\u00e7\u00e3o induzida.\n",
      "    * Quais s\u00e3o os valores de inclina\u00e7\u00e3o e declina\u00e7\u00e3o que devem ser usados para fazer a redu\u00e7\u00e3o ao polo? Por que?\n",
      "    * O que acontece com a propor\u00e7\u00e3o entre a \u00e1rea positiva e negativa quando voc\u00ea utiliza a dire\u00e7\u00e3o de magnetiza\u00e7\u00e3o (inclina\u00e7\u00e3o e declina\u00e7\u00e3o) correta?\n",
      "2. Aumente a magnetiza\u00e7\u00e3o remanente.\n",
      "    * O que acontece com a redu\u00e7\u00e3o ao polo se utilizar os valores de inc e dec determinados na quest\u00e3o anterior?\n",
      "    * Como fica a propor\u00e7\u00e3o entre a \u00e1rea positiva e negativa conforme a magnetiza\u00e7\u00e3o remanente aumenta?\n",
      "    * Como poderiamos utilizar a  redu\u00e7\u00e3o ao polo para verificar se h\u00e1 magnetiza\u00e7\u00e3o remanente ou n\u00e3o?\n",
      "3. Determine a inclina\u00e7\u00e3o e declina\u00e7\u00e3o da magnetiza\u00e7\u00e3o do corpo. Lembre-se que quando h\u00e1 magnetiza\u00e7\u00e3o remanente a anomalia \u00e9 causada pela magnetiza\u00e7\u00e3o total = induzida + remanente.\n",
      "    * Dica: Utilize as observa\u00e7\u00f5es feitas nas quest\u00f5es 1 e 2.\n",
      "    * Descreva seu m\u00e9todo para determinar a dire\u00e7\u00e3o de magnetiza\u00e7\u00e3o total."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Rode a c\u00e9lula abaixo para conferir se o valor de magnetiza\u00e7\u00e3o que voc\u00ea determinou est\u00e1 correto."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print(u'Magnetiza\u00e7\u00e3o total verdadeira: inc={:.2f} dec={:.2f}'.format(w.result[1], w.result[2]))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}