unit werteunit;

{$mode objfpc}{$H+}

interface

uses
  Classes, SysUtils, typenunit, math, process, lowlevelunit, matheunit, fftunit;

type
// tLLWerte ********************************************************************
  pTLLWerteSingle = ^tLLWerteSingle;
  pTLLWerteDouble = ^tLLWerteDouble;
  pTLLWerteExtended = ^tLLWerteExtended;
  generic tLLWerte<wgen> = class(tObject)
{
    Diese Klasse stellt nur die Berechnungsroutinen und ähnliches bereit,
    nicht jedoch ein (wie auch immer geartetes) Userinterface.
    Auch die korrekte Parallelisierung obliegt dem übergeordneten Programmteil.
}
  private
  public
    werte:                                       array of wgen;
    params:                                      tExtrainfos;
    constructor create(ps: tExtrainfos); overload;
    constructor create(original: pTLLWerteSingle; ps: tExtrainfos; xmin,xmax: longint); overload;
    constructor create(original: pTLLWerteDouble; ps: tExtrainfos; xmin,xmax: longint); overload;
    constructor create(original: pTLLWerteExtended; ps: tExtrainfos; xmin,xmax: longint); overload;
    procedure kopiereVon(st: boolean; original: pTLLWerteSingle); overload;
    procedure kopiereVon(st: boolean; original: pTLLWerteDouble); overload;
    procedure kopiereVon(st: boolean; original: pTLLWerteExtended); overload;
    procedure kopiereVon(st: boolean; original: pTLLWerteSingle; xmin,xmax: longint); overload;
    procedure kopiereVon(st: boolean; original: pTLLWerteDouble; xmin,xmax: longint); overload;
    procedure kopiereVon(st: boolean; original: pTLLWerteExtended; xmin,xmax: longint); overload;
    procedure kopiereVon(st: boolean; original: pTLLWerteSingle; xmin,xmax,tmin,tmax: longint); overload;
    procedure kopiereVon(st: boolean; original: pTLLWerteDouble; xmin,xmax,tmin,tmax: longint); overload;
    procedure kopiereVon(st: boolean; original: pTLLWerteExtended; xmin,xmax,tmin,tmax: longint); overload;
    procedure kopiereVonNach(original: pTLLWerteSingle; qxmin,qxmax,qtmin,qtmax,zxmin,ztmin: longint); overload;
    procedure kopiereVonNach(original: pTLLWerteDouble; qxmin,qxmax,qtmin,qtmax,zxmin,ztmin: longint); overload;
    procedure kopiereVonNach(original: pTLLWerteExtended; qxmin,qxmax,qtmin,qtmax,zxmin,ztmin: longint); overload;
    procedure kopiereVerzerrt(original: pTLLWerteSingle; ZPs: tIntPointArray; ZGs: tExtPointArray; ZAs: tExtendedArray; xmin,xmax,tmin,tmax: longint; vb,nb: tTransformation; va,na: longint); overload;
    procedure kopiereVerzerrt(original: pTLLWerteDouble; ZPs: tIntPointArray; ZGs: tExtPointArray; ZAs: tExtendedArray; xmin,xmax,tmin,tmax: longint; vb,nb: tTransformation; va,na: longint); overload;
    procedure kopiereVerzerrt(original: pTLLWerteExtended; ZPs: tIntPointArray; ZGs: tExtPointArray; ZAs: tExtendedArray; xmin,xmax,tmin,tmax: longint; vb,nb: tTransformation; va,na: longint); overload;
    destructor destroy; override;
    function liesDateien(dateien: tGenerischeInputDateiInfoArray): boolean;
    function fft(senkrecht,invers: boolean; algo: tFFTAlgorithmus; fen: tFenster; out pvFehler: extended): boolean; overload; inline;
    function fft(smin,smax: longint; senkrecht,invers: boolean; algo: tFFTAlgorithmus; fen: tFenster; out pvFehler: extended): boolean; overload;
    procedure spiegle; overload;
    procedure spiegle(tmin,tmax: longint); overload;
    procedure fft2dNachbearbeitungA(nb: tFFTDatenordnung);
    procedure fft2dNachbearbeitungB(xmin,xmax: longint; nb: tFFTDatenordnung);
    procedure schreibeWert(var f: textfile; x,y: longint); overload;
    procedure schreibeWert(var f: textfile; x,y,wert: extended); overload;
    procedure schreibeWertIntegriert(var f: textfile; i: longint; hor: boolean);
    procedure erzeugeBinning(senkrecht,linien: boolean; x0,dx: extended; s: string);
    procedure holeRam; inline; overload;
    procedure holeRam(ausgaben: byte); inline; overload;
    procedure holeRam(ausgaben: byte; gemaeszTXMinMax: boolean); inline; overload;
    procedure gibMinMaxDichten(out wMi,wMa: extended; xmin,xmax,tmin,tmax: longint);
    function zuPixelWerten(whoehe,wbreite,xpmi,xmi,tmi: longint; xz,yz: extended; pPWerte: pTExtendedArray; pPAnzahlen: pTLongintArray): boolean;
    function findeSchwellwerte(xmi,xma,tmi,tma: longint; Schw: extended): tExtPointArray;
    procedure integriereSingle(qu: pTLLWerteSingle; xmi,xma,tmi,tma,xof,tof: longint; richtung: tIntegrationsRichtung);
    procedure integriereDouble(qu: pTLLWerteDouble; xmi,xma,tmi,tma,xof,tof: longint; richtung: tIntegrationsRichtung);
    procedure integriereExtended(qu: pTLLWerteDouble; xmi,xma,tmi,tma,xof,tof: longint; richtung: tIntegrationsRichtung);
    procedure gauszFit(amplituden,breiten,positionen,ueberlappe,hintergruende: pTLLWerteExtended; von,bis: longint; senkrecht: boolean; fensterBreite,maxBreite,maxVerschiebung: extended; positionsMitten: tExtendedArray);
  end;
  tLLWerteSingle = specialize tLLWerte<single>;
  tLLWerteDouble = specialize tLLWerte<double>;
  tLLWerteExtended = specialize tLLWerte<extended>;
// andere Typen ****************************************************************
  tWavelet = class(tObject)
    freq,tfwhm,pvFehler: extended;
    hlen:                longint;
    typ:                 tWaveletTyp;
    werte:               tLLWerteDouble;
    mitFFT:              boolean;
    constructor create;
    destructor destroy; override;
    function setzeTyp(s: string): boolean;
    function berechneWerte: boolean;
  end;

const
  anzSergeyFelder = 12;

implementation

uses systemunit;

constructor tLLWerte.create(ps: tExtrainfos);
begin
  inherited create;
  params:=ps;
  setlength(werte,0);
end;

constructor tLLWerte.create(original: pTLLWerteSingle; ps: tExtrainfos; xmin,xmax: longint);
begin
  inherited create;
  params:=ps;
  kopiereVon(false,original,xmin,xmax);
end;

constructor tLLWerte.create(original: pTLLWerteDouble; ps: tExtrainfos; xmin,xmax: longint);
begin
  inherited create;
  params:=ps;
  kopiereVon(false,original,xmin,xmax);
end;

constructor tLLWerte.create(original: pTLLWerteExtended; ps: tExtrainfos; xmin,xmax: longint);
begin
  inherited create;
  params:=ps;
  kopiereVon(false,original,xmin,xmax);
end;

procedure tLLWerte.kopiereVon(st: boolean; original: pTLLWerteSingle);
begin
  kopiereVon(st,original,0,original^.params.xsteps-1);
end;

procedure tLLWerte.kopiereVon(st: boolean; original: pTLLWerteDouble);
begin
  kopiereVon(st,original,0,original^.params.xsteps-1);
end;

procedure tLLWerte.kopiereVon(st: boolean; original: pTLLWerteExtended);
begin
  kopiereVon(st,original,0,original^.params.xsteps-1);
end;

procedure tLLWerte.kopiereVon(st: boolean; original: pTLLWerteSingle; xmin,xmax: longint);
begin
  kopiereVon(st,original,xmin,xmax,0,original^.params.tsiz-1);
end;

procedure tLLWerte.kopiereVon(st: boolean; original: pTLLWerteDouble; xmin,xmax: longint);
begin
  kopiereVon(st,original,xmin,xmax,0,original^.params.tsiz-1);
end;

procedure tLLWerte.kopiereVon(st: boolean; original: pTLLWerteExtended; xmin,xmax: longint);
begin
  kopiereVon(st,original,xmin,xmax,0,original^.params.tsiz-1);
end;

procedure tLLWerte.kopiereVon(st: boolean; original: pTLLWerteSingle; xmin,xmax,tmin,tmax: longint);
var
  i,j: longint;
begin
  inherited create;
  tmax:=min(tmax,original^.params.tsiz-1);
  tmin:=max(tmin,0);
  params.tsiz:=tmax+1-tmin;
  xmax:=min(xmax,original^.params.xsteps-1);
  xmin:=max(xmin,0);
  params.xsteps:=xmax+1-xmin;
  zerstoereTransformationWennObsolet(params.transformationen);
  params.transformationen:=tKoordinatenAusschnitt.create(original^.params.transformationen,xmin,xmax,tmin,tmax);
  params.maxW:=0;
  params.minW:=0;
  params.np:=original^.params.np;
  params.beta:=original^.params.beta;
  params.refreshKnownValues;
  if not st then begin
    holeRam(0);
    for i:=xmin to xmax do
      for j:=tmin to tmax do
        werte[i-xmin+(j-tmin)*params.xsteps]:=original^.werte[i+j*original^.params.xsteps];
  end;
