# How to draw force lines onto an equipotential contour plot using TikZ?

I'm really new to TikZ and need to reproduce the following figure: which corresponds to the gravitational field produced by two dissimilar masses (being the mass on the left twice the mass on the right).

There's a special need to show the contour plot of equipotential surfaces (in light blue) and field lines (in red).

I've already managed to generate the contour plot using gnuplot but don't get a clue regarding how to superimpose the field lines onto it.

My TikZ code for the contour plot is as follows:

\documentclass{article}
\usepackage{tikz,pgfplots}

\begin{document}
\thispagestyle{empty}
\begin{tikzpicture}
\begin{axis}[
view={0}{90},
clip = false,
xmin = -4,
xmax = 4,
ymin = -2.3,
ymax = 2.3,
point meta min = -4,
point meta max = 3,
y axis line style={draw opacity=0},
x axis line style={draw opacity=0},
xtick=\empty,
ytick=\empty,
unit vector ratio=1 1 1,
width=\linewidth,
colormap name=viridis
]
no markers,
raw gnuplot,
contour prepared,
contour/labels=false
] gnuplot {set samples 50, 50;
set isosamples 55, 55;
set contour base;
set cntrparam levels incremental -4,0.22,4;
set style data lines;
splot [-4:4] [-2.3:2.3] (-2/sqrt((x+2)**2+y**2)-1/sqrt((x-2)**2+y**2));
};
\end{axis}
\end{tikzpicture}
\end{document}


The result of this code is: Any help will be very welcome.

• Not a real answer, just an idea: since field lines are always tangent to the equipotential lines, they have to satisfy the equation dr/ds = F(r(s)), where r is the current vector and F is your vector field (s is a curve coordinate) . You can solve this differential equation (which is not a real TeX problem) and then do a loop of initial conditions (i.e. vectors pointing in different directions) and plot it as parametric plot. – jaytar Sep 24 '17 at 7:11
• You can use a quiver plot but that will only give you straight arrows at every sampling point and not a smooth field line. pst-electricfield can do what you're looking for but uses PSTricks rather than TikZ :-/ – Henri Menke Sep 24 '17 at 7:29

(The answer is inspired by the jaytar's comment)

You can use pythontex for that. Note that if your potential function is U(x) then the field lines are orthogonal to its level curves, so they are solutions of the differential equation dx/dt=DU(x) (or -DU(x), it doesn't matter). In the following code I use the scpiy.integrate.odeint function to solve the equation numerically, and then just use the obtained coordinate pairs to plot the trajectory. Also, notice some clipping, otherwise it'd be a bit harder to find when the curve should hit the figure boundary.

To run the code below you'll have to install pythontex (and probably scipy). In Debian I just do

apt-get install texlive-extra-utils python-scipy


and then run

pdflatex --shell-escape file.tex
pythontex file.tex
pdflatex --shell-escape file.tex


The code is:

\documentclass{article}
\usepackage{tikz,pgfplots}
\pgfplotsset{compat=1.14}
\usepackage{pythontex}

\begin{pycode}
import numpy as np
from scipy import integrate

def func(Y, t):
x, y = Y
return [(2*(x+2))/((x+2)**2+y**2)**1.5 + (x-2)/((x-2)**2+y**2)**1.5, (2*y)/((x+2)**2+y**2)**1.5 + y/((x-2)**2+y**2)**1.5]

def parenthesize(x):
return '(%f,%f)' % (x[0], x[1])

t = np.arange(0, 7.0, 0.01)
sol = integrate.odeint(func, [x, y], t)
print('\\addplot[red] coordinates { ' + '\n'.join([parenthesize(x) for x in sol]) + ' };')
\end{pycode}

\begin{document}
\thispagestyle{empty}
\begin{tikzpicture}
\begin{axis}[
view={0}{90},
clip = false,
xmin = -4,
xmax = 4,
ymin = -2.3,
ymax = 2.3,
point meta min = -4,
point meta max = 3,
y axis line style={draw opacity=0},
x axis line style={draw opacity=0},
xtick=\empty,
ytick=\empty,
unit vector ratio=1 1 1,
width=\linewidth,
colormap name=viridis
]
no markers,
raw gnuplot,
contour prepared,
contour/labels=false
] gnuplot {set samples 50, 50;
set isosamples 55, 55;
set contour base;
set cntrparam levels incremental -4,0.22,4;
set style data lines;
splot [-4:4] [-2.3:2.3] (-2/sqrt((x+2)**2+y**2)-1/sqrt((x-2)**2+y**2));
};
\begin{scope}
\path[clip] (-4,-2.3)--(-4,2.3)--(4,2.3)--(4,-2.3)--cycle;
\pyc{for i in range(16):