def estimateDeltaStaeckel(pot,R,z, no_median=False):
    """
    NAME:
       estimateDeltaStaeckel
    PURPOSE:
       Estimate a good value for delta using eqn. (9) in Sanders (2012)
    INPUT:
       pot - Potential instance or list thereof
       R,z- coordinates (if these are arrays, the median estimated delta is returned, i.e., if this is an orbit)
       no_median - (False) if True, and input is array, return all calculated values of delta (useful for quickly 
       estimating delta for many phase space points)
    OUTPUT:
       delta
    HISTORY:
       2013-08-28 - Written - Bovy (IAS)
       2016-02-20 - Changed input order to allow physical conversions - Bovy (UofT)
    """
    if isinstance(R,nu.ndarray):
        delta2= nu.array([(z[ii]**2.-R[ii]**2. #eqn. (9) has a sign error
                           +(3.*R[ii]*_evaluatezforces(pot,R[ii],z[ii])
                             -3.*z[ii]*_evaluateRforces(pot,R[ii],z[ii])
                             +R[ii]*z[ii]*(evaluateR2derivs(pot,R[ii],z[ii],
                                                            use_physical=False)
                                           -evaluatez2derivs(pot,R[ii],z[ii],
                                                             use_physical=False)))/evaluateRzderivs(pot,R[ii],z[ii],use_physical=False)) for ii in range(len(R))])
        indx= (delta2 < 0.)*(delta2 > -10.**-10.)
        delta2[indx]= 0.
        if not no_median:
        	delta2= nu.median(delta2[True^nu.isnan(delta2)])
    else:
        delta2= (z**2.-R**2. #eqn. (9) has a sign error
                 +(3.*R*_evaluatezforces(pot,R,z)
                   -3.*z*_evaluateRforces(pot,R,z)
                   +R*z*(evaluateR2derivs(pot,R,z,use_physical=False)
                         -evaluatez2derivs(pot,R,z,use_physical=False)))/evaluateRzderivs(pot,R,z,use_physical=False))
        if delta2 < 0. and delta2 > -10.**-10.: delta2= 0.
    return nu.sqrt(delta2)
NAME:
       estimateDeltaStaeckel
    PURPOSE:
       Estimate a good value for delta using eqn. (9) in Sanders (2012)
    INPUT:
       pot - Potential instance or list thereof
       R,z- coordinates (if these are arrays, the median estimated delta is returned, i.e., if this is an orbit)
       no_median - (False) if True, and input is array, return all calculated values of delta (useful for quickly 
       estimating delta for many phase space points)
    OUTPUT:
       delta
    HISTORY:
       2013-08-28 - Written - Bovy (IAS)
       2016-02-20 - Changed input order to allow physical conversions - Bovy (UofT)