Gravitační pole prstence - osa kolmá k rovině prstence

výpočet intenzity a potenciálu gravitačního pole nekonečně tenkého homogenního prstence o poloměru $R$ a hmotnosti $m$

import numpy as np
import matplotlib.pyplot as plt

Prstenec leží v rovine $xy$

Zvolme takové jednotky, že $G = 1$

Parametry prstence:

• poloměr $R = 1$
• hmotnost $m = 1$

Funkce počítající potenciál gravitačního pole v ose prstence ležící v rovině prstence

#potencial
def pot(x,y,z):
    N=1000
    dalpha=2*np.pi/N
    dm=1/N
    sum=0.0
    R=1
    alpha=0
    for i in range(N):
        alpha=alpha+dalpha
        r=np.sqrt((R*np.cos(alpha)-x)**2+(R*np.sin(alpha)-y)**2+z**2)
        if r>0:
            sum=sum-1/r*dm
    return(sum)

Funkce počítající intenzitu gravitačního pole v rovině prstence

#x-ova slozka intenzity
def Kx(x,y,z):
    N=1000
    dalpha=2*np.pi/N
    dm=1/N
    sum=0.0
    R=1
    alpha=0
    for i in range(N):
        alpha=alpha+dalpha
        r=np.sqrt((R*np.cos(alpha)-x)**2+(R*np.sin(alpha)-y)**2+z**2)
        if r>0:
            sum=sum-1/r**3*(x-R*np.cos(alpha))*dm
    return(sum)

#y-ova slozka intenzity
def Ky(x,y,z):
    N=1000
    dalpha=2*np.pi/N
    dm=1/N
    sum=0.0
    R=1
    alpha=0
    for i in range(N):
        alpha=alpha+dalpha
        r=np.sqrt((R*np.cos(alpha)-x)**2+(R*np.sin(alpha)-y)**2+z**2)
        if r>0:
            sum=sum-1/r**3*(y-np.sin(alpha))*dm
    return(sum)

#z-ova slozka intenzity
def Kz(x,y,z):
    N=1000
    dalpha=2*np.pi/N
    dm=1/N
    sum=0.0
    R=1
    alpha=0
    for i in range(N):
        alpha=alpha+dalpha
        r=np.sqrt((R*np.cos(alpha)-x)**2+(R*np.sin(alpha)-y)**2+z**2)
        if r>0:
            sum=sum-1/r**3*z*dm
    return(sum)

Vektorizujeme funkce pot, Kx, Ky, Kz abychom s nimi mohli počítat s celým polem najednou

vpot=np.vectorize(pot)
vKx=np.vectorize(Kx)
vKy=np.vectorize(Ky)
vKz=np.vectorize(Kz)

Zvolíme body na ose, kde budeme intenzitu a potenciál počítat

z=np.arange(-5,5,0.01)

Provedeme vlastní výpočet pomocí vektorizovaných funkcí vpota vKz
Složky vektoru intenzity $K_x$ a $K_y$ jsou na ose kolmé na rovinu prstence nulové, tj. vektor $\vec{K}$ má tvar $(0,0,K_z)$.

p=vpot(0,0,z)
K=vKz(0,0,z)

Nyní můžeme vykreslit průběh potenciálu.

ig,ax=plt.subplots()
ax.plot(z,p)
ax.plot([-1/np.sqrt(2),-1/np.sqrt(2)],[-10,10],ls='dashed',c='grey')
ax.plot([1/np.sqrt(2),1/np.sqrt(2)],[-10,10],ls='dashed',c='grey')
plt.ylim(-3,0)
plt.xlim(-5,5)
plt.title('potenciál',fontsize=14)
ax.set_xlabel('z',fontsize=14)
ax.set_ylabel('$ \phi $',fontsize=14)
plt.show()

a průběh $z$-tové složky intenzity

fig,ax=plt.subplots()
ax.plot(z,K)
ax.plot([-1/np.sqrt(2),-1/np.sqrt(2)],[-10,10],ls='dashed',c='grey')
ax.plot([1/np.sqrt(2),1/np.sqrt(2)],[-10,10],ls='dashed',c='grey')
ax.plot([-5,5],[-0,0],ls='dashed',c='grey')
plt.ylim(-1,1)
plt.xlim(-5,5)
plt.title('z-tová složka intenzity',fontsize=14)
ax.set_xlabel('z',fontsize=14)
ax.set_ylabel('$ K_z $',fontsize=14)
plt.show()