function handle = arrowhead(x,y,clr,Asize,head)
% arrowhead Draw a arrow head
% This program was written to place arrowheads on trajectories
% calculated using one of MATLAB's ode solvers. However, the program
% does not require the solvers, it only requires two points as shown
% in the figure below.
% (x2,y2) (x2,y2) (x2,y2)
% head=1 o head=2 o head=3 o
% /|\ /|\ /|\
% / | \ / | \ / | \
% / | \ / o \ / | \
% / | \ / / | \ \ / _-o-_ \
% o----o----o o/ | \o o | o
% | | |
% | | |
% (x1,y1) o (x1,y1) o (x1,y1) o
% arrowhead(x,y,color,Asize,head)
% required arguments: the points x=[x1 x2] and y=[y1 y2]
% optional arguments:
% color: a string, same format used in the plot command
% Asize = [Asize(1) Asize(2)]: this rescales the arrowhead
% Asize(1) = the length in changed by this factor
% Asize(2) = the width in changed by this factor
% head: an integer, this specifies the shape of the arrow head
% head = 1 straight back (the default)
% head = 2 slanted back
% head = 3 curved back
% Examples:
% [t,y] = ode45(@rhs,[0 100],[0 1]);
% hold on
% plot(y(:,1),y(:,2))
% i=40; ii=i+3;
% arrowhead([y(i,1) y(ii,1)],[y(i,2) y(ii,2)]);
% arrowhead([y(i,1) y(ii,1)],[y(i,2) y(ii,2)],'r');
% arrowhead([y(i,1) y(ii,1)],[y(i,2) y(ii,2)],'k',[2 1],2);
% arrowhead([y(i,1) y(ii,1)],[y(i,2) y(ii,2)],[],[0.4 1.2]);
% arrowhead([y(i,1) y(ii,1)],[y(i,2) y(ii,2)],[],[],3);
% As indicated in the above three figures, the arrow head is centered on the
% line through (x1,y1), (x2,y2), with the front of the arrow head at (x2,y2).
% In the above examples these points are from the computed solution of an ODE.
% In this context, (x1,y1) is the solution at time t1, while (x2,y2) is the
% solution at time t2, where t1 < t2. depending on how picky you are about
% the arrow head being centered on the solution curve, you can adjust t1 to
% rotate the arrow's direction.
% the default length and width of the arrow head are determined by the height
% and width of the currently active plot, or subplot, window. these values can
% be adjusted using the asize argument.
handle = [];
% errors
if nargin < 2
error('You need to specify two points');
end
if (length(x) ~= 2) || (length(y) ~= 2),
error('x and y vectors must each have two components');
end
if (x(1) == x(2)) && (y(1) == y(2)),
error('The two points are equal - cannot determine direction!');
end
% set the defaults
if nargin <= 2
clr = 'b';
end
if nargin <= 3
Asize = [1 1];
end
if nargin <= 4
head = 1;
end
% check if variables left empty
if isempty(clr)
clr = 'b';
end;
if isempty(Asize)
Asize = [1 1];
end
% determine and remember the hold status, toggle if necessary
if ishold,
WasHold = 1;
else
WasHold = 0;
hold on;
end
x1 = x(1); x2 = x(2); y1 = y(1); y2 = y(2);
% get the axis dimensions
Scales = get(gca,'Position');
% get ranges for axis
TheAxis = axis;
Xlimit = abs(TheAxis(2) - TheAxis(1));
Ylimit = abs(TheAxis(4) - TheAxis(3));
% 2*Whead = width of arrow head
Whead = Asize(2)*min(Scales(3),Scales(4))/50;
% Lhead = length of arrow head
Lhead = 2.8*Whead*Asize(1);
% calculate the two tips of the arrow head
if x1==x2
xa=x1+Xlimit*Whead/Scales(3);
ya=y2-sign(y2-y1)*Ylimit/Scales(4);
xb=x1-Xlimit*Whead/Scales(3);
yb=ya;
else
% calculate slope of arrow head
s0=Scales(4)*Xlimit/(Scales(3)*Ylimit);
m = s0*(y2 - y1)/(x2 - x1);
% calculate center-point of wide end
if x1>x2
x3=x2+Xlimit*Lhead/(Scales(3)*sqrt(1+m^2));
else
x3=x2-Xlimit*Lhead/(Scales(3)*sqrt(1+m^2));
end
y3=y2+m*(x3-x2)/s0;
% calculate the two tips of the arrow head
xa=x3+Xlimit*m*Whead/(Scales(3)*sqrt(1+m^2));
ya=y3-(xa-x3)/(s0*m);
xb=x3-Xlimit*m*Whead/(Scales(3)*sqrt(1+m^2));
yb=y3-(xb-x3)/(s0*m);
end;
% the polygon forming the arrow head
if head==1
xd = [x2 xa xb];
yd = [y2 ya yb];
elseif head==2
LL=0.7*Lhead;
if x1==x2
x4=x2;
y4=y2-sign(y2-y1)*Ylimit*LL/Scales(4);
elseif x1>x2
x4=x2+Xlimit*LL/(Scales(3)*sqrt(1+m^2));
y4=y2+m*(x4-x2)/s0;
else
x4=x2-Xlimit*LL/(Scales(3)*sqrt(1+m^2));
y4=y2+m*(x4-x2)/s0;
end
xd = [x2 xa x4 xb];
yd = [y2 ya y4 yb];
else
z2 = [x2 y2];
z3 = [x3 y3];
zb = [xb yb];
za = [xa ya];
z23=norm(z2-z3);
zb3=norm(zb-z3);
za3=norm(za-z3);
z2b=norm(z2-zb);
z2a=norm(z2-za);
db=(zb-z3);
da=(za-z3);
d=(z2-z3)/z23;
q=(zb-z2)/z2b;
qq=(za-z2)/z2a;
zi = z2 + 0.7*(z2-z3);
z2i=norm(z2-zi);
b=z23/z2i;
a=z2i/zb3^b;
aa=z2i/za3^b;
nn=6;
for i=1:nn
zl=z3+db*(i-1)/(nn-1);
zll=z3+da*(i-1)/(nn-1);
k=z2b*norm(zl-z3)/zb3;
kk=z2a*norm(zll-z3)/za3;
zt=z2+k*q;
ztt=z2+kk*qq;
g=norm(zt-zl)-a*norm(zb-zl)^b;
gg=norm(ztt-zll)-aa*norm(za-zll)^b;
zm=zl+g*d;
zmm=zll+gg*d;
if i==1
xd = [zm(1)]; yd=[zm(2)];
else
xd = [zmm(1) xd zm(1)]; yd=[zmm(2) yd zm(2)];
end;
end;
xd = [x2 xd];
yd = [y2 yd];
end
% draw the arrowhead
fill(xd,yd,clr,'EdgeColor',clr);
% restore original axe ranges and hold status
axis(TheAxis);
if ~WasHold,
hold off
end