function  [E, B, Phi] = Landau_Mathieu(Z, N, L, Neig, nu, lambda, ...
                                        sigma, mu, omega)
% Computing the spectrum of the Schroedinger operator H = (p-eA/c)^2
% on a 2 dimensional flat torus with boundary conditions dictated by
% gauge symmetry.  The magnetic field is given by 
%        B (1+lambda( sin(2*pi*nu* x/L)+sigma sin(4*pi*nu* x/L))
% Dirac condition requires B L^2/(2 pi) = Z = integer.
% This is automatically taken into account by deriving B as a
% function of Z=monopole charge=integer.
% The potential is periodic (differentiable on the torus)
% Momentum variables are computed via FFT.
% Units: hbar = 1, e/c = 1, mass=1
% Torus has equal sides of length L=sqrt(2);
%
% Usage: [E, B, Phi] = Landau_Mathieu(Z, N, L, Neig, nu, lambda,
%  sigma, mu, omega)
%
% All inputs have a default value:
%   mu:             % periodicity of scalar pot.
%   omega: 0.0      % scalar potential coupling
%   sigma: 0.5      % second harmonic of B
%   lambda: 0.5     % fundamental harmonic of B (0.5 means 50%)
%   nu:    4.0      % # of periods of B perturbation
%   Neig: 32        % number of eig.v's
%   L:  sqrt(2)
%   N:    64        % number of grid points in each direction
%   Z:    4         % magnetic charge (flux/2/pi)

%% Load parameters from command line or apply defaults
if nargin<9,  omega=0;     end
if nargin<8,  mu=4;        end
if nargin<7,  sigma=0.5;   end       
if nargin<6,  lambda=0.5;  end       
if nargin<5,  nu=4;        end           
if nargin<4,  Neig=32;     end        
if nargin<3,  L=sqrt(2);     end        
if nargin<2,  N=64;        end        
if nargin<1,  Z=4;         end

B = 2*Z*pi/4/L^2;
lambda=lambda*B;   % input lambda is relative to B

%% -------- build the space Lattice ----------------------
dx     = 2*L/N;
[x,y]  = meshgrid(-L*(1-1/N) : dx : L*(1-1/N)); 

V      = 0.5*omega.^2*(sin(mu*pi*x/L)+...
         sin(2*mu*pi*x/L)/4);
                                         
%% -------- and the momentum lattice --------------------
[kx,ky] = meshgrid(fftshift(pi*(-N/2:(N/2-1))/L));

%% ----------------- GAUGE PHASE FACTOR
chi =- (B*x-L/nu/pi*lambda*(cos(nu*pi*x/L)+sigma*cos(2*nu*pi*x/L)/2)).*y;
gf  = exp(i*chi);             
gfb = exp(-i*chi);             

%%  initial wavefunction

%opts.v0 = init(x/(1-1/N), y/(1-1/N), 2*L, B, N, lambda, nu);


%%------------------ Options for eigs

opts.disp   = 1;    % eigs works quietly - alternatively =1 , =2.
opts.isreal = false;
opts.issym  = true;

%%----------------- Now compute eig.fcts and eig.vs

[Psi, E, flag] = eigs(@Hamiltonian, N^2, Neig, 'SR', opts);
[E,ind] = sort(diag(E));
E       = E+B/2;                        % restore zero p. energy

for j=1:Neig
  Phi{j} = reshape(Psi(:,ind(j)),N,N);
end

if flag
  disp(['error flag=',num2str(flag)]);
  disp('see eigs documentation');
end

subplot(211)
plot(E/B,'.--','Markersize',20)
ylabel('E/B','FontSize',20,'Rotation',0)
%grid
subplot(212)
plot(E(1:min(Z,Neig))/B,'.--','Markersize',20)
xlabel('n','FontSize',20)
ylabel('E/B','FontSize',20,'Rotation',0)
%grid

%          - CAUTION - 
% Embedded function - all variables
% above are seen by the subroutine
function xout=Hamiltonian(xin)

psi = reshape(xin, N, N);  % eigs works with vectors!

% Use fft instead of FFT2! 

psix  = ifft(kx.*fft(psi,N,2), N,2) - (B+lambda*(sin(nu*pi*x/L) + ...
                               sigma*sin(2*nu*pi*x/L))).*y.*psi;
psixx = ifft(kx.*fft(psix, N,2), N,2) - (B+lambda*(sin(nu*pi*x/L)+ ...
                               sigma*sin(2*nu*pi*x/L))).*y.*psix;

psiy  = gfb .* ifft(ky.*fft(gf.*psi,  N,1), N,1)+...
        (B*x-lambda*L/nu/pi*(cos(nu*pi*x/L)+sigma*cos(2*nu*pi*x/L)/2)).*psi;
psiyy = gfb .* ifft(ky.*fft(gf.*psiy, N,1), N,1)+...
        (B*x-lambda*L/nu/pi*(cos(nu*pi*x/L)+...
                             sigma*cos(2*nu*pi*x/L)/2)).*psiy;

Hpsi   = V.*psi + (psixx+psiyy)/2.0 - B/2*psi;

xout = reshape(Hpsi, N^2, 1);   
end
end

%%------ build initial wf with correct b.c.
function psi = init(x, y, L, B, N, lambda, nu)

psi = exp(-B/2*(x.^2 + y.^2));

for n1=--4:4
  for n2=-4:4
    psi = psi + exp(-B/2*((x+n1*L).^2 + (y+n2*L).^2)) .*...
          exp(i*L*n2*(B*x-lambda*L/2/nu/pi*cos(2*nu*pi*x/L)));
  end
end

psi = psi - exp(-B/2*(x.^2 + y.^2));
psi=reshape(psi,N^2,1)/norm(psi);

end