end;

procedure tLLWerte.kopiereVon(st: boolean; original: pTLLWerteDouble; xmin,xmax,tmin,tmax: longint);
var
  i,j: longint;
begin
  inherited create;
  tmax:=min(tmax,original^.params.tsiz-1);
  tmin:=max(tmin,0);
  params.tsiz:=tmax+1-tmin;
  xmax:=min(xmax,original^.params.xsteps-1);
  xmin:=max(xmin,0);
  params.xsteps:=xmax+1-xmin;
  zerstoereTransformationWennObsolet(params.transformationen);
  params.transformationen:=tKoordinatenAusschnitt.create(original^.params.transformationen,xmin,xmax,tmin,tmax);
  params.maxW:=0;
  params.minW:=0;
  params.np:=original^.params.np;
  params.beta:=original^.params.beta;
  params.refreshKnownValues;
  if not st then begin
    holeRam(0);
    for i:=xmin to xmax do
      for j:=tmin to tmax do
        werte[i-xmin+(j-tmin)*params.xsteps]:=original^.werte[i+j*original^.params.xsteps];
  end;
end;

procedure tLLWerte.kopiereVon(st: boolean; original: pTLLWerteExtended; xmin,xmax,tmin,tmax: longint);
var
  i,j: longint;
begin
  inherited create;
  tmax:=min(tmax,original^.params.tsiz-1);
  tmin:=max(tmin,0);
  params.tsiz:=tmax+1-tmin;
  xmax:=min(xmax,original^.params.xsteps-1);
  xmin:=max(xmin,0);
  params.xsteps:=xmax+1-xmin;
  zerstoereTransformationWennObsolet(params.transformationen);
  params.transformationen:=tKoordinatenAusschnitt.create(original^.params.transformationen,xmin,xmax,tmin,tmax);
  params.maxW:=0;
  params.minW:=0;
  params.np:=original^.params.np;
  params.beta:=original^.params.beta;
  params.refreshKnownValues;
  if not st then begin
    holeRam(0);
    for i:=xmin to xmax do
      for j:=tmin to tmax do
        werte[i-xmin+(j-tmin)*params.xsteps]:=original^.werte[i+j*original^.params.xsteps];
  end;
end;

procedure tLLWerte.kopiereVonNach(original: pTLLWerteSingle; qxmin,qxmax,qtmin,qtmax,zxmin,ztmin: longint);
var
  i,j: longint;
begin
  inherited create;
  for i:=qxmin to qxmax do
    for j:=qtmin to qtmax do
      werte[i-qxmin+zxmin + (j-qtmin+ztmin)*params.xsteps]:=
        original^.werte[i+j*original^.params.xsteps];
end;

procedure tLLWerte.kopiereVonNach(original: pTLLWerteDouble; qxmin,qxmax,qtmin,qtmax,zxmin,ztmin: longint);
var
  i,j: longint;
begin
  inherited create;
  for i:=qxmin to qxmax do
    for j:=qtmin to qtmax do
      werte[i-qxmin+zxmin + (j-qtmin+ztmin)*params.xsteps]:=
        original^.werte[i+j*original^.params.xsteps];
end;

procedure tLLWerte.kopiereVonNach(original: pTLLWerteExtended; qxmin,qxmax,qtmin,qtmax,zxmin,ztmin: longint);
var
  i,j: longint;
begin
  inherited create;
  for i:=qxmin to qxmax do
    for j:=qtmin to qtmax do
      werte[i-qxmin+zxmin + (j-qtmin+ztmin)*params.xsteps]:=
        original^.werte[i+j*original^.params.xsteps];
end;

procedure tLLWerte.kopiereVerzerrt(original: pTLLWerteSingle; ZPs: tIntPointArray; ZGs: tExtPointArray; ZAs: tExtendedArray; xmin,xmax,tmin,tmax: longint; vb,nb: tTransformation; va,na: longint);
var
  i,j,k: longint;
  tmp:   extended;
begin
  for i:=tmin to tmax do
    for j:=xmin to xmax do
      werte[j+i*params.xsteps]:=0;
  for i:=0 to length(ZPs)-1 do
    for j:=0 to 1 do
      for k:=0 to 1 do
        if (ZPs[i]['x']+j>=xmin) and (ZPs[i]['x']+j<=xmax) and
            (ZPs[i]['y']+k>=tmin) and (ZPs[i]['y']+k<=tmax) then begin
          tmp:=original^.werte[i];
          if (va>0) or (na>0) then
            tmp:=(tmp-original^.params.minW)/(original^.params.maxW-original^.params.minW);
          if va>0 then
            vb.transformiereWert(tmp,va-1);
          tmp:=tmp * (ZGs[i]['x'] * (2*j-1) + 1-j) * (ZGs[i]['y'] * (2*k-1) + 1-k);
          werte[ZPs[i]['x']+j + (ZPs[i]['y']+k)*params.xsteps]:=
            werte[ZPs[i]['x']+j + (ZPs[i]['y']+k)*params.xsteps] +
            tmp;
        end;
  for i:=tmin to tmax do
    for j:=xmin to xmax do begin
      tmp:=werte[j + i*params.xsteps] / ZAs[j + i*params.xsteps];
      if na>0 then
        tmp:=nb.transformiereWert(tmp,na-1);
      werte[j + i*params.xsteps]:=tmp;
    end;
end;

procedure tLLWerte.kopiereVerzerrt(original: pTLLWerteDouble; ZPs: tIntPointArray; ZGs: tExtPointArray; ZAs: tExtendedArray; xmin,xmax,tmin,tmax: longint; vb,nb: tTransformation; va,na: longint);
var
  i,j,k: longint;
  tmp:   extended;
begin
  for i:=tmin to tmax do
    for j:=xmin to xmax do
      werte[j+i*params.xsteps]:=0;
  for i:=0 to length(ZPs)-1 do
    for j:=0 to 1 do
      for k:=0 to 1 do
        if (ZPs[i]['x']+j>=xmin) and (ZPs[i]['x']+j<=xmax) and
          (ZPs[i]['y']+k>=tmin) and (ZPs[i]['y']+k<=tmax) then begin
          tmp:=original^.werte[i];
          if (va>0) or (na>0) then
            tmp:=(tmp-original^.params.minW)/(original^.params.maxW-original^.params.minW);
          if va>0 then
            vb.transformiereWert(tmp,va-1);
          tmp:=tmp * (ZGs[i]['x'] * (2*j-1) + 1-j) * (ZGs[i]['y'] * (2*k-1) + 1-k);
          werte[ZPs[i]['x']+j + (ZPs[i]['y']+k)*params.xsteps]:=
            werte[ZPs[i]['x']+j + (ZPs[i]['y']+k)*params.xsteps] +
            tmp;
        end;
  for i:=tmin to tmax do
    for j:=xmin to xmax do begin
      tmp:=werte[j + i*params.xsteps] / ZAs[j + i*params.xsteps];
      if na>0 then
        tmp:=nb.transformiereWert(tmp,na-1);
      werte[j + i*params.xsteps]:=tmp;
    end;
end;

procedure tLLWerte.kopiereVerzerrt(original: pTLLWerteExtended; ZPs: tIntPointArray; ZGs: tExtPointArray; ZAs: tExtendedArray; xmin,xmax,tmin,tmax: longint; vb,nb: tTransformation; va,na: longint);
var
  i,j,k: longint;
  tmp:   extended;
begin
  for i:=tmin to tmax do
    for j:=xmin to xmax do
      werte[j+i*params.xsteps]:=0;
  for i:=0 to length(ZPs)-1 do
    for j:=0 to 1 do
      for k:=0 to 1 do
        if (ZPs[i]['x']+j>=xmin) and (ZPs[i]['x']+j<=xmax) and
          (ZPs[i]['y']+k>=tmin) and (ZPs[i]['y']+k<=tmax) then begin
          tmp:=original^.werte[i];
          if (va>0) or (na>0) then
            tmp:=(tmp-original^.params.minW)/(original^.params.maxW-original^.params.minW);
          if va>0 then
            vb.transformiereWert(tmp,va-1);
          tmp:=tmp * (ZGs[i]['x'] * (2*j-1) + 1-j) * (ZGs[i]['y'] * (2*k-1) + 1-k);
          werte[ZPs[i]['x']+j + (ZPs[i]['y']+k)*params.xsteps]:=
            werte[ZPs[i]['x']+j + (ZPs[i]['y']+k)*params.xsteps] +
            tmp;
        end;
  for i:=tmin to tmax do
    for j:=xmin to xmax do begin
      tmp:=werte[j + i*params.xsteps] / ZAs[j + i*params.xsteps];
      if na>0 then
        tmp:=nb.transformiereWert(tmp,na-1);
      werte[j + i*params.xsteps]:=tmp;
    end;
end;

destructor tLLWerte.destroy;
begin
  setlength(werte,0);
  inherited destroy;
end;

function tLLWerte.liesDateien(dateien: tGenerischeInputDateiInfoArray): boolean;
var
  i,j,k,l,tmpi,etsiz,spAnz,br: longint;
  f:                           file;
  tmps:                        single;
  tmpd:                        double;
  tmpe,Zeit:                   extended;
  sa:                          tSingleArray;
  da:                          tDoubleArray;
  ea:                          tExtendedArray;
  ipp:                         tProcess;
  buf:                         tByteArray;
  etwasGelesen:                boolean;
