#file1='tbt_0._bl_bz_100.e-6_br_bz_0.0_1_-7000000._eta_0.0_radial_frac_0.0'
#file1 = 'Bz100ppm_cos_noradial.'; label1='B_{long}=100 ppm, N=1, cos, B_{radial}=0'

#file1='Bz100ppm_sin_noradial.'; label1='B_{long}=100 ppm, N=1, sin, B_{radial}=0'

file1='Bz100ppm_sin_zero7.'; label1='B_{long}=100 ppm, N=1, sin, Maxwell =>  B_{radial}; Order 7'

#file1='Bz0_sin_Brzero7.'; label1='B_{long}=0 ppm, N=1, sin, Maxwell =>  B_{radial}; Order 7'

#file1 ='Bz0Brad100ppb_uniform.'; label1='B_{long}=0, B_{radial}=100ppb; uniform'

#file1= 'Bz0Brad100ppm_cos.'; label1 = 'B_{long}=0, B_{radial}=100ppm, cos'

set term pdf fontscale 0.5 size 6in,4in
outfile = file1.'pdf'
set output outfile

set timestamp 
set label 1 '$d9/g-2/examples/example10/edm/plotting_scripts/spin_vs_long_phase.gnu' at graph 0.5 ,-0.10 font 'Verdana,6' noenhance
!pwd > directory
!sed -n -e 's/^/name="/' -e 's/$/"/p' directory > name
load 'name'
print 'directory='.name
#set label 2 at graph 0.0,1.02 name font 'Verdana,10' noenhance
set label 2 at graph 0.45, -0.12 name font 'Verdana,8' noenhance
set title ' '
set xrange [0:10]
set xlabel 'time [us]'
set ylabel 's_{long}'
set y2label 's_{vert}'
set yrange [-1:1]
set ytics -1,.2
set  y2range [-1.:1.]
set y2tics -1.,0.2
set label 3 at graph 0.05,0.95 '{/Symbol h}=10^{-6}'
set label 5 at graph 0.05,0.9 'd = 4.67 X 10^{-20} e-cm'
factory2=1000.
print log10(factory2)
#set label 1 at graph 1, 0 sprintf("× 10^{%d}",log10(factorX)) offset 0.5,0.3
set label 4 at graph 1.04, .99 sprintf("× 10^{%i}",log10(1./10000.)) offset 0,0.7
plot 'tbt_0.0_bl_bz_0._br_bz_0._0_-7000000._eta_1.e-6' u ($14*1.e6):($11*factory2) w l axes x1y2 t 's_{vert}', '' u ($14*1.e6):12 w l dt 2 t 's_{long}'


stats file1 u 11:12
int = STATS_max_y/STATS_max_x
int = int(log10(int))
print 'STATS_max_x =', STATS_max_x, ' STATS_max_y = ', STATS_max_y, ' int = ', int

factory2=10.**int


set label 3 at graph 0.05,0.95 label1
set label 4 at graph 1.04, .99 sprintf("× 10^{%i}",log10(1./factory2*.99999999)) offset 0,0.7
set label 5 at graph 0.05,0.9 'On closed orbit'

plot file1  u ($14*1.e6):($11*factory2) w l axes x1y2 t 's_{vert}', '' u ($14*1.e6):12 w l dt 2 t 's_{long}'
plot file1.'_Filt' u 2:($6*factory2) w l axes x1y2 t 's_{vert}', '' u 2:5 w l dt 2 t 's_{long}'

#set table 'location.dat'
#set label 6 at graph 0.05,0.04 'Measurement plane at 292.5 deg azimuth'
set title 'Measurement plane at 292.5 deg azimuth'
plot file1   u ($14*1.e6):($2 == 13 ? $11*factory2:1/0) w lp axes x1y2  t 's_{vert}','' u ($14*1.e6):($2 == 13 ? $12:1/0) w lp t 's_{long}'
unset label 6

set key outside
unset y2label
unset y2tics
set ylabel 's_{vert}'
set label 4 at graph -0.04, 1.01 sprintf("× 10^{%i}",log10(1./factory2*.99999999)) offset 0,0.7


set title 'Measurement planes at 0,22.5,45, ..., 292.5 deg azimuth'

set yrange [-0.2:0.2]
plot file1   u ($14*1.e6):($2 == 17 ? $11*factory2:1/0) w lp t '0 deg',\
 ''  u ($14*1.e6):($2 == 1 ? $11*factory2:1/0) w lp t '22.5 deg', ''  u ($14*1.e6):($2 == 2 ? $11*factory2:1/0) w lp t '45 deg',\
 ''  u ($14*1.e6):($2 == 4 ? $11*factory2:1/0) w lp t '90 deg' ,\
 ''  u ($14*1.e6):($2 == 6 ? $11*factory2:1/0) w lp t '135 deg', ''  u ($14*1.e6):($2 == 8 ? $11*factory2:1/0) w lp t '180 deg', ''  u ($14*1.e6):($2 == 10 ? $11*factory2:1/0) w lp t '225 deg',\
 ''  u ($14*1.e6):($2 == 12 ? $11*factory2:1/0) w lp t '270 deg', ''  u ($14*1.e6):($2 == 13 ? $11*factory2:1/0) w lp t '292.5 deg'
unset title


#unset table
set yrange [-1:1]

set key inside
#plot 'location.dat' u 1:2 w l

