function [idos,wf,len]=calc_spectr_Fib_2D(V,n,eigmod_flg);
% calc_spectr_Fib_2D - calculate spectrum and states of quasiperiodic 2D strip
% Usage: [idos,wf,len]=calc_spectr_Fib_2D(V,n,eigmod_flg);
% Input:
%   V - forced scaling of amplitude of the binary potential {default: V=1}
%   n - generation number of the Fibonacci sequence
%   eigmod_flg - eigenmode calculation method flag:
%                0 - eigenmode AMPLITUDE due to the 1st transverse mode is
%                calculated, or
%                1 - sqrt of the 1st TRANSVERSE mode eigenmode INTENSITIES including 
%                   evanescent contributions of the higher UNPERTURBED modes.   
%   NOTE: the following internal parameters exist
%       N_trvs - Number of the transverse modes
%       Use_Correction - flag to use the next mode correction in the 1D solution
%       smth_scale - width profile smoothness scale in inits of the letter length
%       P_grating - number of the numerical grid points per letter (i.e. the resolution)
%       n_ref  - Effective refraction index of the Bragg cavity
%       system_flg - (1,2,3,4) choice of the experimental system
% 
% Output:
%   idos - 2 x len array: 1st row - energies in units of the recoil energy Er
%                       2nd row - IDOS, the number of the states under the energy
%   wf - n x len array of the "eigenmodes" calculated due to the eigmod_flg
%        (see input description above)
%   len  - length of the system
% 
% Detailed description: 
% ====================
% (1) 2D calculcation:
% An exact diagonalization is applied to the trancated Hamiltonian
% matrix, calculated using N_trvs transverse modes with N_sub states in the lowest one. The basis is
% spanned by the product states of sines(cosines) with the transverse (y) component with the x-dependent
% wavevector and the longitudinal (x) component the eigenfunctions of the uniform system.
% 
% (2) 1D calculation:
% This is done by choosing the number of the transverse modes (N_trvs) equal to
% unity. Then, there is another internal parameter "Use_Correction" which
% define whether to apply the correction due to the coupling of the lowest
% mode to the next one. For Use_Correction=0 only the single-mode trancated
% Hamiltonanian is applied (i.e. H_00). For Use_Correction=1, an additional
% contribution (derived due to H_01*G_11*H_10) is added.
% 
% (3) The code below has traces of several different attempts to approximate
% the above correction. Some of them are still selectable (through if's). The last
% version is using the approximate but fully analytic (up to some tunable parameter
% of the order of unity) expression for the mean/typical Green's function G_11.
% Former attempts employed some phenomenological approximations (see the document).
% 
% (4) The polaritonic system parameter are choosable from four diffent
% experimental configurations (Fib2, Fib18 (x2) and Fib15 , - see
% Dimitrii's resumes). The choice is done using the internal parameter
% "system_flg". For each configuration, the corect values of w_a,w_b,LL (letter widths 
% and length) and other parameters of the polaritonic system like Ex, Ec etc.
% are taken(see the code below).
% 
% (5) Other internally set parameters exist (see input description above):
%
% Info: this routine was last modified 4/02/2014 (seems to be the last
%       version for the polariton spectrum calculation)

% Calculational and phenomenological parameters:
%===============================================
% Number of the subspace wave functions for zero transverse mode (basis size)
N_sub=round(2*400);
% Number of the transverse modes
N_trvs=1
if N_trvs==1, % Define flag to use the next mode correction in the 1D solution (only for single-mode)
    Use_Correction=0
else
    Use_Correction=0
end
P_grating=80*(4/8); % Grating period (letter length) in units of TBM lattice (site) spacings (affects numerical accuracy)
en_shift=0; %2.2; % manual energy shift of the 1st band (meV) (to tune/correct the spectral edge)
smth_scale=0.5*(1.35/0.8)*1.07/8; %0.22*1.07/8; %1.3*1.07/8; % Smoothness scale in inits of the letter length (1.07/8 - to match the experiment)


% Experimental Parameters:
%===========================
% Effective refraction index of the Bragg cavity (value recommended by Dimitrii was 3.12-3.15)
n_ref=3.25+0.08; %3.15;  

system_flg=1;
if system_flg==1, % System "Fib2" used initially in the EPL paper
    disp('System "Fib2" used initially in the EPL paper')
    Ex=1581.7+1*0.4;  % Exciton energy E_x [in meV]
    Ec=1575.7-0.5*1;  % Photon 1st band energy (modified from usual "ground" energy!) E_c [in meV]
    E_b1=1570.8;      % Energy of the 1st polariton band (in the experiment) [meV]
    w_a=3.5;  % Transverse width of the A-letter (mkm)
    w_b=2.04; % Transverse width of the A-letter (mkm)
    LL=1.35;  % Letter length(mkm)
elseif system_flg==2; % System "Wire_18" described in the "French" PDF memo describing first measurements
    disp('System "Wire_18" described in the first PDF memo')
    Ex=1607.5; % Exciton energy E_x [in meV]
    Ec=Ex-6; % Photon 1st band energy (modified from usual "ground" energy!) E_c [in meV]
    E_b1=1596.5; % Energy of the 1st polariton band (in the experiment) [meV]
    w_a=3.5;  % Transverse width of the A-letter (mkm)
    w_b=1.86; % Transverse width of the A-letter (mkm)
    LL=0.8;  % Letter length(mkm)
elseif system_flg==3; % System "Wire_18" chosen as the 1st system for the EPL paper v2
    disp('System "Wire_18" chosen as the 1st system for the EPL paper v2')
    Ex=1605; % Exciton energy E_x [in meV]
    Ec=Ex-8.3-0.0; % Photon 1st band energy (modified from usual "ground" energy!) E_c [in meV]
    E_b1=1592.3; % Energy of the 1st polariton band (in the experiment) [meV]
    w_a=3.5;  % Transverse width of the A-letter (mkm)
    w_b=1.86; % Transverse width of the A-letter (mkm)
    LL=0.8;  % Letter length(mkm)
elseif system_flg==4; % System "Wire_15"
    disp('System "Wire_15"')
    Ex=1607; % Exciton energy E_x [in meV]
    Ec=Ex-8.7-0.0; % Photon 1st band energy (modified from usual "ground" energy!) E_c [in meV]
    E_b1=1594; % Energy of the 1st polariton band (in the experiment) [meV]
    w_a=3.5;  % Transverse width of the A-letter (mkm)
    w_b=1.86; % Transverse width of the A-letter (mkm)
    LL=0.5;  % Letter length(mkm)
end

%  Derived parameters:
%============================

% Define recoil energy for photons (not polaritones) [meV]:
E_r_ph=0.5*(1.054*10^(-34))^2*(3*10^8/n_ref)^2*(pi/2/LL*10^6)^2/Ec/(1.6*10^(-19))^2*10^6; 

% link between the Rabi frequency and the Ec energy

E0=Ex-E_b1; % Difference between the excitonic level and the 1st polaritonic band
detuning_0=Ec-Ex; % Detuning at k=0
% Calculate Rabi frequency (exciton-photon coupling), which is not independent of
% Ex and Ec, once E_b1 is fixed:
Omega=-detuning_0*sqrt((1+2*E0/detuning_0)^2-1) % Rabi frequency [meV]
%===================

Hopf_c0=(abs(detuning_0)/sqrt(1*detuning_0^2+Omega^2)+1)/2 % Hopfeld coef at k=0

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                   Define potential
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Letter length (half grating period) in units of TBM lattice (site) spacing.
% The "recoil energy" is E_r=pi^2*hbar^2/(8*M*P_grating^2), and we use units where
% hbar^2/(2*M)=1.

% Express letter widths and length in the units of the numerical lattice spacing:
w_a=w_a/LL*P_grating;
w_b=w_b/LL*P_grating;
LL=P_grating;

% Express energy parameters and potential amplitude in units of hbar^2/(2*M):
Er=(pi^2/P_grating^2/4);
%V=V*Er; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                           Calculate photonic spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Generate quasiperiodic sequence and width profile
% =================================================
seq=[];
for n_gen=n:n,
   seq=[seq gen_fib_seq(n_gen)]; % the letters are mapped as A=1, B=-1
end

if system_flg==1,
    n_pad=7; % Number of the A-letters to pad each side (depends on the system profile)
else
    n_pad=0;
end
seq=[ones(1,n_pad) seq ones(1,n_pad)];
% seq=[repmat([1 -1],1,77) 1]; % periodic

w=w_a+(w_b-w_a)*(1-seq)/2; % map binary code sequence into width profile sequence

w=repmat(w,P_grating,1); % replicate width profile using P_grating points per letter
w=w(:)';

% Smooth width profile according to "smth_scale":
% =====================
ker=gausswin(P_grating,1/smth_scale)'; % smoothing convolution kernel
% Alternative smoothing kernel (half-circle) - used to check generality
% ker=(smth_scale*P_grating).^2-([1:2*P_grating]-P_grating).^2; ker(ker<0)=[];ker=sqrt(ker);

ker=ker/sum(ker); % kernel integral is normalized to unity
% Do Convolution smoothing of the wire width profile:
w=filter2(ker, [fliplr(w(1:5*P_grating)) w fliplr(w(end-5*P_grating+1:end))]);
% Correct (discard) width in the convolution-invalid regions of the "paddings"    
w(1:5*P_grating)=[];w(end-5*P_grating+1:end)=[];
% w(1:P_grating)=w_a; w(end-P_grating:end)=w_a; % correct potential in the "paddings"

% Calculcate the required derivatives of the wire width:
% =====================================================
Dw=([0 diff(w)]+[diff(w) 0])/2; % 1st derivative of the width
DDw=([0 diff(Dw)]+[diff(Dw) 0])/2; % 2nd derivative of the width
DDDw=([0 diff(DDw)]+[diff(DDw) 0])/2; % 3rd derivative of the width

% Prepare Effective potential ingredients:
% =======================================
inv_w2=1./w.^2;
DwOw=Dw./w;
DwOw2=DwOw.^2;
DDwOw=DDw./w;
% DDDwOw=DDDw./w; % not used

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       Calculate the Hamiltonian matrix in the chosen sub-space
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

len=length(w); % System length
x=pi*[1:len]/(2+len); % x-coordinate vector

% Calculate required dimensions of the higher transverse mode subspaces 
% (truncating the subspaces at same cut-off energy):
Ns2=N_sub^2-(len/w_b)^2*((1+2*[0:(N_trvs-1)]).^2-1); % squared dimension per tranverse mode
% Check modes higher in energy then max energy in zero transverse mode subspace:
if any(Ns2<=0), % reduce the # of transv. modes if some Ns2 are negative (those modes are too high in energy)
    Ns2(Ns2<=0)=[];
    Ns=round(sqrt(Ns2));
    N_trvs=length(Ns);
    % Inform about reducing the number of the transverse modes
    disp(['!!! The number of the transverse modes has been reduced to N_trvs=' num2str(N_trvs)]);
else
    Ns=round(sqrt(Ns2));
end
% Total size of the subspace (i.e. including all the transverse modes used):
N_tot=sum(Ns);
h_bl_sizes=cumsum([0 Ns]); % CUMULATIVE sizes of the transverse mode subspaces

% Allocate the total Hamiltonian matrix
h=zeros(N_tot,N_tot); 

% Construct basis functions in x-direction of the size of the ZERO mode only (which is largest):
% ==========================================================================
base=zeros(len,N_sub); % allocate for the basis functions
dbase=zeros(len,N_sub); % allocate for the derivatives of the basis functions
for ii=1:N_sub,
    base(:,ii)=sin(ii*x);
    dbase(:,ii)=(ii*pi/(len+2))*cos(ii*x); % eigenfunction derivative (includes the k-factor from the 1st der)
end
base=base/sqrt((len+2)/2); % normalize
dbase=dbase/sqrt((len+2)/2); % normalize

% Calculate the off-diagonal (due the effective potential) part of the Hamiltonian
% ==============================================================================

for q=1:N_trvs, % Double-Loop through the transverse modes and calculate Hamiltonian blocks
    for n=1:N_trvs,
        if 0, % Same as below in "else", but due to the non-symmetric calc of the matrix elements
            IQN=iqn(q-1,n-1); IYQN=iyqn(q-1,n-1); IYYQN=iyyqn(q-1,n-1);
            z=(1+2*(n-1));
            v1=IQN*pi^2*z^2*inv_w2+...
                (-3/4*IQN+z*pi*(3*IYQN+z*pi*IYYQN))*DwOw2+...
                (IQN-2*IYQN*z*pi)/2*DDwOw;
            v2=(IQN-2*IYQN*z*pi)*DwOw;
            
            h((h_bl_sizes(q)+1):h_bl_sizes(q+1),(h_bl_sizes(n)+1):h_bl_sizes(n+1))=...
                base(:,1:Ns(q))'*...
                (repmat(v1',1,Ns(n)).*base(:,1:Ns(n))+repmat(v2',1,Ns(n)).*dbase(:,1:Ns(n)));
        else % Due to symmetric calculation of the matrix elements
            % Construct terms in the potential:
            H0=h0(q-1,n-1); H1=h1(q-1,n-1); H2=h2(q-1,n-1);
            v1=H0*inv_w2 + H2*DwOw2;
            v2=H1*DwOw;
            
 
            % Calculate the Hamiltonian matrix
            %=================================
            h((h_bl_sizes(q)+1):h_bl_sizes(q+1),(h_bl_sizes(n)+1):h_bl_sizes(n+1))=...
                base(:,1:Ns(q))'*...
                (repmat(v1',1,Ns(n)).*base(:,1:Ns(n))+repmat(v2',1,Ns(n)).*dbase(:,1:Ns(n)))-...
                dbase(:,1:Ns(q))'*(repmat(v2',1,Ns(n)).*base(:,1:Ns(n)));
            
            if Use_Correction, % Define correction terms due to the next-mode-coupling in 1D calculation
                non_local_correction=1
                if non_local_correction, % Apply non-local Green's function correction
                    % Construct the approximate (or effective) Green's function (see document)
                    if 0, % Obsolete - woring, but based on wrongly (typo) calculated theoretical coefficients
                        E11=(11+(3+9*pi^2)/3/pi^2*(1.5/2.7)^2/(0.07*(smth_scale/(0.25*1.07/8)))^2)*Er;
                        % smoothing convolution kernel (effective Green's function (???)):
                        g_ker=(4.75/6.83/sqrt(E11))*exp(-abs([0:6*P_grating]-3*P_grating)*sqrt(E11)*0.9); 
                        %                     g_ker=g_ker/sum(g_ker); % kernel integral is normalized to unity
    
                    else % based on approx but fully analytical construction of the meam Green's function
                        beta=1.05*sqrt(2*pi)*smth_scale/sqrt(2)/2^(1/4); % width of the rect. barrier
                        m=1; % mode index (=1 for the second mode)
                        U=(4*(1+2*m)^2*(2*LL/(w_a+w_b))^2 - 4*(LL/w_a)^2 +...
                            (3+(1+2*m)^2*pi^2)/3/pi^2*(2*(w_a-w_b)/(w_a+w_b))^2/beta^2)*Er;
                        k=sqrt(U);
                        kx0=k*LL*beta;
                        g0=(cosh(kx0)+sinh(kx0)/kx0)/(2*k*sinh(kx0)); % value of the Green's func. at x1=x2
                        % Exponent (it is different from k !):
                        kappa=k*(kx0^2*sinh(kx0)+kx0*cosh(kx0)-sinh(kx0))/(kx0*cosh(kx0)+sinh(kx0))/kx0;
                        
                        g_ker=g0*exp(-abs([0:6*P_grating]-3*P_grating)*kappa); % smoothing convolution kernel
                    end
                    
                    F=(7/2)*DwOw2-DDwOw;
                    Q=-2*DwOw;
                    
                    % Convolution smoothing:
                    GF=filter2(g_ker, F );
                    GQ=filter2(g_ker, Q );
                    % correct width in the convolution-invalid regions of the "paddings"
                    GF(1:P_grating)=F(1:P_grating); GF(end-P_grating:end)=F(end-P_grating:end); % correct potential in the "paddings"
                    GQ(1:P_grating)=Q(1:P_grating); GQ(end-P_grating:end)=Q(end-P_grating:end); % correct potential in the "paddings"
                    
                    ket=-(3/4)*((repmat(GF',1,Ns(n)).*base(:,1:Ns(n))) + (repmat(GQ',1,Ns(n)).*dbase(:,1:Ns(n))));
                    % Alternatively, apply Green func to the wavefucntion as well
                    %    ket=-(3/4)*((repmat(F',1,Ns(n)).*base(:,1:Ns(n))) + (repmat(Q',1,Ns(n)).*dbase(:,1:Ns(n))));
                    %    ket=filter2(g_ker', ket );
                    
                    bra=-(3/4)*((repmat(F',1,Ns(n)).*base(:,1:Ns(n))) + (repmat(Q',1,Ns(n)).*dbase(:,1:Ns(n))));
                    
                    % Correction to the Hamiltonain matrix due to the coupling to the 2nd mode
                    h_corr=-bra'*ket;
                    
                else % local correction (using delta function approximation for the Green's function)
                    % Energy to put into the denominator (this is an empirical formula):
                    E11=(0.6/smth_scale)*(11+19*((0.25*1.07/8)/smth_scale)^2)*Er %(10+20*(0.25*1.07/8)/smth_scale)*Er;
                    mass=(9/4/E11)*DwOw2;
                    v010=(9/16/E11)*(35/4*DwOw2.^2 - 18*DDwOw.*DwOw2 + DDwOw.^2 + 2*DwOw.*DDDwOw);
                    % Correction to the Hamiltonain matrix due to the coupling to the 2nd mode
                    h_corr=-dbase(:,1:Ns(q))'*(repmat(mass',1,Ns(n)).*dbase(:,1:Ns(n)))+...
                        base(:,1:Ns(q))'*(repmat(v010',1,Ns(n)).*base(:,1:Ns(n)));
                end
                % Add correction term to the Hamiltonain matrix:
                h((h_bl_sizes(q)+1):h_bl_sizes(q+1),(h_bl_sizes(n)+1):h_bl_sizes(n+1))=h_corr+...
                    h((h_bl_sizes(q)+1):h_bl_sizes(q+1),(h_bl_sizes(n)+1):h_bl_sizes(n+1));
            end
        end
        % % %         if q~=n, % if off-diagonal block, then add transpose to it
        % % %             h((h_bl_sizes(n)+1):h_bl_sizes(n+1),(h_bl_sizes(q)+1):h_bl_sizes(q+1))=...
        % % %                 h((h_bl_sizes(q)+1):h_bl_sizes(q+1),(h_bl_sizes(n)+1):h_bl_sizes(n+1));
        % % %         else
        % % %             h((h_bl_sizes(q)+1):h_bl_sizes(q+1),(h_bl_sizes(n)+1):h_bl_sizes(n+1))=...
        % % %                 (h((h_bl_sizes(q)+1):h_bl_sizes(q+1),(h_bl_sizes(n)+1):h_bl_sizes(n+1))+...
        % % %                 h((h_bl_sizes(q)+1):h_bl_sizes(q+1),(h_bl_sizes(n)+1):h_bl_sizes(n+1))')/2;
        % % %         end
    end
end

% Symmetrize the Hamiltonian matrix (otherwise the result could be unstable numerically)
% ============================================================
h=(h+h')/2; 

% Add The kinetic (unperturbed diagonal) part of the Hamiltonian
% ==============================================================
% h(:,:)=0;
for q=1:N_trvs,
    for ii=(h_bl_sizes(q)+1):h_bl_sizes(q+1),
        h_kin=(pi*(ii-h_bl_sizes(q))/(2+len))^2; % Kinetic energy
        h(ii,ii)=h(ii,ii)+h_kin;
    end
end

clear w Dw DDw DwOw DDwOw DwOw2 dbase inv_w2 v1 v2 x

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                   Diagonalize the Hamiltonian matrix:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[modes,eig_ens]=eig(h);

% Get eigenmods (or thier intensities) - this is same for photon or polaritons:
wf=zeros(size(base)); % (base*modes(1:Ns(1),1:Ns(1)))';
if eigmod_flg==0, % Calc eigenmode AMPLITUDEs due to the 1st transverse mode
    q=1;
    wf=wf+(base(:,1:Ns(q))*...
        modes((h_bl_sizes(q)+1):h_bl_sizes(q+1),1:Ns(1)));
else
    % Calculate sqrt of the 1st TRANSVERSE mode eigenmode INTENSITIES including evanescent 
    % contributions of higher tr. modes (intra-mode states sum in amplitude, while 
    % inter-mode states sum in intenities)
    for q=1:N_trvs,
        wf=wf+abs(base(:,1:Ns(q))*...
            modes((h_bl_sizes(q)+1):h_bl_sizes(q+1),1:Ns(1))).^2;
    end
    wf=sqrt(wf);
end
wf=wf';

% Get photonic eigenenergies:
% ==========================
eig_ens=diag(eig_ens);
% Shift photonic eigenenergies to get zero for the ground state:
if system_flg==1,
    en_shift= -en_shift*Er/E_r_ph+(80/P_grating)^2*4.4e-04; %min(eig_ens)
    eig_ens=eig_ens-en_shift;
else
    eig_ens=eig_ens-min(eig_ens);
end

% Return to the units of meV:
eig_ens=eig_ens*E_r_ph/Er;

% Calculate polariton energies by coupling photons to excitons:
% E_0=(detuning_0+sqrt(detuning_0^2+Omega^2))/2; % Energy offset to make E(k=0)=0
del=(eig_ens+detuning_0);
eig_ens=(Ex+(del-sqrt(del.^2+Omega^2))/2);


% Calculate IDOS (two points per state define step of the width "step_width")
step_width=10^(-13);
idos(1,:)=reshape([eig_ens'-step_width;eig_ens'],1,2*N_tot);
idos(2,:)=reshape([[0:N_tot-1];[1:N_tot]],1,2*N_tot);

return


% The subroutines below define coefficients in the effective potential
function x=iqn(q,n),
if q==n;
    x=1;
else x=0;
end

function x=iyqn(q,n),
if q==n;
    x=1/(2*pi*(1+2*n));
else
    x=-(-1)^(n+q)*(1+2*q)/(2*pi*(n-q)*(1+n+q));
end

function x=iyyqn(q,n),
if q==n;
    x=1/12-2/(2*pi*(1+2*n))^2;
else
    x=2*(-1)^(n+q)*(1+2*q)*(1+2*n)/(2*pi*(n-q)*(1+n+q))^2;
end

function x=h0(q,n),
if q==n;
    x=pi^2*(1+2*n)^2;
else
    x=0;
end

function x=h1(q,n),
if q==n;
    x=0;
else
    x=(-1)^(n+q)*(1+2*q)*(1+2*n)/(2*(n-q)*(1+n+q));
end

function x=h2(q,n),
if q==n;
    x=(3+pi^2*(1+2*n)^2)/12;
else
    x=(-1)^(n+q)*(1+2*q)*(1+2*n)*(1+2*n*(1+n)+2*q*(1+q))/2/((n-q)*(1+n+q))^2;
end