
%Visualization of chain_growth. Monomers are shown in red, growing ends are
%shown in green, polymers are shown in blue. The polymer chains are frozen
%in place, which is not physically realistic. Chains grow when a monomer
%hits the gorwing end. Set the parameters and run the script. Explainations
%of each parameter is next to the parameter. To cancel the simulation,
%press ctrl-c when focused on the MATLAB command window. Some non-physical
%thngs that currently happen include: 1) the frozen polymers trapping
%monomers that don't bond, 2) instead of requiring very small timesteps,
%particles that overlap with others will be forcefully moved to a random
%part of the simulation. This allows the simulation to be run in real time
%though, so it's probably worth it. 

box_side_length=80; %Side the side length of the box containing particles. Monomers have a diameter of 1 unit. 
number_of_particles=400; %Number of monomers.
step_size=0.1; %Timestep size. Best to keep it at 0.1, or particles will start to go through each other and pop out of eistence. A small bump can be okay on weaker systems.
timesteps=5000; %Number of timesteps in the simulation.
number_of_chains=5; %Number of growing polymer chains. 
termination_chance=1/200; %Chance of a growing chain terminating each timestep. A new chain will initiate at a the location of another random monomer. 
replace_particle_flag=0; %When set to 1, all particles that polymerize will be replaced with a particle in a random location. When set to 0, then will be removed. 

%%Run

sidelength=box_side_length;
balls=number_of_particles;
timestep=step_size;
numsteps1=0;
numsteps2=timesteps;
numinit=number_of_chains;
newinitfrac=termination_chance;
replacement=replace_particle_flag;

% sidelength=80;
% balls=400*4;
% timestep=0.1;
% numsteps1=5;
% numsteps2=50000;
% numinit=5;
% newinitfrac=1/10000;
% replacement=0;


f=figure



pos=rand(1,2)*sidelength;

for i=2:balls
    i
    test=0;
    while test==0
        postest=rand(1,2)*sidelength;
        dsquared=sum((postest-pos).^2,2);
        if sum(dsquared<1.5)>0
            test=0;
        else
            pos(end+1,:)=postest;
            test=1;
        end
    end
end

% pos=rand(balls,2)*sidelength;
vel=-2*rand(balls,2)+1;
f.Units='inches'
f.Position=[0.6528 3.1042 8.7569 7.0694];
axis([0 sidelength 0 sidelength])
f.Units='points';
f.Position(4)=f.Position(3);
a=gca;
a.Units='points'
a.Position(4)=a.Position(3);
diameter=a.Position(3)/sidelength
area=(diameter/2)^2*pi*1;
b=scatter(pos(:,1),pos(:,2),area,'r');
f.Color='w';
hold on
% box on
plot([0-0.5 sidelength+0.5 sidelength+0.5 -0.5 -0.5],[-0.5 -0.5 sidelength+0.5 sidelength+0.5 -0.5],'k-');
set(gca,'visible','off')
axeslimits=[a.XLim a.YLim];
hold off
diameter=a.Position(3)/axeslimits(2);
area=(diameter/2)^2*pi*1;
b=scatter(pos(:,1),pos(:,2),area,'r');
hold on
plot([0-0.5 sidelength+0.5 sidelength+0.5 -0.5 -0.5],[-0.5 -0.5 sidelength+0.5 sidelength+0.5 -0.5],'k-');
hold off
set(gca,'visible','off')
pause(0.1)


