from numpy import *
from scipy import *
from pylab import *
import numpy as np
import pylab as pl
#from matplotlib.ticker import MultipleLocator, FormatStrFormatter
import scipy as sp

Fontsize = 18
SmallFont = 15
SmallerFont = 12
TickPad = 6
MajorTickSize = 10
MajorTickWidth = 3.0
MinorTickSize = 5
MinorTickWidth = 3.0
MinorTickInterval = 2
AxisLabelWidth = 3.0
XLabelPad = 10
YLabelPad = 12.5

# set font properties with dictionary
font = {'fontname'   : 'Helvetica',
        'color'      : 'black',
        'fontweight' : 'bold',
        'fontsize'   : Fontsize}
# math mode
rcParams['mathtext.default']='regular'
# set positions of subplots
LM = 0.15  # the left side of the subplots of the figure
RM = 0.95  # the right side of the subplots of the figure
BM = 0.15  # the bottom of the subplots of the figure
TM = 0.95  # the top of the subplots of the figure
WD = 0.05  # the amount of width reserved for blank space between subplots
HT = 0.10  # the amount of height reserved for white space between subplots
pl.subplots_adjust(left=LM, bottom=BM, right=RM, top=TM, wspace=WD, hspace=HT)
# lower panel, initial size distributions
ax1 = subplot(212)
# make tick labels a bold font
setp(ax1.get_xticklabels(), color='black', fontsize=Fontsize, fontweight="bold" )
setp(ax1.get_yticklabels(), color='black', fontsize=Fontsize, fontweight="bold" )
# pad between tick marks and labels
pl.rc(("xtick.major"), pad=TickPad)
pl.rc(("ytick.major"), pad=TickPad)
# length of major ticks
pl.rc(("xtick.major"), size=MajorTickSize)
pl.rc(("ytick.major"), size=MajorTickSize)
# length of minor ticks
pl.rc(("xtick.minor"), size=MinorTickSize)
pl.rc(("ytick.minor"), size=MinorTickSize)
# pad between tic mark labels and axis labels
#pl.xlabel("...", labelpad=XLabelPad)
#pl.ylabel("...", labelpad=YLabelPad)
# width of axis lines
pl.rc("axes", linewidth=AxisLabelWidth)
# width of tick lines
pl.rc("lines", markeredgewidth=MajorTickWidth)
# tick params
ax1.tick_params(which='major',width=MajorTickWidth,length=MajorTickSize,labelsize=SmallFont)
ax1.tick_params(which='minor',width=MinorTickWidth,length=MinorTickSize)
ax1.tick_params(which='both',direction="in",top=True,right=True)
# turn minor ticks on
ax1.minorticks_on()
# label axes with defined font
ax1.set_xlabel('Radius (km)', fontsize=Fontsize, fontweight="bold")
# import data & plot 
# 
a, r, n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11, n12, n13, n14, n15, n16 = loadtxt('chc1-0111.dat', unpack=True)
rad = pow(10., r)
rp = rad / 1.e5
ratio = 1.016012 
iratio = 1. / ratio
factor = pow(rad, 2.) / (ratio - iratio)
rnum = factor * pow(10., n1)
norm = 1. / rnum[672]
np = norm * rnum
scatter(rp, np, s=20, color="black")
rnum = factor * pow(10., n2)
norm = 1. / rnum[672]
np = norm * rnum
ns = savgol_filter(np, 11, 3) # window size 11, polynomial order 3
scatter(rp, ns, s=20, color="darkmagenta")
rnum = factor * pow(10., n4)
norm = 1. / rnum[672]
np = norm * rnum
ns = savgol_filter(np, 11, 3) # window size 11, polynomial order 3
scatter(rp, ns, s=20, color="lime")
rnum = factor * pow(10., n6)
norm = 1. / rnum[672]
np = norm * rnum
ns = savgol_filter(np, 11, 3) # window size 11, polynomial order 3
scatter(rp, ns, s=20, color="darkorange")
# legend
x = [30., 60.]
y = [500., 500.]
plot(x, y, linestyle="solid", linewidth=4.0, color="darkorange")
y = [3000., 3000.]
plot(x, y, linestyle="solid", linewidth=4.0, color="lime")
y = [18000., 18000.]
plot(x, y, linestyle="solid", linewidth=4.0, color="darkmagenta")
y = [108000., 108000.]
plot(x, y, linestyle="solid", linewidth=4.0, color="black")
# labels
text(75., 64800., '0.0', fontsize=SmallFont, fontweight='bold')
text(75., 10800., '0.1', fontsize=SmallFont, fontweight='bold')
text(75., 1800., '0.5', fontsize=SmallFont, fontweight='bold')
text(75., 300., '3.0', fontsize=SmallFont, fontweight='bold')
text(5.e-4, 1.e11, r'$R$ = $N$ / $[r^{-3} (r_{up}$ - $r_{low})]$', fontsize=Fontsize, fontweight='bold', rotation=90.)
ax1.set_xscale('log')
ax1.set_yscale('log')
xticks = [1.e-4,1.e-3,1.e-2,1.e-1,1.e0,1.e1,1.e2,1.e3,1.e4,1.e5,1.e6,1.e7,1.e8]
ax1.xaxis.set_ticks(xticks)
yticks = [1.e-8,1.e-6,1.e-4,1.e-2,1.e0,1.e2,1.e4,1.e6,1.e8]
ax1.yaxis.set_ticks(yticks)
XMin = 3.e-3
XMax = 3.e2
YMin = 5.e-2
YMax = 2.e6
pl.xlim(XMin,XMax)
pl.ylim(YMin,YMax)
# end of lower plot

