diff --git a/math/boxkde/GausLib.cpp b/math/boxkde/GausLib.cpp
new file mode 100644
index 0000000..8c4b98f
--- /dev/null
+++ b/math/boxkde/GausLib.cpp
@@ -0,0 +1,83 @@
+#include "GausLib.h"
+
+#include "benchmark.h"
+#include "DataStructures.h"
+#include "Image2D.h"
+#include "BoxGaus.h"
+#include "Grid2D.h"
+#include <cstring>
+
+// X sind 2D Samples [x,x,x,.... y,y,y]
+// sizeX = Anzahl Samples
+// hx, hy Bandbreite KDE (default=100)
+// nBins = Anzahl der Bins in die jeweilige Richtung
+// boundingBox ist Größe des Gebäudes (40 m x 60 m)
+// algorithm = define method  
+
+void gausKde2D_simple(float* X, int sizeX, float* weights, float hx, float hy, int nBinsX, int nBinsY, c_bbox* boundingBox, int algorithm, double* out_maxValue, double* out_maxPosX, double* out_maxPosY, double* out_runTimeInNS, float* out_density)
+{
+    BenchResult bench = benchmarkEx("", 1, [&](){
+        // Create bounding box
+        BoundingBox<float> bb(boundingBox->minX, boundingBox->maxX, boundingBox->minY, boundingBox->maxY);
+
+        // Create histogram
+        Grid2D<float> grid(bb, nBinsX, nBinsY);
+        for (size_t i = 0; i < sizeX; i++)
+        {
+            float x = X[i];
+            float y = X[i + sizeX];
+            float w = weights[i];
+
+            grid.add(x, y, w);
+        }
+
+        int nFilt = 3;
+
+        float sigmaX = hx / grid.binSizeX;
+        float sigmaY = hy / grid.binSizeY;
+
+        // Apply gaus filter
+        switch (algorithm)
+        {
+            default:
+            case ALGO_BOX:
+            {
+                BoxGaus<float> boxGaus;
+                boxGaus.approxGaus(grid.image(), sigmaX, sigmaY, nFilt);
+                break;
+            }
+            case ALGO_EXBOX:
+            {
+                //ExtendedBox<float> boxGaus;
+                //boxGaus.approxGaus(grid.image(), sigmaX, sigmaY, nFilt);
+                break;
+            }
+            case ALGO_BOX_SIMD:
+            {
+                //BoxGausSIMD<float> boxGaus;
+                //boxGaus.approxGaus(grid.image(), sigmaX, sigmaY, nFilt);
+                break;
+            }
+            case ALGO_BOX_CL:
+            {
+                //clBox->writeBuffer(grid.image().data().data());
+                //clBox->execute(sigmaX, nFilt);
+                //clBox->readBuffer(grid.image().data().data());
+                break;
+            }
+        }
+
+
+        Point2D<float> maxPos;
+        *out_maxValue = grid.maximum(maxPos);
+        *out_maxPosX = static_cast<double>(maxPos.X);
+        *out_maxPosY = static_cast<double>(maxPos.Y);
+
+        if (out_density)
+        {
+            memcpy(out_density, grid.image().data().data(), sizeof(float)*grid.image().data().size());
+        }
+    });
+
+    *out_runTimeInNS = bench.min;
+}