#Bradial vs bessel zeros
#set term pdf
outfile ='bessel_vs_zero.pdf'
#set output outfile
array xroot[10]
xroot[1] = 3.8317
xroot[2] = 7.0156
xroot[3] = 10.1735
xroot[4] = 13.3237
xroot[5] = 16.4701
xroot[6]=  19.61585851
xroot[7] = 22.76008438
xroot[8] =25.90367209
xroot[9] =29.04682853
xroot[10]=32.18967991
array ratio[10]

B_l0=100.
rho_0=1007.112
rho_bl=rho_0+0.570
set xrange [0:8]

do for [i=1:10]{
k=xroot[i]/rho_bl
dr=0.0001


dJdr(rho)=(besjn(1,k*(rho+dr/2))-besjn(1,k*(rho-dr/2)))/dr

D = -B_l0*rho_0/besjn(1,k*rho_0)
z=0
theta=0
Blong(rho) = B_l0   * cosh(k*z) * besjn(1,k*rho)/besjn(1,k*rho_0)*(rho_0/rho)
Bradial(rho)= D*dJdr(rho)*cos(theta)*cosh(k*z)
set label 1 at graph 0.1,0.97'root '.sprintf("%i",i)
plot Blong(x)  t 'B_{long}',Bradial(x) t 'B_{radial}'

ratio[i]=Bradial(rho_0)/Blong(rho_0)
ratio_rho(x)=Bradial(x)/Blong(x)

set xrange [0:1000]
set samples 1000
plot Bradial(x) t 'B_{radial}', Blong(x) t 'B_{long}' 
#plot Bradial(x)/Blong(x)

pause -1

}
print ratio
unset label 1
set grid
set xlabel 'root'
set ylabel 'B_{radial}/B_{long}'
set xrange [0:12]
set title 'Radial/Longitudinal'
plot ratio
print outfile