set term pdf fontscale 0.7 size 6in,4in
outfile = '2D_histogram.pdf'
set output outfile

load 'fort.16'
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")

set xrange [0.068:0.086]
set yrange [-0.006:0.006]
set ylabel '{/Symbol d}'
set xlabel 'x_{inf} [m]'


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:2:3 w pm3d   #for data generated with bmad
#splot '2D_hist.dat' u 1:2:3 w image

reset


load 'p_vs_x_fullfit.dat' 

set ylabel '<{/Symbol d}>'
set xlabel 'x_{inf} [m]'
load 'labels.gnu'
set yrange [-0.008:0.008]
set xrange [0.068:0.086]
set grid
print 'slope=',slope,' offset=',offset
h(x) = offset + slope * x
j(x) = 1/eta/2. * (x-0.077)

dx_ddelta = 1./slope #[m]
set key Left right


delta_omega = dphi_dx  * dgamma_dt/gamma * dx_ddelta
print ' '
print 'dphi_dx[rad/m] = ',dphi_dx,' dgamma_dt[/us] =', dgamma_dt,' dx_ddelta[m] = ', dx_ddelta
print 'delta_omega [rad/us] = ', delta_omega, ' omega_a[rad/us] = ',omega_a
print 'delta_omega/omega = ', delta_omega/omega_a*1.e9,' [ppb]'
#set label at graph 0.7,0.8 '{/Symbol Dw}/{/Symbol w}_a ='.sprintf("%6.1f",delta_omega/omega_a*1.e9).' ppb' font ',12'
plot 'p_vs_x.dat' u 2:1 w p pt 6 ps .5 t 'data' , h(x) t 'fit-d{/Symbol d}/dx ='.sprintf("%6.4f",slope), j(x) dt 2 t '1/(2{/Symbol h})='.sprintf("%6.4f",1/eta/2.)

reset

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")

set xrange [-0.006:0.006]
set yrange [0.068:0.086]

set xlabel '{/Symbol d}'
set ylabel 'x_{inf} [m]'
splot 'momentum_x_histogram.dat' u 4:3:5 w pm3d
plot 'momentum_x_histogram.dat' u 4:3:5 w image

reset
set xlabel '{/Symbol d}'
set ylabel '<x_{inf}> [m]'

load 'x_vs_momentum_fit.dat'
load 'x_vs_p_fullfit.dat' 
print 'slope=',slope,' offset=',offset
dx_ddelta = slope #m
delta_omega = dphi_dx  * dgamma_dt/gamma * dx_ddelta
print ' '
print 'dphi_dx[rad/m] = ',dphi_dx,' dgamma_dt[/us] =', dgamma_dt,' dx_ddelta[m] = ', dx_ddelta
print 'delta_omega [rad/us] = ', delta_omega, ' omega_a[rad/us] = ',omega_a
print 'delta_omega/omega = ', delta_omega/omega_a*1.e9,' [ppb]'
#set label 13 at graph 0.8,0.1 '{/Symbol Dw}/{/Symbol w}_a ='.sprintf("%6.3f",delta_omega/omega_a*1.e9).' ppb' font ',12'
#slope = 0.5646; offset = 0.075754
f(x) = slope * x + offset
j(x) = 2*eta *  x +0.077
load 'labels.gnu'
set ylabel '<x_{inf}> [m]'
set xrange [-0.006:0.006]
set yrange [0.068:0.086]
set grid

plot 'x_vs_momentum.dat' u 2:1 w p t 'data'
unset label
set label 14 at graph 0.7,0.15 'd<x_{inf}>/d{/Symbol d} ='.sprintf("%5.3f",slope) font ',12'
#set label 13 at graph 0.8,0.1 '{/Symbol Dw}/{/Symbol w}_a ='.sprintf("%6.3f",delta_omega/omega_a*1.e9).' ppb' font ',12'
plot 'x_vs_momentum.dat' u 2:1 w p t 'data', f(x) t 'fit'  #j(x) dt 2 t '2{/Symbol h}'

print outfile