#Bradial vs bessel zeros
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

rho_0=7.112
rho_bl=rho_0+0.570
k=xroot[1]/rho_bl
dr=0.001


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

Bradial(rho)= D*dJdr(rho)*cos(theta)*cosh(k*z)