Trouble understanding Skyfield visualization calculations

I have found this code written in python showing a visualization of the Roadster's orbit around Earth using Skyfield library. I had trouble understanding following piece of code snippets in terms of its physical importance to the visualization.

theta = np.linspace(0, twopi, 201)
cth, sth, zth = [f(theta) for f in (np.cos, np.sin, np.zeros_like)]
lon0 = re*np.vstack((cth, zth, sth))
lons = []
for phi in rads*np.arange(0, 180, 15):
    cph, sph = [f(phi) for f in (np.cos, np.sin)]
    lon = np.vstack((lon0[0]*cph - lon0[1]*sph,
                     lon0[1]*cph + lon0[0]*sph,
                     lon0[2]) )
    lons.append(lon)

lat0 = re*np.vstack((cth, sth, zth))
lats = []
for phi in rads*np.arange(-75, 90, 15):
    cph, sph = [f(phi) for f in (np.cos, np.sin)]
    lat = re*np.vstack((cth*cph, sth*cph, zth+sph))
    lats.append(lat)

Briefly, I am having trouble understanding some "magic" numbers that appear around the code. I would really appreciate if someone can walk me through at least this piece of code if not the full.

Cheers.