%Define number of elements
%Define Initial Positions
clear all;
x=[0 1 2 3 4 5 6 7 8 9 9 9 10 11 11 11 12 13 14 15 16 17 18 19 20];
y=[5 5 5 5 5 5 5 5 5 5 4 3 3 3 4 5 5 5 5 5 5 5 5 5 5];
plot(x,y,'ko-')
%define properties
deltaT=0.0005;
gammaS=1;
gammaB=1;
deltagammaS=0;
deltagammaB=0;
gI=1;gII=1;
m=1;
%start loop
for v=1:1000;
    %plot current coordinates
    close;
    figure;plot(x,y,'ko-');
    axis([0 20 2 6]);
    F(v)=getframe;
	filename=['diff',num2str(v)];
	hgsave(eval(['filename']));
    n=length(x)-1;
    %assembling H matrix
    %combine each node's H matrix into global H matrix
    H=[zeros(2*n+2,2*n+2)];
    f=[zeros(2*n+2,1)];
    %define elements
    %calculate theta(i) and l(i)
    for k=1:n
        l(k)=((y(k+1)-y(k))^2+(x(k+1)-x(k))^2)^(1/2);
        theta(k)=asin((y(k+1)-y(k))/l(k));
        if x(k+1)<x(k)
        theta(k)=pi-theta(k);
        end
    end
   for k=1:n;
        %First Row
        H(2*k-1,2*k-1)=l(k)/6/m*2*(sin(theta(k)))^2+H(2*k-1,2*k-1); 
        H(2*k-1,2*k)=l(k)/6/m*-2*sin(theta(k))*cos(theta(k))+H(2*k-1,2*k);
        H(2*k-1,2*k+1)=l(k)/6/m*(sin(theta(k)))^2+H(2*k-1,2*k+1); 
        H(2*k-1,2*k+2)=l(k)/6/m*-sin(theta(k))*cos(theta(k))+H(2*k-1,2*k+2); 
        %Second Row
        H(2*k,2*k-1)=l(k)/6/m*-2*sin(theta(k))*cos(theta(k))+H(2*k,2*k-1);
        H(2*k,2*k)=l(k)/6/m*2*(cos(theta(k)))^2+H(2*k,2*k);
        H(2*k,2*k+1)=l(k)/6/m*-sin(theta(k))*cos(theta(k))+H(2*k,2*k+1);
        H(2*k,2*k+2)=l(k)/6/m*(cos(theta(k)))^2+H(2*k,2*k+2);
        %Third Row
        H(2*k+1,2*k-1)=l(k)/6/m*(sin(theta(k)))^2+H(2*k+1,2*k-1); 
        H(2*k+1,2*k)=l(k)/6/m*-sin(theta(k))*cos(theta(k))+H(2*k+1,2*k);
        H(2*k+1,2*k+1)=l(k)/6/m*2*(sin(theta(k)))^2+H(2*k+1,2*k+1); 
        H(2*k+1,2*k+2)=l(k)/6/m*-2*sin(theta(k))*cos(theta(k))+H(2*k+1,2*k+2); 
        %Fourth Row
        H(2*k+2,2*k-1)=l(k)/6/m*-sin(theta(k))*cos(theta(k))+H(2*k+2,2*k-1);
        H(2*k+2,2*k)=l(k)/6/m*(cos(theta(k)))^2+H(2*k+2,2*k);
        H(2*k+2,2*k+1)=l(k)/6/m*-2*sin(theta(k))*cos(theta(k))+H(2*k+2,2*k+1);
        H(2*k+2,2*k+2)=l(k)/6/m*2*(cos(theta(k)))^2+H(2*k+2,2*k+2);
        %assembling f column
        f(2*k-1)=gammaS*cos(theta(k))+deltagammaS*(-sin(theta(k)))+l(k)*(gI-gII)/2*(-sin(theta(k)))+f(2*k-1);
        f(2*k)=gammaS*sin(theta(k))+deltagammaS*(cos(theta(k)))+l(k)*(gI-gII)/2*(cos(theta(k)))+f(2*k);
        f(2*k+1)=gammaS*-cos(theta(k))+deltagammaS*(sin(theta(k)))+l(k)*(gI-gII)/2*(-sin(theta(k)))+f(2*k+1);
        f(2*k+2)=gammaS*-sin(theta(k))+deltagammaS*(-cos(theta(k)))+l(k)*(gI-gII)/2*(cos(theta(k)))+f(2*k+2);
    end
        %Add Artificial Finite Mobility (small number to diagonal points in H matrix)
    for m=1:(2*n+2);
        H(m,m)=0.000001+H(m,m);
    end
    %boundary conditions
    %fix end points
	H(2*n+1,2*n+1)=(1+H(2*n+1,2*n+1))*1000000;
	H(2*n+2,2*n+2)=(1+H(2*n+2,2*n+2))*1000000;
    H(1,1)=(1+H(1,1))*1000000;
	H(2,2)=(1+H(2,2))*1000000;
	%solve for velocities
	rate1=H\f;
    %update positions
	for i=1:(n+1);
        x1(i)=x(i)+0.5*rate1(2*i-1)*deltaT;
        y1(i)=y(i)+0.5*rate1(2*i)*deltaT;
	end
    %recalc with new positions
    %assembling H matrix
    %combine each node's H matrix into global H matrix
    H=[zeros(2*n+2,2*n+2)];
    f=[zeros(2*n+2,1)];
    %define elements
    %calculate theta(i) and l(i)
    for k=1:n
        l(k)=((y1(k+1)-y1(k))^2+(x1(k+1)-x1(k))^2)^(1/2);
        theta(k)=asin((y1(k+1)-y1(k))/l(k));
        if x1(k+1)<x1(k)
        theta(k)=pi-theta(k);
        end
    end
   for k=1:n;
        %First Row
        H(2*k-1,2*k-1)=l(k)/6/m*2*(sin(theta(k)))^2+H(2*k-1,2*k-1); 
        H(2*k-1,2*k)=l(k)/6/m*-2*sin(theta(k))*cos(theta(k))+H(2*k-1,2*k);
        H(2*k-1,2*k+1)=l(k)/6/m*(sin(theta(k)))^2+H(2*k-1,2*k+1); 
        H(2*k-1,2*k+2)=l(k)/6/m*-sin(theta(k))*cos(theta(k))+H(2*k-1,2*k+2); 
        %Second Row
        H(2*k,2*k-1)=l(k)/6/m*-2*sin(theta(k))*cos(theta(k))+H(2*k,2*k-1);
        H(2*k,2*k)=l(k)/6/m*2*(cos(theta(k)))^2+H(2*k,2*k);
        H(2*k,2*k+1)=l(k)/6/m*-sin(theta(k))*cos(theta(k))+H(2*k,2*k+1);
        H(2*k,2*k+2)=l(k)/6/m*(cos(theta(k)))^2+H(2*k,2*k+2);
        %Third Row
        H(2*k+1,2*k-1)=l(k)/6/m*(sin(theta(k)))^2+H(2*k+1,2*k-1); 
        H(2*k+1,2*k)=l(k)/6/m*-sin(theta(k))*cos(theta(k))+H(2*k+1,2*k);
        H(2*k+1,2*k+1)=l(k)/6/m*2*(sin(theta(k)))^2+H(2*k+1,2*k+1); 
        H(2*k+1,2*k+2)=l(k)/6/m*-2*sin(theta(k))*cos(theta(k))+H(2*k+1,2*k+2); 
        %Fourth Row
        H(2*k+2,2*k-1)=l(k)/6/m*-sin(theta(k))*cos(theta(k))+H(2*k+2,2*k-1);
        H(2*k+2,2*k)=l(k)/6/m*(cos(theta(k)))^2+H(2*k+2,2*k);
        H(2*k+2,2*k+1)=l(k)/6/m*-2*sin(theta(k))*cos(theta(k))+H(2*k+2,2*k+1);
        H(2*k+2,2*k+2)=l(k)/6/m*2*(cos(theta(k)))^2+H(2*k+2,2*k+2);
        %assembling f column
        f(2*k-1)=gammaS*cos(theta(k))+deltagammaS*(-sin(theta(k)))+l(k)*(gI-gII)/2*(-sin(theta(k)))+f(2*k-1);
        f(2*k)=gammaS*sin(theta(k))+deltagammaS*(cos(theta(k)))+l(k)*(gI-gII)/2*(cos(theta(k)))+f(2*k);
        f(2*k+1)=gammaS*-cos(theta(k))+deltagammaS*(sin(theta(k)))+l(k)*(gI-gII)/2*(-sin(theta(k)))+f(2*k+1);
        f(2*k+2)=gammaS*-sin(theta(k))+deltagammaS*(-cos(theta(k)))+l(k)*(gI-gII)/2*(cos(theta(k)))+f(2*k+2);
    end
        %Add Artificial Finite Mobility (small number to diagonal points in H matrix)
    for m=1:(2*n+2);
        H(m,m)=0.000001+H(m,m);
    end
    %boundary conditions
    %fix end points
	H(2*n+1,2*n+1)=(1+H(2*n+1,2*n+1))*1000000;
	H(2*n+2,2*n+2)=(1+H(2*n+2,2*n+2))*1000000;
    H(1,1)=(1+H(1,1))*1000000;
	H(2,2)=(1+H(2,2))*1000000;
    
    %solve for velocities
	rate=H\f;
    %update positions
	for i=1:(n+1);
        x(i)=x(i)+rate(2*i-1)*deltaT;
        y(i)=y(i)+rate(2*i)*deltaT;
	end
  
