diff --git a/geometry.py b/geometry.py
new file mode 100644
index 0000000..7fdd8cd
--- /dev/null
+++ b/geometry.py
@@ -0,0 +1,73 @@
+import numpy as np
+
+# Define radii
+r1 = 1.0
+r2 = 1.0
+
+# Define centers
+center1 = (-r1, 0)
+center2 = (r2, 0)
+
+def parameterize_circle(center, radius, num_points):
+    """
+    Generates points on the circumference of a circle.
+
+    Args:
+        center (tuple): (x, y) coordinates of the circle's center.
+        radius (float): The radius of the circle.
+        num_points (int): The number of points to generate.
+
+    Returns:
+        tuple: Two NumPy arrays (x_coords, y_coords) representing the points.
+    """
+    angles = np.linspace(0, 2 * np.pi, num_points, endpoint=False)
+    x_coords = center[0] + radius * np.cos(angles)
+    y_coords = center[1] + radius * np.sin(angles)
+    return x_coords, y_coords
+
+def incident_field(x, y):
+    """
+    Defines the incident field.
+    For now, phi_incident(x, y) = x.
+    """
+    return x
+
+def boundary_condition_on_surface(x, y, incident_field_func):
+    """
+    Calculates the desired value for the scattered field on the surface.
+    This is -incident_field_func(x, y).
+    """
+    return -incident_field_func(x, y)
+
+if __name__ == "__main__":
+    num_points_c1 = 100
+    num_points_c2 = 100
+
+    x_c1, y_c1 = parameterize_circle(center1, r1, num_points_c1)
+    x_c2, y_c2 = parameterize_circle(center2, r2, num_points_c2)
+
+    print("First 5 points of C1:")
+    for i in range(5):
+        print(f"({x_c1[i]:.2f}, {y_c1[i]:.2f})")
+
+    print("\nFirst 5 points of C2:")
+    for i in range(5):
+        print(f"({x_c2[i]:.2f}, {y_c2[i]:.2f})")
+
+    # Calculate boundary conditions for C1
+    bc_c1 = np.zeros(num_points_c1)
+    for i in range(num_points_c1):
+        bc_c1[i] = boundary_condition_on_surface(x_c1[i], y_c1[i], incident_field)
+
+    print("\nFirst 5 boundary condition values for C1:")
+    for i in range(5):
+        print(f"{bc_c1[i]:.2f}")
+
+    # Calculate boundary conditions for C2
+    bc_c2 = np.zeros(num_points_c2)
+    for i in range(num_points_c2):
+        bc_c2[i] = boundary_condition_on_surface(x_c2[i], y_c2[i], incident_field)
+
+    print("\nFirst 5 boundary condition values for C2:")
+    for i in range(5):
+        print(f"{bc_c2[i]:.2f}")
diff --git a/solver.py b/solver.py
new file mode 100644
index 0000000..dfcf236
--- /dev/null
+++ b/solver.py
@@ -0,0 +1,118 @@
+import numpy as np
+from geometry import parameterize_circle, incident_field, boundary_condition_on_surface
+from solver_setup import place_poles, construct_system_matrix
+
+# Define problem parameters (can be adjusted)
+r1_val = 1.0
+r2_val = 1.0
+center1_val = (-r1_val, 0)
+center2_val = (r2_val, 0)
+
+num_poles_c = 50  # Number of poles per circle
+pole_sf = 0.95    # Pole scale factor
+num_coll_pts_c = 50 # Number of collocation points per circle
+
+def setup_problem(r1, r2, center1, center2, num_poles_circle, pole_scale_factor, num_collocation_points_circle):
+    """
+    Generates all necessary components for the scattering problem.
+    """
+    # Generate Collocation Points for C1 and C2
+    coll_x1, coll_y1 = parameterize_circle(center1, r1, num_collocation_points_circle)
+    coll_x2, coll_y2 = parameterize_circle(center2, r2, num_collocation_points_circle)
+    all_coll_x = np.concatenate((coll_x1, coll_x2))
+    all_coll_y = np.concatenate((coll_y1, coll_y2))
+    z_collocation_pts = all_coll_x + 1j * all_coll_y
+
+    # Generate Pole Locations for C1 and C2
+    poles_x1, poles_y1 = place_poles(center1, r1, num_poles_circle, pole_scale_factor)
+    poles_x2, poles_y2 = place_poles(center2, r2, num_poles_circle, pole_scale_factor)
+    all_poles_x = np.concatenate((poles_x1, poles_x2))
+    all_poles_y = np.concatenate((poles_y1, poles_y2))
+    z_poles_pts = all_poles_x + 1j * all_poles_y
+
+    # Calculate Boundary Values (b_vector) at Collocation Points
+    b_vec = np.zeros(len(all_coll_x))
+    for i in range(len(all_coll_x)):
+        # Using boundary_condition_on_surface which is -incident_field
+        b_vec[i] = boundary_condition_on_surface(all_coll_x[i], all_coll_y[i], incident_field)
+
+    # Construct System Matrix A
+    A_mat = construct_system_matrix(z_collocation_pts, z_poles_pts)
+
+    return z_collocation_pts, z_poles_pts, A_mat, b_vec
+
+# Solve for Pole Coefficients
+# This will be done in the main block after calling setup_problem
+
+def calculate_scattered_field(x, y, z_p, coeffs):
+    """
+    Calculates the scattered field at a point (x,y) due to poles.
+    """
+    z_eval = x + 1j * y
+    scattered_field_value = 0.0
+    for k in range(len(z_p)):
+        if not np.isclose(z_eval, z_p[k]): # Avoid division by zero if poles are at eval point
+            scattered_field_value += np.real(coeffs[k] / (z_eval - z_p[k]))
+    return scattered_field_value
+
+def calculate_total_field(x, y, z_p, coeffs, inc_field_func):
+    """
+    Calculates the total field (incident + scattered) at a point (x,y).
+    """
+    scattered_val = calculate_scattered_field(x, y, z_p, coeffs)
+    incident_val = inc_field_func(x, y)
+    return scattered_val + incident_val
+
+if __name__ == "__main__":
+    print("Setting up the scattering problem...")
+    z_collocation, z_poles, A_matrix, b_vector = setup_problem(
+        r1_val, r2_val, center1_val, center2_val,
+        num_poles_c, pole_sf, num_coll_pts_c
+    )
+
+    print(f"A_matrix shape: {A_matrix.shape}")
+    print(f"b_vector shape: {b_vector.shape}")
+    print(f"z_poles shape: {z_poles.shape}")
+    print(f"z_collocation shape: {z_collocation.shape}")
+
+    print("\nSolving for pole coefficients...")
+    coefficients = np.linalg.solve(A_matrix, b_vector)
+
+    print("\nFirst 10 coefficients:")
+    for i in range(min(10, len(coefficients))):
+        print(f"Coeff {i}: {coefficients[i]:.4f}")
+
+    test_points = [
+        (2.0, 0.0),
+        (0.0, 2.0),
+        (-2.0, 0.5), # Test point near C1 but outside
+        (2.0, -0.5)  # Test point near C2 but outside
+    ]
+
+    print("\nField calculations at test points (outside circles):")
+    for pt_x, pt_y in test_points:
+        sc_field = calculate_scattered_field(pt_x, pt_y, z_poles, coefficients)
+        tot_field = calculate_total_field(pt_x, pt_y, z_poles, coefficients, incident_field)
+        print(f"Point ({pt_x:.1f}, {pt_y:.1f}): Scattered = {sc_field:.4f}, Total = {tot_field:.4f}, Incident = {incident_field(pt_x, pt_y):.4f}")
+
+    print("\nVerification: Total field at first 3 collocation points (should be ~0):")
+    for i in range(min(3, len(z_collocation))):
+        coll_x = np.real(z_collocation[i])
+        coll_y = np.imag(z_collocation[i])
+        tot_field_coll = calculate_total_field(coll_x, coll_y, z_poles, coefficients, incident_field)
+        # Expected b_vector value at this point (which is -incident_field)
+        expected_boundary_val = b_vector[i]
+        # Total field = scattered + incident. Boundary condition was scattered = -incident. So total = 0.
+        print(f"Collocation Pt ({coll_x:.2f}, {coll_y:.2f}): Total Field = {tot_field_coll:.4e} (Expected on boundary: ~0, b_vector value: {expected_boundary_val:.2f})")
+
+    print("\nVerification: Total field at a point on C1 (idx N_coll_c // 2) and C2 (idx N_coll_c + N_coll_c // 2) (should be ~0):")
+    if len(z_collocation) > num_coll_pts_c + num_coll_pts_c // 2:
+        indices_to_check = [num_coll_pts_c // 2, num_coll_pts_c + num_coll_pts_c // 2 -1] # one from each circle
+        for idx in indices_to_check:
+            coll_x = np.real(z_collocation[idx])
+            coll_y = np.imag(z_collocation[idx])
+            tot_field_coll = calculate_total_field(coll_x, coll_y, z_poles, coefficients, incident_field)
+            expected_boundary_val = b_vector[idx]
+            print(f"Collocation Pt ({coll_x:.2f}, {coll_y:.2f}): Total Field = {tot_field_coll:.4e} (Expected on boundary: ~0, b_vector value: {expected_boundary_val:.2f})")
+    else:
+        print("Not enough collocation points to run the second verification set.")
diff --git a/solver_setup.py b/solver_setup.py
new file mode 100644
index 0000000..7a9f078
--- /dev/null
+++ b/solver_setup.py
@@ -0,0 +1,113 @@
+import numpy as np
+from geometry import r1, r2, center1, center2, parameterize_circle, incident_field, boundary_condition_on_surface
+
+# Define Pole Placement Parameters
+num_poles_circle = 50  # number of poles per circle
+pole_scale_factor = 0.95  # scales poles inward from the boundary
+
+def place_poles(center, radius, num_poles, scale_factor):
+    """
+    Places poles inside a circle, scaled inward from the boundary.
+
+    Args:
+        center (tuple): (x, y) coordinates of the circle's center.
+        radius (float): The radius of the circle.
+        num_poles (int): The number of poles to generate.
+        scale_factor (float): Factor to scale poles inward from the boundary.
+
+    Returns:
+        tuple: Two NumPy arrays (pole_x_coords, pole_y_coords).
+    """
+    # Generate points on the boundary
+    boundary_x, boundary_y = parameterize_circle(center, radius, num_poles)
+
+    # Scale points inward to get pole locations
+    pole_x_coords = center[0] + (boundary_x - center[0]) * scale_factor
+    pole_y_coords = center[1] + (boundary_y - center[1]) * scale_factor
+
+    return pole_x_coords, pole_y_coords
+
+# Generate Collocation Points
+num_collocation_points_circle = 50  # can be same as num_poles for now
+coll_x1, coll_y1 = parameterize_circle(center1, r1, num_collocation_points_circle)
+coll_x2, coll_y2 = parameterize_circle(center2, r2, num_collocation_points_circle)
+
+all_coll_x = np.concatenate((coll_x1, coll_x2))
+all_coll_y = np.concatenate((coll_y1, coll_y2))
+z_collocation = all_coll_x + 1j * all_coll_y
+
+# Generate Pole Locations
+poles_x1, poles_y1 = place_poles(center1, r1, num_poles_circle, pole_scale_factor)
+poles_x2, poles_y2 = place_poles(center2, r2, num_poles_circle, pole_scale_factor)
+
+all_poles_x = np.concatenate((poles_x1, poles_x2))
+all_poles_y = np.concatenate((poles_y1, poles_y2)) # Corrected concatenation
+z_poles = all_poles_x + 1j * all_poles_y
+
+# Define System Matrix Construction Function
+def construct_system_matrix(z_coll, z_p):
+    """
+    Constructs the system matrix A for the method of fundamental solutions.
+    A[i, j] = real(1 / (z_coll[i] - z_p[j]))
+
+    Args:
+        z_coll (np.array): Complex array of collocation points.
+        z_p (np.array): Complex array of pole locations.
+
+    Returns:
+        np.array: The system matrix A.
+    """
+    num_coll_points = len(z_coll)
+    num_poles = len(z_p)
+    A_matrix = np.zeros((num_coll_points, num_poles))
+
+    for i in range(num_coll_points):
+        for j in range(num_poles):
+            A_matrix[i, j] = np.real(1 / (z_coll[i] - z_p[j]))
+    return A_matrix
+
+# Calculate Boundary Values at Collocation Points
+b_vector = np.zeros(len(all_coll_x))
+for i in range(len(all_coll_x)):
+    incident_val = incident_field(all_coll_x[i], all_coll_y[i])
+    b_vector[i] = -incident_val # Using the simplified direct relation for now
+
+if __name__ == "__main__":
+    print("Shapes of generated arrays:")
+    print(f"all_coll_x: {all_coll_x.shape}")
+    print(f"z_collocation: {z_collocation.shape}")
+    print(f"all_poles_x: {all_poles_x.shape}")
+    print(f"z_poles: {z_poles.shape}")
+    print(f"b_vector: {b_vector.shape}")
+
+    # Construct System Matrix
+    A_matrix = construct_system_matrix(z_collocation, z_poles)
+    print(f"A_matrix shape: {A_matrix.shape}")
+
+    print("\nFirst 5 pole coordinates for C1:")
+    for i in range(5):
+        print(f"({poles_x1[i]:.2f}, {poles_y1[i]:.2f})")
+
+    print("\nFirst 5 pole coordinates for C2:")
+    for i in range(5):
+        print(f"({poles_x2[i]:.2f}, {poles_y2[i]:.2f})")
+
+    print("\nFirst 5 collocation points for C1:")
+    for i in range(5):
+        print(f"({coll_x1[i]:.2f}, {coll_y1[i]:.2f})")
+
+    print("\nFirst 5 collocation points for C2:")
+    for i in range(5):
+        print(f"({coll_x2[i]:.2f}, {coll_y2[i]:.2f})")
+
+    print("\nFirst 10 values of b_vector:")
+    for i in range(10):
+        print(f"{b_vector[i]:.2f}")
+
+    print("\nFirst 5x5 submatrix of A:")
+    # Ensure we only print up to 5x5, even if matrix is smaller
+    rows_to_print = min(5, A_matrix.shape[0])
+    cols_to_print = min(5, A_matrix.shape[1])
+    for i in range(rows_to_print):
+        row_str = " ".join([f"{A_matrix[i, j]:9.4f}" for j in range(cols_to_print)])
+        print(f"[{row_str} ]")
diff --git a/visualize_field.py b/visualize_field.py
new file mode 100644
index 0000000..4ad20d0
--- /dev/null
+++ b/visualize_field.py
@@ -0,0 +1,102 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+# Import from our other modules
+from geometry import r1, r2, center1, center2, parameterize_circle, incident_field
+from solver import setup_problem, calculate_total_field
+# solver.py's setup_problem uses its own internal parameters, but we can pass
+# the ones from geometry.py to ensure consistency for this script.
+# The parameters from solver.py's setup_problem are:
+# r1_val, r2_val, center1_val, center2_val, num_poles_c, pole_sf, num_coll_pts_c
+# We will use r1, r2, center1, center2 from geometry.py
+# And define num_poles_c, pole_sf, num_coll_pts_c here for clarity for visualization
+num_poles_vis = 50
+pole_scale_factor_vis = 0.95
+num_collocation_points_vis = 50
+
+
+print("Setting up the problem for visualization...")
+# Use r1, r2, center1, center2 imported from geometry.py
+z_collocation, z_poles, A_matrix, b_vector = setup_problem(
+    r1, r2, center1, center2,
+    num_poles_vis, pole_scale_factor_vis, num_collocation_points_vis
+)
+
+print("Solving for coefficients...")
+coefficients = np.linalg.solve(A_matrix, b_vector)
+
+print("Creating grid for plotting...")
+# Define a grid range (adjust as needed)
+x_min = -2.5 * r1  # Adjusted for better visualization of two circles
+x_max = 2.5 * r1
+y_min = -2.0 * r1
+y_max = 2.0 * r1
+num_grid_points = 150  # Increased for higher resolution
+
+x_grid = np.linspace(x_min, x_max, num_grid_points)
+y_grid = np.linspace(y_min, y_max, num_grid_points)
+X, Y = np.meshgrid(x_grid, y_grid)
+
+print("Calculating total field on the grid...")
+Z_field = np.zeros_like(X)
+
+for i in range(num_grid_points):
+    if i % 10 == 0: # Simple progress indicator
+        print(f"Calculating row {i+1}/{num_grid_points}")
+    for j in range(num_grid_points):
+        current_x, current_y = X[i, j], Y[i, j]
+        
+        dist_to_c1_center_sq = (current_x - center1[0])**2 + (current_y - center1[1])**2
+        dist_to_c2_center_sq = (current_x - center2[0])**2 + (current_y - center2[1])**2
+        
+        # Using squared distance to avoid sqrt, and a small tolerance epsilon
+        epsilon_sq = 1e-4 # To avoid issues exactly on the boundary if parameterize_circle points are used
+                          # and to make the exclusion slightly larger than the circle radius for visual clarity
+
+        if dist_to_c1_center_sq < r1**2 - epsilon_sq or \
+           dist_to_c2_center_sq < r2**2 - epsilon_sq:
+            Z_field[i, j] = np.nan # Mask points inside circles
+        else:
+            Z_field[i, j] = calculate_total_field(current_x, current_y, z_poles, coefficients, incident_field)
+
+print("Plotting the field...")
+plt.figure(figsize=(12, 9)) # Adjusted figure size
+
+# Plot the field using contourf
+contour_levels = np.linspace(np.nanmin(Z_field), np.nanmax(Z_field), 50)
+plt.contourf(X, Y, Z_field, levels=contour_levels, cmap='viridis', extend='both')
+plt.colorbar(label='Total Potential ($\phi_{total}$)')
+
+# Overlay Circle Boundaries
+num_boundary_points = 100
+c1_x_boundary, c1_y_boundary = parameterize_circle(center1, r1, num_boundary_points)
+c2_x_boundary, c2_y_boundary = parameterize_circle(center2, r2, num_boundary_points)
+
+# Plot boundaries (optional, as patches will cover this)
+# plt.plot(c1_x_boundary, c1_y_boundary, color='black', linewidth=1.0)
+# plt.plot(c2_x_boundary, c2_y_boundary, color='black', linewidth=1.0)
+
+# Fill the circles with a solid color
+circle1_patch = plt.Circle(center1, r1, color='lightgray', ec='black', linewidth=1.5, zorder=10)
+circle2_patch = plt.Circle(center2, r2, color='lightgray', ec='black', linewidth=1.5, zorder=10)
+plt.gca().add_patch(circle1_patch)
+plt.gca().add_patch(circle2_patch)
+
+
+plt.title('Total Potential Field around Two Touching Circles (Incident Field $\phi_i = x$)')
+plt.xlabel('x')
+plt.ylabel('y')
+plt.axis('equal')
+plt.grid(True, linestyle='--', alpha=0.5)
+
+# Show plot
+print("Displaying plot...")
+# plt.show() # This might block in some environments, savefig is safer for automated runs
+
+# Save plot
+output_filename = 'total_field_plot.png'
+plt.savefig(output_filename)
+print(f"Plot saved to {output_filename}")
+
+# To ensure script finishes in automated environments if plt.show() blocks
+plt.close()