function schrsolve(potential,fl)

%SCHRSOLVE solve Schroedinger equation
%
%schrsolve(potential,flag)
%
%Solves Schroedingers equation for all quantum boxes and all carrier types.
%The resulting energies and wavefunctions are stored in the global
%structure aquila_subbands.xx (xx=ge,xe,le,hh,lh,so).
%The necessary information concerning material composition etc. are taken from the 
%global variables characterizing the grid.
%
%potential is the electrostatic potential for which Schroedingers equation is to be solved.
%flag=1 indicates, that old solutions should be adapted to a slightly changed
%   potential via inverse vectoriteration instead of a complete recomputation
%   of the energy spectrum.
%
%Result: aquila_subbands.xx(nr).E = vector of energies or matrix for 1D quantization
%        aquila_subbands.xx(nr).psi = wavefunctions of all subbands
%           size=nr. of solutions * size of corresponding QBOX
%           i.e. psi=[1st 2nd 3rd ..] subband-wavefunction
%        nr numbers the QBOXes in order of definition
%(See file 'structures' for a description of aquila_subbands).

%Copyright 1999 Martin Rother
%
%This file is part of AQUILA.
%
%AQUILA is free software; you can redistribute it and/or modify
%it under the terms of the GNU General Public License as published by
%the Free Software Foundation; either version 2 of the License, or
%(at your option) any later version.
%
%AQUILA is distributed in the hope that it will be useful,
%but WITHOUT ANY WARRANTY; without even the implied warranty of
%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%GNU General Public License for more details.
%
%You should have received a copy of the GNU General Public License
%along with AQUILA; if not, write to the Free Software
%Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

global aquila_subbands aquila_control aquila_structure aquila_material

%check for correct execution order
if bitget(aquila_control.progress_check,6)==0
   error('schrsolve: You must run BUILDSTRUCTURE before solving Schroedingers equation !')
end

%define some constants
constants

%tell the user, what we are doing
if aquila_control.verbose>0
   disp('schrsolve: solving Schroedinger equation for all quantum regions')               
end

