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

set timestamp 
set label 1 '$d9/g-2/examples/example10/edm/plotting_scripts/spin_vs_radial_phase_avbrad0.gnu' at graph 0.6 ,-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.55, -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}'
factory2=10000.

set label 3 at graph 0.05,0.95 'B_{radial}= measured'
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'
file = 'tbt_9.1584e-06_bl_bz_0._br_bz_2._0_-7000000._eta_0._radial_frac_2.'

plot file  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 file.'_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 225 deg azimuth'
set title 'Measurement plane at 292.5 deg azimuth'
plot file   u ($14*1.e6):($2 == 10 ? $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'
plot file   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 key inside
#plot 'location.dat' u 1:2 w l

plot file   u ($14*1.e6):($2 == 17 ? $11*factory2:1/0) w lp t '0 deg', '' u ($14*1.e6):($2 == 17 ? $12: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 == 1 ? $12:1/0) w lp t '22.5 deg', ''  u ($14*1.e6):($2 == 3 ? $11*factory2:1/0) w lp t '45 deg',\
 ''  u ($14*1.e6):($2 == 3 ? $12:1/0) w lp t '45 deg'


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) file   u ($14*1.e6):($2 == i ? $11*factory2:1/0) via a[i],b[i],c[i],d[i]
plot file   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]
g(x) = aa[i]*cos(bb[i]*x+cc[i])+dd[i]
fit g(x) file   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 file   u ($14*1.e6):($2 == i ? $12:1/0), g(x)
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 [3*pi/4:5*pi/4]
set ytics ("3pi/4" 3*pi/4,"pi" pi, "5pi/4" 5*pi/4)
plot phase_dif u ($1*22.5):2 w lp t '

set title '    '
set ylabel 'Amplitude'
factory2=10000
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