################################################################################
#                         TAREA 1 METODOS NUMERICOS                            #
################################################################################

'''
Metodos Numericos para la Ciencia e Ingenieria FI3104
Nicolas Troncoso Kurtovic
'''

import matplotlib.pyplot as p  #modulo para graficar
import numpy   as np   #este modulo es para trabajar con matrices como en matlab
import math    as m
from   astropy import constants as cte
from   astropy import units     as units
from   scipy   import integrate as sci

################################################################################

# Extraer datos en dos arreglos

a=np.loadtxt('sun_AM0.dat')  # Se abre el archivo

z = np.loadtxt('sun_AM0.dat')
x = z[:,0]          # Obtiene x de z. Arreglo de lambda en nm
y = z[:,1]          # Obtiene y de z. Arreglo de flujo en W/(nm**2*s)

T=5778.0*units.K    # Temperatura del sol

################################################################################
#                               FLUJO RECIBIDO                                 #
################################################################################
# Calculo del flujo utilizando rectangulos de altura 
# Promedio entre los dos puntos contiguos

flujo=0.0
for i in range(len(x)-1):
  flujo+=((y[i]+y[i+1])/2.0)*(x[i+1]-x[i])

flujo = flujo* units.J/(units.m**2*units.s)
print 'Flujo Recibido SI  = ',flujo

print 'Flujo Recibido cgs = ',flujo.cgs

print '\n'

################################################################################
#                                 LUMINOSIDAD                                  #
################################################################################

lum=4.0*np.pi*(cte.au**2)*flujo

print 'Luminosidad SI  = ',lum

print 'Luminosidad cgs = ',lum.cgs

print '\n'

################################################################################
#                                FLUJO EMITIDO                                 #
################################################################################
# Energia total por unidad de area emitida por el sol

def func(x):
  return m.tan(x)**3/(m.exp(m.tan(x))-1)*(1.0/m.cos(x)**2)

tol = 0.001             # Tolerancia error en los bordes
pas = 0.001             # PASO
ang = np.arange(tol, m.pi/2.0-tol, pas)  # Angulos

integ=0.0
for i in range(len(ang)-1):
  integ += ((func(ang[i])+func(ang[i+1]))/2.0)*(ang[i+1]-ang[i])

# This flux is just the integral of x3/(e(x)-1)
cc    = (2.0*m.pi*cte.h/(cte.c**2))*((cte.k_B*T/cte.h)**4)
flux = integ*cc

print 'Flujo emitido SI  = ',flux

print 'Flujo emitido cgs = ',flux.cgs

print '\n'

print 'pi^4/15 = ',m.pi**4/15.0

print 'pi^4/15 integrado = ',integ

################################################################################
#                                RADIO DEL SOL                                 #
################################################################################
# Estimando el radio del sol

Rsol = np.sqrt( lum / (4.0*np.pi*flux) )

Rsun = cte.R_sun

print 'Radio Sol calculado = ', Rsol

print 'Radio Sol real = ',1*Rsun

print 'Radio Sol cgs = ',1*Rsun.cgs

print '\n'

################################################################################
#                           INTEGRALES USANDO SCIPY                            #
################################################################################
#scipy.integrate.trapz y scipy.integrate.quad

fu = lambda x: m.tan(x)**3/(m.exp(m.tan(x))-1)*(1.0/m.cos(x)**2)  #x va a ser una variable

inttrap_sci = sci.trapz(y, x=x) * units.J/(units.m**2*units.s)
#integral = sci.trapz(variable y, x=variable x) 

intquad_sci = cc * sci.quad(fu, 0.01, np.pi/2.0-0.01)[0]     
#integral = sci.quad(funcion, punto inicial, punto final)[0]
#usar paso 0.001 produce error en esta funcion


print 'Flujo Recibido scipy trapz = ', inttrap_sci

print 'Flujo Emitido scipy quad = ', intquad_sci

################################################################################
#                                  GRAFICOS                                    #
################################################################################
# Grafico loglog

logx=np.zeros(len(x))
logy=np.zeros(len(y))

for i in range(len(x)):
  logx[i]=m.log(x[i])/m.log(10.0)
  logy[i]=m.log(y[i])/m.log(10.0)

x = x*10     # Transformar nm a Angstrom

y = y*10**3

# Plot sin cambios en los ejes

p.plot(x,y)
p.axis([2000, 15000, -0.2, 2500])
p.xlabel("longitud de onda [A]")
p.ylabel("Flujo [erg cm^-2 s^-1]")
p.show()

# Plot en loglog

p.plot(logx,logy)
p.axis([2.0, 4.5, -7, 1])
p.xlabel("log(longitud de onda) [A]")
p.ylabel("log(Flujo) [erg cm^-2 s^-1]")
p.show()


