#!/usr/bin/env python
# coding: utf-8

# # 5 - Sísmica de reflexão: Reflexão e seções Common Mid Point (CMP)
# 
# 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.
# 
# 
# 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.

# ## Objetivos
# 
# * Identificar a reflexão no grafico tempo x distância.
# * Explicar por que a reflexão aparece como uma hiperbole nesse gráfico.
# * Identificar a reflexão na seção CMP.

# ## Questão para entregar
# 
# <div class="alert alert-info" style="font-size:12pt; margin-top:20px">
# <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>
# </div>
# 
# ### Regras para a resposta
# 
# * Coloque **nome, data e o número da prática** em sua resposta. 
# * A resposta pode ter no **máximo 1 página** (não uma folha).
# * **Execute o notebook** antes de responder. As simulações abaixo foram feitas para te ajudar.
# * **Pense e organize** sua resposta andtes de começar a escrever.

# ## Instruções
# 
# 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).
# 
# O notebook te fornecerá exemplos interativos que trabalham os temas abordados no questionário. Utilize esses exemplos para responder as perguntas.
# 
# 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.
# 
# **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.

# ## Setup
# 
# Rode as células abaixo para carregar os módulos necessários para essa prática.

# In[ ]:


get_ipython().run_line_magic('matplotlib', 'inline')
from __future__ import division, print_function
from multiprocessing import Pool, cpu_count
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import ipywidgets as ipw
from fatiando.seismic import RickerWavelet, FDAcoustic2D
from fatiando.vis import anim_to_html
from fatiando.vis.mpl import seismic_image


# ## Reflexão no gráfico tempo x distância
# 
# 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. 
# 
# O nosso modelo será uma interface plana separando dois meios.

# In[ ]:


shape = (300, 1000)
spacing = 10
extent = [0, shape[1]*spacing, shape[0]*spacing, 0]
interface = 100
density = np.zeros(shape, dtype='float32') + 2000
density[interface:,:] = 2500
velocity = np.zeros(shape, dtype='float32') + 4000
velocity[interface:,:] = 5000


# In[ ]:


shot = FDAcoustic2D(velocity, density, spacing=spacing, taper=0.004, padding=60)
fonte = (0, 100)
shot.add_point_source(fonte, RickerWavelet(1, 60))


# In[ ]:


shot.run(3000)


# Rode a célula abaixo para gerar uma animação da aquisição.

# In[ ]:


every = 20
frames = shot.simsize//every
receptores = np.arange(fonte[1] + 50, shape[1], 15, dtype='int')
x = receptores*spacing
dados = shot[:, 0, receptores]
times = np.linspace(0, shot.dt*shot.simsize, shot.simsize)
dados_dummy = np.zeros_like(dados)

fig, axes = plt.subplots(2, 1, sharex=True, sharey=False, figsize=(12, 8), facecolor='white')
ax1 = axes[0]
ax1.set_title('iteration: 0')
ax1.set_ylabel('time (s)')
cutoff = 0.4
section = ax1.imshow(dados_dummy, extent=[x.min(), x.max(), times.max(), times.min()], 
                     aspect=1400, cmap='Greys', vmin=-cutoff, vmax=cutoff,
                     interpolation='nearest')
ax1.set_ylim(0, times.max())
ax2 = axes[1]
ax2.set_xlabel('x (m)')
ax2.set_ylabel('depth (m)')
cutoff = 1
wavefield = ax2.imshow(shot[0, :, :], extent=extent, vmin=-cutoff, vmax=cutoff, cmap='Greys')
ax2.plot(fonte[1]*spacing, 0, '*y', markersize=15)
ax2.plot(x, np.zeros_like(receptores), 'vb', markersize=10)
ax2.hlines(interface*spacing, 0, shape[1]*spacing)
fig.tight_layout(h_pad=0)

def anim_shot(frame):
    ax1.set_title('iteration: {:d}'.format(frame*every))
    u = shot[frame*every]
    wavefield.set_array(u)
    dados_dummy[:frame*every, :] = dados[:frame*every, :]
    section.set_array(dados_dummy)
    return wavefield

anim = FuncAnimation(fig, anim_shot, frames=frames)
anim_to_html(anim, fps=6, dpi=60)


# Na animação acima:
# 
# * Os triângulo azuis são os receptores.
# * A estrela amarela é a fonte.
# * O painel de baixo mostra a simulação da propagação da onda.
# * 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.

# ### Para pensar
# 
# * Tente identificar a onda direta, a refletida e a refratada no ângulo crítico na animação.
# * Repare que a frente de onda refletida se aproxima da direta (no painel de baixo).
# * Isso é condizente com a fórmula para o tempo de chegada da refletida (o que acontece com a fórmula quando x fica grande)?

# ## Simulando um CMP
# 
# 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.

# In[ ]:


