diff --git a/Blender/Python Files/code b/Blender/Python Files/code
deleted file mode 100644
index 8b13789..0000000
--- a/Blender/Python Files/code	
+++ /dev/null
@@ -1 +0,0 @@
-
diff --git a/Blender/rod_sim.py b/Blender/rod_sim.py
deleted file mode 100644
index 010a91e..0000000
--- a/Blender/rod_sim.py
+++ /dev/null
@@ -1,151 +0,0 @@
-# Helper fn
-# Call in main function; Replace vals in euler rotation and location for cylinders
-# POTENTIALLY USE CYLINDER.SCALE[2] for length stretch on z-axis
-
-# Spheres[i].location and spheres[i+1].location = pos1, pos2 pass in
-
-
-import bpy
-import numpy as np
-
-# Can access values one at a time if you return as an array
-
-
-def calc_cyl_orientation(pos1, pos2):
-    pos1 = np.array(pos1)
-    pos2 = np.array(pos2)
-    depth = np.linalg.norm(pos2 - pos1)
-    dz = pos2[2] - pos1[2]
-    dy = pos2[1] - pos1[1]
-    dx = pos2[0] - pos1[0]
-    center = (pos1 + pos2) / 2
-    # Spherical coords (phi and theta); Can be used in cyl. euler rotations (Look at wikipedia diagram)
-    phi = np.arctan2(dy, dx)
-    theta = np.arccos(dz / depth)
-    angles = np.array([phi, theta])
-    # Should it return something?
-    return depth, center, angles
-
-
-def f(x):
-    k = 5
-    y, v = x[0], x[1]
-    return np.array([v, (-1) * (k) * (y)])
-
-
-# Clear existing mesh objects in the scene
-bpy.ops.object.select_all(action="DESELECT")
-bpy.ops.object.select_by_type(type="MESH")
-bpy.ops.object.delete()
-
-
-# Create spheres with different initial velocities
-
-v0_values = [20, 25, 30, 35, 40]
-spheres = (
-    []
-)  # sphere is now a list a tuples containing the sphere object and its z-velocity
-for i, v0 in enumerate(v0_values):
-    bpy.ops.mesh.primitive_uv_sphere_add(radius=0.2, location=(i * 2, 0, 0))
-    sphere = bpy.context.active_object
-    spheres.append([sphere, v0])
-
-# Create cylinders to connect the spheres
-cylinders = []
-for i in range(len(spheres) - 1):
-    bpy.ops.mesh.primitive_cylinder_add(radius=0.1, depth=1)
-    cylinder = bpy.context.active_object
-
-    # connect cylinders to spheres on creation
-    depth, center, angles = calc_cyl_orientation(
-        spheres[i][0].location, spheres[i + 1][0].location
-    )
-    cylinder.location = (center[0], center[1], center[2])
-    cylinder.rotation_euler = (0, angles[1], angles[0])
-    cylinder.scale[2] = depth
-    cylinders.append(cylinder)
-
-
-####### SIMULATION ########
-
-
-# Creating Simulation Parameters
-
-# dt = 10**(-6)
-dt = 10 ** (-3)
-framerate = 25
-simulation_ratio = int(1 / framerate / dt)
-time = np.arange(0, 10, dt)
-
-# Simulate the behavior of all objects over "time"
-for time_index, t in enumerate(time[:-1]):
-    for i, [sphere, vz] in enumerate(spheres):
-        # update 3D-world position of each sphere
-        x = np.array([sphere.location.z, vz])
-        x = x + f(x) * dt
-        sphere.location.z = x[0]
-        sphere.location.y = x[0]
-        spheres[i][1] = x[1]
-
-    if (
-        time_index % simulation_ratio
-    ) == 0:  # this is an index which we want to write to a keyframe
-        # then, we add updated sphere locations to the keyframe
-        for i, (sphere, vz) in enumerate(spheres):
-            sphere.keyframe_insert(
-                data_path="location",
-                frame=int(time_index / simulation_ratio) + 1,
-            )
-
-        # now we update cylinder orientation and then draw those to the keyframe.
-        for i in range(len(cylinders)):
-            depth, center, angles = calc_cyl_orientation(
-                spheres[i][0].location, spheres[i + 1][0].location
-            )
-            cylinders[i].location = (center[0], center[1], center[2])
-            cylinders[i].rotation_euler = (0, angles[1], angles[0])
-            cylinders[i].scale[2] = depth
-
-            # Keyframe the cylinder's location , rotation, and scaling
-            cylinders[i].keyframe_insert(
-                data_path="location",
-                frame=int(time_index / simulation_ratio) + 1,
-            )
-            cylinders[i].keyframe_insert(
-                data_path="rotation_euler",
-                frame=int(time_index / simulation_ratio) + 1,
-            )
-            cylinders[i].keyframe_insert(
-                data_path="scale", frame=int(time_index / simulation_ratio) + 1
-            )
-
-
-"""
-for time_index, t in enumerate(time[:-1]):
-    for i, (sphere, v0) in enumerate(spheres):
-        pos = [0]
-        x = np.array([0, v0])
-        x = x + f(x) * dt
-        pos.append(x[0])
-        if (time_index % simulation_ratio) == 0: # this is an index which we want to write to a keyframe
-            sphere.location.z = pos[time_index]
-            sphere.keyframe_insert(data_path="location", index=2, frame=int(time_index/simulation_ratio) + 1)
-            sphere.keyframe_insert(data_path="location", index=1, frame=int(time_index/simulation_ratio) + 1)
-            sphere.keyframe_insert(data_path="location", index=0, frame=int(time_index/simulation_ratio) + 1)
-
-            # Check if the current index is within the range of cylinders
-            if i < len(cylinders):
-                # Update the cylinder's location and rotation to connect the current sphere to the next one
-                # xloc = sphere.location.x + (spheres[i+1][0].location.x - sphere.location.x) / 2
-                depth, center, angles = calc_cyl_orientation(sphere.location, spheres[i+1][0].location)
-                cylinders[i].location = (center[0], center[1], center[2])
-                cylinders[i].rotation_euler = (0,angles[1], angles[0])
-                cylinders[i].scale[2] = depth
-
-
-                # Keyframe the cylinder's location and rotation
-                cylinders[i].keyframe_insert(data_path="location", index=2, frame=int(time_index/simulation_ratio) + 1)
-                cylinders[i].keyframe_insert(data_path="location", index=0, frame=int(time_index/simulation_ratio) + 1)
-                cylinders[i].keyframe_insert(data_path="rotation_euler", index=2, frame=int(time_index/simulation_ratio) + 1)
-
-"""
diff --git a/examples/single_rigid_rod_spring_action.py b/examples/single_rigid_rod_spring_action.py
new file mode 100644
index 0000000..0539e79
--- /dev/null
+++ b/examples/single_rigid_rod_spring_action.py
@@ -0,0 +1,55 @@
+import numpy as np
+
+import bsr
+
+bsr.clear_mesh_objects()
+
+# Environment configuration
+k = 5  # spring constant
+
+
+def f(x):
+    y, v = x[0], x[1]
+    return np.array([v, (-1) * (k) * (y)])
+
+
+def euler_forward_step(x, f, dt):
+    return x + f(x) * dt
+
+
+# Rod configuration
+N = 5
+velocity = np.arange(20, 45, N)  # initial z-velocities
+position = np.stack(
+    [
+        np.zeros(N),  # x
+        np.zeros(N),  # y
+        np.linspace(0, 10, N),  # z
+    ],
+    axis=0,
+)  # (3, 5)
+radius = 0.2 * np.ones(N)  # radius of the rod
+
+rod = bsr.rod()
+rod.update(keyframe=0, positions=position, radius=radius)
+
+
+####### SIMULATION ########
+# Simulation parameters
+
+dt = 10 ** (-3)
+framerate = 25
+simulation_ratio = int(1 / framerate / dt)
+time = np.arange(0, 10, dt)
+
+# Euler-Forward time stepper
+for time_index, t in enumerate(time[:-1]):
+    state = np.stack([position[2], velocity], axis=0)
+    position[2] = euler_forward_step(state, f, dt)
+
+    if time_index % simulation_ratio == 0:
+        # update the rod
+        keyframe = int(time_index / simulation_ratio) + 1
+        rod.update(keyframe=keyframe, positions=position, radius=radius)
+
+bsr.save("single_rigid_rod_spring_action.blend")
diff --git a/load_write_blender.py b/load_write_blender.py
deleted file mode 100644
index 7594e89..0000000
--- a/load_write_blender.py
+++ /dev/null
@@ -1,174 +0,0 @@
-import os
-import numpy as np
-import bpy
-
-import bsr
-
-class PendulumBlender:
-    def __init__(self, location, ball_radius):
-        self.location = location
-        self.ball_radius = ball_radius
-
-        origin_location = np.array([0, 0, 0])
-        # Add a sphere
-        self.add_sphere(
-            location=self.location, 
-            radius=self.ball_radius
-        )
-        # Add a cylinder
-        center, angle, length = PendulumBlender.calculate_cylinder_pose(
-            loc1=origin_location, 
-            loc2=self.location
-        )
-        self.add_cylinder(
-            location=center, 
-            rotation=angle, 
-            depth=length, 
-            radius=0.1*self.ball_radius
-        )
-
-    def __repr__(self):
-        return f"Pendulum at {self.location} with ball radius {self.ball_radius}"
-    
-    def add_sphere(self, location, radius):
-        self.sphere = bpy.ops.mesh.primitive_uv_sphere_add(
-            radius=radius, 
-            location=location
-        )
-        self.sphere = bpy.context.active_object
-
-    def add_cylinder(self, location, rotation, depth, radius):
-        bpy.ops.mesh.primitive_cylinder_add(
-            location=location,
-            depth=depth,
-            radius=radius,
-        )
-        self.cylinder = bpy.context.active_object
-        self.cylinder.rotation_euler = rotation
-
-    @staticmethod
-    def calculate_cylinder_pose(
-        loc1: np.ndarray, 
-        loc2: np.ndarray
-    ) -> tuple[float, np.ndarray, np.ndarray]:
-        # Calculate the depth of the cylinder
-        length = np.linalg.norm(loc1 - loc2)
-        # Calculate the center of the cylinder
-        center = (loc1 + loc2) / 2
-        # Calculate the difference in z, y, and x
-        delta = loc2 - loc1
-        # Calculate the angles of the cylinder
-        angle = np.array([
-            0, 
-            np.arccos(delta[2] / length), 
-            np.arctan2(delta[1], delta[0])
-        ])
-        return center, angle, length
-    
-    def update(self, position):
-        self.sphere.location = position
-        center, angle, length = PendulumBlender.calculate_cylinder_pose(
-            loc1=np.array([0, 0, 0]), 
-            loc2=position
-        )
-        self.cylinder.location = center
-        self.cylinder.rotation_euler = angle
-        self.cylinder.scale[2] = length
-
-class Pendulum:
-    def __init__(self, length, euler_angles):
-        self.length = length
-        self.euler_angles = euler_angles
-        self.position = self.calculate_position()
-        self.velocity = np.array([0., 0., 0.])
-    
-    def calculate_position(self):
-        x = self.length * np.cos(self.euler_angles[1]) * np.cos(self.euler_angles[0])
-        y = self.length * np.cos(self.euler_angles[1]) * np.sin(self.euler_angles[0])
-        z = self.length * np.sin(self.euler_angles[1])
-        return np.array([x, y, z])
-    
-    def get_position(self):
-        return self.position
-    
-    def update(self, dt):
-        # Update the velocity
-        self.position += self.velocity * dt/2
-        self.velocity += np.array([0, 0, -9.81]) * dt
-        self.position += self.velocity * dt/2
-
-
-def delete_all():
-    bpy.ops.object.select_all(action='SELECT')
-    bpy.ops.object.delete()
-
-
-def main():
-    delete_all()
-    
-    pendulum_length = 0.3
-    pendulum_euler_angles = np.array([0., 0.])
-    pendulum = Pendulum(
-        length=pendulum_length, 
-        euler_angles=pendulum_euler_angles
-    )
-    pendulum_blender = PendulumBlender(
-        location=pendulum.get_position(), 
-        ball_radius=0.2
-    )
-    
-    dt = 10**(-3)
-    framerate = 25
-    simulation_ratio = int(1 / framerate / dt)
-    time = np.arange(0, 10, dt)
-
-    for k, t in enumerate(time):
-        pendulum.update(dt)
-        if k % simulation_ratio == 0:
-            # Update the location of the pendulum
-            pendulum_blender.update(position=pendulum.get_position())
-            
-
-            # # Update the scene
-            bpy.context.view_layer.update()
-    
-    ### Saving the file ###
-    write_filepath = "pendulum.blend"
-    write_filepath = os.path.join(os.path.dirname(os.path.realpath(__file__)), write_filepath)
-    #Saves above script as .blend file; stores on personal device using filepath
-    bpy.ops.wm.save_as_mainfile(filepath=write_filepath)  
-    with bpy.data.libraries.load(load_filepath, link=True) as (
-        data_from,
-        data_to,
-    ):
-        data_to.objects = [
-            name for name in data_from.objects if name.startswith("S")
-        ]
-
-
-    # load_filepath = "Blender/3rodbr2.blend"
-
-    # with bpy.data.libraries.load(load_filepath, link=True) as (data_from, data_to):
-    #     data_to.objects = [name for name in data_from.objects if name.startswith("S")]
-
-    # for obj in data_to.objects:
-    #     assert obj is not None
-    #     print(obj.name)
-
-    # # TODO: The writing part is not working
-    # # data_to.objects["Cube"].select_set(True)
-    # write_filepath = "/Blender/3rodbr2_write.blend"
-    # bpy.data.libraries.write(write_filepath, set(bpy.context.selected_objects), path_remap="RELATIVE")
-
-    
-    
-    # TODO: The writing part is not working
-    # data_to.objects["Cube"].select_set(True)
-    write_filepath = "/Blender/3rodbr2_write.blend"
-    bpy.data.libraries.write(
-        write_filepath, set(bpy.context.selected_objects), path_remap="RELATIVE"
-    )
-
-
-if __name__ == "__main__":
-    main()
diff --git a/practice.py b/practice.py
deleted file mode 100644
index 9ced419..0000000
--- a/practice.py
+++ /dev/null
@@ -1,22 +0,0 @@
-import bpy
-
-
-class SimpleOperator(bpy.types.Operator):
-    bl_idname = "object.simple_operator"
-    bl_label = "Tool Name"
-
-    def execute(self, context):
-        print("Hello World")
-        return {"FINISHED"}
-
-
-def register():
-    bpy.utils.register_class(SimpleOperator)
-
-
-def unregister():
-    bpy.utils.unregister_class(SimpleOperator)
-
-
-if __name__ == "__main__":
-    register()
diff --git a/src/bsr/__init__.py b/src/bsr/__init__.py
index e773e9e..b3b919e 100644
--- a/src/bsr/__init__.py
+++ b/src/bsr/__init__.py
@@ -3,6 +3,7 @@
 
 from .collections import *
 from .rod import *
+from .macros import *
 
 
 def get_version() -> str:
diff --git a/src/bsr/macros.py b/src/bsr/macros.py
new file mode 100644
index 0000000..4730364
--- /dev/null
+++ b/src/bsr/macros.py
@@ -0,0 +1,10 @@
+__all__ = ["clear_mesh_objects"]
+
+import bpy
+
+
+def clear_mesh_objects() -> None:
+    # Clear existing mesh objects in the scene
+    bpy.ops.object.select_all(action="DESELECT")
+    bpy.ops.object.select_by_type(type="MESH")
+    bpy.ops.object.delete()