Take the 2-minute tour ×
TeX - LaTeX Stack Exchange is a question and answer site for users of TeX, LaTeX, ConTeXt, and related typesetting systems. It's 100% free, no registration required.

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: Vector field on torus

share|improve this question
1  
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... –  dingo_d Jul 5 '12 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. –  Håkon Marthinsen Jul 6 '12 at 0:02

2 Answers 2

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)
{
  double rad(R0 + r0*Cos(v));
  return P(rad*Cos(u), rad*Sin(u), r0*Sin(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)
{
  double rad(R0 + r0*Cos(v));
  return P(-rad*Sin(u), rad*Cos(u), 0);
}

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);

  label_mask(White());
  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();
}
share|improve this answer

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

share|improve this answer
    
Not sure the export would be a perspective projection as it is in the provided example. –  pluton Aug 14 '13 at 21:54

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.