diff --git a/code_gen/cgen.py b/code_gen/cgen.py
index f1a913b649402b2791a47686b384875ebab8ee7b..227f15cb932376eec2521a885a86685084083461 100644
--- a/code_gen/cgen.py
+++ b/code_gen/cgen.py
@@ -10,7 +10,7 @@ from ir.lit import Lit
 from ir.loops import For, Iter, ParticleFor, While
 from ir.math import Sqrt
 from ir.memory import Malloc, Realloc
-from ir.properties import Property, PropertyList
+from ir.properties import Property, PropertyList, RegisterProperty
 from ir.select import Select
 from ir.sizeof import Sizeof
 from ir.utils import Print
@@ -42,8 +42,14 @@ class CGen:
         self.print("#include <stdio.h>")
         self.print("#include <stdlib.h>")
         self.print("#include <stdbool.h>")
+        self.print("//---")
+        self.print("#include \"runtime/pairs.hpp\"")
+        self.print("#include \"runtime/read_from_file.hpp\"")
+        self.print("")
+        self.print("using namespace pairs;")
         self.print("")
         self.print("int main() {")
+        self.print("    PairsSim *ps = new PairsSim();")
         self.generate_statement(ast_node)
         self.print("}")
         self.print.end()
@@ -146,6 +152,16 @@ class CGen:
             array_name = ast_node.array.name()
             self.print(f"{array_name} = ({tkw} *) realloc({array_name}, {size});")
 
+        if isinstance(ast_node, RegisterProperty):
+            p = ast_node.property()
+            ptype = "Prop_Integer"  if p.type() == Type_Int else \
+                    "Prop_Float"    if p.type() == Type_Float else \
+                    "Prop_Vector"   if p.type() == Type_Vector else \
+                    "Prop_Invalid"
+
+            assert ptype != "Prop_Invalid", "Invalid property type!"
+            self.print(f"ps->addProperty(Property({p.id()}, \"{p.name()}\", {p.name()}, {ptype}));")
+
         if isinstance(ast_node, Timestep):
             self.generate_statement(ast_node.block)
 
@@ -239,7 +255,7 @@ class CGen:
         if isinstance(ast_node, PropertyList):
             tid = CGen.temp_id
             list_ref = f"prop_list_{tid}"
-            list_def = ", ".join(str(p.id) for p in ast_node)
+            list_def = ", ".join(str(p.id()) for p in ast_node)
             self.print(f"const int {list_ref}[] = {{{list_def}}};")
             CGen.temp_id += 1
             return list_ref
diff --git a/ir/properties.py b/ir/properties.py
index 3b5635c4a1bb9bc3cb97bdfe4444077409b2b07e..17029a69b73ae796130f28ce831c3c7f24ff47ca 100644
--- a/ir/properties.py
+++ b/ir/properties.py
@@ -38,13 +38,16 @@ class Properties:
 
         return None
 
+    def __iter__(self):
+        yield from self.props
+
 
 class Property(ASTNode):
     last_prop_id = 0
 
     def __init__(self, sim, name, dtype, default, volatile, layout=Layout_AoS):
         super().__init__(sim)
-        self.id = Property.last_prop_id
+        self.prop_id = Property.last_prop_id
         self.prop_name = name
         self.prop_type = dtype
         self.prop_layout = layout
@@ -56,7 +59,7 @@ class Property(ASTNode):
         return f"Property<{self.prop_name}>"
 
     def id(self):
-        return self.id
+        return self.prop_id
 
     def name(self):
         return self.prop_name
@@ -94,3 +97,16 @@ class PropertyList(ASTNode):
 
     def length(self):
         return len(self.list)
+
+
+class RegisterProperty(ASTNode):
+    def __init__(self, sim, prop):
+        super().__init__(sim)
+        self.prop = prop
+        self.sim.add_statement(self)
+
+    def property(self):
+        return self.prop
+
+    def __str__(self):
+        return f"Property<{self.prop.name()}>"
diff --git a/pairs.py b/pairs.py
index dab286c5177b26f3d468ebc9da6254341246e2b8..d816433b91e4b1ca5260685c9bd7c971eac7b0b3 100644
--- a/pairs.py
+++ b/pairs.py
@@ -3,4 +3,4 @@ from sim.particle_simulation import ParticleSimulation
 
 
 def simulation(ref, dims=3, timesteps=100):
