clc
clear all
close all




dist=[5];
burst=[500];
dt=[1000];
ext=40;


for dd=1:size(dist,1)
    for bb=1:size(burst,1)
        for tt=1:size(dt,1)
            
            figure
            grid on
            hold on            
            
            assimilations=[];
            
            B=(burst(bb));
            distance=(dist(dd));
            timestep=(dt(tt));
            
            
            R=2500*10^-9;          % RX radius 
            distance=5013*10^-9;	% TX-RX distance
            
            k=1.38*10^-23;	% Boltzmann's constant
            T=310;		% Temperature Kelvin degrees
            eta=0.0011;	% viscosity
            r=1.75E-9;	% Carrier radius
            
            t=0:0.01:80;  						% time          
            D=(k*T/(6*pi*eta*r));					% Diffusion Coefficient
            
            P=R/(R+distance)*erfc(distance./(2*sqrt(D*t)));	% Theoretical curve - Probability of reception
            theory=P*B;  						% Amount of received carriers
            
         
            plot(t,theory,'b')
            
                
            for ee=0:ext
                
                if(ee>0)
                    est=strcat('.',num2str(ee));
                else
                    est='';
                end
                         
                filename=strcat('nanonet_dist_',num2str(dist(dd)),'_burst_',num2str(burst(bb)),'_dt_',num2str(dt(tt)),'.dat',est);
                                
                if(exist(filename,'file')==2)
                    
                    data=load(filename);
                    
                    s1=size(data,1);
                    s2=size(assimilations,1);
                    
                    
                    if(s2==0)
                        assimilations=data;
                        
                    else if(s1>s2)
                            assimilations=[assimilations,data(1:s2,2)];
                        else
                            assimilations=[assimilations(1:s1,:),data(:,2)];
                        end
                    end
                end
                
                
                plot(data(:,1)*dt*10^-9,data(:,2),'-.');
                
            end
            
            
            res=assimilations(:,2:end);
            
            stt=std(res,0,2);
            mean_assimilations=mean(res,2);
            
            plot(assimilations(:,1)*dt*10^-9,mean_assimilations,'g-');  
            
            
            confidence=95;
            vall=1-(1-confidence/100)/2;
            for i=1:size(res,1),
                appo=tinv(vall,size(res,2)-1);
                confa(i)=appo*stt(i)/sqrt(size(res,1));
            end
            
            errorbar(assimilations(:,1)*dt*10^-9, mean_assimilations, confa,'r:')
            axis([0 2.5 0 150]);
            xlabel('time [s]');
            ylabel('Assimilations');

        end
    end
end