Matplotlib plot_surface spurious faces

These equations

Xs = radius * np.outer(np.cos(latitudes), np.sin(longitudes))
Ys = radius * np.outer(np.sin(latitudes), np.sin(longitudes))
Zs = radius * np.outer(np.ones(latitudes.size), np.cos(longitudes))

are computing cartesian X,Y,Z coordinates given spherical coordinates for radius, latitude and longitude. But if that is so, then the longitude should range from 0 to pi, not 0 to 2pi. Therefore, change

longitudes = np.linspace(0, 2*np.pi, N_POINTS)

to

longitudes = np.linspace(0, np.pi, N_POINTS)

import math
import numpy as np
import scipy.special as special
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

" Problem setup. "
# m**3/s**2, from EGM96.
GM = 3986004.415E8 
# m, from EGM96.
a = 6378136.3 
# Number of lattitudes and longitudes used to plot the geoid.
N_POINTS = 50 
# Geocentric     latitudes and longitudes where the geoid will be visualised.
latitudes = np.linspace(0, 2*np.pi, N_POINTS) 
longitudes = np.linspace(0, np.pi, N_POINTS)
# Radius at which the equipotential will be computed, m.
radius = 6378136.3 
# Maximum degree of the geopotential to visualise.
MAX_DEGREE = 2 

" EGM96 coefficients - minimal working example. "
Ccoeffs={2:[-0.000484165371736, -1.86987635955e-10, 2.43914352398e-06]}
Scoeffs={2:[0.0, 1.19528012031e-09, -1.40016683654e-06]}

" Compute the gravitational potential at the desired locations. "
# Gravitational potentials computed with the given geoid. Start with ones and
# just add the geoid corrections.
gravitationalPotentials = np.ones( latitudes.shape ) 

# Go through all the desired orders and compute the geoid corrections to the
# sphere.
for degree in range(2, MAX_DEGREE+1): 
    # Contribution to the potential from the current degree and all
    # corresponding orders.
    temp = 0. 
    # Legendre polynomial coefficients corresponding to the current degree.
    legendreCoeffs = special.legendre(degree) 
    # Go through all the orders corresponding to the currently evaluated degree.
    for order in range(degree): 
        temp += (legendreCoeffs[order] 
                 * np.sin(latitudes) 
                 * (Ccoeffs[degree][order]*np.cos( order*longitudes ) 
                    + Scoeffs[degree][order]*np.sin( order*longitudes )))

    # Add the contribution from the current degree.
    gravitationalPotentials = math.pow(a/radius, degree) * temp 

# Final correction.
gravitationalPotentials *= GM/radius 

" FIGURE THAT SHOWS THE SPHERICAL HARMONICS. "
fig = plt.figure(figsize=(12,8))
ax = Axes3D(fig)
ax.set_aspect("equal")
ax.view_init(elev=45., azim=45.)
ax.set_xlim([-1.5*radius, 1.5*radius])
ax.set_ylim([-1.5*radius, 1.5*radius])
ax.set_zlim([-1.5*radius, 1.5*radius])

# Make sure the shape of the potentials is the same as the points used to plot
# the sphere.

# Don't need the second copy of colours returned by np.meshgrid
gravitationalPotentialsPlot = np.meshgrid(
    gravitationalPotentials, gravitationalPotentials )[0] 
# Normalise to [0 1]
gravitationalPotentialsPlot /= gravitationalPotentialsPlot.max() 

" Plot a sphere. "
Xs = radius * np.outer(np.cos(latitudes), np.sin(longitudes))
Ys = radius * np.outer(np.sin(latitudes), np.sin(longitudes))
Zs = radius * np.outer(np.ones(latitudes.size), np.cos(longitudes))
equipotential = ax.plot_surface(
    Xs, Ys, Zs, facecolors=plt.get_cmap('jet')(gravitationalPotentialsPlot), 
    rstride=1, cstride=1, linewidth=0, antialiased=False, shade=False)

plt.show()

yields