How to use the galpy.actionAngle.estimateDeltaStaeckel function in galpy

### integrate the same orbit (note different calling conventions - cartesian coordinates as input)
    ### using both the orbit integration routine and the potential from Agama - much faster
    dt = time.time()
    times_c, c_orb_car = agama.orbit(ic=[ic[0],0,ic[1],ic[3],ic[5],ic[4]], \
        potential=c_pot, time=inttime, trajsize=numsteps)
    print 'Time to integrate orbit in Agama: %.4g s' % (time.time()-dt)

    ### make it compatible with galpy's convention (native output is in cartesian coordinates)
    c_orb = c_orb_car*1.0
    c_orb[:,0] = (c_orb_car[:,0]**2+c_orb_car[:,1]**2)**0.5
    c_orb[:,3] =  c_orb_car[:,2]

    ### in galpy, this is the only tool that can estimate focal distance,
    ### but it requires the orbit to be computed first
    delta = galpy.actionAngle.estimateDeltaStaeckel(g_pot, g_orb[:,0], g_orb[:,3])
    print "Focal distance estimated from the entire trajectory: Delta=%.4g" % delta

    ### plot the orbit(s) in R,z plane, along with the prolate spheroidal coordinate grid
    plt.figure(figsize=(12,8))
    plt.axes([0.04, 0.54, 0.45, 0.45])
    plotCoords(delta, 1.5)
    plt.plot(g_orb[:,0],g_orb[:,3], 'b', label='galpy')  # R,z
    plt.plot(c_orb[:,0],c_orb[:,3], 'g', label='Agama')
    plt.plot(a_orb[:,0],a_orb[:,3], 'r', label='galpy using Agama potential')
    plt.xlabel("R/8kpc")
    plt.ylabel("z/8kpc")
    plt.xlim(0, 1.2)
    plt.ylim(-1,1)
    plt.legend(loc='lower left')

    ### create galpy action/angle finder for the given value of Delta