# Jean-Baptiste Juin
# Astronomy Python course
# Practice 3

#------------------------
#------------------------
# Title: Plotting with pylab/matplotlib
#------------------------
#------------------------

# If not in ipython (-pylab):
# import pylab

#------------------------
# 1D
#------------------------
# Simple curve
t=arange(0.,20.,20./100)
d=cos(t)
d2=t**2
d3=sin(t)
figure()
plot(t,d)
plot(t,d2,linestyle='-')
plot(t,d3,linestyle='.-',linewidth=5)
plot(t,d3)
plot(t,d2-0.5*d3)
ylim(-1.0,2.0)
# Reading galaxy cluster catalog from fits file
catalog=pyfits.getdata("./Practice3.Catalog.fits")
chdr=pyfits.getheader("./Practice3.Catalog.fits")
print chdr
flux=catalog[:,0]
flux0=catalog[:,1]
figure()
plot(flux,flux0,".")
loglog(flux,flux0,".")
xlabel("Flux")
ylabel("Flux0")
savefig("./Flux-Flux0.ps")
savefig("./Flux-Flux0.png")
# Scatter plotting
xv=arange(0.,10.0,10.0/20.0)
yv=rand(20)
zv=rand(20)*200
figure()
scatter(xv,yv,zv,"b")
#------------------------
# 2D image
#------------------------
# Plotting a map
Map=pyfits.getdata("./Practice3.Map.fits")
maphdr=pyfits.getheader("./Practice3.Map.fits")
print maphdr
figure()
imshow(flipud(Map))
imshow(flipud(Map),interpolation='nearest',vmin=0.0,vmax=5.0)
submap=Map[440:490,330:390]
imshow(flipud(submap),interpolation='nearest',vmin=0.0,vmax=5.0)
# Over plotting the catalog
imshow(flipud(Map),interpolation='nearest',vmin=0.0,vmax=5.0)
xpos=catalog[:,10]
ypos=catalog[:,11]
plot(xpos,ypos,'.r')
#------------------------
# Contours
#------------------------
# A cool surface
delta = 0.1
x = arange(-3.0, 3.0, delta)
y = arange(-2.0, 2.0, delta)
X, Y = meshgrid(x, y)
Z1 = bivariate_normal(X, Y, 1.0, 1.0, 0.0, 0.0)
Z2 = bivariate_normal(X, Y, 1.5, 0.5, 1, 1)
Z = 10.0 * (Z2 - Z1)
figure()
CS = contour(X, Y, Z, 10)
CS = contourf(X, Y, Z, 10, aplha=0.5)
clabel(CS, inline=1, fontsize=10)
title('Contours with labels')