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

# # 2 - Reflexão e refração: Lei de Snell
# 
# Como é o ângulo de reflexão e refração das ondas P e S para diferentes ângulos de incidência? E como isso varia dependendo das velocidades de onda sísmica do meio? A Lei de Snell explica:
# 
# $$\dfrac{\sin \theta_1}{V_1} = \dfrac{\sin \theta_2}{V_2}$$
# 
# Nessa prática, vocês deverão usar a Lei de Snell para prever e explicar o comportamento das ondas P e S ao se depararem com uma interface. 
# 
# 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
# 
# * Entender como funciona a Lei de Snell
# * Observar as previsões da Lei de Snell na propagação de ondas P e S
# * Checar se os valores previstos na Lei de Snell são observados nas simulações

# ## Questão para entregar
# 
# <div class="alert alert-info" style="font-size:12pt; margin-top:20px"> 
# Utilize a Lei de Snell para explicar como serão os ângulos de reflexão e refração das ondas P e S quando uma onda P incide em uma interface. Sua resposta deve explicar como isso varia dependendo das velocidades das ondas sísmicas e possíveis casos extremos e limites.
# </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[1]:


get_ipython().run_line_magic('matplotlib', 'inline')
from __future__ import division, print_function
import numpy as np
import matplotlib.pyplot as plt
import ipywidgets as ipw
from fatiando.seismic import RickerWavelet, GaussianWavelet, FDElasticPSV, FDAcoustic2D


# ## Incidência de ondas P em meios fluidos: V1 < V2
# 
# Vamos simular uma onda P incidindo na interface entre dois meios fluidos com velocidades aumentando com a profundidade.

# In[2]:


shape = (200, 600)
spacing = 5
extent = [0, shape[1]*spacing, shape[0]*spacing, 0]


# In[3]:


vp1 = 1500
vp2 = 2000


# In[4]:


increase_density = np.zeros(shape, dtype='float32') + 1000
increase_velocity = np.zeros(shape, dtype='float32') + vp1
increase_density[100:,:] = 1500
increase_velocity[100:,:] = vp2


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

# In[5]:


