# -*- coding: utf-8 -*-

#AGA0503 - Métodos Numéricos em Astronomia
#Professor: Alex Carciofi
#Aluna: Anne Viegas Rathsam
#No USP: 10873535
#Exercício de Programação 5 - EP5
#Problema 3 - Rotina Simpson

############################## PACOTE UTILIZADO ##############################

import math

################### DETERMINAÇÃO DAS CONSTANTES UTILIZADAS ###################

c = 2.99793e10    # velocidade da luz em m/s
k = 1.38065e-16    # constante de Boltzmann em erg/K
h = 6.62607e-27    # constante de Planck em erg.s

######################## DEFINIÇÃO DA FUNÇÃO PRINCIPAL #######################

def main():
	''' Dado um intervalo de comprimento de onda lm1 e lm2, o programa calcula
	a integral da função de corpo negro para esses comprimentos de onda (a uma
	temperatura de 2500 K) utilizando o método de Simpson, até obter uma
	precisão eps dada.
	Utiliza as funções corpo_negro, trapezios e simpson.
	'''
	
	# Leitura dos dados de entrada
	
	lm1 = float(input("Digite o valor do comprimento de onda do extremo"
				   " inferior da integral, em micrometros: ")) * 1e-4
	
	lm2 = float(input("Digite o valor do comprimento de onda do extremo"
				   " superior da integral, em micrometros: ")) * 1e-4
				   
   # Multiplicando por 1e-4 para termos os valores em cm
	
	eps = float(input("Digite a precisão desejada: "))
	
	# Agora chamamos a função simpson para calcular a integral
	
	simp, inc, n = simpson(corpo_negro, lm1, lm2, eps)
	
	print("\nO valor da integral calculada no intervalo dado é", simp, "com"
	   " precisão de", inc)
	print("\nO número de intervalos necessários para atingir essa precisão foi", n)	
	
#################### DEFINIÇÃO DA FUNÇÃO DE CORPO NEGRO #####################	
	
def corpo_negro(lm):
	''' Essa função calcula o valor da função de corpo negro para um 
	comprimento de onda lm de um corpo à temperatura 2500 K.
	'''
	
	numerador = 2*math.pi*h*(c**2)
	exponencial = math.exp((h*c)/(lm*k*2500))  
	denominador = (lm**5) * (exponencial - 1)
	
	B = numerador/denominador # B é a função de corpo negro
	
	return B

####################### DEFINIÇÃO DA FUNÇÃO TRAPEZIOS ########################
	
def trapezios(f, x1, x2, n, trap_ant):
	''' Essa função calcula a integral de f no intervalo [x1, x2] pelo método
	dos trapézios utilizando n intervalos. A váriável trap_ant é a estimativa
	daintegral calculada na chamada de função anterior.
	'''

	h = (x2 - x1)/n # Distância entre os pontos
	
	if n == 1: 
		# Só existeum intervalo, ou seja, existem apenas os dois pontos
		# do extremo do intervalo
		termo = (f(x1) + f(x2))/2
		integral = termo*h
		
	else:
		termo = 0
		xn = x1 + h
		
		# Queremos calcular o valor de f apenas nos pontos novos, então vamos 
		# calcular para x's variando de (n+1) * h em (n+1) * h, pois os 
		# x's + 2h já foram calculados na chamada anterior da função
		
		while xn < x2:
			termo += f(xn)
			xn += 2*h
			
		integral = (termo*h) + (trap_ant/2)
		
	return integral
	
####################### DEFINIÇÃO DA FUNÇÃO TRAPEZIOS ########################
	
def simpson(f, x1, x2, eps):
	''' Essa função calcula a integral de f no intervalo [x1, x2] pelo método
	de Simpson utilizando n intervalos. 
	Utiliza a função trapezios.
	'''
	
	# Iniciando as chamadas da função trapézio com n = 2
	n = 2
	trap_ant = trapezios(f, x1, x2, 1, 0)
	simp_ant = trap_ant
	inc = trap_ant
	
	while inc > eps and n < 1e11: # Limitando n para não entrar em loop infinito
		trap = trapezios(corpo_negro, x1, x2, n, trap_ant)
		simp = 1/3 * ((4*trap) - trap_ant)
		inc = abs((simp - simp_ant)/simp)
		
		trap_ant = trap
		simp_ant = simp
		n = n*2
		
	return simp, inc, n

main()