# How can I draw a vector field (quiver plot) on a surface in 3D?

I'd like to draw a vector field on a torus using PGFPlots and/or TikZ. I managed to do it using ePiX, based on the file decorate.xp from the example gallery, but then I have to do z-sorting by hand, which is very ugly. How can I get PGFPlots and/or TikZ to give the same result? It seems to me that PGFPlots can plot the vector field using `quiver`, and plot the torus as a `surf`, but not both at once. Is it possible to get around this by chopping up the drawing into smaller pieces? How should I do that?

NB: The code is a hack, and is probably badly written and formatted. Any tips for improving the figure and/or code are appreciated!

``````/* -*-ePiX code based on on the example file decorate.xp-*- */
#include <algorithm>
#include "epix.h"
using namespace ePiX;

const int N1(60);  // latitudes
const int N2(60);  // longitudes

const double du(2*M_PI/N1);
const double dv(2*M_PI/N2);

const double r_0(.999); // minor radius
const double R_0(1);    // major radius

P F(double u, double v)
{
return polar(R_0 + r_0*Cos(u), v) + P(0,0,r_0*Sin(u));
}

namespace ePiX {
private:
P pt1, pt2, pt3, pt4, center;
double distance;

public:
mesh_quad(P f(double u, double v), double u0, double v0)
: pt1(f(u0,v0)), pt2(f(u0+du,v0)), pt3(f(u0+du,v0+dv)), pt4(f(u0,v0+dv)),
center(0.25*(pt1 + pt2 + pt3 + pt4)),
distance(max(max(max(norm(pt1-camera.viewpt()), norm(pt2-camera.viewpt())),
norm(pt3-camera.viewpt())), norm(pt4-camera.viewpt())))
{ }

double how_far() const { return distance; }

void draw() const
{
P direction(center-camera.viewpt());
P normal((pt2 - pt1)*(pt4 - pt1));
normal *= 1/norm(normal);

fill(Gray(normal|(recip(distance)*direction)));

pen(0);

fill(false);
}
};

class mesh_vf {
private:
P pt1, pt2, center;
double distance;

public:
mesh_vf(P f(double u, double v), double u0, double v0)
{
pt1 = f(u0,v0);

double C = 2 * (pt1.x2() - 8);
double Gamma = -C * Sin(u0) * Sin(v0);
double Delta = C * (1 + Cos(u0)) * Cos(v0);

pt2 = f(u0,v0) + .005 * P(Gamma * Sin(u0) * Cos(v0)
+ Delta * (1 + Cos(u0)) * Sin(v0),
Gamma * Sin(u0) * Sin(v0) - Delta * (1 + Cos(u0)) * Cos(v0),
-Gamma * Cos(u0));

center = .5*(pt1 + pt2);
distance = min(norm(pt1-camera.viewpt()),norm(pt2-camera.viewpt()));
//distance = norm(center-camera.viewpt());
}

double how_far() const { return distance; }

void draw() const
{
plain(Black());
arrow(pt1, pt2, .3);
}
};

class by_distance {
public:
{ return arg1.how_far() > arg2.how_far(); }
bool operator() (const mesh_vf& arg1, const mesh_vf& arg2)
{ return arg1.how_far() > arg2.how_far(); }
bool operator() (const mesh_quad& arg1, const mesh_vf& arg2)
{ return arg1.how_far() > arg2.how_far(); }
bool operator() (const mesh_vf& arg1, const mesh_quad& arg2)
{ return arg1.how_far() > arg2.how_far(); }
};
} // end of namespace

int main()
{
picture(P(-2.5,-3.5),P(3.5,3), "12x10cm");

begin();

set_crop();

viewpoint(sph(10, 80 * M_PI / 180, 30 * M_PI / 180));

// background grids
pen(Gray(.5));
grid(P(-2, -2, -2), P(2, 2, -2), 4,4);
grid(P(-2, -2, -2), P(2, -2, 2), 4,4); // n.b. (z,x) divisions
grid(P(-2, -2, -2), P(-2, 2, 2), 4,4);

label(P(2, 0, -2), P(-15, -15), "\$y\$", bl);
label(P(2, -2, -2), P(-2, -2), "\$-2\$", bl);
label(P(2, -1, -2), P(-2, -2), "\$-1\$", bl);
label(P(2, 0, -2), P(-2, -2), "\$0\$", bl);
label(P(2, 1, -2), P(-2, -2), "\$1\$", bl);

label(P(-2, 2, 0), P(15, -15), "\$z\$", br);
label(P(-2, 2, -1), P(2, 0), "\$-1\$", br);
label(P(-2, 2, 0), P(2, 0), "\$0\$", br);
label(P(-2, 2, 1), P(2, 0), "\$1\$", br);
label(P(-2, 2, 2), P(2, 0), "\$2\$", br);

label(P(0, 2, -2), P(0, -20), "\$x\$", b);
label(P(-1, 2, -2), P(0, -3), "\$-1\$", b);
label(P(-2, 2, -2), P(0, -3), "\$-2\$", b);
label(P(0, 2, -2), P(0, -3), "\$0\$", b);
label(P(1, 2, -2), P(0, -3), "\$1\$", b);
label(P(2, 2, -2), P(0, -3), "\$2\$", b);

// build and draw a torus with vector field

for (int i=0; i<N1; ++i)
for (int j=0; j<N2; ++j)

sort(mesh.begin(), mesh.end(), by_distance());

std::vector<mesh_vf> vf;

for (int i=0; i<N1; ++i)
for (int j=0; j<N2; ++j)
vf.push_back(mesh_vf(F, i*du, j*dv));
// vf.push_back(mesh_vf(F, (i+.5)*du, (j+.5)*dv));

sort(vf.begin(), vf.end(), by_distance());

int j = 0;
for (unsigned int i=0; i<mesh.size(); ++i)
{
mesh.at(i).draw();
if (j < vf.size())
{
for ( ; (vf.at(j).how_far() > mesh.at(i).how_far()) && (j < vf.size()); j++)
vf.at(j).draw();
}
}
for ( ; j < vf.size(); j++)
vf.at(j).draw();

//label_color(Red());
//spot(P(0,2,0));

tikz_format();
end();
}
``````

