# Plotting the bimodal galaxy distribution from Blanton et al. 2003
# http://adsabs.harvard.edu/abs/2003ApJ...594..186B
# using a simplified RGB model based on SDSS Skyserver Image colors
# to visualize the g-r values in the histogram
#
# Markus Pössel, Haus der Astronomie, 2013
# http://www.haus-der-astronomie.de
#
# More information about this diagram: http://www.mpia.de/home/poessel/Bimodal/

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.path import Path
from matplotlib.patches import PathPatch

# fig1. dat contains a list of g-r values and frequency values (unnormalized)
# from the g-r histogram in Fig. 7 of Blanton 2003. Loading as list:
galList=loadtxt("bimodal.dat",unpack=True).tolist()

pregMr, prefreq = galList[0:2]

# Extract the bin size from the data:
binSize = (pregMr[1]-pregMr[0])

# Setting the x range within the plot by hand:
plotXrange = 1.25

# Add initial and final values to list of x-values so that pyplot.step will make a nice, closed plot:
gMr = [-binSize,pregMr[0]-binSize] + pregMr + [pregMr[-1]+binSize,plotXrange]

# Normalize the frequency list so it will show percentages; add initial and final values as per previous comment:
sumHist = sum(prefreq)
def normFreq(x): return x/sumHist*100

freq = [0,0] + map(normFreq, prefreq) + [0,0]

# Make the plot area a bit higher than the part filled by the histogram;
# I find these proportions pleasing, using phi from the golden ratio:
figHeight = max(freq)*(1+0.5*0.618)

# Since I want to show the colors in the histogram, I'm creating an image for the purpose:
imSizex = 200
imSizey = 1

# The following are simplified curves for the RGB vs. (g-r) dependence of 45 galaxies I
# sampled from the SDSS [9 (g-r) values within our range; 5 galaxies to average over
# at each value]. RGB colors are from those displayed by SDSS Skyserver.
def Rfunc (px): return np.minimum(177.2304828 +  157.51724123*px/imSizex*plotXrange,255)/255.0

def Gfunc (px): return (1/(1+np.exp(20.0*(plotXrange*px/imSizex-0.4795)))*(218.90963984*(plotXrange*px/imSizex)**1.91215098 + 178.36899012)+ 1/(1+np.exp(-20.0*(plotXrange*px/imSizex-0.4795)))*(260.00036998-59.80533941*px/imSizex*plotXrange))/255.0

def Bfunc (px):return np.minimum(177.873/(1.0+np.exp(5.476*(plotXrange*px/imSizex-0.803)))+84.13511856,255)/255.0

# Prepare an array that will become the image of the gradient:
gradientArray = zeros((imSizey,imSizex,3), dtype=float32);

# Fill this array with the proper colors determined by the functions Rfunc, Gfunc, Bfunc:
for i in range(imSizey):
    for j in range(imSizex):
        gradientArray[i][j] = np.array([Rfunc(j), Gfunc(j), Bfunc(j)]   , dtype=float32)

# The image with the gradient is not supposed to be rectangular. It is supposed to fill the histogram. To do so, it needs to be clipped with a curve that is the same as for the histogram.

# Preparing the clip path. Start with a list of points:
pathList=[]
pathLength = 2*len(gMr)

# For each data points, we will need two path points (defining the horizontal line that caps each bin of the histogram):
for i in range(len(gMr)):
    pathList.append([gMr[i]-0.5*binSize,freq[i]])
    pathList.append([gMr[i]+0.5*binSize,freq[i]])

# Make this list into a path into a patch suitable for clipping:
histoPatch = PathPatch(Path(pathList), facecolor='none', lw=0)

# Start plotting

# I like my tick labels a little larger and don't want the zero point to be too crowded:
plt.rc('xtick.major', size=6, pad=8)
plt.rc('ytick.major', size=6, pad=8)
plt.rc('xtick', color='#000000')
plt.rc('ytick', color='#000000')

# Setting figure size and defining axis ranges; clearing the figure is useful if you're in interactive mode and are trying out different things:
plt.figure(figsize=(8,6))
axis([0.0,plotXrange,0.0,figHeight])

# Here, redefine the y axis tick labels to include a percent sign:
locs,labels = plt.yticks()
loc_list=locs.tolist()

def reFormat(y): return '{0:.0f}%'.format(y)

newlabels = map(reFormat,loc_list[:-1])

plt.yticks(locs[:-1],newlabels)

#axis([0.0,plotXrange,0.0,figHeight])


# Cosmetics: Make the tick labels larger:
#ax.xticks(fontsize=14)
#ax.yticks(fontsize=14)

# Start with the bar plot of the histogram. "Mid" is important since I want the steps centered on the x-axis anchor points.
plt.step(gMr, freq, color="#555555", linewidth=1, where='mid')

# Here, I've introduced vertical lines for a pedagogical version that explains the bins. Uncomment to show those:
#for xVal in gMr: plt.vlines(xVal+0.5*binSize, 0, figHeight, color='#bbbbbb', linestyles='solid', linewidth=1)



# Add the patch to the current axis (=plot), and add the gradient image clipped by that patch:
plt.gca().add_patch(histoPatch)
plt.imshow(gradientArray,extent=[0,plotXrange,0,figHeight],aspect='auto',interpolation='sinc',clip_path=histoPatch, clip_on=True)

# Add the axis labels and title/subtitle to the diagram:
plt.ylabel("Fraction of galaxies per bin\n",horizontalalignment='right',fontsize=14)
plt.xlabel("\nGalaxy color (g-r)",fontsize=14)
plt.title("Data from Blanton et al. 2003: 183,487 galaxies from SDSS",fontsize=12)
plt.suptitle("Galaxy distribution by color\n", fontsize=16)

# Done!

# For saving, use these:

# Save as PDF:
#savefig("bimodal.pdf")

# Chosen for width of 580:
#savefig("fig1a.png", dpi=72.6)

# Higher-res PNG:
savefig("bimodal.png", dpi=300)