%% Mie Theory Fitting for Silver Nanoparticles%%
%% Reference for Optical Constants, n and k, epslion1=n^2-k^2, epslion2=2*n*k: 1. Palik  E. D., Handbook of Optical Constants., Academic Press Inc., San Diego (1985); 2. Johnson  P., Christy  R. W., Optical constants of the noble metals, Phys. Rev. B. 6, (12), 4370 - 4379 (1972)%%
%% This MATLAB code is generated from MATLAB R2013a.%%
%%

clc
clear all


%Step 1: raw data modification%
fprintf('Select the test.xlsx file!\n\n');
FileName = uigetfile('*.xlsx'); %select the 'test.xlsx' file%
num = xlsread(FileName);
xx=num(:,1);
xx=roundn(xx,0);
num1(:,1)=xx;
num1(:,2)=num(:,2);
xlswrite('textdata.xlsx',num1)
fprintf('Step 1: Data modification done!\n\n');

%Step 2: experimental data import%
pause; %press any key to strat%
fprintf('Importing data strats! Select the textdata.xlsx file!\n\n');
FileName = uigetfile('*.xlsx'); %select the 'textdata.xlsx' file%
data=xlsread(FileName);
WL=data(:,1);
AB=data(:,2);
NormalAB=AB/max(AB);
wlIndex = find(AB == max(AB), 1, 'first');
maxwlValue = WL(wlIndex);
fprintf('Maximun extintion is located at %d nm.\n\n',maxwlValue);
Normaldata(:,1)=WL;
Normaldata(:,2)=NormalAB;
xlswrite('normaltextdata.xlsx',Normaldata);
fprintf('Step 2: experimental data import done!\n\n');

%Step 3: generate the constant file for the wavelength used%
pause;%press any key to strat%
fprintf('Generating the consatant file starts! Select the epslion1,epslion2 for Ag.xlsx file!\n\n');
FileName=uigetfile('*.xlsx');
num2=xlsread(FileName);
x2=Normaldata(:,1);
x1=num2(:,1);
y1=num2(:,2);
z1=num2(:,3);
y2=interp1(x1,y1,x2);
z2=interp1(x1,z1,x2);
JCConst=[x2,y2,z2];
xlswrite('J&CAg.xlsx',JCConst);
fprintf('Step 3: constant file generation done!\n\n');

%Step 4: Mie theory fitting for nanoparticle size%
pause;%press any key to strat%
DIConst=1/(3.1*10^-14);%the bulk metal value for Ag, s^-1%
vF=1.39*10^6;%the Fermi speed of Ad, m*s^-1%
ne=5.86*10^28;%Free electron number density for Ag,atom/m^3%
mEFF=0.99*9.1093897*10^-31;%The effective mass of electrons in Ag, kilogram%
e=1.60217733*10^-19;%elementary charge of proton, C%
e0=8.854*10^-12;%Permittivity of free space, F/m%
c=299792458;%Speed of light, m/s%
wp=(ne*e*e/e0/mEFF)^0.5;

nmedium=input('refraction index for the solvent= \n');%
%nmedium=1.34;
coeffA1=input('coeffA1= \n');%modification of parameter A%
%coeffA1=-0.5495;
coeffA2=input('coeffA2= \n');
%coeffA2=0.6905;
coeffA3=input('coeffA3= \n');
%coeffA3=0;
coeffA4=input('coeffA4= \n');
%coeffA4=0;
Dini=input('the diameter of the fitting starts from (nm) \n');
%Dini=1;
Dfinal=input('the diameter of the fitting ends with (nm) \n');
%Dfinal=8;
Dstep=input('the diameter interval of the fitting is (nm) \n');
%Dstep=1;
Lmax=3;%set multipolar oder of the Mie calculation%

target_D=0;
delta1=+inf;
target_spec=[];