shape = (100, 200)
spacing = 10
extent = [0, shape[1]*spacing, shape[0]*spacing, 0]
interface = shape[0]//2
density = np.zeros(shape, dtype='float32') + 2000
density[interface:,:] = 2500
velocity = np.zeros(shape, dtype='float32') + 4000
velocity[interface:,:] = 5000


# In[ ]:


step = 3
fontes = np.array(list(reversed(range(55, shape[1]//2 - step//2, step))))
recep = np.array([shape[1] - s for s in fontes])
offsets = (recep - fontes)*spacing
print("Utilizando {} fontes e {} receptores.".format(len(fontes), len(recep)))
print('Fontes (m):     {}'.format(fontes))
print('Receptores (m): {}'.format(recep*spacing))
print('Offsets (m):    {}'.format(offsets))


# 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.

# In[ ]:


def simul_shot(fonte, its=800, verbose=False):
    shot = FDAcoustic2D(velocity, density, spacing=spacing, taper=0.005, padding=50, verbose=verbose)
    shot.add_point_source((0, fonte), RickerWavelet(1, 100))
    shot.run(its)
    return shot


# In[ ]:


get_ipython().run_cell_magic('time', '', "print('Simulando...')\npool = Pool(processes=cpu_count())\nshots = pool.map(simul_shot, fontes)\npool.close()\nprint('Terminado.')\n")


# Rode a célula abaixo para gerar uma animação da aquisição do CMP.

# In[ ]:


every = 30
frames = shots[0].simsize//every
dt = shots[0].dt
times = np.linspace(0, dt*shots[0].simsize, shots[0].simsize)
CMP = np.empty((shots[0].simsize, len(recep)))
for i, s in enumerate(shots):
    CMP[:, i] = s[:, 0, recep[i]]
CMP_dummy = np.zeros_like(CMP)

fig, axes = plt.subplots(2, 1, sharex=False, sharey=False, figsize=(7, 8), facecolor='white')
ax1 = axes[1]
ax1.set_title('shot: 1')
ax1.set_xlabel('x (m)')
ax1.set_ylabel('depth (m)')
cutoff = 1
wavefield = ax1.imshow(shots[0][500, :, :], extent=extent, vmin=-cutoff, vmax=cutoff, cmap='Greys')
src, = ax1.plot(fontes[0]*spacing, 0, '*y', markersize=15)
rec, = ax1.plot(recep[0]*spacing, 0, 'vb', markersize=10)
ax1.hlines(interface*spacing, 0, shape[1]*spacing)

ax2 = axes[0]
ax2.set_title('CMP: 1')
ax2.set_xlabel('offset (m)')
ax2.set_ylabel('times (s)')
cutoff = 1
section = ax2.imshow(CMP_dummy, extent=[offsets.min(), offsets.max(), times.min(), times.max()], 
                     aspect=690, cmap='Greys', vmin=-cutoff, vmax=cutoff, interpolation='nearest')

fig.tight_layout()

def anim_shot(i):
    shot = i//frames
    frame = i%frames
    ax1.set_title('shot: {:d}'.format(shot + 1))
    u = shots[shot][frame*every]
    wavefield.set_array(u)
    src.set_xdata(fontes[shot]*spacing)
    rec.set_xdata(recep[shot]*spacing)
    ax2.set_title('CMP: {:d}'.format(shot + 1))
    CMP_dummy[:, :shot] = CMP[:, :shot]
    CMP_dummy[:frame*every, shot] = CMP[:frame*every, shot]
    section.set_array(CMP_dummy)
    return wavefield, CMP_dummy

anim = FuncAnimation(fig, anim_shot, frames=frames*len(shots))
anim_to_html(anim, fps=6, dpi=60)


# Na animação acima:
# 
# * O painel de baixo mostra a simulação de cada tiro.
# * A estrela representa a fonte e o triângulo representa o receptor.
# * O painel de cima mostra a seção CMP. Cada coluna nesse gráfico representa os dados medidos por um receptor.

# ### Para pensar
# 
# * Identifique a onda direta e a reflexão no CMP.
# * De quais parâmetros depende a curva do tempo de chegada da onda refletida?
# * O que conseguimos obter com um CMP que não podemos obter com um único par fonte-receptor?
# * Por que utilizar um CMP ao invés de um único tiro com vários receptores, como na primeira simulação?
# * 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?
# * É prático fazer a aquisição de CMPs da forma como simulamos (dê um tiro, mova fonte e receptor, dê um tiro, ...)?

# ## License and information
# 
# **Course website**: https://github.com/leouieda/geofisica2
# 
# **Note**: This notebook is part of the course "Geofísica 2" of Geology program of the 
# [Universidade do Estado do Rio de Janeiro](http://www.uerj.br/). 
# All content can be freely used and adapted under the terms of the 
# [Creative Commons Attribution 4.0 International License](http://creativecommons.org/licenses/by/4.0/).
# 
# ![Creative Commons License](https://i.creativecommons.org/l/by/4.0/88x31.png)