#set term pdf fontscale 0.75 size 6in,4in
outfile = 'ddelta_dx.pdf'
set output outfile
index = 0.108
R0=7.112
eta=R0/(1-index)
beta=R0/sqrt(1-index)
Ap=0.045 
x0=0.077 
theta_0=x0/beta
b=-0.325
bb=1.642
bb=0.56

theta_low = (x0-Ap/2)/beta
theta_low_mid = (x0-Ap/6)/beta

theta_high= (x0+Ap/2)/beta
theta_high_mid= (x0+Ap/6)/beta
A1=(beta*bb)**2+1
theta=theta_low
B1(x)=-2*((bb*beta)**2*x0+x*beta)
C1(x)=-Ap**2+(x*beta)**2+(beta*bb*x0)**2
xmn(x)=(-B1(x)-sqrt(B1(x)**2-4*A1*C1(x)))/2./A1
xmx(x)=(-B1(x)+sqrt(B1(x)**2-4*A1*C1(x)))/2./A1
minh=xmn(theta_low)
maxh=xmx(theta_low)
print theta_low,minh,maxh

min0=xmn(theta_0)
max0=xmx(theta_0)
print theta_0,min0,max0

minl=xmn(theta_high)
maxl=xmx(theta_high)
print theta_high,minl,maxl

reset

set xrange [0.068:0.086]
set samples 10000
set grid
plot (bb*beta*(x-x0))**2
pause -1

B(x,theta)=1./(Ap**2-(x-theta*beta)**2)
D(x,theta) = 1.-(bb*beta*(x-x0))**2*B(x,theta)
plot B(x,theta)
pause -1
plot D(x,theta_low), D(x,theta_0), D(x,theta_high)
pause -1
#f(x,amp,theta) = 1./eta*(1- B(x,theta) *  ( beta_xp2(x,amp)*(1+2*(x-theta*beta)**2*B(x,theta) )-2*(x-x0)*(x-theta*beta) ))*100

flow(x,theta) = xmn(theta_high)<x && x<xmx(theta_high) ? 0.5/eta*(1-beta**2*bb**2*(x-x0)*B(x,theta)*((x-x0)*(1+2*(x-theta*beta)**2*B(x,theta))+2*(x-beta*theta))) :1/0
flow1(x,theta) = xmn(theta_high_mid)<x && x<xmx(theta_high_mid) ? 0.5/eta*(1-beta**2*bb**2*(x-x0)*B(x,theta)*((x-x0)*(1+2*(x-theta*beta)**2*B(x,theta))+2*(x-beta*theta))) :1/0
f(x,theta) = min0<x && x<max0 ? 0.5/eta*(1-beta**2*bb**2*(x-x0)*B(x,theta)*((x-x0)*(1+2*(x-theta*beta)**2*B(x,theta))+2*(x-beta*theta))) :1/0
fhigh1(x,theta) = xmn(theta_low_mid)<x && x<xmx(theta_low_mid) ? 0.5/eta*(1-beta**2*bb**2*(x-x0)*B(x,theta)*((x-x0)*(1+2*(x-theta*beta)**2*B(x,theta))+2*(x-beta*theta))) :1/0
fhigh(x,theta) = minh<x && x<maxh ? 0.5/eta*(1-beta**2*bb**2*(x-x0)*B(x,theta)*((x-x0)*(1+2*(x-theta*beta)**2*B(x,theta))+2*(x-beta*theta))) :1/0

