def cmag(mesh, retinotopy='<STR_LIT>', surface=None, to='<STR_LIT>'):
<EOL>if pimms.is_str(retinotopy):<EOL><INDENT>retino = retinotopy_data(mesh, retinotopy)<EOL><DEDENT>else:<EOL><INDENT>retino = retinotopy<EOL><DEDENT>if isinstance(mesh, geo.Topology): mesh = mesh.white_surface<EOL>vcoords = np.asarray(as_retinotopy(retino, '<STR_LIT>'))<EOL>if surface is None: scoords = mesh.coordinates<EOL>else:<EOL><INDENT>scoords = surface.coordinates<EOL>if scoords.shape[<NUM_LIT:1>] > mesh.vertex_count:<EOL><INDENT>scoords = scoords[:, surface.index(mesh.labels)]<EOL><DEDENT><DEDENT>faces = mesh.tess.indexed_faces<EOL>sx = mesh.face_coordinates<EOL>(s0,s1,s2) = np.sqrt(np.sum((np.roll(sx, -<NUM_LIT:1>, axis=<NUM_LIT:0>) - sx)**<NUM_LIT:2>, axis=<NUM_LIT:1>))<EOL>(s0_2,s1_2,s2_2) = (s0**<NUM_LIT:2>, s1**<NUM_LIT:2>, s2**<NUM_LIT:2>)<EOL>s0_inv = zinv(s0)<EOL>b = <NUM_LIT:0.5> * (s0_2 - s1_2 + s2_2) * s0_inv<EOL>h = <NUM_LIT:0.5> * np.sqrt(<NUM_LIT:2>*s0_2*(s1_2 + s2_2) - s0_2**<NUM_LIT:2> - (s1_2 - s2_2)**<NUM_LIT:2>) * s0_inv<EOL>h_inv = zinv(h)<EOL>vx = np.asarray([vcoords[:,f] for f in faces])<EOL>s_areas = geo.triangle_area(*sx)<EOL>v_areas = geo.triangle_area(*vx)<EOL>arl_mag = s_areas * zinv(v_areas)<EOL>dvx0_dfx = (vx[<NUM_LIT:2>] - vx[<NUM_LIT:1>]) * s0_inv<EOL>dvx1_dfx = (vx[<NUM_LIT:0>] - (vx[<NUM_LIT:1>] + b*dvx0_dfx)) * h_inv<EOL>grad = np.asarray([dvx0_dfx, dvx1_dfx])<EOL>fsgn = np.sign(grad[<NUM_LIT:0>,<NUM_LIT:0>]*grad[<NUM_LIT:1>,<NUM_LIT:1>] - grad[<NUM_LIT:0>,<NUM_LIT:1>]*grad[<NUM_LIT:1>,<NUM_LIT:0>])<EOL>grad_norms_2 = np.sum(grad**<NUM_LIT:2>, axis=<NUM_LIT:1>)<EOL>grad_norms = np.sqrt(grad_norms_2)<EOL>(dvx_norms_inv, dvy_norms_inv) = zinv(grad_norms)<EOL>ngrad = grad * ((dvx_norms_inv, dvx_norms_inv), (dvy_norms_inv, dvy_norms_inv))<EOL>dp = np.clip(np.sum(ngrad[<NUM_LIT:0>] * ngrad[<NUM_LIT:1>], axis=<NUM_LIT:0>), -<NUM_LIT:1>, <NUM_LIT:1>)<EOL>fang = fsgn * np.arccos(dp)<EOL>vx_ctr = np.mean(vx, axis=<NUM_LIT:0>)<EOL>(x0, y0) = vx_ctr<EOL>den_inv = zinv(np.sqrt(x0**<NUM_LIT:2> + y0**<NUM_LIT:2>))<EOL>drad_dvx = np.asarray([x0,  y0]) * den_inv<EOL>dtan_dvx = np.asarray([-y0, x0]) * den_inv<EOL>drad_dfx = np.asarray([np.sum(drad_dvx[i]*grad[:,i], axis=<NUM_LIT:0>) for i in [<NUM_LIT:0>,<NUM_LIT:1>]])<EOL>dtan_dfx = np.asarray([np.sum(dtan_dvx[i]*grad[:,i], axis=<NUM_LIT:0>) for i in [<NUM_LIT:0>,<NUM_LIT:1>]])<EOL>rad_mag = zinv(np.sqrt(np.sum(drad_dfx**<NUM_LIT:2>, axis=<NUM_LIT:0>)))<EOL>tan_mag = zinv(np.sqrt(np.sum(dtan_dfx**<NUM_LIT:2>, axis=<NUM_LIT:0>)))<EOL>if to == '<STR_LIT>':<EOL><INDENT>return {'<STR_LIT>': rad_mag, '<STR_LIT>': tan_mag, '<STR_LIT>': arl_mag, '<STR_LIT>': fsgn}<EOL><DEDENT>mtx = simplex_summation_matrix(mesh.tess.indexed_faces)<EOL>cols = np.asarray(mtx.sum(axis=<NUM_LIT:1>), dtype=np.float)[:,<NUM_LIT:0>]<EOL>cols_inv = zinv(cols)<EOL>s_areas = mtx.dot(s_areas)<EOL>v_areas = mtx.dot(v_areas)<EOL>arl_mag = s_areas * zinv(v_areas)<EOL>(rad_mag, tan_mag, fsgn) = [cols_inv * mtx.dot(x) for x in (rad_mag, tan_mag, fsgn)]<EOL>return {'<STR_LIT>': rad_mag, '<STR_LIT>': tan_mag, '<STR_LIT>': arl_mag, '<STR_LIT>': fsgn}<EOL>
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.