(* ::Package:: *)

function [ax,ay] = QuinticCurve(u,v)
% u and v are column 3-vectors.
% Let
% d12 = distance from (u(1),v(1)) to (u(2),v(2))
% d23 = distance from (u(2),v(2)) to (u(3),v(3))
% d12 = distance from (u(3),v(3)) to (u(1),v(1))
%
% ax and ay are column 5-vectors with the propery that if
%
% px(d) = ax(1) + ax(2)d + ax(3)d^2 + ax(4)d^3 + ax(5)d^4
% py(d) = ay(1) + ay(2)d + ay(3)d^2 + ay(4)d^3 + ay(5)d^4
%
% then
% (px(0),py(0)) = (u(1),v(1))
% (px(d12),py(d12)) = (u(2),v(2))
% (px(d12+d23),py(d12+d23)) = (u(3),v(3))
% (px(d12+d23+d31),py(d12+d23+d31)) = (u(4),v(4))
% (px'(0),py'(0)) = (px'(d12+d23+d31),py'(d12+d23+d31))
% (px "(0),py"(0)) = (px "(d12+d23+d31),py"(d12+d23+d31))

d = zeros(3,1);
d(1) = 0;
d(2) = d(1) + sqrt((u(1)-u(2))^2+(v(1)-v(2))^2);
d(3) = d(2) + sqrt((u(2)-u(3))^2+(v(2)-v(3))^2);
d(4) = d(3) + sqrt((u(3)-u(1))^2+(v(3)-v(1))^2);
ax = Quintic(d,[u;u(1)]);
ay = Quintic(d,[v;v(1)]);

% Display...
t = linspace(0,d(4));
xVals = ((((ax(6)*t + ax(5)).*t + ax(4)).*t + ax(3)).*t + ax(2)).*t + ax(1);
yVals = ((((ay(6)*t + ay(5)).*t + ay(4)).*t + ay(3)).*t + ay(2)).*t + ay(1);
plot(u,v,'ok',u(1),v(1),'or',xVals,yVals,'k')
axis equal off
shg

% If all results are correct: 4 points
% If the results are not correct: 0-3 points, depending on your code

% Displaying the plots correctly: 1 point