%                __________________________________________________________
%               %\                                                       
%    welcome    %%\            
%               %%%\                                                             
%%%%%%%%%%%%%%%%%%%%\______________________________________________________
                                                                  clear all
                                                                  clear all
                                                                        clc
%%% input parameters %%\

N=500;
px_size= 0.0506;                           %%%%%%%%%%%%%%% number of frames
f=0.2041;                                  %%%%%%%%%%% time lapse (s/frame)
cutoff=8;                                  %% please insert only 2^n values
TauLimit=int64(N/10)+1;                    %%%%%%%%%%%%%%% maximum time-lag

Filter=0;                                  %% bkg subtraction (0=no; 1=yes)
av_toll=10000;                              %% threshold for bkg subtraction
bit=16;                                    % image-type (bit=8, bit=16 ...) 

%_________________________________________________________________________%
%                                  import image-stack   \%%%%%%%%%%%%%%%%%%
%%%    string before indexing      \%%%%%%%%%%%%%%%%%
name_start='PhogrinGFP_preCHOL.lif_Series001_t';    
%%%     string after indexing      \%%%%%%%%%%%%%%%%
name_end='_ch00';            
%___________________________________________________________ _____________            
                                                           %\             |
                                                            %\    START   |
                                                             %\___________|
disp('     importing image-stack')
disp(' ')
for k=1:N
if k<11
inputFileName = sprintf([name_start '00%d' name_end '.tif'],k-1);
else
if k<101
inputFileName = sprintf([name_start '0%d' name_end '.tif'],k-1);
else
if k<1001
inputFileName = sprintf([name_start '%d' name_end '.tif'],k-1);
else
end
end
end
A(:,:,k) = imread(inputFileName);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  starting fitting parameters %%%%%%%%%%%

%linear stics [scale D offset]
    %start2=[1 0.005 1];
%parabolic stics [scale D v offset]
    start3=[1 0.005 0.01 1];
%linear imsd
    %startS2=[D_1 omega^2];
    startS2=[0.001 1];
%parabolic imsd
    %startS=[D_1 v_1 omega^2];
    startS=[0.001 0.005 1];

           %%% set rho_\infty = cont_0*omega+fact*vdrift*TauLimit %%%%%%%%%
cont_0=1;  
fact=0;
    
K=double(TauLimit);

 P=size(A,1);
 R=px_size*P;
 A=double(A);
 BB=A;
 off2=mean(mean(mean(A)));
 K1=K*ones(N,1);
 K2=K1;
 for k=1:K
     K1(k)=k-1;
     K2(N-k+1)=k-1;
 end
   A7=A;
   if Filter==0
 else
  A5=A-av_toll;
  A5(find(A5<0))=0;
  A5=A5+av_toll;
  A5(find(A5==av_toll))=0;
  A=A5;
  end
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% plot1      

dim=0:px_size:R-px_size;
orig_h0 = hist (A7(:,:,:),50);
check_H0 = sum(sum(orig_h0,2),1);
orig_H0 = sum(orig_h0,2);
filt_h0 = hist (A(:,:,:),50);
filt_H0 = sum(filt_h0,2);
max_bar=max(orig_H0);

figure (1)

subplot (2,3,1)
pcolor (dim,dim,A7(:,:,1))
axis square
shading interp
colormap default
title ('Raw image','fontsize',14)
%caxis ([0 255])
xlabel ('x (\mum)','fontsize',14)
ylabel ('y (\mum)','fontsize',14)
set (gca,'fontsize',14)

subplot (2,3,3)
bar (2^bit/50*(2:50),orig_H0(2:end))
axis square
xlabel ('Intensisty','fontsize',14)
ylabel ('Counts','fontsize',14)
% set(gca,'YScale','log')
%ylim ([0 max_bar*(1.05)])
xlim ([0 2^bit])
set (gca,'fontsize',14)

subplot (2,3,2)
surf (dim,dim,A7(:,:,1))
axis square
shading interp
xlim ([0 R])
ylim ([0 R])
%zlim ([-10 260])
box on
xlabel ('x (\mum)','fontsize',14)
ylabel ('y (\mum)','fontsize',14)
zlabel ('Intensisty','fontsize',14)
set (gca,'fontsize',14)

subplot (2,3,4)
pcolor (dim,dim,A(:,:,1))
axis square
shading interp
title ('Bkg subtraction','fontsize',14)
%caxis ([0 255])
xlabel ('x (\mum)','fontsize',14)
ylabel ('y (\mum)','fontsize',14)
set (gca,'fontsize',14)

subplot (2,3,6)
bar (2^bit/50*(2:50),filt_H0(2:end))
axis square
xlabel ('Intensisty','fontsize',14)
ylabel ('Counts','fontsize',14)
xlim ([0 2^bit])
set (gca,'fontsize',14)

subplot (2,3,5)
surf (dim,dim,A(:,:,1))
axis square
shading interp
xlim ([0 R])
ylim ([0 R])
xlabel ('x (\mum)','fontsize',14)
ylabel ('y (\mum)','fontsize',14)
zlabel ('Intensisty','fontsize',14)
set (gca,'fontsize',14)
%zlim ([-10 260])
box on
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% spatial autocorrelation function %%%%%%%%%

a = (-R/2:px_size:(R-px_size)/2);
b = (-R/2:px_size:(R-px_size)/2);

