-- plot vector field

-- by Gareth Rees  wgr2@cam.ac.uk
-- 11 March 2007
-- my first MathPad routine so not very elegant

-- Data to plot are in 4 arrays: xp, yp (x and y coordinates), xp2 and yp2 
-- (x- and y-components of field).
-- Vector arrows are centred at (xp, yp) and point in direction (xp2, yp2).
-- You will probably want to set plot limits Xmin, Xmax, Ymin, Ymax, and
-- also to set the Tracecolor (see example calling routine).

-- set up vector rescaling

-- This increases the range of vector lengths that can usefully be plotted
-- by applying a nonlinear (hyperbolic tangent) rescaling to them. If the
-- physical vector length is x, its plotted length is p = s1*s2*tanh(x/s1). 
-- This means that if x<<s1, p is approximately s2*x, while if x>>s1, p is
-- approximately s1*s2. s2 can be thought of as the scale factor, s1 as the
-- largest value of x that is correctly represented. Set s1 to be very large
-- to turn this nonlinear rescaling off, or modify the lines cp=f(xp2) and sp=f(yp2).

s1=0.3
s2=0.5
tanh(x)=(exp(x)-exp(-x))/(exp(x)+exp(-x))
f(x)=s1*s2*tanh(x/s1)

-- make the plot data

cp=f(xp2)
sp=f(yp2)

A[i]=v1[trunc((i-1)/6)] when (i-6*trunc((i-1)/6)=1),
     v2[trunc((i-1)/6)] when (i-6*trunc((i-1)/6)=2),
     v3[trunc((i-1)/6)] when (i-6*trunc((i-1)/6)=3),
     v4[trunc((i-1)/6)] when (i-6*trunc((i-1)/6)=4),
     v5[trunc((i-1)/6)] when (i-6*trunc((i-1)/6)=5) dim[count(xp)*6]

v1[i]={xp[i]-cp[i]/2,yp[i]-sp[i]/2}
v2[i]={xp[i]+cp[i]/2,yp[i]+sp[i]/2}
v3[i]={xp[i]+0.3*cp[i]-0.2*sp[i],yp[i]+0.3*sp[i]+0.2*cp[i]}
v4[i]={xp[i]+0.3*cp[i]+0.2*sp[i],yp[i]+0.3*sp[i]-0.2*cp[i]}
v5[i]={xp[i]+cp[i]/2,yp[i]+sp[i]/2}