n=size(aquila_structure.qbox);
for count=1:n(1) %for all QBOXes   
   if bitand(aquila_structure.qbox(count,9),GE)>0 %the QBOX works on gamma electrons
      %find the indices of the corresponding nodes in the grid
      [ix,iy]=boxindex(aquila_structure.xpos,aquila_structure.ypos,...
         aquila_structure.qbox(count,[1:4]));
      %get the masses
      mx=aquila_material.meg(iy,ix);
      my=mx; 
      %correct the indices
      ix=[ix ix(end)+1];
      if aquila_control.mode==2
         iy=[iy iy(end)+1];
      end
      %compute the potential = sum of conduction band edge and electric potential
      pot=aquila_material.eg(iy,ix)-potential(iy,ix);
      %some output for the user
      if aquila_control.verbose>1
         os=sprintf('qbox nr %d, Gamma electrons, type %d',count,aquila_structure.qbox(count,8));
         disp(os)               
      end
      %call the driver routine for the solver
      if length(aquila_subbands.ge)>=count %we already have subbands -> use them in the solver
         [E,ev]=solveit(pot,ix,iy,mx,my,aquila_structure.qbox(count,7),...
            aquila_structure.qbox(count,8),aquila_subbands.ge(count),fl);
      else %we have no subbands computed before
         [E,ev]=solveit(pot,ix,iy,mx,my,aquila_structure.qbox(count,7),...
            aquila_structure.qbox(count,8),[],fl);
      end
      %store the result
      aquila_subbands.ge(count).E=E;
      aquila_subbands.ge(count).psi=ev;
   end
   
   %here comes the same stuff for X-point electrons
   if bitand(aquila_structure.qbox(count,9),XE)>0
      [ix,iy]=boxindex(aquila_structure.xpos,aquila_structure.ypos,...
         aquila_structure.qbox(count,[1:4]));
      mx=aquila_material.mex(iy,ix);
      my=mx;
      ix=[ix ix(end)+1];
      if aquila_control.mode==2
         iy=[iy iy(end)+1];
      end
      pot=aquila_material.ex(iy,ix)-potential(iy,ix);
      if aquila_control.verbose>1
         os=sprintf('qbox nr %d, X electrons, type %d',count,aquila_structure.qbox(count,8));
         disp(os)               
      end
      if length(aquila_subbands.xe)>=count
         [E,ev]=solveit(pot,ix,iy,mx,my,aquila_structure.qbox(count,7),...
            aquila_structure.qbox(count,8),aquila_subbands.xe(count),fl);
      else
         [E,ev]=solveit(pot,ix,iy,mx,my,aquila_structure.qbox(count,7),...
            aquila_structure.qbox(count,8),[],fl);
      end
      aquila_subbands.xe(count).E=E;
      aquila_subbands.xe(count).psi=ev;
   end
   
   %L-point electrons
   if bitand(aquila_structure.qbox(count,9),LE)>0
      [ix,iy]=boxindex(aquila_structure.xpos,aquila_structure.ypos,...
         aquila_structure.qbox(count,[1:4]));
      mx=aquila_material.mel(iy,ix);
      my=mx;
      ix=[ix ix(end)+1];
      if aquila_control.mode==2
         iy=[iy iy(end)+1];
      end
      pot=aquila_material.el(iy,ix)-potential(iy,ix);
      if aquila_control.verbose>1
         os=sprintf('qbox nr %d, L electrons, type %d',count,aquila_structure.qbox(count,8));
         disp(os)               
      end
      if length(aquila_subbands.le)>=count
         [E,ev]=solveit(pot,ix,iy,mx,my,aquila_structure.qbox(count,7),...
            aquila_structure.qbox(count,8),aquila_subbands.le(count),fl);
      else
         [E,ev]=solveit(pot,ix,iy,mx,my,aquila_structure.qbox(count,7),...
            aquila_structure.qbox(count,8),[],fl);
      end
      aquila_subbands.le(count).E=E;
      aquila_subbands.le(count).psi=ev;
   end
   
   %light holes
   if bitand(aquila_structure.qbox(count,9),LH)>0
      [ix,iy]=boxindex(aquila_structure.xpos,aquila_structure.ypos,...
         aquila_structure.qbox(count,[1:4]));
      mx=aquila_material.mlh001(iy,ix);
      my=aquila_material.mlh110(iy,ix);
      ix=[ix ix(end)+1];
      if aquila_control.mode==2
         iy=[iy iy(end)+1];
      end
      pot=-(aquila_material.ev(iy,ix)-potential(iy,ix));
      if aquila_control.verbose>1
         os=sprintf('qbox nr %d, light holes, type %d',count,aquila_structure.qbox(count,8));
         disp(os)               
      end
      if length(aquila_subbands.lh)>=count
         [E,ev]=solveit(pot,ix,iy,mx,my,aquila_structure.qbox(count,7),...
            aquila_structure.qbox(count,8),aquila_subbands.lh(count),fl);
      else
         [E,ev]=solveit(pot,ix,iy,mx,my,aquila_structure.qbox(count,7),...
            aquila_structure.qbox(count,8),[],fl);
      end
      aquila_subbands.lh(count).E=-E;
      aquila_subbands.lh(count).psi=ev;
   end
   
   %heavy holes
   if bitand(aquila_structure.qbox(count,9),HH)>0
      [ix,iy]=boxindex(aquila_structure.xpos,aquila_structure.ypos,...
         aquila_structure.qbox(count,[1:4]));
      mx=aquila_material.mhh001(iy,ix);
      my=aquila_material.mhh110(iy,ix);
      ix=[ix ix(end)+1];
      if aquila_control.mode==2
         iy=[iy iy(end)+1];
      end
      pot=-(aquila_material.ev(iy,ix)-potential(iy,ix));
      if aquila_control.verbose>1
         os=sprintf('qbox nr %d, heavy holes, type %d',count,aquila_structure.qbox(count,8));
         disp(os)               
      end
      if length(aquila_subbands.hh)>=count
         [E,ev]=solveit(pot,ix,iy,mx,my,aquila_structure.qbox(count,7),...
            aquila_structure.qbox(count,8),aquila_subbands.hh(count),fl);
      else
         [E,ev]=solveit(pot,ix,iy,mx,my,aquila_structure.qbox(count,7),...
            aquila_structure.qbox(count,8),[],fl);
      end
      aquila_subbands.hh(count).E=-E;
      aquila_subbands.hh(count).psi=ev;
   end
   
   %split-off holes
   if bitand(aquila_structure.qbox(count,9),SO)>0
      [ix,iy]=boxindex(aquila_structure.xpos,aquila_structure.ypos,...
         aquila_structure.qbox(count,[1:4]));
      mx=aquila_material.mso(iy,ix);%DOS mass probably not correct here
      my=mx;
      ix=[ix ix(end)+1];
      if aquila_control.mode==2
         iy=[iy iy(end)+1];
      end
      pot=-(aquila_material.eso(iy,ix)-potential(iy,ix));
      if aquila_control.verbose>1
         os=sprintf('qbox nr %d, split-off holes, type %d',count,aquila_structure.qbox(count,8));
         disp(os)               
      end
      if length(aquila_subbands.so)>=count
         [E,ev]=solveit(pot,ix,iy,mx,my,aquila_structure.qbox(count,7),...
            aquila_structure.qbox(count,8),aquila_subbands.so(count),fl);
      else
         [E,ev]=solveit(pot,ix,iy,mx,my,aquila_structure.qbox(count,7),...
            aquila_structure.qbox(count,8),[],fl);
      end
      aquila_subbands.so(count).E=-E;
      aquila_subbands.so(count).psi=ev;
   end
