maxflux=max(max(xy(2).psi));
V1=[0:(maxflux/40):maxflux];

X1=reshape(xy(1).rz(1,:,:),mxbl(1)+1,mybl(1)+1);
Y1=reshape(xy(1).rz(2,:,:),mxbl(1)+1,mybl(1)+1);
X2=reshape(xy(2).rz(1,:,:),mxbl(2)+1,mybl(2)+1);
Y2=reshape(xy(2).rz(2,:,:),mxbl(2)+1,mybl(2)+1);
X3=reshape(xy(3).rz(1,:,:),mxbl(3)+1,mybl(3)+1);
Y3=reshape(xy(3).rz(2,:,:),mxbl(3)+1,mybl(3)+1);

[c1,h1]=contour(X1,Y1,xy(1).psi,V1);
hold on
[c2,h2]=contour(X2,Y2,xy(2).psi,V1);
[c3,h3]=contour(X3,Y3,xy(3).psi,V1);

bR=[0 .02 .25 .4 .52 .58 .58 .52 .4 .25 .02 0];
bZ=[-.06 -.06 -.29 -.29 -.03 -.03 .03 .03 .29 .29 .06 .06];
plot(bR,bZ,'k--')

axis image
[ix,iy]=find(xy(2).psi==max(max(xy(2).psi)));
r=X2(ix,iy);
z=Y2(ix,iy);
magnetic_axis=[r,z]

colorbar