Initialer Commit

8 files changed, 1583 insertions, 0 deletions

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..92b4f6f
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,13 @@
+*.bmp
+*.png
+*.bak
+*.ppu
+*.o
+*.zip
+*.tar.gz
+*.dat
+Ergebnisse
+lib
+ROM
+*~
+Log*
diff --git a/Bild.pdf b/Bild.pdf
new file mode 100644
index 0000000..52cdcb7
--- /dev/null
+++ b/Bild.pdf
diff --git a/ROM.lpi b/ROM.lpi
new file mode 100644
index 0000000..d57d7be
--- /dev/null
+++ b/ROM.lpi
@@ -0,0 +1,96 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<CONFIG>
+  <ProjectOptions>
+    <Version Value="9"/>
+    <General>
+      <Flags>
+        <MainUnitHasCreateFormStatements Value="False"/>
+        <MainUnitHasTitleStatement Value="False"/>
+      </Flags>
+      <SessionStorage Value="InProjectDir"/>
+      <MainUnit Value="0"/>
+      <Title Value="ROM"/>
+      <UseAppBundle Value="False"/>
+      <ResourceType Value="res"/>
+    </General>
+    <i18n>
+      <EnableI18N LFM="False"/>
+    </i18n>
+    <VersionInfo>
+      <StringTable ProductVersion=""/>
+    </VersionInfo>
+    <BuildModes Count="1">
+      <Item1 Name="Default" Default="True"/>
+    </BuildModes>
+    <PublishOptions>
+      <Version Value="2"/>
+    </PublishOptions>
+    <RunParams>
+      <local>
+        <FormatVersion Value="1"/>
+      </local>
+    </RunParams>
+    <Units Count="3">
+      <Unit0>
+        <Filename Value="ROM.lpr"/>
+        <IsPartOfProject Value="True"/>
+        <UnitName Value="ROM"/>
+      </Unit0>
+      <Unit1>
+        <Filename Value="romunit.pas"/>
+        <IsPartOfProject Value="True"/>
+        <UnitName Value="romunit"/>
+      </Unit1>
+      <Unit2>
+        <Filename Value="mathunit.pas"/>
+        <IsPartOfProject Value="True"/>
+        <UnitName Value="mathunit"/>
+      </Unit2>
+    </Units>
+  </ProjectOptions>
+  <CompilerOptions>
+    <Version Value="11"/>
+    <Target>
+      <Filename Value="ROM"/>
+    </Target>
+    <SearchPaths>
+      <IncludeFiles Value="$(ProjOutDir)"/>
+      <UnitOutputDirectory Value="lib/$(TargetCPU)-$(TargetOS)"/>
+    </SearchPaths>
+    <CodeGeneration>
+      <Checks>
+        <IOChecks Value="True"/>
+        <RangeChecks Value="True"/>
+        <OverflowChecks Value="True"/>
+        <StackChecks Value="True"/>
+      </Checks>
+      <VerifyObjMethodCallValidity Value="True"/>
+    </CodeGeneration>
+    <Linking>
+      <Options>
+        <Win32>
+          <GraphicApplication Value="True"/>
+        </Win32>
+      </Options>
+    </Linking>
+    <Other>
+      <CompilerMessages>
+        <MsgFileName Value=""/>
+      </CompilerMessages>
+      <CompilerPath Value="$(CompPath)"/>
+    </Other>
+  </CompilerOptions>
+  <Debugging>
+    <Exceptions Count="3">
+      <Item1>
+        <Name Value="EAbort"/>
+      </Item1>
+      <Item2>
+        <Name Value="ECodetoolError"/>
+      </Item2>
+      <Item3>
+        <Name Value="EFOpenError"/>
+      </Item3>
+    </Exceptions>
+  </Debugging>
+</CONFIG>
diff --git a/ROM.lpr b/ROM.lpr
new file mode 100644
index 0000000..5863de0
--- /dev/null
+++ b/ROM.lpr
@@ -0,0 +1,150 @@
+program ROM;
+
+{$DEFINE UseCThreads}
+
+{$mode objfpc}{$H+}
+
+uses
+  {$IFDEF UNIX}{$IFDEF UseCThreads}
+  cthreads,
+  {$ENDIF}{$ENDIF}
+  Classes
+  { you can add units after this },
+  SysUtils,ROMunit, mathunit, Math;
+
+var inPuls,refPuls,surTraj: tExtPointArray;
+    i,smooth:               longint;
+    tmax,wmax,absShift:     extended;
+    force,fourier:          boolean;
+
+const Verwendung='Verwendung: ROM ($Einfallspuls_Datei $Ausfallspuls_Datei)/(- $trace-Datei-Prefix) $output_inPuls $output_refPuls $output_Trajektorie '+
+   '[-s/--smooth $n] [-f/--force] [-t/--tmax $t] [-w/--wmax $w] [-F/--FFT] [-d/--dt $dt]';
+
+begin
+  if paramcount<5 then Fehler(Verwendung);
+  i:=6;
+  force:=false;
+  smooth:=1;
+  tmax:=-1;
+  wmax:=-1;
+  fourier:=false;
+  absShift:=-1e9;
+  while i<=paramcount do begin
+    if (paramstr(i)='--force') or (paramstr(i)='-f') then begin
+      force:=true;
+      inc(i);
+      continue;
+    end;
+    if (paramstr(i)='--smooth') or (paramstr(i)='-s') then begin
+      inc(i);
+      smooth:=strtoint(paramstr(i));
+      inc(i);
+      continue;
+    end;
+    if (paramstr(i)='--tmax') or (paramstr(i)='-t') then begin
+      inc(i);
+      tmax:=strtofloat(paramstr(i));
+      inc(i);
+      continue;
+    end;
+    if (paramstr(i)='--wmax') or (paramstr(i)='-w') then begin
+      inc(i);
+      wmax:=strtofloat(paramstr(i));
+      fourier:=true;
+      inc(i);
+      continue;
+    end;
+    if (paramstr(i)='--dt') or (paramstr(i)='-d') then begin
+      inc(i);
+      if paramstr(i)='auto' then begin
+        absShift:=-1e9;
+        inc(i);
+        continue;
+      end;
+      if paramstr(i)='input' then begin
+        if paramstr(1)<>'-' then Fehler('Ich brauche zur Bestimmung der Gesamtverschiebung die Inputdatei vom LPIC!');
+        absShift:=-2e9;
+        inc(i);
+        continue;
+      end;
+      absShift:=strtofloat(paramstr(i));
+      inc(i);
+      continue;
+    end;
+    if (paramstr(i)='--FFT') or (paramstr(i)='-F') then begin
+      fourier:=true;
+      inc(i);
+      continue;
+    end;
+    Fehler(Verwendung);
+  end;
+  if paramstr(1)='-' then
+    readRawInputs(paramstr(2),inPuls,refPuls,absShift)
+  else begin
+    readTextInput(paramstr(1),inPuls);
+    readTextInput(paramstr(2),refPuls);
+  end;
+  write(stderr,'Input sortieren ...');
+  sort(inPuls);
+  sort(refPuls);
+  writeln(stderr,' fertig');
+  uniq(inPuls,false);
+  uniq(refPuls,false);
+  write(stderr,'Input interpolieren ...');
+  interpoliere(inPuls);
+  interpoliere(refPuls);
+  writeln(stderr,' fertig');
+  flip(refPuls);
+  integrate(inPuls);
+  integrate(refPuls);
+  removeLinearOffset(inPuls);
+  removeLinearOffset(refPuls);
+  if smooth>1 then begin
+    write(stderr,'glätten ...');
+    smoothen(inPuls,smooth);
+    smoothen(refPuls,smooth);
+    writeln(stderr,' fertig');
+  end;
+  cut(inPuls,tmax);
+  cut(refPuls,tmax);
+  for i:=3 to 5 do
+    if (paramstr(i)<>'-') and fileexists(paramstr(i)) and not force then Fehler('Die Ausgabedatei '''+paramstr(i)+''' existiert bereits!');
+  gesamtverschiebung(inPuls,refPuls,absShift);
+  write(stderr,'Trajektorie berechnen ...');
+  berechneTrajektorie(inPuls,refPuls,surTraj,absShift*byte(not fourier));
+  writeln(stderr,' fertig');
+  write(stderr,'Ergebnis sortieren ...');
+  sort(surTraj);
+  writeln(stderr,' fertig');
+  uniq(surTraj,false);
+  if fourier then begin
+    write(stderr,'Ergebnis interpolieren ...');
+    interpoliere(surTraj);
+    writeln(stderr,' fertig');
+    fft(inPuls);
+    fft(refPuls);
+    fft(surTraj);
+    inPuls[0].y:=0;
+    refPuls[0].y:=0;
+    surTraj[0].y:=0;
+    if wmax<0 then begin
+      cut(surTraj,min(refPuls[length(refPuls)-1].x,inPuls[length(inPuls)-1].x)/2);
+      cut(inPuls,surTraj[length(surTraj)-1].x);
+      cut(refPuls,surTraj[length(surTraj)-1].x);
+    end
+    else begin
+      cut(surTraj,wmax);
+      cut(inPuls,wmax);
+      cut(refPuls,wmax);
+    end;
+    write(stderr,'alles normieren ...');
+    normiere(inPuls);
+    normiere(refPuls);
+    normiere(surTraj);
+    writeln(stderr,' fertig');
+  end;
+  if paramstr(3)<>'-' then writeOutput(paramstr(3),inPuls);
+  if paramstr(4)<>'-' then writeOutput(paramstr(4),refPuls);
+  if paramstr(5)<>'-' then writeOutput(paramstr(5),surTraj);
+end.
+
diff --git a/ROM.lps b/ROM.lps
new file mode 100644
index 0000000..b6e2d2d
--- /dev/null
+++ b/ROM.lps
@@ -0,0 +1,166 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<CONFIG>
+  <ProjectSession>
+    <Version Value="9"/>
+    <BuildModes Active="Default"/>
+    <Units Count="3">
+      <Unit0>
+        <Filename Value="ROM.lpr"/>
+        <IsPartOfProject Value="True"/>
+        <UnitName Value="ROM"/>
+        <IsVisibleTab Value="True"/>
+        <EditorIndex Value="0"/>
+        <WindowIndex Value="0"/>
+        <TopLine Value="1"/>
+        <CursorPos X="151" Y="20"/>
+        <UsageCount Value="90"/>
+        <Loaded Value="True"/>
+      </Unit0>
+      <Unit1>
+        <Filename Value="romunit.pas"/>
+        <IsPartOfProject Value="True"/>
+        <UnitName Value="romunit"/>
+        <EditorIndex Value="1"/>
+        <WindowIndex Value="0"/>
+        <TopLine Value="1004"/>
+        <CursorPos X="1" Y="998"/>
+        <FoldState Value=" T3iA041 pkRkZ0Y113A"/>
+        <UsageCount Value="90"/>
+        <Loaded Value="True"/>
+      </Unit1>
+      <Unit2>
+        <Filename Value="mathunit.pas"/>
+        <IsPartOfProject Value="True"/>
+        <UnitName Value="mathunit"/>
+        <EditorIndex Value="2"/>
+        <WindowIndex Value="0"/>
+        <TopLine Value="4"/>
+        <CursorPos X="1" Y="43"/>
+        <UsageCount Value="89"/>
+        <Loaded Value="True"/>
+      </Unit2>
+    </Units>
+    <General>
+      <ActiveWindowIndexAtStart Value="0"/>
+    </General>
+    <JumpHistory Count="29" HistoryIndex="28">
+      <Position1>
+        <Filename Value="ROM.lpr"/>
+        <Caret Line="93" Column="58" TopLine="77"/>
+      </Position1>
+      <Position2>
+        <Filename Value="ROM.lpr"/>
+        <Caret Line="123" Column="66" TopLine="90"/>
+      </Position2>
+      <Position3>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="577" Column="69" TopLine="393"/>
+      </Position3>
+      <Position4>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="212" Column="2" TopLine="38"/>
+      </Position4>
+      <Position5>
+        <Filename Value="mathunit.pas"/>
+        <Caret Line="22" Column="1" TopLine="1"/>
+      </Position5>
+      <Position6>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="528" Column="14" TopLine="476"/>
+      </Position6>
+      <Position7>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="519" Column="37" TopLine="475"/>
+      </Position7>
+      <Position8>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="542" Column="1" TopLine="512"/>
+      </Position8>
+      <Position9>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="536" Column="34" TopLine="514"/>
+      </Position9>
+      <Position10>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="537" Column="25" TopLine="508"/>
+      </Position10>
+      <Position11>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="679" Column="28" TopLine="649"/>
+      </Position11>
+      <Position12>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="333" Column="150" TopLine="305"/>
+      </Position12>
+      <Position13>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="310" Column="3" TopLine="53"/>
+      </Position13>
+      <Position14>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="26" Column="95" TopLine="1"/>
+      </Position14>
+      <Position15>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="727" Column="6" TopLine="694"/>
+      </Position15>
+      <Position16>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="26" Column="1" TopLine="1"/>
+      </Position16>
+      <Position17>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="334" Column="31" TopLine="328"/>
+      </Position17>
+      <Position18>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="311" Column="41" TopLine="288"/>
+      </Position18>
+      <Position19>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="310" Column="3" TopLine="291"/>
+      </Position19>
+      <Position20>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="22" Column="1" TopLine="16"/>
+      </Position20>
+      <Position21>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="1018" Column="1" TopLine="891"/>
+      </Position21>
+      <Position22>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="669" Column="16" TopLine="625"/>
+      </Position22>
+      <Position23>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="621" Column="38" TopLine="527"/>
+      </Position23>
+      <Position24>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="312" Column="1" TopLine="295"/>
+      </Position24>
+      <Position25>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="657" Column="15" TopLine="626"/>
+      </Position25>
+      <Position26>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="655" Column="16" TopLine="637"/>
+      </Position26>
+      <Position27>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="285" Column="9" TopLine="254"/>
+      </Position27>
+      <Position28>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="312" Column="49" TopLine="301"/>
+      </Position28>
+      <Position29>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="21" Column="69" TopLine="1"/>
+      </Position29>
+    </JumpHistory>
+  </ProjectSession>
+  <EditorMacros Count="0"/>
+</CONFIG>
diff --git a/make.gnu b/make.gnu
new file mode 100644
index 0000000..9a85259
--- /dev/null
+++ b/make.gnu
@@ -0,0 +1,26 @@
+set terminal pdf
+set decimalsign ","
+
+set output "~/Dokumente/Prograemmchen/ROM/Bild.pdf"
+set title "Oberflächentrajektorie"
+set size .9,.9
+set format x "%g"
+set format y "%g"
+
+if (system("echo $logscale")=="1") {
+  set xlabel "omega in omega_L"
+  set xrange [0:40]
+  set yrange [10**-11:*]
+  set logscale y
+}
+else {
+  set xlabel "t in T"
+}
+
+#  set ylabel sprintf("$\\left|\\Big(\\mathcal{FT}\\left(%s\\right)\\Big)(\\omega)\\right|$ in $%s$",groesze,einheit)
+#    set xtics out
+#    set ytics 10**-2 out
+plot "~/Dokumente/Prograemmchen/ROM/surface.dat" using 1:2 with lines lt 1 title "surface", \
+     "~/Dokumente/Prograemmchen/ROM/input.dat" using 1:2 with lines lt 2 title "input", \
+     "~/Dokumente/Prograemmchen/ROM/output.dat" using 1:2 with lines lt 3 title "output"
+
diff --git a/mathunit.pas b/mathunit.pas
new file mode 100644
index 0000000..3d5e111
--- /dev/null
+++ b/mathunit.pas
@@ -0,0 +1,46 @@
+unit mathunit;
+
+{$mode objfpc}{$H+}
+
+interface
+
+uses
+  Classes, SysUtils;
+
+type
+  tExtPoint = record
+    x,y: extended;
+  end;
+  tExtPointArray = array of tExtPoint;
+  pTExtPointArray = ^tExtPointArray;
+  tLongintArray = array of longint;
+  pTLongintArray = ^tLongintArray;
+  tExtendedArray = array of extended;
+
+function Plus(a,b: tExtPoint): tExtPoint;
+function Durch(a: tExtPoint; b: extended): tExtPoint;
+function myFrac(x: extended): extended;
+
+implementation
+
+function Plus(a,b: tExtPoint): tExtPoint;
+begin
+  result.x:=a.x+b.x;
+  result.y:=a.y+b.y;
+end;
+
+function Durch(a: tExtPoint; b: extended): tExtPoint;
+begin
+  result.x:=a.x/b;
+  result.y:=a.y/b;
+end;
+
+function myFrac(x: extended): extended;
+begin
+  result:=frac(x);
+  while result<0 do
+    result:=result+1;
+end;
+
+end.
+
diff --git a/romunit.pas b/romunit.pas
new file mode 100644
index 0000000..aaab446
--- /dev/null
+++ b/romunit.pas
@@ -0,0 +1,1086 @@
+unit romunit;
+
+{$mode objfpc}{$H+}
+
+interface
+
+uses
+  Classes, SysUtils, Math, mathunit;
+
+procedure Fehler(s: string);
+procedure readRawInputs(nam: string; out d1,d2: tExtPointArray; var absShift: extended);
+procedure readTextInput(nam: string; out dat: tExtPointArray);
+procedure writeOutput(nam: string; const dat: tExtPointArray);
+procedure berechneTrajektorie(const inPuls,outPuls: tExtPointArray; out surTraj: tExtPointArray; absShift: extended);
+procedure smoothen(var dat: tExtPointArray; width: longint);
+procedure sort(var dat: tExtPointArray); overload;
+procedure sort(var dat: tExtPointArray; start, stop, threads: longint); overload;
+procedure cut(var dat: tExtPointArray; tmax: extended);
+function  findeExtrema(const dat: tExtPointArray; maxima: boolean): tLongintArray;
+procedure filtereExtrema(const dat: tExtPointArray; var extrema: tLongintArray; sollAbst, toleranz: extended);
+procedure monotonieHerstellen(const dat: tExtPointArray; var minima,maxima: tLongintArray);
+procedure gesamtverschiebung(var inPuls,outPuls: tExtPointArray; var absShift: extended);
+procedure integrate(var dat: tExtPointArray);
+procedure removeLinearOffset(var dat: tExtPointArray);
+procedure flip(var dat: tExtPointArray);
+procedure vereineExtrema(const dat1,dat2: tExtPointArray; var el1,el2: tLongintArray; toleranz: extended);
+procedure uniq(var dat: tExtPointArray; streng: boolean);
+procedure fft(var dat: tExtPointArray);
+procedure interpoliere(var dat: tExtPointArray);
+procedure normiere(var dat: tExtPointArray);
+
+type
+  tSortThread = class(tThread)
+    _pdat:        pTExtPointArray;
+    _start,_stop: longint;
+    fertig:       boolean;
+    constructor create(pd: pTExtPointArray; sta, sto: longint);
+    destructor destroy; override;
+    procedure execute; override;
+  end;
+  tPhaseShiftThread = class(tThread)
+    _dats:          array[boolean] of tExtPointArray;
+    _ergebnis:      tExtPointArray;
+    _starts,_stops: array[boolean] of longint;
+    _absShift:      extended;
+    fertig:         boolean;
+    constructor create(p1d,p2d: pTExtPointArray; sta1, sto1, sta2, sto2: longint; absShift: extended);
+    destructor destroy; override;
+    procedure execute; override;
+  end;
+
+implementation
+
+procedure Fehler(s: string);
+begin
+  writeln(s);
+  halt(1);
+end;
+
+procedure readRawInputs(nam: string; out d1,d2: tExtPointArray; var absShift: extended);
+var f:                             file;
+    tf:                            textfile;
+    s:                             string;
+    i,j,k,start,stop,traces,steps: longint;
+    i32:                           int32;
+    sr:                            tSearchRec;
+    fl:                            single;
+    buff:                          array of single;
+    factor,dx,cells_left,maxAmp:   extended;
+begin
+  if not fileexists(nam) then Fehler('Die Datei '''+nam+''' existiert nicht!');
+  assignfile(tf,nam);
+  i:=0;
+  reset(tf);
+  factor:=1;
+  j:=0;
+  dx:=0;
+  cells_left:=0;
+  maxAmp:=0;
+  while not eof(tf) do begin
+    readln(tf,s);
+    if pos('output path        :',s)=1 then begin
+      delete(s,1,pos(':',s));
+      while (length(s)>0) and (s[1]=' ') do
+        delete(s,1,1);
+      while (length(s)>0) and (s[length(s)]=' ') do
+        delete(s,length(s),1);
+      nam:=extractfilepath(nam)+s+'/trace-1-';
+      i:=i or 1;
+    end;
+    if s='pulse # 1' then
+      i:=i or 2;
+    if s='pulse # 2' then
+      i:=i and not $7E;
+    if s='pulse component # 0' then
+      i:=(i or $4)  and not $78;
+    if s='pulse component # 1' then
+      i:=(i or $8)  and not $74;
+    if s='pulse component # 2' then
+      i:=(i or $10) and not $7C;
+    if s='pulse component # 3' then
+      i:=(i or $20) and not $5C;
+    if odd(i shr 1) and (pos('.a0                 :',s)=1) then begin
+      delete(s,1,pos(':',s));
+      while (length(s)>0) and (s[1]=' ') do
+        delete(s,1,1);
+      if strtofloat(s)>maxAmp then begin
+        i:=i or $40;
+        maxAmp:=strtofloat(s);
+      end;
+    end;
+    if odd(i shr 6) and (pos('.frequency          :',s)=1) then begin
+      delete(s,1,pos(':',s));
+      while (length(s)>0) and (s[1]=' ') do
+        delete(s,1,1);      writeln('bla'+s+'blu');
+      factor:=strtofloat(s);
+    end;
+    case j of
+      0:
+        if s='domain: plasma' then
+          j:=1;
+      1: begin
+        if pos('cells_per_wl       :',s)=1 then begin
+          delete(s,1,pos(':',s));
+          while (length(s)>0) and (s[1]=' ') do
+            delete(s,1,1);
+          dx:=1/strtofloat(s);
+        end;
+        if pos('cells_left         :',s)=1 then begin
+          delete(s,1,pos(':',s));
+          while (length(s)>0) and (s[1]=' ') do
+            delete(s,1,1);
+          cells_left:=strtofloat(s);
+        end;
+        if s='domain: particles' then
+          j:=2;
+      end;
+    end{of case};
+  end;
+  closefile(tf);
+  if (not odd(i)) or ((absShift<-1.5e9) and ((dx=0) or (cells_left=0))) then Fehler('Unerwartetes Dateiende in '''+nam+'''!');
+
+  setlength(d1,0);
+  setlength(d2,0);
+  i32:=0;
+  fl:=0;
+  start:=-1;
+  stop:=-1;
+  i:=findFirst(nam+'*',$3F,sr);
+  while i=0 do begin
+    s:=sr.Name;
+    while pos('-',s)>0 do
+      delete(s,1,pos('-',s));
+    i:=strtoint(s);
+    if (start=-1) or (start>i) then
+      start:=i;
+    if (stop=-1) or (stop<i) then
+      stop:=i;
+    i:=findNext(sr);
+  end;
+  findClose(sr);
+  if start=-1 then
+    Fehler('Keine Dateien '''+nam+'*'' gefunden!');
+  for i:=start to stop do
+    if not fileexists(nam+inttostr(i)) then
+      Fehler('Die Dateien beginnend mit '''+nam+''' bilden keine vollständige Abfolge! ('+inttostr(i)+' fehlt)');
+
+  for i:=start to stop do begin
+    write(stderr,'Datei '''+nam+inttostr(i)+''' einlesen ');
+    assignfile(f,nam+inttostr(i));
+    reset(f,1);
+    blockread(f,i32,sizeof(int32));              // Dateinummer
+    if i32<>i then begin
+      closefile(f);
+      Fehler('Die Datei '''+nam+inttostr(i)+''' behauptet, sie sei die '+inttostr(i32)+'-te!');
+    end;
+
+    blockread(f,i32,sizeof(int32));              // Anzahl x-Positionen
+    traces:=i32;
+    blockread(f,i32,sizeof(int32));              // Anzahl t-Positionen
+    steps:=i32;
+
+    setlength(buff,steps);
+    setlength(d1,length(d1)+steps);
+    setlength(d2,length(d2)+steps);
+
+    for j:=0 to traces-1 do begin
+      blockread(f,fl,sizeof(single));            // x-Position
+      if j=0 then begin
+        if i=start then cells_left:=cells_left-fl;
+        blockread(f,buff[0],sizeof(single)*steps); // fp
+        for k:=0 to steps-1 do begin
+          d1[k+length(d1)-steps].x:=(i-1+k/steps)*factor;
+          d1[k+length(d1)-steps].y:=buff[k];
+        end;
+        blockread(f,buff[0],sizeof(single)*steps); // fm
+        for k:=0 to steps-1 do begin
+          d2[k+length(d2)-steps].x:=(i-1+k/steps)*factor;
+          d2[k+length(d2)-steps].y:=buff[k];
+        end;
+      end;
+      for k:=0 to 7+2*byte(j<>0) do                // (fp,fm,)gp,gm,ex,de,di,jx,jy,jz
+        blockread(f,buff[0],sizeof(single)*steps);
+    end;
+
+    if not eof(f) then
+      Fehler('Die Datei '''+nam+inttostr(i)+''' ist zu groß!');
+
+    closefile(f);
+    writeln(stderr,'fertig');
+  end;
+  if absShift<-1.5e9 then begin
+    absShift:=dx*cells_left*factor*2;
+    writeln(stderr,'Aus Inputdatei ausgelesene konstante Verschiebung des Outputpulses: '+floattostr(absShift)+' T');
+  end;
+end;
+
+procedure readTextInput(nam: string; out dat: tExtPointArray);
+var f:  textfile;
+    i:  longint;
+    s:  string;
+begin
+  if not fileexists(nam) then
+    Fehler('Datei '''+nam+''' existiert nicht!');
+  write(stderr,'Datei '''+nam+''' einlesen ');
+  assignfile(f,nam);
+  reset(f);
+  i:=0;
+  while not eof(f) do begin
+    readln(f,s);
+    if pos('#',s)>0 then delete(s,pos('#',s),length(s));
+    while (length(s)>0) and (s[1]=' ') do delete(s,1,1);
+    if s='' then continue;
+    inc(i);
+  end;
+  closefile(f);
+  reset(f);
+  setlength(dat,i);
+  for i:=0 to length(dat)-1 do begin
+    readln(f,s);
+    if pos('#',s)>0 then delete(s,pos('#',s),length(s));
+    while (length(s)>0) and (s[1]=' ') do delete(s,1,1);
+    if s='' then continue;
+    dat[i].x:=strtofloat(copy(s,1,pos(' ',s)-1));
+    delete(s,1,pos(' ',s));
+    while (length(s)>0) and (s[1]=' ') do delete(s,1,1);
+    dat[i].y:=strtofloat(copy(s,1,pos(' ',s+' ')-1));
+    if i and 1023 = 0 then write(stderr,'.');
+  end;
+  closefile(f);
+  writeln(stderr,'fertig');
+end;
+
+procedure writeOutput(nam: string; const dat: tExtPointArray);
+var f:       textfile;
+    i,ml,mr: longint;
+    leer:    string;
+    strings: array of array[boolean] of string;
+begin
+  if length(dat)=0 then exit;
+  write(stderr,'Datei '''+nam+''' schreiben ');
+  assignfile(f,nam);
+  rewrite(f);
+  setlength(strings,length(dat));
+  ml:=0;
+  mr:=0;
+  for i:=0 to length(dat)-1 do begin
+    strings[i,false]:=floattostr(dat[i].x);
+    strings[i,true]:=floattostr(dat[i].y);
+    ml:=max(ml,length(strings[i,false]));
+    mr:=max(mr,pos('.',strings[i,true]+'.'));
+    if i and 65535 = 0 then write(stderr,'.');
+  end;
+  setlength(leer,ml+mr);
+  fillchar(leer[1],length(leer),' ');
+  for i:=0 to length(dat)-1 do begin
+    writeln(f,strings[i,false]+copy(leer,1,1+ml+mr-length(strings[i,false])-pos('.',strings[i,true]+'.'))+strings[i,true]);
+    if i and 65535 = 0 then write(stderr,'.');
+  end;
+  closefile(f);
+  writeln(stderr,' fertig');
+end;
+
+procedure berechneTrajektorie(const inPuls,outPuls: tExtPointArray; out surTraj: tExtPointArray; absShift: extended);
+var i,j:       longint;
+    exList:    array[0..3] of tLongintArray;
+    psThreads: array of tPhaseShiftThread;
+    stuecke:   array of array[boolean,boolean] of longint;
+    erledigt:  array of boolean;
+    b1,b2:     boolean;
+begin
+  for i:=0 to length(inPuls)-2 do
+    if inPuls[i+1].x<=inPuls[i].x then
+      writeln('Der Inpuls ist nicht sortiert!');
+  for i:=0 to length(outPuls)-2 do
+    if outPuls[i+1].x<=outPuls[i].x then
+      writeln('Der Outpuls ist nicht sortiert!');
+
+  // Extrema identifizieren
+  exList[0]:=findeExtrema(inPuls,false);
+  exList[1]:=findeExtrema(inPuls,true);
+  exList[2]:=findeExtrema(outPuls,false);
+  exList[3]:=findeExtrema(outPuls,true);
+
+  // und filtern
+  filtereExtrema(inPuls,exList[0],1,2);
+  filtereExtrema(inPuls,exList[1],1,2);
+  filtereExtrema(outPuls,exList[2],1,2);
+  filtereExtrema(outPuls,exList[3],1,2);
+
+  monotonieHerstellen(inPuls,exList[0],exList[1]);
+  monotonieHerstellen(outPuls,exList[2],exList[3]);
+
+  writeln;
+  for i:=0 to 3 do begin
+    for j:=0 to length(exList[i])-1 do
+      write(exList[i,j],'(',j,') ');
+    writeln;
+  end;
+
+  vereineExtrema(inPuls,outPuls,exList[0],exList[2],0.5);
+  vereineExtrema(inPuls,outPuls,exList[1],exList[3],0.5);
+
+writeln(length(exList[0]),' ',length(exList[1]),' ',length(exList[2]),' ',length(exList[3]));
+
+{  for k:=0 to 1 do begin
+    i:=0;
+    while (i<length(exList[k])-1) and (exList[k,i+1]<exList[1-k,0]) do
+      inc(i);
+    for j:=i to length(exList[k])-1 do
+      exList[k,j-i]:=exList[k,j];
+    setlength(exList[k],length(exList[k])-i);
+    for j:=i to length(exList[2+k])-1 do
+      exList[2+k,j-i]:=exList[2+k,j];
+    setlength(exList[2+k],length(exList[2+k])-i);
+
+    i:=length(exList[k])-1;
+    while (i>0) and (exList[k,i-1]>exList[1-k,length(exList[1-k])-1]) do
+      dec(i);
+    setlength(exList[k],i+1);
+    setlength(exList[2+k],i+1);
+  end; }
+
+  writeln;
+  for i:=0 to 3 do begin
+    for j:=0 to length(exList[i])-1 do
+      write(exList[i,j],'(',j,') ');
+    writeln;
+  end;
+
+  setlength(stuecke,length(exList[0])+length(exList[1])-1);
+    writeln(length(stuecke),' ',length(exList[0]),' ',length(exList[1]),' ',length(exList[2]),' ',length(exList[3]));
+  for i:=0 to length(stuecke)-1 do
+    for b1:=false to true do    // Anfang / Ende
+      for b2:=false to true do  // input / output
+        stuecke[i,b1,b2]:=
+          exList[(i and 1) xor byte(exList[0,0]>exList[1,0]) xor (2*byte(b2)) xor byte(b1),(i+byte(b1)) div 2];
+
+  writeln;
+  for i:=0 to length(stuecke)-1 do
+    writeln(i,': ',stuecke[i,false,false],' ',stuecke[i,false,true],' ',stuecke[i,true,false],' ',stuecke[i,true,true]);
+  writeln;
+{  for i:=0 to length(stuecke)-1 do
+    writeln(inPuls[stuecke[i,false,false]].x,' ',outPuls[stuecke[i,false,true]].x,' ',inPuls[stuecke[i,true,false]].x,' ',outPuls[stuecke[i,true,true]].x);
+  writeln; }
+
+  setlength(erledigt,length(stuecke));
+  for i:=0 to length(erledigt)-1 do
+    erledigt[i]:=false;
+  setlength(surTraj,0);
+  setlength(psThreads,8);
+  for i:=0 to length(psThreads)-1 do
+    psThreads[i]:=nil;
+  b1:=true;
+  repeat
+    b1:=true;
+    for i:=0 to length(psThreads)-1 do begin
+      if (psThreads[i]<>nil) and (psThreads[i].fertig) then begin
+        setlength(surTraj,length(surTraj)+length(psThreads[i]._ergebnis));
+        for j:=0 to length(psThreads[i]._ergebnis)-1 do
+          surTraj[length(surTraj)-1-j]:=psThreads[i]._ergebnis[j];
+        psThreads[i].free;
+        psThreads[i]:=nil;
+      end;
+      if psThreads[i]=nil then begin
+        j:=0;
+        while (j<length(erledigt)) and erledigt[j] do
+          inc(j);
+        if j>=length(erledigt) then continue;
+        psThreads[i]:=
+          tPhaseShiftThread.create(
+            @inPuls,
+            @outPuls,
+            stuecke[j,false,false],
+            stuecke[j,true,false],
+            stuecke[j,false,true],
+            stuecke[j,true,true],
+            absShift);
+        psThreads[i].Suspended:=false;
+        erledigt[j]:=true;
+      end;
+      b1:=false;
+    end;
+    if not b1 then sleep(100);
+    for i:=0 to length(erledigt)-1 do
+      b1:=b1 and erledigt[i];
+  until b1;
+
+
+(*  setlength(surTraj,min(length(inPuls),length(outPuls)));
+  for i:=0 to length(surTraj)-1 do begin
+    surTraj[i].x:=(inPuls[i].x+outPuls[i].x)/2;
+    surTraj[i].y:=inPuls[i].y+outPuls[i].y;
+  end;
+  for i:=0 to length(exList[0])-1 do
+    surTraj[exList[0,i]].y:=-1;
+  for i:=0 to length(exList[1])-1 do
+    surTraj[exList[1,i]].y:=1;
+  for i:=0 to length(exList[2])-1 do
+    surTraj[exList[2,i]].y:=0;
+  for i:=0 to length(exList[3])-1 do
+    surTraj[exList[3,i]].y:=0; *)
+end;
+
+procedure smoothen(var dat: tExtPointArray; width: longint);
+var i,j: longint;
+begin
+  for i:=0 to length(dat)-width do begin
+    for j:=1 to width-1 do
+      dat[i]:=Plus(dat[i],dat[i+j]);
+    dat[i]:=Durch(dat[i],width);
+  end;
+  setlength(dat,length(dat)-width+1);
+end;
+
+procedure sort(var dat: tExtPointArray);
+begin
+  sort(dat,0,length(dat)-1,8);
+end;
+
+procedure sort(var dat: tExtPointArray; start, stop, threads: longint);
+var sortThreads: array of tSortThread;
+    i,j,k:       longint;
+    fertig:      boolean;
+    apos,epos:   array of longint;
+    tmp:         tExtPointArray;
+begin
+  if start>=stop then exit;
+  if threads=1 then begin
+    sort(dat,start,(start+stop) div 2,threads);
+    sort(dat,(start+stop) div 2 + 1,stop,threads);
+    setlength(apos,2);
+    setlength(epos,2);
+    apos[0]:=start;
+    apos[1]:=(start+stop) div 2 + 1;
+    epos[0]:=(start+stop) div 2;
+    epos[1]:=stop;
+  end
+  else begin
+    threads:=min(threads,stop-start+1);
+    setlength(sortThreads,threads);
+    for i:=0 to length(sortThreads)-1 do
+      sortThreads[i]:=
+        tSortThread.create(
+          @dat,
+          start+round(i  *  (stop+1-start)/threads),
+          start+round((i+1)*(stop+1-start)/threads)-1);
+    for i:=0 to length(sortThreads)-1 do
+      sortThreads[i].Suspended:=false;
+    repeat
+      fertig:=true;
+      sleep(100);
+      for i:=0 to length(sortThreads)-1 do
+        fertig:=fertig and sortThreads[i].fertig;
+    until fertig;
+    for i:=0 to length(sortThreads)-1 do
+      sortThreads[i].free;
+    setlength(sortThreads,0);
+    setlength(apos,threads);
+    setlength(epos,threads);
+    for i:=0 to threads-1 do begin
+      apos[i]:=start+round(i  *  (stop+1-start)/threads);
+      epos[i]:=start+round((i+1)*(stop+1-start)/threads)-1;
+    end;
+  end;
+  // und jetzt nur noch mergen
+  setlength(tmp,stop-start+1);
+  k:=0;
+  repeat
+    fertig:=true;
+    j:=-1;
+    for i:=0 to length(apos)-1 do begin
+      if apos[i]>epos[i] then continue;
+      fertig:=false;
+      if (j=-1) or (dat[apos[i]].x<dat[apos[j]].x) then j:=i;
+    end;
+    if not fertig then begin
+      tmp[k]:=dat[apos[j]];
+      inc(k);
+      inc(apos[j]);
+    end;
+  until fertig;
+  for i:=start to stop do
+    dat[i]:=tmp[i-start];
+  setlength(tmp,0);
+end;
+
+procedure cut(var dat: tExtPointArray; tmax: extended);
+var i: longint;
+begin
+  if tmax<0 then exit;
+  i:=0;
+  while (i<length(dat)) and (dat[i].x<tmax) do
+    inc(i);
+  setlength(dat,i);
+end;
+
+function  findeExtrema(const dat: tExtPointArray; maxima: boolean): tLongintArray;
+var i,j: longint;
+begin
+  j:=0;
+  for i:=1 to length(dat)-2 do
+    if (maxima    and    (dat[i].y>dat[i-1].y) and (dat[i].y>=dat[i+1].y) and (dat[i].y>0)) or
+       ((not maxima) and (dat[i].y<dat[i-1].y) and (dat[i].y<=dat[i+1].y) and (dat[i].y<0)) then
+      inc(j);
+  setlength(result,j);
+  j:=0;
+  for i:=1 to length(dat)-2 do
+    if (maxima    and    (dat[i].y>dat[i-1].y) and (dat[i].y>=dat[i+1].y) and (dat[i].y>0)) or
+       ((not maxima) and (dat[i].y<dat[i-1].y) and (dat[i].y<=dat[i+1].y) and (dat[i].y<0)) then begin
+      result[j]:=i;
+      inc(j);
+    end;
+end;
+
+procedure filtereExtrema(const dat: tExtPointArray; var extrema: tLongintArray; sollAbst, toleranz: extended);
+var mitte,i,j,k: longint;
+    mx,dist:     extended;
+    behalten:    array of boolean;
+    Toepfe:      array of longint;
+begin
+  setlength(behalten,length(extrema));
+  for i:=0 to length(behalten)-1 do
+    behalten[i]:=false;
+
+  setlength(toepfe,round(sqrt(length(extrema))));
+  for i:=0 to length(toepfe)-1 do
+    toepfe[i]:=0;
+
+  for i:=0 to length(extrema)-1 do
+    inc(toepfe[round(dat[extrema[i]].x/sollAbst*length(toepfe)) mod length(toepfe)]);
+
+  mitte:=0;
+  for i:=1 to length(toepfe)-1 do
+    if toepfe[i]>toepfe[mitte] then
+      mitte:=i;
+
+  mx:=0;
+  for i:=0 to length(extrema)-1 do
+    mx:=mx + myfrac(dat[extrema[i]].x - mitte/length(toepfe) + 0.5);
+
+  mx:=(mx/length(extrema) + mitte/length(toepfe) - 0.5)*sollAbst;
+
+  for i:=round((dat[extrema[0]].x - mx)/sollAbst) to round((dat[extrema[length(extrema)-1]].x - mx)/sollAbst) do begin
+    dist:=abs(dat[extrema[0]].x - i*sollAbst - mx);
+    k:=0;
+    for j:=1 to length(behalten)-1 do
+      if abs(dat[extrema[j]].x - i*sollAbst - mx) < dist then begin
+        dist:=abs(dat[extrema[j]].x - i*sollAbst - mx);
+        k:=j;
+      end;
+    if 2*dist<toleranz then
+      behalten[k]:=true;
+  end;
+
+(*// Variante 1
+  mx:=(dat[0].x+dat[length(dat)-1].x)/2;
+  mitte:=0;
+  for i:=1 to length(extrema)-1 do
+    if abs(dat[extrema[mitte]].x-mx)>abs(dat[extrema[i]].x-mx) then
+      mitte:=i;
+  repeat
+    j:=mitte;
+    repeat
+      k:=-1;
+      for i:=j-1 downto 0 do
+        if (abs(dat[extrema[j]].x-dat[extrema[i]].x-sollAbst)<=toleranz) and
+           ((k=-1) or (abs(dat[extrema[j]].x-dat[extrema[i]].x-sollAbst)<abs(dat[extrema[j]].x-dat[extrema[k]].x-sollAbst))) then
+          k:=i;
+      if k<>-1 then begin
+        behalten[k]:=true;
+        behalten[j]:=true;
+        j:=k;
+      end;
+    until k=-1;
+    if j=mitte then
+      inc(mitte);
+  until (j<>mitte) or (mitte>length(extrema));
+  if j=mitte then Fehler('So, wie sich der Programmierer das gedacht hat, findet man die richtigen Maxima nicht ...');
+  j:=mitte;
+  repeat
+    k:=-1;
+    for i:=j+1 to length(extrema)-1 do
+      if (abs(dat[extrema[i]].x-dat[extrema[j]].x-sollAbst)<=toleranz) and
+         ((k=-1) or (abs(dat[extrema[i]].x-dat[extrema[j]].x-sollAbst)<abs(dat[extrema[k]].x-dat[extrema[j]].x-sollAbst))) then
+        k:=i;
+    if k<>-1 then begin
+      behalten[k]:=true;
+      behalten[j]:=true;
+      j:=k;
+    end;
+  until k=-1;
+  if j=mitte then Fehler('So, wie sich der Programmierer das gedacht hat, findet man die richtigen Maxima nicht ...');
+*)
+
+  j:=0;
+  for i:=0 to length(behalten)-1 do
+    if behalten[i] then begin
+      extrema[j]:=extrema[i];
+//      write(extrema[j],' ');
+      inc(j);
+    end;
+  setlength(extrema,j);
+//  writeln;
+end;
+
+procedure monotonieHerstellen(const dat: tExtPointArray; var minima,maxima: tLongintArray);
+var i,j,k:    longint;
+    b:        boolean;
+    behalten: array[boolean] of array of Boolean;
+begin
+  setlength(behalten[false],length(maxima));
+  setlength(behalten[true],length(minima));
+  for b:=false to true do
+    for i:=0 to length(behalten[b])-1 do
+      behalten[b,i]:=false;
+
+  i:=0;
+  j:=0;
+  b:=maxima[0]>minima[0];
+  repeat
+    if b then begin
+      k:=j;
+      repeat
+        if (j<length(minima)) and (dat[minima[j]].y<dat[minima[k]].y) then
+          k:=j;
+        inc(j);
+      until (j>=length(minima)) or ((i<length(maxima)) and (maxima[i]<minima[j]));
+      behalten[true,k]:=true;
+      b:=false;
+    end
+    else begin
+      k:=i;
+      repeat
+        if (i<length(maxima)) and (dat[maxima[i]].y>dat[maxima[k]].y) then
+          k:=i;
+        inc(i);
+      until (i>=length(maxima)) or ((j<length(minima)) and (maxima[i]>=minima[j]));
+      behalten[false,k]:=true;
+      b:=true;
+    end;
+  until (i>=length(maxima)) and (j>=length(minima));
+
+  i:=0;
+  j:=0;
+  while i<length(maxima) do begin
+    if behalten[false,i] then begin
+      maxima[j]:=maxima[i];
+      inc(j);
+    end;
+    inc(i);
+  end;
+  setlength(maxima,j);
+
+  i:=0;
+  j:=0;
+  while i<length(minima) do begin
+    if behalten[true,i] then begin
+      minima[j]:=minima[i];
+      inc(j);
+    end;
+    inc(i);
+  end;
+  setlength(minima,j);  writeln(length(maxima),' ',length(minima));
+end;
+
+procedure gesamtverschiebung(var inPuls,outPuls: tExtPointArray; var absShift: extended);
+var iMax,oMax,i: longint;
+begin
+  iMax:=0;
+  for i:=1 to length(inPuls)-1 do
+    if inPuls[i].y>inPuls[iMax].y then iMax:=i;
+  oMax:=0;
+  for i:=1 to length(outPuls)-1 do
+    if outPuls[i].y>outPuls[oMax].y then oMax:=i;
+  if absShift<-0.9e9 then absShift:=outPuls[oMax].x-inPuls[iMax].x
+  else absShift:=(outPuls[oMax].x-inPuls[iMax].x)-round((outPuls[oMax].x-inPuls[iMax].x)-absShift);
+
+  for i:=0 to length(outPuls)-1 do
+    outPuls[i].x:=outPuls[i].x-absShift;
+  oMax:=0;
+  while (oMax<length(outPuls)) and (outPuls[oMax].x<inPuls[0].x) do
+    inc(oMax);
+  for i:=oMax to length(outPuls)-1 do
+    outPuls[i-oMax]:=outPuls[i];
+  setlength(outPuls,length(outPuls)-oMax);
+  cut(inPuls,outPuls[length(outPuls)-1].x);
+  writeln(stderr,'Die konstante Verschiebung wurde auf '+floattostr(absShift)+' T optimiert.');
+end;
+
+procedure integrate(var dat: tExtPointArray);
+var i: longint;
+begin
+  for i:=1 to length(dat)-1 do
+    dat[i].y:=dat[i].y * (dat[i].x-dat[i-1].x)+dat[i-1].y;
+end;
+
+procedure removeLinearOffset(var dat: tExtPointArray);
+var i:     longint;
+    dx,dy: extended;
+begin
+  dx:=0;
+  dy:=0;
+  for i:=0 to 99 do begin
+    dy:=dy + dat[length(dat)-1-i].y - dat[i].y;
+    dx:=dx + dat[length(dat)-1-i].x - dat[i].x;
+  end;
+  for i:=0 to length(dat)-1 do
+    dat[i].y:=dat[i].y - dy/dx*(dat[i].x-dat[0].x);
+  dy:=0;
+  for i:=0 to 99 do
+    dy:=dy + dat[length(dat)-1-i].y + dat[i].y;
+  dy:=dy/200;
+  for i:=0 to length(dat)-1 do
+    dat[i].y:=dat[i].y - dy;
+end;
+
+procedure flip(var dat: tExtPointArray);
+var i: longint;
+begin
+  for i:=0 to length(dat)-1 do
+    dat[i].y:=-dat[i].y;
+end;
+
+procedure vereineExtrema(const dat1,dat2: tExtPointArray; var el1,el2: tLongintArray; toleranz: extended);
+var i,j,k:    longint;
+    b:        boolean;
+    behalten: array[boolean] of array of boolean;
+begin
+  setlength(behalten[false],length(el1));
+  setlength(behalten[true],length(el2));
+  for b:=false to true do
+    for i:=0 to length(behalten[b])-1 do
+      behalten[b,i]:=false;
+
+(*  i:=0;
+  j:=0;
+  b:=dat1[el1[0]].x > dat2[el2[0]].x;
+  repeat
+    if b then begin
+      behalten[true,i]:=true;
+      while (i+1<length(el1)) and ((j>=length(el2)) or (dat1[el1[i+1]].x <= dat2[el2[j]].x)) do
+        inc(i);
+      inc(i);
+      b:=false;
+    end
+    else begin
+      behalten[false,j]:=true;
+      while (j+1<length(el2)) and ((i>=length(el1)) or (dat1[el1[i]].x > dat2[el2[j+1]].x)) do
+        inc(j);
+      inc(j);
+      b:=true;
+    end;
+  until (i>=length(el1)) and (j>=length(el2)); *)
+
+  for i:=0 to length(el1)-1 do begin
+    k:=-1;
+    for j:=0 to length(el2)-1 do
+      if (abs(dat1[el1[i]].x-dat2[el2[j]].x)<=toleranz) and
+         ((k=-1) or (abs(dat1[el1[i]].x-dat2[el2[j]].x)<abs(dat1[el1[i]].x-dat2[el2[k]].x))) then
+        k:=j;
+    if (k<>-1) and not behalten[true,k] then begin
+      behalten[false,i]:=true;
+      behalten[true,k]:=true;
+    end;
+  end;
+
+  j:=0;
+  for i:=0 to length(behalten[false])-1 do
+    if behalten[false,i] then begin
+      el1[j]:=el1[i];
+      inc(j);
+    end;
+  setlength(el1,j);
+  j:=0;
+  for i:=0 to length(behalten[true])-1 do
+    if behalten[true,i] then begin
+      el2[j]:=el2[i];
+      inc(j);
+    end;
+  setlength(el2,j);
+end;
+
+procedure uniq(var dat: tExtPointArray; streng: boolean);
+var i,j,k: longint;
+    tmp:   extended;
+    s:     string;
+begin
+  write(stderr,'Dopplungen entfernen ...');
+  j:=0;
+  k:=0;
+  tmp:=0;
+  for i:=0 to length(dat)-1 do begin
+    if (i<length(dat)-1) and (dat[i].x=dat[i+1].x) and ((dat[i].y=dat[i+1].y) or not streng) then begin
+      inc(k);
+      tmp:=tmp+dat[i].y;
+      continue;
+    end;
+    if (k>0) then begin
+      dat[j].y:=(tmp+dat[i].y)/(k+1);
+      dat[j].x:=dat[i].x;
+      k:=0;
+      tmp:=0;
+    end
+    else dat[j]:=dat[i];
+    inc(j);
+  end;
+  s:=floattostr(round(100*(length(dat)-j)/length(dat)*10)/10);
+  if pos('.',s)>0 then delete(s,pos('.',s)+2,length(s));
+  writeln(stderr,' fertig (es gab '+inttostr(length(dat)-j)+' Dopplungen - etwa '+s+'%)');
+  setlength(dat,j);
+end;
+
+procedure fft(var dat: tExtPointArray);
+var i,j,k,n,dist,absch,wnum,wstep,haL:   longint;
+    in0,out0:                            boolean;
+    wRe,wIm:                             tExtendedArray;
+    t1,t2,vorher,nachher,fstep,pvFehler: extended;
+    umsortierung:                        tLongintArray;
+begin
+  write(stderr,'FFT ... ');
+  fstep:=dat[1].x-dat[0].x;
+  j:=length(dat);
+  setlength(dat,2*round(power(2,ceil(ln(length(dat))/ln(2)))));
+  k:=(length(dat)-j) div 2;
+  for i:=j-1 downto 0 do
+    dat[i+k]:=dat[i];
+  for i:=0 to k-1 do                     // weich auf 0 abklingen lassen
+    dat[i].y:=dat[k].y   *   sqr(sin(pi/2*i/k));
+  for i:=j+k to length(dat)-1 do         // dito
+    dat[i].y:=dat[j+k-1].y * sqr(sin(pi/2*(length(dat)-1-i)/(length(dat)-j-k)));
+  for i:=0 to length(dat)-1 do begin
+    dat[i].x:=dat[i].y;
+    dat[i].y:=0;
+  end;
+
+  haL:=length(dat) div 2;
+  n:=round(ln(2*hal)/ln(2));
+
+  fstep:=1/(2*haL)/fstep;
+
+  vorher:=0;
+  for i:=0 to 2*haL-1 do
+    vorher:=vorher + sqr(dat[i].x);
+
+  setlength(umsortierung,2*haL);
+  for j:=0 to 2*haL-1 do begin
+    k:=0;
+    for i:=0 to n-1 do             // Bits umkehren
+      k:=2*k + byte(odd(j shr i));
+    umsortierung[j]:=k;
+  end;
+  setlength(wRe,haL);
+  setlength(wIm,haL);
+  for i:=0 to haL-1 do begin
+    wRe[i]:=cos(pi*i/haL);
+    wIm[i]:=sin(pi*i/haL);
+  end;
+
+  in0:=true;
+  out0:=true;
+  for i:=0 to 2*haL-1 do begin
+    if umsortierung[i]>i then begin
+      t1:=dat[i].x;
+      dat[i].x:=dat[umsortierung[i]].x;
+      dat[umsortierung[i]].x:=t1;
+    end;
+    in0:=in0 and (dat[i].x=0);
+  end;
+
+  dist:=1;
+  wstep:=haL;
+  while dist<2*haL do begin
+    absch:=0;
+    while absch<2*haL do begin
+      wnum:=0;
+      for i:=0 to dist-1 do begin
+        // x_links:  [absch+j]
+        // x_rechts: [absch+j+dist]
+        t1:=wRe[wnum]*dat[absch+i+dist].x - wIm[wnum]*dat[absch+i+dist].y;
+        t2:=wRe[wnum]*dat[absch+i+dist].y + wIm[wnum]*dat[absch+i+dist].x;
+
+        dat[absch+i+dist].x:=dat[absch+i].x-t1;
+        dat[absch+i+dist].y:=dat[absch+i].y-t2;
+        dat[absch+i].x:=dat[absch+i].x+t1;
+        dat[absch+i].y:=dat[absch+i].y+t2;
+        wnum:=wnum+wstep;
+      end;
+      absch:=absch+2*dist;
+    end;
+    wstep:=wstep div 2;
+    dist:=dist*2;
+  end;
+
+  for i:=0 to 2*haL-1 do begin
+    dat[i].y:=(sqr(dat[i].x)+sqr(dat[i].y))/(2*haL);
+    dat[i].x:=fStep*i;
+  end;
+  for j:=0 to 2*haL-1 do
+    out0:=out0 and (dat[i].y=0);
+
+  nachher:=0;
+  for i:=0 to 2*haL-1 do
+    nachher:=nachher + dat[i].y;
+
+  if (nachher=0) and (vorher=0) then pvFehler:=0
+  else pvFehler:=abs(nachher-vorher)/(nachher+vorher);
+
+  writeln(stderr,' fertig (Parseval-Fehler = '+floattostr(pvFehler)+' ... nämlich von '+floattostr(vorher)+' zu '+floattostr(nachher)+')');
+  if in0 then writeln(stderr,'Nur Nullen im Input der FFT!');
+  if out0 then writeln(stderr,'Nur Nullen im Output der FFT!');
+end;
+
+procedure interpoliere(var dat: tExtPointArray);
+var i,j:  longint;
+    tmp:  extended;
+    tdat: tExtPointArray;
+begin
+  tmp:=dat[1].x-dat[0].x;
+  for i:=2 to length(dat)-1 do
+    tmp:=min(tmp,dat[i].x-dat[i-1].x);
+  if tmp<=0 then Fehler('Die Daten müssen sortiert sein und dürfen keine doppelten x-Werte enthalten! ('+floattostr(tmp)+')');
+  setlength(tdat,length(dat));
+  for i:=0 to length(dat)-1 do
+    tdat[i]:=dat[i];
+  setlength(dat,max(2*length(tdat),round(min(power(length(tdat),1.3),(tdat[length(tdat)-1].x-tdat[0].x)/tmp+1))));
+  j:=0;
+  for i:=0 to length(dat)-1 do begin
+    dat[i].x:=tdat[0].x+(tdat[length(tdat)-1].x-tdat[0].x)*i/(length(dat)-1);
+    while (j<length(tdat)-1) and (tdat[j+1].x <= dat[i].x) do
+      inc(j);
+    if j=length(tdat)-1 then tmp:=tdat[j].y
+    else begin
+      tmp:=(dat[i].x-tdat[j].x)/(tdat[j+1].x-tdat[j].x);
+      tmp:=tdat[j].y*(1-tmp) + tdat[j+1].y*tmp;
+    end;
+    dat[i].y:=tmp;
+  end;
+end;
+
+procedure normiere(var dat: tExtPointArray);
+var i: longint;
+    m: extended;
+begin
+  if length(dat)=0 then exit;
+  m:=dat[0].y;
+  for i:=1 to length(dat)-1 do
+    m:=max(m,dat[i].y);
+  if m=0 then exit;
+  for i:=0 to length(dat)-1 do
+    dat[i].y:=dat[i].y/m;
+end;
+
+// tSortThread *****************************************************************
+
+constructor tSortThread.create(pd: pTExtPointArray; sta, sto: longint);
+begin
+  inherited create(true);
+  _pdat:=pd;
+  _start:=sta;
+  _stop:=sto;
+  fertig:=false;
+end;
+
+destructor tSortThread.destroy;
+begin
+  inherited destroy;
+end;
+
+procedure tSortThread.execute;
+begin
+  sort(_pdat^,_start,_stop,1);
+  fertig:=true;
+end;
+
+// tPhaseShiftThread ***********************************************************
+
+constructor tPhaseShiftThread.create(p1d,p2d: pTExtPointArray; sta1, sto1, sta2, sto2: longint; absShift: extended);
+var i: longint;
+    b: boolean;
+begin
+  inherited create(true);
+  _starts[false]:=sta1;
+  _stops[false]:=sto1;
+  _starts[true]:=sta2;
+  _stops[true]:=sto2;
+  _absShift:=absShift;
+  for b:=false to true do
+    setlength(_dats[b],_stops[b]-_starts[b]+1);
+  for i:=0 to length(_dats[false])-1 do
+    _dats[false,i]:=p1d^[i+_starts[false]];
+  for i:=0 to length(_dats[true])-1 do
+    _dats[true,i]:=p2d^[i+_starts[true]];
+  setlength(_ergebnis,2+sto1+sto2-sta1-sta2);
+  fertig:=false;
+end;
+
+destructor tPhaseShiftThread.destroy;
+begin
+  setlength(_ergebnis,0);
+  inherited destroy;
+end;
+
+procedure tPhaseShiftThread.execute;
+var i,j,offset:               longint;
+    b:                        boolean;
+    t1,t2,t3,bestErr,bestPos: extended;
+begin
+  for b:=false to true do begin // normieren (auf 0..1)
+    t1:=0;
+    t2:=0;
+    for i:=0 to length(_dats[b])-1 do begin
+      t1:=max(t1,_dats[b,i].y);
+      t2:=min(t2,_dats[b,i].y);
+    end;
+    for i:=0 to length(_dats[b])-1 do
+      _dats[b,i].y:=(_dats[b,i].y-t2)/(t1-t2);
+  end;
+
+  for b:=false to true do begin
+    offset:=byte(b)*(_stops[false]-_starts[false]+1);
+    for i:=0 to _stops[b]-_starts[b] do begin
+      _ergebnis[offset + i].x:=_dats[b,i].x;
+      if (i=0) then bestPos:=_dats[not b,0].x
+      else begin
+        t1:=_dats[b,i].y;
+        j:=0; // floor(_ergebnis[offset + i - 1].y); <- das war Murks
+(*        while (j<_stops[not b]-_starts[not b]) and
+              (_dats[not b,j].x<_ergebnis[offset + i - 1].y) do
+          inc(j);*)
+        bestPos:=-1;
+        bestErr:=-1;
+        repeat
+          while (j<_stops[not b]-_starts[not b]) and
+                ((_dats[not b,j].y - t1) * (_dats[not b,j + 1].y - t1) > 0) do
+            inc(j);
+          t2:=_dats[not b,j].x;
+          if j<_stops[not b]-_starts[not b] then // x1 = dx * y1 / (y1 - y2) ... Berechnung des gebrochenen Anteils
+            t2:=t2 + (_dats[not b,j+1].x-_dats[not b,j].x) * (_dats[not b,j].y - t1) / (_dats[not b,j].y - _dats[not b,j+1].y);
+          t3:=abs(t2-_ergebnis[offset + i].x-_ergebnis[offset + i - 1].y+_ergebnis[offset + i - 1].x);
+          if (bestPos<0) or (bestErr>t3) then begin
+            bestErr:=t3;
+            bestPos:=t2;
+          end;
+          inc(j);
+        until j>_stops[not b]-_starts[not b];
+      end;
+      _ergebnis[offset + i].y:=bestPos;
+    end;
+  end;
+  for i:=_stops[false]-_starts[false]+1 to length(_ergebnis)-1 do begin
+    t1:=_ergebnis[i].x;
+    _ergebnis[i].x:=_ergebnis[i].y;
+    _ergebnis[i].y:=t1;
+  end;
+  for i:=0 to length(_ergebnis)-1 do begin
+    _ergebnis[i].y:=_ergebnis[i].y+_absShift;
+    t1:=(_ergebnis[i].y+_ergebnis[i].x)/2;
+    _ergebnis[i].y:=(_ergebnis[i].y-_ergebnis[i].x)/2;
+    _ergebnis[i].x:=t1;
+  end;
+  fertig:=true;
+end;
+
+end.
+