disp ('     calculating g(\xi, eta, 0)_n ...')


for k=1:N
    prc = sprintf('%d',k);

    disp (['     computing g(\xi, eta, 0)_' prc])
    scale(k)=mean(mean(A(:,:,k)))*P;
    %Z(:,:,k)=abs(fftshift(ifft2(fft2(A(:,:,k)).*conj(fft2(A(:,:,k))))))./(scale(k)^2*double(TauLimit));
    Z(:,:,k)=(abs(fftshift(ifft2(fft2(A(:,:,k)).*conj(fft2(A(:,:,k))))))-mean(mean(A(:,:,k))))./(scale(k)^2*double(TauLimit));
    ZB(:,:,k)=(abs(fftshift(ifft2(fft2(BB(:,:,k)).*conj(fft2(BB(:,:,k))))))-mean(mean(BB(:,:,k))))./(scale(k)^2*double(TauLimit));
end

Z_mean=mean(Z,3);

a1 = (-P/cutoff*px_size:px_size:(P-1)/cutoff*px_size);
b1 = a1;
%Z_cut=Z(P/2+1+int64(a1/px_size),P/2+1+int64(b1/px_size),:);
%Z1=Z_mean(P/2+1+int64(a1/px_size),P/2+1+int64(b1/px_size));
Z2=Z(P/2+1+int64(a1/px_size),P/2+1+int64(a1/px_size),:);
A2=A(P/2+1+int64(a1/px_size),P/2+1+int64(b1/px_size),:);
ZB2=ZB(P/2+1+int64(a1/px_size),P/2+1+int64(a1/px_size),:);

for i=1:N
startGxy=[10 0.85 0.0];
lGxy=[0 0 -10];
uGxy=[inf inf inf];
[X Y]=meshgrid(a1,a1);
xy(:,:,1)=X;
xy(:,:,2)=Y;
options=optimset('display','off');
modelFunGxy=@(pgxy,xy)(pgxy(1)*exp(-((xy(:,:,1)).^2+(xy(:,:,2)).^2)./pgxy(2).^2)+pgxy(3));
%options = optimset('Display','off');
 Gxy=lsqcurvefit(modelFunGxy,startGxy,xy,Z2(:,:,i),lGxy,uGxy,options);
 Gaussxy(:,:,i)=modelFunGxy(Gxy,xy);
 s_Gxy(i,1)=Gxy(2);
 o_Gxy(i,1)=Gxy(3);
 %s_G(1,1)=omega^2;
 Gxy_x(:,i)= Gaussxy(:,length(a1)/2+1,i);
 Gxy_y(:,i)= Gaussxy(length(a1)/2+1,:,i);
 
end
Z1=mean(Z2,3);
for i=1:N
    Gx(:,i)= Z2(:,length(a1)/2+1,i);
    Gy(:,i)= Z2(length(a1)/2+1,:,i);
end

omega2=mean(s_Gxy,1);
err_omega2=std(s_Gxy,1);
omega=omega2;
err_omega=err_omega2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% temporal correlation function %%%%

t_lg=TauLimit;

[X,Y,num_frame] = size(A);

STCorr = zeros(X,Y,TauLimit);  
ST_std = zeros(X,Y,TauLimit);
SeriesMean = squeeze(mean(mean(A)));
imgser_fluct = zeros(X,Y,num_frame);
for i = 1:num_frame
    imgser_fluct(:,:,i) = A(:,:,i) - SeriesMean(i);
end

for tau = 0 : TauLimit-1
clc
disp(['                                               iMSD   '])
disp ('     input parameters')
disp (['         image size (px)     = ' sprintf('%d*%d',P,P)])
disp (['         pixel size (\mum)   = ' sprintf('%.3f',px_size)])
disp (['         time lapse (s)      = ' sprintf('%.3f',f)])
disp (['         roi_stics size (px) = ' sprintf('%d*%d',2*P/cutoff,2*P/cutoff)])
disp (['         lag_time (s)        = ' sprintf('%.2f',(TauLimit-1).*f)]);
disp (['         n_frames            = ' sprintf('%d',N)])
disp (' ')
prc = sprintf('%.d',(100*(tau)/(TauLimit-1)));
disp ('     fitting g(\xi, \eta, 0)                   done')
disp (['         omega^2 (\mum^2)    = ' sprintf('%.6f', omega^2) ' +- ' sprintf('%.6f',2*err_omega*omega)])
disp (' ')
disp (['     computing g(\xi, \eta, \tau)              ' prc '%'])
    
    lagcorr = zeros(X,Y,(num_frame-tau));
    for pair=1:(num_frame-tau)
        lagcorr(:,:,pair) = fftshift(real(ifft2(fft2(imgser_fluct(:,:,pair)).*conj(fft2(imgser_fluct(:,:,(pair+tau)))))));
        lagcorr(:,:,pair) = lagcorr(:,:,pair)/SeriesMean(pair)/SeriesMean(pair+tau);
    end
    STCorr(:,:,(tau+1)) = mean(lagcorr,3);
    ST_std(:,:,(tau+1)) = std(lagcorr,0,3);
end

STCorr_a=STCorr(P/2+1+int64(a1/px_size),P/2+1+int64(b1/px_size),:);

