%matplotlib inline
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
import ephem # para las conversiones de coordenadas
from IPython.core.display import HTML # Para incluir imágenes como HTML
# 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, ephem

r = 90
theta = np.arange(0,2*np.pi,0.1)
x = np.array([0,np.pi/3,np.pi/3,0,0])
y = np.array([0,0,np.pi/4,np.pi/4,0])
fig = plt.figure(figsize=(10, 5))
ax = fig.add_subplot(111, projection="mollweide", axisbg ='LightCyan')
ax.grid(True)
ax.plot(x,y)
ax.plot(x+np.pi/6,y+np.pi/6)
ax.plot(x-np.pi/2,y-np.pi/2);

def plot_mwd(RA,Dec,org=0,title=u'Proyección de Mollweide', projection='mollweide'):
    ''' RA, Dec son arrays de la misma longitud.
    RA toma valores en [0,360), Dec en [-90,90],
    los cuales representan ángulos en grados.
    org  es el origen del gráfico, 0 o multipo de 30, en [0,360)
    title el título del gráfico
    projection es el tipo de proyección: 'mollweide', 'aitoff', 'hammer', 'lambert'
    '''
    x = np.remainder(RA+360-org,360) # desplazamos los valores de RA
    ind = x>180
    x[ind] -=360    # Conversión de escala a [-180, 180]
    x=-x    # invertir la escala: el Este a la izquierda
    tick_labels = np.array([150, 120, 90, 60, 30, 0, 330, 300, 270, 240, 210])
    tick_labels = np.remainder(tick_labels+360+org,360)
    fig = plt.figure(figsize=(10, 5))
    ax = fig.add_subplot(111, projection=projection, axisbg ='LightCyan')
    ax.scatter(np.radians(x),np.radians(Dec))  # convertir grados a radianes
    ax.set_xticklabels(tick_labels)     # Añadimos escala en el eje x
    ax.set_title(title)
    ax.title.set_fontsize(15)
    ax.set_xlabel("RA")
    ax.xaxis.label.set_fontsize(12)
    ax.set_ylabel("Dec")
    ax.yaxis.label.set_fontsize(12)
    ax.grid(True)

coord = np.array([(0,30), (60,-45), (240,15), (150,-75)])
plot_mwd(coord[:,0],coord[:,1], org=90, title = u'Test función plot_mwd')

np.random.seed(0)
nrand = np.random.rand(100,2)    # Array (100,2) de valores aleatorios en [0,1)
nrand *= np.array([360.,180.]) # array de valores en [0,360)x[0,180)
nrand -=np.array([0.,90.]) # array de valores en [0,360)x[-90,90)
nrand = nrand[(-86 < nrand[:,1]) & (nrand[:,1] < 86)] # Evitar Matplotlib Runtime Warning,
                               # manteneniendonos lejos de las singularidades en -90º y 90º
RA = nrand[:,0]
Dec = nrand[:,1]
plot_mwd(RA,Dec)

lon_array = np.arange(0,360)
lat = 0.
eq_array = np.zeros((360,2))
for lon in lon_array:
    ga = ephem.Galactic(np.radians(lon), np.radians(lat))
    eq = ephem.Equatorial(ga)
    eq_array[lon] = np.degrees(eq.get())
RA = eq_array[:,0]
Dec = eq_array[:,1]
plot_mwd(RA, Dec, 180, title = u'Plano galáctico en coordenadas ecuatoriales \n (Proyección de Mollweide)')

r = 90
theta = np.arange(0,2*np.pi,0.1)
x = np.array([0,np.pi/3,np.pi/3,0,0])
y = np.array([0,0,np.pi/4,np.pi/4,0])
fig = plt.figure(figsize=(10, 5))
ax = fig.add_subplot(111, projection="aitoff", axisbg ='LightCyan')
ax.grid(True)
ax.plot(x,y)
ax.plot(x+np.pi/6,y+np.pi/6)
ax.plot(x-np.pi/2,y-np.pi/2);

lon_array = np.arange(0,360)
lat = 0.
eq_array = np.zeros((360,2))
for lon in lon_array:
    ga = ephem.Galactic(np.radians(lon), np.radians(lat))
    eq = ephem.Equatorial(ga)
    eq_array[lon] = np.degrees(eq.get())
RA = eq_array[:,0]
Dec = eq_array[:,1]
plot_mwd(RA, Dec, 180, title = u'Plano galáctico en coordenadas ecuatoriales \n (Proyección de Aitoff)', projection='aitoff')

HTML('<img src="http://classic.sdss.org/dr7/dr7photo_big.gif" WIDTH=500>')

theta_1 = 30./180 * np.pi # 15 grados expresados en radianes

# Creamos un array con valores de theta (ángulo medido desde el polo N)
# Que delimitan regiones entre paralelos celestes de igual área
theta_a = np.zeros(7)
for n in range(0,7):
    theta_n = np.arccos(n*np.cos(theta_1)-n+1)
    theta_a[n] = theta_n
theta_a

# Los convertimos a valores de declinaciones
# medidas desde el ecuador celeste
dec_a = np.pi/2 - theta_a
dec_a

# Y finalmente representamos 6 regiones rectangulares de igual área

fig = plt.figure(figsize=(10, 5))
ax = fig.add_subplot(111, projection="mollweide", axisbg ='LightCyan')
ax.set_title(u'Ejemplo de áreas iguales \n (proyección de Mollweide)')
ax.grid(True)

for n in range(6):
    rec_n = np.array([(n* np.pi/6, dec_a[n+1]),
                      ((n +1) * np.pi/6, dec_a[n+1]),
                      ((n +1) * np.pi/6, dec_a[n]),
                      (n* np.pi/6, dec_a[n]),
                      (n* np.pi/6, dec_a[n+1])])
    x = rec_n[:,0]
    y = rec_n[:,1]
    ax.plot(x,y)