% 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);