Plot an Analema visible from Madrid

AnalemaMadrid

#!/usr/bin/python

"""
    Traza la Analema para un Madrid
    
    @precondition: se necesita tener el paquete Ephem
    
"""

import ephem
from pylab import *

obs = ephem.Observer()
obs.long = '-3:44:16'
obs.lat = '40:29:4'
hour = 12 - ephem.degrees(obs.long)/pi*12
obs.date = (2012, 1, 1, hour)

sun = ephem.Sun()
sun_alt = []
sun_az = []

for i in range(365):   
    sun.compute(obs)   
    sun_az.append(sun.az*180/pi)   
    sun_alt.append(sun.alt*180/pi)   
    obs.date = obs.date + 1

plot(sun_az, sun_alt, '.')
axis('equal')
grid()
show()