% SuperGaussian_Fit.m
% read one image from files.
% fits a super Gaussian funtion on the projection, get the rms values.
% then pass them to fit parabola
%
% April-2011
% Sadiq Setiniyaz

clc;
clear

me = 0.511; % electron rest mass in MeV
%En = 14; %enter beam energy here in MeV unit. 14 MeV energy beam used at Mar 14, 15, 16 and 17th, 2011.
En = 14; %12.74 MeV energy at Mar 18th, 2011.
L = 0.08; % pole face is 8 cm long in z direction.
R_Bore = 0.0254; %Radius of the Bore aperture = center to pole face = 1 inch.
%L_ef = L + R_Bore;
%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=-3.0; % Scan starting current
%I_end=-2.0; % Scan ending current
I_increment=0.2;
scan_current_number=31; % do 13 scan.
scan_times=7;

%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.04327*0.001; %m/px
calibration_y=0.04204*0.001; %m/px
er_calibration_x=0.00016*0.001;
er_calibration_y=0.00018*0.001;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%enter calibration here:
% calibration_x=0.04327; %mm/px
% calibration_y=0.04204; %mm/px
% er_calibration_x=0.00016;
% er_calibration_y=0.00018;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% now we do the action......................................
disp('Read image ...')
% get image

for ii=1:scan_times
    switch ii

        case 1
            addpath(genpath('C:\Documents and Settings\setisadi\My Documents\MATLAB\analysis\Mar_2011\Mar_17_2011\Q1_Scan\Q1_Scan_by_Sadiq_at_40mA_Peak_Current\Scan_1'));
        case 2
            addpath(genpath('C:\Documents and Settings\setisadi\My Documents\MATLAB\analysis\Mar_2011\Mar_17_2011\Q1_Scan\Q1_Scan_by_Sadiq_at_40mA_Peak_Current\Scan_2'));
        case 3
            addpath(genpath('C:\Documents and Settings\setisadi\My Documents\MATLAB\analysis\Mar_2011\Mar_17_2011\Q1_Scan\Q1_Scan_by_Sadiq_at_40mA_Peak_Current\Scan_3'));
        case 4
            addpath(genpath('C:\Documents and Settings\setisadi\My Documents\MATLAB\analysis\Mar_2011\Mar_17_2011\Q1_Scan\Q1_Scan_by_Sadiq_at_40mA_Peak_Current\Scan_4'));
        case 5
            addpath(genpath('C:\Documents and Settings\setisadi\My Documents\MATLAB\analysis\Mar_2011\Mar_17_2011\Q1_Scan\Q1_Scan_by_Sadiq_at_40mA_Peak_Current\Scan_5'));
        case 6
            addpath(genpath('C:\Documents and Settings\setisadi\My Documents\MATLAB\analysis\Mar_2011\Mar_17_2011\Q1_Scan\Q1_Scan_by_Sadiq_at_40mA_Peak_Current\Scan_6'));
        case 7
            addpath(genpath('C:\Documents and Settings\setisadi\My Documents\MATLAB\analysis\Mar_2011\Mar_17_2011\Q1_Scan\Q1_Scan_by_Sadiq_at_40mA_Peak_Current\Scan_7'));
    end

    bg = imread('RF_on_Gun_off_bg.bmp');
    imd_bg = double(bg(:,:,1));

    for scan = 1:scan_current_number;
        thisfilename=ii;
        imNo = ii;
        I(scan) = I_start + (scan-1)*I_increment;
        switch scan
            case 1

