#input 'number = ?'
set term pdf fontscale 0.5 size 6in,4in
outfile = '2D_histogram.pdf'
set output outfile

set timestamp offset graph 1.135,0.7 rotate font 'Verdana,6'
set label 'plotting_scripts/2D_histogram.gnunu' at graph 1.01,0.02 rotate left font 'Verdana,6' noenhance
!pwd > directory
!sed -n -e 's/^/name="/' -e 's/$/"/p' directory > name
load 'name'
print 'directory='.name
set label at graph 0.0,1.02 name font 'Verdana,6' noenhance


#x-px phase space

#spin - amplitude
ix_bins=40;iy_bins=40
xmax=0.016;xmin=0.0
ymax=0.5;ymin=-0.5
xconvert=4.e-4
yconvert=2.5e-2
x0=0.3

if(number == 7){
load '2D_hist_fit_coef_err_7.dat'
set xlabel 'displacement [m]'
set ylabel 'spin angle [rad]'
}

    xconvert = (xmax-xmin)/ix_bins
    xoffset =  (xmax+xmin)/2.
    yconvert = (ymax-ymin)/iy_bins
    yoffset = (ymax+ymin)/2.
    x0=(xmax+xmin)/2

#y-y phase space
#xmax=0.030;xmin=-0.030
#ymax=0.006;ymin=-0.006

#x_e - Amplitude phase space
#ymax=0.4e-3;ymin=0
#xmin=-0.02;xmax=0.04

#ix_bins=100;iy_bins=100
#xmax=0.045;xmin=-0.045;ix_bins=100
#ymax=0.006;ymin=-0.006;iy_bins=100
set key off
set view map
 set pm3d at s hidden3d 100
 set surface
 set dgrid3d 100,100,3
 set view map  scale 1
 set style increment default
 set style line 100 lt 5 lw 0.5
 unset hidden3d
 unset surf
 set pm3d interpolate 0,0
# set autoscale fix
# set contour base
# unset dgrid
set pm3d noborder
set palette defined ( 0 0.05 0.05 0.2, 0.1 0 0 1, 0.25 0.7 0.85 0.9,\
   0.4 0 0.75 0, 0.5 1 1 0, 0.7 1 0 0, 0.9 0.6 0.6 0.6,\
  1 0.95 0.95 0.95 )
    # set palette defined (0 0 0 0, 1 0 0 1, 3 0 1 0, 4 1 0 0, 6 1 1 1)
     set palette color
#     set palette defined ( 0 "white", 1 "black", 2 "blue" , 3"green")
s=1
fx=(xmax-xmin)/ix_bins*s
#x0=(xmax+xmin)/2*s
#x0=0
fy=(ymax-ymin)/iy_bins*s
y0=(ymax+ymin)/2*s
#y0=0
#set ylabel 'spin angle [rad]'
#set xlabel 'amplitude [cm-rad]'
#set ylabel "y' [mrad]"
#set ylabel "Amplitude ((x-D_x{/Symbol d})^2/{/Symbol b}+x'^2{/Symbol b}) [mm-mrad]"
#set xrange [xmin*s:xmax*s]

#set xtics xmin*s,2
#set yrange [ymin*s:ymax*s]
#set zrange [0.00001:]
#set logscale z
#splot '2D_hist.dat' u (($1*fx+x0)):(($2*fy+y0)):($3) w pm3d
#xconvert = yconvert*2.5/1.5
set xrange [xmin: xmax]
set yrange [ymin:ymax]
splot '2D_hist.dat' u (($1)*xconvert+x0):($2*yconvert+y0):($2>-iy_bins/2 && $2 < iy_bins/2 ? $3:1/0) w pm3d
set autoscale
#splot '2D_hist.dat' u 1:2:($2>-iy_bins/2 && $2 < iy_bins/2 ? $3:1/0) w pm3d
set xrange [-ix_bins/2:ix_bins/2]
set yrange [-iy_bins/2:iy_bins/2]
splot '2D_hist.dat' u 1:2:3 w pm3d
f(x,beta)= a+asin(b*x/sqrt(beta))
g(x,beta) =a-asin(b*x/sqrt(beta))
a=2.38
mpi= 139.57
mmu=105.658
ppi=3100.
beta0=20
b = ppi*(2*mmu)/(mpi**2-mmu**2)/100
print 'ppi=',ppi,' b=',b   
#set multiplot

set xrange [xmin: xmax]
set yrange [ymin:ymax]
splot '2D_hist.dat' u (($1)*xconvert):($2*yconvert+y0):($2>-iy_bins/2 && $2 < iy_bins/2 ? $3:1/0) w pm3d
reset
set xrange [0:1.2]
set yrange [2.1:2.7]
set size 0.825,0.7925
unset xtics
unset ytics
unset xlabel
unset ylabel
plot f(x,9) lc 1 not, g(x,9) lc 1 not, f(x,1.) lc 2 dt 2 not, g(x,1.) lc 2 dt 2 not '{/Symbol b}=1', f(x,81) lc 3 dt 2 not, g(x,81) lc 3 dt 2 not '{/Symbol b}=81'
set key below
plot f(x,9) lc 1 not, g(x,9) lc 1 t '{/Symbol b}_0=9', f(x,1) lc 2 not, g(x,1) lc 2 t '{/Symbol b}_0=1', f(x,81) lc 3 not, g(x,81) lc 3 t '{/Symbol b}_0=81'
#unset multiplot
#set autoscale
#set xrange [0:1]
#plot f(x),g(x)

print outfile