begin
  result:=false;
  gibAus('... Dateien einlesen ...',1);
  zeit:=now;
  tmpi:=0;
  tmps:=0;
  tmpd:=0;
  etsiz:=0;
  spAnz:=-1;
  for i:=0 to length(dateien)-1 do begin
    gibAus(' '+dateien[i].Name,1);
    etwasGelesen:=false;
    if dateien[i] is tPipeInputDateiInfo then begin
      if ((dateien[i] as tPipeInputDateiInfo).bytesPerSample<>4) or
        ((dateien[i] as tPipeInputDateiInfo).Kodierung<>k32BitSignedInteger) then begin
        gibAus('Ich kann nur vier Bytes mit einem mal als Integer interpretiert aus einer Pipe einlesen!',3);
        exit;
      end;
      tmpe:=power(2,-31);
      ipp:=tProcess.create(nil);
      ipp.Options:=ipp.Options + [poUsePipes];
      ipp.Executable:=(dateien[i] as tPipeInputDateiInfo).Executable;
      ipp.Parameters.Text:=(dateien[i] as tPipeInputDateiInfo).ParametersText;
      ipp.execute;
      setlength(buf,0);
      br:=0;
      while ipp.running or (ipp.Output.NumBytesAvailable>0) do begin
        if ipp.Output.NumBytesAvailable > 0 then begin
          if br+ipp.Output.NumBytesAvailable>=length(buf) then
            setlength(buf,br+ipp.Output.NumBytesAvailable+65536*1024);
          tmpi:=ipp.Output.Read(buf[br],min(ipp.Output.NumBytesAvailable,length(buf)-br));
          if ((br+tmpi) shr 24) > (br shr 24) then gibAus(inttostr(br+tmpi)+' von '+inttostr(dateien[i].tsiz*dateien[i].xsteps*4)+' Bytes bisher gelesen ('+floattostrtrunc((br+tmpi)/(dateien[i].tsiz*dateien[i].xsteps*4)*100,2,true)+'%).',1);
          br:=br+tmpi;
        end
        else
          sleep(10);
      end;
      ipp.free;
      gibAus('insgesamt '+inttostr(br div 1024 div 1024)+' MB gelesen.',1);
      setlength(buf,br);
      if dateien[i].tsiz*dateien[i].xsteps*4>length(buf) then begin
        gibAus('Ich habe '+inttostr(length(buf))+' Bytes aus der Pipe gelesen, anstelle von wenigstens '+inttostr(dateien[i].tsiz*dateien[i].xsteps*4)+', wie erwartet!',3);
        setlength(buf,0);
        exit;
      end;
      tmpi:=length(buf)-4*dateien[i].tsiz*dateien[i].xsteps;  // der Offset des ersten Daten-Wortes
      for j:=dateien[i].tmin to dateien[i].tmax do
        for k:=dateien[i].xmin to dateien[i].xmax do
          werte[dateien[i].t0abs+ k-dateien[i].xmin + (j-dateien[i].tmin)*(dateien[i].xmax-dateien[i].xmin+1)]:=
            int32((((((buf[tmpi+3+4*(k+j*dateien[i].xsteps)] shl 8) or
                      buf[tmpi+2+4*(k+j*dateien[i].xsteps)]) shl 8) or
                    buf[tmpi+1+4*(k+j*dateien[i].xsteps)]) shl 8) or
                    buf[tmpi+4*(k+j*dateien[i].xsteps)]) * tmpe;
      setlength(buf,0);
      if etwasGelesen then begin
        gibAus('Ich habe diese Runde schon Daten gelesen!',3);
        exit;
      end;
      etwasGelesen:=true;
    end;
    if (dateien[i] is tSpaceTimeInputDateiInfo) or
      (dateien[i] is tTraceInputDateiInfo) then begin
      assign(f,dateien[i].Name);
      reset(f,1);
      blockread(f,tmpi,sizeof(integer));
      dec(tmpi);
      if tmpi-round(params.tstart/dateien[i].groeszenFaktor)<>i then begin
        gibAus('Datei '''+dateien[i].Name+''' kommt nicht an '+inttostr(i)+'-ter Stelle, wie sie sollte, sondern an '+inttostr(tmpi-round(params.tstart/dateien[i].groeszenFaktor))+'-ter.',3);
        close(f);
        exit;
      end;
      if dateien[i] is tTraceInputDateiInfo then begin
        blockread(f,tmpi,sizeof(integer));  // #Traces
        spAnz:=tmpi;
      end;
      blockread(f,etsiz,sizeof(integer));
      if dateien[i] is tSpaceTimeInputDateiInfo then begin
        for j:=0 to etsiz-1 do begin
          case Dateien[i].Genauigkeit of
            gSingle: begin
              blockread(f,tmps,sizeof(single)); // xstart
              tmpe:=tmps;
            end;
            gDouble: begin
              blockread(f,tmpd,sizeof(double)); // xstart
              tmpe:=tmpd;
            end;
            gExtended:
              blockread(f,tmpe,sizeof(extended)); // xstart
          end{of Case};
          tmpe:=tmpe*dateien[i].groeszenFaktor;
          if j=0 then params.transformationen.xstart:=tmpe;
          if tmpe<>params.transformationen.xstart then begin
            gibAus('Falscher linker Rand in '''+dateien[i].Name+''' im Schritt '+inttostr(j)+', nämlich '+myfloattostr(tmpe)+' statt '+myfloattostr(params.transformationen.xstart)+'!',3);
            close(f);
            exit;
          end;
          case Dateien[i].Genauigkeit of
            gSingle: begin
              blockread(f,tmps,sizeof(single)); // xstop
              tmpe:=tmps;
            end;
            gDouble: begin
              blockread(f,tmpd,sizeof(double)); // xstop
              tmpe:=tmpd;
            end;
            gExtended:
              blockread(f,tmpe,sizeof(extended)); // xstop
          end{of Case};
          tmpe:=tmpe*dateien[i].groeszenFaktor;
          if j=0 then params.transformationen.xstop:=tmpe;
          if tmpe<>params.transformationen.xstop then begin
            gibAus('Falscher rechter Rand in '''+dateien[i].Name+''' im Schritt '+inttostr(j)+', nämlich '+myfloattostr(tmpe)+' statt '+myfloattostr(params.transformationen.xstop)+'!',3);
            close(f);
            exit;
          end;
          blockread(f,tmpi,sizeof(integer)); // xsteps
          if tmpi<>params.xsteps then begin
            gibAus('Falsche Anzahl an x-Schritten in '''+dateien[i].Name+''' im Schritt '+inttostr(j)+', nämlich '+inttostr(tmpi)+' statt '+myfloattostr(params.xsteps)+'!',3);
            close(f);
            exit;
          end;
          if sizeof(wgen) = wertGroesze(Dateien[i].Genauigkeit) then
            blockread(f,werte[(j+Dateien[i].t0abs)*params.xsteps],params.xsteps*sizeof(wgen))
          else begin
            setlength(sa,params.xsteps);
            blockread(f,sa[0],params.xsteps*sizeof(single));
            for k:=0 to params.xsteps-1 do
              werte[(j+Dateien[i].t0abs)*params.xsteps+k]:=sa[k];
          end;
          if power(dateien[i].gamma,3)/sqr(dateien[i].groeszenFaktor)<>1 then    // gamma^3 als Skalierungsfaktor für Dichten ?
            for k:=0 to params.xsteps-1 do
              werte[(j+Dateien[i].t0abs)*params.xsteps+k]:=werte[(j+Dateien[i].t0abs)*params.xsteps+k]*power(dateien[i].gamma,3)/sqr(dateien[i].groeszenFaktor);
        end;
        if etwasGelesen then begin
          gibAus('Ich habe diese Runde schon Daten gelesen!',3);
          exit;
        end;
        etwasGelesen:=true;
      end;
      if dateien[i] is tTraceInputDateiInfo then begin
        case Dateien[i].Genauigkeit of
          gSingle: begin
            setlength(sa,etsiz);
            setlength(da,0);
            setlength(ea,0);
          end;
          gDouble: begin
            setlength(sa,0);
            setlength(da,etsiz);
            setlength(ea,0);
          end;
          gExtended: begin
            setlength(sa,0);
            setlength(da,0);
            setlength(ea,etsiz);
          end;
        end{of case};
        for j:=0 to spAnz-1 do
          case Dateien[i].Genauigkeit of
            gSingle: begin
              blockread(f,tmps,sizeof(single));   // x
              if j=(Dateien[i] as tTraceInputDateiInfo).Spurnummer then begin
                params.transformationen.xstop:=tmps;
                params.transformationen.xstart:=params.xstop;
              end;
              for k:=0 to length(FeldgroeszenNamen)-1 do begin
                blockread(f,sa[0],sizeof(single)*length(sa));
                if (j=(Dateien[i] as tTraceInputDateiInfo).Spurnummer) and
                  (k=(Dateien[i] as tTraceInputDateiInfo).Feldnummer) then begin
                  for l:=0 to length(sa)-1 do
                    werte[l+Dateien[i].t0abs]:=sa[l]*sqr(dateien[i].gamma)/sqr(dateien[i].groeszenFaktor);
                  if etwasGelesen then begin
                    gibAus('Ich habe diese Runde schon Daten gelesen!',3);
                    exit;
                  end;
                  etwasGelesen:=true;
                end;
              end;
            end;
            gDouble: begin
              blockread(f,tmpd,sizeof(double)); // x
              if j=(Dateien[i] as tTraceInputDateiInfo).Spurnummer then begin
                params.transformationen.xstop:=tmpd;
                params.transformationen.xstart:=params.xstop;
              end;
              for k:=0 to length(FeldgroeszenNamen)-1 do begin
                blockread(f,da[0],sizeof(double)*length(da));
                if (j=(Dateien[i] as tTraceInputDateiInfo).Spurnummer) and
                  (k=(Dateien[i] as tTraceInputDateiInfo).Feldnummer) then begin
                  for l:=0 to length(da)-1 do
                    werte[l+dateien[i].t0abs]:=da[l]*sqr(dateien[i].gamma)/sqr(dateien[i].groeszenFaktor);
                  if etwasGelesen then begin
                    gibAus('Ich habe diese Runde schon Daten gelesen!',3);
                    exit;
                  end;
                  etwasGelesen:=true;
                end;
              end;
            end;
            gExtended: begin
              blockread(f,tmpe,sizeof(extended)); // x
              if j=(Dateien[i] as tTraceInputDateiInfo).Spurnummer then begin
                params.transformationen.xstop:=tmpe;
                params.transformationen.xstart:=params.xstop;
              end;
              for k:=0 to length(FeldgroeszenNamen)-1 do begin
                blockread(f,ea[0],sizeof(extended)*length(ea));
                if (j=(Dateien[i] as tTraceInputDateiInfo).Spurnummer) and
                  (k=(Dateien[i] as tTraceInputDateiInfo).Feldnummer) then begin
                  for l:=0 to length(ea)-1 do
                    werte[l+dateien[i].t0abs]:=ea[l]*sqr(dateien[i].gamma)/sqr(dateien[i].groeszenFaktor);
                  if etwasGelesen then begin
                    gibAus('Ich habe diese Runde schon Daten gelesen!',3);
                    exit;
                  end;
                  etwasGelesen:=true;
                end;
              end;
            end;
          end{of Case};
      end;
      if not eof(f) then begin
        gibAus('Zu viele Daten in '''+dateien[i].Name+'''!',3);
        close(f);
        exit;
      end;
      close(f);
    end;
    if dateien[i] is tPhaseSpaceInputDateiInfo then begin
      if i<>0 then begin
        gibAus('Ich kann Phasenraumdateien nicht kaskadieren!',3);
        close(f);
        exit;
      end;
      assign(f,dateien[i].Name);
      reset(f,1);
      seek(f,filepos(f)
        + 4*wertGroesze(dateien[i].genauigkeit) // xstart,xstop,tstart,tstop
        + 2*sizeof(longint));                   // xsteps,tsiz
      blockread(f,werte[0],params.xsteps*params.tsiz*wertGroesze(dateien[i].genauigkeit));
      close(f);
      etwasGelesen:=true;
    end;
    if dateien[i] is tSergeyInputDateiInfo then begin
      etsiz:=Dateien[i].tsiz;

      setlength(buf,sizeof(extended)*anzSergeyFelder);
      tmpi:=wertGroesze(Dateien[i].Genauigkeit);

      assign(f,dateien[i].Name+'traces/traces.dat');
      reset(f,1);

      setlength(buf,tmpi*anzSergeyFelder);

      for j:=0 to etsiz-1 do begin
        if (Dateien[i] as tSergeyInputDateiInfo).Feldnummer>0 then
          blockread(f,buf[0],(Dateien[i] as tSergeyInputDateiInfo).Feldnummer*tmpi);
        case dateien[i].genauigkeit of
          gSingle: begin
            blockread(f,tmps,tmpi);
            werte[j]:=tmps;
          end;
          gDouble: begin
            blockread(f,tmpd,tmpi);
            werte[j]:=tmpd;
          end;
          gExtended: begin
            blockread(f,tmpe,tmpi);
            werte[j]:=tmpe;
          end;
        end{of Case};
        if (Dateien[i] as tSergeyInputDateiInfo).Feldnummer<anzSergeyFelder-1 then
          blockread(f,buf[0],(anzSergeyFelder-1-(Dateien[i] as tSergeyInputDateiInfo).Feldnummer)*tmpi);
      end;

      setlength(buf,0);
      if not eof(f) then begin
        gibAus('Zu viele Daten in '''+dateien[i].Name+'traces/traces.dat''!',3);
        close(f);
        exit;
      end;
      close(f);
      etwasGelesen:=true;
    end;
    if not etwasGelesen then begin
      gibAus('Ich habe diese Runde keine Daten gelesen!',3);
      exit;
    end;
  end;
  params.refreshKnownValues;
  gibAus('... fertig '+timetostr(now-Zeit),1);
  result:=true;