timestepcounter=0;
%warmup
for j=1:numsteps1
    j

    pos=pos+vel*timestep;
    d1=pdist(pos)<1 & pdist(pos)>0;
    d2=triu(squareform(d1));
    %bouncing code
    for i=find(d2)'
        [row,col]=ind2sub(size(d2),i);
        if row~=col
            v1=vel(row,:);
            v2=vel(col,:);
            x1=pos(row,:);
            x2=pos(col,:);

            v1new=v1-(dot(v1-v2,x1-x2)/(norm(x1-x2)^2))*(x1-x2);
            v2new=v2-(dot(v2-v1,x2-x1)/(norm(x2-x1)^2))*(x2-x1);
            vel(row,:)=v1new;
            vel(col,:)=v2new;

            if isnan(v1new) | isnan(v2new)
                error('vel is nan')
            end
        end
    end
    
    %x outside box
    for i=find(pos(:,1)<0)
        vel(i,1)=-vel(i,1);
        pos(i,1)=-pos(i,1);
    end
    for i=find(pos(:,1)>sidelength)
        vel(i,1)=-vel(i,1);
        pos(i,1)=2*sidelength-pos(i,1);
    end
    
    %y outside box
    for i=find(pos(:,2)<0)
        vel(i,2)=-vel(i,2);
        pos(i,2)=-pos(i,2);
    end
    for i=find(pos(:,2)>sidelength)
        vel(i,2)=-vel(i,2);
        pos(i,2)=2*sidelength-pos(i,2);
    end
%     pos=pos+vel*timestep;
    if sum(sum(isnan(pos)))~=0
    error('particle Nan position')
    end
    
    
    scatter(pos(:,1),pos(:,2),area,'r')
