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

Question

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 {
class mesh_quad {
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);
        quad(pt1, pt2, pt3, pt4);

        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:
    bool operator() (const mesh_quad& arg1, const mesh_quad& arg2)
    { 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
    std::vector<mesh_quad> mesh;

    for (int i=0; i<N1; ++i)
    for (int j=0; j<N2; ++j)
    mesh.push_back(mesh_quad(F, i*du, j*dv));

    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:

Answer 1

What about directly exporting your matlab figure to Tikz?

Then you can do some further formatting in the generated tikz file.

Take a look at http://www.mathworks.com/matlabcentral/fileexchange/22022

Cheers