# 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[0])/((x-a[0])**2+(y-a[1])**2)**(1.5), \
q*(y-a[1])/((x-a[0])**2+(y-a[1])**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[0]
Ey=Ey+E[1]
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[0] - x[0])**2 + (C.pos[1] - x[1])**2)
c = c*cutv(r)
print c
print C.pos[0],C.pos[1]
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[0], C.pos[1], 'bo', ms=8*sqrt(C.q))
if C.q<0:
plot(C.pos[0], C.pos[1], '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[0], C.pos[1], 'bo', ms=8*sqrt(C.q))
if C.q<0:
plot(C.pos[0], C.pos[1], '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? Commented Jan 28, 2015 at 8:47