function [] = wavepack_expan(idos,ff,let_len,x_0,width0,tf,k0,f_name_pref)
% wavepack_expan - calculate wave packet (WP) expansion in a 1D complex system
%
% Usage: [] = wavepack_expan(idos,ff,let_len,x_0,width0,tf,k0,f_name_pref);
%
% Input:
%   idos - integrated density of states (use output of 'calc_spectr_Fib_PL3.m')
%   ff   - list of eigenmodes (use output of 'calc_spectr_Fib_PL3.m')
%   let_len - a single letter length: the # of the numerical lattice points per letter
%   x_0   - (a scalar) center of initial wave packet (in units of letters)
%   width0 - width of initial wave packet (= twice the STD) in units of letters 
%            [multiple values can be given as a vector of max length=5];
%   tf=10^6.5 - final time in units of meV^(-1) {default: tf=10^6.5}
%   k0 - a row-array of k-vectors of the initial WP modulation [in units of 2*pi/letter] {0}
%   f_name_pref - preffix of the names of the saved figure files {'WPE_Fib_<date number>'}
%   NOTE: the following internal parameters exist:
%       plot_WP_spectroscopy - flag(0/1) to plot initial WP real space-energy map
%       dx=let_len/2; - the return probability is averaged over 2*dx+1 sites
%       animation_flag - flag(0/1) to plot animation of the WP expansion
%       close_flg - flag (0/1) to close the figure after saving them
% Output:
%   none (the results are saved in the figures)
%
% Detailed Description:
% =====================
%   1. Initial WP is defined by a real-valued Gaussian with a given center
% and width (x_0 and width_0). Moreover, to shift the mean energy upwards,
% modulation by sum of exp(i*k0*x)) can be added, where k0 may have multiple values.
%   2. To solve for the time evolution, initial WP is then decomposed into
% the eigenmodes provided in the input (ff).
%   3. The time evolution of the state is calculated using the above spectral
% decomposition and the energies provided in the input (as idos)
%   4. The calculated features are the (mean) return probability, the RMS
% displacement and the mean participation ratio.
%   5. The results of the calculation are plotted in three figures, which
% are saved in the current directory (presently, the spread and the part. ratio are
% plotted in same figure).
%   6. Steps 2-4 can be repeated for up to 5 different initial WP widthes
% (this limit is just the number of the plot colors, and can be modified in the
% code), and the results for different widths are plotted in same figures.

if ~exist('f_name_pref','var') || isempty(f_name_pref);
    t=clock;[num2str(t(2)) num2str(t(3)) num2str(t(4)) num2str(t(5))];
    f_name_pref=['WPE_Fib_' [num2str(t(2)) num2str(t(3)) num2str(t(4)) num2str(t(5))]]
end
if ~exist('k0','var') || isempty(k0);k0=0; end
if ~exist('tf','var') || isempty(tf);tf=10^6.5; end % default final time in meV^(-1)
len=size(ff,2); % deduce the system length from the # of point in the eigenmodes

%=============================================================================
% Initial parameters:
plot_WP_spectroscopy=0; % flag(0/1) to plot initial WP real space-energy map
dx=let_len/2; % the return probability is averaged over 2*dx+1 sites
animation_flag=1; % flag(0/1) to plot animation of the WP expansion (0/1)
close_flg=0; % flag (0/1) to close the figure after saving them
%=============================================================================

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [idos,ff,len]=calc_spectr_Fib_PL3(1,10); % calculcate idos (when not using input)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Create required graphical objects:
%===================================
% Plot IDOS and the (later) WP spectral decomposition coefficients:
h_ocp=figure; % Open figure and get handle
plot(idos(1,:),idos(2,:),'k');
xlabel('E [meV]','Interpreter','LaTex','FontSize',14)
title('Coefficients of the WP spectral decomposition and IDOS')
hold on
mrks={'og','sr','.b','*m','dk'};
clrs='brgmk';
h_rp=figure; % Open figure and get handle
% h_pr=figure; % Open figure and get handle
h_dx=figure; % Open figure and get handle
if animation_flag; h_anim=figure; end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     Calculate (1) the WP spectral decomposition and (2) WP time Evolution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for kk=1:length(width0), % Loop through different intial widths
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %         (1) Calculate the WP spectral decomposition
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    w0=4*(width0(kk)*let_len); %4*80; %20*16*3 % initial width
    x0=round(x_0*let_len); %1420; %1580; %round(len/2); 9300/4; % initial center
    if x0>len; error('x0 exceeds system length');end
