Python: график плотности тепла в диске

Question

Моя цель - сделать график плотности плотности тепла сферы в 2D. Код построения под строкой работает, когда я использую прямоугольные домены. Тем не менее, я пытаюсь использовать код для круговой области. Радиус сферы 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 времени мощности полинома Лежандра где аргумент из полинома legendre равен cos(theta).

В python: integrating a piecewise function я научился использовать полиномы Лежандра в суммировании, но этот метод несколько отличается, и для построения графика мне нужна функция 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 для справки о том, что моя конечная цель

Answer 1

Если установить значения вне области диска (или в зависимости от того домена вы хотите), чтобы плавать («нан»), эти точки будут игнорироваться при построении (оставляя их белым цветом).