end;

procedure tLLWerte.gibMinMaxDichten(out wMi,wMa: extended; xmin,xmax,tmin,tmax: longint);
var
  i,j:    longint;
begin
  wMi:=werte[xmin+tmin*params.xsteps];
  wMa:=Werte[xmin+tmin*params.xsteps];
  for i:=xmin to xmax do
    for j:=tmin to tmax do begin
      wMi:=min(wMi,werte[i+j*params.xsteps]);
      wMa:=max(wMa,werte[i+j*params.xsteps]);
    end;
end;

function tLLWerte.fft(senkrecht,invers: boolean; algo: tFFTAlgorithmus; fen: tFenster; out pvFehler: extended): boolean;
begin
  if senkrecht then
    result:=fft(0,params.xsteps-1,senkrecht,invers,algo,fen,pvFehler)
  else
    result:=fft(0,params.tsiz-1,senkrecht,invers,algo,fen,pvFehler);
end;

function tLLWerte.fft(smin,smax: longint; senkrecht,invers: boolean; algo: tFFTAlgorithmus; fen: tFenster; out pvFehler: extended): boolean;
var
  i,j,pmax,pstep,sstep:  longint;
  in0,out0:              boolean;
  vorher,nachher:        extended;
begin
  result:=false;
  if senkrecht then begin
    pmax:=params.tsiz-1;
    pstep:=params.xsteps;
    sstep:=1;
  end
  else begin
    pmax:=params.xsteps-1;
    pstep:=1;
    sstep:=params.xsteps;
  end;

  if assigned(fen) and fen.aktiv then begin
    if length(fen.Werte)<>pmax+1 then
      fen.berechneWerte(pmax+1);
    for i:=smin to smax do  // Werte fenstern
      for j:=0 to pmax do
        werte[i*sstep+j*pstep]:=
          werte[i*sstep+j*pstep]*fen.Werte[j];
  end;

  vorher:=0;
  nachher:=0;
  in0:=true;
  out0:=true;

  case sizeof(wgen) of
    sizeof(single):
      for i:=smin to smax do begin
        algo.laden(invers,pSingle(@(werte[i*sstep])),pstep);
        algo.summen(true,vorher,in0);
        algo.ausfuehren;
        algo.summen(false,nachher,out0);
        algo.speichern(invers,pSingle(@(werte[i*sstep])),pstep);
      end;
    sizeof(double):
      for i:=smin to smax do begin
        algo.laden(invers,pDouble(@(werte[i*sstep])),pstep);
        algo.summen(true,vorher,in0);
        algo.ausfuehren;
        algo.summen(false,nachher,out0);
        algo.speichern(invers,pDouble(@(werte[i*sstep])),pstep);
      end;
    sizeof(extended):
      for i:=smin to smax do begin
        algo.laden(invers,pExtended(@(werte[i*sstep])),pstep);
        algo.summen(true,vorher,in0);
        algo.ausfuehren;
        algo.summen(false,nachher,out0);
        algo.speichern(invers,pExtended(@(werte[i*sstep])),pstep);
      end;
  end{of case};

  if (nachher=0) and (vorher=0) then pvFehler:=0
  else pvFehler:=abs(nachher-vorher)/(nachher+vorher);

  gibAus(inttostr(byte(senkrecht))+' '+inttostr(byte(invers))+' (Parseval-Fehler = '+floattostr(pvFehler)+') ... nämlich von '+floattostr(vorher)+' zu '+floattostr(nachher),1);
  if in0 then gibAus('Nur Nullen im Input der FFT!',1);
  if out0 then gibAus('Nur Nullen im Output der FFT!',1);
  result:=not (in0 or out0);
end;

procedure tLLWerte.schreibeWert(var f: textfile; x,y: longint);
begin
  schreibeWert(f,x,y,werte[x+y*params.xsteps]);
end;

procedure tLLWerte.schreibeWert(var f: textfile; x,y,wert: extended);
var
  xa,ta: extended;
begin
  if params.xstop=params.xstart then
    xa:=params.xstop
  else
    xa:=params.xstart+x/(params.xsteps-1)*(params.xstop-params.xstart);
  if params.tstop=params.tstart then
    ta:=params.tstop
  else
    ta:=params.tstart+y/(params.tsiz-1)*(params.tstop-params.tstart);
  writeln(f,floattostr(xa)+' '+floattostr(ta)+' '+floattostr(wert));
end;

procedure tLLWerte.schreibeWertIntegriert(var f: textfile; i: longint; hor: boolean);
var
  j:   longint;
  tmp: extended;
begin
  tmp:=0;
  if hor then begin
    for j:=0 to params.xsteps-1 do
      tmp:=tmp+werte[j+i*params.xsteps];
    schreibeWert(f,(params.xsteps-1)/2,i,tmp);
  end
  else begin
    for j:=0 to params.tsiz-1 do
      tmp:=tmp+werte[i+j*params.xsteps];
    schreibeWert(f,i,(params.tsiz-1)/2,tmp);
  end;
