%matplotlib inline

from __future__ import division

import quantities as pq
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

# Generar un cuadro con versiones de las librerías utilizadas en este notebook
#https://github.com/jrjohansson/version_information
%load_ext version_information
%version_information numpy, matplotlib, quantities, pandas

def B(wl,T):
    '''wl es un array de longitudes de onda con unidades de longitud
    T es una temperatura expresada en Kelvin
    el resultado es un array de r.e. con unidades W/(m**2 * nm * sr)
    '''
    I = 2 * pq.constants.h * (pq.c)**2 / wl**5 *  \
        1 / (np.exp((pq.constants.h*pq.c \
        / (wl*pq.constants.k*T)).simplified)-1)
    return I.rescale(pq.watt/(pq.m**2 * pq.nm *pq.sr))

TCMB = 2.725 *pq.K
lmax = pq.constants.b / TCMB
print lmax.simplified

wl = np.arange(0.1,10,0.1) * pq.mm
I = B(wl,TCMB).rescale(pq.watt/(pq.m**2 * pq.mm * pq.sr))

fig, ax = plt.subplots(figsize=(10, 8))
ax.plot(wl, I*1e7)
ax.set_title('Espectro del cuerpo negro para T = 2.725 K \
\n (Longitudes de onda en mm)')
ax.title.set_fontsize(20)
ax.set_xlabel('Longitud de onda en mm')
ax.xaxis.label.set_fontsize(15)
ax.set_ylabel('Radiancia espectral ($10^{-7} W m^{-2} mm^{-1} sr^{-1}$)')
ax.yaxis.label.set_fontsize(15)
ax.grid()

def B_f(wf,T):
    '''wf es un array de frecuencias con unidades en GHz
    T es una temperatura expresada en Kelvin
    el resultado es un array de radiancias espectrales
    con unidades W/(m**2 x GHz x sr)
    '''
    wf = wf.rescale(pq.hertz)
    I = 2 * pq.constants.h * wf**3 / pq.c**2 * 1 \
        / (np.exp((pq.constants.h*wf \
        / (pq.constants.k*T)).simplified)-1)
    return I.rescale(pq.watt/(pq.m**2 * pq.gigahertz *pq.sr))

B_f(10**13*pq.hertz, 7000*pq.K)

wf = np.arange(0.1,1000,1)* pq.gigahertz
I = B_f(wf,TCMB)

fig, ax = plt.subplots(figsize=(10, 8))
ax.plot(wf, I*10**9) # Escalamos las unidades
ax.set_title('Espectro del cuerpo negro  para T = 2.725K \
 \n (Frecuencias en GHz)')
ax.title.set_fontsize(20)
ax.set_xlabel('Frecuencia ($GHz)$')
ax.xaxis.label.set_fontsize(15)
ax.set_ylabel('Radiancia espectral   $(10^{-9} W m^{-2} GHz^{-1} sr^{-1})$')
ax.yaxis.label.set_fontsize(15)
ax.grid()

# Conversión de 160 GHz en cm**(-1)

frec = 160 * pq.gigahertz
wl = pq.c  / frec
print 'longitud de onda en cm = ', wl.rescale(pq.cm)
k = 1/( wl.rescale(pq.cm))
print u'número de onda en ciclos por cm = ', k

# La conversión se puede hacer de forma más abreviada
# simplemente dividiendo por la velocidad de la luz:
fec_en_wave_number = (frec/pq.c).rescale(1/pq.cm)
print fec_en_wave_number

def B_k(wk,T):
    '''wk es un array de frecuencias con unidades en cm^-1
    T es una temperatura expresada en Kelvin
    el resultado es un array de radiancias espectrales
    con unidades W/(m**2 x cm**(-1) x sr)
    '''
    wk = wk.rescale(1/pq.m)
    I = 2 * pq.constants.h * pq.c**2 * wk**3  * 1 \
        / (np.exp((pq.constants.h*pq.c*wk \
        / (pq.constants.k*T)).simplified)-1)
    return I.rescale(pq.watt/(pq.m**2 * pq.cm**(-1)*pq.sr))

f = 5.35 * pq.cm**(-1)
B_k(f, TCMB)

wk = np.arange(0.1,20,0.1) / pq.cm
I = B_k(wk,TCMB)
I = I*10**7 # Expresamos las unidades en múltiplos de 10^(-7)

fig, ax = plt.subplots(figsize=(10, 8))
ax.plot(wk, I)
ax.set_title('Espectro del cuerpo negro  para T = 2.725K \
 \n (Frecuencias en $cm^{-1}$)')
ax.title.set_fontsize(20)
ax.set_xlabel('Frecuencia (ciclos / cm )')
ax.xaxis.label.set_fontsize(15)
ax.set_ylabel('Radiancia espectral $(10^{-7}W m^{-2} (cm^{-1})^{-1}sr^{-1})$')
ax.yaxis.label.set_fontsize(15)
ax.grid()

Jy = pq.UnitQuantity('jansky', 1e-26*pq.watt/(pq.m**2*pq.hertz), symbol = 'Jy')
MJy = pq.UnitQuantity('megajansky', 10**6 * Jy, symbol = 'MJy')

I = 3.34 * 10**(-18) * pq.W * pq.m**(-2) * pq.hertz**(-1) * pq.sr**(-1)
I.rescale(MJy)

wf = np.arange(0.1,1000,1)* pq.gigahertz
I = B_f(wf,TCMB)
I = I.rescale(MJy*pq.sr**(-1)) # convertimos a MJy * sr**(-1)

wk = (wf/pq.c).rescale(1/pq.cm) # convertimos frecuencias a 1/cm
fig, ax = plt.subplots(figsize=(10, 8))
ax.plot(wk, I)
ax.set_title('Espectro del cuerpo negro  para T = 2.725K \
 \n radiancia espectral en MJy / sr')
ax.title.set_fontsize(20)
ax.set_xlabel('Frecuencia (ciclos/cm)')
ax.xaxis.label.set_fontsize(15)
ax.set_ylabel('Radiancia espectral    ($MJy \, sr^{-1}$)')
ax.yaxis.label.set_fontsize(15)
ax.grid()

#En este caso he descargado el fichero en el directorio ../datos
df = pd.read_table('../datos/firas_monopole_spec_v1.txt', 
    skiprows=18, sep='\s+', header=None, 
    names =['freq', 'I', 'residual', 'uncert', 'poles'])

df[0:10]

wf = np.arange(0.1,1000,1)* pq.gigahertz
I = B_f(wf,TCMB)
I = I.rescale(MJy*pq.sr**(-1)) # convertimos a MJy * sr**(-1)

wk = (wf/pq.c).rescale(1/pq.cm) # convertimos frecuencias a 1/cm

fig, ax = plt.subplots(figsize=(10, 8))
ax.plot(wk, I)
ax.set_title(u'Espectro del cuerpo negro  para T = 2.725K \
\n con datos del espectro del CMB del satélite COBE')
ax.title.set_fontsize(20)
ax.set_xlabel('Frecuencia (ciclos/cm)')
ax.xaxis.label.set_fontsize(15)
ax.set_ylabel('Radiancia espectral ($MJy \, sr^{-1}$)')
ax.yaxis.label.set_fontsize(15)
ax.set_xlim(0,25)
ax.set_ylim(0,400)

ax.scatter(df['freq'], df['I'],c='red', s= 50)
ax.grid()