%check element size
    %Check for elements that are too small
	for j=1:2
	n=length(x);
	tr=0;
	for k=1:(n-2);
        l(k)=((y(k+1)-y(k))^2+(x(k+1)-x(k))^2)^(1/2);
        if l(k)<0.2
            tr=tr+1;   
            for j=(k+1):(length(x)-1);
                x(j)=x(j+1);
                y(j)=y(j+1);
            end;
        end;
        if (n-2-tr)<=k, break, end
	end;
        
	xr=x;yr=y;
	x=0;y=0;
	for k=1:(length(xr)-tr)
        x(k)=xr(k); y(k)=yr(k);
	end
	end
	
	%check last point
	lpc=0;
	for m=(length(x)-1)
        l(m)=((y(m+1)-y(m))^2+(x(m+1)-x(m))^2)^(1/2);
		if l(m)<0.2
            lpc=1;
            x(m)=x(m+1);
            y(m)=y(m+1);
		end
	end
	
	xlp=x; ylp=y;
	x=0; y=0;
	for m=1:(length(xlp)-lpc)
        x(m)=xlp(m); y(m)=ylp(m);
	end
	
	%Check for elements that are too large
	sn=length(x)-1;
	nn=sn;
	for k=1:(length(x)-1);
        l(k)=((y(k+1)-y(k))^2+(x(k+1)-x(k))^2)^(1/2);
        if l(k)>5;
            nn=nn+1;
        end;
	end;
	for k=1:sn;
        l(k)=((y(k+1)-y(k))^2+(x(k+1)-x(k))^2)^(1/2);
        if l(k)>5;
            for c=(length(x)):-1:(k+1);
                y(c+1)=y(c); x(c+1)=x(c);
            end;
        x(k+1)=(x(k)+x(k+1))/2; y(k+1)=(y(k)+y(k+1))/2;
        end;
	end;
	for k=sn:nn;
        l(k)=((y(k+1)-y(k))^2+(x(k+1)-x(k))^2)^(1/2);
        if l(k)>5;
            for c=(length(x)):-1:(k+1);
                y(c+1)=y(c); x(c+1)=x(c);
            end;
        x(k+1)=(x(k)+x(k+1))/2; y(k+1)=(y(k)+y(k+1))/2;
        end;
	end;
end
movie2avi(F,'webhole','fps',40,'quality',90);


 

