epostunit.pas: faktorisieren neu

1 files changed, 211 insertions, 0 deletions

diff --git a/epostunit.pas b/epostunit.pas
index 0d7f829..7f190c8 100644
--- a/epostunit.pas
+++ b/epostunit.pas
@@ -166,6 +166,7 @@ type
     procedure initFuerGauszFit(sT: boolean; daten: tWerte; senkrecht: boolean; adLaenge: longint; adStart,adStop: extended);
     function fitteGausze(sT: boolean; f: tMyStringList; threads: longint): boolean;
     function fitte2dGausze(sT: boolean; f: tMyStringList): boolean;
+    function faktorisiere(sT: boolean; f: tMyStringList): boolean;
     function berechneZeitfrequenzanalyse(sT: boolean; f: tMyStringList; threads: longint; quelle: tWerte; warn: tWarnStufe): boolean;
     function berechneVerzerrung(sT: boolean; f: tMyStringList; threads: longint; quelle: tWerte; warn: tWarnStufe): boolean;
     function berechneLambdaZuOmegaVerzerrung(sT: boolean; f: tMyStringList; threads: longint; quelle: tWerte): boolean;
@@ -3940,6 +3941,216 @@ begin
   result:=true;
 end;
 
+function tWerte.faktorisiere(sT: boolean; f: tMyStringList): boolean;
+var
+  bekannteBefehle:                     tMyStringList;
+  Zeit,tmpE:                           extended;
+  anzahl,i,j:                          longint;
+  s,datei:                             string;
+  pDatei:                              textFile;
+  lineOut,tmpEA:                       tExtendedArray;
+  maxima,maxima2:                      tLongintArray;
+  toepfe:                              tIntPointArray;
+  formen:                              array of tExtendedArray;
+  residuenBerechnen,auszenHorizontal:  boolean;
+begin
+  result:=false;
+  if not sT then
+    gibAus('faktorisiere ...',3);
+  bekannteBefehle:=tMyStringList.create;
+  Zeit:=now;
+  anzahl:=1;
+  datei:='';
+  residuenBerechnen:=false;
+  auszenHorizontal:=false;
+  repeat
+    if not f.metaReadln(s,true) then begin
+      gibAus('Unerwartetes Dateiende!',3);
+      exit;
+    end;
+    bekannteBefehle.clear;
+    if istDasBefehl('Ende',s,bekannteBefehle,false) then break;
+    if istDasBefehl('Anzahl:',s,bekannteBefehle,true) then begin
+      anzahl:=round(exprToFloat(sT,s));
+      continue;
+    end;
+    if istDasBefehl('Residuen berechnen',s,bekannteBefehle,false) then begin
+      residuenBerechnen:=true;
+      continue;
+    end;
+    if istDasBefehl('Datei:',s,bekannteBefehle,true) then begin
+      if datei<>'' then begin
+        gibAus('Habe bereits eine Zieldatei beim Fitten eines 2d-Gaußes!',3);
+        bekannteBefehle.free;
+        exit;
+      end;
+      datei:=s;
+      continue;
+    end;
+    if istDasBefehl('äußere Dimension:',s,bekannteBefehle,true) then begin
+      if s='horizontal' then
+        auszenHorizontal:=true
+      else if s='vertikal' then
+        auszenHorizontal:=false
+      else begin
+        gibAus(''''+s+''' ist keine gültige äußere Dimension beim Faktorisieren - ich kenne nur ''horizontal'' und ''vertikal''!',3);
+        bekannteBefehle.free;
+        exit
+      end;
+      continue;
+    end;
+    bekannteBefehle.sort;
+    gibAus('Verstehe Option '''+s+''' nicht beim Faktorisieren!'#10'Ich kenne:'#10+bekannteBefehle.text,3);
+    bekannteBefehle.free;
+    exit;
+  until false;
+  bekannteBefehle.free;
+
+  if (datei='') and not residuenBerechnen then begin
+    gibAus('Keine Ausgaben beim Faktorisieren!',3);
+    exit;
+  end;
+
+  if sT then begin
+    result:=true;
+    exit;
+  end;
+
+  case genauigkeit of
+    gSingle:
+      lineOut:=sWerte.integriereZuLineOut(not auszenHorizontal);
+    gDouble:
+      lineOut:=dWerte.integriereZuLineOut(not auszenHorizontal);
+    gExtended:
+      lineOut:=eWerte.integriereZuLineOut(not auszenHorizontal);
+  end{of case};
+
+  setLength(maxima,0);
+  i:=0;
+  for j:=0 to length(lineOut)-1 do
+    if ((j=0) or (lineOut[j]>lineOut[j-1])) and
+      ((j=length(lineOut)-1) or (lineOut[j]>lineOut[j+1])) then begin
+      if length(maxima)<=i then
+        setLength(maxima,i+128);
+      maxima[i]:=j;
+      inc(i);
+    end;
+  setLength(maxima,i);
+
+  sortiereNachDominanz(maxima,lineOut);
+  setLength(maxima,anzahl);
+
+  setLength(maxima2,anzahl);
+  setLength(tmpEA,anzahl);
+  for i:=0 to length(tmpEA)-1 do begin
+    tmpEA[i]:=maxima[i];
+    maxima2[i]:=maxima[i];
+    maxima[i]:=i;
+  end;
+  sortiereNachWert(maxima,tmpEA,false);
+  setLength(tmpEA,0);
+  setLength(toepfe,anzahl);
+  for i:=0 to length(toepfe)-1 do begin
+    if i>0 then
+      toepfe[length(toepfe)-1-i]['y']:=(maxima2[maxima[i-1]]+maxima2[maxima[i]]) div 2
+    else if auszenHorizontal then
+      toepfe[length(toepfe)-1-i]['y']:=_xSteps-1
+    else
+      toepfe[length(toepfe)-1-i]['y']:=_tSiz-1;
+    if i<length(toepfe)-1 then
+      toepfe[length(toepfe)-1-i]['x']:=(maxima2[maxima[i]]+maxima2[maxima[i+1]]) div 2 + 1
+    else
+      toepfe[length(toepfe)-1-i]['x']:=0;
+  end;
+
+  case genauigkeit of
+    gSingle:
+      sWerte.integriereTopfweise(auszenHorizontal,toepfe,formen);
+    gDouble:
+      dWerte.integriereTopfweise(auszenHorizontal,toepfe,formen);
+    gExtended:
+      eWerte.integriereTopfweise(auszenHorizontal,toepfe,formen);
+  end{of case};
+
+  // normieren
+  for i:=0 to length(toepfe)-1 do begin
+    tmpE:=0;
+    for j:=toepfe[i]['x'] to toepfe[i]['y'] do
+      tmpE:=tmpE + lineOut[j];
+    if tmpE = 0 then
+      continue;
+    tmpE:=1/tmpE;
+    for j:=0 to length(formen[i])-1 do
+      formen[i,j]:=formen[i,j] * tmpE;
+  end;
+
+  if datei<>'' then begin
+    gibAus('... speichern ...',3);
+    assignFile(pDatei,datei);
+    rewrite(pDatei);
+    if auszenHorizontal then begin
+      writeln(pDatei,'Töpfe');
+      for i:=0 to length(toepfe)-1 do
+        for j:=toepfe[i]['x'] to toepfe[i]['y'] do
+          writeln(
+            pDatei,
+            intToStr(i) + #9 +
+            myFloatToStr(transformationen.positionAufAchseZuWert(lOben,j/_xSteps)) + #9 +
+            myFloatToStr(lineOut[j])
+          );
+      writeln(pDatei);
+      for i:=0 to length(toepfe)-1 do begin
+        writeln(pDatei,'Topf '+intToStr(i));
+        for j:=0 to length(formen[i])-1 do
+          writeln(
+            pDatei,
+            myFloatToStr(transformationen.positionAufAchseZuWert(lLinks,j/_tSiz)) + #9 +
+            myFloatToStr(formen[i][j])
+          );
+        writeln(pDatei);
+      end;
+    end
+    else begin
+      writeln(pDatei,'Töpfe');
+      for i:=0 to length(toepfe)-1 do
+        for j:=toepfe[i]['x'] to toepfe[i]['y'] do
+          writeln(
+            pDatei,
+            intToStr(i) + #9 +
+            myFloatToStr(transformationen.positionAufAchseZuWert(lLinks,j/_tSiz)) + #9 +
+            myFloatToStr(lineOut[j])
+          );
+      writeln(pDatei);
+      for i:=0 to length(toepfe)-1 do begin
+        writeln(pDatei,'Topf '+intToStr(i));
+        for j:=0 to length(formen[i])-1 do
+          writeln(
+            pDatei,
+            myFloatToStr(transformationen.positionAufAchseZuWert(lOben,j/_xSteps)) + #9 +
+            myFloatToStr(formen[i][j])
+          );
+        writeln(pDatei);
+      end;
+    end;
+    closeFile(pDatei);
+  end;
+
+  if residuenBerechnen then begin
+    gibAus('... Residuen berechnen ...',3);
+    case genauigkeit of
+      gSingle:
+        sWerte.produktSubtrahieren(auszenHorizontal,toepfe,lineOut,formen);
+      gDouble:
+        dWerte.produktSubtrahieren(auszenHorizontal,toepfe,lineOut,formen);
+      gExtended:
+        eWerte.produktSubtrahieren(auszenHorizontal,toepfe,lineOut,formen);
+    end{of case};
+  end;
+
+  gibAus('... fertig '+timetostr(now-Zeit),3);
+  result:=true;
+end;
+
 function tWerte.berechneZeitfrequenzanalyse(sT: boolean; f: tMyStringList; threads: longint; quelle: tWerte; warn: tWarnStufe): boolean;
 var
   i,tMin,tMax,qlen:    longint;