[Tex/LaTex] How to draw a vector field (quiver plot) on a surface in 3D

3dpgfplotstikz-pgf

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

Best Answer

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

Related Solutions

[Tex/LaTex] How to draw a quiver plot with varying line thickness

Pgfplots allows to employ point meta as additional input coordinate. Typically, point meta is used as color data, but it can also be used for different purposes.

In this case, you can use it as line thickness parameter for a quiver plot:

\documentclass{article}

\usepackage{pgfplots}

\begin{document}

\begin{tikzpicture}
\def\U{1}
\def\V{2*x}
\def\LEN{(sqrt((\U)^2 + (\V)^2)}
\begin{axis}[axis equal]
    \addplot[blue,
        point meta={\LEN},
        quiver={
            u={(\U)/\LEN},v={(\V)/\LEN},
            scale arrows=2,
            every arrow/.append style={line width=\pgfplotspointmetatransformed/1000 * 2pt},
        },
        -stealth,samples=15,
    ] {x^2};
\end{axis}
\end{tikzpicture}

\end{document}

enter image description here

The example defines some macros \U, \V, and \LEN where \LEN is the length of these vectors. Then, a quiver plot is drawn where each vector is normalized to unit length - and the length parameter is assigned as point meta. This is a simple and standard way of communicating extra data to plots (you could also input extra values in a table).

Finally, \pgfplotspointmetatransformed is a scalar value in the range [0,1000] such that 0 corresponds to the minimum point meta value and 1000 corresponds to the maximum point meta value.

The expression for the line width scales this linearly such that the minimum line thickness is 0 and the maximum is 2pt.

Note that axis equal is necessary if we want to see that the vectors have the same sizes.

[Tex/LaTex] quiver plot on hyperboloid surface

Asymptote is probably better for this, since it allows for hiding the arrows behind the hyperboloid surface, but here's how you can draw the arrows using PGFPlots.

To calculate the tangent vector, you can simply evaluate the y and z values at a location a small distance along the x axis.

\documentclass[12pt]{book}
\usepackage{tikz}
\usepackage{pgfplots}
\pgfplotsset{compat=1.9}

\begin{document}
\begin{tikzpicture}
  \begin{axis}[view={110}{20}, %
      scale = 1.2, y post scale = 1.5,
    xlabel = $x$, ylabel = $y$, zlabel = $z$]
    \addplot3[surf, samples=8, variable = \u, variable y = \v, z buffer = sort,
      y domain = 0:2*pi,
      quiver={
        u={(sqrt(1+(u+0.01)^2)*cos(deg(v)))-x},
        v={0.01},
        w={(sqrt(1+(u+0.01)^2)*sin(deg(v)))-z},
        scale arrows=75
      },
      -stealth, thick]
      ({sqrt(1+u^2)*cos(deg(v))},
      {u},
      {sqrt(1+u^2)*sin(deg(v))}); 
\end{axis}
\end{tikzpicture}
\end{document}

This approach works for other functions as well. You need to make sure to explicitly assign the independent variables of your parametric plot variable names other than x and y, however, otherwise it's not clear whether x refers to the independent variable or to the x coordinate:

\documentclass[border=5mm]{standalone}
\usepackage{tikz}
\usepackage{pgfplots}
\pgfplotsset{compat=1.9}

\begin{document}
\begin{tikzpicture}
  \begin{axis}[view={110}{20}, %
      scale = 1.2, y post scale = 1.5,
    xlabel = $x$, ylabel = $y$, zlabel = $z$]
 \addplot3[surf,domain=1:2, y domain = 0:2*pi, z buffer=sort, samples = 5, samples y=10,
  variable = \s, variable y=\t,
    quiver = {
      u = {(s+0.01)*cos(deg(t)) - x},
      v = {(s+0.01)*sin(deg(t)) - y},
      w = {1/(s+0.01) - z},
      scale arrows=15
    },
    -stealth, thick
  ]
  ({s*cos(deg(t))}, {s*sin(deg(t))}, {1/s});
\end{axis}
\end{tikzpicture}
\end{document}

Best Answer

Related Solutions

[Tex/LaTex] How to draw a quiver plot with varying line thickness

[Tex/LaTex] quiver plot on hyperboloid surface

Related Question