set term pdf fontscale 0.5 size 6in,4in
#x-px phase space
set timestamp offset graph 1.14,0.7 rotate font ',6'
set label 'plotting scripts/2D_hist_test.gnu' at graph 1.01,0.02 rotate left font 'Verdana,6' noenhance
#the following 1-create a file named directory. 2-modify the file adding name=" at the beginning and " at the end. Then loading. Now the character string <name> is the path
!pwd > directory
!sed -n -e 's/^/name="/' -e 's/$/"/p' directory > name
load 'name'
print 'directory='.name
set label at graph 0.,1.02 name font 'Verdana,6'noenhance 
load '2D_hist_input.dat'

outfile = '2D_histogram_'.sprintf("%1i",number).'.pdf'
print outfile
set output outfile

yf=1
xf=1
if(number == 6){
ymin=ymin*1.e6
ymax=ymax*1.e6
yf=1.e6
filename = 'bin_average_rms_6.dat'
D_filename = '2D_hist_6.dat'
}
if(number == 3){
xmin=xmin*1.e6
xmax=xmax*1.e6
xf=1.e6
filename = 'bin_average_rms_3.dat'
D_filename = '2D_hist_3.dat'
}
if(number == 4){
filename = 'bin_average_rms_4.dat'
D_filename = '2D_hist_4.dat'
}
if(number == 5){
filename = 'bin_average_rms_5.dat'
D_filename = '2D_hist_5.dat'
}
print filename

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

if(number == 6){
set ylabel 'time [{/Symbol m}s]'
set xlabel '{/Symbol D}p/p [%]'
}
if(number == 3){
set xlabel 'time [{/Symbol m}s]'
set ylabel '{/Symbol D}p/p [%]'
}
if(number == 4){
set xlabel 'momentum [{/Symbol D}p/p]'
set ylabel 'x_{inf} [m]'
}
if(number == 5){
set xlabel 'momentum [{/Symbol D}p/p]'
set ylabel "x'_{inf} [rad]"
}


 xscale = (xmax-xmin)/ix_bins
 yscale=(ymax-ymin)/iy_bins
 yconvert(x) = yscale * x + (ymax+ymin)/2.
#set xrange [-0.3:0.5]
set yrange [ymin:ymax]
print D_filename
splot D_filename u ($1*xscale):(yconvert($2)):3 w pm3d
#pause -1
reset

set ylabel '<time> [{/Symbol m}s]'
set xlabel '{/Symbol D}p/p'

f(x) = offset + slope*x + q*x*x
q=0
reset
set label at graph 0.,1.02 name font 'Verdana,6'noenhance 
set timestamp offset graph 1.03,0.7 rotate font 'Verdana,6'
set label 'plotting scripts/2D_hist_test.gnu' at graph 1.015,0.02 rotate left font 'Verdana,6' noenhance

if(number == 6){
set ylabel '<time> [{/Symbol m}s]'
set xlabel '{/Symbol D}p/p'
legend  = 'd<t>/d{/Symbol d}'
units = '[{/Symbol m}s]'
}
if(number == 3){
set xlabel 'time [{/Symbol m}s]'
set ylabel '<{/Symbol D}p/p>'
legend  = 'd<t>/d{/Symbol d} = '
units = '[{/Symbol m}s]'
}
if(number == 4){
set ylabel '<x_{inf}> [m]'
set xlabel '{/Symbol D}p/p'
legend  = 'd<x_{inf}>/d{/Symbol d} = '
units='[m]'
}
if(number == 5){
set ylabel "<x'_{inf}> [rad]"
set xlabel '{/Symbol D}p/p'
legend  = "d<x'_{inf}>/d{/Symbol d} = "
units = "[rad]"
}
fit f(x) filename u ($1*xf):($2*yf):($5 > 2 ? (1./$5**.5):1/0) yerror via slope, offset
 
legends = legend.sprintf("%5.4f",slope).'{/Symbol \261}'.sprintf("%5.4f",slope_err).units
#plot filename u ($1*xf):($2*yf):($4*yf) w err not, f(x) t  legend.' ='.sprintf("%6.4f",slope).' [{/Symbol m}s]'
plot filename u ($1*xf):($2*yf):($4*yf) w err not, f(x) t  legends 


#time_to_phase = t/T*2*pi = t*omega_a = t *gamma/(gamma^2-1)* omega_c = t*gamma/(gamma^2-1)*2*pi/.1492
gamma=29.3
time_to_phase_factor = gamma/(gamma**2-1)*2*pi/.1492
ttp = time_to_phase_factor

if(number == 6){
reset
set label at graph 0.,1.02 name font 'Verdana,6'noenhance 
set ylabel '<{/Symbol f}_a> [rad]'
set xlabel '{/Symbol D}p/p'


f(x) = offset + slope*x +q*x*x
q=0
 fit f(x) filename u ($1):($2*yf*ttp):($5>2 ? (1./$5**.5):1/0) yerror via slope, offset


plot filename u ($1):($2*yf*ttp):($4*yf*ttp) w err not, f(x) t  'fit- d<t>d{/Symbol d} ='.sprintf("%6.4f",slope).' [rad]'
}

print outfile