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

# # 3 - Refração crítica
# 
# Hoje vamos estudar o fenômeno da refração no ângulo crítico (ângulo de incidência que gera refração a 90°). Essa onda possui propriedades distintas que possibilitam seu uso para inferir propriedades da subsuperfície (o método da sísmica de refração). Veremos quais são essas propriedades e possíveis fatores limitantes para a aplicação do método.
# 
# 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
# 
# * Observar como a refração crítica ocorre durante a propagação de ondas P.
# * Verificar as condições para que haja a refração crítica.
# * Treinar como calcular o ângulo crítico e a distância mínima para haver refração crítica.

# ## Questão para entregar
# 
# <div class="alert alert-info" style="font-size:12pt; margin-top:20px">
# <b>Explique o fenômeno da refração no ângulo crítico.</b>
# <p>
# Sua resposta deve conter no mínimo:
# </p>
# <ul>
#     <li>As condições para ela existir.</li>
#     <li>Condições para ela ser observada da superfície.</li>
#     <li>O caminho que a onda percorre da fonte à superfície.</li>
#     <li>Como é o tempo de chegada da onda na superfície e como ele compara com outras ondas.</li>
#     <li>Como calcular o tempo de chegada.</li>
# </ul>
# </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, FDAcoustic2D


# ## Onda P gerada por uma fonte na superfície
# 
# Vamos simular uma onda P incidindo na interface entre dois meios gerada por uma explosão na superfície. 
# 
# Dessa vez, vamos ignorar as ondas S para não complicar a vida.

# In[2]:


shape = (250, 900)
spacing = 30
extent = [0, shape[1]*spacing, shape[0]*spacing, 0]


# In[3]:


vp1 = 3000
vp2 = 4000
print("Vp1 = {} m/s".format(vp1))
print("Vp2 = {} m/s".format(vp2))


# In[4]:


density = np.zeros(shape, dtype='float32') + 2000
velocity = np.zeros(shape, dtype='float32') + vp1
interface = 80
density[interface:,:] = 2300
velocity[interface:,:] = vp2
print("Interface entre os dois meios a {} m de profundidade.".format(interface*spacing))


# A fonte será uma explosão na superfície com as seguintes coordenadas:

# In[5]:


fonte = (0, shape[1]//5)
print('Coordenadas da fonte:')
print('  x = {} m'.format(fonte[1]*spacing))
print('  z = {} m'.format(fonte[0]*spacing))


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

# In[6]:


simul = FDAcoustic2D(velocity, density, spacing=spacing)
simul.add_point_source(fonte, RickerWavelet(1, 30))


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

# In[7]:


simul.run(2000)


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

# In[8]:


simul.animate(every=30, embed=True, dpi=60, cutoff=2, fps=6, cmap='Greys')


# Note que na simulação existem ondas refletindo nas bordas da área. **Ignorem essas ondas**. Elas são artefatos da simulação que são difíceis de evitar.

# Rode a próxima célula para gerar uma figura interativa com:
# 
# * As frentes de onda para cada etapa da simulação, uma de cada vez.
# * Os raios das ondas incidente e refratada para um ângulo de incidência que você poderá controlar.
# * Um marcador (azul) que você poderá controlar a distância e o ângulo (use para achar a coordenada x de algum evento ou seu ângulo de propagação).

# In[9]:


def plot_with_rays(tempo, incidencia, x, angulo):
    fig = plt.figure()
    ax = plt.subplot(111)
    simul.snapshot(frame=tempo, ax=ax, cutoff=2, cmap='Greys')
    fig.set_size_inches(14, 5.4)
    ax.set_title('Tempo em segundos: {} s'.format(tempo*simul.dt))
    y_bottom = shape[0]*spacing
    y_interface = interface*spacing
    y_source = fonte[0]*spacing
    x_source = fonte[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)
    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)    
    tmp = np.tan(np.radians(angulo))
    x_marc_topo = (y_interface - 0)*tmp + x
    x_marc_base = (y_interface - y_bottom)*tmp + x
    ax.plot([x_marc_base, x_marc_topo], [y_bottom, 0], '--b', linewidth=2)
    ax.hlines(y_interface, 0, shape[1]*spacing, colors='grey')
    
ipw.interactive(plot_with_rays, 
                tempo=ipw.IntSlider(value=0, min=0, max=simul.it, step=50),
                incidencia=ipw.FloatSlider(value=45, min=0, max=90, step=0.5),
                x=ipw.FloatSlider(value=10e3, min=0, max=shape[1]*spacing, step=500),
                angulo=ipw.FloatSlider(value=0, min=0, max=90, step=0.5))


# ### Para pensar
# 
# * Use a figura acima para descobrir o caminho que a onda refratada no ângulo crítico percorre até um ponto na superfície. **Dica**: tente achar o caminho que é perpendicular a frente de onda que volta para a superfície.
# * Qual é o ângulo com o qual a onda refratada no ângulo crítico retorna a superfície?
# * Qual é a distância mínima da fonte para que haja refração no ângulo crítico? Tente calcular essa distância?
# * Sabendo a distância acima, qual seria a distância mínima para que a refração seja observada na superfície?
# * Tente observar quando a onda refratada no ângulo crítico passa pelas outras **na superfície**.
# * O que acontece quando a onda refletida na superfície incide na interface com ângulo crítico? 

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