syms V1 V2 L
XL1=V1*(pi/2/V1)^0.5*besselj(L+0.5,V1);
XL2=V2*(pi/2/V2)^0.5*besselj(L+0.5,V2);
YL1=V1*(pi/2/V1)^0.5*(besselj(L+0.5,V1)+i*bessely(L+0.5,V1));
XL11=diff(V1*(pi/2/V1)^0.5*besselj(L+0.5,V1),V1);
XL22=diff(V2*(pi/2/V2)^0.5*besselj(L+0.5,V2),V2);
YL11=diff(V1*(pi/2/V1)^0.5*(besselj(L+0.5,V1)+i*bessely(L+0.5,V1)),V1);


ATAT=[];

for D=Dini:Dstep:Dfinal
    delta=0;
    A=coeffA1+coeffA2*(D/2)+coeffA3*(D/2)^2+coeffA4*(D/2)^3;
    DI=DIConst+2*A*vF/D*10^9;
    calresult=[];
    for n=1:1:length(x2)
        lambda=x2(n);
        k=2*pi*10^9/lambda;
        w=c*k;
        enne=(y2(n)+wp^2*(1/(w^2+DIConst^2)-1/(w^2+DI^2))+i*(z2(n)+wp^2/w*(DI/(w^2+DI^2)-DIConst/(w^2+DIConst^2))))^0.5;
        x=k*D/2*enne*10^-9;
        m=nmedium/enne;
        extcal=0;
        for L=1:Lmax%Lmax=3%

        V1=m*x;
        V2=x;

        XL1_0=eval(XL1);
        XL2_0=eval(XL2);
        YL1_0=eval(YL1);
        XL11_0=eval(XL11);
        XL22_0=eval(XL22);
        YL11_0=eval(YL11);
        aL=(-(m*XL1_0*XL22_0-XL11_0*XL2_0)/(m*YL1_0*XL22_0-YL11_0*XL2_0));
        bL=(-(XL1_0*XL22_0-m*XL11_0*XL2_0)/(YL1_0*XL22_0-m*YL11_0*XL2_0));
        extcal=extcal+(-2*pi)/(nmedium*k)^2*(2*L+1)*real(aL+bL);
        end
        calresult(n,:)=[x2(n),extcal];
        
        
    end
    ATAT=[ATAT calresult];
    calresult_nor(:,1)=calresult(:,1);
    calresult_nor(:,2)=calresult(:,2)/max(calresult(:,2));
    NormalAB_cal=calresult_nor(:,2);
    WL_cal=calresult_nor(:,1);
    wlIndex_cal = find(NormalAB_cal == max(NormalAB_cal), 1, 'first');
    maxwlValue_cal = WL_cal(wlIndex_cal);
    shift=maxwlValue_cal-maxwlValue;
    if shift~=0
       calresult_nor(:,1)=calresult_nor(:,1)-shift;
    end
    step=shift/(x2(1)-x2(2));
    nini=1;
    nfinal=length(x2);
    if step<0
        nini=1;
        nfinal=length(x2)+step;
    elseif step>0
                nini=1+step;
                nfinal=length(x2);
    end
    for n1=nini:nfinal
        delta=delta+(NormalAB(n1)-NormalAB_cal(n1-step))^2;
    end
    if delta<delta1
        delta1=delta;
        target_D=D;
        target_spec=calresult_nor;
        else delta1=delta1;
        target_D=target_D;
        target_spec=target_spec;
    end
    fprintf('D = %3.2f nm.\n',D);
    fprintf('Delta is %6f.\n',delta);
    fprintf('Best fitting is %3.2f nm.\n\n',target_D);
    
    
end
fprintf('The best fitting of Uv-vis spectra is %3.2f nm.\n',target_D);
fprintf('The best calculated spectra has been shifted of %d nm.\n',shift);
fprintf('The delta of best fitting is %6f.\n\n',delta1);
figure
scatter(target_spec(:,1),target_spec(:,2),'r')
hold on
scatter(Normaldata(:,1),Normaldata(:,2),'b')
hold off
legend('calculated','experimental')
axis([300 800 0 1])
fprintf('Step 4: Mie fitting done!\n\n');

