diff --git a/src/geometry/nerf.py b/src/geometry/nerf.py
new file mode 100644
index 0000000..ee53fd2
--- /dev/null
+++ b/src/geometry/nerf.py
@@ -0,0 +1,50 @@
+import numpy as np
+
+# "Practical conversion from torsion space to Cartesian space for in silico protein synthesis"
+#     https://pubmed.ncbi.nlm.nih.gov/15898109/
+# Natural Extension Reference Frame algorithm
+def nerf(a, b, c, bond_length, bond_angle_deg, torsion_angle_deg):
+    """
+    Parameters:
+    - a, b, c: (3,) numpy arrays of xyz coordinates of the previous atoms.
+    - bond_length: distance between c and d (Angstroms).
+    - bond_angle_deg: angle b-c-d (Degrees).
+    - torsion_angle_deg: dihedral angle a-b-c-d (Degrees).
+    
+    Returns:
+    - d: (3,) numpy array of coordinates for the new atom.
+    """
+    # convert to radians
+    bond_angle_rad = np.radians(bond_angle_deg)
+    torsion_angle_rad = np.radians(torsion_angle_deg)
+    
+    # firstly calculate the position of D in the local reference frame
+    # the bond b-c is the local x axis
+   
+
+    bond_angle_complement = np.pi - bond_angle_rad
+    cd_local = np.array([
+        bond_length * np.cos(bond_angle_complement),
+        bond_length * np.sin(bond_angle_complement) * np.cos(torsion_angle_rad),
+        bond_length * np.sin(bond_angle_complement) * np.sin(torsion_angle_rad)
+    ])
+    
+    # next few steps are converting from local frame to global frame 
+    bc = c - b
+    bc_u = bc / np.linalg.norm(bc)
+    ab = b - a
+    
+    # normal to the plane a-b-c is the local z axis
+    n = np.cross(ab, bc_u)
+    n_u = n / np.linalg.norm(n)
+    
+    # the local y axis is perpendicular to both the bond bc and the normal
+    m_u = np.cross(n_u, bc_u)
+    
+    # rotation matrix columns: [x_axis, y_axis, z_axis] (3,3)
+    rotation_matrix = np.column_stack((bc_u, m_u, n_u))
+    
+    # Global_Pos = Origin(C) + (Rotation * Local_Pos)
+    d_global = c + np.dot(rotation_matrix, cd_local)
+    
+    return d_global
diff --git a/src/legacy/backbone.py b/src/legacy/backbone.py
index 22c1102..470e15f 100644
--- a/src/legacy/backbone.py
+++ b/src/legacy/backbone.py
@@ -17,7 +17,7 @@ GEO = {
     'C_N_H_angle': 119.0, 
 }
 
-# Legacy function checking if two atoms (implementation as CAs) are within threshold.
+# Legacy function checking if two atoms are within threshold.
 def check_clashes(new_coord, existing_coords, threshold=3.0):
     if len(existing_coords) == 0: return False
     diff = existing_coords - new_coord
