Visualize singular vector fields with TikZ or PsTricks and friends

I am trying to find a good way to make really visually appealing plots of singular vector fields. For example electric fields of point like charges. I already tried this in python without too much succes:

from pylab import *
from scipy.integrate import odeint
from matplotlib import animation
from matplotlib import cm
import numpy as np

rc('font', **{'family': 'serif', 'serif': ['Computer Modern']})
rc('text', usetex=True)

## Set up charges

class charge:
def __init__(self, q, pos):
self.q=q
self.pos=pos

chargesPlus=[]
chargesMinus=[]

#for i in arange(0,1,1):
chargesPlus.append(charge(1,[3,0]))
chargesMinus.append(charge(-1,[-3,0]))
charges = chargesPlus + chargesMinus

def E_point_charge(q, a, x, y,r):
return q*(x-a)/((x-a)**2+(y-a)**2)**(1.5), \
q*(y-a)/((x-a)**2+(y-a)**2)**(1.5)

def E_total(x, y, charges):
Ex, Ey=0, 0
for C in charges:
E=E_point_charge(C.q, C.pos, x, y,1)
Ex=Ex+E
Ey=Ey+E
return [Ex, Ey]

domain =2

## Cut Quiver plot
def cut(r):
if r < domain:
return 0
else:
return 1

cutv = np.vectorize(cut)

def cut_total(charges,x):
c = 1
for C in charges:
r = sqrt((C.pos - x)**2 + (C.pos - x)**2)
c = c*cutv(r)
print c
print C.pos,C.pos
return c

fig = figure()

xMin,xMax=-15,15
yMin,yMax=-10,10

#ax.plot(x,y)
ax.axis('tight')
xlim([xMin,xMax])
ylim([yMin,yMax])

# plot point charges
for C in charges:
if C.q>0:
plot(C.pos, C.pos, 'bo', ms=8*sqrt(C.q))
if C.q<0:
plot(C.pos, C.pos, 'ro', ms=8*sqrt(-C.q))

xG,yG = meshgrid(linspace(xMin,xMax,25),linspace(yMin,yMax,25))

# plot vector field
E_totalX,E_totalY = E_total(xG,yG,charges)

EAbs = (E_totalX**2 + E_totalY**2)**(0.5)
E_XX = E_totalX/EAbs
E_YY = E_totalY/EAbs
#EAbs = np.nan_to_num(EAbs)

#ax.streamplot(xG,yG,E_XX,E_YY,color=EAbs,cmap=cm.autumn)
ax.quiver(xG,yG,E_XX,E_YY,EAbs,cmap=cm.GnBu)

xlabel('x')
ylabel('y')
ax.set_aspect(1)

plt.savefig('fig1.png')

E_totalX = E_totalX*cut_total(charges,[xG,yG])
E_totalY = E_totalY*cut_total(charges,[xG,yG])

ax.cla()

# plot point charges
for C in charges:
if C.q>0:
plot(C.pos, C.pos, 'bo', ms=8*sqrt(C.q))
if C.q<0:
plot(C.pos, C.pos, 'ro', ms=8*sqrt(-C.q))

xlabel('x')
ylabel('y')
ax.set_aspect(1)

plt.savefig('fig2.png')
#show()  The problem with the first example is how to choose colormap. I played with different maps from http://matplotlib.org/examples/color/colormaps_reference.html but none gave me a really satisfying result.

The problem with the first example is the cut function which doesn't work automatically. You need to cut out a specific region manually.

I know that there are powerful LaTeX packages for drawing and plotting and I am curious if this problem could be solved better, more automatically and more visual appealing using some LaTeX drawing or plotting techniques as my python approach above. Perhaps it would be also a good idea to use sagetex or pythontex to make the calculations and use the TeX packages for plotting. So, how would a TikZ/PSTricks ... guru solve this problem?

Please convert this question to community wiki if possible.

Run with xelatex(takes some time) or latex->dvips->ps2pdf:

\documentclass{article}
\usepackage{pst-electricfield}
\begin{document}

\begin{pspicture*}(-6,-6)(6,6)
\psframe*[linecolor=lightgray!50](-6,-6)(6,6)
\psgrid[subgriddiv=0,gridcolor=gray,griddots=10]
\psElectricfield[Q={[-1 -2 2][1 2 2][-1 2 -2][1 -2 -2]},linecolor=red]
\psEquipotential[Q={[-1 -2 2][1 2 2][-1 2 -2][1 -2 -2]},
linecolor=blue](-6.1,-6.1)(6.1,6.1)
\psEquipotential[Q={[-1 -2 2][1 2 2][-1 2 -2][1 -2 -2]},
linecolor=green,linewidth=2\pslinewidth,Vmax=0,Vmin=0](-6.1,-6.1)(6.1,6.1)
\end{pspicture*}

\end{document} and for a stream density plot:

\documentclass[pstricks,border=5pt]{standalone}
\usepackage{pst-magneticfield}
\begin{document}

\begin{pspicture}(-6,-4)(6,4)
\psmagneticfield[N=3,R=2,L=2,StreamDensityPlot](-6,-4)(6,4)
\end{pspicture}

\end{document} • Nice. What's missing is some type of color map for the field lines and (or) the equipotential lines to indicate the field strength. Is it also possible to get the vector field plot and not the field lines? – Julia Jan 28 '15 at 8:47
• see edited answer – user2478 Jan 28 '15 at 9:25
• Hi. I have tried to run the first code but there is an error: ! Emergency stop. <*> linee_campo_el_4_cariche.tex *** (job aborted, no legal \end found) Here is how much of TeX's memory you used: 7441 strings out of 493727 135849 string characters out of 6147890 186550 words of memory out of 5000000 10729 multiletter control sequences out of 15000+600000 3640 words of font info for 14 fonts, out of 8000000 for 9000 1328 hyphenation exceptions out of 8191 43i,1n,51p,485b,203s stack positions out of 5000i,500n,10000p,200000b,80000s No pages of output. – ryuk May 4 '16 at 5:12