%     k0=0; %pi*6/160; % initial wave vector
    x=[1:len]; % create x-coordinate vector
    
    % Create initial wavepacket (function f0)
    f0=zeros(1,len); % allocation
    f0(x0)=1;f0=filter2(gausswin(w0,4)', f0 ); % Build a real-valued Gaussian
    if any(k0~=0),
        f0=f0.*sum(exp(-1i*kron(x,k0'/let_len))) ; % Add harmonic modulation as a sum of plane waves with equal weight
    end
    f0=f0/sqrt(f0*f0'); % normalize to unity
    
    % f0(81:end)=0;
    % Calculate coefficients of spectral expansion of the initial WP
    c=f0*ff';
%     for ii=1:length(idos(1,2:2:end)); % Loop over the eigenfunction subspace [REMOVE THIS LOOP!]
%         c(ii)=f0*ff(ii,:)';
%     end
    
    % Plot the expansion coefficients on top of the (plotted previously) IDOS
    figure(h_ocp);
    h=semilogy(idos(1,2:2:end),abs(c).^2*1000,mrks{kk});
    legend(h,['width_0=' num2str(width0(kk))]);
    pause(0.001)
    
    if plot_WP_spectroscopy,
        % Visualize occupation in space-energy plane:
        figure
        [map,ens] = eigmode_map(ff(1:180,:).*repmat(c(1:180)',1,len),idos(1,2:2:360),3*10^(-4)/1);
        imagesc([1:len]/20*1.35,ens/(pi/2/20)^2,map); axis xy
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %         (2) Calculate WP Evoltion
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % Define time vector with the time points spaced in the logartihmic scale
    tt=[0 exp([[-2:log(1.04):0],[log(1.001):log(1.001):log(tf)]])];
    dt=diff([tt]); % time increaments (used for time averaging below)
    
    % Allocation:
    px0=zeros(1,length(tt));xm=px0;var_x=px0;ipr=px0;
    count=0;
    for ii=1:length(tt); % Loop through the time points and calc wavefunction for them
        % WP wavefunction at the current time tt(ii)
        f_t=ff'*(exp(-1i*tt(ii)*idos(1,2:2:end)').*c'); 
        f_t_int=f_t.*conj(f_t); % WP wavefunction intensity
        f_t_int=f_t_int/sum(f_t_int); % correct the normalizaion to be sure
        
        % Calculate the required WP features:
        %====================================
        % The return probability smeared over 2*dx+1 sites:
        px0(ii)=sum(f_t_int([x0-dx:x0+dx])); %f_t_int(x0); 
        % The Spatial spread variance:
        xm(ii)=sum(x*f_t_int); % mean coordinate
        var_x(ii)=sum((x-x0).^2*f_t_int); %sum((x-xm(ii)).^2.*f_t_int); % spatial spread
        ipr(ii)=sum(f_t_int.^2); % Inverse participation ratio (IPR)
        count=count+1;
        
        % Plot animation of the WP expansion
        if animation_flag & count==100,
            figure(h_anim)
            title(['tt=' num2str(tt(ii))])
            pause(0.00000001);plot(f_t_int);pause(0.0000001)
            count=0;
        end
    end
    
    % Calculate time-average of the Return probability
    rpx0=[1 cumsum([px0(1:end-1)+px0(2:end)].*dt/2)./tt(2:end)/px0(1)]; 
    % Calculate time-average of the Participation Ration (NOT IPR):
    mpr=[1/ipr(1) cumsum([1./ipr(1:end-1)+1./ipr(2:end)].*dt/2)./tt(2:end)];
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %                PLOTTING
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %     Plot time-averaged return probability
    figure(h_rp);
    h=loglog(tt,[rpx0',px0'/px0(1)],clrs(kk));
    set(h(1),'LineWidth',2)
    legend(h,{['C(t), w_0=' num2str(width0(kk))],...
        ['p(t), w_0=' num2str(width0(kk))]});
    
    xlabel('$t [meV^{-1}]$','Interpreter','LaTex','FontSize',14)
    ylabel('$p(t)$ and $C(t)$','Interpreter','LaTex','FontSize',14)
    hold on
    
    %     Plot spatial spread and participation ratio in same plot
    figure(h_dx);
   % h=loglog(tt',[var_x'.^0.5 mpr'],clrs(kk));
    h=loglog(tt',[var_x' mpr'],clrs(kk));
    legend(h,{['\Delta x, w_0=' num2str(width0(kk))],['PR, w_0=' num2str(width0(kk))]});
    
    xlabel('$t [meV^{-1}]$','Interpreter','LaTex','FontSize',14)
    ylabel('$\sqrt{<(x(t)-x_0)^2>}$ and PR','Interpreter','LaTex','FontSize',14)
    set(h(2),'LineWidth',2)
    hold on
    
    
end

% Save figures:
figure(h_ocp);legend('-DynamicLegend');
saveas(h_ocp,[f_name_pref '_c_x0_' strrep(num2str(x_0),'.','p')],'fig');

figure(h_rp);legend('-DynamicLegend');
saveas(h_rp,[f_name_pref '_RetPr_x0_' strrep(num2str(x_0),'.','p')],'fig');

figure(h_dx);legend('-DynamicLegend');
saveas(h_dx,[f_name_pref '_PartRat_x0_' strrep(num2str(x_0),'.','p')],'fig');
if close_flg,
    close(h_ocp)
    close(h_rp)
    close(h_dx)
end