docs: def cmag(mesh, retinotopy='any', surface=None, to='vertices'):
      '''
      cmag(mesh) yields the neighborhood-based cortical magnification for the given mesh.
      cmag(mesh, retinotopy) uses the given retinotopy argument; this must be interpretable by
        the as_retinotopy function, or should be the name of a source (such as 'empirical' or
        'any').
  
      The neighborhood-based cortical magnification data is yielded as a map whose keys are 'radial',
      'tangential', 'areal', and 'field_sign'; the units of 'radial' and 'tangential' magnifications 
      are cortical-distance/degree and the units on the 'areal' magnification is
      (cortical-distance/degree)^2; the field sign has no unit.
  
      Note that if the retinotopy source is not given, this function will by default search for any
      source using the retinotopy_data function.
  
      The option surface (default None) can be provided to indicate that while the retinotopy and
      results should be formatted for the given mesh (i.e., the result should have a value for each
      vertex in mesh), the surface coordinates used to calculate areas on the cortical surface should
      come from the given surface. The surface option may be a super-mesh of mesh.
  
      The option to='faces' or to='vertices' (the default) specifies whether the return-values should
      be for the vertices or the faces of the given mesh. Vertex data are calculated from the face
      data by summing and averaging.
      '''
      # First, find the retino data
      if pimms.is_str(retinotopy):
          retino = retinotopy_data(mesh, retinotopy)
      else:
          retino = retinotopy
      # If this is a topology, we want to change to white surface
      if isinstance(mesh, geo.Topology): mesh = mesh.white_surface
      # Convert from polar angle/eccen to longitude/latitude
      vcoords = np.asarray(as_retinotopy(retino, 'geographical'))
      # note the surface coordinates
      if surface is None: scoords = mesh.coordinates
      else:
          scoords = surface.coordinates
          if scoords.shape[1] > mesh.vertex_count:
              scoords = scoords[:, surface.index(mesh.labels)]
      faces = mesh.tess.indexed_faces
      sx = mesh.face_coordinates
      # to understand this calculation, see this stack exchange question:
      # https://math.stackexchange.com/questions/2431913/gradient-of-angle-between-scalar-fields
      # each face has a directional magnification; we need to start with the face side lengths
      (s0,s1,s2) = np.sqrt(np.sum((np.roll(sx, -1, axis=0) - sx)**2, axis=1))
      # we want a couple other handy values:
      (s0_2,s1_2,s2_2) = (s0**2, s1**2, s2**2)
      s0_inv = zinv(s0)
      b = 0.5 * (s0_2 - s1_2 + s2_2) * s0_inv
      h = 0.5 * np.sqrt(2*s0_2*(s1_2 + s2_2) - s0_2**2 - (s1_2 - s2_2)**2) * s0_inv
      h_inv = zinv(h)
      # get the visual coordinates at each face also
      vx = np.asarray([vcoords[:,f] for f in faces])
      # we already have enough data to calculate areal magnification
      s_areas = geo.triangle_area(*sx)
      v_areas = geo.triangle_area(*vx)
      arl_mag = s_areas * zinv(v_areas)
      # calculate the gradient at each triangle; this array is dimension 2 x 2 x m where m is the
      # number of triangles; the first dimension is (vx,vy) and the second dimension is (fx,fy); fx
      # and fy are the coordinates in an arbitrary coordinate system built for each face.
      # So, to reiterate, grad is ((d(vx0)/d(fx0), d(vx0)/d(fx1)) (d(vx1)/d(fx0), d(vx1)/d(fx1)))
      dvx0_dfx = (vx[2] - vx[1]) * s0_inv
      dvx1_dfx = (vx[0] - (vx[1] + b*dvx0_dfx)) * h_inv
      grad = np.asarray([dvx0_dfx, dvx1_dfx])
      # Okay, we want to know the field signs; this is just whether the cross product of the two grad
      # vectors (dvx0/dfx and dvx1/dfx) has a positive z
      fsgn = np.sign(grad[0,0]*grad[1,1] - grad[0,1]*grad[1,0])
      # We can calculate the angle too, which is just the arccos of the normalized dot-product
      grad_norms_2 = np.sum(grad**2, axis=1)
      grad_norms = np.sqrt(grad_norms_2)
      (dvx_norms_inv, dvy_norms_inv) = zinv(grad_norms)
      ngrad = grad * ((dvx_norms_inv, dvx_norms_inv), (dvy_norms_inv, dvy_norms_inv))
      dp = np.clip(np.sum(ngrad[0] * ngrad[1], axis=0), -1, 1)
      fang = fsgn * np.arccos(dp)
      # Great; now we can calculate the drad and dtan; we have dx and dy, so we just need to
      # calculate the jacobian of ((drad/dvx, drad/dvy), (dtan/dvx, dtan/dvy))
      vx_ctr = np.mean(vx, axis=0)
      (x0, y0) = vx_ctr
      den_inv = zinv(np.sqrt(x0**2 + y0**2))
      drad_dvx = np.asarray([x0,  y0]) * den_inv
      dtan_dvx = np.asarray([-y0, x0]) * den_inv
      # get dtan and drad
      drad_dfx = np.asarray([np.sum(drad_dvx[i]*grad[:,i], axis=0) for i in [0,1]])
      dtan_dfx = np.asarray([np.sum(dtan_dvx[i]*grad[:,i], axis=0) for i in [0,1]])
      # we can now turn these into the magnitudes plus the field sign
      rad_mag = zinv(np.sqrt(np.sum(drad_dfx**2, axis=0)))
      tan_mag = zinv(np.sqrt(np.sum(dtan_dfx**2, axis=0)))
      # this is the entire result if we are doing faces only
      if to == 'faces':
          return {'radial': rad_mag, 'tangential': tan_mag, 'areal': arl_mag, 'field_sign': fsgn}
      # okay, we need to do some averaging!
      mtx = simplex_summation_matrix(mesh.tess.indexed_faces)
      cols = np.asarray(mtx.sum(axis=1), dtype=np.float)[:,0]
      cols_inv = zinv(cols)
      # for areal magnification, we want to do summation over the s and v areas then divide
      s_areas = mtx.dot(s_areas)
      v_areas = mtx.dot(v_areas)
      arl_mag = s_areas * zinv(v_areas)
      # for the others, we just average
      (rad_mag, tan_mag, fsgn) = [cols_inv * mtx.dot(x) for x in (rad_mag, tan_mag, fsgn)]
      return {'radial': rad_mag, 'tangential': tan_mag, 'areal': arl_mag, 'field_sign': fsgn}