Matlab 画箭头矢量场图

贡献者： addis

未完成：这里只有单个箭头，需要画出矢量场。

　　这里提供一个在平面上画箭头的 Matlab 函数 arrow，该函数的特点是可以自动检测当前 plot box 的长宽比，使箭头的形状不因此失真。这里是一个演示程序

figure; axis([0, 1000, 0, 1]);
arrow(0, 0, 800, 0.7, 50);

运行结果如下，可以看到虽然该图横轴单位距离比纵轴小很多，但画出的箭头仍不失真。

图 1：运行结果

代码 1：arrow.m

% plot arrow without distortion by getting current aspect ratio
% will 'hold on' and 'axis manual'
function arrow(x, y, vx, vy, headWidth, color, lineWidth)
hold on; axis manual;
A = norm([vx,vy]);
if ~exist('headWidth', 'var') || isempty(headWidth)
    headWidth = 0.2 * A;
end
if ~exist('color', 'var') || isempty(color)
    color = 'k';
end
if ~exist('lineWidth', 'var') || isempty(lineWidth)
    lineWidth = 1;
end
[x1, y1] = xy2equi(x, y);
[vx1, vy1] = xy2equi(vx, vy);
A1 = norm([vx1, vy1]);
headWidth1 = headWidth / A * A1;
Ptri1 = arrowTri(x1, y1, vx1, vy1, headWidth1);
tri_x1 = Ptri1(:,1); tri_y1 = Ptri1(:,2);
[tri_x2, tri_y2] = equi2xy(tri_x1, tri_y1);
vx = vx * (1-headWidth/A*0.3); vy = vy * (1-headWidth/A*0.3);
% plot line
plot([x, x + vx], [y, y + vy], 'Color', color, 'LineWidth', lineWidth);
hold on;
% plot arrow head (triangle)
fill(tri_x2, tri_y2, 'w', 'FaceColor', color, 'EdgeAlpha', 0);
end

% calculate the arrow triangle
function Ptri = arrowTri(x, y, vx, vy, headWidth)
A = norm([vx,vy]);
u = [vx, vy, 0] / A;
Ptri = [0, -1; 0.5, -1; 0, 0; -0.5, -1; 0, -1]; % triangle
Ptri = Ptri * headWidth;
Ptri = [Ptri, zeros(5, 1)];
Ptri = turn(Ptri, [0, 1, 0], u);
Ptri = Ptri(:, 1:2) + [x + vx, y + vy];
end

代码 2：xy2equi.m

% map x, y coordinates to x1, y1 with equal axis
% x range is [0 1], y range [0, ymax] depending on aspect ratial
function [x1, y1] = xy2equi(x, y)
xrange = get(gca, 'XLim');
yrange = get(gca, 'YLim');
aspect = get(gca, 'PlotBoxAspectRatio');
aspect = aspect(2)/aspect(1);

x1 = x / (xrange(2) - xrange(1));
y1 = y / (yrange(2) - yrange(1)) * aspect;
end

代码 3：equi2xy.m

% the inverse of function xy2equi
function [x1, y1] = equi2xy(x, y)
xrange = get(gca, 'XLim');
yrange = get(gca, 'YLim');
aspect = get(gca, 'PlotBoxAspectRatio');
aspect = aspect(2)/aspect(1);

x1 = x * (xrange(2) - xrange(1));
y1 = y * ((yrange(2) - yrange(1))/aspect);
end

代码 4：xy2equi_demo.m

figure; axis([-1000,1000,-1,1]);
set(gca, 'PlotBoxAspectRatio', [0.5, 1, 0.5]);
% plot a square
x = [0 1 1 0 0]/3-0.5/3;
y = [0 0 1 1 0]/3-0.5/3;
[x1, y1] = equi2xy(x, y);
hold on; plot(x1, y1);
[x2, y2] = xy2equi(x1, y1);
% check x2 == x, y2 == y
if norm(x2 - x) + norm(y2 - y) > 1e-14
    error('something wrong');
end

致读者：小时百科一直以来坚持所有内容免费，这导致我们处于严重的亏损状态。长此以往很可能会最终导致我们不得不选择大量广告以及内容付费等。因此，我们请求广大读者热心打赏 ，使网站得以健康发展。如果看到这条信息的每位读者能慷慨打赏 10 元，我们一个星期内就能脱离亏损，并保证在接下来的一整年里向所有读者继续免费提供优质内容。但遗憾的是只有不到 1% 的读者愿意捐款，他们的付出帮助了 99% 的读者免费获取知识，我们在此表示感谢。