Surface area of an oblate spheroid using gaussian quadrature

Question

I want to compute the surface area of an oblate spheroid using gaussian quadrature, the parametrization of the oblate spheroid is given by:

$$x = a \cdot \sin\theta \cdot \cos \phi \\ y = a \cdot \sin\theta \cdot \sin \phi \\ z = b \cdot \cos\theta $$

Where $b<a$, in order to compute this integral I am using a Gauss-legendre quadrature to compute the points and the weights for the integral for the first octant on the spheroid so $\theta = [0,\pi/2]$, $\phi = [0,\pi/2]$. So using the weights and the nodes of a Gauss-Legendre quadrature of order $n$ I can define the weights and the nodes in $\theta$ and $\phi$ by:

$$X_{\theta/\phi} = \frac12 \cdot(X_{\text{Gauss-Leg}} + 1)\cdot \frac{\pi}2 \\ W_{\theta/\phi} = \frac12\cdot W_{\text{Gauss-Leg}} \cdot\frac{\pi}2$$

So having this defined, I can compute the surface integral as:

$$\int_{0}^{\pi/2}\int_{0}^{\pi/2} dS = \sum\sum W_{\theta}\cdot W_{\phi}\cdot dS $$

Where I think $dS = a^2\cdot b \cdot\rho^2 \cdot \sin\theta \cdot d\theta \cdot d\phi$. I find a bit confusing the $\rho$ here but I have tried using the standard definition of $\rho$ as:

$$\rho = \sqrt{x^2 + y^2 + z^2} = \sqrt{a^2\cdot \sin^2\theta + b^2 \cdot \cos^2\theta} $$

So the expresion that I am using for solving this integral is:

$$\int_{0}^{\pi/2}\int_{0}^{\pi/2} dS = \sum\sum W_{\theta}\cdot W_{\phi}\cdot a^2\cdot b \cdot(a^2\cdot \sin^2\theta + b^2 \cdot \cos^2\theta) \cdot \sin\theta $$

So I am doing something wrong since the surface of the oblate spheroid can be computed using:

$$S = 2\pi \cdot \left(a^2 + \frac{b^2}{\sin(ae)} \ln\Bigl(\frac{1 + \sin(ae)}{\cos(ae)} \Bigr) \right)$$

with $ae = \arccos(b/a)$. And the obtained result is not the same than the analytical surface. I have made an small python script that computes both results and prints them:

import numpy as np
from scipy.special import roots_legendre

#Define a and b
b = 2.
a = 100.

#Compute the Weights and nodes

x_phi, w_phi = roots_legendre(150)
x_theta, w_theta = roots_legendre(100)

#Translate them
x_phi = 0.5 * (x_phi + 1.) * np.pi/2.
x_theta = 0.5 * (x_theta + 1.) * np.pi/2.
w_phi = 0.5 * w_phi * np.pi/2.
w_theta = 0.5 * w_theta * np.pi/2.

#Compute the integral

integral = 0
for i in xrange(len(x_phi)):
    for j in xrange(len(x_theta)):
        integral += w_phi[i] * w_theta[j] * a**2 * b  * (a**2 * np.sin(x_theta[j])**2 + b**2 * np.cos(x_theta[j])**2) * np.sin(x_theta[j])


print("Estimated int: %f" %(8*integral))

ae = np.arccos(b/a)
surface = 2*np.pi*(a**2 + b**2/np.sin(ae) * np.log((1+np.sin(ae))/np.cos(ae)))

print("Real int: %f" %(surface))

So what am I doing wrong? (I have to mention that this is just a simple test, what I really want to do is to compute the surface integral of any arbitrary function on this spheroid)

Answer 1

I found the error, it was on the element of surface dS, here is the working code:

import numpy as np
from scipy.special import roots_legendre

#Define a and b
b = 2.
a = 100.

#Compute the Weights and nodes

x_phi, w_phi = roots_legendre(150)
x_theta, w_theta = roots_legendre(100)

#Translate them
x_phi = 0.5 * (x_phi + 1.) * np.pi/2.
x_theta = 0.5 * (x_theta + 1.) * np.pi/2.
w_phi = 0.5 * w_phi * np.pi/2.
w_theta = 0.5 * w_theta * np.pi/2.

#Compute the integral

integral = 0
for i in xrange(len(x_phi)):
    for j in xrange(len(x_theta)):
        integral += w_phi[i] * w_theta[j] * np.sqrt(b**2/a**2 * np.tan(x_theta[j])**2 + 1) * a**2 * np.sin(x_theta[j]) * np.cos(x_theta[j])

print("Estimated int: %f" %(8*integral))

ae = np.arccos(b/a)
surface = 2*np.pi*(a**2 + b**2/np.sin(ae) * np.log((1+np.sin(ae))/np.cos(ae)))

print("Real int: %f" %(surface))
```

Answer 2

Even though you eventually got the right areal element you seem to have an unsteady method of deriving it. Since on the surface we have parameterized $$\vec r=\langle x,y,z\rangle=\langle a\sin\theta\cos\phi,a\sin\theta\sin\phi,b\cos\theta\rangle$$ We can take the differential to get $$d\vec r=\langle a\cos\theta\cos\phi,a\cos\theta\sin\phi,-b\sin\theta\rangle\,d\theta+\langle-a\sin\theta\sin\phi,a\sin\theta\cos\phi,0\rangle\,d\phi$$ We can form the cross product to get the areal element $$\begin{align}d^2\vec A&=\pm\langle a\cos\theta\cos\phi,a\cos\theta\sin\phi,-b\sin\theta\rangle\,d\theta\times\langle-a\sin\theta\sin\phi,a\sin\theta\cos\phi,0\rangle\,d\phi\\ &=\pm\langle ab\sin^2\theta\cos\phi,ab\sin^2\theta\sin\phi,a^2\sin\theta\cos\theta\rangle\,d\theta\,d\phi\end{align}$$ And we could take the $+$ sign to get the outward normal. Then we can take its magnitude to get the scalar areal element: $$d^2A=\lVert d^2\vec A\rVert=\sqrt{b^2\sin^2\theta+a^2\cos^2\theta}\,a\sin\theta\,d\theta\,d\phi$$ It's equivalent to the areal element you eventually arrived at but I don't see how you could have gotten such a convoluted expression for it.