In [10]:
#P1 tarea 1
import numpy as np
import matplotlib.pyplot as plt
from pylab import *
datos=np.loadtxt('sun_AM0.dat')
londa=datos[:,0]*0.001 #lambda en micrometros
fsol=datos[:,1]*1000000 #en sistema cgs (W*m-2*nm-1->sist cgs=>*10000000/(10000*(1/1000))=*1000000=*10^6)
fig=plt.figure(1)
fig.set_size_inches(16, 8, forward=True)
plt.xlim(0,3) # xlim(xmin,xmax)
plt.ylim(0,2500000)
plt.xlabel(r'$Longitud \ de \ onda \ (\mu m)$') #r' es para poder usar latex
plt.ylabel(r'$Flujo \ (\frac{erg}{s \ cm^{2} \ \mu m })$')
plt.title(r'$Flujo \ v/s \ longitud \ de \ onda$')
plt.plot(londa,fsol,'r-')
grid(True)
plt.savefig('fl')
show()
#P2 tarea 1
#metodo compuesto por simpson y trapezoides
def intpar(n,x0,xf): #integral numerica simpson para numero de intervalos par y x0 y xf indice valores inicial y final respec. en sun_AM0.dat desde e 0 al 1697
spar=0 #suma pares
simp=0 #suma impares
for i in range(xf-x0+1):
if i%2==0 and i!=0:
spar+=fsol[x0+i]
elif i%2!=0 and i!=0:
simp+=fsol[x0+i]
Ipar=((n*0.001)/3.0)*(fsol[x0]+4*simp+2*spar+fsol[xf])
return Ipar #como esta integral es solo para intervalos par, entonces falta definir otra, para intervalos impar
def intimpar(n,x0,xf): #integral numerica para numero de intervalos impares
finter=0
for i in range(xf-x0-1):
finter+=fsol[x0+1+i]
Iimpar=((n*0.001)/2.0)*(fsol[x0]+2*finter+fsol[xf])
return Iimpar
#Fsolipar intervalos i par, definidos y descritos en informe
Fsol1par=intpar(1,1,510)
Fsol2par=intpar(2,511,1446)
Fsol3par=intpar(20,1446,1571)
Fsol4par=intpar(50,1572,1671)
Fsol5par=intpar(1000,1671,1681)
Fsol6par=intpar(5000,1681,1685)
Fsol7par=intpar(10000,1685,1687)
Fsol8par=intpar(20000,1688,1690)
Fsol9par=intpar(50000,1691,1693)
Fsolintpar=sum(Fsol1par+Fsol2par+Fsol3par+Fsol4par+Fsol5par+Fsol6par+Fsol7par+Fsol8par+Fsol9par)
#Fsoliimpar intervalos i impart, definidos y descritos en informe
Fsol1impar=intimpar(1,0,1)
Fsol2impar=intimpar(1.5,510,511)
Fsol3impar=intimpar(50,1571,1572)
Fsol4impar=intimpar(20000,1687,1688)
Fsol5impar=intimpar(30000,1690,1691)
Fsol6impar=intimpar(50000,1693,1694)
Fsol7impar=intimpar(100000,1694,1695)
Fsol8impar=intimpar(600000,1695,1696)
Fsolintimpar=sum(Fsol1impar+Fsol2impar+Fsol3impar+Fsol4impar+Fsol5impar+Fsol6impar+Fsol7impar+Fsol8impar)
#Fsol: flujo total del sol
Fsol=Fsolintpar+Fsolintimpar #Flujo del sol a la distancia en que se encuentra orbitando la tierra
UA=15000000000000 #unidad astronomica en cm
Lsol=Fsol*4*pi*(UA**2)
print "Flujo del sol a distancia de 1 UA c/datos (algoritmo propio)=",Fsol
print "Luminosidad del sol=",Lsol
#metodo neto de trapezoides
def trap(): #integral numerica para numero de intervalos impares
Itraptot=0
for j in range(1696):
xf=j+1
x0=j
Itrap=(((londa[xf]-londa[x0]))/2.0)*(fsol[x0]+fsol[xf])
Itraptot+=Itrap
return Itraptot
print "Flujo del sol a distancia de 1 UA c/datos (met. del trapecio)=",trap()
#P3 Tarea 1
from astropy import constants as cte
def intB(y): #int de B sin ctes, donde y=arctan(x)=>x->0=>y->0 ^ x->inf=>y->pi/2 => los limites de la nueva integral son y=0 e y=pi/2
return ((tan(y)**3)*((1/cos(y))**2))/(exp(tan(y))-1)
def intsimpson(n): #n:numero de intervalos, debe ser par
a=0
b=pi/2
sumpar=0
sumimpar=0
for i in range(n+1):
if i%2==0 and i!=0 and i!=n:
sumpar+=intB(i*((b-a)/n))
elif i%2!=0 and i!=0 and i!=n:
sumimpar+=intB(i*((b-a)/n))
simps=(((b-a)/n)/3)*(6+2*sumpar+4*sumimpar+0) #simps=((b-a)/n)*(intB(a)+sumpar+sumimpar+int(b)) con intB(a)=6 e intB(b)=0
return simps
I=intsimpson(1000)
T=5778 # °K
c=cte.c.cgs #cte. de veloc. de la luz en cm/s
h=cte.h.cgs #cte. de planck en erg*s
kb=cte.k_B.cgs #cte. de boltzman en erg/K
Fsolteo=((2*pi*h)/(c**2))*(((kb*T)/h)**4)*I #Flujo del sol en su superficie o flujo "emitido"
Reff=(Lsol/(Fsolteo*4*pi))**(0.5)
print "Flujo del sol en la superficie teorico (algoritmo propio)=",Fsolteo #en cgs
print "Radio efectivo del sol=",Reff #en cgs, dado que las ctes no tienen todas unidades, dan unidades distintas a cm, pero el valor corresponde a cm
#P4 Tarea 1
from scipy import integrate as intg
Itrap=intg.trapz(fsol,x=londa) #trapezoides
Y=lambda y:((tan(y)**3)*((1/cos(y))**2))/(exp(tan(y))-1) #definiendo funcion para integrar con intg.quad
Idef=intg.quad(Y,0,pi/2) #integral definida en un intervalo
print "Flujo sol a 1 UA c/datos (trapezoides)=",Itrap #P2
print "Flujo del sol en la superficie teorico (int. definida)=",Idef #P3
mgc=get_ipython().magic #funcion que entrega tiempo de ejecucion
print "Tiempo ejecución Flujo sol a 1 UA c/datos (algoritmo propio trapezoides)=",mgc('%timeit trap()') #P2
print "Tiempo ejecución Flujo sol a 1 UA c/datos (algoritmo propio)=",mgc('%timeit Fsol') #P2
print "Tiempo ejecución Flujo sol a 1 UA c/datos (trapezoides)=",mgc('%timeit Itrap') #P2
print "Tiempo ejecución Flujo sol teorico (algoritmo propio)=",mgc('%timeit Fsolteo') #P3
print "Tiempo ejecución Flujo sol teorico (int. definida)=",mgc('%timeit Idef') #P3
In [ ]: