unit romunit;

{$mode objfpc}{$H+}

interface

uses
  Classes, SysUtils, Math, matheunit, lowlevelunit;

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 ableiten(indat: tExtPointArray; out outdat: tExtPointArray; dist: longint; bound: extended);
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 interpoliere(var dat: tExtPointArray);
procedure normiere(var dat: tExtPointArray);
procedure berechneRefPuls(inPuls,surTraj: tExtPointArray; betaGlaette: longint; betaBound: extended; 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);
      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]/factor;
        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]/factor;
        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
  i: longint;
  f: textfile;
begin
  if length(dat)=0 then
    exit;
  write(stderr,'Datei '''+nam+''' schreiben ');
  assignfile(f,nam);
  rewrite(f);
  for i:=0 to length(dat)-1 do
    writeln(f,floatToStr(dat[i]['x'])+' '+floatToStr(dat[i]['y']));
  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]);

  vereineExtrema(inPuls,outPuls,exList[0],exList[2],0.5);
  vereineExtrema(inPuls,outPuls,exList[1],exList[3],0.5);

  setlength(stuecke,length(exList[0])+length(exList[1])-1);

  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];

  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;
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;
begin
  setlength(behalten,length(extrema));
  for i:=0 to length(behalten)-1 do
    behalten[i]:=false;

  mitte:=0;
  mx:=abs(dat[extrema[0]]['y']);
  for i:=1 to length(extrema)-1 do
    if abs(dat[extrema[i]]['y'])>mx then begin
      mx:=abs(dat[extrema[i]]['y']);;
      mitte:=i;
    end;
  mx:=dat[extrema[mitte]]['x'];

  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;

  j:=0;
  for i:=0 to length(behalten)-1 do
    if behalten[i] then begin
      extrema[j]:=extrema[i];
      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);
end;

procedure gesamtverschiebung(var inPuls,outPuls: tExtPointArray; var absShift: extended);
var
  oMin,i: longint;
begin
  for i:=0 to length(outPuls)-1 do
    outPuls[i]['x']:=outPuls[i]['x']-absShift;
  oMin:=0;
  while (oMin<length(outPuls)) and (outPuls[oMin]['x']<inPuls[0]['x']) do
    inc(oMin);
  if oMin>0 then
    dec(oMin);
  for i:=oMin to length(outPuls)-1 do
    outPuls[i-oMin]:=outPuls[i];
  setlength(outPuls,length(outPuls)-oMin);
  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 ableiten(indat: tExtPointArray; out outdat: tExtPointArray; dist: longint; bound: extended);
var
  i: longint;
begin
  bound:=abs(bound);
  if dist<1 then
    dist:=1;
  setlength(outdat,length(indat)-dist);
  for i:=0 to length(outdat)-1 do begin
    outdat[i]['x']:=(indat[i+dist]['x']+indat[i]['x'])/2;
//    outdat[i]['y']:=min(max((indat[i+dist]['y']-indat[i]['y'])/(indat[i+dist]['x']-indat[i]['x']),-bound),bound);
    outdat[i]['y']:=(indat[i+dist]['y']-indat[i]['y'])/(indat[i+dist]['x']-indat[i]['x']);
    if outdat[i]['y']>bound then
      outdat[i]['y']:=bound;
    if outdat[i]['y']<-bound then
      outdat[i]['y']:=-bound;
  end;
end;

procedure removeLinearOffset(var dat: tExtPointArray);
var
  i: longint;
  d: tExtPoint;
begin
  d:=extPoint(0,0);
  for i:=0 to 99 do
    d:=d + dat[length(dat)-1-i] - dat[i];
  for i:=0 to length(dat)-1 do
    dat[i]['y']:=dat[i]['y'] - d['y']/d['x']*(dat[i]['x']-dat[0]['x']);
  d['y']:=0;
  for i:=0 to 99 do
    d['y']:=d['y'] + dat[length(dat)-1-i]['y'] + dat[i]['y'];
  d['y']:=d['y']/200;
  for i:=0 to length(dat)-1 do
    dat[i]['y']:=dat[i]['y'] - d['y'];
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;

  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 interpoliere(var dat: tExtPointArray);
var
  i,j:      longint;
  tmp,xLen: extended;
  tdat:     tExtPointArray;
begin
  tmp:=0;
  xLen:=dat[length(dat)-1]['x']-dat[0]['x'];
  for i:=1 to length(dat)-1 do begin
    if dat[i]['x']<=dat[i-1]['x'] then
      Fehler('Die Daten müssen sortiert sein und dürfen keine doppelten x-Werte enthalten! ('+floattostr(tmp)+')');
    tmp:=tmp + 1/(dat[i]['x']-dat[i-1]['x']);
  end;
  tmp:=length(dat)/tmp;
  setlength(tdat,length(dat));
  for i:=0 to length(dat)-1 do
    tdat[i]:=dat[i];
  setlength(dat,min(3*length(dat),round(xLen/tmp+1)));
  j:=0;
  for i:=0 to length(dat)-1 do begin
    dat[i]['x']:=tdat[0]['x'] + xLen*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; betaGlaette: longint; betaBound: extended; out cRefPuls: tExtPointArray);
var
  iSur,iIn,anz: longint;
  beta,frac:    extended;
begin
  if betaGlaette<1 then
    betaGlaette:=1;
  setlength(cRefPuls,0);
  anz:=0;
  iIn:=0;
  betaBound:=abs(betaBound);
  for iSur:=0 to length(surTraj)-betaGlaette-1 do begin
    if inPuls[0]['x']>surTraj[iSur]['x']-surTraj[iSur]['y'] then
      continue;
    while (iIn<length(inPuls)-1) and (inPuls[iIn+1]['x']<=surTraj[iSur]['x']-surTraj[iSur]['y']) do
      inc(iIn);
    if iIn>=length(inPuls)-1 then
      break;

    if anz>=length(cRefPuls) then
      setlength(cRefPuls,anz+32768);
    beta:=min(max((surTraj[iSur+betaGlaette]['y']-surTraj[iSur]['y'])/(surTraj[iSur+betaGlaette]['x']-surTraj[iSur]['x']),-betaBound),betaBound);
   // beta:=0;
    frac:=(surTraj[iSur]['x']-surTraj[iSur]['y'] - inPuls[iIn]['x'])/(inPuls[iIn+1]['x'] - inPuls[iIn]['x']);
    cRefPuls[anz]['x']:=surTraj[iSur]['x']+surTraj[iSur]['y'];
    cRefPuls[anz]['y']:=
      (inPuls[iIn]['y']  *  (1-frac) +
       inPuls[iIn+1]['y'] * frac) *
      (1-beta)/(1+beta);
    inc(anz);
  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.