使用 Legendre 多项式在极坐标中建立和绘制 Matplotlib 2d 直方图答案

【问题标题】：Build up and plot Matplotlib 2d histogram in polar coordinates using the Legendre Polynomials使用 Legendre 多项式在极坐标中建立和绘制 Matplotlib 2d 直方图
【发布时间】：2013-10-09 01:20:20
【问题描述】：

我正在尝试绘制分布：

这是半径为 (a) 的球体内的温度分布，其上半球保持在 T=1，下半球保持在 T=0（忽略两个半球边界处的不连续性），P_l 是第一类勒让德多项式。

import pylab as pl
from scipy.special import eval_legendre as Leg
import math,sys

def sumTerm(a,r,theta,l):
    """ 
    Compute term of sum given radius of sphere (a),
    y and z coordinates, and the current index of the 
    Legendre polynomials (l) over the entire range
    where these polynomials are orthogonal [-1,1].
    """
    xRange = pl.arange(-0.99,1.0,0.01)
    x = pl.cos(theta)
    # correct for scipy handling negative indices incorrectly
    lLow = l-1
    lHigh = l+1
    if lLow < 0:
        lLow = -lLow-1
    return 0.5*((r/a)**l)*Leg(l,x)*(Leg(lLow,0)-Leg(lHigh,0))

def main():

    n = 10      # number of l terms to expand to
    a = 1.0     # radius of sphere

    # generate r, theta values
    aBins = pl.linspace(0, 2*pl.pi, 360)      # 0 to 360 in steps of 360/N.
    rBins = pl.linspace(0, 1, 50)
    theta,r = pl.meshgrid(aBins, rBins)

    tempProfile = pl.zeros([50,360])
    for nr,ri in enumerate(rBins):
        for nt,ti in enumerate(aBins):
            temp = 0.0
            for l in range(n):
                temp += sumTerm(a, ri, ti, l)
            tempProfile[nr,nt] = temp

    # plot the Temperature profile
    pl.imshow(tempProfile)
    pl.colorbar()
    pl.axes().set_aspect('equal')
    pl.show()

if __name__=='__main__':
    main()

这会产生以下情节：

这看起来不错，但是如何在极坐标中显示呢？

【问题讨论】：

您的代码没有运行。您需要 import numpy as np 并在某处定义 gs
是的，对不起，这是用于在笛卡尔坐标系中绘制事物的代码的合并（不是我想要的）以及我在 SO 上可以找到的在极坐标中绘制二维直方图的代码。绘图调用也会引发错误。我为草率的代码道歉（令人尴尬），但我认为有人可能马上就知道如何做类似的事情。我也很乐意接受有人建立任何其他发行版的任何示例。

标签： python matplotlib histogram2d

【解决方案1】：

好的，所以我想通了。这是我的解决方案（我觉得自己的解决方案很奇怪）。

# =============================================================================
# Plot central cross-section of sphere under steady-state conditions
# where the temperature on upper hemisphere is T=T_0 and the lower 
# hemisphere is held at T=0.  This is an expansion in Legendre polynomials.
#
# Author:           Max Graves
# Last Revised:     8-OCT-2013
# =============================================================================

import pylab as pl

from scipy.special import eval_legendre as Leg
import math,sys

def sumTerm(a,r,theta,l):
    """ 
    Compute term of sum given radius of sphere (a),
    y and z coordinates, and the current index of the 
    Legendre polynomials (l) over the entire range
    where these polynomials are orthogonal [-1,1].
    """
    xRange = pl.arange(-0.99,1.0,0.01)
    x = pl.cos(theta)
    # correct for scipy handling negative indices incorrectly
    lLow = l-1
    lHigh = l+1
    if lLow < 0:
        lLow = -lLow-1
    return 0.5*((r/a)**l)*Leg(l,x)*(Leg(lLow,0)-Leg(lHigh,0))

def main():

    n = 20      # number of l terms to expand to
    a = 1.0     # radius of sphere

    # generate r, theta values
    aBins = pl.linspace(0, 2*pl.pi, 360)      # 0 to 360 in steps of 360/N.
    rBins = pl.linspace(0, 1, 50)
    theta,r = pl.meshgrid(aBins, rBins)

    tempProfile = pl.zeros([50,360])
    for nr,ri in enumerate(rBins):
        print nr
        for nt,ti in enumerate(aBins):
            temp = 0.0
            for l in range(n):
                temp += sumTerm(a, ri, ti, l)
            tempProfile[nr,nt] = temp

    # plot the Temperature profile
    fig, ax = pl.subplots(subplot_kw=dict(projection='polar'))
    pax = ax.pcolormesh(theta, r, tempProfile)
    ax.set_theta_zero_location("N") # 'north' location for theta=0
    ax.set_theta_direction(-1)      # angles increase clockwise
    fig.colorbar(pax)

    pl.show()

if __name__=='__main__':
    main()

产生以下情节：

【讨论】：