1 files changed, 218 insertions, 0 deletions
diff --git a/harmoneff/myMath.pas b/harmoneff/myMath.pas
new file mode 100644
index 0000000..04681d1
--- /dev/null
+++ b/harmoneff/myMath.pas
@@ -0,0 +1,218 @@
+unit myMath;
+
+interface
+
+uses math, sysutils;
+
+type TExtPoint = record
+                  case byte of
+                   0: (x,y: extended;);     // coordinates or
+                   1: (val,err: extended;); // value and relative error
+                 end;
+
+function myfloattostr(f: double): string;
+
+function multiply(a,b: TExtPoint): TExtPoint;
+function divide(a,b: TExtPoint): TExtPoint;
+function valuetostr(a: TExtPoint; separator: String): String;
+function valuetostr(a: TExtPoint): String;
+function valuetostr(a: TExtPoint; separator: String; mitFehler: boolean): String;
+function valuetostr(a: TExtPoint; mitFehler: boolean): String;
+
+procedure GaussFitDownhillSimplex(W: array of TExtPoint; var y0,dt,delta: extended);
+procedure GaussFitGradientenverfahren(W: array of TExtPoint; var y0,dt,delta: extended);
+
+var Komma: char = ',';
+
+implementation
+
+function myfloattostr(f: double): string;
+var s,t: string;
+begin
+ s:=floattostr(f);
+ t:='';
+ while pos('.',s)>0 do
+  begin
+   t:=t+copy(s,1,pos('.',s)-1)+komma;
+   delete(s,1,pos('.',s));
+  end;
+ myfloattostr:=t+s;
+end;
+
+//*****************************************************************************
+
+function multiply(a,b: TExtPoint): TExtPoint;
+var erg: TExtPoint;
+begin
+ erg.val:=a.val*b.val;
+ erg.err:=a.err+b.err;
+ multiply:=erg;
+end;
+
+function divide(a,b: TExtPoint): TExtPoint;
+var erg: TExtPoint;
+begin
+ if b.val=0 then writeln('Division durch Null!');
+ erg.val:=a.val/b.val;
+ erg.err:=a.err+b.err;
+ divide:=erg;
+end;
+
+function valuetostr(a: TExtPoint; separator: String): String;
+begin
+ valuetostr:=valuetostr(a,separator,true);
+end;
+
+function valuetostr(a: TExtPoint): String;
+begin
+ valuetostr:=valuetostr(a,#9);
+end;
+
+function valuetostr(a: TExtPoint; separator: String; mitFehler: boolean): String;
+begin
+ if mitFehler then valuetostr:=myfloattostr(a.val)+separator+myfloattostr(a.val*a.err)+separator+myfloattostr(a.err)
+              else valuetostr:=myfloattostr(a.val);
+end;
+
+function valuetostr(a: TExtPoint; mitFehler: boolean): String;
+begin
+ valuetostr:=valuetostr(a,#9,mitFehler);
+end;
+
+//*****************************************************************************
+
+function qerr(W: array of TExtPoint; Param: TExtPoint): extended;
+var I:   Integer;
+    erg: extended;
+begin
+ erg:=0;
+ For I:=0 to length(W)-1 do
+  erg:=erg +
+       sqr(W[I].y - Param.y*exp(-sqr(W[I].x/Param.x)));
+ qerr:=erg;
+end;
+
+procedure GaussFitDownhillSimplex(W: array of TExtPoint; var y0,dt,delta: extended);
+var Simplex:    array[0..2] of TExtPoint;
+    Errs:       array[0..2] of Extended;
+    mitt,np:    TExtPoint;
+    ma,mi,nerr: extended;
+    I,J:        Integer;
+const gutFaktor = 2.5;
+      schlechtFaktor = 0.5;
+      mittelFaktor = 1.0;
+begin
+ Simplex[0].x:=dt;
+ Simplex[0].y:=y0;
+ Simplex[1].x:=dt+delta*1000;
+ Simplex[1].y:=y0;
+ Simplex[2].x:=dt;
+ Simplex[2].y:=y0+delta*1000;
+ For I:=0 to 2 do
+  Errs[I]:=qerr(W,Simplex[I]);
+ repeat
+  J:=0;
+  For I:=1 to 2 do
+   if Errs[I]>Errs[J] then
+    J:=I;
+  mitt.x:=0;
+  mitt.y:=0;
+  Ma:=Errs[Byte(J=0)];
+  Mi:=Errs[Byte(J=0)];
+  For I:=0 to 2 do
+   if J<>I then
+    begin
+     Ma:=Max(Ma,Errs[I]);
+     Mi:=Min(Mi,Errs[I]);
+     mitt.x:=mitt.x+Simplex[I].x;
+     mitt.y:=mitt.y+Simplex[I].y;
+    end;
+  mitt.x:=mitt.x/2;
+  mitt.y:=mitt.y/2;
+  np.x:=2*mitt.x-Simplex[J].x;
+  np.y:=2*mitt.y-Simplex[J].y;
+  nerr:=qerr(W,np);
+  np.x:=mitt.x-Simplex[J].x;
+  np.y:=mitt.y-Simplex[J].y;
+  if nerr<Mi then // Minimum !
+   begin
+    Simplex[J].x:=mitt.x+gutFaktor*np.x;
+    Simplex[J].y:=mitt.y+gutFaktor*np.y;
+   end
+    else
+   if nerr>Ma then // Maximum !
+    begin
+     Simplex[J].x:=mitt.x+schlechtFaktor*np.x;
+     Simplex[J].y:=mitt.y+schlechtFaktor*np.y;
+    end
+     else
+    begin
+     Simplex[J].x:=mitt.x+mittelFaktor*np.x;
+     Simplex[J].y:=mitt.y+mittelFaktor*np.y;
+    end;
+  Errs[J]:=qerr(W,Simplex[J]);  
+ until (abs((Simplex[0].x-Simplex[1].x)*(Simplex[0].y-Simplex[2].y) -
+            (Simplex[0].x-Simplex[2].x)*(Simplex[0].y-Simplex[1].y)) <= delta);
+ J:=0;
+ For I:=1 to 2 do
+  if Errs[I]<Errs[J] then
+   J:=I;
+ delta:=Errs[J];
+ dt:=Simplex[J].x;
+ y0:=Simplex[J].y;
+end;
+
+procedure GaussFitGradientenverfahren(W: array of TExtPoint; var y0,dt,delta: extended);
+var ex,tmp,err: extended;
+    i,vz,vzz:   integer;
+    dy0,
+    ddt,fac:    extended;
+    weiter:     boolean;
+begin
+ dt:=1/dt;
+ err:=0;
+ fac:=2;
+ vz:=0;
+ vzz:=0;
+ for i:=0 to length(W)-1 do
+  err:=err+sqr(W[i].y-y0*exp(-sqr(W[i].x*dt)));
+ repeat
+  dy0:=0;
+  ddt:=0;
+  for i:=0 to length(W)-1 do
+   begin
+    ex:=exp(-sqr(W[i].x*dt));
+    tmp:=W[i].y*ex-y0*sqr(ex);
+    dy0:=dy0 - tmp;
+    ddt:=ddt + sqr(W[i].x)*tmp;
+   end;
+  ddt:=ddt*dt*y0*2;
+  weiter:=sqr(ddt)+sqr(dy0)>1;
+  if weiter then tmp:=1/sqrt(sqr(ddt)+sqr(dy0))
+            else tmp:=1;
+  if (vz<>2*byte(ddt>0)+byte(dy0>0)) or
+     (ddt=0) or (dy0=0) then
+   begin
+    fac:=fac/2;
+    vzz:=0;
+   end
+    else
+   begin
+    inc(vzz);
+    for i:=5 to vzz do
+     fac:=min(1,fac*1.2);
+   end;
+  vz:=2*byte(ddt>0)+byte(dy0>0);
+  tmp:=tmp*fac;
+  dt:=dt-ddt*tmp;
+  y0:=y0-dy0*tmp;
+  tmp:=err;
+  err:=0;
+  for i:=0 to length(W)-1 do
+   err:=err+sqr(W[i].y-y0*exp(-sqr(W[i].x*dt)));
+ until (tmp-err<delta*err) and (fac<delta);
+ dt:=1/dt;
+ delta:=err;
+end;
+
+end.