diff --git a/analysis.m b/analysis.m
new file mode 100644
index 0000000000000000000000000000000000000000..e7eb1a2d40db6e60d44553f34eb1927067f86039
--- /dev/null
+++ b/analysis.m
@@ -0,0 +1,80 @@
+
+%% distribution of cluster sizes N(n,t)= mean number of clusters of size n at time t
+clear N;
+N=zeros(700,size(c,2));
+T=1:size(c,2);
+edges=0.5:1:700.5;
+for I=T
+            for i=1:size(c,1)
+                N(:,I)=N(:,I)+histcounts(c{i,I},edges)';
+            end
+end
+    
+I=find(sum(N,2),1,'last');
+N=N(1:I,:)/size(c,1);
+
+    %%
+    figure(1)
+    clf;
+    I=2500;
+    n=1:size(N,1);
+    loglog(n,N(:,I));
+    %image(N)
+
+
+
+%% mean number of clusters greater than k, number of free particles and size of largest cluster
+
+k=10;
+figure(2)
+clf;
+subplot(3,1,1)
+    loglog(t,sum(N(k+1:end,:)))
+    hold on;
+    loglog(t(1500:2500),1000*(t(1500:2500)).^(-0.25))
+    hold off;
+
+subplot(3,1,2)
+    semilogx(t,500-[2:size(N,1)]*N(2:end,:))
+    ylim([0,500])
+    
+subplot(3,1,3)
+    
+    for i=T
+    largest(i)=find(N(:,i),1,'last');
+    end
+    loglog(t,largest)
+    
+%% mean cluster size bigger than k
+k=10;
+n=k+1:size(N,1);
+tmp=sqrt(n/pi)*N(k+1:end,:);
+%tmp=n*N(k+1:end,:);
+R=tmp./sum(N(k+1:end,:));
+
+figure(4)
+loglog(t,R)
+hold on
+range=1600:length(t);
+power=(4/15);
+a=R(range(1))/t(range(1))^(power);
+loglog(t(range),a*t(range).^(power))
+hold off
+%xlim([1,1000])
+%ylim([1,10])
+xlabel('times (MCS)')
+ylabel('linear cluster size')
+
+%%
+%alternative method
+%mean_size=cellfun(@(x) nanmean(x(x>k)),c);
+%figure(5)
+%loglog(t,nanmean(mean_size,1))
+
+%xlim([1,100])
+%ylim([100,1000])
+
+%%
+
+
+    
\ No newline at end of file
diff --git a/include/Particles.h b/include/Particles.h
index 3b4a5a0a3e301a844e077aa603b787caf902904e..4e7ab2cbdb11a4e7eb22df329d753078c877cac5 100644
--- a/include/Particles.h
+++ b/include/Particles.h
@@ -129,6 +129,12 @@ public:
         return ori==1;
 
     }
+
+    inline bool is_diffuse_ind(const int ind) const
+    {
+        return particles[ind]->ori==1;
+    }
+
     inline bool is_bound(const int pos) const
     {
         int ori=0;
@@ -139,8 +145,13 @@ public:
             ori = p-> ori;
         }
         return ori>1;
+    }
 
+        inline bool is_bound_ind(const int ind) const
+    {
+        return particles[ind]->ori>1;
     }
+
     inline int get_orientation(const int pos) const
     {
         int ori=0;
@@ -150,12 +161,24 @@ public:
         return ori-3;
     }
 
+    inline int get_orientation_ind(const int ind) const
+    {
+        return particles[ind]->ori-3;
+    }
+
     inline void set_orientation(const int pos, const int ori)
     {
         assert(!grid1[pos].expired());
         auto p(grid1[pos].lock());
         p-> ori = ori;
     }
