There are several other questions on this issue, for instance:
As pointed out by Matthew Leingang, this is due to a limitation of pgfplots.
Unlike all workarounds so far, I have found one (a hack, really) that allows to draw more than one, in fact an arbitrary number of, potentially intersecting, surfaces in a single addplot3 command, with automatic z buffering, without doing anything manually.
We use addplot3 table instead of addplot3 coordinates, and we generate data externally. A single surface needs three matrices with x, y, z cordinates. For two surfaces, we can stack together [x; x], [y; y] and [z1; z2]. To make the two surfaces disconnected, we can insert a vector, say n, of NaNs of appropriate size between the stacked matrices, e.g. [x; n; x], [y; n; y] and [z1; n; z2], together with option unbounded coords=jump. Finally, we save the three matrices as three stacked columns representing (x,y,z) triplets as pgfplots expects. This also requires specifying the number of columns of the matrices with mesh/cols.
To implement this idea, I define some macros that allow calling arbitrary python code via addplot shell, saving the data to a text file in tabular form, and then loading for display. This requires the -shell-escape flag in pdflatex.
Unfortunately, because this is a single plot, I cannot see how to specify different properties (e.g. color or opacity) for each individual surfaces. Well, maybe by adding a fourth column in the data combined with point meta option as in scatter plots, but I haven't tried that.
Also, by trying more complex examples, one realizes that, although patch visibility is computed correctly, we don't really get patch intersection. So, to get the feeling of a smooth curve at the surface intersection, one needs to increase the resolution. I do not intend to use this; I am just sharing because I found it interesting.
Below, I am giving an example of two intersecting planes, but really one could compute anything with the same 'method'. It is in beamer, because this is what I was trying already.
\documentclass{beamer}
\usefonttheme[onlymath]{serif}
\setbeamersize{text margin left=10pt}
\setbeamersize{text margin right=10pt}
\usepackage{pgfplots}
\pgfplotsset{
every axis/.append style={font=\scriptsize},
plain/.style={every axis plot/.append style={mark=none},enlargelimits=false,grid=none},
z-sort/.style={z buffer=sort,unbounded coords=jump},
}
\newcommand{\python}[1]{python -c "%
import math, sys; import numpy as np;%
#1 np.savetxt(sys.stdout, data)%
"}
\newcommand<>{\pyplot}[3][]%
{\only#4{\addplot[#1] shell[prefix=fig/data/,id=#2,] {\python{#3}};}}
\newcommand<>{\pyplott}[3][]%
{\only#4{\addplot3[z-sort,#1] shell[prefix=fig/data/,id=#2,] {\python{#3}};}}
\newcommand<>{\pyload}[3][]%
{\only#4{\addplot[#1] table[x index=0,y index=#2] {fig/data/#3.out};}}
\newcommand<>{\pyloadt}[2][]%
{\only#3{\addplot3[z-sort,#1] table {fig/data/#2.out};}}
\newcommand{\pysave}[2]{
\begin{tikzpicture}[overlay,opacity=0]
\begin{axis} \pyplot{#1}{#2} \end{axis}
\end{tikzpicture}
}
\begin{document}
\begin{frame}
\pysave{surf}{
n = 11; x = np.linspace(0,1,n); y = x;
X, Y = np.meshgrid(x,y);
Z1 = X + Y;
Z2 = 1 - X + Y;
N = np.ones([1, n]) * np.NaN;
X = np.r_[X, N, X ].reshape([-1, 1]);
Y = np.r_[Y, N, Y ].reshape([-1, 1]);
Z = np.r_[Z1, N, Z2].reshape([-1, 1]);
data = np.c_[X, Y, Z];
}
\begin{center}
\begin{tikzpicture}
\begin{axis}[plain,width=\textwidth,height=.8\textwidth]
\pyloadt[surf,mesh/cols=11]{surf};
\end{axis}
\end{tikzpicture}
\end{center}
\end{frame}
\end{document}
The result looks like this:
EDIT
It is possible, eventually, to color each surface differently. I couldn't make point meta=explicit symbolic or point meta=explicit work, but what did work is point meta=\thisrowno{3}. Here is the code:
\documentclass{beamer}
\usefonttheme[onlymath]{serif}
\setbeamersize{text margin left=10pt}
\setbeamersize{text margin right=10pt}
\usepackage{pgfplots}
\pgfplotsset{
every axis/.append style={font=\scriptsize},
plain/.style={every axis plot/.append style={mark=none},enlargelimits=false,grid=none},
z-sort/.style={z buffer=sort,unbounded coords=jump},
}
\newcommand{\python}[1]{python -c "%
import math, sys; import numpy as np;%
#1 np.savetxt(sys.stdout, data)%
"}
\newcommand<>{\pyplot}[3][]%
{\only#4{\addplot[#1] shell[prefix=fig/data/,id=#2,] {\python{#3}};}}
\newcommand<>{\pyplott}[3][]%
{\only#4{\addplot3[z-sort,#1] shell[prefix=fig/data/,id=#2,] {\python{#3}};}}
\newcommand<>{\pyload}[3][]%
{\only#4{\addplot[#1] table[x index=0,y index=#2] {fig/data/#3.out};}}
\newcommand<>{\pyloadt}[2][]%
{\only#3{\addplot3[z-sort,#1] table {fig/data/#2.out};}}
\newcommand{\pysave}[2]{
\begin{tikzpicture}[overlay,opacity=0]
\begin{axis} \pyplot{#1}{#2} \end{axis}
\end{tikzpicture}
}
\begin{document}
\begin{frame}
\pysave{surf}{
n = 31; x = np.linspace(0,1,n); y = x;
X, Y = np.meshgrid(x,y);
Z1 = X + Y;
Z2 = 1 - X + Y;
Z3 = 1- X + 1 - Y;
M1 = np.ones([n, n]);
M2 = 2 * M1;
M3 = 3 * M1;
N = np.ones([1, n]) * np.NaN;
X = np.r_[X, N, X, N, X ].reshape([-1, 1]);
Y = np.r_[Y, N, Y, N, Y ].reshape([-1, 1]);
Z = np.r_[Z1, N, Z2, N, Z3].reshape([-1, 1]);
M = np.r_[M1, N, M2, N, M3].reshape([-1, 1]);
data = np.c_[X, Y, Z, M];
}
\begin{center}
\begin{tikzpicture}
\begin{axis}[
plain,width=\textwidth,height=.8\textwidth,
colormap={summap}{color=(green);color=(red);color=(yellow);},
]
\pyloadt[surf,opacity=.7,mesh/cols=31,point meta=\thisrowno{3}]{surf};
\end{axis}
\end{tikzpicture}
\end{center}
\end{frame}
\end{document}
In this example I am showing three planes colored in green, red, yellow. A little transparency helps seeing what is going on. It would be very complex in this case to compute intersections manually with min and max as in previous workarounds.
However, the missing patch intersections are now evident between the green and yellow planes, so I increased the resolution to 31x31. Further increasing to 41x41 gives TeX capacity exceeded, which is very sad. Anyhow, here is the result:
mlf2pdf.m.