diff --git a/stardis/radiation_field/base.py b/stardis/radiation_field/base.py
index 6464b13d..22de6b33 100644
--- a/stardis/radiation_field/base.py
+++ b/stardis/radiation_field/base.py
@@ -43,9 +43,18 @@ def __init__(self, frequencies, source_function, stellar_model, num_of_thetas):
         self.opacities = Opacities(frequencies, stellar_model)
         self.F_nu = np.zeros((stellar_model.no_of_depth_points, len(frequencies)))
 
-        thetas, weights = np.polynomial.legendre.leggauss(num_of_thetas)
-        self.thetas = (thetas / 2) + 0.5 * np.pi / 2
-        self.I_nus_weights = weights * np.pi / 2
+        # thetas, weights = np.polynomial.legendre.leggauss(num_of_thetas) # This block migrates to Gauss-Legendre quadrature
+        # self.thetas = (thetas / 2) + 0.5 * np.pi / 2
+        # self.I_nus_weights = weights * np.pi / 2
+
+        dtheta = (np.pi / 2) / num_of_thetas  # Korg uses Gauss-Legendre quadrature here
+        start_theta = dtheta / 2
+        end_theta = (np.pi / 2) - (dtheta / 2)
+        self.thetas = np.linspace(start_theta, end_theta, num_of_thetas)
+        self.I_nus_weights = (
+            2 * np.pi * dtheta * np.sin(self.thetas) * np.cos(self.thetas)
+        )
+
         self.I_nus = np.zeros(
             (stellar_model.no_of_depth_points, len(frequencies), len(self.thetas))
         )
diff --git a/stardis/radiation_field/radiation_field_solvers/base.py b/stardis/radiation_field/radiation_field_solvers/base.py
index 6d31700d..ef15dbf6 100644
--- a/stardis/radiation_field/radiation_field_solvers/base.py
+++ b/stardis/radiation_field/radiation_field_solvers/base.py
@@ -468,14 +468,14 @@ def calculate_spherical_ray(thetas, depth_points_radii):
     ray_distance_through_layer_by_impact_parameter = np.zeros(
         (len(depth_points_radii) - 1, len(thetas))
     )
-    
+
     dr = np.diff(depth_points_radii)
     for theta_index, theta in enumerate(thetas):
         b = depth_points_radii[-1] * np.sin(theta)  # impact parameter of the ray
         ray_z_coordinate_grid = np.sqrt(
             depth_points_radii**2 - b**2
         )  # Rays that don't go deeper than a layer will have a nan here
-        
+
         ray_distance = dr * depth_points_radii[1:] / ray_z_coordinate_grid[1:]
         # ray_distance = np.diff(ray_z_coordinate_grid)
         ray_distance_through_layer_by_impact_parameter[