Matlab 画箭头矢量场图

                     

贡献者： addis

　　

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

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

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

% 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

% 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

% 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

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