finv(x,theta) = 1./(0.5/eta*(1-beta**2*bb**2*(x-x0)*B(x,theta)*((x-x0)*(1+2*(x-theta*beta)**2*B(x,theta))-2*(x-beta*theta))))
g(x,theta) = 1/f(x,theta)*0.45e-7*1.e6/42.112e6*1.e9
hlow(x,theta)=xmn(theta_low)<x && x<xmx(theta_low) ? 0.5/eta*(x-theta*beta)*(1-(beta*bb*(x-x0))**2*B(x,theta)) : 1/0
hlow1(x,theta)=xmn(theta_low_mid)<x && x<xmx(theta_low_mid) ? 0.5/eta*(x-theta*beta)*(1-(beta*bb*(x-x0))**2*B(x,theta)) : 1/0
h(x,theta)=xmn(theta_0)<x && x<xmx(theta_0) ?   0.5/eta*(x-theta*beta)*(1-(beta*bb*(x-x0))**2*B(x,theta)): 1/0
hhigh1(x,theta)= xmn(theta_high_mid)<x && x<xmx(theta_high_mid) ? 0.5/eta*(x-theta*beta)*(1-(beta*bb*(x-x0))**2*B(x,theta)) : 1/0
hhigh(x,theta)= xmn(theta_high)<x && x<xmx(theta_high) ? 0.5/eta*(x-theta*beta)*(1-(beta*bb*(x-x0))**2*B(x,theta)) : 1/0
set xrange [0.068:0.086]
set grid
set yrange [-1.6:1]
set ytics -1.6,0.2
set ylabel 'd{/Symbol d}/dx_{inf} [m^{-1}]'
set xlabel 'displacement at inflector exit [m]'
#set key center Left right
set term x11 0
plot flow(x,theta_low) t '{/Symbol q}_k ='.sprintf("%3.1f",theta_low*1000).' mrad',\
flow1(x,theta_high_mid) t '{/Symbol q}_k ='.sprintf("%3.1f",theta_high_mid*1000).' mrad',\
f(x,theta_0) lt -1 t '{/Symbol q}_k ='.sprintf("%3.1f",theta_0*1000).' mrad',\
fhigh1(x,theta_low_mid) t '{/Symbol q}_k ='.sprintf("%3.1f",theta_low_mid*1000).' mrad',\
fhigh(x,theta_high) t '{/Symbol q}_k ='.sprintf("%3.1f",theta_high*1000).' mrad'
print outfile
pause -1
outfile = 'delta_omega.pdf'
set output outfile

set ylabel '{/Symbol Dw}/{/Symbol w} [ppb]'
set yrange [-4:64]
set ytics -8,8
set term x11 1
plot g(x,theta_low) t '{/Symbol q}_k ='.sprintf("%3.1f",theta_low*1000).' mrad',g(x,theta_0) lt -1 t '{/Symbol q}_k ='.sprintf("%3.1f",theta_0*1000).' mrad',g(x,theta_high) t '{/Symbol q}_k ='.sprintf("%3.1f",theta_high*1000).' mrad'
print outfile
pause -1
set autoscale y
#set key at 0.076,0.002
outfile = 'delta_vs_x.pdf'
set output outfile
set samples 10000
set yrange [-0.002:0.002]
set ytics -0.002,0.0005
set ylabel '{/Symbol d}'
set term x11 2
plot hlow(x,theta_low) t '{/Symbol q}_k ='.sprintf("%3.1f",theta_low*1000).' mrad',\
hlow1(x,theta_low_mid) t '{/Symbol q}_k ='.sprintf("%3.1f",theta_low_mid*1000).' mrad',\
h(x,theta_0) lt -1 t '{/Symbol q}_k ='.sprintf("%3.1f",theta_0*1000).' mrad', \
hhigh1(x,theta_high_mid) t '{/Symbol q}_k ='.sprintf("%3.1f",theta_high_mid*1000).' mrad',\
hhigh(x,theta_high) t '{/Symbol q}_k ='.sprintf("%3.1f",theta_high*1000).' mrad'
print outfile

set term x11 4
set xrange [(x0-Ap)/beta:(x0+Ap)/beta]
set xlabel 'kick [mrad]'
set autoscale y
set ylabel '<d{/Symbol d}/dx_{inf}> m^{-1}'
set yrange [-0.02:0.02]
set ytics -0.02, 0.002
set multiplot

