function varargout = CltvGUI(varargin)

% CLTVGUI M-file for CltvGUI.fig

%      CLTVGUI, by itself, creates a new CLTVGUI or raises the existing

%      singleton*.

%

%      H = CLTVGUI returns the handle to a new CLTVGUI or the handle to

%      the existing singleton*.

%

%      CLTVGUI('CALLBACK',hObject,eventData,handles,...) calls the local

%      function named CALLBACK in CLTVGUI.M with the given input arguments.

%

%      CLTVGUI('Property','Value',...) creates a new CLTVGUI or raises the

%      existing singleton*.  Starting from the left, property value pairs are

%      applied to the GUI before CltvGUI_OpeningFunction gets called.  An

%      unrecognized property name or invalid value makes property application

%      stop.  All inputs are passed to CltvGUI_OpeningFcn via varargin.

%

%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one

%      instance to run (singleton)".

%

% See also: GUIDE, GUIDATA, GUIHANDLES



% Edit the above text to modify the response to help CltvGUI



% Last Modified by GUIDE v2.5 10-Dec-2003 00:40:06


% Begin initialization code - DO NOT EDIT

gui_Singleton = 1;

gui_State = struct('gui_Name',       mfilename, ...

                   'gui_Singleton',  gui_Singleton, ...

                   'gui_OpeningFcn', @CltvGUI_OpeningFcn, ...

                   'gui_OutputFcn',  @CltvGUI_OutputFcn, ...

                   'gui_LayoutFcn',  [] , ...

                   'gui_Callback',   []);

if nargin & isstr(varargin{1})

    gui_State.gui_Callback = str2func(varargin{1});

end



if nargout

    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});

else

    gui_mainfcn(gui_State, varargin{:});

end

% End initialization code - DO NOT EDIT





% --- Executes just before CltvGUI is made visible.

function CltvGUI_OpeningFcn(hObject, eventdata, handles, varargin)

% This function has no output args, see OutputFcn.

% hObject    handle to figure

% eventdata  reserved - to be defined in a future version of MATLAB

% handles    structure with handles and user data (see GUIDATA)

% varargin   command line arguments to CltvGUI (see VARARGIN)



% Choose default command line output for CltvGUI

handles.output = hObject;



% Initialize the input;

handles.Current     =       500;                %Current in A

handles.Radius      =       5*10^(-6);          %Beam radius in m

handles.Energy      =       10;                 %Beam energy in MeV

handles.PulseTime   =       300*10^(-15);        %Pulse duration in s

handles.Density     =       1.33*10^(23);       %Electron's number density in #/m^3



handles.EnergySpread    =   0;                  %Beam energy spread

handles.Emittance   =       0;                  %Beam transverse emittance in m.rad



handles.a           =       1;                  %Normalized vector potential

handles.lambda      =       1*10^(-6);          %Laser wave length

handles.Cycles      =       100;                %The initial interaction length in terms of the optical cycle.



handles.seed        =       7;                 %Assume a noise seeding;

handles.NumofRun    =       0;

% Update handles structure

guidata(hObject, handles);



% UIWAIT makes CltvGUI wait for user response (see UIRESUME)

% uiwait(handles.figure1);





% --- Outputs from this function are returned to the command line.

function varargout = CltvGUI_OutputFcn(hObject, eventdata, handles)

% varargout  cell array for returning output args (see VARARGOUT);

% hObject    handle to figure

% eventdata  reserved - to be defined in a future version of MATLAB

% handles    structure with handles and user data (see GUIDATA)



% Get default command line output from handles structure

varargout{1} = handles.output;



% --- Executes during object creation, after setting all properties.
function InputCurrent_CreateFcn(hObject, eventdata, handles)
% hObject    handle to InputCurrent (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end



function InputCurrent_Callback(hObject, eventdata, handles)
% hObject    handle to InputCurrent (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of InputCurrent as text
%        str2double(get(hObject,'String')) returns contents of InputCurrent as a double
handles.Current     =       str2double(get(hObject, 'String'));

handles.Density     =       CalNumDensity(handles);

guidata(hObject, handles);

h       =       findobj('Tag', 'NumberDensity');

set(h, 'String', num2str(handles.Density));




% --- Executes during object creation, after setting all properties.
function InputVoltage_CreateFcn(hObject, eventdata, handles)
% hObject    handle to InputVoltage (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end



function InputVoltage_Callback(hObject, eventdata, handles)
% hObject    handle to InputVoltage (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of InputVoltage as text
%        str2double(get(hObject,'String')) returns contents of InputVoltage as a double
handles.Energy      =       str2double(get(hObject, 'String'));

handles.Density     =       CalNumDensity(handles);

guidata(hObject, handles);

h       =       findobj('Tag', 'NumberDensity');

set(h, 'String', num2str(handles.Density));



% --- Executes during object creation, after setting all properties.
function InputRadius_CreateFcn(hObject, eventdata, handles)
% hObject    handle to InputRadius (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end



function InputRadius_Callback(hObject, eventdata, handles)
% hObject    handle to InputRadius (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of InputRadius as text
%        str2double(get(hObject,'String')) returns contents of InputRadius as a double
handles.Radius      =       str2double(get(hObject, 'String'))*10^(-6);

handles.Density     =       CalNumDensity(handles);

guidata(hObject, handles);

h       =       findobj('Tag', 'NumberDensity');

set(h, 'String', num2str(handles.Density));



% --- Executes during object creation, after setting all properties.
function InputPulseTime_CreateFcn(hObject, eventdata, handles)
% hObject    handle to InputPulseTime (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end



function InputPulseTime_Callback(hObject, eventdata, handles)
% hObject    handle to InputPulseTime (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of InputPulseTime as text
%        str2double(get(hObject,'String')) returns contents of InputPulseTime as a double
handles.PulseTime   =       str2double(get(hObject, 'String'))*10^(-15);

handles.Density     =       CalNumDensity(handles);

guidata(hObject, handles);

h       =       findobj('Tag', 'NumberDensity');

set(h, 'String', num2str(handles.Density));



% --- Executes during object creation, after setting all properties.
function InputEnergySpread_CreateFcn(hObject, eventdata, handles)
% hObject    handle to InputEnergySpread (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end



function InputEnergySpread_Callback(hObject, eventdata, handles)
% hObject    handle to InputEnergySpread (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of InputEnergySpread as text
%        str2double(get(hObject,'String')) returns contents of InputEnergySpread as a double
handles.EnergySpread    =       str2double(get(hObject, 'String'));

guidata(hObject, handles);



% --- Executes during object creation, after setting all properties.
function InputEmittance_CreateFcn(hObject, eventdata, handles)
% hObject    handle to InputEmittance (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end



function InputEmittance_Callback(hObject, eventdata, handles)
% hObject    handle to InputEmittance (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of InputEmittance as text
%        str2double(get(hObject,'String')) returns contents of InputEmittance as a double
handles.Emittance   =       str2double(get(hObject, 'String'));

guidata(hObject, handles);



% --- Executes during object creation, after setting all properties.
function Input_a_CreateFcn(hObject, eventdata, handles)
% hObject    handle to Input_a (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end



function Input_a_Callback(hObject, eventdata, handles)
% hObject    handle to Input_a (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of Input_a as text
%        str2double(get(hObject,'String')) returns contents of Input_a as a double
handles.a       =       str2double(get(hObject, 'String'));

guidata(hObject, handles);



% --- Executes during object creation, after setting all properties.
function InputWaveLength_CreateFcn(hObject, eventdata, handles)
% hObject    handle to InputWaveLength (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end



function InputWaveLength_Callback(hObject, eventdata, handles)
% hObject    handle to InputWaveLength (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of InputWaveLength as text
%        str2double(get(hObject,'String')) returns contents of InputWaveLength as a double
handles.lambda      =       str2double(get(hObject, 'String'));

guidata(hObject, handles);



% --- Executes during object creation, after setting all properties.
function ResultList_CreateFcn(hObject, eventdata, handles)
% hObject    handle to ResultList (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: listbox controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end


% --- Executes on selection change in ResultList.
function ResultList_Callback(hObject, eventdata, handles)
% hObject    handle to ResultList (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns ResultList contents as cell array
%        contents{get(hObject,'Value')} returns selected item from ResultList


% --- Executes on button press in Run.
function Run_Callback(hObject, eventdata, handles)
% hObject    handle to Run (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)



% Initialize a progress bar

info.color  =       [0, 1, 0];

info.title  =       'Progress Bar';

info.msg    =       'Simulation Running, Please Wait';

pb          =       progbar(info);





handles.Cycles  =   500;

guidata(hObject, handles);

L1          =       handles.Cycles*handles.lambda;
Rslt1       =       CalCollective(handles);

xi          =       min(Rslt1.omega):(max(Rslt1.omega)-min(Rslt1.omega))/3000:max(Rslt1.omega);

gain1       =       interp1(Rslt1.omega, Rslt1.gain, xi, 'cubic', 0);

eta1        =       interp1(Rslt1.omega, Rslt1.eta, xi, 'cubic', 0);

[Maxeta, index2]   =       max(eta1);

[MaxGain, index]   =       max(gain1);



%update progess bar

figure(pb);

progbar(pb, 30);



%handles.Cycles  =   15;

%guidata(hObject, handles);

%L2          =       handles.Cycles*handles.lambda;

%Rslt2       =       CalCollective(handles);

%gain2       =       interp1(Rslt2.omega, Rslt2.gain, xi, 'cubic', 0);
%h1          =       axes('Tag', 'GainPlot');

%h2          =       axes('Tag', 'ExpPlot');

%plot(Results.omega, Results.gain, 'k', 'LineWidth', 1.5);

%xlabel('\omega/\omega_0', 'FontSize', 10);

%ylabel('P(z)/P(0)', 'FontSize', 10);

%set(h, 'FontSize', 10, 'LineWidth', 1.5);

%legend(Results.str);



TotalGain   =       eta1(index)*handles.Current*handles.Energy*10^6/handles.seed;

%Lgain       =       (L1-L2)/log(gain1(index)/gain2(index));

%Ls          =       L1 + log(TotalGain/gain1(index))*Lgain;

i           =       1;

gain(1)     =       1;

flag        =       0;

while (gain(i) <= TotalGain)

    if flag == 0

        L           =       handles.lambda;

        handles.Cycles  =   1;

        guidata(hObject, handles);

        TempRslt    =       CalCollective(handles, xi(index));

        gain(1)     =       TempRslt.gain;

        distance(1) =       handles.lambda;

        eta(1)      =       TempRslt.eta;

        flag        =       1;



    else

        i           =       i + 1;

        L           =       handles.lambda + handles.lambda*i;

        handles.Cycles  =   L/handles.lambda;

        guidata(hObject, handles);

        TempRslt    =       CalCollective(handles, xi(index));

        gain(i)     =       TempRslt.gain;

        distance(i) =       L;

        eta(i)      =       TempRslt.eta;

        

%        axes(handles.ExpPlot);

%        plot(distance/handles.lambda, log10(gain), 'k', 'LineWidth', 1.5);

%        drawnow;



        figure(pb);

        progress    =       30 + log(gain(i-1))/log(TotalGain)*30;

        progbar(pb, progress);

    end



    %axes(handles.ExpPlot);

    %plot(distance, log10(gain), 'rs');

    %drawnow;

end



L1          =       (TempRslt.Cycles-1)*handles.lambda;

handles.Cycles  =   TempRslt.Cycles - 1;

guidata(hObject, handles);



Rslt1       =       CalCollective(handles);

xi          =       min(Rslt1.omega):(max(Rslt1.omega)-min(Rslt1.omega))/3000:max(Rslt1.omega);

gain1       =       interp1(Rslt1.omega, Rslt1.gain, xi, 'cubic', 0);

eta1        =       interp1(Rslt1.omega, Rslt1.eta, xi, 'cubic', 0);

[Maxeta, index2]   =       max(eta1);

[MaxGain, index]   =       max(gain1);



axes(handles.GainPlot);

%plot(xi, gain1, 'k', 'LineWidth', 1.5, xi(index), gain1(index), '*');

%plot(xi, gain1, 'k', xi(index), gain1(index), 'r*', 'LineWidth', 1.5, 'MarkerSize', 8);

plot(xi, gain1, 'k', 'LineWidth', 1.5);

xlabel('\omega/\omega_0', 'FontSize', 10);

ylabel('P(z)/P(0)', 'FontSize', 10);

set(gca, 'FontSize', 10, 'LineWidth', 1.5);

legend(Rslt1.str);

drawnow;



% Update Progress;

figure(pb);

progbar(pb, 90);



div         =       (max(distance) - min(distance))/handles.lambda + 500;

Intpxi      =       min(distance):(max(distance)-min(distance))/div:max(distance);

IntpGain    =       interp1(distance, gain, Intpxi, 'cubic');

axes(handles.ExpPlot);

%plot(distance/handles.lambda, log10(gain), 'rs',  xi/handles.lambda, log10(IntpGain), 'k', 'LineWidth', 1.5, 'MarkerEdgeColor', 'k', 'MarkerFaceColor', 'g', 'MarkerSize', 5);

plot(Intpxi/handles.lambda, log10(IntpGain), 'k', 'LineWidth', 1.5);

xlabel('L/\lambda_0', 'FontSize', 10);

ylabel('log(Gain)', 'FontSize', 10);

set(gca, 'FontSize', 10, 'LineWidth', 1.5);



figure(5);

%plot(distance/handles.lambda, log10(gain), 'rs',  xi/handles.lambda, log10(IntpGain), 'k', 'LineWidth', 1.5, 'MarkerEdgeColor', 'k', 'MarkerFaceColor', 'g', 'MarkerSize', 5);

plot(Intpxi/handles.lambda, log10(IntpGain), 'k', 'LineWidth', 1.5);

xlabel('L/\lambda_0', 'FontSize', 10);

ylabel('log(Gain)', 'FontSize', 10);

set(gca, 'FontSize', 10, 'LineWidth', 1.5);



figure(6);

%plot(xi, gain1, 'k', 'LineWidth', 1.5, xi(index), gain1(index), '*');

plot(xi, gain1, 'k', 'LineWidth', 1.5);

[FWHM_m, FWHM_n] = size(xi)

for i=1:FWHM_n

    difference(i) = abs(gain1(i)-gain1(index)/2);

end

[FWHM_diff, FWHM_index] =   min(difference);

FWHM    =   abs(xi(index) - xi(FWHM_index))/xi(index);

xlabel('\omega/\omega_0', 'FontSize', 10);

ylabel('P(z)/P(0)', 'FontSize', 10);

set(gca, 'FontSize', 10, 'LineWidth', 1.5);

legend(Rslt1.str);

FWHM_str=num2str(FWHM);

legend(FWHM_str)



% Update Progress;

figure(pb);

progbar(pb, 95);



[m, n]      =       size(eta);

handles.NumofRun    =   handles.NumofRun + 1;

guidata(hObject, handles);

strName     =       strvcat('List of Results: Run', 'a', 'EnergySpread', 'LaserWaveLength(um)', 'NumberDensity(#/m^3)', 'Emittance', 'Efficiency', 'SaturationGain', 'SaturationPower(W)', 'SaturationLength(m)', 'NormalizedFreq');

str1        =       num2str(handles.NumofRun);

str2        =       num2str(handles.a);

str3        =       num2str(handles.EnergySpread);

str4        =       num2str(handles.lambda);

str5        =       num2str(handles.Density);

str6        =       num2str(handles.Emittance);

str7        =       num2str(eta(n-1));

str8        =       num2str(TotalGain);

str9        =       num2str(eta(n-1)*handles.Current*handles.Energy*10^6);

str10       =       num2str(L1);

str11       =       num2str(xi(index));



strValue    =       strvcat(str1, str2, str3, str4, str5, str6, str7, str8, str9, str10, str11);

str         =       strvcat(strName, strValue);

str         =       cat(2, str(1:11, :), str(12:22, :));



set(handles.ResultList, 'String', str);

%axes(handles.ExpPlot);

%plot(distance, log10(gain));



%close progress bar;

figure(pb);

progbar(pb, 100);

progbar(pb, -1);   

    

% --------------------------------------------------------------------
function Untitled_1_Callback(hObject, eventdata, handles)
% hObject    handle to Untitled_1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)


function Density = CalNumDensity(handles)

% Get new beam parameters and re-calculate the electron's number density;



gamma               =       handles.Energy/0.511 + 1;

velocity            =       sqrt(1-1/gamma^2)*3.0*10^8;

volume              =       pi*handles.Radius^2*velocity*handles.PulseTime;

Density             =       handles.Current*handles.PulseTime/(1.602*10^(-19))/volume;





function Results = CalCollective(handles, OmegaOfIntrst)

%Calculate the collective gain and efficiency;



% Initialization

a           =       handles.a;

lambda      =       handles.lambda;

n_0         =       handles.Density;

gamma_0     =       handles.Energy/0.511 + 1;  %(MeV)

E_0         =       handles.Energy; %(MeV)

E_spread    =       handles.EnergySpread; 

r_b0        =       handles.Radius;

l           =       handles.Cycles*lambda;

%transvers emittance

epsilon_tr  =       handles.Emittance;







% Physics constants

e           =       1.602*10^(-19);

epsilon     =       8.8542*10^(-12);

m           =       9.1095*10^(-31);

c           =       3*10^8;



% Electron dynamics in laser field

gamma_0     =       E_0/0.511 + 1;  

E_1         =       E_0 - E_spread*E_0;

E_2         =       E_0 + E_spread*E_0;

gamma_low   =       E_1/0.511 + 1;   %low end of the thermal spread;

gamma_high  =       E_2/0.511 + 1;   %high end of the thermal spread;



beta_z0     =       sqrt(1 - 1/gamma_0^2);

gamma_avg   =       gamma_0 + a^2/(4*gamma_0*(1+beta_z0));

beta_z_avg  =       abs((gamma_avg - gamma_0*(1 + beta_z0))/gamma_avg);  %transvers avg velocity;

gamma_0z    =       gamma_avg/sqrt(1 + a^2);

momentum_avg=       abs(gamma_avg - gamma_0*(1 + beta_z0));



beta_low    =       sqrt(1 - 1/gamma_low^2);

gamma_avg_l =       gamma_low + a^2/(4*gamma_low*(1+beta_low));

momentum_l  =       abs(gamma_avg_l - gamma_low*(1 + beta_low));



beta_high   =       sqrt(1 - 1/gamma_high^2);

gamma_avg_h =       gamma_high + a^2/(4*gamma_high*(1+beta_high));

momentum_h  =       abs(gamma_avg_h - gamma_high*(1 + beta_high));



if momentum_h==momentum_l

    x1      =       -1/2;

    x2      =       1/2;

else

    x1          =       (momentum_l - momentum_avg)/(momentum_h - momentum_l);

    x2          =       (momentum_h - momentum_avg)/(momentum_h - momentum_l);

end



% Longitudinal velocity spread;

v_zth       =       c*(momentum_h - momentum_l)/gamma_avg/gamma_0z^2;

r           =       v_zth/c/beta_z_avg;

k_0         =       2*pi/lambda;

% Tansverse velocity spread;

v_zth_epsilon   =   (epsilon_tr/pi)^2/2/r_b0^2*c;



omega_p     =       sqrt(e^2*n_0/gamma_0z^2/gamma_avg/m/epsilon);

theta_p     =       omega_p/(beta_z_avg*c);



if nargin==1

    for j   =   1:2001;

        theta_tilde(j)  =       (j-1001)/10;

        omega_rad       =       (theta_tilde(j)/l + k_0 + k_0/beta_z_avg)/(1/beta_z_avg/c - 1/c);

        lambda_rad      =       c/omega_rad*2*pi;

        %   With the nonlinear correction, i.e. JJ term or KK term. 

        kappa           =       1/2/pi*(a*2*pi)^2/lambda^2*(1 + beta_z_avg)^2*gamma_0z^4/gamma_avg^2*lambda_rad*(besselj(0, a^2/2/(1+a^2)) - besselj(1, a^2/2/(1+a^2))); 

        % Without the nonlinear correction

        %    kappa           =       1/2/pi*(a*2*pi)^2/lambda^2*(1 + beta_z_avg)^2*gamma_0z^4/gamma_avg^2*lambda_rad;

        Q               =       kappa*theta_p^2;

    

        theta_p_tilde   =   theta_p*l;

        kappa_tilde     =   kappa*l;

        %theta_p_tilde   =  1;      % test for Gover's

        %kappa_tilde     =  1;      % test for Gover's   

        %theta_th_tilde  =  10;     % test

        Q_tilde         =       kappa_tilde*theta_p_tilde^2;

        Q_tilde1(j)     =       Q_tilde;

        %velocity spread due to the energy spread

        theta_th_tilde_E=       (theta_tilde(j) + l/k_0 + 2*pi*l/lambda_rad)*r;

        %spread due to transverse emittance

        theta_th_tilde_epsilon  =   epsilon_tr^2*l/pi/r_b0^2/lambda_rad;

    

        %sum up the spread due to different sources linearly;

        %    theta_th_tilde  =   theta_th_tilde_E + theta_th_tilde_epsilon;

    

        %sum up the spread due to different sources in quadratures

        theta_th_tilde  =   sqrt(theta_th_tilde_E^2 + theta_th_tilde_epsilon^2);

        

        b_up        =       [0 1+x1*x2*r^2+(x1+x2)*r -2*theta_tilde(j)+(x1+x2)*theta_th_tilde-2*theta_tilde(j)*r*(x1+x2)+2*r*theta_th_tilde*x1*x2-2*theta_tilde(j)*r^2*x1*x2 theta_tilde(j)^2+x1*x2*theta_th_tilde^2-2*r*theta_th_tilde*theta_tilde(j)*x1*x2+theta_tilde(j)^2*r*(x1+x2)-theta_tilde(j)*theta_th_tilde*(x1+x2)+x1*x2*theta_tilde(j)^2*r^2-theta_p_tilde^2];

        d_down      =       i*[1+x1*x2*r^2+(x1+x2)*r -2*theta_tilde(j)+(x1+x2)*theta_th_tilde-2*theta_tilde(j)*r*(x1+x2)+2*r*theta_th_tilde*x1*x2-2*theta_tilde(j)*r^2*x1*x2 theta_tilde(j)^2+x1*x2*theta_th_tilde^2-2*r*theta_th_tilde*theta_tilde(j)*x1*x2+theta_tilde(j)^2*r*(x1+x2)-theta_tilde(j)*theta_th_tilde*(x1+x2)+x1*x2*theta_tilde(j)^2*r^2-theta_p_tilde^2 Q_tilde];

        [r2, p, k]  =       residue(b_up, d_down);

        sz          =       (p/l + 2*pi/lambda_rad)*l;    

        Magnitude   =       sum(r2.*exp(i*sz));

        sj_1(j)     =       p(1)/l;

        sj_2(j)     =       p(2)/l;

        sj_3(j)     =       p(3)/l;

        remains_1(j)=       r2(1);

        remains_2(j)=       r2(2);

        remains_3(j)=       r2(3);

    

        gain(j)     =       Magnitude*conj(Magnitude);

    end



    %omega           =       (theta_tilde(index)/l + k_0 + k_0/beta_z_avg)/(1/beta_z_avg/c - 1/c)/(c*k_0);

    omega_tilde     =       (theta_tilde/l + k_0 + k_0/beta_z_avg)/(1/beta_z_avg/c - 1/c)/(c*k_0);

    eta_tilde       =       real(sj_2) - theta_tilde/l;

    eta             =       beta_z_avg^3/(1-1/gamma_avg^2)*gamma_0z^2*(1+1/gamma_avg)*eta_tilde*lambda./omega_tilde/pi;

    [Maxgain, index] =       max(gain);



    str             =       sprintf('a=%g, n_0=%g/m^3\nE=%gMeV, \\DeltaE/E=%g\nL=%g', a, n_0, E_0,E_spread, l);



    Results.omega   =       omega_tilde;

    Results.gain    =       gain;

    Results.eta     =       eta;

    Results.str     =       str;



elseif nargin==2

    for j   =   1:1;

        theta_tilde(j)  =       (OmegaOfIntrst*(1/beta_z_avg/c - 1/c)*c*k_0 - k_0 - k_0/beta_z_avg)*l;

        omega_rad       =       (theta_tilde(j)/l + k_0 + k_0/beta_z_avg)/(1/beta_z_avg/c - 1/c);

        lambda_rad      =       c/omega_rad*2*pi;

        kappa           =       1/2/pi*(a*2*pi)^2/lambda^2*(1 + beta_z_avg)^2*gamma_0z^4/gamma_avg^2*lambda_rad*(besselj(0, a^2/2/(1+a^2)) - besselj(1, a^2/2/(1+a^2))); 

        Q               =       kappa*theta_p^2;

    

        theta_p_tilde   =       theta_p*l;

        kappa_tilde     =       kappa*l;

        Q_tilde         =       kappa_tilde*theta_p_tilde^2;

        Q_tilde1(j)     =       Q_tilde;

        %velocity spread due to the energy spread

        theta_th_tilde_E=       (theta_tilde(j) + l/k_0 + 2*pi*l/lambda_rad)*r;

        %spread due to transverse emittance

        theta_th_tilde_epsilon  =   epsilon_tr^2*l/pi/r_b0^2/lambda_rad;

    

        %sum up the spread due to different sources linearly;

        %    theta_th_tilde  =   theta_th_tilde_E + theta_th_tilde_epsilon;

    

        %sum up the spread due to different sources in quadratures

        theta_th_tilde  =   sqrt(theta_th_tilde_E^2 + theta_th_tilde_epsilon^2);

        

        b_up        =       [0 1+x1*x2*r^2+(x1+x2)*r -2*theta_tilde(j)+(x1+x2)*theta_th_tilde-2*theta_tilde(j)*r*(x1+x2)+2*r*theta_th_tilde*x1*x2-2*theta_tilde(j)*r^2*x1*x2 theta_tilde(j)^2+x1*x2*theta_th_tilde^2-2*r*theta_th_tilde*theta_tilde(j)*x1*x2+theta_tilde(j)^2*r*(x1+x2)-theta_tilde(j)*theta_th_tilde*(x1+x2)+x1*x2*theta_tilde(j)^2*r^2-theta_p_tilde^2];

        d_down      =       i*[1+x1*x2*r^2+(x1+x2)*r -2*theta_tilde(j)+(x1+x2)*theta_th_tilde-2*theta_tilde(j)*r*(x1+x2)+2*r*theta_th_tilde*x1*x2-2*theta_tilde(j)*r^2*x1*x2 theta_tilde(j)^2+x1*x2*theta_th_tilde^2-2*r*theta_th_tilde*theta_tilde(j)*x1*x2+theta_tilde(j)^2*r*(x1+x2)-theta_tilde(j)*theta_th_tilde*(x1+x2)+x1*x2*theta_tilde(j)^2*r^2-theta_p_tilde^2 Q_tilde];

        [r2, p, k]  =       residue(b_up, d_down);

        sz          =       (p/l + 2*pi/lambda_rad)*l;    

        Magnitude   =       sum(r2.*exp(i*sz));

        sj_1(j)     =       p(1)/l;

        sj_2(j)     =       p(2)/l;

        sj_3(j)     =       p(3)/l;

        remains_1(j)=       r2(1);

        remains_2(j)=       r2(2);

        remains_3(j)=       r2(3);

    

        gain(j)     =       Magnitude*conj(Magnitude);

    end        

    omega_tilde     =       (theta_tilde/l + k_0 + k_0/beta_z_avg)/(1/beta_z_avg/c - 1/c)/(c*k_0);

    eta_tilde       =       real(sj_2) - theta_tilde/l;

    eta             =       beta_z_avg^3/(1-1/gamma_avg^2)*gamma_0z^2*(1+1/gamma_avg)*eta_tilde*lambda./omega_tilde/pi;



    Results.gain    =       gain;

    Results.eta     =       eta;

    Results.Cycles  =       l/handles.lambda;

end



function [handle]=progbar(varargin)



% progbar.m: general purpose progress bar for MATLAB

% type 'progbar demo' to get started

%

% pb=progbar(info)             initialises a new progbar

%    progbar(pb,progress)      updates an existing progbar

%    progbar(pb,progress,msg)  update the progress and message on an existing progbar

%    progbar(pb,-1)            closes an existing progbar (no error if closed by user)

%    close(pb)                 alternative way to close (error if already closed)

%

% info, if supplied, should be a structure with some

% or all of the following fields (listed as =<default>)

%

%    title='Progress'  title

%    msg=''            message

%    size=1            small integer dictating size of bar

%    period=0.1        a time in seconds dictating the minimum update

%                      rate of the progress bar

%    pos='centre'      any of 'topleft,top,top right,centreleft,centre,

%                        centreright,bottomleft,bottom,bottomright'. note

%                        american spelling accepted, but frowned upon.

%    color='blue'      RGB triple in MATLAB format, [R G B]. if you have

%                        the color tool 'rgb.m' on your path, you may use

%                        any format that rgb accepts. rgb.m is available

%                        at www.mathworks.com/matlabcentral

%    clearance=0       number between 0 and 1, affects appearance



% VERSION:  06/06/2002

% AUTHOR:   ben mitch

% CONTACT:  footballpitch@hotmail.com

% LOCATION: interaction\guis\



if nargin==0 info.dummy='dummy'; end

if nargin==1 info=varargin{1}; end



if nargin<2

	% check type

	if isa(info,'char') & strcmp(info,'demo')

		progbar_demo

		return

	end

	

	if ~isstruct(info)

		error('Single argument should be an info structure or ''demo''');

	end

	

	% extract

	if ~isfield(info,'title') info.title='Progress'; end

	if ~isfield(info,'msg') info.msg=''; end

	if ~isfield(info,'size') info.size=1; end

	if ~isfield(info,'period') info.period=0.1; end

	if ~isfield(info,'pos') info.pos='centre'; end

	if ~isfield(info,'color') info.color=[0 0 1]; end

	if ~isfield(info,'clearance') info.clearance=0; end

	

	% check values

	if floor(info.size)~=info.size | info.size<1 | info.size>5

		error('size should be a small integer (1-5)');

	end

	if info.period<0

		error('period should not be negative');

	end

	if info.clearance<0 | info.clearance>1

		error('clearance should be between 0 and 1');

	end

	

	spars=get(0,'screensize');

	sl=spars(1);

	sb=spars(2);

	sw=spars(3);

	sh=spars(4);

	

	pw=200*info.size;			% progress bar width

	ph=20*info.size;			% progress bar height

	mh=~isempty(info.msg)*30;	% message bar height

	th=ph+mh;					% total height

	

	% check position

	switch info.pos

	case {'center','centre'}

		x=sl+sw/2;

		y=sb+sh/2;

	case {'centerleft','centreleft'}

		x=sl+sw/4;

		y=sb+sh/2;

	case {'centerright','centreright'}

		x=sl+3*sw/4;

		y=sb+sh/2;

	case {'top'}

		x=sl+sw/2;

		y=sb+3*sh/4;

	case {'topleft'}

		x=sl+sw/4;

		y=sb+3*sh/4;

	case {'topright'}

		x=sl+3*sw/4;

		y=sb+3*sh/4;

	case {'bottom'}

		x=sl+sw/2;

		y=sb+sh/4;

	case {'bottomleft'}

		x=sl+sw/4;

		y=sb+sh/4;

	case {'bottomright'}

		x=sl+3*sw/4;

		y=sb+sh/4;

	otherwise

		error(['pos ''' info.pos ''' was not recognised']);

	end

	

	info.lastclock=clock;

	info.id='progbar';

	handle=figure(...

	'MenuBar','none',...

	'numbertitle','off',...

	'userdata',info,...

	'name',info.title);

	set(handle,'position',[x-pw/2 y-th/2 pw th]);

	set(handle,'resize','off');

	

	% get color

	if ~isempty(which('rgb.m')) & 0

		color=rgb(info.color);

	else

		color=info.color;

	end

	

	% check color

	if ~isa(color,'double') | size(color,1)~=1 | size(color,2)~=3

		close(handle);

		error('Unrecognised color');

	end

	

	rectangle('position',[0 0 0.001 1],'facecolor',color,'edgecolor',color,'visible','off')

	set(gca,'position',[0 0 1 ph/(ph+mh)])

	axis([0 100 -info.clearance 1+info.clearance])

	set(gca,'Xtick',0:10:100);

	set(gca,'Ytick',[-1 2]);

	set(gca,'box','on');

	

	title(info.msg);

	drawnow

elseif nargin==2

	pb=varargin{1};

	progress=varargin{2};

	

	% check types

	if ~isnumeric(pb)

		error('pb should be a progress bar handle');

	end

	if ~isnumeric(progress)

		error('progress should be numeric');

	end

	if ~ishandle(pb) return; end % ignore if closed already

	

	% close

	if(progress==-1) close(pb); return; end



	% update

	info=get(pb,'userdata');

	if (info.period>0)

		if etime(clock,info.lastclock)<info.period return; end

		info.lastclock=clock;

		set(pb,'userdata',info);

	end

	

	if(progress<=0) progress=0.001; end

	if(progress>100) progress=100; end

	

	BarHandle=get(get(pb,'Children'),'Children');

	set(BarHandle,'position',[0 0 progress 1],'visible','on');

	progress=ceil(progress-0.5);

	set(pb,'name',[int2str(progress) '% - ' info.title]);

	drawnow

elseif nargin==3

	pb=varargin{1};

	progress=varargin{2};

	msg=varargin{3};

	

	% check types

	if ~isnumeric(pb)

		error('pb should be a progress bar handle');

	end

	if ~isnumeric(progress)

		error('progress should be numeric');

	end

	if ~ischar(msg)

		error('msg should be a string');

	end

	if ~ishandle(pb) return; end % ignore if closed already

	

	% close

	if(progress==-1) close(pb); return; end



	% update

	info=get(pb,'userdata');

	if (info.period>0)

		if etime(clock,info.lastclock)<info.period return; end

		info.lastclock=clock;

		set(pb,'userdata',info);

	end

	

	if(progress<=0) progress=0.001; end

	if(progress>100) progress=100; end

	

	BarHandle=get(get(pb,'Children'),'Children');

	set(BarHandle,'position',[0 0 progress 1],'visible','on');

	progress=ceil(progress-0.5);

	set(pb,'name',[int2str(progress) '% - ' info.title]);

	

	ax=get(pb,'children');

	ti=get(ax,'title');

	set(ti,'string',msg);

	

	drawnow

else

	error('Wrong number of input args');

end