W(:,1)=STCorr(P/2+1,P/2+1,:)/(P^2*double(t_lg));
t(:,1)=double(0:t_lg-1);
 start2=[1 0.0025 1 0.05 1];                                                                          t=f*t;
    
    t1=t(2:end);
    W1=W(2:end);
 
STCorr=STCorr/(P^2*double(t_lg));
STCorr_a=STCorr(P/2+1+int64(a1/px_size),P/2+1+int64(b1/px_size),:);

                                                                           %  %  %  %  % % iMSD
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:t_lg
    B=STCorr_a(:,:,i);
    [value, location] = max(B(:));
    [R,C] = ind2sub(size(B),location);
    RR(i,1)=R;
    CC(i,1)=C;
end

o2=0;
o1=0;
ws = warning('off','all');

for i=1:t_lg
Iy(:,i)= STCorr_a(:,CC(i,1),i);
Ix(:,i)= STCorr_a(RR(i,1),:,i);   
modelFunyy = @(py,a1)(py(1).*exp(-(a1-py(4)).^2./py(2).^2)+o2);
    options = statset('nlinfit');
    %options.Weigths=[zeros(T) s0];
    options.RobustWgtFun= 'bisquare';
    options.Display= 'Off';
   
startG=[1 0 0 1 0];
startG1=[1 0 0 1 0];
startG12=[1 0 0 1 0 0];
lG=[0 -100 -100 0 0 -10];
uG=[100 100 100 100 100 100];
[X Y]=meshgrid(a1,a1);
xy(:,:,1)=X;
xy(:,:,2)=Y;
options=optimset('display','off');

 if i==1
 modelFunG=@(pg,xy)(pg(1)*exp(-((xy(:,:,1)-pg(2)).^2+(xy(:,:,2)-pg(3)).^2)./pg(4).^2)+pg(5));
 %options = optimset('Display','off');
 G=lsqcurvefit(modelFunG,startG,xy,STCorr_a(:,:,i),[],[],options);
 Gaussb(:,:,1)=modelFunG(G,xy);
 x_G(i,1)=G(2);
 y_G(i,1)=G(3);
 s_G(i,1)=G(4).^2;
 %s_G(1,1)=omega^2;
 o_G(i,1)=G(5);
 O_G=o_G;
 else
 modelFunG=@(pg,xy)(pg(1)*exp(-((xy(:,:,1)-pg(2)).^2+(xy(:,:,2)-pg(3)).^2)./pg(4).^2)+O_G);
   G=lsqcurvefit(modelFunG,startG12,xy,STCorr_a(:,:,i),[lG uG],[],options);
  Gaussb(:,:,i)=modelFunG([G(1) -G(2) G(3) G(4)],xy);
 x_G(i,1)=-G(2);
 y_G(i,1)=G(3);
 s_G(i,1)=G(4).^2;
 end
Gauss_y(:,i)=Gaussb(:,size(STCorr_a,1)/2+1+int64(y_G(i,1)*px_size),i);
Gauss_x(:,i)=Gaussb(size(STCorr_a,1)/2+1+int64(x_G(i,1)*px_size),:,i);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% drift fitting %%%%%%%%%%%%%%