diff --git a/src/legacy/nerf.py b/src/legacy/nerf.py
deleted file mode 100644
index f4f3a0c..0000000
--- a/src/legacy/nerf.py
+++ /dev/null
@@ -1,77 +0,0 @@
-import numpy as np
-
-def place_atom_nerf(a, b, c, bond_length, bond_angle_deg, torsion_angle_deg):
-    """
-    Calculates the coordinates of atom D given atoms A, B, C and internal coords.
-    
-    Parameters:
-    - a, b, c: (3,) numpy arrays of xyz coordinates of the previous atoms.
-    - bond_length: distance between c and d (Angstroms).
-    - bond_angle_deg: angle b-c-d (Degrees).
-    - torsion_angle_deg: dihedral angle a-b-c-d (Degrees).
-    
-    Returns:
-    - d: (3,) numpy array of coordinates for the new atom.
-    """
-    # 1. Convert to rad 
-    bond_angle_rad = np.radians(bond_angle_deg)
-    torsion_angle_rad = np.radians(torsion_angle_deg)
-    
-    # 2. Calculate the position of D in the local reference frame
-    # We align the local X-axis with the bond B->C.
-    # Note: The geometric calculation uses the complement of the bond angle 
-    #  so we use pi - bond_angle.
-    
-    # D_local represents the vector C->D in the local frame
-    d_local = np.array([
-        bond_length * np.cos(np.pi - bond_angle_rad),
-        bond_length * np.sin(np.pi - bond_angle_rad) * np.cos(torsion_angle_rad),
-        bond_length * np.sin(np.pi - bond_angle_rad) * np.sin(torsion_angle_rad)
-    ])
-    
-    # 3. Build the transformation matrix from Local to Global Frame
-    
-    # Vector from B to C (normalized) - this is our local X-axis
-    bc = c - b
-    bc_u = bc / np.linalg.norm(bc)
-    
-    # Vector from A to B
-    ab = b - a
-    
-    # Normal to the plane A-B-C (normalized) - this is our local Z-axis
-    # We use Cross Product to find the perpendicular vector
-    n = np.cross(ab, bc_u)
-    n_u = n / np.linalg.norm(n)
-    
-    # The "Up" vector in the plane (normalized) - this is our local Y-axis
-    # It is perpendicular to both the bond BC and the normal N
-    m_u = np.cross(n_u, bc_u)
-    
-    # Create the Rotation Matrix columns [x_axis, y_axis, z_axis]
-    # Shape: (3, 3)
-    rotation_matrix = np.column_stack((bc_u, m_u, n_u))
-    
-    # 4. Transform D_local to Global Coordinates
-    # Global_Pos = Origin(C) + (Rotation * Local_Pos)
-    d_global = c + np.dot(rotation_matrix, d_local)
-    
-    return d_global
-
-# --- Simple Test Case ---
-if __name__ == "__main__":
-    # Define 3 arbitrary starting atoms (A, B, C)
-    atom_a = np.array([0.0, 1.0, 0.0])
-    atom_b = np.array([0.0, 0.0, 0.0]) # Origin
-    atom_c = np.array([1.5, 0.0, 0.0]) # Along X-axis
-    
-    # Define internal coordinates for the next atom (D)
-    length = 1.5       # Angstroms
-    b_angle = 120.0    # Degrees
-    torsion = 90.0     # Degrees (Should point "out" of the screen/plane)
-
-    new_atom = place_atom_nerf(atom_a, atom_b, atom_c, length, b_angle, torsion)
-    
-    print(f"Atom A: {atom_a}")
-    print(f"Atom B: {atom_b}")
-    print(f"Atom C: {atom_c}")
-    print(f"Atom D (Calculated): {np.round(new_atom, 3)}")
diff --git a/src/legacy/ramachandran.py b/src/legacy/ramachandran.py
deleted file mode 100644
index 58655c3..0000000
--- a/src/legacy/ramachandran.py
+++ /dev/null
@@ -1,74 +0,0 @@
-import numpy as np
-
-class RamachandranSampler:
-    def __init__(self):
-        
-        # Generalised for 18 amino acids. Weighted for beta strand region
-        self.general_dist = [
-            # beta strand (top left)
-            {'phi_m': -120, 'phi_s': 20, 'psi_m': 140, 'psi_s': 20, 'w': 0.70},
-            # right handed helix (bottom left)
-            {'phi_m': -65,  'phi_s': 15, 'psi_m': -40, 'psi_s': 15, 'w': 0.25},
-            # left handed helix (top right)
-            {'phi_m': 60,   'phi_s': 15, 'psi_m': 40,  'psi_s': 15, 'w': 0.05}
-        ]
-
-        # Glycine more flexible
-        self.glycine_dist = [
-            # beta strand (top left)
-            {'phi_m': -100, 'phi_s': 30, 'psi_m': 140, 'psi_s': 30, 'w': 0.3},
-            # right handed helix (bottom left)
-            {'phi_m': -60,  'phi_s': 30, 'psi_m': -30, 'psi_s': 30, 'w': 0.2},
-            # left handed helix (top right)
-            {'phi_m': 60,   'phi_s': 30, 'psi_m': 30,  'psi_s': 30, 'w': 0.2},
-            # bottom right for glycine
-            {'phi_m': 100,  'phi_s': 30, 'psi_m': -140,'psi_s': 30, 'w': 0.3}
-        ]
-
-        # Proline no left handed helix
-        self.proline_dist = [
-            # beta strand (top left)
-            {'phi_m': -63,  'phi_s': 5,  'psi_m': 150, 'psi_s': 20, 'w': 0.8},
-            # right handed helix (bottom left)
-            {'phi_m': -63,  'phi_s': 5,  'psi_m': -35, 'psi_s': 20, 'w': 0.2}
-        ]
-
-    # returns the correct distribution
-    def _get_distribution(self, res_name):
-        if res_name == 'GLY':
-            return self.glycine_dist
-        elif res_name == 'PRO':
-            return self.proline_dist
-        else:
-            return self.general_dist
-
-    # returns (phi, psi) for the residue (units degrees)
-    def sample(self, res_name): 
-        # the possible regions to sample 
-        dist_options = self._get_distribution(res_name)
-        
-        weights = [d['w'] for d in dist_options]
-        # normalise the weights to ensure sum is exactly 1.0 (no float errors)
-        weights = np.array(weights) / np.sum(weights)
-        # choose a region based on the weights 
-        choice_idx = np.random.choice(len(dist_options), p=weights)
-        selected = dist_options[choice_idx]
-        
-        # do a gaussian sample to get some variation
-        phi = np.random.normal(selected['phi_m'], selected['phi_s'])
-        psi = np.random.normal(selected['psi_m'], selected['psi_s'])
-        
-        # wrap angles to +-180 degrees
-        phi = (phi + 180) % 360 - 180
-        psi = (psi + 180) % 360 - 180
-        
-        return phi, psi
-
-# example
-if __name__ == "__main__":
-    sampler = RamachandranSampler()
-    
-    print("Sampling 5 residues...")
-    for res in ['ALA', 'GLY', 'PRO', 'TRP', 'VAL']:
-        phi, psi = sampler.sample(res)
-        print(f"{res}: Phi={phi:.1f}, Psi={psi:.1f}")
diff --git a/src/main.py b/src/main.py
index 209c53c..f39665b 100644
--- a/src/main.py
+++ b/src/main.py
@@ -1,15 +1,4 @@
-from backbone.ramachandran import RamachandranSampler
-from rotamers.dunbrack import DunbrackRotamerLibrary
-
-rs = RamachandranSampler()
-rl = DunbrackRotamerLibrary()
-
-res_name = "PHE"
-
-print(f"Backbone torsion angles for {res_name}")
-phi, psi = rs.sample(res_name)
-print(f"phi: {phi}, psi: {psi}")
-
-print(f"Sidechain rotamers for {res_name}")
-for rotamer in rl.rotamer_params(res_name, 100, 100):
-    print(rotamer.p, rotamer.chis)
+# Get primary sequence i.e. "PHE", "SER", "THR", "PRO", "VAL"...
+# Build first residue including sidechain
+# Add ghost C' for first NeRF iteration
+# Using C 
diff --git a/src/test.py b/src/test.py
new file mode 100644
index 0000000..6c7edd3
--- /dev/null
+++ b/src/test.py
@@ -0,0 +1,17 @@
+from backbone.ramachandran import RamachandranSampler
+from rotamers.dunbrack import DunbrackRotamerLibrary
+
+rs = RamachandranSampler()
+rl = DunbrackRotamerLibrary()
+
+res_name = "PHE"
+
+print(f"Backbone torsion angles for {res_name}:")
+phi, psi = rs.sample(res_name)
+print(f"phi: {phi}, psi: {psi}")
+print(" ")
+
+print(f"Sidechain rotamers for {res_name} (top 5):")
+for rotamer in rl.rotamer_params(res_name, 100, 100)[:5]:
+    print(rotamer.p, rotamer.chis)
+print(" ")