☰

如何在 3D 圆环上绘制矢量场（箭头图）？

tag-icon • tag-icon tikz-pgf pgfplots 3d
如何在 3D 圆环上绘制矢量场（箭头图）？

我想使用 PGFPlots 和/或 TikZ 在圆环上绘制矢量场。我设法使用派克斯，基于示例库中的 decorate.xp 文件，但随后我必须手动进行 z 排序，这非常难看。如何让 PGFPlots 和/或 TikZ 给出相同的结果？在我看来，PGFPlots 可以使用绘制矢量场quiver，并将圆环绘制为surf，但不能同时进行。是否可以通过将绘图分成更小的部分来解决这个问题？我该怎么做？

注意：代码是黑客编写的，可能写得和格式都很糟糕。任何改进图表和/或代码的建议都值得赞赏！

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

以下是结果的截图：

环面上的矢量场

答案1

回复不及时，但是为了帮助后人，下面的代码使用了较少的 ePiX 内部构件：

该策略是先绘制表面，然后手动叠加矢量；vec_field()绘制单个箭头，draw_surface()然后叠加。场景被手动切成三个块，以微调隐藏物体的移除（见注释）。

向量场由 指定my_F()，并且可以根据需要更改。返回值应为 u 和 v 的周期对值函数。

安迪

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

答案2

那么将你的 matlab 图形直接导出到 Tikz 怎么样？

然后您可以在生成的 tikz 文件中进行进一步的格式化。

看一眼http://www.mathworks.com/matlabcentral/fileexchange/22022

干杯

相关内容