Here is a screenshot of the result:

• Wow, this looks impressive, but wouldn't it be easier to use maybe Mathematica or similar software to generate an image and then export it as .eps or whatever format you need? Because this looks kinda daunting... Jul 5, 2012 at 10:52
• Yes, it's very easy to do it in Matlab (that was actually the starting point), but the resulting figure gets messed up when you save it as .eps or .pdf. I don't have any experience with Mathematica. However, I prefer using free software. Jul 6, 2012 at 0:02

Not a timely reply, but in case it's helpful to posterity, the code below uses fewer ePiX internals:

The strategy is to plot the surface first, then to overlay the vectors manually; `vec_field()` draws individual arrows, `draw_surface()` does the overlaying. The scene is manually chopped into three blocks to fine-tune the hidden object removal (see comments).

The vector field is specified by `my_F()`, and may be changed as needed. The return value should be a periodic, pair-valued function of u and v.

Andy

``````/* -*-ePiX-*- */
#include "epix.h"
using namespace ePiX;

const double MAX(2); // outside radius of torus

const int N1(40);  // longitudes
const int N2(20);  // latitudes

const double r0(0.95); // minor radius
const double R0(MAX - r0);    // major radius

domain dom(P(0, 0), P(360, 360), mesh(N1, N2), mesh(3*N1, 3*N2));

// torus parametrization
P Phi(double u, double v)
{
}

// vector field to plot (in plane coordinates)
P my_F(double u, double v)
{
P loc(Phi(u, v));
double x(loc.x1()), y(loc.x2()), z(loc.x3());
return 0.05*P((x - y)*(x + y), 2*x*y); // adjust as desired
}

// auxiliary mappings
// partial derivatives (manually computed)
P Phi_u(double u, double v)
{
}

P Phi_v(double u, double v)
{
return P(-r0*Cos(u)*Sin(v), -r0*Sin(u)*Sin(v), r0*Cos(v));
}

// draw vector if surface element faces camera
const double du(360/N1);
const double dv(360/N2);

void vec_field(double u, double v)
{
double u1(u + 0.5*du), v1(v + 0.5*dv); // middle of mesh element
P loc(Phi(u1, v1)),
Xu(Phi_u(u1, v1)), Xv(Phi_v(u1, v1)), dir(Xu*Xv); // cross product

if (0 < ((camera.viewpt() - loc)|dir))
{
P tmp(my_F(u1, v1)), // vector field in coordinates
vec(tmp.x1()*Xu +  tmp.x2()*Xv); // field in R^3

dart(loc, loc + vec);
}
}

// will chop scene manually; group repeated drawing commands
void draw_surface()
{
pen(White(), 0);
surface(Phi, dom, -1); // cull back-facing elements

plain(Black());
for (int i = 0; i < 2*N1; ++i)
for (int j = 0; j < 2*N2; ++j)
vec_field(0.5*i*du, 0.5*j*dv);
}

int main()
{
picture(P(-3, -2.5),P(3, 2.5), "12x10cm");

begin();
degrees();
camera.at(sph(10, 80, 30));

P O(0, 0, 0), // origin
back(-MAX, -MAX, -MAX); // back corner

// axes and labels
arrow(back, back + (2*MAX + 1)*E_1);
label(back + (2*MAX + 1)*E_1, P(-4, 0), "\$x\$", l);

arrow(back, back + (2*MAX + 1)*E_2);
label(back + (2*MAX + 1)*E_2, P(0, -4), "\$y\$", b);

arrow(back, back + (2*MAX + 1)*E_3);
label(back + (2*MAX + 1)*E_3, P(0, 4), "\$z\$", t);

axis aX(P(-MAX,  MAX, -MAX), P( MAX, MAX, -MAX), 2*MAX, P(0, -4), b);
axis aY(P( MAX, -MAX, -MAX), P( MAX, MAX - 1, -MAX), 2*MAX - 1, P(-2, 0), l);
axis aZ(P(-MAX,  MAX, -MAX + 1), P(-MAX, MAX,  MAX), 2*MAX - 1, P(0, -2), br);

aX.draw_labels();
aY.draw_labels();
aZ.draw_labels();

// background grids
black(0.5);
grid(P(-MAX, -MAX, -MAX), P(MAX,  MAX, -MAX), 2*MAX, 2*MAX);
grid(P(-MAX, -MAX, -MAX), P(MAX, -MAX,  MAX), 2*MAX, 2*MAX);
grid(P(-MAX, -MAX, -MAX), P(-MAX, MAX,  MAX), 2*MAX, 2*MAX);

black();

fill(White());
// manually chop the scene
clip_face(O, -E_3); // bottom half
draw_surface();
clip_restore();

clip_face(O, E_3); // top back quarter
clip_face(O, -E_2);
draw_surface();
clip_restore();

clip_face(O, E_3); // top front quarter
clip_face(O, E_2);
draw_surface();
clip_restore();

tikz_format();
end();
}
``````