pomega2=xx^2+yy^2
pomega=sqrt(pomega2)
;
phi=atan(yy,xx)
sinp=sin(phi)
cosp=cos(phi)
;
B0=.0
bamp=1.
@parameters
rmax=1.
mu=3.
mask=make_array(size=size(phi))
mask(where(pomega le rmax))=1.
Bz=Bamp*beselj(mu*pomega,0)*mask
Bp=Bamp*beselj(mu*pomega,1)*mask
Bx=-Bp*sinp
By=+Bp*cosp
;
;  visualize B field
;
m=(m1+m2)/2
n=(n1+n2)/2
print
print,'Visualize B field constructed from Bessel functions'
contour,Bz(*,*,n),x,y,nlev=30,/fill
velovect,Bx(*,*,n),By(*,*,n),x,y,/over
;
;  B = curlA
;  curlB = curlcurlA = -del^2 A
;
bb=aa
bb(*,*,3,0)=Bx(*,*,3)
bb(*,*,3,1)=By(*,*,3)
bb(*,*,3,2)=Bz(*,*,3)
jj=curl(bb)
print
print,'Visualize J after taking the curl'
contour,jj(*,*,n,2),x,y,nlev=30,/fill
velovect,jj(*,*,n,0),jj(*,*,n,1),x,y,/over
;
aa(*,*,3,0)=-poisson_2d(reform(jj(*,*,3,0)))
aa(*,*,3,1)=-poisson_2d(reform(jj(*,*,3,1)))
aa(*,*,3,2)=-poisson_2d(reform(jj(*,*,3,2)))
;
bb=curl(aa)
bb(*,*,*,2)=bb(*,*,*,2)+B0
;jj=curl2(aa)
print
print,'Visualize B after having calculated A=-del^-2(J)'
contour,bb(*,*,n,2),x,y,nlev=30,/fill
velovect,bb(*,*,n,0),bb(*,*,n,1),x,y,/over
;stop
print
print,'Visualize radial profiles of B'
plot,x(l1:l2),bb(l1:l2,m,n,2),yr=[-1,1]*1.
oplot,x(l1:l2),bb(l1:l2,m,n,1),li=2
oplot,x(l1:l2),x(l1:l2)*0.,li=2,col=122
;oplot,[1,1]*2.4,[-1,1]
;
;  overwrite data file
;
openw, 1, datadir+'/'+varfile,/F77
writeu, 1, aa,alpm
if (lshear) then writeu, 1, t, x,y,z, dx,dy,dz, deltay else writeu, 1, t, x,y,z, dx,dy,dz
close, 1
print, 't=',t
;
END