set term pdf fontscale 0.75 size 6in,4in enhanced
outfile='spin_vs_energy_at_time.pdf'
set output outfile 
set timestamp font 'Verdana,6' #offset graph 1.13,0.7 rotate font 'Verdana,6'
set label 4 '$ps/spin_vs_energy_at_time.gnu' at graph 1.02,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 1 at graph 0.0, 1.02 name font 'Verdana,6' noenhance


set autoscale
unset label 2
set xlabel '{/Symbol f}'
set ylabel '{/Symbol Dw}/{/Symbol w} [ppm]'
stats 'all_omega_a.dat' u (abs($3-1.44e6)<0.02e6 ? $3:1/0) 
plot 'all_omega_a.dat' u 17:(abs($3-1.44e6)<0.02e6 ? (($3-STATS_mean)/STATS_mean)*1.e6:1/0) pt 6 ps 0.2 t 'all fitted {/Symbol w}_a'

f(x) = a+b*x

set xlabel '{/Symbol D}p/p'
set ylabel '{/Symbol f} [rad]'


#plot 2
set xrange [-0.05:0.05]
plot '../EndofM5_Valetov_withInit.dat_trimmed' u 9:($15 <98 && $15>0 && abs($9)<0.05 ? $15:1/0) pt 6 ps .2 lc 1 not,'' u 9:($15<98 && $15<0 && abs($9)<0.05 ?($15+2*pi):1/0) pt 6 ps .2 lc 1 t 'Valetov Distribution'

system("rm junk.dat")
set table 'junk.dat'
plot '../EndofM5_Valetov_withInit.dat_trimmed' u 9:($15 <98 && $15>0  ? $15:1/0) pt 6 ps .2 lc 1,'' u 9:($15<98 && $15<0  ?($15+2*pi):1/0) pt 6 ps .2 lc 1  
unset table

fit f(x) 'junk.dat' u 1:2 via a,b

system("$pb/bin_data_from_two_columns junk.dat 1 2 50 -0.03 0.03")
set autoscale #unset xrange #[-0.03:0.03]
set label 2 at graph 0.6,0.9 'd{/Symbol f}/d{/Symbol d}= '.sprintf("%4.2f",b).'{/Symbol \261}'.sprintf("%4.2f",b_err).' [rad]' font ',10'

#plot 3
plot 'junk.dat_binned' u 2:3:(1/sqrt($4)) w e t  'Valetov Distribution', f(x) not

set autoscale


do for [j=1:5]{
print 'j=',j
if(j==1){n=81;m=86;p=82; t1='#4';t2='#5'}
if(j==2){n=9;m=14;p=4; t1='#6';t2='#7'}
if(j==3){n=99;m=104;p=100; t1='#8';t2='#9'}
if(j==4){n=93;m=98;p=94; t1='#10';t2='#11'}
if(j==5){n=105;m=110;p=106; t1='#12'; t2='#13'}

fit f(x) 'all_Energy_spin_time.dat' u n:m via a,b

set label 2 at graph 0.6,0.9 'd{/Symbol f}/d{/Symbol d}= '.sprintf("%4.2f",b).'{/Symbol \261}'.sprintf("%4.2f",b_err).' [rad]' font ',10'
set label 10 at graph 0.8,-0.12 'Entries='.sprintf("%6i",FIT_NDF) font',10' 


stats 'all_Energy_spin_time.dat' u p
set label 3 at graph 0.6,0.85 'time= '.sprintf("%5.2f",STATS_mean*1.e6).'{/Symbol m}s' font ',10'


#plots 4,6,8,10

plot 'all_Energy_spin_time.dat' u n:m pt 6 ps .2 not,f(x) ls -1 lw 2 not t1
system("rm junk.dat")
set table 'junk.dat'
plot 'all_Energy_spin_time.dat' u n:m pt 6 ps .2 not,f(x) ls -1 lw 2 not
unset table

system("$pb/bin_data_from_two_columns junk.dat 1 2 50 -0.004 0.004")

plot 'junk.dat_binned' u 2:3:(1/sqrt($4)) w e not, f(x) not t2 #5,7,9,11




unset label 2
unset label 3
}

#quit

#plot 14

set xlabel '{/Symbol D}p/p'
set ylabel '{/Symbol f} [rad]'

 fit f(x) 'all_MARK_INFLECTOR_DS_phase_space.dat_trimmed' u 9:($15 <10 ? $15:1/0) via a,b
set label 2 at graph 0.6,0.9 'd{/Symbol f}/d{/Symbol d}= '.sprintf("%4.2f",b).'{/Symbol \261}'.sprintf("%4.2f",b_err).' [rad]' font ',10'
set label 10 at graph 0.8,-0.12 'Entries='.sprintf("%6i",FIT_NDF) font',10' 
stats 'all_MARK_INFLECTOR_DS_phase_space.dat_trimmed' u 10
#set label 3 at graph 0.4,0.95 'time= '.sprintf("%5.2f",STATS_mean*1.e6).'{/Symbol m}s' font ',10'
plot 'all_MARK_INFLECTOR_DS_phase_space.dat_trimmed' u 9:($15 < 10. ? $15:1/0) ps .2 pt 6 not, f(x) ls -1 lw 2 not '#14'

