Intersection of 2 surface plots in pgfplots

My question is pretty much a duplicate of this one. I'd like pgfplots to display the intersection between tho surface plots properly. Unfortunately, my function is quite different from the one in the mentioned question, so I don't know, where the surfaces will meet. Would there be any kind of workaround for this kind of function?

\documentclass{scrartcl}
\usepackage{pgfplots}
\begin{document}
\begin{tikzpicture}
\begin{axis}[domain=0.01:30]
\end{axis}
\end{tikzpicture}
\end{document}


• Unfortunately, I think this is a job for PSTricks... – Jake Sep 24 '12 at 19:05
• As Jake already said, this is beyond the capabilities of pgfplots. The golatex forum has a related post in which the OP managed to generate the graphics as such with matlab and imported it via \addplot3 graphics: golatex.de/… (in german only - but pictures speak in their own language, I guess) – Christian Feuersänger Sep 25 '12 at 18:40
• I just discovered this website tlhiv.org/mpgraph that render 3D surfaces in MetaPost (and other formats, such as PDF), with surface intersections. Unfortunately there is no TikZ output. – alfC Oct 3 '17 at 8:05

Another hackish solution:

\documentclass{scrartcl}
\usepackage{pgfplots}
\begin{document}
\begin{tikzpicture}
\begin{axis}[domain=0.01:30]
\end{axis}
\end{tikzpicture}
\end{document}


Note that this kind of solution is less flexible, because the correct hidden surface removal depends on the position of the camera, and also on special properties of shape of the functions. If one could know the point of view internally, one could generalize the max and min function to make it camera dependent and in this way simulate hidden surfaces.

• Thanks a lot for both your solutions. Would you be so kind as to explain the rationale behind the value 118.89 and especially behind the second solution? I won't need different camera angles, but I will need to plot a lot of similar functions. – meep.meep Sep 25 '12 at 8:20
• Second solution meaning this one, the order changed because of my accepting this answer. – meep.meep Sep 25 '12 at 9:20
• I got the solution by myself. I just was missing the rule: ln(e^{f(x)}) = f(x). Using that, I now understand where the 118.89 come from. Thanks again. – meep.meep Sep 25 '12 at 14:06
• Missing a parentheses ) at the end of the first \addplot3 – juanuni Sep 5 '16 at 1:58
• I have found another hack that allows to plot an arbitrary number of surfaces correctly, without doing anything manually: tex.stackexchange.com/a/394066/38641. Moreover, it does not depend on the viewpoint. – iavr Oct 1 '17 at 10:58

A combination of surface colors, opacities and parametric plots can get you close to the desired result:

Code follows:

\documentclass{scrartcl}
\usepackage{pgfplots}
\begin{document}
\begin{tikzpicture}
\begin{axis}[domain=0.01:30]
\end{axis}
\end{tikzpicture}
\end{document}

• Could you please explain how this works? I am trying to modify it for the x=z plane and a different function in place of the exponential but failing. In fact as soon as I replace {0} with {x} in the first addplot3 I just get two planes. – Arrow Jan 17 '18 at 16:54
• @Arrow what do you expect? You have a problem of scale. The surface looks a plane because of the scale. – alfC Jan 18 '18 at 2:01

Yet another hackish solution, based on this answer. I am adding it here in response to the comment by alfC:

\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; import scipy.linalg as la;%
#1 np.savetxt(sys.stdout, data)%
"}

\newcommand<>{\pyplot}[3][]%

\newcommand<>{\pyplott}[3][]%

\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,30,n); y = x;
X, Y = np.meshgrid(x, y);
Z1 = np.zeros([n, n]);
Z2 = (1-0.3)*np.exp(-X*(Y/100)*(1-0.3))-np.exp(-X*(Y/100));
M1 = np.ones([n, n]);
M2 = 2 * M1;
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]);
M = np.r_[M1, N, M2].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);},
]
\end{axis}
\end{tikzpicture}
\end{center}

\end{frame}

\end{document}


Here is the result:

So, we only need to provide one expression for each surface, and nothing manual for the intersection(s). And it works for any viewpoint. Such automated solution is needed when surfaces are very complex.

However, the intersection looks ugly because each patch is either drawn or not - while two intersecting patches themselves should be partially visible. I guess a single self-intersecting surface would have the same problem.

• Not that bad. you can cheat and add an intersection line to hide the rough intersection. With lualatex maybe you can have a more sense grid? – alfC Oct 2 '17 at 13:23
• Maybe, but for now I am experimenting with matplotlib PGF backend instead: matplotlib.org/users/pgf.html – iavr Oct 2 '17 at 13:34
• If we are going to the external tools realm, many possibilities open up! tex.stackexchange.com/a/394337/1871 – alfC Oct 3 '17 at 7:48

Use asymptote, this is a more appropriate tool for plotting in 3D. With some additional work you can make it look like a pgfplots graph. Below there are three ways of doing this with asymptote.

Transparent surfaces

settings.outformat="png";
settings.render=16;

import graph3;
import palette;

size3(150,IgnoreAspect);

currentprojection=orthographic(50,-60,0.4);
currentlight=(5,-3,1);

real f(pair z){return (1.-0.3)*exp(-z.x*(z.y/100.)*(1.-0.3))-exp(-z.x*(z.y/100.));}
real cero(pair z){return 0.;}

