def calc_potential_c(pot,R,z,rforce=False,zforce=False):
    """
    NAME:
       calc_potential_c
    PURPOSE:
       Use C to calculate the potential on a grid
    INPUT:
       pot - Potential or list of such instances
       R - grid in R
       z - grid in z
       rforce=, zforce= if either of these is True, calculate the radial or vertical force instead
    OUTPUT:
       potential on the grid (2D array)
    HISTORY:
       2013-01-24 - Written - Bovy (IAS)
       2013-01-29 - Added forces - Bovy (IAS)
    """
    from galpy.orbit.integrateFullOrbit import _parse_pot #here bc otherwise there is an infinite loop
    #Parse the potential
    npot, pot_type, pot_args= _parse_pot(pot)

    #Set up result arrays
    out= numpy.empty((len(R),len(z)))
    err= ctypes.c_int(0)

    #Set up the C code
    ndarrayFlags= ('C_CONTIGUOUS','WRITEABLE')
    if rforce:
        interppotential_calc_potentialFunc= _lib.calc_rforce
    elif zforce:
        interppotential_calc_potentialFunc= _lib.calc_zforce
    else:
        interppotential_calc_potentialFunc= _lib.calc_potential
    interppotential_calc_potentialFunc.argtypes= [ctypes.c_int,
                                                  ndpointer(dtype=numpy.float64,flags=ndarrayFlags),
                                                  ctypes.c_int,
                                                  ndpointer(dtype=numpy.float64,flags=ndarrayFlags),
                                                  ctypes.c_int,
                                                  ndpointer(dtype=numpy.int32,flags=ndarrayFlags),
                                                  ndpointer(dtype=numpy.float64,flags=ndarrayFlags),
                                                  ndpointer(dtype=numpy.float64,flags=ndarrayFlags),
                                                  ctypes.POINTER(ctypes.c_int)]

    #Array requirements, first store old order
    f_cont= [R.flags['F_CONTIGUOUS'],
             z.flags['F_CONTIGUOUS']]
    R= numpy.require(R,dtype=numpy.float64,requirements=['C','W'])
    z= numpy.require(z,dtype=numpy.float64,requirements=['C','W'])
    out= numpy.require(out,dtype=numpy.float64,requirements=['C','W'])

    #Run the C code
    interppotential_calc_potentialFunc(len(R),
                                       R,
                                       len(z),
                                       z,
                                       ctypes.c_int(npot),
                                       pot_type,
                                       pot_args,
                                       out,
                                       ctypes.byref(err))
    
    #Reset input arrays
    if f_cont[0]: R= numpy.asfortranarray(R)
    if f_cont[1]: z= numpy.asfortranarray(z)

    return (out,err.value)