startX=[0 0];
modelFunX = @(pX,t)(pX(1).*t+pX(2));
options = statset('nlinfit');
options.RobustWgtFun= 'bisquare';
t_1=t(1:t_lg);
[wnlmX, RX,JX,CovX] = nlinfit(t_1,x_G,modelFunX,startX,options);
errX=sqrt(diag(CovX));
rsqX = 1-(sum(RX.^2,1) ./ sum((x_G - mean(x_G)).^2,1));
resultsX=[wnlmX' errX];
diff=omega^2/(4*wnlmX(2));
fitFunX=modelFunX(wnlmX,t_1);
v_driftX=resultsX(1,1);
err_v_driftX=resultsX(1,2);

startY=[0 0];
modelFunY = @(pY,t)(pY(1).*t+pY(2));
options = statset('nlinfit');
options.RobustWgtFun= 'bisquare';
t_1=t(1:t_lg);
[wnlmY, RY,JY,CovY] = nlinfit(t_1,y_G,modelFunY,startY,options);
errY=sqrt(diag(CovY));
rsqY = 1-(sum(RY.^2,1) ./ sum((y_G - mean(y_G)).^2,1));
resultsY=[wnlmY' errY];
diff=omega^2/(4*wnlmY(2));
fitFunY=modelFunX(wnlmY,t_1);
v_driftY=resultsY(1,1);
err_v_driftY=resultsY(1,2);

v_drift1=sqrt(v_driftX^2+v_driftY^2);
err_v_drift1=(2*err_v_driftX*v_driftX+2*err_v_driftY*v_driftY)/(2*v_drift1);


t_1=t(2:end);
v_drift=sqrt(x_G(end)^2/(double(t_lg))^2+y_G(end)^2/(double(t_lg))^2);

v_drift_x=x_G(end)/(double(t_lg));
v_drift_y=y_G(end)/(double(t_lg));

if v_driftY>0
    ALPHA=acos(v_driftX/v_drift1);
else
    ALPHA=-acos(v_driftX/v_drift1);
end
   
Sigma=s_G;
Sigma_1=Sigma(2:end);
 y_th(1:t_lg,1)=px_size;
 square = [-px_size -px_size; -px_size px_size; px_size px_size; px_size -px_size; -px_size -px_size];
  
disp(' ')
disp (['         v_\phi (\mum/s)    =  ' sprintf('%.6f', v_drift1) ' +- ' sprintf('%.6f',err_v_drift1)])
disp (['         \phi (rad)         =  ' sprintf('%.3f', ALPHA)])
disp ('     _____________')
disp ('     iMSD analysis')
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% fitting Sigma
 
  options = statset('nlinfit');
  options.Robust= 'Off ';
 
lb=[0 0.00005 0];
ub=[10 0.1 1];

t_2=t(2:end);
Sigma2=Sigma(2:end);
lb=[0 0 0 ];
ub=[1 1 10 ];
modelFunS = @(pS)(pS(3)+4*pS(1).*t_2+pS(2).^2.*t_2.^2-Sigma2);
 [SS,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(modelFunS,startS,lb,ub,options);
 dfe = length(t_2)- 3;
 [Q,R] = qr(jacobian,0);
 Rinv = inv(R);
 se = sqrt(sum(Rinv.*Rinv,2)*resnorm/dfe);
rsqS=1-(sum(residual.^2,1) ./ sum((Sigma2 - mean(Sigma2)).^2,1));
fitFunS=modelFunS(SS)+Sigma2;

par=zeros(3,1);
par(1)=se(1);
par(2)=se(2);
par(3)=se(3);
D_Sigma=SS(1);
v_Sigma=SS(2);
err_D_Sigma=par(1);
err_v_Sigma=par(2);

lb=[0.00005 0];
ub=[0.1 10];
modelFunS2 = @(pS2)(pS2(2)+4*pS2(1).*t_2-Sigma2);
 [SS2,resnorm2,residual2,exitflag,output,lambda,jacobian2]=lsqnonlin(modelFunS2,startS2,lb,ub,options);
 dfe2 = length(t_2)- 3;
 [Q2,R2] = qr(jacobian2,0);
 R2inv = inv(R2);
 se2 = sqrt(sum(R2inv.*R2inv,2)*resnorm2/dfe2);
rsqS2=1-(sum(residual2.^2,1) ./ sum((Sigma2 - mean(Sigma2)).^2,1));
fitFunS2=modelFunS2(SS2)+Sigma2;

par2=zeros(2,1);
par2(1)=se2(1);
par2(2)=se2(2);
D_Sigma2=SS2(1);
err_D_Sigma2=par2(1);

v_tot_Sigma=sqrt(v_Sigma^2+v_drift1^2);
err_v_tot_Sigma=1/(2*v_tot_Sigma)*(2*v_drift1*err_v_drift1+2*v_Sigma*err_v_Sigma);
A_simm=(v_Sigma/v_tot_Sigma)^2;
A_asimm=(v_drift1/v_tot_Sigma)^2;
a_simm=sqrt(A_simm);
a_asimm=sqrt(A_asimm);
factors=[A_asimm A_simm; a_asimm a_simm];

theta1=(0:0.005:2*pi);
ftheta1=(2./theta1.*sin(theta1./2)).^2;
g=A_asimm;
vec=find(g>=ftheta1);
beta=theta1(vec(1));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% anomalous diffusion
lbc=[10^-4 10^-4 0];
ubc=[Inf Inf 2];
startc=[0 10^-3 0.8];

modelFun_conf = @(pS)(pS(1)+4*pS(2).*t_2.^pS(3)-Sigma2);
 
 [SSc,resnormc,residualc,exitflagc,outputc,lambdac,jacobianc]=lsqnonlin(modelFun_conf,startc,lbc,ubc,options);
 dfec = length(t_2)-3;
 [Qc,Rc] = qr(jacobianc,0);
 Rinvc = inv(Rc);
 sec = sqrt(sum(Rinvc.*Rinvc,2)*resnormc/dfec);
rsqS_conf=1-(sum(residualc.^2,1) ./ sum((Sigma2 - mean(Sigma2)).^2,1));
fitFun_conf=modelFun_conf(SSc)+Sigma2;
par_c=zeros(3,1);
par_c(1)=sec(1);
par_c(2)=sec(2);
par_c(3)=sec(3);
pow=SSc(3);
err_pow=par_c(3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% confined diffusion 
lbc2=[0 10^-4 0 0];
ubc2=[Inf 10^-1 Inf 8];
startc2=[0 10^-3 0.5 1];

modelFun_conf2 = @(pC)(pC(1)+4*pC(2).*t_2+pC(3)./3.*(1-exp(-t_2./pC(4)))-Sigma2);
[SSc2,resnormc2,residualc2,exitflagc2,outputc2,lambdac2,jacobianc2]=lsqnonlin(modelFun_conf2,startc2,lbc2,ubc2,options);
 dfec2 = length(t_2)-4;
 [Qc2,Rc2] = qr(jacobianc2,0);
 Rinvc2 = inv(Rc2);
 sec2 = sqrt(sum(Rinvc2.*Rinvc2,2)*resnormc2/dfec2);
rsqS_conf2=1-(sum(residualc2.^2,1) ./ sum((Sigma2 - mean(Sigma2)).^2,1));
fitFun_conf2=modelFun_conf2(SSc2)+Sigma2;
par_c2=zeros(4,1);
par_c2(1)=sec2(1);
par_c2(2)=sec2(2);
par_c2(3)=sec2(3);
par_c2(4)=sec2(4);
D_M=SSc2(2);
err_D_M=par_c2(2);
D_m=SSc2(3)./(12.*SSc2(4))+D_M;
err_D_m=err_D_M+(1/12)*SSc2(3)/SSc2(4)*(par_c2(3)/SSc2(3)+par_c2(3)/SSc2(4));
Lquadro=SSc2(3);
tauc=SSc2(4);

disp(' ')
disp ('       -  Brownian diffusion')
disp (['         offset (\mum^2)    =  ' sprintf('%.6f', SS2(2)) ' +- ' sprintf('%.6f',par2(2))])
disp (['         D (\mum^2/s)       =  ' sprintf('%.6f', D_Sigma2) ' +- ' sprintf('%.6f',err_D_Sigma2)])
disp (['         R^2                =  ' sprintf('%.6f', rsqS2)])
disp(' ')
disp ('       -  Anomalous diffusion')
disp (['         offset (\mum^2)    =  ' sprintf('%.6f', SSc(1)) ' +- ' sprintf('%.6f',par_c(1))])
disp (['         \alpha             =  ' sprintf('%.6f', pow) ' +- ' sprintf('%.6f',err_pow)])
disp (['         R^2                =  ' sprintf('%.6f', rsqS_conf)])
disp(' ')
if pow<1
disp ('       -  Confined diffusion')
disp (['         offset (\mum^2)    =  ' sprintf('%.6f', SSc2(1)) ' +- ' sprintf('%.6f',par_c2(1))])
disp (['         D_M (\mum^2/s)     =  ' sprintf('%.6f', D_M) ' +- ' sprintf('%.6f',err_D_M)])
disp (['         D_m (\mum^2/s)     =  ' sprintf('%.6f', D_m) ' +- ' sprintf('%.6f',err_D_m)])
disp (['         L^2 (\mum^2)       =  ' sprintf('%.6f', Lquadro) ' +- ' sprintf('%.6f',par_c2(3))])
disp (['         \tau_c (s)         =  ' sprintf('%.6f', tauc) ' +- ' sprintf('%.6f',par_c2(4))])
disp (['         R^2                =  ' sprintf('%.6f', rsqS_conf2)])
%disp(' ')
else
disp ('       -  Brownian diffusion + flow motion')
disp (['         offset (\mum^2)    =  ' sprintf('%.6f', SS(3)) ' +- ' sprintf('%.6f',par(3))])
disp (['         D (\mum^2/s)       =  ' sprintf('%.6f', D_Sigma) ' +- ' sprintf('%.6f',err_D_Sigma)])
disp (['         v_\sigma (\mum/s)  =  ' sprintf('%.6f', v_Sigma) ' +- ' sprintf('%.6f',err_v_Sigma)])
disp (['         R^2                =  ' sprintf('%.6f', rsqS)])
disp(' ')
%disp (['         v_\sigma (\mum/s)  =  ' sprintf('%.6f', SS(2)) ' +- ' sprintf('%.6f',par(2))])
disp (['         v_tot    (\mum/s)  =  ' sprintf('%.6f', v_tot_Sigma) ' +- ' sprintf('%.6f',err_v_tot_Sigma)])
%disp (' ')
%disp (['         {a_\phi}^2         =  ' sprintf('%.3f', A_asimm)])
%disp (['         {a_\sigma}^2       =  ' sprintf('%.3f', A_simm)])
%disp ('')
disp (['         \psi (rad)         =  ' sprintf('%.3f', beta)])
%disp('')
end
disp ('     ______________________________________________')
disp('')
disp('                                               DONE')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% g(\rho, \theta, \tau) %%%%%%%%

Cont=STCorr_a(size(STCorr_a,1)/2+1+int64((cont_0*omega+0*v_drift*double(t_lg))/px_size),...
  (size(STCorr_a,1)/2)+1+int64((cont_0*omega+0*v_drift*double(t_lg))/px_size),end);

figure (10)
[c,h]=contour(1:length(a1),1:length(a1),STCorr_a(:,:,1),Cont);
[d,h1]=contour(a1,a1,STCorr_a(:,:,1),Cont);
close (10)
cc=c(:,2:end);
dd_=d(:,2:end);
dd_=dd_';
tht=2*pi/50;
theta_=-pi:tht:pi-tht;
%rho=[dd_(1,1) dd_(1,2)];
%RHO=sqrt(rho(1)^2+rho(2)^2);
RHO=cont_0*omega+fact*v_drift*double(t_lg);
dd=[-RHO*cos(theta_') RHO*sin(theta_')];
DD_= int64([-RHO*cos(theta_') RHO*sin(theta_')]/px_size);
DD=int64([-RHO*cos(theta_') RHO*sin(theta_')]/px_size)+size(STCorr_a,1)/2+1;

theta(1:length(theta_)/2)=-acos(double(DD_(length(theta_)/2+1:end,1))/(RHO/px_size));
theta(length(theta_)/2+1:length(theta_))=acos(double(DD_(1:length(theta_)/2,1))/(RHO/px_size));

theta(1)=theta_(1);
for j=2:length(theta_)
    if theta_(j)>=theta_(j-1)
        theta(j)=theta_(j);
    end
end

cc=int64(cc');
mid=mean(double(cc(:,1)));
dc=double(cc-mid);
tht=2*pi/length(cc);

d_theta=zeros(length(DD),t_lg);
t2=t/f;
tau=t2+1;

for k=1:t_lg
        for i=1:length(theta);
                d_theta(i,k)=STCorr_a(DD(i,2),DD(i,1),int64(tau(k)));
    end
    end

F_theta=d_theta(:,end);
Ftheta=F_theta';
lbFt=[-1 -pi 0 0 -1 -pi 0];
ubFt=[1 pi 10 10 1 pi 10];
startFt=[0.1 ALPHA 0.4 0.2 0.1 -ALPHA 0.3];

options.Robust='off';
options.MaxIter=500;
 ref1=int64(mean(find(int64(theta-ALPHA)==0)));
 Gthetapm=[Ftheta Ftheta Ftheta];
 Gtheta=Gthetapm(length(theta)+ref1-length(theta)/2:length(theta)+ref1+length(theta)/2-1);
 Htheta(1:length(theta)/2)=Gtheta(length(theta)/2+1:end);
 Htheta(length(theta)/2+1:length(theta))=Gtheta(1:length(theta)/2);
 
%%%%%%%%%%%%%%%%%%%%%%%%% 
lGt=[0 -pi 0 -10];
uGt=[10 pi 2*pi 10];
startGt=[0.1 0 0.1 0.1]; 
modelFunGt = @(pG)(pG(1).*((sinc(pG(3).*(theta-pG(2)))).^2)+pG(4)-Gtheta);
[SSGt,resnormGt,residualGt,exitflagGt,outputGt,lambdaGt,jacobianGt]=lsqnonlin(modelFunGt,startGt,lGt,uGt,options);
 dfeGt = length(theta)- 4;
 [QGt,RGt] = qr(jacobianGt,0);
 RGtinv = inv(RGt);
 seGt = sqrt(sum(RGtinv.*RGtinv,2)*resnormGt/dfeGt);
rsqGt=1-(sum(residualGt.^2,2) ./ sum((Gtheta - mean(Gtheta)).^2,2));
fitFunGt=modelFunGt(SSGt)+Gtheta;
%%%%%%%%%%%%%%%%%%%%%%%%%%% 
lHt=[0 -pi 0 -10];
uHt=[10 pi 2*pi 10];
startHt=[0.1 0 0.1 0.1]; 
modelFunHt = @(pH)(pH(1).*(1-(sinc(pH(3).*((theta+pi)-pH(2)))).^2)+pH(4)-Htheta);
[SSHt,resnormHt,residualHt,exitflagHt,outputHt,lambdaHt,jacobianHt]=lsqnonlin(modelFunHt,startHt,lHt,uHt,options);
 dfeHt = length(theta)- 4;
 [QHt,RHt] = qr(jacobianHt,0);
 RHtinv = inv(RHt);
 seHt = sqrt(sum(RHtinv.*RHtinv,2)*resnormHt/dfeHt);
rsqHt=1-(sum(residualHt.^2,2) ./ sum((Htheta - mean(Htheta)).^2,2));
fitFunHt=modelFunHt(SSHt)+Htheta;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lbFt=[-1 -pi 0 0 -1 -pi 0];
ubFt=[1 pi 1 1 1 pi 1];
startFt=[0.1 ALPHA 0.4 0.2 0.1 pi-ALPHA 0.3];

modelFunFt = @(pF)(pF(1).*((sinc(pF(3).*(theta-pF(2)))).^2)+pF(5).*(1-(sinc(pF(7).*(theta-pF(2)))).^2)+pF(4)-Ftheta);
 [SSFt,resnormFt,residualFt,exitflagFt,outputFt,lambdaFt,jacobianFt]=lsqnonlin(modelFunFt,startFt,lbFt,ubFt,options);
 dfeFt = length(theta)- 7;
 [QFt,RFt] = qr(jacobianFt,0);
 RFtinv = inv(RFt);
 seFt = sqrt(sum(RFtinv.*RFtinv,2)*resnormFt/dfeFt);
rsqFt=1-(sum(residualFt.^2,2) ./ sum((Ftheta - mean(Ftheta)).^2,2));
fitFunFt=modelFunFt(SSFt)+Ftheta;
par2=zeros(2,1);
par2(1)=seFt(2);
par2(2)=seFt(3);
phi_Ft=SSFt(2);
psi_Ft=2*SSFt(3)/(pi);
err_phi_Ft=par2(1);
err_psi_Ft=par2(2);
ang_distr=[phi_Ft err_phi_Ft psi_Ft err_psi_Ft rsqFt];

thr1=zeros(length(t))+px_size/(2*sqrt(2));
thr2=zeros(length(t))+px_size/(2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  plot2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure (2)
subplot (2,3,1)
pcolor (dim,dim,A(:,:,1))
axis square
shading interp
title ('image','fontsize',14)
xlabel ('x (\mum)','fontsize',14)
ylabel ('y (\mum)','fontsize',14)
set (gca,'fontsize',14)

subplot (2,3,2)
pcolor (-a1,a1,STCorr_a(:,:,end))
shading interp
colormap jet
hold on
plot (x_G,y_G,'-b',x_G(1),y_G(1),'black+',x_G(end),y_G(end),'blacko','linewidth',1,'markersize',8)
axis square
plot (dd(:,1),dd(:,2),'--black','linewidth',1.1)
%plot (cc(:,1),cc(:,2))
title ('g -> top view','fontsize',14)
xlabel ('\xi (\mum)','fontsize',14)
ylabel ('\eta (\mum)','fontsize',14)
 set (gca,'fontsize',14)
hold off

subplot (2,3,3)
p3=zeros(length(a1),1)+max(max(a1))+px_size/2;
p4=zeros(length(a1),1)-min(min(a1))-px_size/2;
 surf (-a1-0*px_size/2,a1,STCorr_a(:,:,end));
 hold on
  plot3 (a1,p4,Gauss_x,'black','linewidth',0.8)
  plot3 (p3,a1,Gauss_y,'black','linewidth',0.8);
 hold off
 xlim([min(min(a1)) max(max(a1))+px_size/2])
 ylim([min(min(a1)) max(max(a1))])
 xlabel '\xi'
 ylabel '\eta'
 zlim ([min(min(min(STCorr_a))) max(max(max(STCorr_a)))])
 grid on
 axis square
 zlabel ('g(\xi,\eta,\tau)','fontsize',14)
 xlabel ('\xi (\mum)','fontsize',14)
 ylabel ('\eta(\mum)','fontsize',14)
 title ('g -> 3d view','fontsize',14)
 set (gca,'fontsize',14)
 box on
 
 subplot (2,3,4)
plot (t,x_G,'blackx',t,y_G,'black+','linewidth',1.5)
legend ('\xi-shift','\eta-shift','location','southwest')
legend ('boxon')
hold on
plot (t,thr2,'-.r',t,-thr2,'-.r','linewidth',1.3)
plot (t,thr1,':r',t,-thr1,':r','linewidth',1.1)
hold off
xlabel ('\tau (s)','fontsize',14)
ylabel ('peak position (\mum)','fontsize',14)
xlim ([0-f t_2(end)+f])
set (gca,'fontsize',14)
title ('net flow','fontsize',14)
axis square

subplot (2,3,5)
plot(t_2,fitFunS2,'--r','linewidth',1.5)
hold on
%plot(t_2,fitFunS,'--b','linewidth',1.5)
plot(t_2,fitFun_conf,'-c','linewidth',1.5)
if pow<1 plot(t_2,fitFun_conf2,'--g','linewidth',1.5)
else plot (t_2,fitFunS,'--b','linewidth',1.5)
end
plot(t_2,Sigma2,'blacko','linewidth',1.5)
hold off
if pow<1 legend ('Brownian','Anomalous','Confined','location','southeast')
else  legend ('Brownian','Anomalous','Flow','location','southeast')
end
legend ('boxon')
xlabel ('\tau (s)','fontsize',14)
ylabel ('\sigma^2 (\mum^2)','fontsize',14)
xlim ([0 t_2(end)+f])
set (gca,'fontsize',14)
title ('iMSD curve','fontsize',14)
axis square
 
subplot (2,3,6)
h=surf(t,theta,d_theta);
hold on
view([-37.5+90 25])
shading interp
xlabel ('\tau (s)','FontSize',14)
ylabel ('\theta (rad)','FontSize',14)
zlabel ('g(\rho_\infty, \theta, \tau)','fontsize',14)
grid on
colormap jet
hold off
ylim ([-pi pi])
xlim ([0 t_2(end)+f])
set (gca,'FontSize',14)
set(h','edgecolor','k')
axis square
title ('g -> circular section','fontsize',14)
box on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
figure (3)
subplot (1,3,1)
plot(t_2,Sigma2,'blacko','linewidth',1.5)
hold on
plot(t_2,fitFunS2,'r','linewidth',1.5)
title ('Brownian Diffusion','FontSize',14)
legend ('iMSD data',sprintf('R^2=%.5f',rsqS2),'location','northwest')
legend ('boxoff')
hold off
xlabel ('\tau (s)','fontsize',14)
ylabel ('\sigma^2 (\mum^2)','fontsize',14)
xlim ([0 t_2(end)+f])
set (gca,'fontsize',14)
axis square

subplot (1,3,2)
plot(t_2,Sigma2,'blacko','linewidth',1.5)
hold on
plot(t_2,fitFun_conf,'c','linewidth',1.5)
legend ('iMSD data',sprintf('R^2=%.5f',rsqS_conf),'location','northwest')
legend ('boxoff')
title (sprintf('Anomalous Diffusion \\alpha = %.3f',pow),'FontSize',14)
hold off
xlabel ('\tau (s)','fontsize',14)
ylabel ('\sigma^2 (\mum^2)','fontsize',14)
xlim ([0 t_2(end)+f])
set (gca,'fontsize',14)
axis square

if pow<1
subplot (1,3,3)
plot(t_2,Sigma2,'blacko','linewidth',1.5)
hold on
plot(t_2,fitFun_conf2,'g','linewidth',1.5)
title ('Confined Diffusion','FontSize',14)
legend ('iMSD data',sprintf('R^2=%.5f',rsqS_conf2),'location','northwest')
legend ('boxoff')
hold off
xlabel ('\tau (s)','fontsize',14)
ylabel ('\sigma^2 (\mum^2)','fontsize',14)
xlim ([0 t_2(end)+f])
set (gca,'fontsize',14)
axis square
else
subplot (1,3,3)
plot(t_2,Sigma2,'blacko','linewidth',1.5)
hold on
plot(t_2,fitFunS,'b','linewidth',1.5)
legend ('iMSD data',sprintf('R^2=%.5f',rsqS),'location','northwest')
legend ('boxoff')
title ('Brownian Diffusion + Flow','FontSize',14)
hold off
xlabel ('\tau (s)','fontsize',14)
ylabel ('\sigma^2 (\mum^2)','fontsize',14)
xlim ([0 t_2(end)+f])
set (gca,'fontsize',14)
axis square
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% END %%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  export %%%%
pre=[name_start '.xlsx'];
delete (pre);

head1=cellstr('t');
head2=cellstr('\xi_0');
head3=cellstr('\eta_0');
head4=cellstr('Sigma^2');
head=[head1 head2 head3 head4];
filename=[name_start '.xlsx'];
xlswrite(filename,head,'Foglio1')
xlswrite(filename,[t_2 x_G(2:end) y_G(2:end) Sigma_1],'Foglio1','A2')

txt1=cellstr('INPUT');
txt2=cellstr('sqrt(image size) (px)');
txt3=cellstr('pixel size (\mum)');
txt4=cellstr('sqrt(fitting-ROI size) (px)');
txt5=cellstr('time lapse (s)');
txt6=cellstr('max lag-time (s)');
txt7=cellstr('number of frames');

text1=cellstr('omega (\mum)');
text2=cellstr('D Brownian (\mum^2/s)');
text3=cellstr('offset (\mum^2)');
text4=cellstr('R^2');
text5=cellstr('\alpha Anomalous');
text6=cellstr('offset (\mum)^2');
text7=cellstr('R^2');
text8=cellstr('D_M (\mum^2/s)');
text9=cellstr('D_m (\mum^2/s)');
text10=cellstr('L^2 (\mum^2)');
text11=cellstr('\tau_c (s)');
text12=cellstr('offset (\mum^2)');
text13=cellstr('R^2');
text14=cellstr('D Br.+Flow (\mum^2/s)');
text15=cellstr('v_\phi (\mum/s)');
text16=cellstr('v_\sigma (\mum/s)');
text17=cellstr('offset (\mum^2)');
text18=cellstr('R^2');

txt0=[txt1;txt2;txt3;txt4;txt5;txt6;txt7;cellstr('bkg-corr')];
nm0=[P;px_size;2*P/cutoff;f;t_2(end);double(N);Filter];
text0=[text1;text15;cellstr('\phi');text2;text3;text4;text5;text6;text7;text8;text9;text10;text11;text12;text13;text14;text16;text17;text18];
head2= [cellstr('OUTPUT') cellstr('mean') cellstr('error')];
nm2=[omega^2;v_drift1;ALPHA;D_Sigma2;SS2(2);rsqS2;pow;SSc(1);rsqS_conf;D_M;D_m;Lquadro;tauc;SSc2(1);rsqS_conf2;D_Sigma;v_Sigma;SS(3);rsqS];
err2=[2*omega*err_omega;err_v_drift1;NaN; err_D_Sigma2;par2(2);NaN;err_pow;par_c(1);NaN;err_D_M;err_D_m;par_c2(3);par_c2(4);par_c2(1);NaN;err_D_Sigma;err_v_Sigma;par(3);NaN];

xlswrite(filename,txt0,'Foglio2','A1')
xlswrite(filename,nm0,'Foglio2','B2')
xlswrite(filename,head2,'Foglio2','A10')
xlswrite(filename, text0,'Foglio2','A11')
xlswrite(filename, [nm2 err2],'Foglio2','B11')
xlswrite(filename, cellstr(sprintf('x%d',P)),'Foglio2','C2')
xlswrite(filename, cellstr(sprintf('x%d',2*P/cutoff)),'Foglio2','C4')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% animation
return
figure

for i=1:t_lg
subplot (1,2,2)
pcolor (-a1, b1,STCorr_a(:,:,i))
shading interp;
axis square
T=sprintf('%d',i-1);
title (['r(\xi,\eta, ' T ')'],'FontSize',12)
%set (xticks, ([min(a1):max(a1)]))
xlabel ('\xi (\mum)','FontSize',12)
ylabel ('\eta (\mum)','FontSize',12')
set(gca,'FontSize',12)

subplot (1,2,1)
surf (-a1, b1,STCorr_a(:,:,i))
%shading interp;
axis square
T=sprintf('%d',i-1);
zlabel (['r(\xi,\eta, ' T ')'],'FontSize',12)
zlim ([0 max(max(max(STCorr)))])
xlabel ('\xi (\mum)','FontSize',12)
ylabel ('\eta (\mum)','FontSize',12)
%xticks (min(a1):max(a1))
xlim ([min(a1) max(a1)])
ylim ([min(a1) max(a1)])
set(gca,'FontSize',12)
box on

pause(0.5)
%saveas(gcf,sprintf('export_AUTOC%d.png',i-1))
colormap jet
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
subplot (2,3,6)
plot (theta,F_theta,':blacko','linewidth',1)
xlabel ('\theta (rad)','FontSize',12)
ylabel ('r(\rho_\infty, \theta, \tau_{M})','FontSize',12)
hold on
plot (theta,fitFunFt,'black','linewidth',1.5')
hold off
set (gca,'FontSize',12)
xlim ([-pi-0.05 pi+0.05])