p_increase = FDAcoustic2D(increase_velocity, increase_density, spacing=spacing)
fonte = (0, shape[1]//4)
p_increase.add_point_source(fonte, RickerWavelet(1, 60))


# Agora que temos nosso simulador pronto, rode a célcula abaixo para avançar a simulação no tempo.

# In[6]:


p_increase.run(1500)


# Essa simulação pode demorar um pouco para terminar.

# In[7]:


p_increase.animate(every=20, embed=True, dpi=60, cutoff=0.4, fps=7)


# Rode a próxima célula para explorar fotos de cada etapa da simulação, uma de cada vez. A figura também mostrará 3 raios (onda incidente, refletida e refratada) para um determinado ângulo de incidência que você poderá controlar. Os ângulos de refração e reflexão são calculados com a Lei de Snell.

# In[8]:


def plot_with_p_rays_increasing(tempo, incidencia):    
    fig = plt.figure()
    ax = plt.subplot(111)
    p_increase.snapshot(frame=tempo, ax=ax, cutoff=0.2, cmap='Greys')
    fig.set_size_inches(14, 5.5)
    y_bottom = shape[0]*spacing
    y_interface = 100*spacing
    y_source = fonte[0]*spacing
    x_source = fonte[1]*spacing
    x_incidence = (np.tan(np.radians(incidencia))*(y_interface - y_source)
                   + x_source)
    x_reflect = 2*x_incidence - x_source
    arg = (vp2/vp1)*np.sin(np.radians(incidencia))
    if arg <= 1:
        refract = np.arcsin(arg)
        x_refract = (np.tan(refract)*(y_bottom - y_interface) + x_incidence)
        ax.plot([x_incidence, x_refract], [y_interface, y_bottom], '-r', linewidth=2)
    ax.plot([x_source, x_incidence], [y_source, y_interface], '-k', linewidth=2)
    ax.plot([x_incidence, x_reflect], [y_interface, 0], '-b', linewidth=2)
    ax.hlines(y_interface, 0, shape[1]*spacing, colors='grey')
    
ipw.interactive(plot_with_p_rays_increasing, 
                tempo=ipw.IntSlider(value=0, min=0, max=p_increase.it, step=50),
                incidencia=ipw.FloatSlider(value=45, min=0, max=90, step=0.5))


# ### Para pensar
# 
# * Como será a proporção entre o ângulo de incidência e o ângulo de reflexão?
# * E com o ângulo de refração? 
# * Calcule os valores com base na Lei de Snell e cheque com o que você vê no gráfico.
# * Como esses ângulos se relacionam com as curvaturas das frentes de onda?
# * Por que no final da simulação a onda refratada se separa da onda refletida e incidente?
# * Por que não há refração acima de um determinado ângulo?
# * Qual ângulo é esse?

# ## Incidência de ondas P em meios fluidos: V1 > V2
# 
# Agora vamos ver o que acontece quando a velocidade diminui com a profundidade.

# In[9]:


vp1_diminui = 1500
vp2_diminui = 1000


# In[10]:


decrease_density = np.zeros(shape, dtype='float32') + 1500
decrease_velocity = np.zeros(shape, dtype='float32') + vp1_diminui
decrease_density[100:,:] = 2000
decrease_velocity[100:,:] = vp2_diminui


# Novamente, crie a simulação com uma fonte explosiva no topo e a avance no tempo.

# In[11]:


p_decrease = FDAcoustic2D(decrease_velocity, decrease_density, spacing=spacing)
p_decrease.add_point_source(fonte, RickerWavelet(1, 60))


# In[12]:


p_decrease.run(1500)


# In[13]:


p_decrease.animate(every=20, embed=True, dpi=60, cutoff=0.4, fps=7)


# Rode abaixo para gerar a mesma figura interativa da simulação anterior.

# In[14]:


def plot_with_p_rays_decreasing(tempo, incidencia):  
    fig = plt.figure()
    ax = plt.subplot(111)
    p_decrease.snapshot(frame=tempo, ax=ax, cutoff=0.4, cmap='Greys')
    fig.set_size_inches(14, 5.5)
    y_bottom = shape[0]*spacing
    y_interface = 100*spacing
    y_source = fonte[0]*spacing
    x_source = fonte[1]*spacing
    x_incidence = (np.tan(np.radians(incidencia))*(y_interface - y_source)
                   + x_source)
    x_reflect = 2*x_incidence - x_source
    arg = (vp2_diminui/vp1_diminui)*np.sin(np.radians(incidencia))
    if arg <= 1:
        refract = np.arcsin(arg)
        x_refract = (np.tan(refract)*(y_bottom - y_interface) + x_incidence)
        ax.plot([x_incidence, x_refract], [y_interface, y_bottom], '-r', linewidth=2)
    ax.plot([x_source, x_incidence], [y_source, y_interface], '-k', linewidth=2)
    ax.plot([x_incidence, x_reflect], [y_interface, 0], '-b', linewidth=2)
    ax.hlines(y_interface, 0, shape[1]*spacing, colors='grey')
    
ipw.interactive(plot_with_p_rays_decreasing, 
                tempo=ipw.IntSlider(value=0, min=0, max=p_decrease.it, step=50),
                incidencia=ipw.FloatSlider(value=45, min=0, max=90, step=0.5))


# ### Para pensar
# 
# * Como será a proporção entre os ângulos de incidência, reflexão e refração?
# * Calcule os valores com base na Lei de Snell e cheque com o que você vê no gráfico.
# * O que mudou em relação a simulação anterior? Isso é condizente com a Lei de Snell?

# ## Incidência de ondas P em meios sólidos
# 
# Vamos fazer agora uma simulação que inclui tanto ondas P quanto ondas S propagando em um meio sólido. A onda P irá passar para um meio com velocidades maiores.

# In[15]:


vp1_solido = 3000
vp2_solido = 4000
vs1_solido = 2000
vs2_solido = 3000


# In[16]:


shape2 = (300, 600)
l = 150
twosolid_density = np.zeros(shape2, dtype='float32') + 1800
twosolid_density[l:, :] = 2100
twosolid_vs = np.zeros(shape2, dtype='float32') + vs1_solido
twosolid_vs[l:, :] = vs2_solido
twosolid_vp = np.zeros(shape2, dtype='float32') + vp1_solido
twosolid_vp[l:, :] = vp2_solido


# Rode as células abaixo para criar e rodar a simulação utilizando uma fonte explosiva na superfície.

# In[17]:


fonte_explosiva = (70, shape[1]//4)


# In[18]:


ps_twosolid = FDElasticPSV(twosolid_vp, twosolid_vs, twosolid_density, spacing=spacing)
ps_twosolid.add_blast_source(fonte_explosiva, GaussianWavelet(1, 100))


# In[19]:


ps_twosolid.run(1200)


# In[20]:


ps_twosolid.animate(every=20, plottype=['vectors', 'wavefield'], 
                    cutoff=1e-6, scale=1e7, every_particle=10,
                    dpi=60, fps=6, embed=True)


# Novamente, use a célula abaixo para ver cada etapa da simulação separadamente e os raios correspondentes as reflexões e refrações.

# In[21]:


def plot_with_ps_rays(tempo, incidencia):    
    fig = plt.figure()
    ax = plt.subplot(111)
    ps_twosolid.snapshot(frame=tempo, ax=ax, plottype=['wavefield'], 
                         cutoff=1e-6, every_particle=10, cmap='Greys')
    fig.set_size_inches(14, 6.5)
    y_bottom = shape2[0]*spacing
    y_interface = l*spacing
    y_source = fonte_explosiva[0]*spacing
    x_source = fonte_explosiva[1]*spacing
    x_incidence = (np.tan(np.radians(incidencia))*(y_interface - y_source)
                   + x_source)    
    ax.plot([x_source, x_incidence], [y_source, y_interface], '-k', linewidth=2)
    # P reflection
    x_preflect = (np.tan(np.radians(incidencia))*(y_interface - 0) + x_incidence)
    ax.plot([x_incidence, x_preflect], [y_interface, 0], '-b', linewidth=2)
    # S reflection
    sreflect = np.arcsin((vs1_solido/vp1_solido)*np.sin(np.radians(incidencia)))
    x_sreflect = (np.tan(sreflect)*(y_interface - 0) + x_incidence)
    ax.plot([x_incidence, x_sreflect], [y_interface, 0], '--b', linewidth=2)
    # P refraction
    arg = (vp2_solido/vp1_solido)*np.sin(np.radians(incidencia))
    if arg <= 1:
        prefract = np.arcsin(arg)
        x_prefract = (np.tan(prefract)*(y_bottom - y_interface) + x_incidence)
        ax.plot([x_incidence, x_prefract], [y_interface, y_bottom], '-r', linewidth=2)
    # S refraction
    arg = (vs2_solido/vp1_solido)*np.sin(np.radians(incidencia))
    if arg <= 1:
        srefract = np.arcsin(arg)
        x_srefract = (np.tan(srefract)*(y_bottom - y_interface) + x_incidence)
        ax.plot([x_incidence, x_srefract], [y_interface, y_bottom], '--r', linewidth=2)        
    ax.hlines(y_interface, 0, shape[1]*spacing, colors='grey')
    
ipw.interactive(plot_with_ps_rays, 
                tempo=ipw.IntSlider(value=0, min=0, max=ps_twosolid.it, step=50),
                incidencia=ipw.FloatSlider(value=45, min=0, max=90, step=0.5))


# ### Para pensar
# 
# * Como é a proporção do ângulo de reflexão das ondas P e S?
# * E o ângulo de refração?
# * Como seria se a velocidade diminuísse com a profunidade?

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