clear
clc
%load Try_data.txt;
%I = Try_data(:,2);
%En  = Try_data(:,3);

load En_vs_I_26Apr2012.txt;
I = En_vs_I_26Apr2012(:,2);
En  = En_vs_I_26Apr2012(:,3);
En_er = En*0+0.01;
I_er = I*0 + 0.1;
%En_er = En*0.01;
%I_er = I*0.1 + 0.1;

[NumPoints, row] =  size(En);

a(1) = max(I);
a(2) = mean(En);
a(3) = std(En);
a(4) = std(En);

%-------------------------- Fitting --------------------------
opts = optimset('TolX',1e-4,'MaxFunEvals',10000,'MaxIter',10000,'Display','on');
[par,fval,flag,err] = fminsearch(@Assym_Gau_devsum,a,opts,En,I);
%-------------------------------------------------------------
sigma = (par(3)+ par(4))/2;
skew = ( par(4) -par(3) )/( par(3)+par(4) );
%plot(En,I), hold on;

Ene=4:0.01:15;

[row,col] = size(Ene);
%for i=1:row
%signe(i) = sign(Ene(i)-par(2));
%end
signe = sign(Ene-par(2));
%signe
size(signe)

set(axes,'FontSize',16)
errorbar(En,I,I_er,'MarkerSize',30,'Marker','.','LineStyle','none',...
    'LineWidth',3,...
    'DisplayName','Data Points',...
    'Color',[0 0 1]); 
 hold on

 for i=1:col
 Gau_fit(i) = par(1)*exp( -0.5*((Ene(i)-par(2))/(sigma*(1+signe(i)*skew))).^2 );
 end
 %Gau_fit'
 %Gau_fit = par(1)*exp( -(Ene-par(2)).^2 / (2*((par(3)+par(4))/2*(1+ (par(3)-par(4))/(par(3)+par(4)))).^2) );% Assymetric Gaussian
 %Gau = a(1)*exp( -(X-a(2)).^2 / (2*((a(3)+a(4))/2*(1+ (a(3)-a(4))/(a(3)+a(4)))).^2) );% Assymetric Gaussian
 %Gau_fit = par(3)*exp( -(Ene-par(2)).^2 / (2*par(1).^2) );% Gaussian
 %Gau_fit = 1/(sqrt(2*pi)*par(1))*exp( -(Ene-par(2)).^2 / (2*par(1).^2) );% Gaussian
 %Gau_fit = par(3) + 1/(sqrt(2*pi)*par(1))*exp( -(Ene-par(2)).^2 / (2*par(1).^2) );% Gaussian
 plot(Ene,Gau_fit,'r','LineWidth',2); %plot fit done by matrix inversion method  

%set(axes,'FontSize',40)
% pardat = errorbar(En,Gau_fit,En_er,'MarkerSize',30,'Marker','.','LineStyle','none',...
%     'LineWidth',3,...
%     'DisplayName','Data Points',...
%     'Color',[0 0 1]); 

 legend('Data', 'Fit');
 %set(legend,'Position',[0.34 0.68 0.1 0.1]);
 set(legend,'FontSize',20);
 set(legend,'FontName','Times New Roman');

 hy = ylabel('I_{FC} (mA)','FontSize',24,'FontName','Times New Roman');
 hx = xlabel('En (MeV)','FontSize',24,'FontName','Times New Roman');
% 
fit=title({'Electron Beam Energy Spread';'at 12 MeV'});
set(fit,'FontSize',24,'FontWeight','bold','FontName','Times New Roman');
% axis([0 0.55 2.6 3.65]); %axis([xmin xmax ymin ymax])
% 
% delete('Gaussian_fit.png')
% saveas(fit,'Gaussian_fit.pdf','pdf');
 saveas(fit,'Gaussian_fit.eps','epsc');
 saveas(fit,'Gaussian_fit.png','png');
 drawnow


 
fprintf(' Amplitude = %.5f, mean = %.5f, sigma_L = %.5f, sigma_R = %.5f \n\n',  par(1), par(2), par(3), par(4));
fprintf(' Sigma = %.5f, Skewness = %.5f \n\n',  sigma, skew);

%fprintf(' Gaussian Fit for Electron Beam Energy: \n En = %.5fexp(-(En-%.5f)^2) / (2*%.5f^2) \n\n', par(3), par(2),par(1));
%fprintf(' Gaussian Energy peak = %.5f \n\n',  par(2));
%fprintf(' Energy spread = %.5f \n rms Energy spread = %.5f \n\n', par(1)*6, par(1));
%fprintf(' Relative Energy spread = %.5f \n Relative rms Energy spread = %.5f \n\n', par(1)*6/par(2),par(1)/par(2));


%fprintf(' mean = %.5f, sigma = %.5f \n\n',  par(2), par(1));
%fprintf('Gaussian Fit for Electron Beam Energy: \n En = %.5fexp(-(En-%.5f)^2) / (2*%.5f^2) \n\n', 1/(sqrt(2*pi)*par(1)), par(2),par(1));

%fprintf('off-set = %.5f, mean = %.5f, sigma = %.5f \n\n', par(3), par(2), par(1));
%Gau_fit = par(3)+1/(sqrt(2*pi)*par(1))*exp( -(En-par(2)).^2/(2*par(1).^2) );% Gaussian
%fprintf('Gaussian Fit for Electron Beam Energy: \n En = %.5f + %.5fexp(-(En-%.5f)^2) / (2*%.5f^2) \n\n', par(3), 1/(sqrt(2*pi)*par(1)), par(2),par(1));