#plot 15
system("rm junk.dat")
set table 'junk.dat'
plot 'all_MARK_INFLECTOR_DS_phase_space.dat_trimmed' u 9:($15 < 10. ? $15:1/0) ps .2 pt 6 not, f(x) ls -1 lw 2
unset table
system("$pb/bin_data_from_two_columns junk.dat 1 2 50 -0.01 0.01")
plot 'junk.dat_binned' u 2:3:(1/sqrt($4)) w e not, f(x) not '#15'

if(1 == 2){
unset label 2

#plot 14
set xlabel '{/Symbol D}p/p'
set ylabel '{/Symbol f} [rad]'


#set xrange [-0.0025:0.003]
fit f(x) 'MARK_INFLECTOR_DS_phase_space.dat_trimmed' u 9:((($9>-0.0025 && $9<0.003 && $4<0.007 && $4>-0.007 && $5 > -0.009 && $5 < 0.012 && $7 > -0.0009 && $7 < 0.0007)) ?$15:1/0) via a,b

set label 2 at graph 0.6,0.9 'd{/Symbol f}/d{/Symbol d}= '.sprintf("%4.2f",b).'{/Symbol \261}'.sprintf("%4.2f",b_err).' [rad]' font ',10'
set label 10 at graph 0.8,-0.12 'Entries='.sprintf("%6i",FIT_NDF) font',10' 

#set label 3 at graph 0.02,0.9 '-0.0025 < {/Symbol D}p/p < 0.003' font ',10'
#set label 5 at graph 0.02,0.85 '-7mm < x < 7mm' font ',10'
#set label 6 at graph 0.02,0.8 "-9 mrad < x' < 12 mrad" font ',10'
#set label 7 at graph 0.02,0.75 "-0.9 mrad < y' < 0.7 mrad" font ',10'

set xrange [-0.004:0.004]
plot 'MARK_INFLECTOR_DS_phase_space.dat_trimmed' u 9:((($9>-0.0025 && $9<0.003 && $4<0.007 && $4>-0.007 && $5 > -0.009 && $5 < 0.012 && $7 > -0.0009 && $7 < 0.0007)) ?$15:1/0) pt 6 ps .2 not, f(x) lt -1 lw 2 not  #14

#plot #15

system("rm junk.dat")
set table 'junk.dat'
plot 'MARK_INFLECTOR_DS_phase_space.dat_trimmed' u 9:((($9>-0.0025 && $9<0.003 && $4<0.007 && $4>-0.007 && $5 > -0.009 && $5 < 0.012 && $7 > -0.0009 && $7 < 0.0007)) ?$15:1/0) pt 6 ps .2 not #, f(x) lt -1 lw 2 not
unset table
system("$pb/bin_data_from_two_columns junk.dat 1 2 70 -0.004 0.004")
plot 'junk.dat_binned' u 2:3:(1/sqrt($4)) w e not, f(x) not #15

unset label 3

unset label 5
unset label 6
unset label 7

#plot 16
	h(hx,hxp) = hx**2/betax + hxp**2*betax
	v(vy,vyp) = vy**2/betay + vyp**2*betay


betax=10.2;   betay=19.9; ax=4e-05; ay=1e-05; pmin=-0.003; pmax=0.004

fit f(x) 'MARK_INFLECTOR_DS_phase_space.dat_trimmed' u 9:(h($4,$5) + v($6,$7)<ax+ay &&   pmin < $9 && $9 < pmax ? $15:1/0) via a,b
set label 2 at graph 0.6,0.9 'd{/Symbol f}/d{/Symbol d}= '.sprintf("%4.2f",b).'{/Symbol \261}'.sprintf("%4.2f",b_err).' [rad]' font ',10'
set label 10 at graph 0.8,-0.12 'Entries='.sprintf("%6i",FIT_NDF) font',10' 
#set label 8 at graph 0.02,0.89 "x'^2{/Symbol b}_x+x^2/{/Symbol b}_x + y'^2{/Symbol b}_y+y^2/{/Symbol b}_y < 50 mm-mrad" font ',10'
#set label 5 at graph 0.02,0.83 '-0.003 < {/Symbol D}p/p < 0.004' font ',10'

set xrange [-0.004:0.005]
plot 'MARK_INFLECTOR_DS_phase_space.dat_trimmed' u 9:(h($4,$5) + v($6,$7)<ax+ay &&   pmin < $9 && $9 < pmax ? $15:1/0) pt 6 ps .2 not, f(x) lt -1 lw 2 not #16

set autoscale x
unset label 5
	he(hx,hxp,delta) = (hx+eta*delta)**2/betax + (hxp+etap*delta)**2*betax +(hx+eta*delta)*(hxp+etap*delta)*alphax
	ve(vy,vyp,delta) = (vy+etay*delta)**2/betay + (vyp+etayp*delta)**2*betay + (vy+etay*delta)*(vyp+etayp*delta)*alphay
betax=10.5; betay=19.6; alphax=0.6; alphay=0.0; eta=-5.5; etap=-0.08; etay=1.7;etayp=-0.07; ax=4e-05; ay=1e-05
pmin=-0.01; pmax=0.01
fit f(x) 'MARK_INFLECTOR_DS_phase_space.dat_trimmed' u 9:(he($4,$5,$9) + ve($6,$7,$9)<ax+ay &&   pmin < $9 && $9 < pmax ? $15:1/0) via a,b

set label 2 at graph 0.6,0.9 'd{/Symbol f}/d{/Symbol d}= '.sprintf("%4.2f",b).'{/Symbol \261}'.sprintf("%4.2f",b_err).' [rad]' font ',10'
set label 10 at graph 0.8,-0.12 'Entries='.sprintf("%6i",FIT_NDF) font',10' 

setax='{/Symbol h}_x' ; setaxp= "{/Symbol h}'_x"; setay='{/Symbol h}_y'; setayp="{/Symbol h}'_y"; sbetax='{/Symbol b}_x'; sbetay='{/Symbol b}_y'; salphax='{/Symbol a}_x'; salphay='{/Symbol a}_y'; sdelta = '{/Symbol d}'; spm='{/Symbol \261}';sgammax='{/Symbol g}_x'; sgammay='{/Symbol g}_y'

#plot #17

#set label 8 at graph 0.02,0.89 "a_x=(x'+".setaxp.sdelta.')^2'.sbetax.' + (x+'.setax.sdelta.')^2'.sgammax." + (x'+".setaxp.sdelta.')(x+'.setax.sdelta.')'.salphax.' ' font ',10'
#set label 9 at graph 0.02,0.83 "a_y=(y'+".setayp.sdelta.')^2'.sbetay.' + (y+'.setay.sdelta.')^2'.sgammay." + (y'+".setayp.sdelta.')(y+'.setay.sdelta.')'.salphay.' ' font ',10'
set xrange [-0.006:0.006]
plot 'MARK_INFLECTOR_DS_phase_space.dat_trimmed' u 9:(he($4,$5,$9) + ve($6,$7,$9)<ax+ay &&   pmin < $9 && $9 < pmax ? $15:1/0) pt 6 ps .2 not, f(x) lt -1 lw 2 not #17


unset label 8
unset label 9

#plot 18
	eta = -6.3
	betax = 6
	etap = 0.08
	betay = 19.4
	ax = 4.4e-05
	ay = 1e-05
	pmin = -0.01
	pmax = 0.01
	etay = -2.3
	etayp = 0.1
	etapy = 4.0
	alphay = 0.02
	alphax = 0.6
	gammax = 0.1
	gammay = 0.07

	he(hx,hxp,delta) = (hx+eta*delta)**2*gammax + (hxp+etap*delta)**2*betax +(hx+eta*delta)*(hxp+etap*delta)*alphax
	ve(vy,vyp,delta) = (vy+etay*delta)**2*gammay + (vyp+etayp*delta)**2*betay + (vy+etay*delta)*(vyp+etayp*delta)*alphay

 fit f(x) 'MARK_INFLECTOR_DS_phase_space.dat_trimmed' u 9:(he($4,$5,$9) + ve($6,$7,$9)<ax+ay &&   pmin < $9 && $9 < pmax ? $15:1/0) via a,b
set label 2 at graph 0.6,0.9 'd{/Symbol f}/d{/Symbol d}= '.sprintf("%4.2f",b).'{/Symbol \261}'.sprintf("%4.2f",b_err).' [rad]' font ',10'
set label 10 at graph 0.8,-0.12 'Entries='.sprintf("%6i",FIT_NDF) font',10' 

 plot 'MARK_INFLECTOR_DS_phase_space.dat_trimmed' u 9:(he($4,$5,$9) + ve($6,$7,$9)<ax+ay &&   pmin < $9 && $9 < pmax ? $15:1/0) pt 6 ps .2 not, f(x) not #18

#plot 19

system("rm junk.dat")
set table 'junk.dat'
 plot 'MARK_INFLECTOR_DS_phase_space.dat_trimmed' u 9:(he($4,$5,$9) + ve($6,$7,$9)<ax+ay &&   pmin < $9 && $9 < pmax ? $15:1/0) pt 6 ps .2 not
unset table

system("$pb/bin_data_from_two_columns junk.dat 1 2 80 -0.006 0.006")

plot 'junk.dat_binned' u 2:3:(1/sqrt($4)) w e not, f(x) not  #19


set autoscale
unset label 2
set xlabel '{/Symbol f}'
set ylabel '{/Symbol Dw}/{/Symbol w} [ppm]'
stats 'all_omega_a.dat' u 3
plot 'all_omega_a.dat' u 17:(abs($3-1.44e6)<0.001e6 ? (($3-STATS_mean)/STATS_mean)*1.e6:1/0) pt 5 ps 0.45 not #20
}  #if(1 == 2)
print outfile