Here are the parameters to my equations, with descriptions:
vshift: The gridline that starts at u=0,v=v0 does not connect to itself-- it connects to u=0,v=v0+vshift. This helps prevent gridline collisions at the circle of self-intersection, and it adds a pleasing twist to the grid. This twist is especially visible looking in the "mouth" of the bottle.
side,infrac: These are the number of units the grid should be thickened along the normal and along the original surface, respectively.
skew: This exponent warps the spine of the bottle in the x-axis so help prevent gridline collisions.
perp: This is the coefficient to an adjustment that bends the cross-section such that, at the intersection area, the section with the small radius tilts to become more perpendicular with the section of larger radius.
perp start/end angle: These are the starting and ending u values for the adjustment described above. A quartic function is used to make the effect start and stop smoothly.
baserad, radamp: The cross-section radius varies from baserad-radamp to baserad+radamp.
zpush: This is the coefficient to an expression that pushes the small neck of the bottle out of the plane of symmetry along the same section that is twisted with 'perp'. (zpush=0 in the final shape)
The trickiest part to writing the program was handling the wrap at the end. I tried to write special cases to connect the "outside" to the "inside", etc, to make the grid close correctly. Still, even though I corrected all visible defects, Sara's analyzer told me there were flipped normals and non-manifold edges. Finally, I took advantage of the analyzer's vertex-merging feature: my program does not imagine the bottle to be connected at all, but of course my parameterization returns identical values at u=0 and u=2*pi. Removing all the end-connection code and extending the evaluation along u all the way to 2*pi created a visually AND geometrically correct part. The analyzer had to fix 2 points where rounding errors outputted slightly different numbers. (View the analyzer's output).
To prevent gridline intersections where the shape self-intersects, I used my interactive geometry generator and the Geomview viewer. Usually, I adjusted Geomview's near and far clipping planes to enclose a very small slice of my shape, so I could see the critical area without any foreground or background geometry in the way. Then I orbited the area over and over so that I could see "sky" between each pair of gridlines that might be intersecting. I never checked the region analytically, but I think I was able to see a gap everywhere there should be one.