clc;
clear

me = 0.511; % electron rest mass in MeV
%En = 14; %enter beam energy here in MeV unit.
En = 12.74; %enter beam energy here in MeV unit.
%I_max_neg = -10; %enter smallest negative scan current. Here -7 Amps
%I_max_pos = 7; %enter biggest positive scan current. Here 10 Amps
I_start=-8; % Scan starting current
%I_end=-2.0; % Scan ending current
I_increment=0.5;
scan_step=13; % do 41 scan.


for ii = 1:scan_step
    I(ii) = I_start + (ii-1)*I_increment;
end


p = sqrt(En*En - me*me)/1000;% e- momentum in GeV/c

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%enter the distance from the center of the quad to the screen
S12=3.1;%(in m unit) S12(q4)=1.875 cm. S12(q1)=3.10 cm
erS12 = 0.005; % assume 51 mm error.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%enter calibration here:
calibration_x=0.04308*0.001; %m/px
calibration_y=0.04228*0.001; %m/px
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% now we do the action......................................
disp('Read image ...')
% get image
%rehash path
rmpath('C:\Documents and Settings\setisadi\My Documents\MATLAB\analysis\Mar_18_2011\Q1_Scan\42mA\positive\2_Amp');
%addpath(genpath('C:\Documents and Settings\setisadi\My Documents\MATLAB\analysis\Mar_18_2011\Q1_Scan\42mA\negatvie\2.4'));
addpath('C:\Documents and Settings\setisadi\My Documents\MATLAB\analysis\Mar_18_2011\Q1_Scan\42mA\positive\1.8_Amp');
%cd('C:\Documents and Settings\setisadi\My Documents\MATLAB\analysis\Mar_18_2011\Q1_Scan\42mA\positive\2_Amp');


bg1 = imread('RF_on_Gun_off_for_bg_1.bmp');
bg2 = imread('RF_on_Gun_off_for_bg_2.bmp');
bg3 = imread('RF_on_Gun_off_for_bg_3.bmp');
bg4 = imread('RF_on_Gun_off_for_bg_4.bmp');
bg5 = imread('RF_on_Gun_off_for_bg_5.bmp');
bg6 = imread('RF_on_Gun_off_for_bg_6.bmp');
bg7 = imread('RF_on_Gun_off_for_bg_7.bmp');
bg8 = imread('RF_on_Gun_off_for_bg_8.bmp');
bg9 = imread('RF_on_Gun_off_for_bg_9.bmp');
bg10 = imread('RF_on_Gun_off_for_bg_10.bmp');
bg = imread('RF_on_Gun_off_for_bg_1.bmp');
%bg=(bg1+bg2+bg3+bg4+bg5+bg6+bg7+bg8+bg9+bg10)/10;
imd_bg = double(bg(:,:,1));

%im = imread('C:\Documents and Settings\setisadi\My Documents\MATLAB\analysis\Mar_18_2011\Q1_Scan\42mA\Matlab_Files\1.jpg'); % reading 1.png
sg = imread('2.bmp'); % reading 1.png

imd_sg = double(sg(:,:,1));
pic = (imd_sg-imd_bg);
%imd = imd_sg;
%imd = imd_bg;

% determine image size in pixels
dim=size(pic);
%dim=size(imd);
dim_x = dim(1);
dim_y = dim(2);


xunit='m';
yunit='m';

% profiles
profx = sum(pic,1);%how to make projection
profy = sum(pic,2)';%how to make projection

xbins=size(profx);
ybins=size(profy);
x=(1:1:xbins(2));
y=(1:1:ybins(2));

% gauss fit with narrowed range
%xl = 150;
%xr = 550;
%yl = 100;
%yr = 480;

xl = 1;
xr = 766;
yl = 1;
yr = 576;

%xl = 100;
%xr = 650;
%yl = 1;
%yr = 576;

print = 1;
%gau_x = gaussfit(x(xl:xr),profx(xl:xr),print);%gaussion fitting in x
%gau_y = gaussfit(y(yl:yr),profy(yl:yr),print);%gaussion fitting in y

% h is the threshold which is the fraction from the maximum y height that the data is been taken from.
%h = 0.2;

X = x;
Y = profx;
%X = y;
%Y = profy;
plot(X',Y','b')
hold on;

% a(1) = base
% a(2) = A; amplitude
% a(3) = x0; center
% a(4) = D; sigma_0
% a(5) = N; 

a1 = 100;
a2 = 7000;
a3 = 280;
a4 = 1.25;
a5 = 1.5;

%SupGau = a(1)+a(2)*exp(-abs(X-a(3))/(a(4))).^a(5); % super gaussian 

a0 = [a1,a2,a3,a4,a5];

opts = optimset('TolX',1e-4,'MaxFunEvals',10000,'MaxIter',10000,'Display','on');
[fitpara,fval,flag,err] = fminsearch(@SupGau_devsum,a0,opts,X,Y);

base = fitpara(1)
amplitude = fitpara(2)
peak_center = fitpara(3)
sig0 = fitpara(4)
N = fitpara(5)
sig = sig0*(pi/2).^(2/fitpara(5)-1)


%y_fit = 1/(sqrt(2*pi)*fitpara(1))*exp(-abs(X-fitpara(3)).^fitpara(2)/(2*fitpara(1).^fitpara(2)));
y_fit = fitpara(1)+fitpara(2)*exp( -0.5*(abs(X-fitpara(3))/(fitpara(4))).^fitpara(5) );

plot(X',y_fit','r');

%The function fminsearch can be tuned by additional options, one can demand a certain accuracy, restrict the number of optimization steps etc. 
%For a detailed description consult Help