pcolormeshでオフセット、contour（Python matplotlib）

Question

私はpcolormeshで2D配列をプロットしようとしています。これは、ブラックホール磁気圏の2D軸対称球面プロットです。問題はなぜシータにオフセットがあるのか​​？...通常、磁力線は赤道のブラックホールによって垂直になり、わずかに湾曲していなければなりません。地図も相殺されているようだ。面白いかもしれない私のコードのここ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で、オフセット減少シータより多くのポイントを取るが、まだここBlack hole magnetosphere (256 points)なります。理由の説明ができる人がいれば、非常に役に立ちます。事前にありがとう。

あなたが統合を行うためにintegrate.simps()を使用しているジェレミー

Answer 1

、それは正確に取得するために多くのポイントを必要としています。

あなたがここにコードがある、統合を行うために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") 

と出力：