For me it seems to work to avoid the zero value for theta, I’m not sure why this breaks the normals, maybe it’s a bug in the automatic normal mesh generation, because the other pole doesn’t have the problem.
If I make that zero value a small positive value instead, I get this:
using GLMakie
GLMakie.activate!()
GLMakie.closeall()
N = 16
nlat = N + 1
nlon = 2nlat
theta = range(0.001, π, length=nlat)
phi = range(0, 2π, length=nlon+1)
Xi = Array{Float64}(undef, nlat, nlon+1, 3)
# Compute the RBC shape in spherical coordinates
alpha = 1.3858189
Xi[:,:,1] .= alpha .* (sin.(theta) .* cos.(phi)')
Xi[:,:,2] .= alpha .* (sin.(theta) .* sin.(phi)')
Xi[:,:,3] .= (0.5 * alpha) .* (0.207 .+ 2.003 .* sin.(theta).^2 .- 1.123 .* sin.(theta).^4) .* cos.(theta) .* ones(1, nlon+1)
fig = Figure(size = (1200, 800))
axis = Axis3(fig[1, 1], aspect = :data)
surface!(axis, Xi[:,:,1], Xi[:,:,2], Xi[:,:,3]; colormap=[:red, :red], transparency=true, backlight=1.0f0, shading=MultiLightShading)
wireframe!(axis, Xi[:,:,1], Xi[:,:,2], Xi[:,:,3], color=:black, linewidth=1, transparency=true)
fig
Another option would be to generate and specify the normals yourself and use a mesh instead of a surface.