%Devide the effetive length by 100 to convert to meter unit


                im = imread('negative\11.bmp');
                tifname='11';
                L_ef(scan) = 9.95/100; %#ok<AGROW>
                Er_L_ef(scan) = 0.17/100; %#ok<AGROW>
            case 2
                im = imread('negative\12.bmp'); % reading 1.png
                tifname='12';
                L_ef(scan) = 9.91/100;
                Er_L_ef(scan) = 0.16/100;
            case 3
                im = imread('negative\13.bmp');
                tifname='13';
                L_ef(scan) = 9.91/100;
                Er_L_ef(scan) = 0.16/100;
            case 4
                im = imread('negative\14.bmp');
                tifname='14';
                L_ef(scan) = 9.85/100;
                Er_L_ef(scan) = 0.14/100;
            case 5
                im = imread('negative\15.bmp');
                tifname='15';
                L_ef(scan) = 9.85/100;
                Er_L_ef(scan) = 0.16/100;
            case 6
                im = imread('negative\16.bmp');
                tifname='16';
                L_ef(scan) = 9.87/100;
                Er_L_ef(scan) = 0.19/100;
            case 7
                im = imread('negative\17.bmp');
                tifname='17';
                L_ef(scan) = 9.82/100;
                Er_L_ef(scan) = 0.19/100;
            case 8
                im = imread('negative\18.bmp');
                tifname='18';
                L_ef(scan) = 9.86/100;
                Er_L_ef(scan) = 0.17/100;
            case 9
                im = imread('negative\19.bmp');
                tifname='19';
                L_ef(scan) = 9.73/100;
                Er_L_ef(scan) = 0.18/100;
            case 10
                im = imread('negative\20.bmp');
                tifname='20';
                L_ef(scan) = 9.65/100;
                Er_L_ef(scan) = 0.22/100;
            case 11
                im = imread('negative\21.bmp');
                tifname='21';
                L_ef(scan) = 9.71/100;
                Er_L_ef(scan) = 0.22/100;
            case 12
                im = imread('negative\22.bmp');
                tifname='22';
                L_ef(scan) = 9.58/100;
                Er_L_ef(scan) = 0.23/100;
            case 13
                im = imread('negative\23.bmp'); % reading 1.png
                tifname='23';
                L_ef(scan) = 9.45/100;
                Er_L_ef(scan) = 0.25/100;
            case 14
                im = imread('negative\24.bmp');
                tifname='24';
                L_ef(scan) = 9.16/100;
                Er_L_ef(scan) = 0.27/100;
            case 15
                im = imread('negative\25.bmp');
                tifname='25';
                L_ef(scan) = 8.80/100;
                Er_L_ef(scan) = 0.26/100;
            case 16
                im = imread('negative\26.bmp');
                tifname='26';
                L_ef(scan) = 8.26/100;
                Er_L_ef(scan) = 0.38/100;
            case 17
                im = imread('positive\27.bmp');
                tifname='27';
                L_ef(scan) = 8.88/100;
                Er_L_ef(scan) = 0.25/100;
            case 18
                im = imread('positive\28.bmp');
                tifname='28';
                L_ef(scan) = 9.12/100;
                Er_L_ef(scan) = 0.22/100;
            case 19
                im = imread('positive\29.bmp');
                tifname='29';
                L_ef(scan) = 9.28/100;
                Er_L_ef(scan) = 0.19/100;
            case 20
                im = imread('positive\30.bmp');
                tifname='30';
                L_ef(scan) = 9.47/100;
                Er_L_ef(scan) = 0.27/100;
            case 21
                im = imread('positive\31.bmp');
                tifname='31';
                L_ef(scan) = 9.55/100;
                Er_L_ef(scan) = 0.20/100;

            case 22
                im = imread('positive\32.bmp'); % reading 1.png
                tifname='32';
                L_ef(scan) = 9.66/100;
                Er_L_ef(scan) = 0.29/100;
            case 23
                im = imread('positive\33.bmp'); % reading 1.png
                tifname='33';
                L_ef(scan) = 9.73/100;
                Er_L_ef(scan) = 0.22/100;
            case 24
                im = imread('positive\34.bmp');
                tifname='34';
                L_ef(scan) = 9.80/100;
                Er_L_ef(scan) = 0.18/100;
            case 25
                im = imread('positive\35.bmp');
                tifname='35';
                L_ef(scan) = 9.72/100;
                Er_L_ef(scan) = 0.17/100;
            case 26
                im = imread('positive\36.bmp');
                tifname='36';
                L_ef(scan) = 9.81/100;
                Er_L_ef(scan) = 0.12/100;
            case 27
                im = imread('positive\37.bmp');
                tifname='37';
                L_ef(scan) = 9.84/100;
                Er_L_ef(scan) = 0.10/100;
            case 28
                im = imread('positive\38.bmp');
                tifname='38';
                L_ef(scan) = 9.90/100;
                Er_L_ef(scan) = 0.11/100;
            case 29
                im = imread('positive\39.bmp');
                tifname='39';
                L_ef(scan) = 9.91/100;
                Er_L_ef(scan) = 0.12/100;
            case 30
                im = imread('positive\40.bmp');
                tifname='40';
                L_ef(scan) = 9.92/100;
                Er_L_ef(scan) = 0.16/100;
            case 31
                im = imread('positive\41.bmp');
                tifname='41';
                L_ef(scan) = 9.94/100;
                Er_L_ef(scan) = 0.18/100;
        end
        sg = im; % 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;

        %plot(x,profx,'k');
        %hold on;
        %plot(y,profy,'k');
        %hold on;

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

        a1 = 100;
        a2 = 7700;
        a3 = 280;
        a4 = 39;
        a5 = 1;

        %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_x,fval_x,flag_x,err_x] = fminsearch(@SupGau_devsum,a0,opts,x,profx);
        [fitpara_y,fval_y,flag_y,err_y] = fminsearch(@SupGau_devsum,a0,opts,y,profy);

        base_x = fitpara_x(1);
        amplitude_x = fitpara_x(2);
        peak_center_x = fitpara_x(3);
        sig0_x = fitpara_x(4);
        N_x = fitpara_x(5);
        sig_x(scan,imNo) = sig0_x*(pi/2).^(2/fitpara_x(5)-1);

        base_y = fitpara_y(1);
        amplitude_y = fitpara_y(2);
        peak_center_y = fitpara_y(3);
        sig0_y = fitpara_y(4);
        N_y = fitpara_y(5);
        sig_y(scan,imNo) = sig0_y*(pi/2).^(2/fitpara_y(5)-1);

        %display('I(A)  Image    Base    Amplitude     PeakCenter     N       Sigma_0     Sigma');
        string_x = [
            num2str(I(scan),'  %2.2f'),'    ', ...
            num2str(imNo,'%d'),'      ', ...
            num2str(base_x,'%12.2f'),'      ', ...
            num2str(amplitude_x,'%12.f'),'         ', ...
            num2str(peak_center_x,'%12.2f'),'      ', ...
            num2str(N_x,'%12.4f'),'    ', ...
            num2str(sig0_x,'%12.2f'),'      ', ...
            num2str(sig_x(scan,imNo),'%12.2f'),'    ', ...
            ];
        %disp(string_x);
        string_y = [
            num2str(I(scan),'  %2.2f'),'    ', ...
            num2str(imNo,'%d'),'      ', ...
            num2str(base_y,'%12.2f'),'      ', ...
            num2str(amplitude_y,'%12.f'),'         ', ...
            num2str(peak_center_y,'%12.2f'),'      ', ...
            num2str(N_y,'%12.4f'),'    ', ...
            num2str(sig0_y,'%12.2f'),'      ', ...
            num2str(sig_y(scan,imNo),'%12.2f'),'    ', ...
            ];
        %disp(string_y);
        %sig

        %plot(x,profx,'k');
        %hold on;
        %plot(y,profy,'k');
        %hold on;

        y_fit_x = fitpara_x(1)+fitpara_x(2)*exp( -0.5*(abs(x-fitpara_x(3))/(fitpara_x(4))).^fitpara_x(5) );
        y_fit_y = fitpara_y(1)+fitpara_y(2)*exp( -0.5*(abs(y-fitpara_y(3))/(fitpara_y(4))).^fitpara_y(5) );
        %y_fit =   fitpara(1)  +fitpara(2)  *exp( -0.5*(abs(X-fitpara(3))/(fitpara(4))).^fitpara(5) );
        
        %plot(x,y_fit_x,'r');
        
        %plot(y,y_fit_y,'r');
        drawnow
    end