end
%set the flag, that we now have subbands
aquila_control.progress_check=bitset(aquila_control.progress_check,7);   

%the driver routine for the solvers
function [E,ev]=solveit(pot,xp,yp,mx,my,subb,tp,subbands,fl);
global aquila_control aquila_structure
constants
%for every type of QBOX call the corresponding solver
%according to the flag 'fl' use Arnoldi (schrsolv..) or inverse vectoriteration (schrtrack..)
switch tp
case QWR %x and y quantized, a quantum wire
   if (bitget(aquila_control.progress_check,7)>0)&(fl==1)&(~isempty(subbands))
      [E,ev]=schrtrack2D(pot,xp,yp,mx,my,subbands.E,subbands.psi);
   else
      [E,ev]=schrsolv2D(pot,xp,yp,mx,my,subb);
   end
case QWX %y motion quantized, a quantum well in y-direction
   if (bitget(aquila_control.progress_check,7)>0)&(fl==1)&(~isempty(subbands))
      [E,ev]=schrtrack1D(pot,aquila_structure.xpos(xp),aquila_structure.ypos(yp),...
         my,subbands.E,subbands.psi);
   else
      [E,ev]=schrsolv1D(pot,aquila_structure.xpos(xp),aquila_structure.ypos(yp),my,subb);
   end
case QWY %x motion quantized, a quantum well in x-direction
   %this uses the same routine as in the QWX-case
   %so we must transpose the whole scenario to exchange x- and y-direction
   %and thereby froming a QWX scenario
   if (fl==1)&(~isempty(subbands))
      if aquila_control.mode==2
         psi=[];
         n=length(xp);
         for i_count=0:length(subbands.E(:,1))-1
            psi=[psi subbands.psi(:,i_count*n+1:(i_count+1)*n)'];
         end
      else
         psi=reshape(subbands.psi',length(xp),length(subbands.E));
      end
   end
   if (bitget(aquila_control.progress_check,7)>0)&(fl==1)&(~isempty(subbands))
      [E,ev]=schrtrack1D(pot',aquila_structure.ypos(yp),aquila_structure.xpos(xp),...
         mx',subbands.E,psi);
   else
      [E,ev]=schrsolv1D(pot',aquila_structure.ypos(yp),aquila_structure.xpos(xp),mx',subb);
   end
   %now we transpose the results to get back to the original QWY scenario
   if aquila_control.mode==2
      psi=[];
      n=length(yp);
      for i_count=0:length(E(:,1))-1
         psi=[psi ev(:,i_count*n+1:(i_count+1)*n)'];
      end
      ev=psi;
   else
      ev=ev(:)';
   end
end