end;

procedure tLLWerte.erzeugeBinning(senkrecht,linien: boolean; x0,dx: extended; s: string);
var
  f:     textfile;
  i:     longint;
  sum,x: extended;
begin
  assignfile(f,s);
  rewrite(f);
  sum:=0;
  while x0<0 do
    x0:=x0+dx;
  for i:=0 to params.xsteps*params.tsiz-1 do
    if i+1>x0 then begin
      sum:=sum+werte[i]*(x0-i);
      if senkrecht then
        x:=x0/(params.tsiz-1)*(params.tstop-params.tstart)+params.tstart
      else
        x:=x0/(params.xsteps-1)*(params.xstop-params.xstart)+params.xstart;
      writeln(f,floattostr(x)+' '+floattostr(sum/dx));
      if linien then begin
        writeln(f,floattostr(x)+' 0');
        writeln(f)
      end;
      sum:=werte[i]*(i+1-x0);
      x0:=x0+dx;
    end
    else
      sum:=sum+werte[i];
  closefile(f);
end;

procedure tLLWerte.spiegle;
begin
  spiegle(0,params.tsiz-1);
end;

procedure tLLWerte.spiegle(tmin,tmax: longint);
var
  i,j:    longint;
  tmp:    wgen;
begin
  for i:=tmin to tmax do
    for j:=0 to params.xsteps div 2 -1 do begin
      tmp:=werte[j+i*params.xsteps];
      werte[j+i*params.xsteps]:=werte[params.xsteps-j-1+i*params.xsteps];
      werte[params.xsteps-j-1+i*params.xsteps]:=tmp;
    end;
end;

procedure tLLWerte.fft2dNachbearbeitungA(nb: tFFTDatenordnung);
var
  i: longint;
begin
  case NB of
    doResIms,doResSmi: ;
    doBetr: begin
      werte[0]:=abs(extended(werte[0]));                                       // rein reell
      for i:=1 to params.xsteps div 2 -1 do begin                             // unterste Zeile ist reell in t
        werte[i]:=sqrt(sqr(extended(werte[i]))+sqr(extended(werte[params.xsteps-i])));
        werte[params.xsteps-i]:=werte[i];
      end;
      for i:=1 to params.xsteps div 2 -1 do begin                             // mittlere Zeile ist reell in t
        werte[i+(params.tsiz div 2)*params.xsteps]:=sqrt(sqr(extended(werte[i+(params.tsiz div 2)*params.xsteps]))+sqr(extended(werte[params.xsteps-i+(params.tsiz div 2)*params.xsteps])));
        werte[params.xsteps-i+(params.tsiz div 2)*params.xsteps]:=werte[i+(params.tsiz div 2)*params.xsteps];
      end;
      werte[params.xsteps div 2]:=abs(extended(werte[params.xsteps div 2])); // rein reell
      werte[params.tsiz div 2]:=abs(extended(werte[params.tsiz div 2]));     // rein reell
      for i:=1 to params.tsiz div 2 -1 do begin                               // linkeste Spalte ist reell in x
        werte[i*params.xsteps]:=sqrt(sqr(extended(werte[i*params.xsteps]))+sqr(extended(werte[(params.tsiz-i)*params.xsteps])));
        werte[(params.tsiz-i)*params.xsteps]:=werte[i*params.xsteps];
      end;
      for i:=1 to params.tsiz div 2 -1 do begin                               // mittlere Spalte ist reell in x
        werte[params.xsteps div 2 + i*params.xsteps]:=sqrt(sqr(extended(werte[params.xsteps div 2 + i*params.xsteps]))+sqr(extended(werte[params.xsteps div 2 + (params.tsiz-i)*params.xsteps])));
        werte[params.xsteps div 2 + (params.tsiz-i)*params.xsteps]:=werte[params.xsteps div 2 + i*params.xsteps];
      end;
    end;
    doBetrQdr: begin
      werte[0]:=sqr(extended(werte[0]));                                       // rein reell
      for i:=1 to params.xsteps div 2 -1 do begin                             // unterste Zeile ist reell in t
        werte[i]:=sqr(extended(werte[i]))+sqr(extended(werte[params.xsteps-i]));
        werte[params.xsteps-i]:=werte[i];
      end;
      for i:=1 to params.xsteps div 2 -1 do begin                             // mittlere Zeile ist reell in t
        werte[i+(params.tsiz div 2)*params.xsteps]:=sqr(extended(werte[i+(params.tsiz div 2)*params.xsteps]))+sqr(extended(werte[params.xsteps-i+(params.tsiz div 2)*params.xsteps]));
        werte[params.xsteps-i+(params.tsiz div 2)*params.xsteps]:=werte[i+(params.tsiz div 2)*params.xsteps];
      end;
      werte[params.xsteps div 2]:=sqr(extended(werte[params.xsteps div 2])); // rein reell
      werte[params.tsiz div 2]:=sqr(extended(werte[params.tsiz div 2]));     // rein reell
      for i:=1 to params.tsiz div 2 -1 do begin                               // linkeste Spalte ist reell in x
        werte[i*params.xsteps]:=sqr(extended(werte[i*params.xsteps]))+sqr(extended(werte[(params.tsiz-i)*params.xsteps]));
        werte[(params.tsiz-i)*params.xsteps]:=werte[i*params.xsteps];
      end;
      for i:=1 to params.tsiz div 2 -1 do begin                               // mittlere Spalte ist reell in x
        werte[params.xsteps div 2 + i*params.xsteps]:=sqr(extended(werte[params.xsteps div 2 + i*params.xsteps]))+sqr(extended(werte[params.xsteps div 2 + (params.tsiz-i)*params.xsteps]));
        werte[params.xsteps div 2 + (params.tsiz-i)*params.xsteps]:=werte[params.xsteps div 2 + i*params.xsteps];
      end;
    end;
  end{of case};
end;

procedure tLLWerte.fft2dNachbearbeitungB(xmin,xmax: longint; nb: tFFTDatenordnung);
var
  i,j:    longint;
begin // bearbeitet nur den Hauptteil (außer erster und mittlerer Zeile/Spalte) nach!
  case nb of
    doBetr: begin
      for i:=xmin to xmax do
        for j:=1 to params.tsiz div 2 -1 do begin
          werte[i+j*params.xsteps]:=
            sqrt(sqr(extended(werte[i+j*params.xsteps]-werte[params.xsteps-i+(params.tsiz-j)*params.xsteps]))   // Re^2
                +sqr(extended(werte[params.xsteps-i+j*params.xsteps]+werte[i+(params.tsiz-j)*params.xsteps]))); // Im^2
          werte[params.xsteps-i+j*params.xsteps]:=werte[i+j*params.xsteps];
          werte[i+(params.tsiz-j)*params.xsteps]:=werte[i+j*params.xsteps];
          werte[params.xsteps-i+(params.tsiz-j)*params.xsteps]:=werte[i+j*params.xsteps];
        end;
    end;
    doBetrQdr: begin
      for i:=xmin to xmax do
        for j:=1 to params.tsiz div 2 -1 do begin
          werte[i+j*params.xsteps]:=
              sqr(extended(werte[i+j*params.xsteps]-werte[params.xsteps-i+(params.tsiz-j)*params.xsteps]))  // Re^2
            +sqr(extended(werte[params.xsteps-i+j*params.xsteps]+werte[i+(params.tsiz-j)*params.xsteps])); // Im^2
          werte[params.xsteps-i+j*params.xsteps]:=werte[i+j*params.xsteps];
          werte[i+(params.tsiz-j)*params.xsteps]:=werte[i+j*params.xsteps];
          werte[params.xsteps-i+(params.tsiz-j)*params.xsteps]:=werte[i+j*params.xsteps];
        end;
    end;
  end{of case}
end;

procedure tLLWerte.holeRam;
begin
  holeRam(1);
end;

procedure tLLWerte.holeRam(ausgaben: byte);
begin
  holeRam(ausgaben,false);
end;

procedure tLLWerte.holeRam(ausgaben: byte; gemaeszTXMinMax: boolean);
var
  Zeit:  extended;
  br,ho: longint;
begin
  Zeit:=now;
  if gemaeszTXMinMax then begin
    br:=params.xsteps_;
    ho:=params.tsiz_;
  end
  else begin
    br:=params.xsteps;
    ho:=params.tsiz;
  end;
  if (ausgaben and __ausgabenMaske) <> 0 then
    gibAus('Fordere '+inttostr(floor(ho*br*sizeof(wgen)/1024/1024))+' MB RAM an ('+inttostr(br)+' x-Schritte mal '+inttostr(ho)+' t-Schritte; bisher '+inttostr(belegterSpeicher div 1024)+' MB belegt). ...',ausgaben);

  setlength(werte,br*ho);
  if (ausgaben and __ausgabenMaske) <> 0 then
    gibAus('... fertig '+timetostr(now-Zeit),ausgaben);
end;

function tLLWerte.zuPixelWerten(whoehe,wbreite,xpmi,xmi,tmi: longint; xz,yz: extended; pPWerte: pTExtendedArray; pPAnzahlen: pTLongintArray): boolean;
var
  i,j,k,l,
    xv,xb,tv,tb: longint;
  b:             boolean;
