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

# # 1 - Ondas sísmicas
# 
# Vamos explorar o que são as ondas sísmicas P e S e e como se comportam, tanto em meios homogêneos como quando há uma interface entre dois meios. Para isso, vamos utilizar 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
# 
# * Visualizar como é a propagação de ondas elásticas P e S em meios homogênios
# * Observar o que acontece quando ondas P e S incidem sobre uma interface entre dois meios
# * Entender os conceitos de raios e frentes de onda

# ## Questão para entregar
# 
# <div class="alert alert-info" style="font-size:12pt; margin-top:20px"> 
# Descreva como é a propagação das ondas P e S e os fatores que diferenciam uma da outra. Fale também sobre o que acontece quando uma onda incide em uma interface que separa dois meios.
# </div>
# 
# Lembre-se de colocar **seu nome** e a **data** em sua resposta. A resposta pode ter no **máximo 1 página** (não uma folha).
# 
# **As simulações abaixo foram feitas para te ajudar a responder a questão.**

# ## 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
import numpy as np
import matplotlib.pyplot as plt
from fatiando.seismic import RickerWavelet, GaussianWavelet, FDElasticPSV, FDAcoustic2D


# ## Propagração de uma onda P
# 
# Vamos simular uma onda P (esqueçam as ondas S por enquanto). Primeiro vamos ver como ela se comporta em um meio homogêneo, sem camadas. 
# 
# Rode a célula abaixo para criar um modelo da subsuperfície homogêneo. Para não termos ondas S, vamos considerar que estamos em um fluido, como a água.

# In[ ]:


shape = (300, 400)
spacing = 5
extent = [0, shape[1]*spacing, shape[0]*spacing, 0]


# In[ ]:


velocity = np.zeros(shape, dtype='float32') + 1500  # m/s
density = np.zeros(shape, dtype='float32') + 1000  # kg/m³


# Agora vamos criar o nosso simulador de ondas P com uma fonte explosiva na superfície do nosso modelo.

# In[ ]:


