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); overload;
procedure integrate(indat: tExtPointArray; out outdat: tExtPointArray); overload;
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);
procedure berechneRefPuls(inPuls,surTraj: tExtPointArray; out cRefPuls: 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 integrate(indat: tExtPointArray; out outdat: tExtPointArray);
var i: longint;
begin
  setlength(outdat,length(indat));
  outdat[0]:=indat[0];
  for i:=1 to length(indat)-1 do begin
    outdat[i].x:=indat[i].x;
    outdat[i].y:=indat[i].y * (indat[i].x-indat[i-1].x)+outdat[i-1].y;
  end;
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;

procedure berechneRefPuls(inPuls,surTraj: tExtPointArray; out cRefPuls: tExtPointArray);
var i,anz: longint;
    beta:  extended;
const step=107;
begin
  i:=0;
  setlength(cRefPuls,0);
  anz:=0;
  while (i<length(surTraj)-step) and (surTraj[i].x-surTraj[i].y<inPuls[length(inPuls)-1].x) do begin
    while (i<length(surTraj)-step) and (inPuls[0].x>surTraj[i].x-surTraj[i].y) do
      inc(i);
    if i>=length(surTraj)-step then break;

    if anz>=length(cRefPuls) then
      setlength(cRefPuls,anz+32768);
    beta:=max((surTraj[i+step].y-surTraj[i].y)/(surTraj[i+step].x-surTraj[i].x),-0.5);
    cRefPuls[anz].x:=surTraj[i].x+surTraj[i].y;
    cRefPuls[anz].y:=
      (inPuls[anz].y    *   (inPuls[anz+step].x - surTraj[i].x+surTraj[i].y) +
       inPuls[anz+step].y * (- inPuls[anz].x + surTraj[i].x-surTraj[i].y)) /
      (inPuls[anz+step].x - inPuls[anz].x) *
      (1-beta)/(1+beta);
    inc(anz);
    inc(i);
  end;
  setlength(cRefPuls,anz);
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.