Generate an stl file from the Ackley function

import numpy as np
from scipy.spatial import Delaunay

x=np.linspace(-5, 5, 100)
y=np.linspace(-5, 5, 100)
x,y=np.meshgrid(x,y)
x=x.flatten()
y=y.flatten()

# evaluate the parameterization at the flattened x and y
# Ackley function
z = -20*np.exp(-0.2*np.sqrt(0.5*(x**2+y**2)))-np.exp(0.5*(np.cos(2*np.pi*x)+np.cos(2*np.pi*y)))+np.e+20

# define 2D points, as input data for the Delaunay triangulation of U
points2D=np.vstack([x, y]).T
tri = Delaunay(points2D)

points3D = np.hstack((points2D, z.reshape(len(z), -1)))

name = 'ackley'
with open(name+'.stl', 'w') as stl:
    stl.write('solid '+name+'\n')
    for tri in tri.simplices:

        # points that make the triangle
        pts = [points3D[i, :] for i in tri]

        # calculate the normal vector
        n = np.cross( pts[1] - pts[0], pts[2] - pts[0])
        nlen = np.sqrt(n[0]**2 + n[1]**2 + n[2]**2)
        n /= nlen
        stl.write('facet normal {} {} {}\n'.format(*n))

        # write the vertices
        stl.write('    outer loop\n')
        for pt in pts:
            stl.write('        vertex {} {} {}\n'.format(*pt))
        stl.write('    endloop\n')
        stl.write('endfacet\n')

    stl.write('endsolid '+name+'\n')