September 22, 2017, Friday
University of Colorado at Boulder Search A to Z Campus Map University of Colorado at Boulder CU 
Search Links


AnimatedPlot.m

From MathBio

Jump to: navigation, search

AnimatedPlot.m

AnimatedPlot.m loads a given data set of points and runs a differential equation solver to find out the position of the bacteria points over time. This script plots (and saves) the final result.

   clear all
   clear, clc
   close all
   load X1.mat
   X = X(1:20,:);
   %X = [1, 3; 1, 5; 1, 7; 1, 9; 1,11; 1,13; 2,7'];
   global restlengths adjacency springconstants
   %Normalizing the matrix (min Y value becomes zero) --------
   lengthX = length(X);
   minX = min(X(:,2));
   Xnormal = ones(lengthX,1)*minX;
   X = [X(:,1),X(:,2) - Xnormal];
   %-----------------
   %Adjacency Solving----------
   TRI = delaunay(X(:,1),X(:,2)); 
   dt = DelaunayTri(X(:,1),X(:,2));
   dt.Triangulation;
               nvert = max(max(dt.Triangulation));            
               nface = length(dt.Triangulation);
              adjacency = zeros(nvert);
        for i=1:nface
              adjacency(dt.Triangulation(i,1),dt.Triangulation(i,2)) = 1; adjacency(dt.Triangulation(i,2),dt.Triangulation(i,3)) = 1;
              adjacency(dt.Triangulation(i,3),dt.Triangulation(i,1)) = 1;
             % make sure that all edges are symmetric
              adjacency(dt.Triangulation(i,2),dt.Triangulation(i,1)) = 1; adjacency(dt.Triangulation(i,3),dt.Triangulation(i,2)) = 1;
              adjacency(dt.Triangulation(i,1),dt.Triangulation(i,3)) = 1;
        end
        add = adjacency;
   %End Adjacency solving ----------------
   %Find the rest lengths between connections
   le_adj = length(adjacency);
   restlengths = zeros(le_adj);
   for i = 1:le_adj
       for j = 1:le_adj
           if adjacency(i,j) == 1 && adjacency(j,i) == 1
               rl = restlength([X(i,1), X(i,2)], [X(j,1),X(j,2)]);
               restlengths(i,j) = rl;
               restlengths(j,i) = rl;
           end
       end
   end
   %-------------------------------------------
   %spring length matrix (normalized).
   % kmax = 2.2; %<<<--------------------------------
   % maxlength = max(max(restlengths));
   % rest = restlengths;
   % rest(~rest)=inf;
   % minlength = min(min(rest));
   % knormal = kmax*minlength;
   % springconstants = zeros(le_adj);
   % for i = 1:le_adj
   %     for j = 1:le_adj
   %         if adjacency(i,j) == 1 && adjacency(j,i) == 1
   %             springconstants(i,j) = knormal/restlengths(i,j);
   %         end
   %     end
   % end
   %----------------------------------
   %Add initial velocities (0 m/s)
   X = [X(:,1), X(:,2), zeros(lengthX,1), zeros(lengthX,1)];
   %Using ODE45
   X = reshape(X', lengthX*4, 1);
   [t,z] = ode45(@differential, [0,1], X');
   %--------------------------------------
   le_t = length(t);
   [le_z,wi_z] = size(z);
   for i = 1:1:le_z
       for j = 1:4:wi_z
       x_pos(j - (3/4)*(j-1), i) = z(i, j);
       y_pos(j - (3/4)*(j-1), i) = z(i, j+1);
       end
   end
   maxval = round(le_t/2);
   min_x = round(min(min(x_pos(:,1:maxval))));
   min_y = round(min(min(y_pos(:,1:maxval))));
   max_x = round(max(max(x_pos(:,1:maxval))));
   max_y = round(max(max(y_pos(:,1:maxval))));
   % min_x = min(min(x_pos));
   % min_y = min(min(y_pos));
   % max_x = max(max(x_pos));
   % max_y = max(max(y_pos));
   le_add = length(add);
   % subplot(2,2,l)%
   fig=figure;
   afile = avifile('case5_higherktest');
       for k = 1:maxval
           plot(x_pos(:,k),y_pos(:,k), 'o')
           %axis([min(x_pos(:,k)), max(x_pos(:,k)), min(y_pos(:,k)), max(y_pos(:,k))])
           axis([min_x max_x min_y max_y]) 
           for i = 1:lengthX
               for j = 1:lengthX
                   if add(i,j)&&add(j,i) ==1
                       if restlength([z(k,4*i-3), z(k,4*i-2)], [z(k,4*j-3),z(k,4*j-2)]) >= 3*restlengths(i,j)
                           add(i,j)=0;
                           add(j,i)=0;
                       else
                           line([x_pos(i,k) x_pos(j,k)],[y_pos(i,k) y_pos(j,k)])
                       end
                   end  
               end
           end
           M(k) = getframe;
           afile = addframe(afile,M(k));
       end
   close(fig)
   xlabel('x position')
   ylabel('y position')
   movie(M,10)
   afile = close(afile);
   % k_orig = kval;%
   % end %