auch reflektierten Puls anhand von Trajektorie rekonstruieren

... bisher unvollständig

5 files changed, 114 insertions, 61 deletions

diff --git a/Bild.pdf b/Bild.pdf
index 52cdcb7..ae01189 100644
--- a/Bild.pdf
+++ b/Bild.pdf
diff --git a/ROM.lpr b/ROM.lpr
index 5863de0..26aebb3 100644
--- a/ROM.lpr
+++ b/ROM.lpr
@@ -12,17 +12,17 @@ uses
   { 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;
+var inPulsO,inPuls,refPuls,surTraj,cRefPuls: 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 '+
+const Verwendung='Verwendung: ROM ($Einfallspuls_Datei $Ausfallspuls_Datei)/(- $trace-Datei-Prefix) $output_inPuls $output_refPuls $output_Trajektorie $output_cRefPuls '+
    '[-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;
+  if paramcount<6 then Fehler(Verwendung);
+  i:=7;
   force:=false;
   smooth:=1;
   tmax:=-1;
@@ -79,23 +79,23 @@ begin
     Fehler(Verwendung);
   end;
   if paramstr(1)='-' then
-    readRawInputs(paramstr(2),inPuls,refPuls,absShift)
+    readRawInputs(paramstr(2),inPulsO,refPuls,absShift)
   else begin
-    readTextInput(paramstr(1),inPuls);
+    readTextInput(paramstr(1),inPulsO);
     readTextInput(paramstr(2),refPuls);
   end;
   write(stderr,'Input sortieren ...');
-  sort(inPuls);
+  sort(inPulsO);
   sort(refPuls);
   writeln(stderr,' fertig');
-  uniq(inPuls,false);
+  uniq(inPulsO,false);
   uniq(refPuls,false);
   write(stderr,'Input interpolieren ...');
-  interpoliere(inPuls);
+  interpoliere(inPulsO);
   interpoliere(refPuls);
   writeln(stderr,' fertig');
   flip(refPuls);
-  integrate(inPuls);
+  integrate(inPulsO,inPuls);
   integrate(refPuls);
   removeLinearOffset(inPuls);
   removeLinearOffset(refPuls);
@@ -142,9 +142,16 @@ begin
     normiere(refPuls);
     normiere(surTraj);
     writeln(stderr,' fertig');
+  end
+  else begin
+    write(stderr,'Reflektierten Puls berechnen ...');
+    berechneRefPuls(inPulsO,surTraj,cRefPuls);
+    integrate(cRefPuls);
+    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);
+  if (paramstr(6)<>'-') and not fourier then writeOutput(paramstr(6),cRefPuls);
 end.
 
diff --git a/ROM.lps b/ROM.lps
index b6e2d2d..99565a0 100644
--- a/ROM.lps
+++ b/ROM.lps
@@ -8,24 +8,24 @@
         <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"/>
+        <TopLine Value="114"/>
+        <CursorPos X="25" Y="149"/>
+        <UsageCount Value="91"/>
         <Loaded Value="True"/>
       </Unit0>
       <Unit1>
         <Filename Value="romunit.pas"/>
         <IsPartOfProject Value="True"/>
         <UnitName Value="romunit"/>
+        <IsVisibleTab Value="True"/>
         <EditorIndex Value="1"/>
         <WindowIndex Value="0"/>
-        <TopLine Value="1004"/>
-        <CursorPos X="1" Y="998"/>
-        <FoldState Value=" T3iA041 pkRkZ0Y113A"/>
-        <UsageCount Value="90"/>
+        <TopLine Value="974"/>
+        <CursorPos X="11" Y="1001"/>
+        <FoldState Value=" T3iC041 pkRkZ0Y2]9clZ093 piXj70G]Rj3jN0w6]9UlG0S[I4Ck10P116k"/>
+        <UsageCount Value="91"/>
         <Loaded Value="True"/>
       </Unit1>
       <Unit2>
