function nullst=intsch(fct,anf,aend,h,gen,varargin)
%
% intsch.m ermittelt Nullstellen einer durch fct gegebenen Funktion 
% mittels Intervallschachtelung nachdem es zuerst eine Grobsuche
% mit der Schrittweite h durchgefuehrt hat.

% Achtung: 1. fct muss die x-Werte als Vektor entgegennehmen koennen und
%             die errechneten y-Werte als gleichlangen Vektor zurueckgeben.
%
% INPUT-Parameter:   fct       Name der zu untersuchenden Funktion (s.u.)
%                    anf,aend  Beginn und Ende des Grobsuchbereiches
%                    h         Schrittweite der Grobsuche
%                    gen       Fehlergrenze
%		     varargin  weitere Parameter für die Funktion muessen 
%			       in einem Array uebergeben werden
%			       Die Funktion fct muss die zusaetzlichen 
%                              Parameter auch geschlossen als Array aufnehmen 
%                              -> fct(x,param)
% 
% OUTPUT-Parameter:  nullst    Vektor der berechneten Nullstellen
%
%
% HINWEIS: Falls die Funktion fct in einem separatem File mit dem Namen
%          "aktuell" liegt, kann sie wie folgt uebergeben werden:
%	      intsch(@aktuell,anf,aend,h,gen);

%          In diesem Fall muss die Funktion im File  aktuell.m  definiert
%          sein.

%          Weitere Parameter der Funktion koennen im Vektor "varargin"
%          uebergeben werden (s.o.)., und/oder man kann die Moeglichkeit
%          "globaler Variablen" nutzen.
%                    
%
% AUTOR: Peter Lichtenberger
% DATUM: 3.Oktober 2002
%
% Erweitert: Martin Ratschek
% Datum:     10. Oktober 2006


eps=1.e-8;
err=gen+1.0;

%Hier werden die zusaetzlichen Argumente fuer die Funktion fct() abgefragt
%Diese Parameter muessen als Zahl oder als Zahlen-Array uebergeben werden 
%und wird auch dann 1-zu-1 an die Function weitergegeben
if (nargin==5)
  param=nan;
elseif (nargin == 6)
  param=varargin{1};
elseif (nargin<5)
  disp('zu wenig Parameter');
  return;
else
  disp('Zu viele Parameter fuer die Funktion "intschm"');
  disp('Hinweis: Parameter fuer die uebergebene Funktion (fct)  muessen als Array uebergeben werden');
  return;
end

x=anf:h:aend;
if (isnan(param))
  y=fct(x);
else
  try,
    y=fct(x,param);
  catch,
    disp('Falsche Anzahl an Argumenten fuer die uebergebene Funktion (fct)');
    return;
  end
end

% Begin der Grobsuche

yl=y(1:end-1);
yr=y(2:end);
comp=yl.*yr>=0;
ind=find(comp==0);

nullstgrob = x(y==0); %Falls schon in der Grobsuche eine Nullstelle getroffen wird

%Grobsuche abgeschlossen

xl=x(ind);
xr=x(ind+1);

yl=y(ind);
yr=y(ind+1);

while any(err>gen)
% eigentliche Intervallschachtelung
% wird fuer alle Nullstellen parallel durchgefuehrt
    
    x=(xl+xr)/2.0;
    if (isnan(param))
      y=fct(x);
    else
      try,
        y=fct(x,param);
      catch,
	disp('Falsche Anzahl an Argumenten');
	return;
      end
    end
    
    xl=(((yl.*y)>0.0)&(yr~=0)).*x+((yl.*y)<0.0).*xl + (yl==0).*xl + ((y==0)&(y~=yl)).*x + ((yr==0)&(yr~=yl)&(yr~=y)).*xr;
    xr=((yl.*y)<0.0).*x+((yl.*y)>0.0).*xr + (yl==0).*xl + ((y==0)&(y~=yl)).*x;
    
    yl=((yl.*y)>0.0).*y+((yl.*y)<0.0).*yl;
    yr=((yl.*y)<0.0).*y+((yl.*y)>0.0).*yr;
    
    err=(abs(x)<eps).*(xr-xl)+(abs(x)>=eps).*(xr-xl)./abs(x);
end

nullst=(xr+xl)/2.0;

nullst=sort([nullst,nullstgrob]); %ergaenzt die Nullstellen mit denen aus der Grobsuche

%Filtert doppelt vorhandene Nullstellen heraus
if(length(nullst)>1)
  hlp=nullst(2:end);
  nullst=[nullst(~(nullst(1:end-1)==hlp)),nullst(end)];
end