-    return ParticleSimulation(CGen(f"{ref}.c"), dims, timesteps)
+    return ParticleSimulation(CGen(f"{ref}.cpp"), dims, timesteps)
diff --git a/runtime/pairs.hpp b/runtime/pairs.hpp
index 57c346cb65ab24d5e9e62c891a29298245f721c6..802eac1dc602d94591f7b8b4d4bfe4548e762a93 100644
--- a/runtime/pairs.hpp
+++ b/runtime/pairs.hpp
@@ -1,6 +1,10 @@
+#include <algorithm>
 #include <iostream>
+#include <memory>
 #include <vector>
 
+#pragma once
+
 #define PAIRS_ASSERT(a)
 #define PAIRS_ERROR(a)
 
@@ -14,45 +18,35 @@ enum PropertyType {
     Prop_Integer = 0,
     Prop_Float,
     Prop_Vector
-}
+};
 
 enum DataLayout {
     Invalid = -1,
     AoS = 0,
     SoA
-}
+};
 
-class PairsSim {
-    PairsSim() : {}
+template<typename real_t>
+class Vector3 {
+private:
+    real_t x, y, z;
 
 public:
-    void addProperty(Property prop) { properties.push_back(prop); }
-
-    Property *getProperty(property_t id) {
-        auto p = std::find_if(properties.begin(), properties.end(), [id](Property p) {
-            return p.getId() == id;
-        });
-
-        return (p != std::end(properties)) ? p : nullptr;
-    }
-
-    Property *getPropertyByName(std::string name) {
-        auto p = std::find_if(properties.begin(), properties.end(), [name](Property p) {
-            return p.getName() == name;
-        });
-
-        return (p != std::end(properties)) ? p : nullptr;
-    }
-
-    VectorPtr getVectorProperty(property_t property) { return static_cast<VectorPtr>(getProperty(property)->getPointer()); }
-    int *getIntegerPropertyPtr(property_t property) { return static_cast<int *>(getProperty(property)->getPointer()); }
-    double *getFloatPropertyPtr(property_t property) { return static_cast<double *>(getProperty(property)->getPointer()); }
-
-private:
-    std::vector<Property> properties;
+    Vector3() : x((real_t) 0.0), y((real_t) 0.0), z((real_t) 0.0) {}
+    Vector3(real_t v) : x(v), y(v), z(v) {}
+    Vector3(real_t x_, real_t y_, real_t z_) : x(x_), y(y_), z(z_) {}
 };
 
 class Property {
+protected:
+    property_t id;
+    std::string name;
+    void *ptr;
+    PropertyType type;
+    layout_t layout;
+    size_t sx, sy;
+
+public:
     Property(property_t id_, std::string name_, void *ptr_, PropertyType type_) :
         id(id_),
         name(name_),
@@ -60,59 +54,80 @@ class Property {
         type(type_),
         layout(Invalid) {}
 
-    Property(property_t id_, std::string name_, void *ptr_, PropertyType type_, layout_t layout_, size_t sx_ size_t sy) :
+    Property(property_t id_, std::string name_, void *ptr_, PropertyType type_, layout_t layout_, size_t sx_, size_t sy_) :
         id(id_),
         name(name_),
+        ptr(ptr_),
         type(type_),
         layout(layout_),
         sx(sx_),
         sy(sy_) {
 
         PAIRS_ASSERT(type != Prop_Invalid && layout_ != Invalid && sx_ > 0 && sy_ > 0);
-        ptr = static_cast<void *>(new VectorPtr(ptr_, layout_, sx_, sy_));
     }
 
-public:
     property_t getId() { return id; }
     std::string getName() { return name; }
-    void *getPointer() { return ptr.get(); }
-    void setPointer(void *ptr_) { ptr = (layout == Invalid) ? ptr_ : static_cast<void *>(new VectorPtr(ptr_, layout, sx, sy)); }
+    void *getPointer() { return ptr; }
+    void setPointer(void *ptr_) { ptr = ptr_; }
+    PropertyType getType() { return type; }
+    layout_t getLayout() { return layout; }
+};
 
-private:
-    property_t id;
-    std::string name;
-    std::shared_ptr<void *> ptr;
-    PropertyType type;
-    layout_t layout;
-    size_t sx, sy;
+class IntProperty : public Property {
+public:
+    inline int &operator()(int i) { return static_cast<int *>(ptr)[i]; }
 };
 
-class VectorPtr {
-    VectorPtr(void *ptr_, layout_t layout_, size_t sx_, size_t sy_) : ptr(ptr_), layout(layout_), sx(sx_), sy(sy_) {}
+class FloatProperty : public Property {
+public:
+    inline double &operator()(int i) { return static_cast<double *>(ptr)[i]; }
+};
 
+class VectorProperty : public Property {
 public:
-    double &operator[](int i, int j) {
-        double *dptr = static_cast<double *>(ptr.get());
+    double &operator()(int i, int j) {
+        PAIRS_ASSERT(type != Prop_Invalid && layout_ != Invalid && sx_ > 0 && sy_ > 0);
+        double *dptr = static_cast<double *>(ptr);
         if(layout == AoS) { return dptr[i * sy + j]; }
         if(layout == SoA) { return dptr[j * sx + i]; }
-        PAIRS_ERROR("VectorPtr::operator[]: Invalid data layout!");
-        return 0.0;
+        PAIRS_ERROR("VectorProperty::operator[]: Invalid data layout!");
+        return dptr[0];
     }
+};
 
+class PairsSim {
 private:
-    std::shared_ptr<void *> ptr;
-    layout_t layout;
-    size_t sx, sy;
-};
+    std::vector<Property> properties;
 
-template<typename real_t>
-class Vector3 {
-    Vector3() : x((real_t) 0.0), y((real_t) 0.0), z((real_t) 0.0) : {}
-    Vector3(real_t v) : x(v), y(v), z(v) : {}
-    Vector3(real_t x_, real_t y_, real_t z_) : x(x_), y(y_), z(z_) {}
+public:
+    PairsSim() {}
+    void addProperty(Property prop) { properties.push_back(prop); }
 
-private:
-    real_t x, y, z;
+    Property &getProperty(property_t id) {
+        auto p = std::find_if(properties.begin(), properties.end(), [id](Property p) {
+            return p.getId() == id;
+        });
+
+        PAIRS_ASSERT(p != std::end(properties));
+        return *p;
+    }
+
+    Property &getPropertyByName(std::string name) {
+        auto p = std::find_if(properties.begin(), properties.end(), [name](Property p) {
+            return p.getName() == name;
+        });
+
+        PAIRS_ASSERT(p != std::end(properties));
+        return *p;
+    }
+
+    inline IntProperty &getAsIntegerProperty(Property &prop) { return static_cast<IntProperty&>(prop); }
+    inline FloatProperty &getAsFloatProperty(Property &prop) { return static_cast<FloatProperty&>(prop); }
+    inline VectorProperty &getAsVectorProperty(Property &prop) { return static_cast<VectorProperty&>(prop); }
+    inline IntProperty &getIntegerProperty(property_t property) { return static_cast<IntProperty&>(getProperty(property)); }
+    inline FloatProperty &getFloatProperty(property_t property) { return static_cast<FloatProperty&>(getProperty(property)); }
+    inline VectorProperty &getVectorProperty(property_t property) { return static_cast<VectorProperty&>(getProperty(property)); }
 };
 
 }
diff --git a/runtime/read_from_file.hpp b/runtime/read_from_file.hpp
index d523479020495c894bd7b5825b921d1c2314ae9c..85d1c2f18a58e5d2288cb86eabde76f58a9018b5 100644
--- a/runtime/read_from_file.hpp
+++ b/runtime/read_from_file.hpp
@@ -1,54 +1,57 @@
-#include <stdio.h>
+#include <iostream>
+#include <string.h>
+#include <fstream>
+#include <sstream>
 //---
 #include "pairs.hpp"
 
+#pragma once
+
 namespace pairs {
 
-size_t read_particle_data(PairsSim *ps, const char *filename, double *grid_buffer, property_t properties[], size_t nprops) {
-    std::ifstream in_file(filename);
+size_t read_particle_data(PairsSim *ps, const char *filename, double *grid_buffer, const property_t properties[], size_t nprops) {
+    std::ifstream in_file(filename, std::ifstream::in);
     std::string line;
     size_t n = 0;
+    int read_grid_data = 0;
 
     if(in_file.is_open()) {
         while(getline(in_file, line)) {
+            std::stringstream line_stream(line);
+            std::string in0;
             int i = 0;
-            char *in0 = strtok(line, ",");
-            while(in0 != NULL) {
-                if(grid_data_unread) {
+
+            while(std::getline(line_stream, in0, ',')) {
+                if(!read_grid_data) {
                     PAIRS_ASSERT(i < ps->getNumDims() * 2);
                     grid_buffer[i] = std::stod(in0);
+                    read_grid_data = 1;
                 } else {
                     PAIRS_ASSERT(i < nprops);
-                    property_t p = properties[i];
-                    auto prop_type = ps->getProperty(p)->getType();
-
-                    switch(prop_type) {
-                        case Prop_Vector:
-                            auto vector_ptr = ps->getVectorPropertyPtr(p);
-                            char *in1 = strtok(NULL, ",");
-                            PAIRS_ASSERT(in1 != NULL);
-                            char *in2 = strtok(NULL, ",");
-                            PAIRS_ASSERT(in2 != NULL);
-                            vector_ptr[n] = Vector3(std::stod(in0), std::stod(in1), std::stod(in2));
-                            break;
-
-                        case Prop_Integer:
-                            auto int_ptr = ps->getIntegerPropertyPtr(p);
-                            int_ptr[n] = std::stoi(in0);
-                            break;
-
-                        case Prop_Float:
-                            auto float_ptr = ps->getFloatPropertyPtr(p);
-                            float_ptr[n] = std::stod(in0);
-                            break;
+                    property_t p_id = properties[i];
+                    auto prop = ps->getProperty(p_id);
+                    auto prop_type = prop.getType();
 
-                        default:
-                            fprintf(stderr, "read_particle_data(): Invalid property type!");
-                            return -1;
+                    if(prop_type == Prop_Vector) {
+                        auto vector_ptr = ps->getAsVectorProperty(prop);
+                        std::string in1, in2;
+                        std::getline(line_stream, in1, ',');
+                        std::getline(line_stream, in2, ',');
+                        vector_ptr(n, 0) = std::stod(in0);
+                        vector_ptr(n, 1) = std::stod(in1);
+                        vector_ptr(n, 2) = std::stod(in2);
+                    } else if(prop_type == Prop_Integer) {
+                        auto int_ptr = ps->getAsIntegerProperty(prop);
+                        int_ptr(n) = std::stoi(in0);
+                    } else if(prop_type == Prop_Float) {
+                        auto float_ptr = ps->getAsFloatProperty(prop);
+                        float_ptr(n) = std::stod(in0);
+                    } else {
+                        std::cerr << "read_particle_data(): Invalid property type!" << std::endl;
+                        return -1;
                     }
                 }
 
-                in0 = strtok(NULL, ",");
                 i++;
             }
 
diff --git a/sim/properties.py b/sim/properties.py
index f99d85778cad1d499e02600e880c843f299977f1..2313be162def600277a0e5a7ed2ad1cba50a60ec 100644
--- a/sim/properties.py
+++ b/sim/properties.py
@@ -1,6 +1,7 @@
 from ir.data_types import Type_Float, Type_Vector
 from ir.loops import ParticleFor
 from ir.memory import Malloc, Realloc
+from ir.properties import RegisterProperty
 from ir.utils import Print
 
 
@@ -26,6 +27,7 @@ class PropertiesAlloc:
                 Realloc(self.sim, p, sizes)
             else:
                 Malloc(self.sim, p, sizes, True)
+                RegisterProperty(self.sim, p)
 
         return self.sim.block