%Step 5: Mie theory fitting with particle size distribution for nanoparticle size%
pause;
target_Xc=0;
target_SD=0;
delta1=+inf;
target_spec=[];

for Xc=target_D/2-2*Dstep:Dstep/10:target_D/2%the setp can be changed to get different resolution%
    for SD=0.05:0.01:1
        BTBT=[];
        LogNorm_sum=0;
        delta=0;
        for g=Dini/2:Dstep/2:Dfinal/2
           
            LogNorm=exp(-log(g/Xc)^2/(2*SD^2))/((2*pi)^0.5*SD*g);%possibility distribution of each Diameter
            LogNorm_sum=LogNorm_sum+LogNorm;%the sum of possibility
            BTBT=[BTBT LogNorm];            
          
            
        end
        BTBT=BTBT/LogNorm_sum;
        CTCT=0;
        for nn=2:2:size(ATAT,2)
            CTCT=CTCT+BTBT(:,nn/2)'*ATAT(:,nn);
        end
        CTCT=CTCT/max(CTCT);
        DTDT=[ATAT(:,1) CTCT];%below the fitting starts%
                
        NormalAB_cal=DTDT(:,2);
    WL_cal=DTDT(:,1);
    wlIndex_cal = find(NormalAB_cal == max(NormalAB_cal), 1, 'first');
    maxwlValue_cal = WL_cal(wlIndex_cal);
    shift=maxwlValue_cal-maxwlValue;
    if shift~=0
       DTDT(:,1)=DTDT(:,1)-shift;
    end
    step=shift/(x2(1)-x2(2));
    nini=1;
    nfinal=length(x2);
    if step<0
        nini=1;
        nfinal=length(x2)+step;
    elseif step>0
                nini=1+step;
                nfinal=length(x2);
    end
    for n1=nini:nfinal
        delta=delta+(NormalAB(n1)-NormalAB_cal(n1-step))^2;
    end
    if delta<delta1
        delta1=delta;
        target_Xc=Xc;
        target_spec=DTDT;
        target_SD=SD;
        else delta1=delta1;
        target_Xc=target_Xc;
        target_spec=target_spec;
        target_SD=target_SD;
    end
    end
end
LogNorm_sum_target=0;
sum=0;
for g=Dini/2:Dstep/2:Dfinal/2
    LogNorm_target=exp(-log(g/target_Xc)^2/(2*target_SD^2))/((2*pi)^0.5*target_SD*g);
    sum=sum+g*LogNorm_target;
    LogNorm_sum_target=LogNorm_sum_target+LogNorm_target;
end
averageR=sum/LogNorm_sum_target;
LogNorm_sum_target=0;
sum=0;
for g=Dini/2:Dstep/2:Dfinal/2
    LogNorm_target=exp(-log(g/target_Xc)^2/(2*target_SD^2))/((2*pi)^0.5*target_SD*g);
    sum=sum+(g-averageR)^2*LogNorm_target;
    LogNorm_sum_target=LogNorm_sum_target+LogNorm_target;    
end
standDEV=(sum/LogNorm_sum_target)^0.5;
fprintf('Best Xc is %3.2f nm.\n',2*target_Xc);
fprintf('Best SD is %3.2f.\n',target_SD);
fprintf('Delta is %6f.\n\n',delta1);
fprintf('Average diameter is %3.2f nm.\n',2*averageR);
fprintf('Standard deviation is %3.2f nm.\n',2*standDEV);
figure
scatter(target_spec(:,1),target_spec(:,2),'r')
hold on
scatter(Normaldata(:,1),Normaldata(:,2),'b')
hold off
legend('calculated','experimental')
axis([300 800 0 1])
fprintf('Step 5: Mie fitting with size distribution done!\n\n');