diff --git a/.gitmodules b/.gitmodules
new file mode 100644
index 0000000..fbb0078
--- /dev/null
+++ b/.gitmodules
@@ -0,0 +1,3 @@
+[submodule "deps/libqasm"]
+	path = deps/libqasm
+	url = https://github.com/QE-Lab/libqasm.git
diff --git a/Makefile b/Makefile
index e65d2e4..7b68273 100644
--- a/Makefile
+++ b/Makefile
@@ -34,11 +34,19 @@ TWOQGATETIMES 		= $(BENCHMARKS_DIR)twoQGateTimes
 THREEQGATETIMES 	= $(BENCHMARKS_DIR)threeQGateTimes
 EPR 				= $(BENCHMARKS_DIR)epr
 
+# cQASM Interface $(LIBQASM_INCL:%=-I%)
+INTERFACE_DEPS 		= $(SRC_DIR)qInterface.hpp
+RQXINTERFACE_DEPS	= $(SRC_DIR)rQXInterface.hpp
+INTERFACE	 		= $(SRC_DIR)qInterface
+RQXINTERFACE 		= $(SRC_DIR)rQXInterface
+LIBQASM_INCL = deps/libqasm/src/cqasm/include deps/libqasm/src/cqasm/tree-gen/include deps/libqasm/pybuild/cbuild/src/cqasm/include deps/libqasm/pybuild/cbuild/src/cqasm/tree-gen deps/libqasm/pybuild/cbuild/src/cqasm/func-gen 
+
+
 # list of object files
-objectFiles = $(TARGET).o $(QUBITLAYER).o $(EXAMPLES).o $(SINGLEQGATETIMES).o $(TWOQGATETIMES).o $(THREEQGATETIMES).o $(EPR).o $(TESTS).o
+objectFiles = $(TARGET).o $(QUBITLAYER).o $(EXAMPLES).o $(SINGLEQGATETIMES).o $(TWOQGATETIMES).o $(THREEQGATETIMES).o $(EPR).o $(TESTS).o $(RQXINTERFACE).o
 
 #list of executables
-executables = $(TARGET) $(SINGLEQGATETIMES) $(TWOQGATETIMES) $(THREEQGATETIMES) $(EPR) $(TESTS)
+executables = $(TARGET) $(SINGLEQGATETIMES) $(TWOQGATETIMES) $(THREEQGATETIMES) $(EPR) $(TESTS) 
 
 # debug directory
 DEBUG = main.dSYM
@@ -62,13 +70,13 @@ BENCHMARKS_STRING = "Running\ benchmarks..."
 
 all: $(TARGET)
 
-$(TARGET): $(TARGET).o $(QUBITLAYER).o $(EXAMPLES).o
-	@printf "%b" "$(CYAN)$(LINK_STRING)   $(NO_COLOR)$(TARGET).o $(QUBITLAYER).o $(EXAMPLES).o  			"
-	@$(CXX) $(CXXFLAGS) -o $(TARGET) $(TARGET).o $(QUBITLAYER).o $(EXAMPLES).o
+$(TARGET): $(TARGET).o $(QUBITLAYER).o $(EXAMPLES).o $(INTERFACE).o $(RQXINTERFACE).o deps/libqasm/lib/libcqasm.so deps/libqasm/lib/libtree-lib.so
+	@printf "%b" "$(CYAN)$(LINK_STRING)   $(NO_COLOR)$(TARGET).o $(QUBITLAYER).o $(EXAMPLES).o $(RQXINTERFACE).o \n"
+	@$(CXX) $(CXXFLAGS) -o $(TARGET) $(TARGET).o $(QUBITLAYER).o $(EXAMPLES).o $(RQXINTERFACE).o $(INTERFACE).o deps/libqasm/lib/libcqasm.so deps/libqasm/lib/libtree-lib.so
 	@printf "%b" "$(GREEN)$(OK_STRING)\n"
 	@printf "%b" "$(GREEN)$(SUCCESS_STRING)$(NO_COLOR)\n";
 
-$(TARGET).o: $(TARGET).cpp $(TARGET_DEPS) $(QLAYER_DEPS)
+$(TARGET).o: $(TARGET).cpp $(TARGET_DEPS) $(QLAYER_DEPS) $(RQXINTERFACE_DEPS) $(INTERFACE_DEPS)
 	@printf "%b" "$(BLUE)$(COM_STRING) $(NO_COLOR)$(@)                             				"
 	@$(CXX) $(CXXFLAGS) -c $(TARGET).cpp -o $(TARGET).o
 	@printf "%b" "$(GREEN)$(OK_STRING)\n"
