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


Differential.m

From MathBio

Jump to: navigation, search

Differential.m

differential.m is a function that calculated the spring forces given a set of bacteria positions and rest lengths. It is solved in AnimatedPlot.m.


   function dxdt = differential(~, X1)
   global restlengths adjacency springconstants
   le = length(X1);
   le4 = le/4;
   %Reordering the matrix to become a 4 column matrix-------
   for i = 1:1:le4
          X(i,1) = X1(4*i-3);
          X(i,2) = X1(4*i-2);
          X(i,3) = X1(4*i-1);
          X(i,4) = X1(4*i);
   end
   %--------------------------
   %Reformating the adjacency matrix for breakage-------------
   for i = 1:1:le4
       for j = 1:1:le4
           if adjacency(i,j) == 1
               rl = restlength([X(i,1), X(i,2)], [X(j,1),X(j,2)]);
                   if rl > 3*restlengths(i,j)
                       adjacency(i,j) = 0;
                   end
           end
       end
   end
   %-----------------------
   %Finding the spring forces for each value and adding them-------------
   springforce  = zeros(le4,2); 
   kval = 2.2e-5;
   k = 1;
   for i = 1:1:le4
       force = [0,0];
       for j = 1:1:le4
           if adjacency(i,j) == 1
   %           forcecomp = sf([X(i,1), X(i,2)], [X(j,1),X(j,2)],restlengths(i,j),springconstants(i,j));
               forcecomp = sf([X(i,1), X(i,2)], [X(j,1),X(j,2)],restlengths(i,j),kval);
               force(1) = force(1) + forcecomp(1);
               force(2) = force(2) + forcecomp(2);
           end
       end
       springforce(k,1) = force(1);
       springforce(k,2) = force(2);
       k = k+1;
   end
   %-------------------------------
   %Using the drag relationship (constants)-----------------
   mass = 1e-12; %mass of the bacteria in g  
   density = .001; %g/mm^3 
   Cd = 0.48; %drag coefficient for a sphere
   max_velocity = 4250; %mm/s maximum velocity. 4.25e3 mm/s
   radius_tube = 200; %radius of aorta = 12.5mm. radius of capillaries = .004 mm
   radius_bacteria = 1e-3; %radius of bacteria
   Area = pi*radius_bacteria^2;
   c = max_velocity/(radius_tube^2); %Constant representing - 1/4n * deltaP/deltax
   %------------------------------
   %Creating a boundary condition (cutoff length = wall attachment)
   dxdt = zeros(le4,4);
   for i = 1:le4
       if X(i,2) > 2
           dxdt(i,1) = X(i,3);
           dxdt(i,2) = X(i,4);
           %dxdt(i,3) = ((0.5*density*((c*(radius_tube^2 - (radius_tube - X(i,2))^2))^2)*Cd*Area) + springforce(i,1))/mass;
           dxdt(i,3) = (0.5*density*Cd*Area*(max_velocity/(radius_tube^2)*(2*radius_tube*X(i,2)-X(i,2)^2))^2 + springforce(i,1))/mass;
           %dxdt(i,3) = (1e-5 + springforce(i,1))/mass;
           dxdt(i,4) = springforce(i,2)/mass;
       end
   end
   dxdt = reshape(dxdt', le, 1);
   end