與pcolormesh，等高線（Python matplotlib）偏移

Question

我試圖用pcolormesh繪製二維數組。這是一個黑洞磁層的二維徑向球面圖。問題是θ在θ上有偏移的原因是什麼？通常情況下，磁場線必須是垂直的，並且在赤道處由黑洞稍微彎曲。地圖也似乎被抵消了。 Black hole magnetosphere (64 points)這部分我的代碼可能是有趣：

import numpy as np, matplotlib.pyplot as plt 
import matplotlib.animation as animation 
import scipy.integrate as integrate 

##### Natural unities 
M = 1.0 
G = 1.0 
c = 1.0 

##### Gravitational radius 
rg = (G * M)/(c*c) 

##### spin 
a = 0.0 

##### Horizon radius 
rh = rg + np.sqrt(rg*rg - a*a) 

##### r, theta parameters 
rmin = 0.9*rh 

rmax = 5.0 

thmin = 0.001*np.pi 
thmax = 0.999*np.pi 

Nr = 64 
Nth = 64 

r = np.logspace(np.log10(rmin),np.log10(rmax),Nr) 

th = np.linspace(thmin,thmax,Nth) 

r_grid, th_grid = np.meshgrid(r,th) 


x = r_grid*np.cos(th_grid) 
y = r_grid*np.sin(th_grid) 

##### ergosphere 

rerg = 1.0 + np.sqrt(1.0-a*a*np.cos(th)*np.cos(th)) 

xerg = rerg*np.cos(th) 
yerg = rerg*np.sin(th) 

##### Horizon 
xrh = rh*np.cos(th) 
yrh = rh*np.sin(th) 


############################# schwarzschild's metric 


Alpha_sch = 1.0/np.sqrt(1.0+(2.0/r_grid)) 

sqr_det_gamma_sch = np.sqrt((1.0+2.0/r_grid)*r_grid*r_grid*r_grid*r_grid*np.sin(th_grid)*np.sin(th_grid)) 


Brg = np.zeros((Nth,Nr)) 


############################# Flux function 

Psy = np.zeros((Nth,Nr)) 

V = [2,4,6,8,10,12] 



##### Wald's solution 
Brg = Alpha_sch*np.cos(th_grid) 


for i in range (0,Nr): 
    for j in range(1,Nth): 
     th3 = th[0:j] 
     Psy[j,i]=integrate.simps(sqr_det_gamma_sch[0:j,i]*Brg[0:j,i],th3) 



plt.figure(figsize=(8,8)) 
ax = plt.subplot(111) 
plt.title('$B^r$') 
circle1 = plt.Circle((0,0),rh,color = 'k') 
ax.add_artist(circle1) 
plt.contour(y,x,Psy,V,colors = 'k') 
plt.pcolormesh(y,x,Brg,cmap='bwr',vmin=-1,vmax=1) 
cbar = plt.colorbar() 
cbar.set_label('Intensité', rotation=270) 
plt.xlim(0,rmax) 
plt.xlabel('$r_g$') 
plt.ylabel('$r_g$') 
plt.legend() 
plt.show() 

當然，如果我在R和THETA偏移減少花費更多點，但還是在這裏Black hole magnetosphere (256 points)。如果你們中的一些人能解釋我爲什麼，這將是非常有益的。先謝謝你。

傑里米

Answer 1

您正在使用integrate.simps()做整合，它需要很多點來獲得精確。

您可以使用integrate.quad()做整合，這裏是代碼：

import numpy as np, matplotlib.pyplot as plt 
import scipy.integrate as integrate 
from math import sqrt, sin, cos 

##### Natural unities 
M = 1.0 
G = 1.0 
c = 1.0 

##### Gravitational radius 
rg = (G * M)/(c*c) 

##### spin 
a = 0.0 

##### Horizon radius 
rh = rg + np.sqrt(rg*rg - a*a) 

##### r, theta parameters 
rmin = 0.9*rh 

rmax = 5.0 

thmin = 0 
thmax = np.pi 

Nr = 64 
Nth = 64 

r = np.logspace(np.log10(rmin),np.log10(rmax),Nr) 

th = np.linspace(thmin,thmax,Nth, endpoint=True) 

r_grid, th_grid = np.meshgrid(r,th) 

x = r_grid*np.cos(th_grid) 
y = r_grid*np.sin(th_grid) 

def f(th, r): 
    th = float(th) 
    r = float(r) 
    Alpha_sch = 1.0/sqrt(1.0+(2.0/r)) 
    sqr_det_gamma_sch = sqrt((1.0+2.0/r)*r*r*r*r*sin(th)*sin(th)) 
    Brg = Alpha_sch*cos(th) 
    return sqr_det_gamma_sch * Brg 

Psy = np.zeros_like(x) 
for i in range (0,Nr): 
    for j in range(Nth): 
     Psy[j, i] = integrate.quad(f, 0, th[j], args=(r[i],))[0] 

fig, ax = plt.subplots(figsize=(8, 8)) 
ax.set_aspect("equal") 
ax.pcolormesh(y, x, Psy) 
ax.contour(y, x, Psy, colors="k") 

和輸出：