2

我的目标是制作二维球体的密度热图。当我使用矩形域时，该行下方的绘图代码有效。但是，我正在尝试将代码用于循环域。球体半径为 1。到目前为止我的代码是：

from pylab import *
import numpy as np
from matplotlib.colors import LightSource
from numpy.polynomial.legendre import leggauss, legval

xi = 0.0
xf = 1.0
numx = 500

yi = 0.0
yf = 1.0
numy = 500


def f(x):
    if 0 <= x <= 1:
        return 100
    if -1 <= x <= 0:
        return 0


deg = 1000
xx, w = leggauss(deg)
L = np.polynomial.legendre.legval(xx, np.identity(deg))
integral = (L * (f(x) * w)[None,:]).sum(axis = 1)

c = (np.arange(1, 500) + 0.5) * integral[1:500]


def r(x, y):
    return np.sqrt(x ** 2 + y ** 2)


theta = np.arctan2(y, x)
x, y = np.linspace(0, 1, 500000)


def T(x, y):
    return (sum(r(x, y) ** l * c[:,None] *
                np.polynomial.legendre.legval(xx, identity(deg)) for l in range(1, 500)))

T(x, y)应该等于c系数的总和乘以作为函数的半径x和幂乘以勒让德多项式y的l函数，其中参数是勒让德多项式cos(theta)。

在python:integration a piecewise function中，我学习了如何在求和中使用 Legendre 多项式，但该方法略有不同，对于绘图，我需要一个函数T(x, y)。

这是绘图代码。

densityinterpolation = 'bilinear'
densitycolormap = cm.jet
densityshadedflag = False
densitybarflag = True
gridflag = True
plotfilename = 'laplacesphere.eps'

x = arange(xi, xf, (xf - xi) / (numx - 1))
y = arange(yi, yf, (yf - yi) / (numy - 1))
X, Y = meshgrid(x, y)
z = T(X, Y)


if densityshadedflag:
    ls = LightSource(azdeg = 120, altdeg = 65)
    rgb = ls.shade(z, densitycolormap)
    im = imshow(rgb, extent = [xi, xf, yi, yf], cmap = densitycolormap)
else:
    im = imshow(z, extent = [xi, xf, yi, yf], cmap = densitycolormap)
im.set_interpolation(densityinterpolation)
if densitybarflag:
    colorbar(im)


grid(gridflag)

show()

我在 Mathematica 中制作了情节以供参考我的最终目标是什么 在此处输入图像描述

4

1 回答 1

0

如果您将磁盘域（或您想要的任何域）之外的值设置为 float('nan')，则在绘图时将忽略这些点（将它们保留为白色）。

于 2015-09-09T20:50:06.760 回答