问三角形等值线网格划分

EN

#!/usr/bin/env python
from __future__ import division
import matplotlib
matplotlib.use("TkAgg")
matplotlib.rc('text', usetex=True)
matplotlib.rcParams['text.latex.preamble']=r"""\usepackage{amsmath}
"""
import math
import scipy.special

root_three_over_two = np.sqrt(3) / 2

def new_figure():
    # 1.45
    plt.figure(figsize = [2.6, 2.6 * root_three_over_two], dpi = 1200)
    plt.axes([0.05, 0.10, 0.90, 0.90], frameon = False)

xsize = 1.0
ysize = root_three_over_two * xsize
plt.axis([0, xsize, 0, ysize])
resolution = 0.05

R = inclusive_arange(0.0, 1.0, resolution)
x, y = np.meshgrid(inclusive_arange(0.0, 1.0, resolution),
                   inclusive_arange(0.0, 1.0, resolution))
# UNFORTUNATELY x, and y include a lot of points where x+y>1
x = []
y = []
for yy in R:
    x.append(list(inclusive_arange(0.0, 1.0 - yy, resolution)))
    y.append([yy for xx in R])
print x
print y
z = 1 - x - y

# We can use these to convert to and from the equilateral triangle.
M = [[1, 0.5], [0, root_three_over_two]]
Mi = np.linalg.inv(M)

def dirichlet(x, y, z, a, b, c):
    if z < 0:
        return 0
    return x ** (a - 1) * y ** (b - 1) * z ** (c - 1) \
            * math.gamma(a + b + c) \
            / (math.gamma(a) * math.gamma(b) * math.gamma(c))

dirichlet = np.frompyfunc(dirichlet, 6, 1)

for (dirichlet_parm, filename) in [((5.0, 1.5, 2.5), "dir_small.pdf")]:
    new_figure()
    height = dirichlet(x, y, z, *dirichlet_parm)
    M = np.max(height)
    cs = plt.contourf(x, y, height, 50)
    S = sum(dirichlet_parm)
    plt.savefig(filename)