@@ -34,132 +34,136 @@
         <UnitName Value="mathunit"/>
         <EditorIndex Value="2"/>
         <WindowIndex Value="0"/>
-        <TopLine Value="4"/>
+        <TopLine Value="3"/>
         <CursorPos X="1" Y="43"/>
-        <UsageCount Value="89"/>
+        <UsageCount Value="90"/>
         <Loaded Value="True"/>
       </Unit2>
     </Units>
     <General>
       <ActiveWindowIndexAtStart Value="0"/>
     </General>
-    <JumpHistory Count="29" HistoryIndex="28">
+    <JumpHistory Count="30" HistoryIndex="29">
       <Position1>
-        <Filename Value="ROM.lpr"/>
-        <Caret Line="93" Column="58" TopLine="77"/>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="26" Column="95" TopLine="1"/>
       </Position1>
       <Position2>
-        <Filename Value="ROM.lpr"/>
-        <Caret Line="123" Column="66" TopLine="90"/>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="727" Column="6" TopLine="694"/>
       </Position2>
       <Position3>
         <Filename Value="romunit.pas"/>
-        <Caret Line="577" Column="69" TopLine="393"/>
+        <Caret Line="26" Column="1" TopLine="1"/>
       </Position3>
       <Position4>
         <Filename Value="romunit.pas"/>
-        <Caret Line="212" Column="2" TopLine="38"/>
+        <Caret Line="334" Column="31" TopLine="328"/>
       </Position4>
       <Position5>
-        <Filename Value="mathunit.pas"/>
-        <Caret Line="22" Column="1" TopLine="1"/>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="311" Column="41" TopLine="288"/>
       </Position5>
       <Position6>
         <Filename Value="romunit.pas"/>
-        <Caret Line="528" Column="14" TopLine="476"/>
+        <Caret Line="310" Column="3" TopLine="291"/>
       </Position6>
       <Position7>
         <Filename Value="romunit.pas"/>
-        <Caret Line="519" Column="37" TopLine="475"/>
+        <Caret Line="22" Column="1" TopLine="16"/>
       </Position7>
       <Position8>
         <Filename Value="romunit.pas"/>
-        <Caret Line="542" Column="1" TopLine="512"/>
+        <Caret Line="1018" Column="1" TopLine="891"/>
       </Position8>
       <Position9>
         <Filename Value="romunit.pas"/>
-        <Caret Line="536" Column="34" TopLine="514"/>
+        <Caret Line="669" Column="16" TopLine="625"/>
       </Position9>
       <Position10>
         <Filename Value="romunit.pas"/>
-        <Caret Line="537" Column="25" TopLine="508"/>
+        <Caret Line="621" Column="38" TopLine="527"/>
       </Position10>
       <Position11>
         <Filename Value="romunit.pas"/>
-        <Caret Line="679" Column="28" TopLine="649"/>
+        <Caret Line="312" Column="1" TopLine="295"/>
       </Position11>
       <Position12>
         <Filename Value="romunit.pas"/>
-        <Caret Line="333" Column="150" TopLine="305"/>
+        <Caret Line="657" Column="15" TopLine="626"/>
       </Position12>
       <Position13>
         <Filename Value="romunit.pas"/>
-        <Caret Line="310" Column="3" TopLine="53"/>
+        <Caret Line="655" Column="16" TopLine="637"/>
       </Position13>
       <Position14>
         <Filename Value="romunit.pas"/>
-        <Caret Line="26" Column="95" TopLine="1"/>
+        <Caret Line="285" Column="9" TopLine="254"/>
       </Position14>
       <Position15>
         <Filename Value="romunit.pas"/>
-        <Caret Line="727" Column="6" TopLine="694"/>
+        <Caret Line="312" Column="49" TopLine="301"/>
       </Position15>
       <Position16>
         <Filename Value="romunit.pas"/>
-        <Caret Line="26" Column="1" TopLine="1"/>
+        <Caret Line="21" Column="69" TopLine="1"/>
       </Position16>
       <Position17>
