% Author: Izaak Neutelings (July 2018)
% page 8 https://archive.org/details/StaticAndDynamicElectricity
% https://stackoverflow.com/questions/4919827/gnuplot-with-texshop-in-osx --shell-escape
% https://tex.stackexchange.com/questions/56353/extract-x-y-coordinate-of-an-arbitrary-point-on-curve-in-tikz
% https://tex.stackexchange.com/questions/412899/tikz-calculate-and-store-the-euclidian-distance-between-two-coordinates

\documentclass[border=3pt,tikz]{standalone}
\usepackage{amsmath} % for \dfrac
\usepackage{bm}
\usepackage{tikz,pgfplots}
\tikzset{>=latex} % for LaTeX arrow head
\pgfplotsset{compat=1.13}
\usetikzlibrary{decorations.markings,intersections,calc}
\usepackage{ifthen}
%\usepackage{etoolbox}
\usepackage{xcolor}
\colorlet{Ecol}{orange!90!black}
\colorlet{EcolFL}{orange!80!black}
\colorlet{Bcol}{blue!90!black}
\tikzstyle{EcolEP}=[blue!80!white]
\tikzstyle{charge+}=[very thin,top color=red!50,bottom color=red!90!black,shading angle=20]
\tikzstyle{charge-}=[very thin,top color=blue!50,bottom color=blue!80,shading angle=20]
%\tikzstyle{EFielLineArrow2}=[
%    EcolFL,decoration={markings,
%          mark=at position 0.5 with {\arrow[rotate=\angle]{latex}}},
%          postaction={decorate}]
\tikzset{
   EFielLineArrow/.style args = {#1}{EcolFL,decoration={markings,
          mark=at position 0.5 with {\arrow[rotate=#1]{latex}}},
          postaction={decorate}}
}

\makeatletter
  \newcommand{\xy}[3]{% % FIND X, Y
    \tikz@scan@one@point\pgfutil@firstofone#1\relax
    \edef#2{\the\pgf@x}%
    \edef#3{\the\pgf@y}%
  }
\makeatother
\newcommand{\EFielLineArrow}[2]{ % ELECTRIC FIELD LINE ARROW
  \pgfkeys{/pgf/fpu,/pgf/fpu/output format=fixed} % for calculation between -1*10^324 and +1*10^324
  \pgfmathsetmacro{\x}{#1/28.45pt}
  \pgfmathsetmacro{\y}{#2/28.45pt}
  \pgfmathsetmacro{\U}{\Q*((\x+\a)^2+(\y)^2)^(3/2)}
  \pgfmathsetmacro{\V}{\q*((\x-\a)^2+(\y)^2)^(3/2)}
  \pgfkeys{/pgf/fpu=false}
  \pgfmathparse{
    atan2(((\y)*\V + (\y)*\U),((\x+\a)*\V + (\x-\a)*\U))
  }
  \edef\angle{\pgfmathresult}
  \pgfmathsetmacro{\D}{int(1000*\q*(\x+\a)/sqrt((\x+\a)^2+\y*\y) + 1000*\Q*(\x-\a)/sqrt((\x-\a)^2+\y*\y))/1000}
  \draw[EFielLineArrow={\angle}] (\xpt,\ypt);
}
\newcommand{\EVector}[2]{ % ELECTRIC FIELD VECTOR
  \pgfmathsetmacro{\x}{#1/28.45pt}
  \pgfmathsetmacro{\y}{#2/28.45pt}
  \pgfmathsetmacro{\U}{((\x+\a)^2+(\y)^2)^(3/2)}
  \pgfmathsetmacro{\V}{((\x-\a)^2+(\y)^2)^(3/2)}
  \pgfmathsetmacro{\Ex}{\q*3*(\x+\a)/\U + \Q*3*(\x-\a)/\V}
  \pgfmathsetmacro{\Ey}{\q*3*(\y)/\U + \Q*3*(\y)/\V}
  \fill[Ecol] (\x,\y) circle (0.8pt);
  \draw[->,Ecol,thick] (\x,\y) --++ (\Ex,\Ey)
    node[Ecol,right,scale=0.8] {$\bm{E}$};
  \pgfmathsetmacro{\D}{int(1000*\q*(\x+\a)/sqrt((\x+\a)^2+\y*\y) + 1000*\Q*(\x-\a)/sqrt((\x-\a)^2+\y*\y))/1000}
}
\newcommand{\EFieldLines}{ % ELECTRIC FIELD LINES
  \message{^^JField lines (\q,\Q) with contours range ^^J\range^^J}
  
  % PATHS for intersections
  \path[name path=ellipse1] (-\a,0) ellipse ({0.90*\R} and {1.5*\R});
  \path[name path=ellipse2] (+\a,0) ellipse ({0.75*\R} and {1.5*\R});
  \path[name path=ellipse3] (  0,0) ellipse ({\a+\R} and 1.5*\R);
  
  % FIELD LINES
  \draw[EcolFL,name path=Elines] plot[id=plot, raw gnuplot, smooth]
    function{
       f(x,y) = \q*(x+\a)/sqrt((x+\a)**2+y**2) + \Q*(x-\a)/sqrt((x-\a)**2+y**2);
       set xrange [\xmin:\xmax];
       set yrange [-\ymax:\ymax];
       set view 0,0;
       set isosample 400,400;
       set cont base;
       set cntrparam levels discrete \range;
       unset surface;
       splot f(x,y)
    };
  
  % ELLIPSE INTERSECTIONS
  \pgfmathsetmacro{\oppositesign}{\q*\Q<0 ? 1 : 0}
  \ifthenelse{\oppositesign>0}{
    % OPPOSITE SIGN
    \foreach \c in {1,2}{
      \message{Intersections \c...}
      \path[name intersections={of=Elines and ellipse\c,total=\t}]
        \pgfextra{\xdef\Nb{\t}};
      \message{ found \Nb ^^J}
      \foreach \i in {1,...,\Nb}{
        \message{  \i}
        \xy{(intersection-\i)}{\xpt}{\ypt}
        \EFielLineArrow{\xpt}{\ypt}
        \message{ (\D,\x,\y)^^J}
      }
    }
  }{
    % SAME SIGN
    \message{Intersections...}
    \path[name intersections={of=Elines and ellipse3,total=\t}]
        \pgfextra{\xdef\Nb{\t}}; 
    \message{ found \Nb ^^J}
    \foreach \i in {1,...,\Nb}{
      \message{  \i}
      \xy{(intersection-\i)}{\xpt}{\ypt}
      \EFielLineArrow{\xpt}{\ypt}
      \message{ (\D,\x,\y)^^J}
    }
  }
}
\newcommand{\EEquipot}{ % EQUIPOTENTIAL SURFACE
  \message{^^JEquipotential surface (\q,\Q) with contours range ^^J\rangeEP}
  
  % FIELD LINES
  \draw[EcolEP] plot[id=plot, raw gnuplot, smooth] %,dashed
    function{
       f(x,y) = \q/sqrt((x+\a)**2+y**2) + \Q/sqrt((x-\a)**2+y**2);
       set xrange [\xmin:\xmax];
       set yrange [-\ymax:\ymax];
       set view 0,0;
       set isosample 400,400;
       set cont base;
       set cntrparam levels discrete \rangeEP;
       unset surface;
       splot f(x,y)
    };
}
\newcommand{\EVectorOnFieldLines}{ % ELECTRIC FIELD VECTOR
  \message{^^JVector on field lines (\q,\Q) with contours range ^^J\C^^J}
  
  % FIELD LINES
  \path[name path=xline] (\Cx,-\ymax) -- (\Cx,-\a) (\Cx,\a) -- (\Cx,\ymax);
  \path[name path=Elines] plot[id=plot, raw gnuplot, smooth]
    function{
       f(x,y) = \q*(x+\a)/sqrt((x+\a)**2+y**2) + \Q*(x-\a)/sqrt((x-\a)**2+y**2);
       set xrange [\xmin:\xmax];
       set yrange [-\ymax:\ymax];
       set view 0,0;
       set isosample 100,100;
       set cont base;
       set cntrparam levels discrete \C;
       unset surface;
       splot f(x,y)
    };
  
    \message{Intersections...}
    \path[name intersections={of=Elines and xline,total=\t}]
      \pgfextra{\xdef\Nb{\t}};
    \message{ found \Nb ^^J}
    \foreach \i in {1,...,\Nb}{
      \message{  \i}
      \xy{(intersection-\i)}{\xpt}{\ypt}
      \EVector{\xpt}{\ypt}
      \message{ (\D,\x,\y)^^J}
    }
}



\begin{document}


% OPPOSITE, +1 -1
\begin{tikzpicture}
  \def\xmin{-3}
  \def\xmax{3}
  \def\ymax{3}
  \def\a{1}
  \def\q{+1}
  \def\Q{-1}
  \def\R{1.0}
  \def\range{0.05,0.2,0.6,1.0,1.4,1.8,1.98}

  % LINES
  \EFieldLines
  
  % VECTORS
  \def\Cx{-0.4}
  \def\C{0.6,1.0}
  \EVectorOnFieldLines
  \def\Cx{-1.26}
  \EVectorOnFieldLines
  \def\Cx{1.26}
  \EVectorOnFieldLines
  
  % CHARGES
  \draw[charge+] (-\a,0) circle (9pt) node[black,scale=1.0] {$+q$};
  \draw[charge-] (+\a,0) circle (9pt) node[black,scale=1.0] {$-q$};
   
\end{tikzpicture}


% OPPOSITE, +2 -1
\begin{tikzpicture}
  \def\xmin{-3}
  \def\xmax{3}
  \def\ymax{3}
  \def\a{1}
  \def\q{+2}
  \def\Q{-1}
  \def\R{1.0}
  %\def\range{-2.9,-2.8,-2.6,-2.2,-1.8,-1.4,
  %           -1.01,-0.6,-0.2,0.2,0.6,0.8,1.0,1.2}
  \def\range{-0.95,-0.8,-0.6,-0.2,0.2,0.6,
             1.01,1.4,1.8,2.2,2.6,2.8,2.94}

  % LINES
  \EFieldLines
  
  % CHARGES
  \draw[charge+] (-\a,0) circle (12pt) node[black,scale=1.0] {$+2q$};
  \draw[charge-] (+\a,0) circle ( 9pt) node[black,scale=1.0] {$-q$};
  
\end{tikzpicture}


% OPPOSITE, +1 +1
\begin{tikzpicture}
  \def\xmin{-3}
  \def\xmax{3}
  \def\ymax{3}
  \def\a{1}
  \def\q{+1}
  \def\Q{+1}
  \def\R{1.2}
  \def\range{-1.99,-1.8,-1.4,-1.0,-0.6,-0.2,
             0.2,0.6,1.0,1.4,1.8,1.99}
  
  % LINES
  \EFieldLines
  
  % VECTORS
  \def\Cx{-1.42}
  \def\C{-1.0}
  \EVectorOnFieldLines
  \def\Cx{-0.85}
  \def\C{-0.6}
  \EVectorOnFieldLines
  \fill[Ecol] (0,0) circle (1pt) node[Ecol,above,scale=0.7] {$\mathbf{E}=0$};
  
  % CHARGES
  \draw[charge+] (-\a,0) circle (9pt) node[black,scale=1.0] {$+q$};
  \draw[charge+] (+\a,0) circle (9pt) node[black,scale=1.0] {$+q$};
  
\end{tikzpicture}


% OPPOSITE, +2 +1 (equipotential surface)
\begin{tikzpicture}
  \def\xmin{-3}
  \def\xmax{3}
  \def\ymax{3}
  \def\a{1}
  \def\q{+2}
  \def\Q{+1}
  \def\R{1.2}
  \def\range{-2.98,-2.7,-2.1,-1.5,-0.9,-0.3,
              0.3,0.9,1.5,2.1,2.7,2.98}

  % LINES
  \EFieldLines
  \fill[Ecol] ({(3-2*sqrt(2))*\a},0) circle (1pt) node[Ecol,above,scale=0.6] {$\mathbf{E}=0$};
  
  % CHARGES
  \draw[charge+] (-\a,0) circle (12pt) node[black,scale=1.0] {$+2q$};
  \draw[charge+] (+\a,0) circle ( 9pt) node[black,scale=1.0] {$+q$};
   
\end{tikzpicture}


% OPPOSITE, +1 -1 (equipotential surface)
\begin{tikzpicture}
  \def\xmin{-3}
  \def\xmax{3}
  \def\ymax{3}
  \def\a{1}
  \def\q{+1}
  \def\Q{-1}
  \def\R{1.0}
  \def\range{0.05,0.2,0.6,1.0,1.4,1.8,1.98}
  \def\rangeEP{-1.8,-1.4,-1.0,-0.6,-0.2,0.0,0.2,0.6,1.0,1.4,1.8}
  
  % LINES
  \EFieldLines
  \EEquipot
  
  % CHARGES
  \draw[charge+] (-\a,0) circle (9pt) node[black,scale=1.0] {$+q$};
  \draw[charge-] (+\a,0) circle (9pt) node[black,scale=1.0] {$-q$};
   
\end{tikzpicture}


% OPPOSITE, +1 +1 (equipotential surface)
\begin{tikzpicture}
  \def\xmin{-3}
  \def\xmax{3}
  \def\ymax{3}
  \def\a{1}
  \def\q{+1}
  \def\Q{+1}
  \def\R{1.2}
  \def\range{-1.99,-1.8,-1.4,-1.0,-0.6,-0.2,
             0.2,0.6,1.0,1.4,1.8,1.99}
  \def\rangeEP{0.0,0.1,0.2,0.4,0.6,0.8,1.0,1.2,1.4,1.6,1.8,2.002,2.2,2.4,2.6}
  %-1.8,-1.4,-1.0,-0.6,-0.2,
  
  % LINES
  \EFieldLines
  \EEquipot
  
  % CHARGES
  \draw[charge+] (-\a,0) circle (9pt) node[black,scale=1.0] {$+q$};
  \draw[charge+] (+\a,0) circle (9pt) node[black,scale=1.0] {$+q$};
   
\end{tikzpicture}



\end{document}
