diff --git a/stardis/radiation_field/radiation_field_solvers/base.py b/stardis/radiation_field/radiation_field_solvers/base.py index 11c4f855..6d31700d 100644 --- a/stardis/radiation_field/radiation_field_solvers/base.py +++ b/stardis/radiation_field/radiation_field_solvers/base.py @@ -468,14 +468,16 @@ 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 = np.diff(ray_z_coordinate_grid) + + 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[ ~np.isnan(ray_distance), theta_index ] = ray_distance[~np.isnan(ray_distance)]