end


for scan=1:scan_current_number
    
    % find mean and standar deviation
    sig_mean_x(scan) = mean(sig_x(scan,:));
    sig_x_er(scan) = std(sig_x(scan,:));
    sig_mean_y(scan) = mean(sig_y(scan,:));
    sig_y_er(scan) = std(sig_y(scan,:));
    
    % conver to length (in m unit) from pixel
    sig_mean_x(scan)=sig_mean_x(scan)*calibration_x; %
    sig_mean_y(scan)=sig_mean_y(scan)*calibration_y; %
    sig_x_er(scan)= sqrt((sig_x_er(scan)*calibration_x).^2+(sig_mean_x(scan)*er_calibration_x).^2); %
    sig_y_er(scan)= sqrt((sig_y_er(scan)*calibration_y).^2+(sig_mean_y(scan)*er_calibration_y).^2); %
    
    % find sigma squared
    sig_sqr_x(scan) = sig_mean_x(scan)*sig_mean_x(scan);
    sig_sqr_x_er(scan) = 2*sig_mean_x(scan)*sig_x_er(scan);
    sig_sqr_y(scan) = sig_mean_y(scan)*sig_mean_y(scan);
    sig_sqr_y_er(scan) = 2*sig_mean_y(scan)*sig_y_er(scan);
    
       
    %I(scan) = I_start + (scan-1)*I_increment;
    %plot(I,sig_mean,'r');
    k1(scan) =0.2998*(3.6*0.0001+1945*0.000001*I(scan))/(0.0254*p);

    k1_Er(1) =0.2998*((3.6+1.3)*0.0001+(1945+2)*0.000001*I(scan))/(0.0254*p);
    k1_Er(2) =0.2998*((3.6-1.3)*0.0001+(1945+2)*0.000001*I(scan))/(0.0254*p);
    k1_Er(3) =0.2998*((3.6+1.3)*0.0001+(1945-2)*0.000001*I(scan))/(0.0254*p);
    k1_Er(4) =0.2998*((3.6-1.3)*0.0001+(1945-2)*0.000001*I(scan))/(0.0254*p);
    er_k1(1)=abs(k1(scan)-k1_Er(1));
    er_k1(2)=abs(k1(scan)-k1_Er(2));
    er_k1(3)=abs(k1(scan)-k1_Er(3));
    er_k1(4)=abs(k1(scan)-k1_Er(4));

    k1_er(scan)=max(er_k1);
    
    k1L(scan) = L_ef(scan)*k1(scan); %k1*L
    er_k1L(scan) = sqrt(  (L_ef(scan)*k1_er(scan) )^2  + (Er_L_ef(scan)*k1(scan))^2  );
    %er_k1L(scan) = L_ef*k1_er(scan);
    
    FitDatX = [k1L',er_k1L',sig_sqr_x',sig_sqr_x_er']
    FitDatY = [k1L',er_k1L',sig_sqr_y',sig_sqr_y_er']
end

%Emit_Parabola_Fit_kl_XProjection(k1L,er_k1L,sig_sqr_x,sig_sqr_x_er,me,En,S12,erS12,scan_current_number);
Emit_Parabola_Fit_kl_XProjection(k1L(16:31),er_k1L(16:31),sig_sqr_x(16:31),sig_sqr_x_er(16:31),me,En,S12,erS12,16);
%Emit_Parabola_Fit_kl_YProjection(k1L,er_k1L,sig_sqr_y,sig_sqr_y_er,me,En,S12,erS12,scan_current_number);
Emit_Parabola_Fit_kl_YProjection(k1L(1:16),er_k1L(1:16),sig_sqr_y(1:16),sig_sqr_y_er(1:16),me,En,S12,erS12,16);