def METHOD_NAME(plot=False):
  #print "roche contact"
  q = 0.5
  F = 0.5
  d = 1
  Omega0 = 2.65
  choice = 2
  # determine area and volume 
  r_av = libphoebe.roche_area_volume(q, F, d, Omega0, choice, larea=True, lvolume=True)
  #print r_av
  # calculate delta
  ntriangles=2000
  delta = np.sqrt(r_av["larea"]/(np.sqrt(3)*ntriangles/4))
  max_triangles = int(1.5*ntriangles)
  # generate mesh  
  r_mesh = libphoebe.roche_marching_mesh(q, F, d, Omega0, delta, choice, max_triangles, 
                            vertices=True, vnormals=True, triangles=True,
                            tnormals=True, areas=True, area=True, volume=True, full=True)  
  # check the number of triagles of the mesh
  #print len(r_mesh["triangles"])
  assert(0.75*ntriangles < len(r_mesh["triangles"]) < 1.35*ntriangles)
  # check is the vertices are really on isosurface
  #print np.max(np.abs([potential_roche(r, q, F, d) - Omega0 for r in r_mesh["vertices"]]))
  assert(np.max(np.abs([potential_roche(r, q, F, d) - Omega0 for r in r_mesh["vertices"]])) < 1e-12)
  # check the area
  #print (r_mesh["area"]-r_av["larea"])/r_av["larea"]
  assert (np.abs(r_mesh["area"]-r_av["larea"])/r_av["larea"] < 1e-2)
  # check the volume
  #print (r_mesh["volume"]- r_av["lvolume"])/r_av["lvolume"]
  assert (np.abs(r_mesh["volume"]- r_av["lvolume"])/r_av["lvolume"] < 1e-2)