+
+    inline void set_orientation_ind(const int ind, const int ori)
+    {
+        particles[ind]-> ori = ori;
+    }
+
+
     inline void set_pos(const int old_pos,const int new_pos, const int ind)
     {
         assert(particles[ind]->pos == old_pos);
@@ -177,7 +200,7 @@ public:
 
         for(unsigned int i=0; i<n.size(); i++)
         {
-            if(int raw_ori=get_ori(n[i].position) && is_interaction_allowed(ori,raw_ori-3,n[i].slope))
+            if(int raw_ori=get_ori(n[i].position); raw_ori>1 && is_interaction_allowed(ori,raw_ori-3,n[i].slope))
             {
                 count_bound++;
             }
@@ -381,7 +404,7 @@ public:
                 interactions.orientations.emplace_back(ori2);
                 interactions.num_bonds++;
                 interactions.num_diffuse++;
-                int bound_neigh_of_neigh = count_interacting_neighbors(n.position,ori);
+                int bound_neigh_of_neigh = count_interacting_neighbors(n.position,ori2);
                 if(bound_neigh_of_neigh>1)
                 {
                     interactions.num_bonds=interactions.num_bonds+(bound_neigh_of_neigh-1);
@@ -424,7 +447,7 @@ public:
             {
 //                std::cout<<"binding"<<"\n";
                 binding_succ++;
-                set_orientation(pos,ori+3);
+                set_orientation_ind(ind,ori+3);
                 for(unsigned int i=0; i<interactions.possible_interaction_pos.size(); i++)
                 {
                     set_orientation(interactions.possible_interaction_pos[i],interactions.orientations[i]+3);
@@ -441,14 +464,13 @@ public:
         // get bound particle
         unbinding_attempt++;
         int particle_pos = get_pos(ind);
-        int ori = get_orientation(particle_pos);
+        int ori = get_orientation_ind(ind);
 
         Interactions interactions;
         interactions.num_bonds=0;
         interactions.num_diffuse=1;
 
-        int pos = get_pos(ind);
-        std::vector<Neighbour>& n(lattice.neighbours[pos]);
+        std::vector<Neighbour>& n(lattice.neighbours[particle_pos]);
 
         auto _is_bound = [this](const Neighbour &n)
         {
@@ -479,7 +501,7 @@ public:
         if (rand<delta_E)
         {
 //            std::cout<<"unbinding"<<"\n";
-            set_orientation(particle_pos,1);
+            set_orientation_ind(ind,1);
             unbinding_succ++;
             if(interactions.possible_interaction_pos.size()>0)
             {
@@ -492,9 +514,11 @@ public:
     }
 
 
-    void label(int ind, int i,std::vector<int> &labels,int &num_bonds) const
+    void label(int ind, int i,std::vector<int> &labels,int &num_bonds,int &num_particles) const
     {
+
         labels[ind]=i;
+        num_particles++;
 
         int pos = get_pos(ind);
         int ori = get_orientation(pos);
@@ -522,10 +546,11 @@ public:
                 int inds = it-particles.begin();
                 return labels[inds]!=0;
             }
+            return false;
         };
 
 
-        auto _label = [this,&labels,i,pos,&num_bonds](const Neighbour &n)
+        auto _label = [this,&labels,i,pos,&num_bonds,&num_particles](const Neighbour &n)
         {
             auto it=(std::find_if(begin(particles),end(particles), [n](std::shared_ptr<Particle> q)
             {
@@ -537,19 +562,21 @@ public:
                 int inds = it-particles.begin();
                 if(labels[inds]==0)
                 {
-                    label(inds,i,labels,num_bonds);
+                    assert(particles[inds]->ori >1);
+                    label(inds,i,labels,num_bonds,num_particles);
                 }
             }
         };
 
 
 
-        auto s1= n | std::views::filter(_is_bound) | std::views::filter(_is_allowed) | std::views::filter(_is_labelled);
-        auto cnt = std::ranges::distance(s1);
+        auto s1= n | std::views::filter(_is_bound) | std::views::filter(_is_allowed);
+
+        auto s2=s1 | std::views::filter(_is_labelled);
+        auto cnt = std::ranges::distance(s2);
         num_bonds=num_bonds+cnt;
 
-        auto s2= n | std::views::filter(_is_bound) | std::views::filter(_is_allowed);
-        std::ranges::for_each(s2,_label);
+        std::ranges::for_each(s1,_label);
     }
 
     std::vector<int> orientations_vec() const
diff --git a/main.cpp b/main.cpp
index 473a53ecc763999030f68f7abbf33d5bb4831ede..aeae9a3c90ef0cbbea4faa1ae625dab724d62d25 100644
--- a/main.cpp
+++ b/main.cpp
@@ -31,14 +31,14 @@ int main(int argc,char *argv[])
     }
     else
     {
-        alpha=0.1;
         J=2.6;
+        alpha=0.1;
         exponent=16;
-        slurm_index = 21;
+        slurm_index = 1001;
         std::cout << "Using default parameters." << '\n';
     }
     density = 0.1;
-    const int MC_steps = 5*pow(10,4); // number of Monte Carlo Steps
+    const int MC_steps = 5*pow(10,5); // number of Monte Carlo Steps
 //    const int MC_steps =500;
     int MC_counter = 0;
 //    long double rand;
@@ -63,13 +63,17 @@ int main(int argc,char *argv[])
 //    std::ofstream out2;
 //    out2.open(fn2.str());
     std::ostringstream fn;
-    fn << slurm_index << "_growth_dFrzB_bondspercluster_J_" << J << "_alpha_" << alpha << "_exp_" << exponent << ".txt";//k_un << "_" << k << ".txt";
+    //fn << slurm_index << "_growth_dFrzB_bondspercluster_J_" << J << "_alpha_" << alpha << "_exp_" << exponent << ".txt";//k_un << "_" << k << ".txt";
+
+    fn << slurm_index << ".txt";//k_un << "_" << k << ".txt";
     std::ofstream out;
     out.open(fn.str());
-    std::ostringstream fn2;
-    fn2 << slurm_index << "_growth_dFrzB_labels_J_" << J << "_alpha_" << alpha << "_exp_" << exponent << ".txt";// k_un << "_" << k << ".txt";
-    std::ofstream out2;
-    out2.open(fn2.str());
+
+//    std::ostringstream fn2;
+//    //fn2 << slurm_index << "_growth_dFrzB_labels_J_" << J << "_alpha_" << alpha << "_exp_" << exponent << ".txt";// k_un << "_" << k << ".txt";
+//    fn2 << slurm_index << "_bonds.txt";//k_un << "_" << k << ".txt";
+//    std::ofstream out2;
+//    out2.open(fn2.str());
 
 //    out << Nx << '\t' << Ny << '\t' <<'\n';
 //    out2 << Nx << '\t' << Ny << '\t' <<'\n';
@@ -127,14 +131,10 @@ int main(int argc,char *argv[])
 //            std::cout<<"ind \t"<<ind<<"\n";
             rand=rand_size-ind;
             assert(rand<1);
-            if(particles.is_diffuse(particles.get_pos(ind)))
+            if(particles.is_diffuse_ind(ind))
             {
 
 //Move diffusive particles
-                if(rand<0)
-                {
-
-                }
                 particles.attempt_diffusion(ind, rand);
 //            print_container(particles.positions);
 
@@ -142,16 +142,26 @@ int main(int argc,char *argv[])
 //BINDING
                 particles.attempt_binding(alpha,J,ind,rand);
 //            print_container(particles.positions);
+            if(particles.is_bound(particles.get_pos(ind))){
 
-
+            assert(particles.get_orientation(particles.get_pos(ind))>-2);
+            assert(particles.get_ori(particles.get_pos(ind))==particles.particles[ind]->ori);
+            assert(particles.count_interacting_neighbors(particles.get_pos(ind),particles.get_orientation(particles.get_pos(ind))));
+            }
 
             }
 
-            else if(particles.is_bound(particles.get_pos(ind)))
+            else if(particles.is_bound_ind(ind))
             {
 // UNBINDING ATTEMPT
+            assert(particles.get_ori(particles.get_pos(ind))==particles.particles[ind]->ori);
+            assert(particles.count_interacting_neighbors(particles.get_pos(ind),particles.get_orientation(particles.get_pos(ind))));
                 particles.attempt_unbinding(alpha,J,ind,rand);
 //            print_container(particles.positions);
+            if(particles.is_bound(particles.get_pos(ind))){
+            assert(particles.get_ori(particles.get_pos(ind))==particles.particles[ind]->ori);
+            assert(particles.count_interacting_neighbors(particles.get_pos(ind),particles.get_orientation(particles.get_pos(ind))));
+            }
 
             }
 
@@ -168,31 +178,42 @@ int main(int argc,char *argv[])
 
 
 
-        if(MC_counter%10000==0)
+        if(MC_counter%(   MC_counter<100 ? 1:  (int) pow(10,floor(log10(MC_counter))-2))==0)
         {
-            std::cout<<Nx<<"num of particles \t"<<particles.particles.size()<<'\n';
+
+            out << MC_counter << '\t';
+            //out2 << MC_counter << '\t';
+
+            //std::cout<<Nx<<"num of particles \t"<<particles.particles.size()<<'\n';
 
 //                std::cout<<"Number of bonds: ";
             std::vector<int> labels(particles.particles.size(),0);
             int label_i = 1;
             for (unsigned int label_index=0; label_index < particles.particles.size(); label_index++)
             {
-                if(particles.is_bound(particles.get_pos(label_index)) && labels[label_index]==0)
+                if(particles.is_bound_ind(label_index) && labels[label_index]==0)
                 {
+                    assert(particles.is_bound(particles.get_pos(label_index)));
+                    assert(particles.get_ori(particles.get_pos(label_index))==particles.particles[label_index]->ori);
                     int num_bonds=0;
-                    particles.label(label_index,label_i,labels,num_bonds);
+                    int num_particles=0;
+                    particles.label(label_index,label_i,labels,num_bonds,num_particles);
                     label_i++;
 //                        std::cout << num_bonds << '\t';
-                    out << num_bonds<< '\t';
+                    assert(num_particles>1);
+                    out << num_particles<< '\t';
+                    //out2 << num_bonds<< '\t';
                 }
             }
 
             //std::cout<<labels.size()<<" "<<particles.positions.size()<<'\n';
-            std::cout<<"Number of clusters: " << std::ranges::max(labels) << '\n';
+            //std::cout<<"Number of clusters: " << std::ranges::max(labels) << '\n';
 //            print_container(labels);
             out << '\n';
-            particles.print_labels(out2,labels);
-            particles.print(out2);
+            //out2 << '\n';
+
+            //particles.print_labels(out2,labels);
+            //particles.print(out2);
         }
         MC_counter++;
     }
diff --git a/main_script.m b/main_script.m
new file mode 100644
index 0000000000000000000000000000000000000000..83516d3bb2b23b3d8960a173157529122932b897
--- /dev/null
+++ b/main_script.m
@@ -0,0 +1,37 @@
+clear;
+
+
+
+
+
+% parfor i=1:36
+%     pause(0.0001*i)
+%     cmd=['set path=%path:C:\Program Files\MATLAB\R2020b\bin\win64;=% & bin\Release\frz_lattice_model.exe 2.6 0.1 16 ',num2str(i)];
+%     system(cmd,'-echo');%outputs to current matlab directory
+% end
+
+
+%%
+for i=1:1000
+    j=1;
+    file=fopen([num2str(i),'.txt']);
+    tline=fgetl(file);
+    while ~isempty(tline) && ischar(tline)
+        z=textscan(tline,'%f','Delimiter','\t');%outputs cell array
+        t(j)=z{1}(1);
+        c{i,j}=z{1}(2:end);
+        tline=fgetl(file);
+        j=j+1;
+    end
+
+fclose(file);
+
+end
+
+save('results.mat','c','t');
+
+
+%%
+
+
+