[Tex/LaTex] How to draw force lines onto an equipotential contour plot using TikZ

Question

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
    ]
\addplot+[
    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.

Answer 1

(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])

def addplot(x, y):
    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
    ]
\addplot+[
    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):
          addplot(-2+0.5*np.cos(2*np.pi*i/16), 0.5*np.sin(2*np.pi*i/16))}
  \pyc{for i in range(16):
          addplot(2+0.3*np.cos(2*np.pi*i/16), 0.3*np.sin(2*np.pi*i/16))}
\end{scope}
\end{axis}
\end{tikzpicture}
\end{document}

And the output is: