%  Daniel Kawano, Rose-Hulman Institute of Technology
%  Last modified:  Feb 23, 2016

function animate_directed_curve(psi, theta, phi, r0, dxi)

%  Solve for how the directors d1 and d2 vary with the convected coordinate
%  xi.  Also, because d3 = d1 x d2 = dr/dxi, after solving for the 
%  evolution of d3, use it to compute the material curve r(xi).  Lastly, 
%  keep track of the tips of the basis vectors d1, d2, and d3 as we move 
%  along the curve:

for k = 1:length(psi);
    P1 = [cos(psi(k)), sin(psi(k)), 0;              %  3-2-3 set of 
          -sin(psi(k)), cos(psi(k)), 0;             %  Euler angles
          0, 0, 1];
    P2 = [cos(theta(k)), 0, -sin(theta(k));
          0, 1, 0;
          sin(theta(k)), 0, cos(theta(k))];   
    P3 = [cos(phi(k)), sin(phi(k)), 0;
          -sin(phi(k)), cos(phi(k)), 0;
          0, 0, 1];
    d1(:,k) = ([1, 0, 0]*(P3*P2*P1))';              %  d1
    d2(:,k) = ([0, 1, 0]*(P3*P2*P1))';              %  d2
    d3(:,k) = ([0, 0, 1]*(P3*P2*P1))';              %  d3
    r(:,k) = r0 + dxi*trapz(d3(:,1:k), 2);          %  m
    r1(:,k) = r(:,k) + d1(:,k);                     %  m
    r2(:,k) = r(:,k) + d2(:,k);                     %  m
    r3(:,k) = r(:,k) + d3(:,k);                     %  m
end

%  Set up the figure window and draw the material curve:

figure
set(gcf, 'color', 'w')
plot3(r(1,:), r(2,:), r(3,:), '-k', 'linewidth', 5)
xlabel('\itx\rm_{1} (m)')
set(gca, 'xdir', 'reverse')
ylabel('\itx\rm_{2} (m)')
set(gca, 'ydir', 'reverse')
zlabel('\itx\rm_{3} (m)            ', 'rotation', 0)
axis equal
xlim([1.2*min(min(min([r1(1,:); r2(1,:); r3(1,:)])), -0.4), 1.2*max(max(max([r1(1,:); r2(1,:); r3(1,:)])), 0.4)])
ylim([1.2*min(min(min([r1(2,:); r2(2,:); r3(2,:)])), -0.4), 1.2*max(max(max([r1(2,:); r2(2,:); r3(2,:)])), 0.4)])
zlim([1.2*min(min(min([r1(3,:); r2(3,:); r3(3,:)])), -0.4), 1.2*max(max(max([r1(3,:); r2(3,:); r3(3,:)])), 0.4)])
grid on

%  Indicate the material curve's starting point r(0), and draw the triad
%  {d1,d2,d3}:

point = line('xdata', r(1,1), 'ydata', r(2,1), 'zdata', r(3,1), 'marker', 'o', 'color', 'k', 'linewidth', 3);
d1vec = line('xdata', [r(1,1), r1(1,1)], 'ydata', [r(2,1), r1(2,1)], 'zdata', [r(3,1), r1(3,1)], 'color', 'b', 'linewidth', 3);
d2vec = line('xdata', [r(1,1), r2(1,1)], 'ydata', [r(2,1), r2(2,1)], 'zdata', [r(3,1), r2(3,1)], 'color', 'r', 'linewidth', 3);
d3vec = line('xdata', [r(1,1), r3(1,1)], 'ydata', [r(2,1), r3(2,1)], 'zdata', [r(3,1), r3(3,1)], 'color', 'g', 'linewidth', 3);

%  Animate the motion of the triad {d1,d2,d3} by updating the figure with 
%  its current location and orientation:

pause

% animation = VideoWriter('directed-curve.avi');
% animation.FrameRate = 1/dxi;
% open(animation);

for k = 1:length(psi)
    set(d1vec, 'xdata', [r(1,k), r1(1,k)], 'ydata', [r(2,k), r1(2,k)], 'zdata', [r(3,k), r1(3,k)]);
    set(d2vec, 'xdata', [r(1,k), r2(1,k)], 'ydata', [r(2,k), r2(2,k)], 'zdata', [r(3,k), r2(3,k)]);
    set(d3vec, 'xdata', [r(1,k), r3(1,k)], 'ydata', [r(2,k), r3(2,k)], 'zdata', [r(3,k), r3(3,k)]);
    drawnow
    % writeVideo(animation, getframe(gcf));
end

% close(animation);