This is a known issue with mplot3d. The issue is that the polygons are essentially abstracted 2-D objects that are using the same backend architecture as the core 2-D plotting software. I think your example script is probably the *best* example of the problem that I have seen and does a very good job of illustrating the problem.
The crux of the problem is that each polygon is assigned a z-order for layering, and this z-order is determined by (i think) the center-of-mass location of the polygon and how it relates the location of the centers of mass of the other objects. In most graphing situations, this is a good enough hack, but there are too many use-cases where it falls apart.
Because we are limited to a single z-order value for each object, and compositing is done through a 2-D rendering engine, this problem will likely not get solved any time soon, unfortunately. The fix would probably require a complete rewrite of the mpl drawing engine or maybe the utilization of OpenGL?
Sorry I can't be more helpful in fixing your problem, it has been a aggravating issue for me as well. Thank you for your excellent example script.