@@ -83,6 +91,16 @@ $(EXAMPLES).o: $(EXAMPLES).cpp $(TARGET_DEPS) $(QLAYER_DEPS) $(EXAMPLES_DEPS)
 	@$(CXX) $(CXXFLAGS) -c $(EXAMPLES).cpp -o $(EXAMPLES).o
 	@printf "%b" "$(GREEN)$(OK_STRING)\n"
 
+$(INTERFACE).o: $(INTERFACE).cpp $(TARGET_DEPS) $(QLAYER_DEPS) $(RQXINTERFACE_DEPS) $(INTERFACE_DEPS)
+	@printf "%b" "$(BLUE)$(COM_STRING) $(NO_COLOR)$(@)                       				"
+	@$(CXX) $(CXXFLAGS) -o $(INTERFACE).o  -c $(INTERFACE).cpp 
+	@printf "%b" "$(GREEN)$(OK_STRING)\n"
+
+$(RQXINTERFACE).o: $(RQXINTERFACE).cpp $(TARGET_DEPS) $(QLAYER_DEPS) $(RQXINTERFACE_DEPS) $(INTERFACE_DEPS)
+	@printf "%b" "$(BLUE)$(COM_STRING) $(NO_COLOR)$(@)                       				"
+	@$(CXX) $(CXXFLAGS) -lcqasm -o $(RQXINTERFACE).o -c $(RQXINTERFACE).cpp 
+	@printf "%b" "$(GREEN)$(OK_STRING)\n"
+
 benchmark: 
 	@make singleQBenchmark 
 	@make twoQBenchmark
diff --git a/deps/libqasm b/deps/libqasm
new file mode 160000
index 0000000..95d346f
--- /dev/null
+++ b/deps/libqasm
@@ -0,0 +1 @@
+Subproject commit 95d346fe79b53b40fe47609643ce559eb3f4dbf6
diff --git a/src/QubitLayer.cpp b/src/QubitLayer.cpp
index 000daea..5eb6fcf 100644
--- a/src/QubitLayer.cpp
+++ b/src/QubitLayer.cpp
@@ -1,6 +1,7 @@
 #include <complex>
 #include <iostream>
 #include <cmath>