pwave = FDAcoustic2D(velocity, density, spacing=spacing)
pwave.add_point_source((0, shape[1]//2), RickerWavelet(1, 60))


# Agora que temos nossa simulação pronta, rode a célcula abaixo para avançar a simulação do tempo 0 até 400 passos no tempo.

# In[ ]:


pwave.run(400)


# Quando a barra acima encher, rode a célula abaixo para gerar uma animação da propagação da onda. Isso pode demorar um pouco. A parte de cima da figura é a superfície.

# In[ ]:


pwave.animate(every=10, embed=True, dpi=50, cutoff=0.5)


# Rode a próxima célula para explorar fotos de cada etapa da simulação, uma de cada vez.

# In[ ]:


pwave.explore(every=20, cutoff=0.5)


# ### Para pensar
# 
# * Como varia a amplitude da onda ao longo do tempo? 
# * Por que isso acontece?
# * Se nossa simulação fosse num espaço infinito e a rodarmos por um tempo infinito, a onda continuaria existindo e se propagando para sempre? 

# ### Nota
# 
# Se quiser continuar a simulação e avançá-la mais no tempo, rode novamente a célula `pwave.run(400)` para avançar mais 400 iterações. Depois, rode novamente a célula que cria a animação e ela será atualizada.

# ## Onda P incidindo em uma interface (fluido)
# 
# Agora que vocês visualizaram como uma onda P se propaga, vamos ver o que acontece quando a subsuperfície contem duas camadas com velocidades e densidades diferentes. No nosso caso de fluidos, poderiam ser duas camadas de água com temperatuda e salinidade diferentes, fazendo com que a velocidade da onda P e a densidade sejam diferentes.
# 
# Vamos fazer um modelo com a camada de baixo mais densa e com velocidade maior. A interface está localizada a 500m de profundidade.

# In[ ]:


twolayer_density = np.zeros(shape, dtype='float32') + 1000
twolayer_velocity = np.zeros(shape, dtype='float32') + 1500
twolayer_density[100:,:] = 1500
twolayer_velocity[100:,:] = 2000


# Novamente, crie a simulação adicione uma fonte explosiva no topo.

# In[ ]:


pwave2 = FDAcoustic2D(twolayer_velocity, twolayer_density, spacing=spacing, taper=0.005)
pwave2.add_point_source((0, shape[1]//2), RickerWavelet(1, 60))


# Rode a simulação por 600 iterações.

# In[ ]:


pwave2.run(700)


# Gere a animação abaixo. O que acontece quando a frente de onda chega na interface?

# In[ ]:


pwave2.animate(every=10, embed=True, dpi=50, cutoff=0.5)


# Rode abaixo para ver cada tempo separado, com calma.

# In[ ]:


pwave2.explore(every=20, cutoff=0.5)


# ### Para pensar
# 
# * A proporção de energia refletida e refratada é igual?
# * A curvatura da frente de onda refratada é a mesma que a original?
# * E a da refletida?

# ## Ondas P e S em meios sólidos
# 
# Vamos fazer agora uma simulação mais realista que inclui tanto ondas P quanto ondas S propagando em um meio sólido.
# 
# Para isso, precisamos definir qual é a velocidade de cada um dos tipos de onda no meio. Note que a onda P é sempre mais rápida que a S.

# In[ ]:


density_solid = np.zeros(shape, dtype='float32') + 2200
vs_solid = np.zeros(shape, dtype='float32') + 3000
vp_solid = np.zeros(shape, dtype='float32') + 4000


# Rode a célula abaixo para criar a simulação de ondas P e S. Dessa vez, **não vamos usar uma fonte explosiva**. A fonte será uma **falha reversa com mergulho de 45 graus**. O epicentro será localizada em profundidade a 750m.

# In[ ]:


ps_wave = FDElasticPSV(vp_solid, vs_solid, density_solid, spacing=spacing)
ps_wave.add_point_source((shape[0]//2, shape[1]//2), dip=45, wavelet=GaussianWavelet(-1, 100))
ps_wave.add_point_source((shape[0]//2 + 1, shape[1]//2 - 1), dip=45, wavelet=GaussianWavelet(1, 100))


# Agora, avance a simulação 400 iterações no tempo. 

# In[ ]:


ps_wave.run(400)


# E gere a animação.

# In[ ]:


ps_wave.animate(every=10, plottype=['vectors', 'wavefield'], cutoff=5e-7, scale=5e7, 
                dpi=80, fps=6, embed=True)


# Note que essa é diferente da outra. Agora, há vetores na imagem além das frentes de onda coloridas. **Os vetores representam a direção e intensidade do deslocamento de cada ponto no meio**. Os pontinhos pretos são vetores com tamanho zero. Note que no começo da animação aparecem vetorres bem grande no centro representando a fonte da onda.

# Novamente, use a célula abaixo para ver cada etapa da simulação separadamente.

# In[ ]:


ps_wave.explore(every=20, plottype=['vectors', 'wavefield'], scale=5e7, cutoff=5e-7)


# ### Para pensar
# 
# * Tente achar as ondas P e as S na simulação. Dica: use os vetores de deslocamento.
# * Qual frente de onda é mais rápida? 
# * Como é a direção de vibração das ondas P? Ela é uniforme para todas as frentes de onda? Tente explicar esse comportamento.

# ## Onda P incidindo em uma interface (sólido)
# 
# Agora que vocês acharam as ondas P e S, vamos fazer a simulação completa do que acontece quando uma onda P encontra uma interface em um meio sólido.
# 
# Rode a célula abaixo para criar um modelo com dois meios com velocidades diferentes. A interface entre os dois meios está a 600m.

# In[ ]:


density_twosolid = np.zeros(shape, dtype='float32') + 1800
vs_twosolid = np.zeros(shape, dtype='float32') + 2000
vp_twosolid = np.zeros(shape, dtype='float32') + 3000
l = 120
density_twosolid[l:, :] = 2300
vs_twosolid[l:, :] = 3000
vp_twosolid[l:, :] = 4000


# A simulação que vamos fazer agora usa uma **fonte explosiva**.  Esse tipo de fonte gera **somente ondas P**. Vamos ver o que acontece quando essa onda encontra a interface. Para variar um pouco, vamos colocar a fonte no fundo do nosso modelo.

# In[ ]:


ps_twolayer = FDElasticPSV(vp_twosolid, vs_twosolid, density_twosolid, spacing, taper=0.004)
ps_twolayer.add_blast_source((3*shape[0]//5, shape[1]//2), GaussianWavelet(1, 150))


# Agora, avance a simulação.

# In[ ]:


ps_twolayer.run(600)


# E gere a animação. Note que no começo há somente ondas P (pergunta, como você sabe que são só ondas P?).

# In[ ]:


ps_twolayer.animate(every=10, plottype=['vectors', 'wavefield'], every_particle=6, 
                    cutoff=5e-7, scale=3e7, dpi=80, fps=6, embed=True)


# In[ ]:


ps_twolayer.explore(every=20, plottype=['vectors', 'wavefield'], 
                    cutoff=5e-7, scale=3e7, every_particle=6)


# ### Para pensar
# 
# * O que acontece quando a onda P incide na interface?
# * A curvatura das frentes de onda refletidas e refratadas são iguais?
# * Avance a simulação no tempo novamente e atualize a animação. O que acontece quando a onda P reflete na superfície? E quando a onda S reflete?

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