statistics.pas new

3 files changed, 170 insertions, 0 deletions

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..5c9870e
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,2 @@
+lib
+statistics
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..8efab47
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,11 @@
+statistics: statistics.pas
+	mkdir -p lib
+	fpc -Mobjfpc -FUlib -Fu/usr/lib/fpc/src/packages/* -Fu../units "$<"
+
+run: statistics
+	find ~/tmp/ -name '*.gpx' \
+	| sort \
+	| tr '\n' '\0' \
+	| xargs -0r cat \
+	| ./gpx-to-tsv \
+	| ./statistics
diff --git a/statistics.pas b/statistics.pas
new file mode 100644
index 0000000..525b353
--- /dev/null
+++ b/statistics.pas
@@ -0,0 +1,157 @@
+program statistics;
+
+uses
+  lowlevelunit, SysUtils, math;
+
+type
+  tPosition = class
+    breite,laenge,hoehe: extended;
+    zeit:                int64;
+    constructor create;
+    destructor destroy; override;
+    procedure load(s: ansiString);
+    function toString: ansiString;
+    procedure copyFrom(p: tPosition);
+  end;
+  tMesswert = class(tPosition)
+    wert: extended;
+  end;
+
+constructor tPosition.create;
+begin
+  inherited create;
+end;
+
+destructor tPosition.destroy;
+begin
+  inherited destroy;
+end;
+
+procedure tPosition.load(s: ansiString);
+begin
+  breite:=strToFloat(erstesArgument(s,#9));
+  laenge:=strToFloat(erstesArgument(s,#9));
+  hoehe:=strToFloat(erstesArgument(s,#9));
+  zeit:=round(strToDate(erstesArgument(s,'T')));
+  zeit:=round(24*60*60*(zeit + strToTime(erstesArgument(s,'Z'))));
+end;
+
+function tPosition.toString: ansiString;
+begin
+  result:=
+    '"'+floatToStr(breite)+'" "'+
+    floatToStr(laenge)+'" "'+
+    floatToStr(hoehe)+'" "'+
+    floatToStr(zeit)+'"';
+end;
+
+procedure tPosition.copyFrom(p: tPosition);
+begin
+  breite:=p.breite;
+  laenge:=p.laenge;
+  hoehe:=p.hoehe;
+  zeit:=p.zeit;
+end;
+
+function abstand(p1,p2: tPosition): extended;
+begin
+  result:=sqrt(
+    sqr(p1.hoehe-p2.hoehe) +
+    sqr((12756320 + p1.hoehe + p2.hoehe)/2) * (
+      sqr((p1.breite-p2.breite)/180*pi) +
+      sqr(
+        sin((p1.breite+p2.breite)/2) *
+        (p1.laenge - p2.laenge)
+      )
+    )
+  )
+end;
+
+var
+  b:                   boolean;
+  s:                   ansiString;
+  i,x,y:               int64;
+  ps:                  array[boolean] of tPosition;
+  ws:                  array of tMesswert;
+  minB,maxB,minL,maxL: extended;
+  zoom:                extended;
+  wL,wU,wBr,wHo:       int64;
+  anz:                 array of int64;
+  sum,qSum:            array of extended;
+
+begin
+  i:=0;
+  ShortDateFormat:='y-m-d';
+  zoom:=500;
+  for b:=false to true do
+    ps[b]:=tPosition.create;
+  setLength(ws,0);
+  minB:=0;
+  maxB:=0;
+  minL:=0;
+  maxL:=0;
+  while not eof() do begin
+    readln(s);
+    ps[true].copyFrom(ps[false]);
+    ps[false].load(s);
+    if i>0 then begin
+      if ps[false].zeit<ps[true].zeit then
+        writeln(stderr, 'ERROR: ' + ps[false].toString + ' < ' + ps[true].toString);
+      if ps[false].zeit-ps[true].zeit < 30 then begin
+        setLength(ws,length(ws)+1);
+        ws[length(ws)-1]:=tMesswert.create;
+        with ws[length(ws)-1] do begin
+          copyFrom(ps[true]);
+          wert:=abstand(ps[true],ps[false])/(ps[false].zeit-ps[true].zeit);
+          if (length(ws)=1) or (minB>breite) then
+            minB:=breite;
+          if (length(ws)=1) or (maxB<breite) then
+            maxB:=breite;
+          if (length(ws)=1) or (minL>laenge) then
+            minL:=laenge;
+          if (length(ws)=1) or (maxL<laenge) then
+            maxL:=laenge;
+        end;
+      end;
+    end;
+    inc(i);
+    if i>10000 then
+      break;
+  end;
+  for b:=false to true do
+    ps[b].free;
+  writeln(minB,' .. ',maxB,' x ',minL,' .. ',maxL);
+  wU:=round(floor(minB*zoom));
+  wHo:=round(floor(maxB*zoom)+1)-wU;
+  wL:=round(floor(minL*cos(minB+maxB)*zoom));
+  wBr:=round(floor(maxL*cos(minB+maxB)*zoom)+1)-wL;
+  writeln(wBr,' x ',wHo);
+  setLength(anz,wBr*wHo);
+  setLength(sum,wBr*wHo);
+  setLength(qSum,wBr*wHo);
+  fillchar(anz[0],sizeof(anz[0])*length(anz),0);
+  fillchar(sum[0],sizeof(sum[0])*length(sum),0);
+  fillchar(qSum[0],sizeof(qSum[0])*length(qSum),0);
+  for i:=0 to length(ws)-1 do begin
+    x:=round(ws[i].laenge*cos(minB+maxB)*zoom)-wL;
+    y:=round(ws[i].breite*zoom)-wU;
+    inc(anz[x+wBr*y]);
+    sum[x+wBr*y]:= sum[x+wBr*y] + ws[i].wert;
+    qSum[x+wBr*y]:= qSum[x+wBr*y] + sqr(ws[i].wert);
+  end;
+  for y:=0 to wHo-1 do begin
+    for x:=0 to wBr-1 do begin
+      if x>0 then
+        write(',');
+      write(intToStr(anz[x+y*wBr])+':');
+      if anz[x+y*wBr]>0 then begin
+        write(floatToStr(sum[x+y*wBr]/anz[x+y*wBr]));
+        write('±');
+        write(floatToStr(sqrt(qSum[x+y*wBr]/anz[x+y*wBr]-sqr(sum[x+y*wBr]/anz[x+y*wBr]))));
+      end
+      else
+        write('NaN±NaN');
+    end;
+    writeln;
+  end;
+end.