+#include <random>
 #include "QubitLayer.hpp"
 
 QubitLayer::QubitLayer(unsigned int numQubits, qubitLayer *qL){
@@ -12,11 +13,13 @@ QubitLayer::QubitLayer(unsigned int numQubits, qubitLayer *qL){
     qEven_ = new qubitLayer[numStates];
     qOdd_ = new qubitLayer[numStates];
     // if input is provided then use that to fill the input qubit state
-    if (!(qL == nullptr))
+    if (!(qL == nullptr)){
         for (int row = 0; row < numStates; row++)
             qEven_[row] = qL[row];
-    else
+    }
+    else{
         qEven_[0] = {1,0};
+        }
 }
 
 QubitLayer::~QubitLayer(){
@@ -212,6 +215,98 @@ void QubitLayer::mcphase(int *controls, int numControls, int target){
     updateLayer();
 }
 
+// void QubitLayer::prepz(int* target, int numTargets){
+//     // Set probability of target bit being |0> to 1 and of being |1> to 0
+//     for (int i = 0; i < numTargets; i++)
+//     {
+//         parity ? qOdd_[target[i]] = -qEven_[i] : qEven_[target[i]] = -qOdd_[i];
+         
+//     }
+//     updateLayer();
+// }
+
+
+void QubitLayer::measure(int target){
+    // Calculate total probability of a bit being zero or one, when measured
+
+    // Generate random measurement result
+
+    // Set probability of all states with bit in other state to 0
+
+    // Add probability of all states with the bit in other state to the probability of the bit in the measurement sate
+
+    // |00> 0.3
+    // |01> 0.2
+    // |10> 0.4
+    // |11> 0.1
+    // measure 1st bit
+    // prob zero: 0.5
+    // prob one: 0.5
+    // -> set 1
+    // |10> 0.7
+    // |11> 0.3
+    std::cout << "Measure bit " << target << std::endl;
+    precision probZero{0};
+    precision probOne{0};
+    for (unsigned long long int i = 0; i < numStates; i++)
+    {
+        if (checkZeroState(i))
+        {
+            std::bitset<maxQubits> state = i;
+            if (state.test(target))
+            {
+                probOne += parity ? abs(qEven_[i]) * abs(qEven_[i]) : abs(qOdd_[i]) * abs(qOdd_[i]);
+            }
+            else
+            {
+                probZero += parity ? abs(qEven_[i]) * abs(qEven_[i]) : abs(qOdd_[i]) * abs(qOdd_[i]);
+            }
+        }
+    }
+    std::cout << "probOne: " << probOne << ", probZero: " << probZero << std::endl;
+    bool result = random_number_gen(probOne);
+    std::cout << "Measurement result of qubit " << target << " is " << result  << std::endl;
+
+    for (unsigned long long int i = 0; i < numStates; i++)
+    {
+        if (checkZeroState(i))
+        {
+            std::bitset<maxQubits> state = i;
+            // If state i contains the bit in the measured state, normalize the probability
+            if (state.test(target) == result)
+            {
+                state.flip(target); // Get the corresponding state where the measured bit has the opposite value
+                double probNormal = parity ? std::sqrt(abs(qEven_[i]) * abs(qEven_[i]) + abs(qEven_[state.to_ulong()]) * abs(qEven_[state.to_ulong()])) : \
+                    std::sqrt(abs(qOdd_[i]) * abs(qOdd_[i]) + abs(qOdd_[state.to_ulong()]) * abs(qOdd_[state.to_ulong()]));
+                // std::cout << "Total probability will be " << probNormal << std::endl;
+                // std::cout << "parity ? qEven_[i] : qOdd_[i] " << (parity ? qEven_[i] : qOdd_[i]) << std::endl;
+                parity ? qOdd_[i] = (qEven_[i])*probNormal/std::sqrt(abs(qEven_[i]) * abs(qEven_[i])) : qEven_[i] = (qOdd_[i])*probNormal/std::sqrt(abs(qOdd_[i]) * abs(qOdd_[i]));
+            
+                std::cout << "State " << state.flip(target) << " will have prob " << (parity ? abs(qOdd_[i]) * abs(qOdd_[i]) : abs(qEven_[i]) * abs(qEven_[i]) ) << std::endl;
+            }
+            // Else the probabilty of this state is zero
+            else
+            {
+                std::cout << "State " << state << " will have prob 0" << std::endl;
+                parity ? qOdd_[i] = 0 : qEven_[i] = 0;
+            }
+        }
+    }
+
+    updateLayer();
+}
+
+// Returns a 1 with the probability dof threshold
+int QubitLayer::random_number_gen(double threshold)
+{
+    std::random_device rd;
+    std::uniform_real_distribution<double> distribution(0.0f, 1.0f);
+    std::mt19937 engine(rd());
+    double value = distribution(engine);
+    std::cout << "[Random number generation] value: " << value << " threshold: " << threshold << std::endl;
+    return (value <  threshold) ? 1 : 0;
+}
+
 qProb QubitLayer::getMaxAmplitude(){
     std::bitset<maxQubits> state;
     qProb result;
diff --git a/src/QubitLayer.hpp b/src/QubitLayer.hpp
index c02c98e..a85a8cd 100644
--- a/src/QubitLayer.hpp
+++ b/src/QubitLayer.hpp
@@ -24,6 +24,9 @@ class QubitLayer{
         void mcnot(int *controls, int numControls, int target);
         void cphase(int control, int target);
         void mcphase(int *controls, int numControls, int target);
+        // void prepz(int *targets, int numTargets);
+        void measure(int target);
+        int random_number_gen(double threshold);
         qProb getMaxAmplitude();
         void printMeasurement();
         void printQubits();
@@ -43,4 +46,5 @@ class QubitLayer{
         bool parity = true;
 };
 
-#endif
\ No newline at end of file
+#endif
+
diff --git a/src/main.cpp b/src/main.cpp
index f344658..0586304 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -2,13 +2,19 @@
 #include <complex>
 #include <chrono>
 #include "QubitLayer.hpp"
-#include "../examples/qAlgorithms.hpp"
 
+#include "rQXInterface.hpp"
+#include "../examples/qAlgorithms.hpp"
 int main(){
-    auto start = std::chrono::steady_clock::now();
-    QubitLayer q = grover(2, 0);
-    auto stop = std::chrono::steady_clock::now();
-    auto duration = std::chrono::duration_cast<std::chrono::microseconds>(stop - start).count();
-    std::cout << "Execution time: " << duration << " µs" << std::endl;
-    q.printMeasurement();
+    // auto start = std::chrono::steady_clock::now();
+    // QubitLayer q = grover(2, 0);
+    // auto stop = std::chrono::steady_clock::now();
+    // auto duration = std::chrono::duration_cast<std::chrono::microseconds>(stop - start).count();
+    // std::cout << "Execution time: " << duration << " µs" << std::endl;
+    // q.printMeasurement();
+    rQXInterface interface;
+    interface.set_qasm("tests/bell.qasm");
+    interface.execute();
+    
+
 }
\ No newline at end of file
diff --git a/src/qInterface.cpp b/src/qInterface.cpp
new file mode 100644
index 0000000..4afda6c
--- /dev/null
+++ b/src/qInterface.cpp
@@ -0,0 +1,7 @@
+#include "qInterface.hpp"
+
+// qInterface::qInterface(){
+// };
+
+qInterface::~qInterface(){
+};
\ No newline at end of file
diff --git a/src/qInterface.hpp b/src/qInterface.hpp
new file mode 100644
index 0000000..8183336
--- /dev/null
+++ b/src/qInterface.hpp
@@ -0,0 +1,17 @@
+#ifndef INTERFACE_H
+#define INTERFACE_H
+#include <complex>
+
+class qInterface{
+    public:
+        virtual int execute() = 0;
+        virtual int set_state(int* state_vec) = 0;
+        virtual int add_circuit() = 0;
+        virtual int* get_state() = 0;
+        virtual int set_qasm(std::string path) = 0;
+        virtual ~qInterface();
+        // qInterface();
+        
+};
+
+#endif
\ No newline at end of file
diff --git a/src/rQXInterface.cpp b/src/rQXInterface.cpp
new file mode 100644
index 0000000..6aa8676
--- /dev/null
+++ b/src/rQXInterface.cpp
@@ -0,0 +1,128 @@
+#include <complex>
+#include "rQXInterface.hpp"
+#include "QubitLayer.hpp"
+#include <cqasm.hpp>
+
+rQXInterface::rQXInterface(){
+ 
+};
+
+// Calls individual gates and keeps track of qubit state and measurement results etc.
+int rQXInterface::execute()
+{
+    QubitLayer q(_numQubits);
+    // DEBUG
+    // std::cout << "printQubits in execute()" << std::endl;
+    // q.printQubits();
+    for (auto subcircuit : _program->subcircuits)
+    {
+        // subcircuit name
+        // subcircuit iterations
+        std::cout << subcircuit->name << "\n";
+
+        execute_circuit(subcircuit, q);
+    }
+
+    q.printMeasurement();
+
+    return 0;
+};
+
+int rQXInterface::set_state(int *state_vec)
+{
+    return 0;
+};
+int rQXInterface::add_circuit()
+{
+    return 0;
+};
+
+int rQXInterface::execute_circuit(tree::base::One<cqasm::semantic::Subcircuit> circuit, QubitLayer qubitLayer)
+{
+    // DEBUG
+    // std::cout << "printQubits in execute_circuit()" << std::endl;
+    // qubitLayer.printQubits();
+    for (auto line : circuit->bundles)
+    {
+        // operation size: parallel or single
+        // operation name
+        // operation parameters: qubits (as_qubit_refs(), size,
+        
+        if (line.size() == 1) // a single operation
+        {
+            auto operation = line->items[0];
+            std::string name = operation->name;
+            std::cout << "gate: " << name << "\n";
+            // std::cout << "operands: " << operation->operands << "\n";
+
+            // auto parameter_0 = operation->operands[0]->as_qubit_refs();
+            // // auto parameter_0 = operation->operands[0];
+            // std::cout << "parameter_0: " << parameter_0->index[0]->value << "\n";
+
+            if (name == "hadamard" || name == "h")
+            {
+                auto parameter_0 = operation->operands[0]->as_qubit_refs();
+                qubitLayer.hadamard(parameter_0->index[0]->value);
+            } else if (name == "identity" || name == "i")
+            { 
+
+            } else if (name == "x"){
+                auto parameter_0 = operation->operands[0]->as_qubit_refs();
+                qubitLayer.pauliX(parameter_0->index[0]->value);
+            } else if (name == "y"){
+                auto parameter_0 = operation->operands[0]->as_qubit_refs();
+                qubitLayer.pauliY(parameter_0->index[0]->value);
+            } else if (name == "z"){
+                auto parameter_0 = operation->operands[0]->as_qubit_refs();
+                qubitLayer.pauliZ(parameter_0->index[0]->value);
+            } else if (name == "rx"){
+                auto parameter_0 = operation->operands[0]->as_qubit_refs();
+                auto parameter_1 = operation->operands[1]->as_const_real();
+                qubitLayer.rx(parameter_0->index[0]->value, parameter_1->value);
+            } else if (name == "ry"){
+                auto parameter_0 = operation->operands[0]->as_qubit_refs();
+                auto parameter_1 = operation->operands[1]->as_const_real();
+                qubitLayer.ry(parameter_0->index[0]->value, parameter_1->value);
+            } else if (name == "rz"){
+                auto parameter_0 = operation->operands[0]->as_qubit_refs();
+                auto parameter_1 = operation->operands[1]->as_const_real();
+                qubitLayer.rz(parameter_0->index[0]->value, parameter_1->value);
+            } else if (name == "cnot"){
+                auto parameter_0 = operation->operands[0]->as_qubit_refs();
+                auto parameter_1 = operation->operands[1]->as_qubit_refs();
+                qubitLayer.rx(parameter_0->index[0]->value, parameter_1->index[0]->value);
+            } else if (name == "toffoli"){
+                auto parameter_0 = operation->operands[0]->as_qubit_refs();
+                auto parameter_1 = operation->operands[1]->as_qubit_refs();
+                auto parameter_2 = operation->operands[2]->as_qubit_refs();
+                qubitLayer.toffoli(parameter_0->index[0]->value, parameter_1->index[0]->value, parameter_2->index[0]->value);
+            } else if (name == "measure"){
+                auto parameter_0 = operation->operands[0]->as_qubit_refs();
+                qubitLayer.measure(parameter_0->index[0]->value);
+            } else {
+                std::cout << "gate name: " << name << " unknown" << std::endl;
+            }
+            // DEBUG
+            qubitLayer.printQubits();
+        };
+    };
+    return 0;
+};
+
+int *rQXInterface::get_state()
+{
+    return 0;
+};
+int rQXInterface::set_qasm(std::string path)
+{
+    auto result = cqasm::analyze(path, "1.0");
+    std::cout << "Analysis successful"
+              << "\n";
+    _numQubits = result->num_qubits;
+    std::cout << "Number of qubits:" << _numQubits << "\n";
+
+    _program = result;
+    return 0;
+};
+
+rQXInterface::~rQXInterface(){};
\ No newline at end of file
diff --git a/src/rQXInterface.hpp b/src/rQXInterface.hpp
new file mode 100644
index 0000000..d60d529
--- /dev/null
+++ b/src/rQXInterface.hpp
@@ -0,0 +1,25 @@
+
+#ifndef RQXINTERFACE_H
+#define RQXINTERFACE_H
+#include "qInterface.hpp"
+#include "QubitLayer.hpp"
+#include <cqasm.hpp>
+
+class rQXInterface : public qInterface {
+    public:
+        int execute();
+        int set_state(int* state_vec);
+        int add_circuit();
+        int execute_circuit(tree::base::One<cqasm::semantic::Subcircuit> circuit, QubitLayer qubitLayer);
+        int* get_state();
+        int set_qasm(std::string path);
+        rQXInterface();
+        virtual ~rQXInterface();
+    private:
+        
+        int _numQubits;
+        cqasm::tree::One<cqasm::semantic::Program> _program;
+        
+};
+
+#endif
\ No newline at end of file
diff --git a/tests/bell.qasm b/tests/bell.qasm
new file mode 100644
index 0000000..0fb805c
--- /dev/null
+++ b/tests/bell.qasm
@@ -0,0 +1,18 @@
+version 1.0
+#define a quantum register of 2 qubits
+qubits 2
+
+#Prep
+#{prep_z q[0] | prep_z q[1]}
+
+# create a Bell pair via a Hadamard rotation
+#h q[0]
+x q[0]
+#y q[0]
+cnot q[0],q[1]
+#rx q[1], 3.787868
+
+# measure both qubits to test correlations
+measure q[0]
+measure q[1]
+