%Define number of elements
%Define Initial Positions
clear all;

%Use this code to solve for a range of Pressures and Voltages. 
sign = [-1 1]
V = [100 500 1000];
lim = size(V);
lim = lim(2);
P = [1 5 10];
%pp = 1;

for pp = 1:2
    
for q = 1:lim

%%%%%%%%%%%%%%%%%Define initial shape from catalyst particle%%%%%%%%%%%%
x=[0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20];
y=sin(x/20*pi)/100;
w=zeros(1,21);  %initialize w
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


scale = 1;

deltaT=0.00005*scale;  %time step
gammaS=36; %36;

deltagammaS=0;
deltagammaB=0;
gI=1.01*P(q);gII=0;
m=1;

%start loop
for v=1:150/scale;  %Run for 150 time steps
    v
    clear H;
    %plot current coordinates
    close;
%     figure;plot(x,y,'ko-');
%     axis([-50 70 0 4]);
%     F(v)=getframe;
	filename=['diff',num2str(v)];
	hgsave(eval('filename'));
    %n=length(x)-1;
    n=20;
    %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
    
    %Calling function "FieldFile" to calculate the energy densities at each node
    % due to the electric field i.e to find w1 and w2 for each node.
    xs=x/20;
    ys=y;
    
    %structure for saving all of data - use later to generate plots
    if pp == 1
    s(q,v).x = x;
    s(q,v).y = y;
    else 
    s(q+lim,v).x = x;
    s(q+lim,v).y = y; %+3 = number of voltages
    end
    
     %ys(21) = 0;
     %xs(21) = 1;

% %Check to make sure x not < than 0
%      for jj = 1:21
%          if y(jj) < 0 
%                y(jj) = 0+0.000001*jj;
%          else 
%          end
%      end
    

     %Define w - Call Function to run FEMLAB -> generate w. 
     wout = FieldFile(xs,ys,V(q));
     w = wout{1}/10^3; %scale to 10nm/21pt
     %w = zeros(1,n+1); 
     if pp == 1
     s(q,v).w = w;
     else 
     s(q+lim,v).w = w; %+3 = number of voltages
     end
    
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  %Assemble the H Matrix
    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 with the 'w' term
        
        f(2*k-1)=gammaS*cos(theta(k)) +deltagammaS*(-sin(theta(k)))+l(k)*(gI-gII)/2*(-sin(theta(k)))+f(2*k-1)+sign(pp)*(l(k)/2)*(2/3*w(k)+1/3*w(k+1))*sin(theta(k));
        f(2*k)  =gammaS*sin(theta(k)) +deltagammaS*(cos(theta(k)))+l(k)*(gI-gII)/2*(cos(theta(k)))  +f(2*k)  -sign(pp)*(l(k)/2)*(2/3*w(k)+1/3*w(k+1))*cos(theta(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)+sign(pp)*(l(k)/2)*(1/3*w(k)+2/3*w(k+1))*sin(theta(k));
        f(2*k+2)=gammaS*-sin(theta(k))+deltagammaS*(-cos(theta(k)))+l(k)*(gI-gII)/2*(cos(theta(k))) +f(2*k+2)-sign(pp)*(l(k)/2)*(1/3*w(k)+2/3*w(k+1))*cos(theta(k));
    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    
    H(n+1,n+1)=(1+H(n+1,n+1))*1000000;   
    %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);
        x(i)=x(i)+0.5*rate1(2*i-1)*deltaT;
        y(i)=y(i)+0.5*rate1(2*i)*deltaT;
    end  

end

end
end

save('filename','s') %stores data as structure 's' for later use


% %loading data from file to compile movies
% for i  = 1:150
% load('minus1000_volts_signs')
% ss(1,i).x = s(1,i).x;
% ss(1,i).y = s(1,i).y;
% load('1000_volts_new')
% ss(2,i).x = s(1,i).x;
% ss(2,i).y = s(1,i).y;
% load('100_volts_new')
% ss(3,i).x = s(1,i).x;
% ss(3,i).y = s(1,i).y;
% end
% 

    %plotting nanotube growths side-by-side
    lim = 3;  
    for r = 1:150
    plot(s(1+lim,r).x/2,s(1+lim,r).y,'-o', s(2+lim,r).x/2+20,s(2+lim,r).y,'-o',s(3+lim,r).x/2+40,s(3+lim,r).y,'-o');
    axis([-10 60 0 2]);
    F(r)=getframe;
    title('CNT Growth with Field')
    xlabel('x (nm)')
    ylabel('height (nm)')
    text(0,1.5, 'V = 100V')
    text(20,1.5, 'V = 500V')
    text(40,1.5, 'V = 1000V')
    text(20,1.8, strcat('t =  ', num2str(r*deltaT), 's'))
    F(r)=getframe;
    end
    
    %plotting w-force with nanotube growth
    for r = 1:140
    plot(s(6,r).x/2,1-s(6,r).w/4,s(6,r).x/2,s(6,r).y,'-o',s(1+lim,r).x/2+20,s(3,r).w/4,s(3+lim,r).x/2+20,s(3+lim,r).y,'-o')
    axis([-10 40 0 1.4]);
    text(0,1.2, 'V = -1000V')
    text(20,1.2, 'V = +1000V')
    F(r)=getframe;
    end
    
    %generating move from plots
    movie(F,1,15);% Play the MATLAB movie 
    save movie_field.mat F;% save the MATLAB movie to a file 
    %mpgwrite(F,jet,'movie.mpg'); % Convert the movie to MPEG format 
    movie2avi(F,'movie_field_W_FORCE2.avi');