function datint=mlrint(dat,IVAR,PVAR,BEREICH,Nmulti,NurStat,varargin)
%
% Eingabe:
% dat
%  Datenarray
% IVAR
%  Spalten der zu interpolierenden Parameter
% PVAR
%  Spalten der Parameter, die zur Interpolation benutzt werden sollen
% BEREICH
%  Definition der "proximity", Angabe von Variablen(Datenspalte),Wert,Prozent
%  als 4'ter Wert kann man den Multiplikationsfaktor angeben (sonst=2)
%  Dieser Wert muss>1 sein, damit der 'vergrösserte' Bereich auch
%  groesser ist. Bei Werten<0 gibt der Absolutwert an, um welchen Anteil
%  der Bereich bei jeder Vergrösserung waechst.
%   Ein Beispiel der Bereichsgrösse mit 1 als Anfangsgrösse
%  Bei +2  : 1 ; 2  ; 4 ; 8
%  Bei +1  : 1 ; 1  ; 1 ; 1
%  Bei -0.5: 1 ; 1.5; 2 ; 2.5 
%  Bei -1  : 1 ; 2  ; 3 ; 4 
% Nmulti
%  maximale Anzahl der Bereichsmultiplikationen
% NurStat
%   1=Es werden nur vorhandene Werte für Statistik berechnet
%  -1=Es werden nur nicht vorhandene Werte berechnet
%   0=Es werden alle Werte berechnet, auch für Statistik
MAXDATA=0;
ii=0;
while ii+7< nargin
  ii=ii+1;
  if strcmp(varargin{ii},'MAXDATA')
    ii=ii+1;
    MAXDATA=varargin{ii};
  else
    error('unknown variable argument')
  end
end

  
[ii,jj]=size(BEREICH);
if ii>=1;
 if jj<=2; BEREICH=[BEREICH zeros(ii,1)]; end
 if jj==3; BEREICH=[BEREICH ones(ii,1)*2]; end
 jj=BEREICH(:,4)<=0; BEREICH(jj,4)=1;
 if max(BEREICH(:,3))>1;   % alle Werte<1 <> schon in anteile 
   BEREICH(:,3)=BEREICH(:,3)/100.;
 end
end
% Angabe der maximalen Anzahl von Multiplikationen
% bei =1 wird also gegebenenfalls einmal multipliziert
if nargin<5; Nmulti=1; end
if Nmulti<0; Nmulti=0; end

if nargin<6; NurStat=-1; end

[IDAT,JDAT]=size(dat);
if MAXDATA==0; MAXDATA=IDAT; end
% Ausgabearray erstmal alles auf NaN
datint=zeros(MAXDATA,JDAT)+NaN;


%% alle Bereichsvariablen duerfen keine NaNs haben
%for ii=BEREICH(:,1)
%  ij=isnan(dat(:,ii));
%  if sum(ij)>0;
%      dat(ij,:)=NaN;
%%     fprintf('Bereichsvariable, %i\n',ii);
%%     error('Bereichsvariable mit NaNs')
%  end
%end

for ivar=IVAR
 for idat=1:MAXDATA
   ij=0;   % Zyklen mit einer Bereichsvariable=NaN werden nicht
   for ii=BEREICH(:,1)   % interpoliert
       if isnan(dat(idat,ii)); ij=1; break; end
   end
   if ij==1; continue; end
   
   pvar=[];  % vorhandene Interpolationsparameter
   for ii=PVAR
       if ~isnan(dat(idat,ii)) & ii~=ivar ;
          pvar=[pvar ii];
       end    
   end
   jj=zeros(IDAT,1); Mmulti=0;
   
   if isnan(dat(idat,ivar)) & NurStat==1; Mmulti=Nmulti+1; end
   if ~isnan(dat(idat,ivar)) & NurStat==-1; Mmulti=Nmulti+1; end

   while Mmulti<=Nmulti
    jj = ~isnan(dat(:,ivar));       % zu interpolierender Parameter muss da sein
    jj(idat) = 0;                   % aktueller Wert wird nicht benutzt
    for ii=pvar
      jj= jj &  ~isnan(dat(:,ii));  % Interpolationsparameter müssen da sein
    end
    for ii=1:length(BEREICH(:,1))             % Bereichsabfragen
      bereich=BEREICH(ii,2)+dat(idat,BEREICH(ii,1))*BEREICH(ii,3);
      if BEREICH(ii,4)<0;
         bereich=bereich+bereich*abs(BEREICH(ii,4))*Mmulti;
      else
         bereich=bereich*BEREICH(ii,4)^Mmulti;
      end
      jj= jj & abs(dat(:,BEREICH(ii,1))-dat(idat,BEREICH(ii,1))) <= bereich;
    end
    Mmulti=Mmulti+1;
    if sum(jj)>length(pvar)*1.5;    % Ueberpruefen ob genuegend Zyklen vorhanden
       a=dat(jj,pvar); b=dat(jj,ivar);
       faktoren=pinv(a)*b;
       datint(idat,ivar)=dat(idat,pvar)*faktoren;
       break; 
    end
   end
 end

 % Ausgabe einiger Informationen
 fprintf('Variable: %i\n',ivar);
 fprintf('Anzahl gemessener     Zyklen: %i\n',sum(~isnan(dat(:,ivar))));
 fprintf('Anzahl interpolierter Zyklen: %i\n',sum(~isnan(datint(:,ivar))));
 jj= ~isnan(dat(1:MAXDATA,ivar)) & ~isnan(datint(:,ivar));
 fprintf('Anzahl der Vergleichs-Zyklen: %i\n',sum(jj));
 if sum(jj)>2
    Diff=dat(jj,ivar)-datint(jj,ivar);
    fprintf('min=%12f ,max=%12f ,std=%12f\n',min(Diff),max(Diff),std(Diff));
 end
end
