diff --git a/script/analysis/plot.py b/script/analysis/plot.py
index 87bd736..63c0fb9 100644
--- a/script/analysis/plot.py
+++ b/script/analysis/plot.py
@@ -123,7 +123,8 @@ def plot_X1X3(ax, geom, var, dump, cmap='jet', vmin=None, vmax=None, cbar=True,
   if not show_axes:
       ax.set_xticklabels([]); ax.set_yticklabels([])
 
-def plot_xz(ax, geom, var, dump, cmap='jet', vmin=None, vmax=None, cbar=True, 
+
+def plot_xz(ax, geom, var, dump, cmap='jet', vmin=None, vmax=None, cbar=True,
             label=None, ticks=None, shading='gouraud',
             l_scale = None,
             reverse_x = False,
@@ -172,7 +173,6 @@ def plot_xz(ax, geom, var, dump, cmap='jet', vmin=None, vmax=None, cbar=True,
   ax.set_aspect('equal')
   ax.set_xlabel('x/M'); ax.set_ylabel('z/M')
   return mesh
-  #ax.grid(True, linestyle=':', color='k', alpha=0.5, linewidth=0.5)
 
 def contour_xz(ax, geom, var, dump,
                levels = None,
@@ -326,32 +326,42 @@ def quiver_xy(ax, geom, dump, varx, vary, C=None,
                           coordinates='figure')
 
 def overlay_field(ax, geom, dump, NLEV=20, linestyle='-', linewidth=1,
-  linecolor='k'):
+  linecolor='k', zorder=2):
   from scipy.integrate import trapz
   hdr = dump['hdr']
-  N1 = hdr['N1']; N2 = hdr['N2']
-  x = flatten_xz(geom['x'], hdr).transpose()
-  z = flatten_xz(geom['z'], hdr).transpose()
+  N1 = hdr['N1']//2; N2 = hdr['N2']
+  x = geom['x']
+  y = geom['y']
+  z = geom['z']
+  if dump['hdr']['N3'] > 1. and dump['hdr']['stopx'][3] >= np.pi:
+    x = flatten_xz(x, dump['hdr'], flip=True)
+    y = flatten_xz(y, dump['hdr'], flip=True)
+    z = flatten_xz(z, dump['hdr'])
+    rcyl = np.sqrt(x**2 + y**2)
+    rcyl[np.where(x<0)] *= -1
+  else:
+    x = x[:,:,0]
+    y = y[:,:,0]
+    z = z[:,:,0]
+    rcyl = np.sqrt(x**2 + y**2)
   A_phi = np.zeros([N2, 2*N1])
   gdet = geom['gdet'].transpose()
   B1 = dump['B1'].mean(axis=-1).transpose()
   B2 = dump['B2'].mean(axis=-1).transpose()
-  print(gdet.shape)
   for j in range(N2):
     for i in range(N1):
-      A_phi[j,N1-1-i] = (trapz(gdet[j,:i]*B2[j,:i], dx=hdr['dx'][1]) - 
+      A_phi[j,N1-1-i] = (trapz(gdet[j,:i]*B2[j,:i], dx=hdr['dx'][1]) -
                          trapz(gdet[:j, i]*B1[:j, i], dx=hdr['dx'][2]))
-      A_phi[j,i+N1] = (trapz(gdet[j,:i]*B2[j,:i], dx=hdr['dx'][1]) - 
+      A_phi[j,i+N1] = (trapz(gdet[j,:i]*B2[j,:i], dx=hdr['dx'][1]) -
                          trapz(gdet[:j, i]*B1[:j, i], dx=hdr['dx'][2]))
   A_phi -= (A_phi[N2//2-1,-1] + A_phi[N2//2,-1])/2.
   Apm = np.fabs(A_phi).max()
   if np.fabs(A_phi.min()) > A_phi.max():
     A_phi *= -1.
-  #NLEV = 20
-  levels = np.concatenate((np.linspace(-Apm,0,NLEV)[:-1], 
+  levels = np.concatenate((np.linspace(-Apm,0,NLEV)[:-1],
                            np.linspace(0,Apm,NLEV)[1:]))
-  ax.contour(x, z, A_phi, levels=levels, colors=linecolor, linestyles=linestyle,
-    linewidths=linewidth)
+  ax.contour(rcyl, z, A_phi, levels=levels, colors=linecolor, linestyles=linestyle,
+             linewidths=linewidth, zorder=zorder)
 
 def plot_xy(ax, geom, var, dump, cmap='jet', vmin=None, vmax=None, cbar=True,
   label=None, ticks=None, shading='gouraud', fix_bounds=True):
diff --git a/script/analysis/snap.py b/script/analysis/snap.py
index 725c946..a621ee2 100755
--- a/script/analysis/snap.py
+++ b/script/analysis/snap.py
@@ -27,6 +27,24 @@
 parser.add_argument('-c','--cmap',
                     type=str,default='jet',
                     help='Colormap used')
+parser.add_argument('--A-contours', action="store_true",
+                    default=False,
+                    help='Flag to overlay contours of the vector potential')
+parser.add_argument('--nlev',
+                    type=int, default=20,
+                    help='Number of positive/negative A-contour levels. Needed only if using --A-contours')
+parser.add_argument('--linestyle',
+                    type=str, default='-',
+                    help='Linestyle of A-contours. Needed only if using --A-contours')
+parser.add_argument('--linewidth', 
+                    type=int, default=1,
+                    help='Linewidth of A-contours. Needed only if using --A-contours')
+parser.add_argument('--linecolor',
+                    type=str, default='k',
+                    help='Linecolor of A-contours. Needed only if using --A-contours')
+parser.add_argument('--zorder',
+                    type=int, default=2,
+                    help='zorder for A-contours. Needed only if using --A-contours')
 parser.add_argument('--save',
                     type=str,default=None,
                     help='Figure filename if you want to save the figure')
@@ -101,6 +119,10 @@ def make_snap(dfnam,vnam,coords,size,cmap,logplot,
       bplt.plot_xz(ax, geom, var, dump, cmap=cmap, vmin=vmin, vmax=vmax,
         cbar=False, label=label, ticks=None, shading='gouraud')
       ax.set_xlim([-size,size]); ax.set_ylim([-size,size])
+    if args.A_contours:
+      bplt.overlay_field(ax, geom, dump, NLEV=args.nlev,
+                         linestyle=args.linestyle, linewidth=args.linewidth,
+                         linecolor=args.linecolor, zorder=args.zorder)
     ax = a1
     if coords == 'mks':
       bplt.plot_X1X3(ax, geom, var, dump, cmap=cmap, vmin=vmin, vmax=vmax, 
@@ -120,6 +142,10 @@ def make_snap(dfnam,vnam,coords,size,cmap,logplot,
       bplt.plot_xz(ax, geom, var, dump, cmap=cmap, vmin=vmin, vmax=vmax,
         cbar=True, label=label, ticks=None, shading='gouraud')
       ax.set_xlim([0,size]); ax.set_ylim([-size,size])
+    if args.A_contours:
+      bplt.overlay_field(ax, geom, dump, NLEV=args.nlev,
+                         linestyle=args.linestyle, linewidth=args.linewidth,
+                         linecolor=args.linecolor, zorder=args.zorder)
 
   if savefig == False:
     plt.show()