Surface area of an oblate spheroid using gaussian quadrature 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)

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.
