In [1]:
from sympy import *
s, t = symbols('s, t')

In [2]:
c = Rational(2, 10)
# x = s + c*sin(pi*s)*sin(pi*t)
# y = t + c*sin(pi*s)*sin(pi*t)

# Flat
# x = s
# y = t

#Chebyshev
# x = -cos(pi*(s+1)/2)
# y = -cos(pi*(t+1)/2)

# Polar
# r = s + 2
# x = r*cos(t)
# y = r*sin(t)

# Squeeze
#x = s + t
#y = t

#Square to circle
# c = Rational(10, 10)
# x = s*(1-c) + c*s*sqrt(1-t**2/2)
# y = t*(1-c) + c*t*sqrt(1-s**2/2)

# u = s + 2
# v = t + 2
# x = u**2 + v**3
# y = u**2 - v**3

#x = (s+2)**3
#y = (t+2)**3

# u = s + 2
# v = t
# x = u**2 - v**2
# y = 2*u*v

# Sine Transform
x = s
y = t - sin(pi*s)

In [3]:
# Convert to C code for easy inclusion in the shaders. 
print(printing.ccode(x, assign_to='x'))
print(printing.ccode(y, assign_to='y'))


x = s;
y = t - sin(M_PI*s);

In [4]:
(x.subs(s, 0).subs(t, 0), 
 x.subs(s,-1).subs(t,-1), 
 y.subs(s, 0).subs(t, 0), 
 y.subs(s,-1).subs(t,-1))


Out[4]:
(0, -1, 0, -1)

In [5]:
J = Matrix([[diff(x, s), diff(x, t)], 
            [diff(y, s), diff(y, t)]])
J


Out[5]:
Matrix([
[            1, 0],
[-pi*cos(pi*s), 1]])

In [6]:
simplify(J.det())


Out[6]:
1

In [7]:
g = J.T*J
simplify(g)


Out[7]:
Matrix([
[pi**2*cos(pi*s)**2 + 1, -pi*cos(pi*s)],
[         -pi*cos(pi*s),             1]])

In [8]:
# Calculate the invariants
I   = g.trace()
II  = (g.trace()**2 - (g*g).trace())/2
III = g.det()

In [9]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.collections as mc
%matplotlib inline
from sympy.utilities.lambdify import lambdify

In [10]:
def cartesian_product(a, b):
    return np.vstack([np.repeat(a, len(b)), np.tile(b, len(a))]).T
cartesian_product([1,2,3], ['a', 'b'])


Out[10]:
array([['1', 'a'],
       ['1', 'b'],
       ['2', 'a'],
       ['2', 'b'],
       ['3', 'a'],
       ['3', 'b']], 
      dtype='<U21')

In [11]:
def plot_grid(n, m, sub, diffeomorphism=(s, t)):
    
    xmin, xmax = -1, +1
    ymin, ymax = -1, +1
    nn = sub*n
    mm = sub*m
    x  = np.linspace(xmin, xmax, n)
    y  = np.linspace(ymin, ymax, m)
    xx = np.linspace(xmin, xmax, nn)
    yy = np.linspace(ymin, ymax, mm)

    plt.figure(figsize=(10,10))
    ax = plt.axes()
    xyy = np.vstack([np.repeat(x, mm), np.tile(yy, n)]).T.reshape(n,mm,2)
    xxy = np.vstack([np.tile(xx, m), np.repeat(y, nn)]).T.reshape(m,nn,2)
    
    varphi = lambdify((s,t), diffeomorphism, 'numpy')
    def add_collection(xy):
        xl, yl = xy[:,:,0], xy[:,:,1]
        xy[:,:,0], xy[:,:,1] = varphi(xl, yl)
        lc = mc.LineCollection(xy, colors='black')
        ax.add_collection(lc)
    
    add_collection(xyy)
    add_collection(xxy)

    xmin = min(xxy[:,:,0].min(), xyy[:,:,0].min())
    xmax = max(xxy[:,:,0].max(), xyy[:,:,0].max())
    ymin = min(xxy[:,:,1].min(), xyy[:,:,1].min())
    ymax = max(xxy[:,:,1].max(), xyy[:,:,1].max())
    ax.axes.set_aspect('equal')
    plt.xlim(xmin-.1, xmax+.1)
    plt.ylim(ymin-.1, ymax+.1)
    plt.axis('off')

In [12]:
N, M = 16, 16
#N, M = 32, 32
#N, M = 64, 64
plot_grid(N, M, 100, (x, y))



In [13]:
plot_grid(N, M, 100)



In [14]:
# Plot a symbolic expression
def plot_st(expr, cmap=plt.cm.seismic, center=None):
    X = Y = np.linspace(-1, +1, 1000)
    X, Y = np.meshgrid(X, Y)
    Z = lambdify((s,t), expr, 'numpy')(X, Y)+ 0*X
    kargs = {}
    if center is not None:
        zmin, zmax = Z.min(), Z.max()
        d = max(center - zmin, zmax - center)
        kargs['vmin'] = center-d
        kargs['vmax'] = center+d
    con = plt.contourf(X, Y, Z, cmap=cmap, **kargs)
    cbar = plt.colorbar(con)

Plot the invariants of the metric


In [15]:
plot_st(I, center=2)



In [16]:
plot_st(II, center=1)



In [17]:
plot_st(III, center=1)


/usr/lib/python3.4/site-packages/matplotlib/colorbar.py:839: RuntimeWarning: invalid value encountered in true_divide
  z = np.take(y, i0) + (xn - np.take(b, i0)) * dy / db
/usr/lib/python3.4/site-packages/matplotlib/colorbar.py:588: RuntimeWarning: invalid value encountered in greater
  inrange = (ticks > -0.001) & (ticks < 1.001)
/usr/lib/python3.4/site-packages/matplotlib/colorbar.py:588: RuntimeWarning: invalid value encountered in less
  inrange = (ticks > -0.001) & (ticks < 1.001)

Plot the components of the metric


In [18]:
plot_st(g[0,0])



In [19]:
plot_st(g[1,1])



In [20]:
plot_st(g[1,0])



In [21]:
plot_st(log(abs(J.det()), 10), cmap=plt.cm.bone)