%% Numerical simulation of the standard glutamate diffusion with sigma_Dglut=0.
%%
%Critical geometric dimensions involved in simplified synapse geometry.
h=20;%height of cleft in nm.
eps=5;%height of the local volume above the PSD in nm.
Rpsd=100;%radius of the PSD in nm.
Rcleft=200;%radius of the entire cleft space in nm.
Dist_absrptn=40;%distance of the absorbtion boundary from Rcleft in nm.
Rabs=Rcleft+Dist_absrptn;%resultant radius of absorption in nm.

%Defining time-dimension.
Ti=0;Tf=1000;%time scale is in microseconds. So the diffusion is observed for 1 millisecond interval.
Ntme=10001;%No. of temporal discretisations.
dt=(Tf-Ti)/(Ntme-1);%time-step of brownian simulation is 0.1 microsecond.
t=linspace(Ti,Tf,Ntme);%time vector on microseconds scale.

%Glutamate diffusion parameters.
Ntrns=2000;%no. of transmitter molecules released from one vesicle.
D=200;%glutamate diffusion coefficient with units in nm^2.us^-1.
Sqrt_distance=sqrt(2*D*dt);

%Glutamate diffusion coordinates.
X=zeros(Ntrns,Ntme);%x-coordinate of glutamate molecules.
Y=zeros(Ntrns,Ntme);%y-coordinate of glutamate molecules.
Z=zeros(Ntrns,Ntme);%z-coordinate of glutamate molecules.
R=zeros(Ntrns,Ntme);%radial coordinate of glutamate molecules.

%Release site in terms of radial distance from center.
Dist_release=0;%nm away from the center of cleft disc along the radius.
R(:,1)=Dist_release;
X(:,1)=Dist_release*cos(0);
Y(:,1)=Dist_release*sin(0); 

for i=1:Ntme-1
    
    %Normal random number generation.
    Xrandom=randn(Ntrns,1);
    Yrandom=randn(Ntrns,1);
    Zrandom=randn(Ntrns,1);
    
    %Diffusion step.
    X(:,i+1)=X(:,i)+(Sqrt_distance*Xrandom);
    Y(:,i+1)=Y(:,i)+(Sqrt_distance*Yrandom);
    Z(:,i+1)=Z(:,i)+(Sqrt_distance*Yrandom);
    R(:,i+1)=sqrt((X(:,i+1).^2)+(Y(:,i+1).^2));
    
    %Checking for absorption boundary.
    Index1=find(R(:,i+1)>=Rabs);
    X(Index1,i+1)=10000;Y(Index1,i+1)=10000;Z(Index1,i+1)=0;
    R(Index1,i+1)=sqrt((X(Index1,i+1).^2)+(Y(Index1,i+1).^2));
    
    %Checking for upper reflecting boundary. 
    Index1=find(R(:,i+1)<Rabs);
    Index2=find(Z(:,i+1)>0);
    Index3=intersect(Index1,Index2);
    Z(Index3,i+1)=-Z(Index3,i+1);
    
    %Checking for lower reflecting boundary.
    Index2=find(Z(:,i+1)<-h);
    Index3=intersect(Index1,Index2);
    Z(Index3,i+1)=-(2*h)-Z(Index3,i+1);
end

Glu_count=sum(double(R<=Rpsd).*double(Z<=eps-h),1);%recording the temporal profile of glutamate population within the prescribed
                          %local volume of radius Rpsd=100nm and height eps=5nm
                          %above the PSD.

%Result:                         
Glu_conc_Molar=10*(Glu_count/(pi*5*Rpsd^2))/6.023;%converting the temporal profile of glutamate population
                                                 %into the molar concentration.