diff --git a/CMakeLists.txt b/CMakeLists.txt
new file mode 100755
index 0000000000000000000000000000000000000000..a2e3b0ed3f7107b5457073a12cec43c272eac057
--- /dev/null
+++ b/CMakeLists.txt
@@ -0,0 +1,13 @@
+cmake_minimum_required(VERSION 3.13) # 3.13 is required for target_link_options
+project(CppFrep CXX)
+
+if(NOT CMAKE_BUILD_TYPE)
+    set(CMAKE_BUILD_TYPE "Release")
+endif()
+message(STATUS "Build type: ${CMAKE_BUILD_TYPE}")
+
+find_package (Eigen3 3.3 REQUIRED NO_MODULE)
+
+include(cmake/shared_settings.cmake)
+
+add_subdirectory(cpp_frep)
diff --git a/README.md b/README.md
index db59d47490c55625bdf25d2c405a2931e4fc2ecf..f68799edc481ffc990559b6e085d698f71bb8904 100644
--- a/README.md
+++ b/README.md
@@ -1,2 +1,34 @@
 # cpp_frep
 
+A simple c++ frep library to design custom functional representation of geometries.
+
+## Dependencies
+
+Similar to [cpp_starter](https://gitlab.cba.mit.edu/erik/cpp_starter/tree/master):
+- a modern C++ compiler (i.e. one that supports C++17 features, ex. gcc 7+)
+- cmake 3.13 or newer
+- Eigen3 (installable via apt, yum, or homebrew)
+
+## Building
+
+From the project's root directory,
+
+```
+mkdir build
+cd build
+cmake ..
+make
+```
+
+Tips:
+- you can build using multiple threads at once using e.g. `make -j4`
+- you can specify a specific compiler using `cmake .. -DCMAKE_CXX_COMPILER=g++-8`
+- for debugging symbols use `cmake .. -DCMAKE_BUILD_TYPE=Debug`
+
+## Running
+
+From the build directory,
+
+```
+cpp_frep/cpp_frep
+```
diff --git a/cmake/shared_settings.cmake b/cmake/shared_settings.cmake
new file mode 100755
index 0000000000000000000000000000000000000000..828fac10346e1139ee7b5b180c07a92b996d0ff0
--- /dev/null
+++ b/cmake/shared_settings.cmake
@@ -0,0 +1,19 @@
+# This file defines an interface library used to add common compile flags to all libraries and
+# executables.
+
+#---------------------------------------------------------------------------------------------------
+add_library(shared_settings INTERFACE)
+
+# Speed flags
+target_compile_options(shared_settings INTERFACE -march=native -ffast-math)
+
+# Warning flags
+target_compile_options(shared_settings INTERFACE
+    -Wall
+    -Wcast-align
+    -Wcast-qual
+    -Wextra
+    -Wundef
+    -Wzero-as-null-pointer-constant
+    -pedantic
+)
diff --git a/cpp_frep/CMakeLists.txt b/cpp_frep/CMakeLists.txt
new file mode 100755
index 0000000000000000000000000000000000000000..16627a6f1642a54631cf932fb80ae60903225d97
--- /dev/null
+++ b/cpp_frep/CMakeLists.txt
@@ -0,0 +1,11 @@
+add_library(cpp_frep_lib
+    frep.cpp
+)
+target_include_directories(cpp_frep_lib PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
+target_link_libraries(cpp_frep_lib shared_settings Eigen3::Eigen)
+target_compile_features(cpp_frep_lib PUBLIC cxx_std_17)
+
+add_executable(cpp_frep
+    main.cpp
+)
+target_link_libraries(cpp_frep shared_settings cpp_frep_lib)
diff --git a/cpp_frep/frep.cpp b/cpp_frep/frep.cpp
new file mode 100755
index 0000000000000000000000000000000000000000..da3ac736590e6224fff19beb293e1a3fcb37608e
--- /dev/null
+++ b/cpp_frep/frep.cpp
@@ -0,0 +1,485 @@
+#include <iostream>
+#include <math.h>
+
+namespace cpp_frep {
+    using scalar=float;
+
+    class Frep {
+        public:
+            scalar X;
+            scalar Y;
+            std::string name;
+
+            Frep():X(0.0f),Y(0.0f){};
+            Frep(std::string name_,scalar X_, scalar Y_) {
+                X=X_;
+                Y=Y_;
+                name=name_;
+            };
+
+            Frep get(){
+                return *this;
+            }
+
+            void setXY(scalar X_,scalar Y_){
+                X=X_;
+                Y=Y_;
+            }
+            virtual scalar eval(){
+                return 0.0f;
+            }
+
+        // private:
+        
+    };
+    //..........................SHAPES..................................//
+    class FrepPrimitive: public Frep {
+        public:
+            scalar center_x;
+            scalar center_y;
+
+            FrepPrimitive():center_x(0.0f),center_y(0.0f){};
+            FrepPrimitive(scalar center_x_, scalar center_y_) {
+                center_x=center_x_;
+                center_y=center_y_;
+            };
+
+            FrepPrimitive get(){
+                return *this;
+            }
+
+            virtual scalar eval(){
+                return 0.0f;
+            }
+
+            
+    };
+
+    class Circle: public FrepPrimitive {
+        public:
+            scalar radius;
+
+            Circle():radius(1.0f){};
+            Circle(scalar radius_,scalar center_x_, scalar center_y_) {
+                radius=radius_;
+                center_x=center_x_;
+                center_y=center_y_;
+            };
+
+            scalar eval(){
+                return ((radius)*(radius)-((X-(center_x))*(X-(center_x))+(Y-(center_y))*(Y-(center_y))));
+            }
+
+            Circle get(){
+                return *this;
+            }
+
+
+        // private:
+            
+    };
+
+    class Rectangle: public FrepPrimitive {
+        public:
+            scalar xmin;
+            scalar xmax;
+            scalar ymin;
+            scalar ymax;
+
+            Rectangle():xmin (0.0f), xmax (0.0f), ymin (0.0f), ymax (0.0f) {};
+            Rectangle(scalar xmin_, scalar xmax_,scalar ymin_,scalar ymax_) {
+                xmin=xmin_;
+                xmax=xmax_;
+                ymin=ymin_;
+                ymax=ymax_;
+            };
+
+            scalar eval(){
+                scalar xfn = std::min(X-(xmin),(xmax)-X);
+                scalar yfn = std::min(Y-(ymin),(ymax)-Y);
+                scalar fn = std::min((xfn),(yfn));
+                return fn;
+            }
+
+            Rectangle get(){
+                return *this;
+            }
+        // private:
+            
+    };
+    //..................................................................//
+
+    //.......................OPERATORS..................................//
+    class FrepOperators: public Frep {
+        public:
+            
+            FrepOperators(){};  
+    };
+
+    class Add: public FrepOperators {
+        public:
+            Frep* shape1;
+            Frep* shape2;
+
+            Add(){};
+            Add(Frep* shape1_,Frep* shape2_) {
+                shape1=shape1_;
+                shape2=shape2_;
+            };
+
+            float eval(){
+                //change x and y
+                shape1->setXY(X,Y);
+                shape2->setXY(X,Y);
+                return std::max(shape1->eval(),shape2->eval());
+            }
+
+            Add get(){
+                return *this;
+            }
+
+            
+    };
+
+    class Subtract: public FrepOperators {
+        public:
+            Frep* shape1;
+            Frep* shape2;
+
+            Subtract(){};
+            Subtract(Frep* shape1_,Frep* shape2_) {
+                shape1=shape1_;
+                shape2=shape2_;
+            };
+
+            float eval(){
+                //change x and y
+                shape1->setXY(X,Y);
+                shape2->setXY(X,Y);
+                return std::min(shape1->eval(),-shape2->eval());
+            }
+
+            Subtract get(){
+                return *this;
+            }
+
+            
+    };
+
+    class Intersect: public FrepOperators {
+        public:
+            Frep* shape1;
+            Frep* shape2;
+
+            Intersect(){};
+            Intersect(Frep* shape1_,Frep* shape2_) {
+                shape1=shape1_;
+                shape2=shape2_;
+            };
+
+            float eval(){
+                //change x and y
+                shape1->setXY(X,Y);
+                shape2->setXY(X,Y);
+                return std::min(shape1->eval(),shape2->eval());
+            }
+
+            Intersect get(){
+                return *this;
+            }
+            
+    };
+
+    class Clearance: public FrepOperators {
+        public:
+            Frep* shape1;
+            Frep* shape2;
+            scalar offset;
+
+            Clearance(){};
+            Clearance(Frep* shape1_,Frep* shape2_,scalar offset_) {
+                shape1=shape1_;
+                shape2=shape2_;
+                offset=offset_;
+            };
+
+            float eval(){
+                //change x and y
+                shape1->setXY(X,Y);
+                shape2->setXY(X,Y);
+                return std::min(shape1->eval(),-(offset+shape2->eval()));
+            }
+
+            Clearance get(){
+                return *this;
+            }
+    };
+
+    class Blend: public FrepOperators {
+        public:
+            Frep* shape1;
+            Frep* shape2;
+            scalar offset;
+
+            Blend(){};
+            Blend(Frep* shape1_,Frep* shape2_,scalar offset_) {
+                shape1=shape1_;
+                shape2=shape2_;
+                offset=offset_;
+            };
+
+            Blend get(){
+                return *this;
+            }
+
+            float eval(){
+                //change x and y
+                shape1->setXY(X,Y);
+                shape2->setXY(X,Y);
+                scalar fn =std::max(shape1->eval(),shape2->eval());
+                return std::max(fn,(shape1->eval())+(shape2->eval())+(offset));
+            }
+
+    };
+
+    class Morph: public FrepOperators {
+        public:
+            Frep* shape1;
+            Frep* shape2;
+            scalar weight;
+
+            Morph(){};
+            Morph(Frep* shape1_,Frep* shape2_,scalar weight_) {
+                shape1=shape1_;
+                shape2=shape2_;
+                weight=weight_;
+            };
+
+            Morph get(){
+                return *this;
+            }
+
+            float eval(){
+                //change x and y
+                shape1->setXY(X,Y);
+                shape2->setXY(X,Y);
+                return (weight*shape1->eval()+(1-weight)*shape2->eval());
+            }
+
+    };
+    //..................................................................//
+
+    //.......................TRANSFORMS.................................//
+    class FrepTransforms: public Frep {
+        public:
+            
+            //Frep(std::string const& name): name_(std::move(name)) {}
+            FrepTransforms(){};  
+    };
+
+    class Translate: public FrepTransforms {
+        public:
+            Frep* shape;
+            scalar dx;
+            scalar dy;
+
+            Translate():dx(0.0f),dy(0.0f){};
+            Translate(Frep* shape_,scalar dx_,scalar dy_) {
+                shape= shape_;
+                dx=dx_;
+                dy=dy_;
+            };
+
+            float eval(){
+                //change x and y
+                shape->setXY(X-dx,Y-dy);
+                return shape->eval();
+            }
+
+            Translate get(){
+                return *this;
+            }
+    };
+
+    class Scale: public FrepTransforms {
+        public:
+            Frep* shape;
+            scalar sx;
+            scalar sy;
+            scalar ox;
+            scalar oy;
+
+            Scale(){};
+            Scale(Frep* shape_,scalar ox_, scalar oy_,scalar sx_,scalar sy_) {
+                shape=shape_;
+                sx=sx_;
+                sy=sy_;
+                ox=ox_;
+                oy=oy_;
+            };
+
+            float eval(){
+                //change x and y
+                shape->setXY(((ox)+(X-(ox))/(sx)),((oy)+(Y-(oy))/(sy)));
+                return shape->eval();
+            }
+
+            Scale get(){
+                return *this;
+            }
+    };
+
+    class Rotate: public FrepTransforms {
+        public:
+            Frep* shape;
+            scalar angle;
+            scalar ox;
+            scalar oy;
+            const float  PI_F=3.14159265358979f;
+
+            Rotate(){};
+            Rotate(Frep* shape_,scalar ox_, scalar oy_, scalar angle_) {
+                shape=shape_;
+                angle=angle_;
+                ox=ox_;
+                oy=oy_;
+            };
+
+            Rotate get(){
+                return *this;
+            }
+
+            float eval(){
+                //change x and y
+                
+                scalar r =angle*PI_F/180.0f;
+                scalar newX=((ox)+(X-(ox))*cos(r)+(Y-(oy))*sin(r));
+                scalar newY=((oy)-(X-(ox))*sin(r)+(Y-(oy))*cos(r));
+                shape->setXY(newX,newY);
+                return shape->eval();
+            }
+
+            
+    };
+
+    class LinearArrayX: public FrepTransforms {
+        public:
+            Frep* shape;
+            scalar spacing;
+            scalar number;
+            scalar size;
+
+            const float  PI_F=3.14159265358979f;
+
+            LinearArrayX(){};
+            LinearArrayX(Frep* shape_,scalar spacing_,scalar number_,scalar size_) {
+                shape=shape_;
+                spacing=spacing_;
+                number=number_;
+                size=size_;
+            };
+
+            LinearArrayX get(){
+                return *this;
+            }
+
+            float eval(){
+                //change x and y
+                scalar xmin=-size;
+                scalar xmax=size;
+
+                scalar newX=remainderf(  X ,spacing+xmax-xmin );
+                shape->setXY(newX,Y);
+                scalar fn=shape->eval();
+
+                fn = std::min(X-(xmin),fn);
+                fn = std::min((xmin+number*(spacing+xmax-xmin))-X,fn);
+
+                return fn;
+            }
+
+            
+    };
+
+    class LinearArrayY: public FrepTransforms {
+        public:
+            Frep* shape;
+            scalar spacing;
+            scalar number;
+            scalar size;
+
+            const float  PI_F=3.14159265358979f;
+
+            LinearArrayY(){};
+            LinearArrayY(Frep* shape_,scalar spacing_,scalar number_,scalar size_) {
+                shape=shape_;
+                spacing=spacing_;
+                number=number_;
+                size=size_;
+            };
+
+            LinearArrayY get(){
+                return *this;
+            }
+
+            float eval(){
+                //change x and y
+                scalar ymin=-size;
+                scalar ymax=size;
+
+                scalar newY=remainderf(  Y ,spacing+ymax-ymin );
+                shape->setXY(X,newY);
+                scalar fn=shape->eval();
+
+                fn = std::min(Y-(ymin),fn);
+                fn = std::min((ymin+number*(spacing+ymax-ymin))-Y,fn);
+
+                return fn;
+            }
+
+            
+    };
+
+    class PolarArray: public FrepTransforms {
+        public:
+            Frep* shape;
+            scalar angle;
+            scalar radius;
+            scalar ox;
+            scalar oy;
+            const float  PI_F=3.14159265358979f;
+
+            PolarArray(){};
+            PolarArray(Frep* shape_,scalar ox_, scalar oy_, scalar angle_,scalar radius_) {
+                shape=shape_;
+                angle=angle_;
+                ox=ox_;
+                oy=oy_;
+                radius=radius_;
+            };
+
+            PolarArray get(){
+                return *this;
+            }
+
+            float eval(){
+                //change x and y
+                
+                scalar r = PI_F*angle/180.0f;
+                scalar tmin = -atan2f(-1.0f-oy,-1.0f-ox);
+                // std::cout << "tmin:" <<tmin<<'\n';
+                
+                // scalar tmax = atan2(radius+oy,radius+ox);
+
+                scalar newX=(ox)+sqrt((X-(ox))*(X-(ox))+(Y-(oy))*(Y-(oy)))*cos((tmin)+remainderf((2.0f*PI_F+atan2f(Y-(oy),X-(ox))-(tmin)),(r)));
+                scalar newY=(oy)+sqrt((X-(ox))*(X-(ox))+(Y-(oy))*(Y-(oy)))*sin((tmin)+remainderf((2.0f*PI_F+atan2f(Y-(oy),X-(ox))-(tmin)),(r)));
+
+                shape->setXY(newX,newY);
+                return shape->eval();
+            }
+
+            
+    };
+    //..................................................................//
+
+}
diff --git a/cpp_frep/main.cpp b/cpp_frep/main.cpp
new file mode 100755
index 0000000000000000000000000000000000000000..534f16b6926b95c23d5add79fd45445817dc677a
--- /dev/null
+++ b/cpp_frep/main.cpp
@@ -0,0 +1,70 @@
+#include "frep.cpp"
+#include <iostream>
+
+using namespace cpp_frep;
+
+using scalar=float;
+//--------------------------------------------------------------------------------------------------
+void drawFrep(Frep* frep,scalar minX,scalar maxX,scalar minY,scalar maxY,scalar dx,scalar dy){
+    for(scalar i=minX;i<maxX;i+=dx){
+        for(scalar j=minY;j<maxY;j+=dy){
+            frep->setXY(i,j);
+            scalar val=frep->eval();
+            if(val>0.0f){
+                std::cout << "+" ;
+
+            }else{
+                std::cout << "o" ;
+            }
+        }
+        std::cout  << '\n';
+    }
+    std::cout  << '\n';
+
+}
+
+int main(int const, char const**) {
+   
+    scalar minX=-3.0f;
+    scalar maxX=3.0f;
+    scalar minY=-5.0f;
+    scalar maxY=5.0f;
+    scalar dx=0.1f;
+    scalar dy=0.1f;
+
+    Circle c(1.0f,0.0f,0.0f);
+    drawFrep(&c,minX, maxX, minY, maxY, dx, dy);
+
+    Translate t(&c,0.0f,1.5f);
+    drawFrep(&t,minX, maxX, minY, maxY, dx, dy);
+
+    Scale sc(&t,0.0f,0.0f,2.0f,2.0f);
+    drawFrep(&sc,minX, maxX, minY, maxY, dx, dy);
+
+    Rectangle r(-1.0f,1.0f,-1.0f,1.0f);
+    drawFrep(&r,minX, maxX, minY, maxY, dx, dy);
+
+    Rotate rt(&r,0.0f,0.0f,15.0f);
+    drawFrep(&rt,minX, maxX, minY, maxY, dx, dy);
+
+    LinearArrayY la(&r,0.2f,4.0f,1.0);
+    drawFrep(&la,minX, maxX, minY, maxY, dx, dy);
+
+    Add a(&t,&r);
+    drawFrep(&a,minX, maxX, minY, maxY, dx, dy);
+
+    Subtract s(&r,&t);
+    drawFrep(&s,minX, maxX, minY, maxY, dx, dy);
+
+    Intersect i(&t,&r);
+    drawFrep(&i,minX, maxX, minY, maxY, dx, dy);
+
+    Clearance cl(&t,&r,1.0);
+    drawFrep(&cl,minX, maxX, minY, maxY, dx, dy);
+
+    PolarArray pa(&r,0.0f, 0.0f, 45.0f,-1.0f);
+    drawFrep(&pa,minX, maxX, minY, maxY, dx, dy);
+
+
+    return 0;
+}