No que se segue, serão dadas as linhas gerais para o cálculo do potencial gravítico causado por um paralelepípedo de densidade \(\rho\). Este método seguirá de perto aquele apresentado no artigo Exterior gravitation of a polyhedron derived an compared with harmonic and mascon gravitation representations of Asteroid 4769 Castalia. Não é difícil concluir que o pontencial gravítico externo causado por um paralelepípedo de matéria é dado pelo integral triplo
\[V=G\rho\int_{-l_1}^{l_1}\int_{-l_2}^{l_2}\int_{-l_3}^{l_3}\frac{dz'dy'dx'}{\sqrt{\left(x-x'\right)^2+\left(y-y'\right)^2+\left(z-z'\right)^2}}\]
Aqui \(G\) representa a constante de gravitação e \(\rho\), a densidade de matéria, assumida como constante ao longo do paralelepípedo. Aplica-se a substituição
\[\left\lbrace\begin{array}{l}x-x'=\xi\\ y-y'=\eta\\ z-z'=\zeta\end{array}\right.\]
de modo a que o integral anterior possa ser colocado na forma
\[\int_{x-l_1}^{x+l_1}\int_{y-l_2}^{y+l_2}\int_{z-l_3}^{z+l_3}\frac{d\zeta d\eta d\xi }{\sqrt{\xi^2+\eta^2+\zeta^2}}\]
Denota-se por \(r=\sqrt{\xi^2+\eta^2+\zeta^2}\) e define-se o vector
\[\vec{r}=\frac{1}{r}\left(\xi,\eta,\zeta\right)\]
Não é difícil verificar que
\[\left\lbrace\begin{array}{l}\frac{\partial}{\partial\xi}\left(\frac{\xi}{r}\right)=\frac{1}{r}-\frac{\xi^2}{r^3}\\ \frac{\partial}{\partial\eta}\left(\frac{\eta}{r}\right)=\frac{1}{r}-\frac{\eta^2}{r^3}\\ \frac{\partial}{\partial\zeta}\left(\frac{\zeta}{r}\right)=\frac{1}{r}-\frac{\zeta^2}{r^3}\end{array}\right.\]
e, portanto,
\[\nabla\vec{r}=\frac{\partial}{\partial\xi}\left(\frac{\xi}{r}\right)+\frac{\partial}{\partial\eta}\left(\frac{\eta}{r}\right)+\frac{\partial}{\partial\zeta}\left(\frac{\zeta}{r}\right)=\frac{2}{r}\]
O integral escreve-se como
\[V=\frac{1}{2}G\rho\int_{x-l_1}^{x+l_1}\int_{y-l_2}^{y+l_2}\int_{z-l_3}^{z+l_3}\left(\frac{\partial}{\partial\xi}\left(\frac{\xi}{r}\right)+\frac{\partial}{\partial\eta}\left(\frac{\eta}{r}\right)+\frac{\partial}{\partial\zeta}\left(\frac{\zeta}{r}\right)\right)\]
Poder-se-ia aqui recorrer ao teorema da divergência para reduzir o integral de volume a um integral de superfície. No entanto, será aqui seguida uma abordagem diferente. O integral divide-se em três parcelas, nomeadamente,
\[ I=I_1+I_2+I_3\]
onde
\[\begin{array}{l}I_1=\int_{x-l_1}^{x+l_1}\int_{y-l_2}^{y+l_2}\int_{z-l_3}^{z+l_3}\frac{\partial}{\partial\xi}\left(\frac{\xi}{r}\right)d\xi d\eta d\zeta\\ I_2=\int_{x-l_1}^{x+l_1}\int_{y-l_2}^{y+l_2}\int_{z-l_3}^{z+l_3}\frac{\partial}{\partial\eta}\left(\frac{\eta}{r}\right)d\xi d\eta d\zeta\\ I_3=\int_{x-l_1}^{x+l_1}\int_{y-l_2}^{y+l_2}\int_{z-l_3}^{z+l_3}\frac{\partial}{\partial\zeta}\left(\frac{\zeta}{r}\right)d\xi d\eta d\zeta\end{array}\]
Trocando a ordem de integração e aplicando o teorema fundamental do cálculo, vem
\[\begin{array}{l}I_1=\frac{1}{2}G\rho\int_{y-l_2}^{y+l_2}\int_{z-l_3}^{z+l_3}\left\lbrack \frac{\xi}{r}\right\rbrack_{\xi=x-l_1}^{\xi=x+l_1} d\zeta d\eta\\ I_2=\frac{1}{2}G\rho\int_{x-l_1}^{z+l_1}\int_{z-l_3}^{z+l_3}\left\lbrack \frac{\eta}{r}\right\rbrack_{\eta=y-l_2}^{\eta=y+l_2} d\zeta d\xi\\ I_3=\frac{1}{2}G\rho\int_{x-l_1}^{x+l_1}\int_{y-l_2}^{y+l_2}\left\lbrack \frac{\zeta}{r}\right\rbrack_{\zeta=z-l_3}^{\zeta=z+l_3} d\zeta d\eta\end{array}\]
Observando que
\[\frac{\partial}{\partial\eta}\left(\frac{\eta}{r}\right)+\frac{\partial}{\partial\zeta}\left(\frac{\zeta}{r}\right)=\frac{2}{r}-\frac{\eta^2+\zeta^2}{r}=\frac{1}{r}-\frac{\xi^2}{r^3}\]
o integral \(I_1\) assume a forma
\[\begin{array}{l}I_1=\frac{1}{2}G\rho\int_{y-l_2}^{y+l_2}\int_{z-l_3}^{z+l_3}\left\lbrack\xi\left(\frac{\partial}{\partial\eta}\left(\frac{\eta}{r}\right)+\frac{\partial}{\partial\zeta}\left(\frac{\zeta}{r}\right)\right)\right\rbrack_{\xi=x-l_1}^{\xi=x+l_1}d\zeta d\eta+\\ +\frac{1}{2}G\rho\int_{y-l_2}^{y+l_2}\int_{z-l_3}^{z+l_3}\left\lbrack\frac{\xi^3}{r^3}\right\rbrack_{\xi=x-l_1}^{\xi=x+l_1}d\zeta d\eta\end{array}\]
o qual, após aplicação do teorema fundamental do cálculo, se pode reduzir à soma de três parcelas
\[I_1=I_{11}+I_{12}+I_{13}\]
onde
\[\begin{array}{l}I_{11}=\int_{z-l_3}^{z+l_3}\left\lbrack\frac{\eta}{r}\right\rbrack_{\begin{array}{l}\eta=y-l_2\\ \xi=x-l_1\end{array}}^{\begin{array}{l}\eta=y+l_2\\ \xi=x+l_1\end{array}}d\zeta\\ I_{12}=\int_{y-l_2}^{y+l_2}\left\lbrack\frac{\zeta}{r}\right\rbrack_{\begin{array}{l}\zeta=z-l_3\\ \xi=x-l_1\end{array}}^{\begin{array}{l}\zeta=z+l_3\\ \xi=x+l_1\end{array}}d\eta\\ I_{13}=G\rho\int_{y-l_2}^{y+l_2}\int_{z-l_3}^{z+l_3}\left\lbrack \frac{\xi^3}{r^3}\right\rbrack_{\xi=x-l_1}^{\xi=x+l_1}d\zeta d\eta\end{array}\]
Os integrais \(I_{11}\) e \(I_{12}\) são fáceis de determinar. Por exemplo,
\[I_{11}=\frac{1}{2}G\rho=\left\lbrack\frac{1}{\sqrt{\xi^2+\eta^2}}\arctan{\left(\frac{\zeta}{\sqrt{\xi^2+\eta^2}}\right)}\right\rbrack_{\begin{array}{l}\xi=x-l_1\\ \eta=y-l_2\\ \zeta=z-l_3\end{array}}^{\begin{array}{l}\xi=x+l_1\\ \eta=y+l_2\\ \zeta=z+l_3\end{array}}\]
O integral \(I_{12}\) é, em tudo, semelhante. Resta determinar o integral de superfície
\[I_{13}=\frac{1}{2}G\rho\int_{y-l_2}^{y+l_2}\int_{z-l_3}^{z+l_3}\left\lbrack \frac{\xi^3}{\left(\xi^2+\eta^2+\zeta^2\right)^{\frac{3}{2}}}\right\rbrack_{\xi=x-l_1}^{\xi=x+l_1}d\zeta d\eta\]
Efectua-se a transformação
\[\lambda=\zeta+\sqrt{\xi^2+\eta^2+\zeta^2}\]
de onde,
\[\lambda^2-2\lambda\zeta-\xi^2-\eta^2=0\]
ou
\[\zeta=\frac{\lambda^2-\xi^2-\eta^2}{2\lambda}\]
e também
\[d\zeta=\frac{\lambda^2+\xi^2+\eta^2}{2\lambda^2}d\lambda\]
bem como
\[\sqrt{\xi^2+\eta^2+\zeta^2}=\lambda-\zeta=\frac{\lambda^2+\xi^2+\eta^2}{2\lambda}\]
A substituição no integral permite obter
\[\begin{array}{l}I_{13}=\frac{1}{2}G\rho\int_{y-l_2}^{y+l_2}\int_{z-l_3}^{z+l_3}\left\lbrack \frac{4\xi^3\lambda}{\left(\lambda^2+\xi^2+\eta^2\right)^2}\right\rbrack_{\xi=x-l_1}^{\xi=x+l_1}d\zeta d\eta=\\ =-G\rho\int_{y-l_2}^{y+l_2}\left\lbrack\frac{\xi^3}{\lambda^2+\xi^2+\eta^2}\right\rbrack_{\begin{array}{l}\xi=x-l_1\\ \lambda=l_3'\end{array}}^{\begin{array}{l}\xi=x+l_1\\ \lambda=l_3''\end{array}}d\eta\end{array}\]
A aplicação da transformação inversa permite obter o integral
\[I_{13}=\int_{y-l_2}^{y+l_2}\left\lbrack\frac{\xi^3\zeta}{\left(\xi^2+\eta^2\right)\sqrt{\xi^2+\eta^2+\zeta^2}}\right\rbrack_{\begin{array}{l}\xi=x-l_1\\ \zeta=z-l_3\end{array}}^{\begin{array}{l}\xi=x+l_1\\ \zeta=z+l_3\end{array}}d\eta\]
O mesmo método acima aplicado permite determinar
\[I_{13}=\left\lbrack\xi^2\arctan{\frac{\eta\sqrt{\xi^2+\eta^2+\zeta^2}+\xi^2+\eta^2}{\xi\zeta}}\right\rbrack_{\begin{array}{l}\xi=x-l_1\\ \eta=y-l_2\\ \zeta=z-l_3\end{array}}^{\begin{array}{l}\xi=x+l_1\\ \eta=y+l_2\\ \zeta=z+l_3\end{array}}\]
Os demais integrais são em tudo semelhantes, obtendo-se uma expressão algo elaborada para o valor do potencial do cubo.