begin
  result:=false;
  for i:=0 to length(pPWerte^)-1 do begin
    pPWerte^[i]:=0;
    pPAnzahlen^[i]:=0;
  end;
  b:=false;
  for j:=0 to whoehe-1 do
    for i:=0 to wbreite-1 do begin
      xv:=min(params.xsteps-1,max(0,ceil((i+max(0,xpmi)-1/2)/xz+xmi)));
      xb:=min(params.xsteps-1,max(0,ceil((i+max(0,xpmi)+1/2)/xz+xmi-1)));
      tv:=min(params.tsiz-1,max(0,ceil((j-1/2)/yz+tmi)));
      tb:=min(params.tsiz-1,max(0,ceil((j+1/2)/yz+tmi-1)));
      if xv>xb then begin
        if (i>0) or (xpmi>0) then
          dec(xv)
        else
          inc(xb);
      end;
      if tv>tb then begin
        if j>0 then
          dec(tv)
        else
          inc(tb);
      end;
      if (xv>xb) or (tv>tb) then begin
        gibAus('Keine Inputwerte für Position '+inttostr(i)+':'+inttostr(j)+'!',1);
        exit;
      end;
      pPanzahlen^[j*wbreite+i]:=(xb-xv+1)*(tb-tv+1);
      for k:=xv to xb do
        for l:=tv to tb do begin
          b:=b or (werte[k+l*params.xsteps]<>0);
          pPWerte^[j*wbreite+i]:=pPWerte^[j*wbreite+i]+min(params.maxW,max(params.minW,werte[k+l*params.xsteps]));
        end;
    end;
  if not b then
    gibAus('Habe nur Nullen im Input!',1);
  result:=true;
end;

function tLLWerte.findeSchwellwerte(xmi,xma,tmi,tma: longint; Schw: extended): tExtPointArray;
var
  i,j,k,l,m,vz: longint;
  dx,dy,x0,y0:  extended;
begin
  setlength(result,0);
  gibAus('Schwellwerte finden ('+inttostr(xmi)+'-'+inttostr(xma)+') '+floattostr(Schw)+' ...',1);
  i:=0;
  dx:=(params.xstop-params.xstart)/(params.xsteps-1);
  x0:=params.xstart;
  dy:=(params.tstop-params.tstart)/(params.tsiz-1);
  y0:=params.tstart;
  for j:=xmi to xma do begin
    for k:=tmi to tma do begin
      vz:=0;
      for l:=0 to 1 do
        for m:=0 to 1 do
          vz:=
              vz or (byte(werte[j-l+(k-m)*params.xsteps]>=Schw) shl 1) or
                    byte(werte[j-l+(k-m)*params.xsteps]<=Schw);
      if vz=3 then begin
        if i>=length(result) then
          setlength(result,length(result)+Speicherhappen);
        result[i]['x']:=j*dx+x0;
        result[i]['y']:=k*dy+y0;
        inc(i);
      end;
    end;
    if (xma-j) and $ff = 0 then
      gibAus('x = '+inttostr(j)+' ('+inttostr(xmi)+'-'+inttostr(xma)+')',1);
  end;
  gibAus('... fertig!',1);
  setlength(result,i);
end;

procedure tLLWerte.integriereSingle(qu: pTLLWerteSingle; xmi,xma,tmi,tma,xof,tof: longint; richtung: tIntegrationsRichtung);
var
  i,j:        longint;
  int,faktor: extended;
begin
  case richtung of
    irHorizontal: begin
      faktor:=(qu^.params.xstop-qu^.params.xstart)/(qu^.params.xsteps-1);
      for i:=tmi to tma do begin
        int:=0;
        for j:=0 to xmi-1 do
          int:=int+qu^.werte[j + i*qu^.params.xsteps];
        for j:=xmi to xma do begin
          int:=int+qu^.werte[j + i*qu^.params.xsteps];
          werte[j-xof + (i-tof)*params.xsteps]:=int*faktor;
        end;
      end;
    end;
    irEinfall: begin
      faktor:=
        sqrt(
          sqr((qu^.params.xstop-qu^.params.xstart)/(qu^.params.xsteps-1)) +
          sqr((qu^.params.tstop-qu^.params.tstart)/(qu^.params.tsiz-1)));
      gibAus('dx = '+floattostr((qu^.params.xstop-qu^.params.xstart)/(qu^.params.xsteps-1)),1);
      gibAus('dt = '+floattostr((qu^.params.tstop-qu^.params.tstart)/(qu^.params.tsiz-1)),1);
      for i:=tmi to tma do begin  // von links eintretendes (inkl. Ecke links unten)
        int:=0;
        for j:=1 to min(xmi,i) do
          int:=int+qu^.werte[xmi-j + (i-j)*qu^.params.xsteps];
        for j:=0 to min(tma-i,xma-xmi) do begin
          int:=int+qu^.werte[xmi+j + (i+j)*qu^.params.xsteps];
          werte[j+xmi-xof + (i+j-tof)*params.xsteps]:=int*faktor;
        end;
      end;
      for i:=xmi+1 to xma do begin  // von unten eintretendes (exkl. Ecke links unten)
        int:=0;
        for j:=1 to min(tmi,i) do
          int:=int+qu^.werte[i-j + (tmi-j)*qu^.params.xsteps];
        for j:=0 to min(tma-tmi,xma-i) do begin
          int:=int+qu^.werte[i+j + (tmi+j)*qu^.params.xsteps];
          werte[i+j-xof + (tmi+j-tof)*params.xsteps]:=int*faktor;
        end;
      end;
    end;
    irAusfall: begin
      faktor:=
        sqrt(
          sqr((qu^.params.xstop-qu^.params.xstart)/(qu^.params.xsteps-1)) +
          sqr((qu^.params.tstop-qu^.params.tstart)/(qu^.params.tsiz-1)));
      gibAus('dx = '+floattostr((qu^.params.xstop-qu^.params.xstart)/(qu^.params.xsteps-1)),1);
      gibAus('dt = '+floattostr((qu^.params.tstop-qu^.params.tstart)/(qu^.params.tsiz-1)),1);
      for i:=tmi to tma do begin  // nach links austretendes (inkl. Ecke links oben)
        int:=0;
        for j:=1 to min(xmi,qu^.params.tsiz-1-i) do
          int:=int+qu^.werte[xmi-j + (i+j)*qu^.params.xsteps];
        for j:=0 to min(i-tmi,xma-xmi) do begin
          int:=int+qu^.werte[xmi+j + (i-j)*qu^.params.xsteps];
          werte[j+xmi-xof + (i-j-tof)*params.xsteps]:=int*faktor;
        end;
      end;
      for i:=xmi+1 to xma do begin  // nach oben austretendes (exkl. Ecke links oben)
        int:=0;
        for j:=1 to min(qu^.params.tsiz-1-tma,i) do
          int:=int+qu^.werte[i-j + (tma+j)*qu^.params.xsteps];
        for j:=0 to min(tma-tmi,xma-i) do begin
          int:=int+qu^.werte[i+j + (tma-j)*qu^.params.xsteps];
          werte[i+j-xof + (tma-j-tof)*params.xsteps]:=int*faktor;
        end;
      end;
    end;
  end{of case};
end;

procedure tLLWerte.integriereDouble(qu: pTLLWerteDouble; xmi,xma,tmi,tma,xof,tof: longint; richtung: tIntegrationsRichtung);
var
  i,j:        longint;
  int,faktor: extended;
begin
  case richtung of
    irHorizontal: begin
      faktor:=(qu^.params.xstop-qu^.params.xstart)/(qu^.params.xsteps-1);
      for i:=tmi to tma do begin
        int:=0;
        for j:=0 to xmi-1 do
          int:=int+qu^.werte[j + i*qu^.params.xsteps];
        for j:=xmi to xma do begin
          int:=int+qu^.werte[j + i*qu^.params.xsteps];
          werte[j-xof + (i-tof)*params.xsteps]:=int*faktor;
        end;
      end;
    end;
    irEinfall: begin
      faktor:=
        sqrt(
          sqr((qu^.params.xstop-qu^.params.xstart)/(qu^.params.xsteps-1)) +
          sqr((qu^.params.tstop-qu^.params.tstart)/(qu^.params.tsiz-1)));
      gibAus('dx = '+floattostr((qu^.params.xstop-qu^.params.xstart)/(qu^.params.xsteps-1)),1);
      gibAus('dt = '+floattostr((qu^.params.tstop-qu^.params.tstart)/(qu^.params.tsiz-1)),1);
      for i:=tmi to tma do begin  // von links eintretendes (inkl. Ecke links unten)
        int:=0;
        for j:=1 to min(xmi,i) do
          int:=int+qu^.werte[xmi-j + (i-j)*qu^.params.xsteps];
        for j:=0 to min(tma-i,xma-xmi) do begin
          int:=int+qu^.werte[xmi+j + (i+j)*qu^.params.xsteps];
          werte[j+xmi-xof + (i+j-tof)*params.xsteps]:=int*faktor;
        end;
      end;
      for i:=xmi+1 to xma do begin  // von unten eintretendes (exkl. Ecke links unten)
        int:=0;
        for j:=1 to min(tmi,i) do
          int:=int+qu^.werte[i-j + (tmi-j)*qu^.params.xsteps];
        for j:=0 to min(tma-tmi,xma-i) do begin
          int:=int+qu^.werte[i+j + (tmi+j)*qu^.params.xsteps];
          werte[i+j-xof + (tmi+j-tof)*params.xsteps]:=int*faktor;
        end;
      end;
    end;
    irAusfall: begin
      faktor:=
        sqrt(
          sqr((qu^.params.xstop-qu^.params.xstart)/(qu^.params.xsteps-1)) +
          sqr((qu^.params.tstop-qu^.params.tstart)/(qu^.params.tsiz-1)));
      gibAus('dx = '+floattostr((qu^.params.xstop-qu^.params.xstart)/(qu^.params.xsteps-1)),1);
      gibAus('dt = '+floattostr((qu^.params.tstop-qu^.params.tstart)/(qu^.params.tsiz-1)),1);
      for i:=tmi to tma do begin  // nach links austretendes (inkl. Ecke links oben)
        int:=0;
        for j:=1 to min(xmi,qu^.params.tsiz-1-i) do
          int:=int+qu^.werte[xmi-j + (i+j)*qu^.params.xsteps];
        for j:=0 to min(i-tmi,xma-xmi) do begin
          int:=int+qu^.werte[xmi+j + (i-j)*qu^.params.xsteps];
          werte[j+xmi-xof + (i-j-tof)*params.xsteps]:=int*faktor;
        end;
      end;
      for i:=xmi+1 to xma do begin  // nach oben austretendes (exkl. Ecke links oben)
        int:=0;
        for j:=1 to min(qu^.params.tsiz-1-tma,i) do
          int:=int+qu^.werte[i-j + (tma+j)*qu^.params.xsteps];
        for j:=0 to min(tma-tmi,xma-i) do begin
          int:=int+qu^.werte[i+j + (tma-j)*qu^.params.xsteps];
          werte[i+j-xof + (tma-j-tof)*params.xsteps]:=int*faktor;
        end;
      end;
    end;
  end{of case};