-        <Filename Value="romunit.pas"/>
-        <Caret Line="334" Column="31" TopLine="328"/>
+        <Filename Value="ROM.lpr"/>
+        <Caret Line="84" Column="38" TopLine="15"/>
       </Position17>
       <Position18>
-        <Filename Value="romunit.pas"/>
-        <Caret Line="311" Column="41" TopLine="288"/>
+        <Filename Value="ROM.lpr"/>
+        <Caret Line="98" Column="3" TopLine="23"/>
       </Position18>
       <Position19>
         <Filename Value="romunit.pas"/>
-        <Caret Line="310" Column="3" TopLine="291"/>
+        <Caret Line="25" Column="1" TopLine="6"/>
       </Position19>
       <Position20>
         <Filename Value="romunit.pas"/>
-        <Caret Line="22" Column="1" TopLine="16"/>
+        <Caret Line="726" Column="60" TopLine="285"/>
       </Position20>
       <Position21>
-        <Filename Value="romunit.pas"/>
-        <Caret Line="1018" Column="1" TopLine="891"/>
+        <Filename Value="ROM.lpr"/>
+        <Caret Line="98" Column="3" TopLine="77"/>
       </Position21>
       <Position22>
         <Filename Value="romunit.pas"/>
-        <Caret Line="669" Column="16" TopLine="625"/>
+        <Caret Line="724" Column="24" TopLine="688"/>
       </Position22>
       <Position23>
         <Filename Value="romunit.pas"/>
-        <Caret Line="621" Column="38" TopLine="527"/>
+        <Caret Line="33" Column="1" TopLine="22"/>
       </Position23>
       <Position24>
         <Filename Value="romunit.pas"/>
-        <Caret Line="312" Column="1" TopLine="295"/>
+        <Caret Line="1099" Column="1" TopLine="1057"/>
       </Position24>
       <Position25>
         <Filename Value="romunit.pas"/>
-        <Caret Line="657" Column="15" TopLine="626"/>
+        <Caret Line="1001" Column="69" TopLine="946"/>
       </Position25>
       <Position26>
         <Filename Value="romunit.pas"/>
-        <Caret Line="655" Column="16" TopLine="637"/>
+        <Caret Line="987" Column="11" TopLine="947"/>
       </Position26>
       <Position27>
         <Filename Value="romunit.pas"/>
-        <Caret Line="285" Column="9" TopLine="254"/>
+        <Caret Line="32" Column="33" TopLine="22"/>
       </Position27>
       <Position28>
         <Filename Value="romunit.pas"/>
-        <Caret Line="312" Column="49" TopLine="301"/>
+        <Caret Line="992" Column="42" TopLine="947"/>
       </Position28>
       <Position29>
         <Filename Value="romunit.pas"/>
-        <Caret Line="21" Column="69" TopLine="1"/>
+        <Caret Line="1003" Column="53" TopLine="973"/>
       </Position29>
+      <Position30>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="1001" Column="79" TopLine="719"/>
+      </Position30>
     </JumpHistory>
   </ProjectSession>
   <EditorMacros Count="0"/>
diff --git a/make.gnu b/make.gnu
index 9a85259..f3736e2 100644
--- a/make.gnu
+++ b/make.gnu
@@ -22,5 +22,6 @@ else {
 #    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"
+     "~/Dokumente/Prograemmchen/ROM/output.dat" using 1:2 with lines lt 3 title "output", \
+     "~/Dokumente/Prograemmchen/ROM/outputC.dat" using 1:2 with lines lt 4 title "output reconstructed"
 
diff --git a/romunit.pas b/romunit.pas
index aaab446..9437006 100644
--- a/romunit.pas
+++ b/romunit.pas
@@ -20,7 +20,8 @@ function  findeExtrema(const dat: tExtPointArray; maxima: boolean): tLongintArra
 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 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);
@@ -28,6 +29,7 @@ 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)
@@ -715,6 +717,17 @@ begin
     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;
@@ -971,6 +984,34 @@ begin
     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);