diff --git a/.gitignore b/.gitignore
index bd2d311..e330d48 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,3 +1,4 @@
 .vscode/
 .DS_Store
-__pycache__/
\ No newline at end of file
+__pycache__/
+sim_results/
\ No newline at end of file
diff --git a/MAGLEV_DIGITALTWIN_PYTHON/topSimulate.py b/MAGLEV_DIGITALTWIN_PYTHON/topSimulate.py
index a6894bf..0a0b3f4 100644
--- a/MAGLEV_DIGITALTWIN_PYTHON/topSimulate.py
+++ b/MAGLEV_DIGITALTWIN_PYTHON/topSimulate.py
@@ -10,11 +10,16 @@ from parameters import QuadParams, Constants
 from utils import euler2dcm, fmag2
 from simulate import simulate_maglev_control
 from visualize import visualize_quad
+import os
 
 
 def main():
     """Main simulation script"""
     
+    # Create output directory if it doesn't exist
+    output_dir = 'sim_results'
+    os.makedirs(output_dir, exist_ok=True)
+    
     # SET REFERENCE HERE
     # Total simulation time, in seconds
     Tsim = 2
@@ -127,7 +132,7 @@ def main():
         'eMat': eMat,
         'plotFrequency': 20,
         'makeGifFlag': True,
-        'gifFileName': 'testGif.gif',
+        'gifFileName': 'sim_results/testGif.gif',
         'bounds': [-1, 1, -1, 1, -300e-3, 0.000]
     }
     visualize_quad(S2)
@@ -161,7 +166,7 @@ def main():
     plt.grid(True)
     
     plt.tight_layout()
-    plt.savefig('simulation_results.png', dpi=150)
+    plt.savefig('sim_results/simulation_results.png', dpi=150)
     print("Saved plot to simulation_results.png")
     
     # Commented out horizontal position plot (as in MATLAB)
@@ -189,8 +194,10 @@ def main():
     plt.title("Complex Magnitude of fft Spectrum")
     plt.xlabel("f (Hz)")
     plt.ylabel("|fft(X)|")
+    # Exclude DC component (first element) from scaling to see AC components better
+    plt.ylim([0, np.max(np.abs(Y[1:])) * 1.05])  # Dynamic y-axis excluding DC
     plt.grid(True)
-    plt.savefig('fft_spectrum.png', dpi=150)
+    plt.savefig('sim_results/fft_spectrum.png', dpi=150)
     print("Saved FFT plot to fft_spectrum.png")
     
     # Show all plots
diff --git a/MAGLEV_DIGITALTWIN_PYTHON/visualize.py b/MAGLEV_DIGITALTWIN_PYTHON/visualize.py
index f4554c8..7123c5c 100644
--- a/MAGLEV_DIGITALTWIN_PYTHON/visualize.py
+++ b/MAGLEV_DIGITALTWIN_PYTHON/visualize.py
@@ -101,6 +101,9 @@ def visualize_quad(S):
         # Plot body axes
         plot_axes(ax, r, RIB)
         
+        # Add steel plate
+        plot_steel_plate(ax, S['bounds'])
+        
         ax.set_xlim([S['bounds'][0], S['bounds'][1]])
         ax.set_ylim([S['bounds'][2], S['bounds'][3]])
         ax.set_zlim([S['bounds'][4], S['bounds'][5]])
@@ -152,6 +155,10 @@ def visualize_quad(S):
             ax.set_zlabel('Z')
             ax.view_init(elev=0, azim=0)
             ax.grid(True)
+            
+            # Add steel plate
+            plot_steel_plate(ax, S['bounds'])
+            
             return []
         
         def update(frame):
@@ -219,6 +226,68 @@ def plot_body(ax, r, RIB, body_faces):
     return poly
 
 
+def plot_steel_plate(ax, bounds):
+    """
+    Plot a steel plate at z=0 with 0.25 inch thickness
+    
+    Parameters
+    ----------
+    ax : Axes3D
+        The 3D axis to plot on
+    bounds : list
+        6-element list [xmin, xmax, ymin, ymax, zmin, zmax]
+    
+    Returns
+    -------
+    plate : Poly3DCollection
+        The plate collection object
+    """
+    plate_thickness = 0.25 * 0.0254  # 0.25 inches to meters
+    x_bounds = [bounds[0], bounds[1]]
+    y_bounds = [bounds[2], bounds[3]]
+    
+    # Create plate vertices (top at z=0, bottom at z=-thickness)
+    plate_vertices = [
+        # Top surface
+        [[x_bounds[0], y_bounds[0], plate_thickness],
+         [x_bounds[1], y_bounds[0], plate_thickness],
+         [x_bounds[1], y_bounds[1], plate_thickness],
+         [x_bounds[0], y_bounds[1], plate_thickness]],
+        # Bottom surface
+        [[x_bounds[0], y_bounds[0], 0],
+         [x_bounds[1], y_bounds[0], 0],
+         [x_bounds[1], y_bounds[1], 0],
+         [x_bounds[0], y_bounds[1], 0]],
+        # Side 1
+        [[x_bounds[0], y_bounds[0], plate_thickness],
+         [x_bounds[1], y_bounds[0], plate_thickness],
+         [x_bounds[1], y_bounds[0], 0],
+         [x_bounds[0], y_bounds[0], 0]],
+        # Side 2
+        [[x_bounds[1], y_bounds[0], plate_thickness],
+         [x_bounds[1], y_bounds[1], plate_thickness],
+         [x_bounds[1], y_bounds[1], 0],
+         [x_bounds[1], y_bounds[0], 0]],
+        # Side 3
+        [[x_bounds[1], y_bounds[1], plate_thickness],
+         [x_bounds[0], y_bounds[1], plate_thickness],
+         [x_bounds[0], y_bounds[1], 0],
+         [x_bounds[1], y_bounds[1], 0]],
+        # Side 4
+        [[x_bounds[0], y_bounds[1], plate_thickness],
+         [x_bounds[0], y_bounds[0], plate_thickness],
+         [x_bounds[0], y_bounds[0], 0],
+         [x_bounds[0], y_bounds[1], 0]]
+    ]
+    
+    # Steel gray color
+    steel_color = [0.7, 0.7, 0.75]
+    plate = Poly3DCollection(plate_vertices, alpha=0.3, facecolor=steel_color, 
+                            edgecolor='darkgray', linewidths=0.5)
+    ax.add_collection3d(plate)
+    return plate
+
+
 def plot_axes(ax, r, RIB):
     """Plot body axes"""
     bodyX = 0.5 * RIB @ np.array([1, 0, 0])