xml=0.068
xmh=0.086
theta=theta_low
if(xmn(theta) > 0.068){xml=xmn(theta)+1.e-7}
if(xmx(theta) < 0.086){xmh=xmx(theta)-1.e-7}
print 'xml=',xml,'  xmh=',xmh
y2= hlow(xml,theta)
#print 'B(xmh,theta) =',B(xmh,theta), '  (bb*beta*(xmh-x0))**2= ',(bb*beta*(xmh-x0))**2, '  B*num =', B(xmh,theta_low)*(bb*beta*(xmh-x0))**2
y1= hlow(xmh,theta)
y3=(y1-y2)/(xmh-xml)
print ' '
print 'kick theta[mrad] = ',theta, ' Range= ', xml, xmh,  '  Average <d delta/dx> [1/m] =', y3
set label at theta,y3 point pt 7 ps 1

xml=0.068
xmh=0.086
theta=theta_low_mid
if(xmn(theta) > 0.068){xml=xmn(theta)+1.e-7}
if(xmx(theta) < 0.086){xmh=xmx(theta)-1.e-7}
#print 'xml=',xml,'  xmh=',xmh
y2= hlow1(xml,theta)
#print 'B(xmh,theta) =',B(xmh,theta), '  (bb*beta*(xmh-x0))**2= ',(bb*beta*(xmh-x0))**2, '  B*num =', B(xmh,theta_low)*(bb*beta*(xmh-x0))**2
y1= hlow1(xmh,theta)
y3=(y1-y2)/(xmh-xml)
print ' '
print 'kick theta[mrad] = ',theta, ' Range= ', xml, xmh,  '  Average <d delta/dx> [1/m] =', y3
set label at theta,y3 point pt 7 ps 1

xml=0.068
xmh=0.086
theta=theta_0
if(xmn(theta) > 0.068){xml=xmn(theta)+1.e-7}
if(xmx(theta) < 0.086){xmh=xmx(theta)-1.e-7}
#print 'xml=',xml,'  xmh=',xmh
y2= h(xml,theta)
#print 'B(xmh,theta) =',B(xmh,theta), '  (bb*beta*(xmh-x0))**2= ',(bb*beta*(xmh-x0))**2, '  B*num =', B(xmh,theta_low)*(bb*beta*(xmh-x0))**2
y1= h(xmh,theta)
y3=(y1-y2)/(xmh-xml)
print ' '
print 'kick theta[mrad] = ',theta, ' Range= ', xml, xmh,  '  Average <d delta/dx> [1/m] =', y3
set label at theta,y3 point pt 7 ps 1

xml=0.068
xmh=0.086
theta=theta_high_mid
if(xmn(theta) > 0.068){xml=xmn(theta)+1.e-7}
if(xmx(theta) < 0.086){xmh=xmx(theta)-1.e-7}
#print 'xml=',xml,'  xmh=',xmh
y2= hhigh1(xml,theta)
#print 'B(xmh,theta) =',B(xmh,theta), '  (bb*beta*(xmh-x0))**2= ',(bb*beta*(xmh-x0))**2, '  B*num =', B(xmh,theta_low)*(bb*beta*(xmh-x0))**2
y1= hhigh1(xmh,theta)
y3=(y1-y2)/(xmh-xml)
print ' '
print 'kick theta[mrad] = ',theta, ' Range= ', xml, xmh,  '  Average <d delta/dx> [1/m] =', y3
set label at theta,y3 point pt 7 ps 1

xml=0.068
xmh=0.086
theta=theta_high
if(xmn(theta) > 0.068){xml=xmn(theta)+1.e-7}
if(xmx(theta) < 0.086){xmh=xmx(theta)-1.e-7}
#print 'xml=',xml,'  xmh=',xmh
y2= hhigh(xml,theta)
#print 'B(xmh,theta) =',B(xmh,theta), '  (bb*beta*(xmh-x0))**2= ',(bb*beta*(xmh-x0))**2, '  B*num =', B(xmh,theta_low)*(bb*beta*(xmh-x0))**2
y1= hhigh(xmh,theta)
y3=(y1-y2)/(xmh-xml)
print ' '
print 'kick theta[mrad] = ',theta, ' Range= ', xml, xmh,  '  Average <d delta/dx> [1/m] =', y3
set label at theta,y3 point pt 7 ps 1

plot '-' w p pt 5 ps 0.1 not
0.0102 0.
e
unset multiplot
unset label

set term x11 3
set xrange [0.0072:0.0132]
plot (xmx(x)-xmn(x))/2. w l