%     axis([0 sidelength 0 sidelength])
    axis(axeslimits)
    hold on
    plot([0-0.5 sidelength+0.5 sidelength+0.5 -0.5 -0.5],[-0.5 -0.5 sidelength+0.5 sidelength+0.5 -0.5],'k-');
    set(gca,'visible','off')
    hold off
    
    drawnow
    filename=['video_out_' num2str(zeros(1,5-numel(num2str(timestepcounter)))')' num2str(timestepcounter) '.png'];
%     export_fig filename -png% -r150
%     movefile('filename.png',filename)
    timestepcounter=timestepcounter+1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start

inits=randi(balls,numinit,1);

initpos=pos(inits,:);
% initpos1=initpos2(1,:);
% initpos2=initpos2(1,:);
% initpos3=initpos2(1,:);
% initpos4=initpos2(1,:);
% initpos1=initpos2(1,:);

initcurrentpos=pos(inits,:);
newinitcounttracker=0;
pos(inits,:)=nan;
scatter(pos(:,1),pos(:,2),area,'r')
hold on
scatter(initpos(1,1,:),initpos(1,2,:),area,'b')
% axis([0 sidelength 0 sidelength])
axis(axeslimits)
plot([0-0.5 sidelength+0.5 sidelength+0.5 -0.5 -0.5],[-0.5 -0.5 sidelength+0.5 sidelength+0.5 -0.5],'k-');
set(gca,'visible','off')
hold off
timestepcounter=1;
for j=1:numsteps2
    j
    pos=pos+vel*timestep;
    
    %test for new init
    a=rand;
    if a<newinitfrac
        newinitcounttracker=newinitcounttracker+1;
        if newinitcounttracker>numinit
            newinitcounttracker=1;
        end
        initcurrentpos(newinitcounttracker,:)=pos(randi(balls,1),:);
    end
    
    
    %find next bond
    for i=1:numinit
        dsquared=sum((initcurrentpos(i,:)-pos).^2,2);
        partlist=setdiff(find(dsquared<1),inits);
        if length(partlist)>0
            newinit=partlist(1);
%             hold on
%             scatter(pos(newinit,1),pos(newinit,2),area,'k')
            initpos(end+1,:)=pos(newinit,:);
            inits(i)=newinit;
            initcurrentpos(i,:)=pos(newinit,:);
            
            if replacement==1
                test=0;
                while test==0
                    postest=rand(1,2)*sidelength;
                    dsquared=[sum((postest-pos).^2,2); sum((postest-initpos).^2,2)];
                    if sum(dsquared<1.5)>0
                        test=0;
                    else
    %                     pos(end+1,:)=postest;
                        pos(newinit,:)=postest;
                        vel(newinit,:)=1-2*rand(1,2);
                        test=1;
                    end
                end
            else
                pos(newinit,:)=pos(newinit,:)*nan;
            end
                
%             pos(newinit,:)=rand(1,2)*sidelength*nan;
%             vel(newinit,:)=[0 0];
        end
    end
%vector to add to bounce off polymers
  
    allpos=[pos;initpos];
%     allpos=pos;
    
    d1=pdist(allpos)<1 & pdist(allpos)>0;
    d2=triu(squareform(d1));
    %bouncing code
    for i=find(d2)'
        [row,col]=ind2sub(size(d2),i);
        if row~=col & (row<=balls & col<=balls)
            v1=vel(row,:);
            v2=vel(col,:);
            x1=pos(row,:);
            x2=pos(col,:);

            v1new=v1-(dot(v1-v2,x1-x2)/(norm(x1-x2)^2))*(x1-x2);
            v2new=v2-(dot(v2-v1,x2-x1)/(norm(x2-x1)^2))*(x2-x1);
            vel(row,:)=v1new;
            vel(col,:)=v2new;

            if isnan(v1new) | isnan(v2new)
                error('vel is nan')
            end     
        elseif row~=col & (row<=balls) %account for hitting polymers
            v1=vel(row,:);
            x1=allpos(row,:);
            x2=allpos(col,:);
            v2=[0 0];
            v1new=v1-(dot(v1-v2,x1-x2)/(norm(x1-x2)^2))*(x1-x2);
            v2new=v2-(dot(v2-v1,x2-x1)/(norm(x2-x1)^2))*(x2-x1);
            vel(row,:)=-v2new;
        elseif row~=col & (col<=balls) %account for hitting polymers
            v2=vel(col,:);
            v1=[0 0];
            x1=allpos(row,:);
            x2=allpos(col,:);
            v1new=v1-(dot(v1-v2,x1-x2)/(norm(x1-x2)^2))*(x1-x2);
            v2new=v2-(dot(v2-v1,x2-x1)/(norm(x2-x1)^2))*(x2-x1);
            vel(row,:)=-v2new;
        end
            
    end
    
    %%%%%%%%%%%%%%%%%%%%%% test too close
    d1=pdist(allpos)<0.85 & pdist(allpos)>0;
    d2=triu(squareform(d1));
    %bouncing code
    for i=find(d2)'
        [row,col]=ind2sub(size(d2),i);
        if row~=col & (row<=balls & col<=balls)
            test=0;
            while test==0
                postest=rand(1,2)*sidelength;
                dsquared=[sum((postest-pos).^2,2); sum((postest-initpos).^2,2)];
                if sum(dsquared<1.5)>0
                    test=0;
                else
%                     pos(end+1,:)=postest;
                    pos(row,:)=postest;
                    vel(row,:)=1-2*rand(1,2);
                    test=1;
                end
            end
        elseif row~=col & (row<=balls)
            test=0;
            while test==0
                postest=rand(1,2)*sidelength;
                dsquared=[sum((postest-pos).^2,2); sum((postest-initpos).^2,2)];
                if sum(dsquared<1.5)>0
                    test=0;
                else
%                     pos(end+1,:)=postest;
                    pos(row,:)=postest;
                    vel(row,:)=1-2*rand(1,2);
                    test=1;
                end
            end
        elseif row~=col & (col<=balls)
            test=0;
            while test==0
                postest=rand(1,2)*sidelength;
                dsquared=[sum((postest-pos).^2,2); sum((postest-initpos).^2,2)];
                if sum(dsquared<1.5)>0
                    test=0;
                else
%                     pos(end+1,:)=postest;
                    pos(col,:)=postest;
                    vel(col,:)=1-2*rand(1,2);
                    test=1;
                end
            end
        end
    end
    %%%%%%%%%%%%%%%%%%%%%%%%end too close
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     d1=pdist(allpos)<1 & pdist(allpos)>0;
%     d2=triu(squareform(d1));
% %     d1=pdist(allpos)<1 & pdist(allpos)>0;
% %     d2=triu(squareform(d1));
%     %bouncing code    
%     for i=find(d2)'
%         [row,col]=ind2sub(size(d2),i);
%         if row~=col & (row<=balls | col<=balls)
%             if row<=balls
%                 v1=vel(row,:);
%                 x1=pos(row,:);
%             else
%                 v1=[0,0];
%                 x1=allpos(row,:);
%             end
%             if col<=balls
%                 v2=vel(col,:);
%                 x2=pos(col,:);
%             else
%                 v2=[0,0];
%                 x2=allpos(col,:);
%             end
%             
%             if row<=balls & col<=balls
%                 v1new=v1-(dot(v1-v2,x1-x2)/(norm(x1-x2)^2))*(x1-x2);
%                 v2new=v2-(dot(v2-v1,x2-x1)/(norm(x2-x1)^2))*(x2-x1);
%             elseif row<=balls
%                 v1newa=v1-(dot(v1-v2,x1-x2)/(norm(x1-x2)^2))*(x1-x2);
%                 v2newa=v2-(dot(v2-v1,x2-x1)/(norm(x2-x1)^2))*(x2-x1);
%                 v1new=-v2new;
%             else
%                 v1newa=v1-(dot(v1-v2,x1-x2)/(norm(x1-x2)^2))*(x1-x2);
%                 v2newa=v2-(dot(v2-v1,x2-x1)/(norm(x2-x1)^2))*(x2-x1);
%                 v2new=-v1new;
%             end
% %             if row<=balls
% %                 vel(row,:)=v1new;
% %             end
% %             if col<=balls
% %                 vel(col,:)=v2new;
% %             end
% 
%             if isnan(v1new) | isnan(v2new)
%             error('vel is nan')
%             end
%         end
%     end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %x outside box
    for i=find(pos(:,1)<0)
        vel(i,1)=-vel(i,1);
        pos(i,1)=-pos(i,1);
    end
    for i=find(pos(:,1)>sidelength)
        vel(i,1)=-vel(i,1);
        pos(i,1)=2*sidelength-pos(i,1);
    end
    
    %y outside box
    for i=find(pos(:,2)<0)
        vel(i,2)=-vel(i,2);
        pos(i,2)=-pos(i,2);
    end
    for i=find(pos(:,2)>sidelength)
        vel(i,2)=-vel(i,2);
        pos(i,2)=2*sidelength-pos(i,2);
    end
    
    if sum(sum(isnan(pos)))~=0
%         error('particle Nan position')
    end
    
    scatter(pos(:,1),pos(:,2),area,'r')
    hold on
    scatter(initpos(:,1),initpos(:,2),area,'b')
    scatter(initcurrentpos(:,1),initcurrentpos(:,2),area,'g')
    plot([0-0.5 sidelength+0.5 sidelength+0.5 -0.5 -0.5],[-0.5 -0.5 sidelength+0.5 sidelength+0.5 -0.5],'k-');
    set(gca,'visible','off')
    hold off
%     axis([0 sidelength 0 sidelength])
    axis(axeslimits)
    drawnow
    filename=['video_out_' num2str(zeros(1,5-numel(num2str(timestepcounter)))')' num2str(timestepcounter) '.png'];
%     export_fig filename -png% -r150
%     movefile('filename.png',filename)
%     saveas(f,filename);
    F=getframe(f);
    F2=F.cdata;
    F3=F2(220:end-220,240:end-200,:);
    F.cdata=F3;
    F4(timestepcounter) = F ;
    timestepcounter=timestepcounter+1;
end

% 
% 
writerObj = VideoWriter('myVideo.avi');
writerObj.FrameRate = 10;
% open the video writer
open(writerObj);
% write the frames to the video
for i=1:length(F4)
    % convert the image to a frame
    frame = F4(i) ;    
    writeVideo(writerObj, frame);
end
% close the writer object
close(writerObj);
