diff --git a/.gitignore b/.gitignore
index bdf3094..1f7461e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,4 +4,10 @@ __pycache__
 /build
 /cmake-build-*
 CMakeFiles/
-*.cmake
\ No newline at end of file
+*.cmake
+/rust/target
+/rust/target/*
+rust/Cargo.lock
+out.out
+output.json
+*.png
diff --git a/CMakeLists.txt b/CMakeLists.txt
index b936e89..fd18ac1 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -82,4 +82,4 @@ if(ENABLE_NCURSES)
     find_package(Curses REQUIRED)
     target_compile_definitions(orbital_simulator PRIVATE NCURSES_ENABLED)
     target_link_libraries(orbital_simulator PRIVATE ${CURSES_LIBRARIES})
-endif() 
\ No newline at end of file
+endif() 
diff --git a/animate.py b/animate.py
index 0b57d29..6773435 100644
--- a/animate.py
+++ b/animate.py
@@ -58,6 +58,22 @@ def read_output_file(filename):
     
     return positions, times
 
+def calculate_l2_point(sun_pos, earth_pos):
+    """Calculate L2 Lagrange point position."""
+    # Calculate distance between bodies
+    dx = earth_pos[0] - sun_pos[0]
+    dy = earth_pos[1] - sun_pos[1]
+    r = np.sqrt(dx*dx + dy*dy)
+    
+    # Calculate mass ratio (using approximate values)
+    mu = 5.972e24 / (1.989e30 + 5.972e24)  # Earth mass / (Sun mass + Earth mass)
+    
+    # Calculate L2 position (beyond Earth)
+    L2_x = sun_pos[0] + (1 + mu**(1/3)) * dx
+    L2_y = sun_pos[1] + (1 + mu**(1/3)) * dy
+    
+    return (L2_x, L2_y)
+
 def center_positions(positions, center_body):
     """Center all positions relative to the specified body."""
     if center_body not in positions:
@@ -103,9 +119,7 @@ def create_animation(positions, times, output_file=None, center_body=None):
     # Set up the plot
     ax.set_xlabel('X (m)')
     ax.set_ylabel('Y (m)')
-    title = 'Orbital Simulation'
-    if center_body:
-        title += f' (centered on {center_body})'
+    title = 'Orbital Simulation (L2 centered)'
     ax.set_title(title)
     ax.legend()
     
@@ -121,11 +135,23 @@ def create_animation(positions, times, output_file=None, center_body=None):
     if not all_x or not all_y:
         print("Error: No valid position data found")
         return
+    
+    # Calculate initial L2 point for reference
+    if 'Sun' in positions and 'Earth' in positions:
+        initial_l2 = calculate_l2_point(positions['Sun'][0], positions['Earth'][0])
+        # Set a fixed zoom level around L2 (adjust the factor to change zoom)
+        zoom_factor = 2e9  # 2 million km view
+        ax.set_xlim(-zoom_factor, zoom_factor)
+        ax.set_ylim(-zoom_factor, zoom_factor)
         
-    max_range = max(max(abs(min(all_x)), abs(max(all_x))),
-                   max(abs(min(all_y)), abs(max(all_y))))
-    ax.set_xlim(-max_range, max_range)
-    ax.set_ylim(-max_range, max_range)
+        # Create L2 point scatter plot at origin
+        l2_scatter, = ax.plot([0], [0], 'gx', markersize=10, label='L2')
+        scatters['L2'] = l2_scatter
+    else:
+        max_range = max(max(abs(min(all_x)), abs(max(all_x))),
+                       max(abs(min(all_y)), abs(max(all_y))))
+        ax.set_xlim(-max_range, max_range)
+        ax.set_ylim(-max_range, max_range)
     
     def init():
         """Initialize the animation."""
@@ -135,10 +161,20 @@ def create_animation(positions, times, output_file=None, center_body=None):
     
     def update(frame):
         """Update the animation for each frame."""
-        for name, scatter in scatters.items():
-            if positions[name]:  # Only update if we have positions
-                x, y = zip(*positions[name][:frame+1])
-                scatter.set_data(x, y)
+        if 'Sun' in positions and 'Earth' in positions:
+            # Calculate L2 point for current frame
+            l2_point = calculate_l2_point(positions['Sun'][frame], positions['Earth'][frame])
+            
+            # Update positions relative to L2
+            for name, scatter in scatters.items():
+                if name == 'L2':
+                    scatter.set_data([0], [0])  # Keep L2 at origin
+                elif positions[name]:  # Only update if we have positions
+                    x, y = zip(*positions[name][:frame+1])
+                    # Shift positions relative to L2
+                    x = [pos - l2_point[0] for pos in x]
+                    y = [pos - l2_point[1] for pos in y]
+                    scatter.set_data(x, y)
         return list(scatters.values())
     
     # Create the animation
diff --git a/config/jwst.json b/config/jwst.json
new file mode 100644
index 0000000..6741461
--- /dev/null
+++ b/config/jwst.json
@@ -0,0 +1,22 @@
+{
+    "bodies": [
+        {
+            "name": "Earth",
+            "mass": 5.972e24,
+            "position": [1.47095e11, 0, 0],
+            "velocity": [0, 29290, 0]
+        },
+        {
+            "name": "Sun",
+            "mass": 1.989e30,
+            "position": [0, 0, 0],
+            "velocity": [0, 0, 0]
+        },
+        {
+            "name": "JWST",
+            "mass": 6500,
+            "position": [149217067274.40607, 0, 0],
+            "velocity": [0, 29729.784, 0]
+        }
+    ]
+} 
\ No newline at end of file
diff --git a/config/planets.json b/config/planets.json
index 13c415e..b66ca48 100644
--- a/config/planets.json
+++ b/config/planets.json
@@ -21,9 +21,15 @@
         {
             "name": "Moon",
             "mass": 7.34767309e22,
-            "position": [149982270700, 0, 0],
-            "velocity": [0, 30822, 0]
+            "position": [146699500000, 0, 0],
+            "velocity": [0, 28278, 0]
         },
+	{
+		"name": "JWST",
+		"mass": 6500,
+		"position": [1.49995e11, 0, 0],
+		"velocity": [0,30603,0]
+	},
         {
             "name": "Sun",
             "mass": 1.989e30,
@@ -31,4 +37,4 @@
             "velocity": [0, 0, 0]
         }
     ]
-} 
\ No newline at end of file
+} 
diff --git a/python/plot_potential.py b/python/plot_potential.py
new file mode 100644
index 0000000..394cc46
--- /dev/null
+++ b/python/plot_potential.py
@@ -0,0 +1,158 @@
+import json
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+from matplotlib.colors import LogNorm
+
+# Constants from our Rust implementation
+G = 6.67430e-11  # Gravitational constant
+R_0 = 149597870700     # Earth radius (normalization length)
+M_0 = 5.972e24   # Earth mass (normalization mass)
+T_0 = 5060.0     # Characteristic time
+
+def load_bodies(config_file):
+    """Load bodies from config file."""
+    with open(config_file, 'r') as f:
+        data = json.load(f)
+    return data['bodies']
+
+def calculate_potential(x, y, bodies):
+    """Calculate gravitational potential at point (x,y)."""
+    potential = 0.0
+    for body in bodies:
+        dx = x - body['position'][0]
+        dy = y - body['position'][1]
+        r = np.sqrt(dx*dx + dy*dy)
+        if r > 0:  # Avoid division by zero
+            potential -= G * body['mass'] / r
+    return potential
+
+def calculate_lagrange_points(bodies):
+    """Calculate approximate positions of Lagrange points."""
+    if len(bodies) != 2:
+        return []
+    
+    # Find primary and secondary bodies
+    if bodies[0]['mass'] > bodies[1]['mass']:
+        primary, secondary = bodies[0], bodies[1]
+    else:
+        primary, secondary = bodies[1], bodies[0]
+    
+    # Calculate distance between bodies
+    dx = secondary['position'][0] - primary['position'][0]
+    dy = secondary['position'][1] - primary['position'][1]
+    r = np.sqrt(dx*dx + dy*dy)
+    
+    # Calculate mass ratio
+    mu = secondary['mass'] / (primary['mass'] + secondary['mass'])
+    
+    # Calculate Lagrange points
+    L1 = np.array([primary['position'][0] + (1 - mu**(1/3)) * dx,
+                   primary['position'][1] + (1 - mu**(1/3)) * dy])
+    
+    L2 = np.array([primary['position'][0] + (1 + mu**(1/3)) * dx,
+                   primary['position'][1] + (1 + mu**(1/3)) * dy])
+    
+    L3 = np.array([primary['position'][0] - (1 + 5*mu/12) * dx,
+                   primary['position'][1] - (1 + 5*mu/12) * dy])
+    
+    L4 = np.array([primary['position'][0] + 0.5 * dx - 0.866 * dy,
+                   primary['position'][1] + 0.5 * dy + 0.866 * dx])
+    
+    L5 = np.array([primary['position'][0] + 0.5 * dx + 0.866 * dy,
+                   primary['position'][1] + 0.5 * dy - 0.866 * dx])
+    
+    return [L1, L2, L3, L4, L5]
+
+def plot_potential(config_file, grid_size=200, extent=1.5):
+    """Plot gravitational potential field."""
+    # Load bodies
+    bodies = load_bodies(config_file)
+    
+    # Find Earth and Sun
+    earth = next((body for body in bodies if body['name'] == 'Earth'), None)
+    sun = next((body for body in bodies if body['name'] == 'Sun'), None)
+    
+    if not (earth and sun):
+        print("Error: Earth and Sun must be present in the config file")
+        return
+    
+    # Calculate Lagrange points for Earth-Sun system
+    lagrange_points = calculate_lagrange_points([sun, earth])
+    
+    # Create coordinate grid centered on Earth
+    x = np.linspace(earth['position'][0] - extent*R_0, earth['position'][0] + extent*R_0, grid_size)
+    y = np.linspace(earth['position'][1] - extent*R_0, earth['position'][1] + extent*R_0, grid_size)
+    X, Y = np.meshgrid(x, y)
+    
+    # Calculate potential at each grid point
+    Z = np.zeros_like(X)
+    for i in range(grid_size):
+        for j in range(grid_size):
+            Z[i,j] = calculate_potential(X[i,j], Y[i,j], bodies)
+    
+    # Take absolute value for logarithmic scale
+    Z_abs = np.abs(Z)
+    
+    # Calculate reference potential at Earth's position
+    earth_potential = abs(calculate_potential(earth['position'][0], earth['position'][1], [sun]))
+    
+    # Set scale to focus on Earth-Sun system
+    vmin = earth_potential / 10  # Show potential variations within an order of magnitude
+    vmax = earth_potential * 10
+    
+    # Create plot
+    plt.figure(figsize=(12, 10))
+    ax = plt.gca()
+    
+    # Plot potential field with logarithmic scale
+    im = ax.imshow(Z_abs, extent=[x[0], x[-1], y[0], y[-1]],
+                   origin='lower', cmap='viridis', norm=LogNorm(vmin=vmin, vmax=vmax))
+    
+    # Add colorbar
+    divider = make_axes_locatable(ax)
+    cax = divider.append_axes("right", size="5%", pad=0.05)
+    plt.colorbar(im, cax=cax, label='|Gravitational Potential| (J/kg)')
+    
+    # Plot bodies with different colors
+    body_colors = {
+        'Sun': 'yellow',
+        'Earth': 'blue',
+        'Moon': 'gray',
+        'Mars': 'red',
+        'Venus': 'orange',
+        'Mercury': 'brown',
+        'Jupiter': 'orange',
+        'Saturn': 'gold',
+        'Uranus': 'lightblue',
+        'Neptune': 'blue'
+    }
+    
+    for body in bodies:
+        color = body_colors.get(body['name'], 'white')  # Default to white if name not found
+        ax.plot(body['position'][0], body['position'][1], 'o', 
+                color=color, markersize=15, label=body['name'])
+    
+    # Plot Lagrange points
+    for i, point in enumerate(lagrange_points, 1):
+        ax.plot(point[0], point[1], 'r+', markersize=12, label=f'L{i}')
+    
+    # Customize plot
+    ax.set_xlabel('X (m)')
+    ax.set_ylabel('Y (m)')
+    ax.set_title('Gravitational Potential Field with Lagrange Points (Log Scale)')
+    ax.legend()
+    
+    # Save plot
+    plt.savefig('potential_field.png', dpi=300, bbox_inches='tight')
+    plt.show()
+
+if __name__ == "__main__":
+    import argparse
+    parser = argparse.ArgumentParser(description='Plot gravitational potential field from config file')
+    parser.add_argument('config_file', help='Path to config file')
+    parser.add_argument('--grid-size', type=int, default=200, help='Grid size for potential calculation')
+    parser.add_argument('--extent', type=float, default=1.5, help='Plot extent in Earth radii')
+    args = parser.parse_args()
+    
+    plot_potential(args.config_file, args.grid_size, args.extent) 
\ No newline at end of file