end;

procedure tLLWerte.integriereExtended(qu: pTLLWerteDouble; xmi,xma,tmi,tma,xof,tof: longint; richtung: tIntegrationsRichtung);
var
  i,j:        longint;
  int,faktor: extended;
begin
  case richtung of
    irHorizontal: begin
      faktor:=(qu^.params.xstop-qu^.params.xstart)/(qu^.params.xsteps-1);
      for i:=tmi to tma do begin
        int:=0;
        for j:=0 to xmi-1 do
          int:=int+qu^.werte[j + i*qu^.params.xsteps];
        for j:=xmi to xma do begin
          int:=int+qu^.werte[j + i*qu^.params.xsteps];
          werte[j-xof + (i-tof)*params.xsteps]:=int*faktor;
        end;
      end;
    end;
    irEinfall: begin
      faktor:=
        sqrt(
          sqr((qu^.params.xstop-qu^.params.xstart)/(qu^.params.xsteps-1)) +
          sqr((qu^.params.tstop-qu^.params.tstart)/(qu^.params.tsiz-1)));
      gibAus('dx = '+floattostr((qu^.params.xstop-qu^.params.xstart)/(qu^.params.xsteps-1)),1);
      gibAus('dt = '+floattostr((qu^.params.tstop-qu^.params.tstart)/(qu^.params.tsiz-1)),1);
      for i:=tmi to tma do begin  // von links eintretendes (inkl. Ecke links unten)
        int:=0;
        for j:=1 to min(xmi,i) do
          int:=int+qu^.werte[xmi-j + (i-j)*qu^.params.xsteps];
        for j:=0 to min(tma-i,xma-xmi) do begin
          int:=int+qu^.werte[xmi+j + (i+j)*qu^.params.xsteps];
          werte[j+xmi-xof + (i+j-tof)*params.xsteps]:=int*faktor;
        end;
      end;
      for i:=xmi+1 to xma do begin  // von unten eintretendes (exkl. Ecke links unten)
        int:=0;
        for j:=1 to min(tmi,i) do
          int:=int+qu^.werte[i-j + (tmi-j)*qu^.params.xsteps];
        for j:=0 to min(tma-tmi,xma-i) do begin
          int:=int+qu^.werte[i+j + (tmi+j)*qu^.params.xsteps];
          werte[i+j-xof + (tmi+j-tof)*params.xsteps]:=int*faktor;
        end;
      end;
    end;
    irAusfall: begin
      faktor:=
        sqrt(
          sqr((qu^.params.xstop-qu^.params.xstart)/(qu^.params.xsteps-1)) +
          sqr((qu^.params.tstop-qu^.params.tstart)/(qu^.params.tsiz-1)));
      gibAus('dx = '+floattostr((qu^.params.xstop-qu^.params.xstart)/(qu^.params.xsteps-1)),1);
      gibAus('dt = '+floattostr((qu^.params.tstop-qu^.params.tstart)/(qu^.params.tsiz-1)),1);
      for i:=tmi to tma do begin  // nach links austretendes (inkl. Ecke links oben)
        int:=0;
        for j:=1 to min(xmi,qu^.params.tsiz-1-i) do
          int:=int+qu^.werte[xmi-j + (i+j)*qu^.params.xsteps];
        for j:=0 to min(i-tmi,xma-xmi) do begin
          int:=int+qu^.werte[xmi+j + (i-j)*qu^.params.xsteps];
          werte[j+xmi-xof + (i-j-tof)*params.xsteps]:=int*faktor;
        end;
      end;
      for i:=xmi+1 to xma do begin  // nach oben austretendes (exkl. Ecke links oben)
        int:=0;
        for j:=1 to min(qu^.params.tsiz-1-tma,i) do
          int:=int+qu^.werte[i-j + (tma+j)*qu^.params.xsteps];
        for j:=0 to min(tma-tmi,xma-i) do begin
          int:=int+qu^.werte[i+j + (tma-j)*qu^.params.xsteps];
          werte[i+j-xof + (tma-j-tof)*params.xsteps]:=int*faktor;
        end;
      end;
    end;
  end{of case};
end;

procedure tLLWerte.gauszFit(amplituden,breiten,positionen,ueberlappe,hintergruende: pTLLWerteExtended; von,bis: longint; senkrecht: boolean; fensterBreite,maxBreite,maxVerschiebung: extended; positionsMitten: tExtendedArray);
var
  i,j,ii,zaehl,                    // Zähler
    qpSchritt,qsSchritt,           // Schrittlängen in der Quelle
    zpSchritt,zdSchritt,           // Schrittlängen im Ziel
    pLen,offset,                   // Datenlänge und Offset in der Quelle
    wiMin,wiMax,                   // momentane Fensterränder (für den Zählindex) in der Quelle
    verbesserung,                  // <0: schlechter, 1: besser, 2: viel besser
    aParams,nParams:     longint;  // Index der alten und neuen Parameter
  parameter:             array[0..1,boolean,0..3] of extended;
    // dim0 nummeriert Parametersätze
    // dim1: Ableitung (true) oder Werte (false)
    // dim2: position, 1/breite, amplitude, hintergrund
  fehlers:               array[0..1] of extended;
  t0,t1,t2,t3,                     // Zwischenergebnisse
    qdrSumm,                       // Quadratesumme im betrachteten Fenster
    pMi,pMa,                       // Achsenenden in Datenrichtung
    schrittFaktor:       extended;
  ignBr,ignVersch:       boolean;  // Breite/Verschiebung am Anschlag!
begin
  if senkrecht then begin
    qpSchritt:=params.xsteps; // Schritt in der Quelle parallel zur Fit-Richtung
    qsSchritt:=1;             // Schritt in der Quelle senkrecht zur Fit-Richtung
    zpSchritt:=params.xsteps; // Schritt im Ziel in positionsMitten-Richtung
    zdSchritt:=1;             // Schritt im Ziel in Daten-Richtung (= senkrecht zur Fit-Richtung)
    pMi:=params.tstart;
    pMa:=params.tstop;
    pLen:=params.tsiz;
  end
  else begin
    qsSchritt:=params.xsteps;
    qpSchritt:=1;
    zdSchritt:=length(positionsMitten);
    zpSchritt:=1;
    pMi:=params.xstart;
    pMa:=params.xstop;
    pLen:=params.xsteps;
  end;
  maxBreite:=1/maxBreite;

  for i:=von to bis do
    for j:=0 to length(positionsMitten)-1 do begin
      offset:=i*qsSchritt;
      aParams:=0;

      wiMin:=0;
      wiMax:=pLen-1;
      if fensterBreite>0 then begin
        wiMin:=max(wiMin,round(positionsMitten[j] - fensterBreite / 2));
        wiMax:=min(wiMax,round(positionsMitten[j] + fensterBreite / 2));
      end;

      qdrSumm:=0;
      t0:=0;
      t1:=0;
      t2:=0;
      schrittFaktor:=0;
      for ii:=wiMin to wiMax do begin
        t3:=werte[offset+qpSchritt*ii];
        t3:=abs(t3);
        qdrSumm:=qdrSumm + sqr(t3);
        t0:=t0 + t3;
        t1:=t1 + ii*t3;
        t2:=t2 + sqr(ii)*t3;
        if t3>schrittFaktor then schrittFaktor:=t3;
      end;

      t1:=t1/t0;
      t2:=t2/t0;

      parameter[aParams,false,0]:=t1;                     // Erwartungswert
      parameter[aParams,false,2]:=schrittFaktor;          // Maximalwert
      parameter[aParams,false,1]:=parameter[aParams,false,2]/t0*sqrt(pi);