# upper plot
ax2 = subplot(211)
# make tick labels a bold font
setp(ax2.get_xticklabels(), visible=False)
setp(ax2.get_yticklabels(), color='black', fontsize=Fontsize, fontweight="bold" )

# pad between tick marks and labels
pl.rc(("xtick.major"), pad=TickPad)
pl.rc(("ytick.major"), pad=TickPad)
# length of major ticks
pl.rc(("xtick.major"), size=MajorTickSize)
pl.rc(("ytick.major"), size=MajorTickSize)
# length of minor ticks
pl.rc(("xtick.minor"), size=MinorTickSize)
pl.rc(("ytick.minor"), size=MinorTickSize)
# pad between tic mark labels and axis labels
#pl.xlabel("...", labelpad=XLabelPad)
#pl.ylabel("...", labelpad=YLabelPad)
# width of axis lines
pl.rc("axes", linewidth=AxisLabelWidth)
# width of tick lines
pl.rc("lines", markeredgewidth=MajorTickWidth)
# tick params
ax2.tick_params(which='major',width=MajorTickWidth,length=MajorTickSize,labelsize=SmallFont)
ax2.tick_params(which='minor',width=MinorTickWidth,length=MinorTickSize)
ax2.tick_params(which='both',direction="in",top=True,right=True)
# turn minor ticks on
ax2.minorticks_on()
# import data  &  plot 
# 
a, r, n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11, n12, n13, n14, n15, n16 = loadtxt('chc1-0111.dat', unpack=True)
rad = pow(10., r)
ratio = 1.016012 
iratio = 1. / ratio
factor = pow(rad, 2.) / (ratio - iratio)
rnum = factor * pow(10., n8)
norm = 1. / rnum[672]
np = norm * rnum
ns = savgol_filter(np, 11, 3) # window size 11, polynomial order 3
scatter(rp, np, s=15, color="darkmagenta")
rnum = factor * pow(10., n10)
norm = 1. / rnum[672]
np = norm * rnum
ns = savgol_filter(np, 11, 3) # window size 11, polynomial order 3
scatter(rp, np, s=15, color="royalblue")
rnum = factor * pow(10., n15)
norm = 1. / rnum[672]
np = norm * rnum
ns = savgol_filter(np, 11, 3) # window size 11, polynomial order 3
scatter(rp, np, s=15, color="lime")
rnum = factor * pow(10., n16)
norm = 1. / rnum[672]
np = norm * rnum
ns = savgol_filter(np, 11, 3) # window size 11, polynomial order 3
scatter(rp, np, s=15, color="darkorange")
# legend
x = [30., 60.]
y = [500., 500.]
plot(x, y, linestyle="solid", linewidth=4.0, color="darkorange")
y = [3000., 3000.]
plot(x, y, linestyle="solid", linewidth=4.0, color="lime")
y = [18000., 18000.]
plot(x, y, linestyle="solid", linewidth=4.0, color="royalblue")
y = [108000., 108000.]
plot(x, y, linestyle="solid", linewidth=4.0, color="darkmagenta")
# labels
text(75., 64800., '10', fontsize=SmallFont, fontweight='bold')
text(75., 10800., '50', fontsize=SmallFont, fontweight='bold')
text(75., 1800., '500', fontsize=SmallFont, fontweight='bold')
text(75., 300., '4500', fontsize=SmallFont, fontweight='bold')
ax2.set_xscale('log')
ax2.set_yscale('log')
#xticks = [1.e-2,1.e-1,1.e0,1.e1,1.e2,1.e3,1.e4,1.e5,1.e6,1.e7,1.e8]
yticks = [1.e-8,1.e-6,1.e-4,1.e-2,1.e0,1.e2,1.e4,1.e6,1.e8]
#ax2.xaxis.set_ticks(xticks)
ax2.yaxis.set_ticks(yticks)
XMin = 3.e-3
XMax = 3.e2
YMin = 5.e-2
YMax = 2.e6
pl.xlim(XMin,XMax)
pl.ylim(YMin,YMax)
# end of upper plot

savefig('cc1.jpg', format='jpg')
savefig('cc1.eps', format='eps')
savefig('cc1.pdf', format='pdf')
pl.show()