set term pdf fontscale 0.5 size 5in,3in
outfile = '2D_histogram_sim_p-xp.pdf'
set output outfile

#load 'fort.16'
load 'line_fit.dat'

load '2D_hist_init.dat'
dphi_dx=2.37 #mrad/mm
dgamma_dt = 5.6e-7 #/us
gamma = 29.3
omega_a=42.112 #[rad/us]

set key off
set view map
 set pm3d at s hidden3d 100
 set surface
 set dgrid3d 100,100,1
 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 10,10
# 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")


#xinj(x) = x/ix_bins*(xmax-xmin) + ix_0 !       ixx = (column(ix)-ix_0)/(xmax-xmin) * ix_bins + sign(size,column(ix)-ix_0)
xfact=(xmax-xmin)/ix_bins
#xfact=1; 
ix_0= 0 #ix_0*1000
#delta(y)= y/iy_bins*(ymax-ymin) + iy_0   !      iyy = (column(iy)-iy_0)/(ymax-ymin) * iy_bins + sign(size,column(iy)-iy_0)
yfact = (ymax-ymin)/iy_bins
#yfact=1; 
iy_0=0 #=iy_0*1000
set yrange [-0.009:0.009]
#set yrange -0.006:0.006]
set xlabel '{/Symbol d} '
set ylabel "x' [rad]"



#splot 'momentum_x_histogram.dat' u 3:4:5 w pm3d   #for data generated analytically
#plot 'momentum_x_histogram.dat' u 3:4:5 w image

splot '2D_hist.dat' u ($1*xfact+ix_0):($2*yfact+iy_0):3 w pm3d   #for data generated with bmad

#pause -1
#set autoscale
#set xrange [-50:50]
#set yrange [-50:50]


#splot '2D_hist.dat' u 1:2:3 w pm3d   #for data generated with bmad
#splot '2D_hist.dat' u ($1*xfact+ix_0):($2*yfact+iy_0):3 w image
#pause -1
reset
set xlabel '{/Symbol d} '
set ylabel "<x'> [rad]"

f(x) = offset_xp + slope_xp*x
plot 'bin_average_rms.dat' u ($1):($2):($4) w err not, f(x) t "fit-dx'/d{/Symbol d} =".sprintf("%6.4f",slope_xp).' rad'
print outfile