plot file1   u ($14*1.e6):($2 == 17 ? $11*factory2:1/0) w lp t '0 deg s_y', '' u ($14*1.e6):($2 == 17 ? $12:1/0) w lp t '0 deg s_z',\
 ''  u ($14*1.e6):($2 == 1 ? $11*factory2:1/0) w lp t '22.5 deg s_y', ''  u ($14*1.e6):($2 == 1 ? $12:1/0) w lp t '22.5 deg s_z', ''  u ($14*1.e6):($2 == 3 ? $11*factory2:1/0) w lp t '45 deg z_y',\
 ''  u ($14*1.e6):($2 == 3 ? $12:1/0) w lp t '45 deg s_z'

plot file1   u ($14*1.e6):($2 == 16 ? $11*factory2:1/0) w lp t '0 deg s_y', '' u ($14*1.e6):($2 == 16 ? $12:1/0) w lp t '0 deg s_z'
plot file1   u ($14*1.e6):($2 == 1 ? $11*factory2:1/0) w lp t '22.5 deg s_y', '' u ($14*1.e6):($2 == 1 ? $12:1/0) w lp t '22.5 deg s_z'
plot file1   u ($14*1.e6):($2 == 2 ? $11*factory2:1/0) w lp t '45 deg s_y', '' u ($14*1.e6):($2 == 2 ? $12:1/0) w lp t '45 deg s_z'
plot file1   u ($14*1.e6):($2 == 3 ? $11*factory2:1/0) w lp t '67.5 deg s_y', '' u ($14*1.e6):($2 == 3 ? $12:1/0) w lp t '67.5 deg s_z'
plot file1   u ($14*1.e6):($2 == 4 ? $11*factory2:1/0) w lp t '90 deg s_y', '' u ($14*1.e6):($2 == 4 ? $12:1/0) w lp t '90 deg s_z'
plot file1   u ($14*1.e6):($2 == 5 ? $11*factory2:1/0) w lp t '112.5 deg s_y', '' u ($14*1.e6):($2 == 5 ? $12:1/0) w lp t '112.5 deg s_z'
plot file1   u ($14*1.e6):($2 == 6 ? $11*factory2:1/0) w lp t '135 deg s_y', '' u ($14*1.e6):($2 == 6 ? $12:1/0) w lp t '135 deg s_z'
plot file1   u ($14*1.e6):($2 == 7 ? $11*factory2:1/0) w lp t '157.5 deg s_y', '' u ($14*1.e6):($2 == 7 ? $12:1/0) w lp t '157.5 deg s_z'
plot file1   u ($14*1.e6):($2 == 8 ? $11*factory2:1/0) w lp t '180 deg s_y', '' u ($14*1.e6):($2 == 8 ? $12:1/0) w lp t '180 deg s_z'

array a[17]
array b[17]
array c[17]
array d[17]
array aa[17]
array bb[17]
array cc[17]
array dd[17]
array cosamp[17]
array sinamp[17]
array phase_dif[17]

a[1]=-7.32353740135554e-05*factory2; b[1]= 1.43934189240645; c[1]= -1.88850748100537; d[1]= 9.37674641284535e-06*factory2
aa[1]=1.; bb[1]= 1.43934189240645; cc[1]= -1.88850748100537; dd[1]= 0.1
do for [i=1:16]{
f(x) = a[i]*cos(b[i]*x+c[i])+d[i]
fit f(x) file1   u ($14*1.e6):($2 == i ? $11*factory2:1/0) via a[i],b[i],c[i],d[i]
#plot file1   u ($14*1.e6):($2 == i ? $11*factory2:1/0), f(x)

a[i+1]=a[i];b[i+1]=b[i];c[i+1]=c[i];d[i+1]=d[i]
a[i+1]=a[1];b[i+1]=b[1];c[i+1]=c[1];d[i+1]=d[1]

g(x) = aa[i]*cos(bb[i]*x+cc[i])+dd[i]
fit g(x) file1   u ($14*1.e6):($2 == i ? $12:1/0) via aa[i],bb[i],cc[i],dd[i]
phase_dif[i] = c[i]-cc[i]
cosamp[i] = a[i]*cos(phase_dif[i])
sinamp[i] = -a[i]*sin(phase_dif[i])
#plot file1   u ($14*1.e6):($2 == i ? $12:1/0), g(x)
#pause -1
aa[i+1]=aa[i];bb[i+1]=bb[i];cc[i+1]=cc[i];dd[i+1]=dd[i]
}
print a
print b
print c
print d
reset
set xrange [0:360]
set xtics 0,60
set xlabel 'azimuth [deg]'
set ylabel 'amplitude'
plot a u ($1*22.5):2 w lp t 'amplitude'
set autoscale y
set ylabel 'frequency'
plot b u ($1*22.5):2  t 'frequency'

set autoscale y
set ylabel 'phase'
plot c u ($1*22.5):2 t 's_{vert}',cc u ($1*22.5):2 t  's_{long}'

set  grid
set ylabel '{/Symbol f}_{vert}-{/Symbol f}_{long}'
set yrange [-pi:pi]
set ytics ("-pi" -pi, "-pi/2" -pi/2, 0, "pi/2" pi/2, "pi" pi)
plot phase_dif u ($1*22.5):2 w lp t '

set title '    '
set ylabel 'Amplitude'

factory2=1000
set label 4 at graph -0.04, 1.01 sprintf("× 10^{%i}",log10(1./factory2*.99999999)) offset 0,0.7
set yrange [-1:1]
set ytics -1,.2
plot cosamp u ($1*22.5):2 w lp t 'A_{g-2}', sinamp u ($1*22.5):2 w lp t 'A_{EDM}'

set autoscale y
plot d u ($1*22.5):2 t 'offset'

print outfile