//      parameter[aParams,false,1]:=1/sqrt(2*(t2-sqr(t1))); // beinhalted Standardabweichung  <-- schlechter!!
      parameter[aParams,false,3]:=0;

      if (maxBreite>=0) and (parameter[aParams,false,1]<2*maxBreite) then
        parameter[aParams,false,1]:=maxBreite*2;

      nParams:=aParams;
      {$DEFINE gauszFitBerechneWerte}
      {$I gauszFit.inc} // Werte + Gradienten berechnen
      {$UNDEF}

      schrittFaktor:=1;
      if maxVerschiebung>=0 then
        schrittFaktor:=min(schrittFaktor,maxVerschiebung/max(abs(parameter[nParams,true,0]),1e-50));
      if maxBreite>=0 then
        schrittFaktor:=min(schrittFaktor,abs(parameter[nParams,false,2]-maxBreite)/max(abs(parameter[nParams,true,1]),1e-50));

      zaehl:=0;
      repeat
        nParams:=(aParams+1) mod length(parameter);
        for ii:=0 to length(parameter[aParams,false])-1 do
          parameter[nParams,false,ii]:=
            parameter[aParams,false,ii] - schrittFaktor * parameter[aParams,true,ii]; // * byte(ii<3);

        ignVersch:=false;
        if maxVerschiebung>0 then begin
          if parameter[nParams,false,0]-positionsMitten[j] > maxVerschiebung then begin
            parameter[nParams,false,0]:=positionsMitten[j] + maxVerschiebung;
            ignVersch:=true;
          end;
          if positionsMitten[j]-parameter[nParams,false,0] > maxVerschiebung then begin
            parameter[nParams,false,0]:=positionsMitten[j] - maxVerschiebung;
            ignVersch:=true;
          end;
        end;

        ignBr:=false;
        if maxBreite>0 then
          if parameter[nParams,false,1] < maxBreite then begin
            parameter[nParams,false,1]:=maxBreite;
            ignBr:=true;
          end;

        {$DEFINE gauszFitBerechneWerte}
        {$I gauszFit.inc} // Werte + Gradienten berechnen
        {$UNDEF}

        if fehlers[aParams]>2*fehlers[nParams] then begin
          verbesserung:=2;
          schrittFaktor:=schrittFaktor * 2;
          aParams:=nParams;
        end
        else if fehlers[aParams]>fehlers[nParams] then begin
          verbesserung:=1;
          schrittFaktor:=schrittFaktor * 1.3;
          aParams:=nParams;
        end
        else begin
          if verbesserung>0 then
            verbesserung:=-1
          else
            dec(verbesserung);
          schrittFaktor:=schrittFaktor * 0.5;
        end;

        if schrittFaktor<1e-50 then
          fehler('Sehr kleiner Schrittfaktor ('+floattostr(schrittFaktor)+') beim Fitten!');

        inc(zaehl);
        {$IFDEF gauszFitStatus}
        if (zaehl mod 10000)=0 then
          gibAus(
            floattostr(fehlers[aParams])+' '+
            floattostr(qdrSumm)+' '+
            inttostr(byte((verbesserung<-10) or (fehlers[aParams]*100<qdrSumm)))+' '+
            inttostr(byte(ignVersch or (abs(schrittFaktor*parameter[aParams,true,0])<1)))+' '+
            inttostr(byte(ignBr or (abs(schrittFaktor*parameter[aParams,true,1]/parameter[aParams,false,1])<0.01)))+' '+
            inttostr(byte(abs(schrittFaktor*parameter[aParams,true,2]/parameter[aParams,false,2])<0.01))+' '+
            inttostr(byte(abs(schrittFaktor*parameter[aParams,true,3]/max(abs(parameter[aParams,false,3]),1e-10))<0.01))+' '+
            inttostr(verbesserung),
            3);
        {$ENDIF}

      until ((zaehl>=100000) or (verbesserung<-10) or (fehlers[aParams]*10<qdrSumm)) and
            (ignVersch or (abs(schrittFaktor*parameter[aParams,true,0])<1)) and
            (ignBr or (abs(schrittFaktor*parameter[aParams,true,1]/parameter[aParams,false,1])<0.01)) and
            (abs(schrittFaktor*parameter[aParams,true,2]/parameter[aParams,false,2])<0.01) and
            (abs(schrittFaktor*parameter[aParams,true,3]/max(abs(parameter[aParams,false,3]),1e-10))<0.01);
      {$IFDEF gauszFitStatus}
      gibAus('',3);
      {$ENDIF}

      if assigned(positionen) then
        positionen^.werte[j*zpSchritt + i*zdSchritt]:=pMi+(pMa-pMi)/(pLen-1)*parameter[aParams,false,0];
      if assigned(breiten) then
        breiten^.werte[j*zpSchritt + i*zdSchritt]:=(pMa-pMi)/(pLen-1)/parameter[aParams,false,1];
      if assigned(amplituden) then
        amplituden^.werte[j*zpSchritt + i*zdSchritt]:=parameter[aParams,false,2];
      if assigned(hintergruende) then
        hintergruende^.werte[j*zpSchritt + i*zdSchritt]:=parameter[aParams,false,3];
      if assigned(ueberlappe) then
        ueberlappe^.werte[j*zpSchritt + i*zdSchritt]:=fehlers[aParams]/qdrSumm;
    end;
end;

// tWavelet ********************************************************************

function tWavelet.setzeTyp(s: string): boolean;
begin
  result:=true;
  if s='Sin2' then begin
    typ:=wtSin2;
    exit;
  end;
  if s='Frequenzfenster' then begin
    typ:=wtFrequenzfenster;
    exit;
  end;
  gibAus('Kenne Wavelettyp '''+s+''' nicht!',3);
  result:=false;
end;

function tWavelet.berechneWerte: boolean;
var
  i:          longint;
  tmp:        extended;
  nur0:       boolean;
  tmpFFTAlgo: tFFTAlgorithmus;
begin
  result:=false;
  werte.params.xsteps:=2;
  werte.params.transformationen.xstart:=0;
  werte.params.transformationen.xstop:=1;
  werte.params.transformationen.tstart:=0;
  werte.params.transformationen.tstop:=1;

  if mitFFT then begin
    if round(power(2,round(ln(werte.params.tsiz)/ln(2))))<>werte.params.tsiz then begin
      gibAus('Waveletlänge muss eine Zweierpotenz sein bei Faltung mittels FFT, '+inttostr(werte.params.tsiz)+' ist das nicht!',3);
      exit;
    end;
    werte.holeRam(1);
    case typ of
      wtSin2: begin
        hlen:=round(tfwhm);
        for i:=0 to hlen do begin
          tmp:=sqr(cos(i*pi/2/tfwhm));
          if i>0 then begin
            werte.werte[(werte.params.tsiz-i)*2]:=  tmp*cos(-i*freq*2*pi);
            werte.werte[(werte.params.tsiz-i)*2+1]:=tmp*sin(-i*freq*2*pi);
          end;
          werte.werte[2*i]:=  tmp*cos(i*freq*2*pi);
          werte.werte[2*i+1]:=tmp*sin(i*freq*2*pi);
        end;
        for i:=hlen+1 to werte.params.tsiz-1-hlen do begin
          werte.werte[2*i]:=0;
          werte.werte[2*i+1]:=0;
        end;
        gibAus(inttostr(werte.params.xsteps)+' '+inttostr(werte.params.tsiz)+' '+inttostr(length(werte.werte)),1);
        tmpFFTAlgo:=createFFTAlgorithmus(2*hlen,doRes,doResIms);
        werte.fft(true,false,tmpFFTAlgo,nil,pvFehler);
        tmpFFTAlgo.free;
      end;
      wtFrequenzfenster: begin
        hlen:=werte.params.tsiz div 2;
        for i:=0 to hlen do begin
          werte.werte[2*i]:=Byte(tfwhm*Abs(i/werte.params.tsiz-freq)<=1); // Re=1 <=> |f-f_Mitte| < Laenge/T_fwhm
          werte.werte[2*i+1]:=0; // Dummy-Wavelet
        end;
        for i:=hlen+1 to 2*hlen-1 do begin
          werte.werte[2*i]:=0;   // Ims = 0
          werte.werte[2*i+1]:=0; // Dummy-Wavelet
        end;
      end;
    end{of case};
  end
  else begin
    if typ <> wtSin2 then begin
      gibAus('Ich kann nur das ''Sin2''-Wavelet im Zeitbereich definieren!',3);
      exit;
    end;
    hlen:=round(tfwhm);
    werte.params.tsiz:=2*hlen+1;
    setlength(werte.werte,werte.params.xsteps*werte.params.tsiz);
    for i:=-hlen to hlen do begin
      tmp:=sqr(cos(i*pi/2/tfwhm));
      werte.werte[2*(i+hlen)]:=  tmp*cos(i*freq*2*pi);
      werte.werte[2*(i+hlen)+1]:=tmp*sin(i*freq*2*pi);
    end;
  end;
  nur0:=true;
  for i:=0 to length(werte.werte)-1 do
    nur0:=nur0 and (werte.werte[i]=0);
  if nur0 then gibAus('Das Wavelet hat nur Nullen!',1);
  result:=not nur0;
  werte.params.refreshKnownValues;
end;

constructor tWavelet.create;
var
  ps: tExtrainfos;
begin
  inherited create;
  ps:=tExtrainfos.create;
  werte:=tLLWerteDouble.create(ps);
  pvFehler:=0;
end;

destructor tWavelet.destroy;
begin
  werte.params.Free;
  werte.free;
  inherited destroy;
end;

end.