Gravitační pole prstence - osa ležící v 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

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

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

p=vpot(x,0,0)
K=vKx(x,0,0)

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

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

a průběh $x$-ové složky intenzity

fig,ax=plt.subplots()
ax.plot(x,K)
ax.plot([-1,-1],[-10,10],ls='dashed',c='grey')
ax.plot([1,1],[-10,10],ls='dashed',c='grey')
ax.plot([-5,5],[-0,0],ls='dashed',c='grey')
plt.ylim(-10,10)
plt.xlim(-5,5)
plt.title('x-ová složka intenzity',fontsize=14)
ax.set_xlabel('x',fontsize=14)
ax.set_ylabel('$ K_x $',fontsize=14)
plt.show()