{
 "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 6 - Magnetometria - O IGRF e a anomalia magnética de campo total\n",
    "\n",
    "## Objetivos\n",
    "\n",
    "* Visualizar o campo magnético de referência [IGRF](http://en.wikipedia.org/wiki/International_Geomagnetic_Reference_Field).\n",
    "* Verificar os efeitos da variação secular do campo geomagnético no último século.\n",
    "* Visualizar o campo magnético induzido em um corpo para diferentes inclinações.\n",
    "* Entender a relação entre campo magnético induzido, o campo magnético da Terra (IGRF) e a anomalia magnética de campo total.\n",
    "* Entender o papel da inclinação magnética na forma da anomalia."
   ]
  },
  {
   "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": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from __future__ import division\n",
    "import numpy as np\n",
    "import pandas\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.basemap import Basemap\n",
    "import ipywidgets as widgets\n",
    "from IPython.display import display\n",
    "import seaborn\n",
    "import fatiando\n",
    "from fatiando import gridder, mesher, utils\n",
    "from fatiando.gravmag import prism, sphere\n",
    "from fatiando.vis import mpl"
   ]
  },
  {
   "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": [
    "## International Geomagnetic Reference Field (IGRF)\n",
    "\n",
    "Abaixo, vamos fazer gráficos do IGRF calculado para a Terra toda. \n",
    "Os dados do IGRF foram gerados com o programa [Geomag 7](http://www.ngdc.noaa.gov/geomag/models.shtml).\n",
    "Como o IGRF é um campo magnético, ele possui 3 componentes: inclinação, declinação e intensidade.\n",
    "\n",
    "Rode as células abaixo para gerar mapas dessas 3 componentes para o ano 2000. \n",
    "\n",
    "Extra: para ver o campo de outras datas, mude o código `dado = igrf[igrf.year == 2000]` abaixo para ser o ano que deseja. Os dados disponíveis são de 1900 a 2000 a cada 10 anos."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "igrf = pandas.read_csv('../data/igrf-1900-2010.csv', sep=';')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "bm = Basemap(projection='moll', lon_0=0, resolution='c')\n",
    "shape = (181, 90)\n",
    "x, y = bm(igrf[igrf.year == 1900].longitude.as_matrix().reshape(shape), \n",
    "          igrf[igrf.year == 1900].latitude.as_matrix().reshape(shape))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dado = igrf[igrf.year == 2000]\n",
    "\n",
    "fig, axes = plt.subplots(3, 1, figsize=(14, 16))\n",
    "ax1, ax2, ax3 = axes\n",
    "ax1.set_title(u'Declinação')\n",
    "ax2.set_title(u'Inclinação')\n",
    "ax3.set_title('Intensidade')    \n",
    "plot = bm.contourf(x, y, dado.dec.as_matrix().reshape(shape), 60, vmin=-100, vmax=100, ax=ax1, cmap='RdBu_r')\n",
    "plt.colorbar(plot, orientation='vertical', aspect=30, shrink=0.7, pad=0.01, ax=ax1).set_label('graus')\n",
    "bm.drawcoastlines(ax=ax1)\n",
    "plot = bm.contourf(x, y, dado.inc.as_matrix().reshape(shape), 40, ax=ax2, cmap='RdBu_r')\n",
    "plt.colorbar(plot, orientation='vertical', aspect=30, shrink=0.7, pad=0.01, ax=ax2).set_label('graus')\n",
    "bm.drawcoastlines(ax=ax2)\n",
    "plot = bm.contourf(x, y, dado.int.as_matrix().reshape(shape)*0.001, 40, ax=ax3, cmap='Reds')\n",
    "plt.colorbar(plot, orientation='vertical', aspect=30, shrink=0.7, pad=0.01, ax=ax3).set_label('microTesla')\n",
    "bm.drawcoastlines(ax=ax3)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## O campo magnético induzido em um cubo\n",
    "\n",
    "Vamos ver como varia o campo magnético induzido em cubo para diferentes inclinações do campo magnético indutor.\n",
    "\n",
    "Rode as células abaixo para criar um cubo e produzir uma figura interativa."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "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",
    "xp, yp, zp = gridder.regular(area, (13, 13), z=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def total_cubo(inc, suceptibilidade):\n",
    "    cubo.addprop('magnetization', utils.ang2vec(suceptibilidade, inc, 0))\n",
    "    def mask(area, v):\n",
    "        x1, x2, y1, y2 = area\n",
    "        v[(xp >= x1) & (xp <= x2) & (yp >= y1) & (yp <= y2)] = 0\n",
    "        return v\n",
    "    terra = utils.ang2vec(10000, inc, 0)\n",
    "    bx = mask(areacubo, prism.bx(xp, zp, yp, [cubo]))\n",
    "    bz = mask(areacubo, prism.bz(xp, zp, yp, [cubo]))\n",
    "    tx =  bx + terra[0]\n",
    "    tz =  bz + terra[2]\n",
    "    fig, axes = plt.subplots(1, 2, sharey='all', figsize=(12, 6))\n",
    "    ax1, ax2 = axes\n",
    "    plt.sca(ax1)\n",
    "    ax1.set_aspect('equal')\n",
    "    ax2.set_aspect('equal')\n",
    "    mpl.square(areacubo, fill='grey')\n",
    "    plt.quiver(xp, yp, bx, bz, linewidth=0.5, \n",
    "               scale=50, pivot='middle', angles='xy', scale_units='xy')\n",
    "    mx, my, mz = utils.ang2vec(suceptibilidade, inc, 0)\n",
    "    plt.quiver([0], [800], [mx], [mz], linewidth=0.5, color='r',\n",
    "               scale=1e-1, pivot='middle', angles='xy', scale_units='xy', zorder=1000)\n",
    "    plt.hlines(0, area[0], area[1], linewidth=3)\n",
    "    mpl.set_area(area)\n",
    "    plt.gca().invert_yaxis()\n",
    "    plt.xlabel('x (m)')\n",
    "    plt.ylabel('z (m)')\n",
    "    plt.title('Campo do corpo')    \n",
    "    plt.sca(ax2)\n",
    "    mpl.square(areacubo, fill='grey')\n",
    "    plt.hlines(0, area[0], area[1], linewidth=3, zorder=1)\n",
    "    plt.quiver(xp, yp, tx, tz, linewidth=0.5, color='b',\n",
    "               scale=50, pivot='middle', angles='xy', scale_units='xy')\n",
    "    plt.quiver([0], [800], [mx], [mz], linewidth=0.5, color='r',\n",
    "               scale=1e-1, pivot='middle', angles='xy', scale_units='xy', zorder=1000)\n",
    "    plt.quiver([-1000], [-0], [terra[0]], [terra[-1]], linewidth=0.5, color='y',\n",
    "               scale=50, pivot='middle', angles='xy', scale_units='xy', zorder=1001)\n",
    "    plt.xlabel('x (m)')\n",
    "    plt.title('Campo total = corpo + Terra')    \n",
    "    plt.tight_layout(pad=0)\n",
    "widgets.interactive(total_cubo, inc=widgets.IntSlider(min=-90, max=90, step=5, value=45),\n",
    "                    suceptibilidade=widgets.IntSlider(min=0, max=50, step=10, value=0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sobre a figura\n",
    "\n",
    "* A figura da **esquerda** mostra o cubo (em cinza). A linha preta é a superfície (z = 0). \n",
    "* Na presença de um campo externo indutor, o corpo se magnetiza e gera um campo magnético. Esse campo será mostrado nessa figura em setas pretas.\n",
    "* A figura da **direita** mostra o campo total presente, ou seja, o campo externo indutor (como o campo da Terra, por exemplo) mais o campo do cubo.\n",
    "* A seta amarela na figura da direita mostra o campo da Terra.\n",
    "* As setas vermelhas mostram a magnetização induzida no cubo pelo campo da Terra.\n",
    "* O botão `inc` controla a inclinação do campo magnético indutor (o da Terra). \n",
    "* O botão `suceptibilidade` controla a susceptibilidade magnética do corpo (o quanto ele se magnetiza na presença de um campo indutor).\n",
    "* A **anomalia magnética de campo total** é a diferença entre as **amplitudes** do campo total (azul na figura da direita) e o campo da Terra."
   ]
  }
 ],
 "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
}