limits((0,0,-0.3),(30,30,0.1));
xaxis3(Label("$x$",1),blue,arrow=Arrow3);
yaxis3(Label("$y$",1),blue,arrow=Arrow3);
zaxis3(Label("$z$",1),blue,arrow=Arrow3);

surface s=surface(f,(0,0),(30,30),30,30,Spline);
surface s2=surface(cero,(0,0),(30,30),30,30,Spline);

s.colors(palette(s.map(zpart),Rainbow()+opacity(0.5)));
s2.colors(palette(s2.map(zpart),Rainbow()+opacity(0.5)));

draw(s,meshpen=blue);
draw(s2,meshpen=red);

shipout(bbox(2mm,Fill(white)));


Solid surfaces

settings.outformat="png";
settings.render=16;

import graph3;
import palette;

size3(150,IgnoreAspect);

currentprojection=orthographic(50,-60,0.4);
//currentlight=(-5,3,20);

real f(pair z){return (1.-0.3)*exp(-z.x*(z.y/100.)*(1.-0.3))-exp(-z.x*(z.y/100.));}
real cero(pair z){return 0.;}

limits((0,0,-0.3),(30,30,0.1));
xaxis3(Label("$x$",1),blue,arrow=Arrow3);
yaxis3(Label("$y$",1),blue,arrow=Arrow3);
zaxis3(Label("$z$",1),blue,arrow=Arrow3);

surface s=surface(f,(0,0),(30,30),30,30,Spline);
surface s2=surface(cero,(0,0),(30,30),30,30,Spline);

//s.colors(palette(s.map(zpart),Rainbow+opacity(0.5)));
//s2.colors(palette(s2.map(zpart),Rainbow()+opacity(0.5)));

draw(s, blue,meshpen=blue);
draw(s2, red, meshpen=red);

shipout(bbox(2mm,Fill(white)));


PGFplots style (still with Asymptote)

This is the closest I was able to make it look like a pgfplots plot.

(Things to do, help is welcomed: 1) automatic point of view (to not depend on the scale of the plot, 2) thick border around each surface 3) vector graphics and hidden surface (e.g. via bsd module) 4) details, like small tics and tics spacing)

settings.outformat="png";
settings.render=8;
import grid3;
import graph3;
import palette;

size3(160,125,90,IgnoreAspect);
currentprojection=orthographic(30,-60,0.8);

real f(pair z){return (1.-0.3)*exp(-z.x*(z.y/100.)*(1.-0.3))-exp(-z.x*(z.y/100.));}
real cero(pair z){return 0.;}

limits((0,0,-0.3),(30,30,0.1));

xaxis3(Bounds(Min,Min), InTicks());
xaxis3(Bounds(Max,Max));
xaxis3(Bounds(Max,Min));
yaxis3(Bounds(Max,Min), InTicks(beginlabel=false) );
yaxis3(Bounds(Min,Max));
yaxis3(Bounds(Min,Min));
zaxis3(Bounds(Min,Min), InTicks() );
zaxis3(Bounds(Max,Max));
zaxis3(Bounds(Min,Max));

surface s=surface(f,(0,0),(30,30),24,24);
surface s2=surface(cero,(0,0),(30,30),24,24);

draw(s2, red, meshpen=0.2pt + black, nolight);
draw(s2, red, meshpen=0.2pt + black, nolight);

shipout(bbox(2mm,Fill(white)));


Notes

asymptote can output vector graphics, but the output is not as good as a raster (rendered) output and removal of hidden faces is very difficult (like with pgfplots.) So, this is big drawback of using asymptote, the plot becomes raster. However, it is common that 3D plots quickly become too complex to render with vector graphics.

• Interesting, but I can see the same "zig-zag" issue with patches at the intersection, as in my pgfplots solution. It will be clearer if you choose a fixed different color for each surface. – iavr Oct 3 '17 at 13:21
• As for other external tools, I tried Python matplotlib/mplot3d, but it doesn't do surface intersection at all, which is amazing. Gnuplot does it nicely, but I cannot see how to make it use latex for text without graphics being rasterized. So, for now, I let Matlab win one more battle. – iavr Oct 3 '17 at 13:29
• @iavr, the zig-zag is an issue of transparency (which was not part of the original question). Solid intersections are extremely excellent quality in Asymptote. See my edit. (Of course we are in raster-land, so all is good down to the pixel level). – alfC Oct 3 '17 at 21:15
• Now that's better!! – iavr Oct 3 '17 at 22:51

This question has already tons of excellent answers, but as far as I can see none of them exploits that (1-0.3)*e^(-x*(y/100)*(1-0.3))-e^(-x*(y/100))=0 implies that y=-1000*ln(0.7)/(3*x)=:pft/x (pft is marmot language and means something like "convenient constant", but is hard to translate precisely;-), and that, in order to draw a plane, one really does not need to use \addplot3. All I am doing here is to draw the part x*y>pft first, than the curved surface, and then the x*y<pft part of the plane.

\documentclass[border=3.14mm]{standalone}
\usepackage{pgfplots}
\pgfplotsset{compat=1.16}
\begin{document}
\begin{tikzpicture}
\begin{axis}[domain=0.01:30,xlabel=$x$]
\pgfmathsetmacro{\pft}{-1000*ln(0.7)/3}
\fill[blue] plot[variable=\x,domain={\pft/30}:30] ({\x},{\pft/\x},0) --
(30,30,0) -- cycle;