# -*- coding: iso-8859-15 -*-
# Programa que calcula a forma de uma corda vibrando considerando 
# alguns harmônicos

from matplotlib.pylab import *  # pylab é a maneira fácil de fazer gráficos
import time                     # Vamos animar os gráficos em tempo real
 
ion()   # Ativa o modo interativo (para a animação do gráfico)

# Definições
L_corda=2.5  # Comprimento da corda
Lambda=5.    # Comprimento de onda da frequencia fundamental
freq=0.1     # Frequência
                         
k=2.0*pi/Lambda  # determina o número de onda
w=2.0*pi*freq    # determina a frequência angular

x = arange(0,L_corda,0.01)  # o gráfico terá o comprimento da corda (em x)
corda, = plot(x,x,lw=2, color="black")     # apenas cria o gráfico inicial

# Cria o gráfico de cada harmônico
harmo1, = plot(x,x, color="red") 
harmo2, = plot(x,x, color="blue")
harmo3, = plot(x,x, color="orange")
harmo4, = plot(x,x, color="purple")
harmo5, = plot(x,x, color="green")

corda.axes.set_ylim(-3,3)        # Define os limites do gráfico
 
starttime = time.time()         # O tempo inicial
t = 0                           # Quanto tempo passou
 
while(t < 130.0):            # while(True): para rodar para sempre

    # Calcula por quanto tempo o programa está rondando
    t = time.time() - starttime 

    # Calcula os primeiros harmônicos da corda
    # As amplitudes de cada harmônico dependem da construção do instrumento
    #  e podem ser estimadas pela análise de Fourier do som correspondente
    h1 = 1.0*sin(k*x)*cos(w*t)  # Resultado analítico da onda estacionária
    h2 = 1.0*sin(2.*k*x)*cos(2.*w*t+.00003) # Segundo harmônico
    h3 = 1.0*sin(3.*k*x)*cos(3.*w*t+.00004) # Terceiro
    h4 = 1.0*sin(4.*k*x)*cos(4.*w*t+.00005) # ...
    h5 = 1.0*sin(5.*k*x)*cos(5.*w*t+.00001)
 
    corda.set_ydata(h1+h2+h3+h4+h5)   # A forma da onda resultante

    # apresenta o gráfico de cada harmônico
    harmo1.set_ydata(h1) 
    harmo2.set_ydata(h2)
    harmo3.set_ydata(h3)
    harmo4.set_ydata(h4)
    harmo5.set_ydata(h5)

    draw()                      # Redesenha a tela

    # Aguarda um instante para salvar um pouco o processador 
    time.sleep(0.03)