Added airy disc visualization to spot_diagram.py

[SeckinBerkay]

In the analysis.spot_diagram, added an airy disc to let the user analyze their system further.
Key points,

• The airy radius is calculated using the primary wavelength (like other design software to my knowledge)
• The view function is changed slightly for an optional airy disc visualization.
• If the airy disc is requested, limits are chosen according to the maximum values of the default axis_lim or calculated airy_radius.
• To avoid changing the x-y limit when the airy disc is not requested, an adjusted axis limit is chosen according to the former method.
• If the airy disc is much larger than the points in the spot diagram, the spots are not visible but the user can see the spots clearly by setting the is_airy_disc as False. An interactive plot could resolve this but arises other issues.

[HarrisonKramer]

Thanks a lot for adding this! This is a great addition.

First, there is already an issue for this feature (as you pointed out), so please assign yourself to that to avoid someone else implementing this too.

A few critical remarks:

• The Airy disk should be centered on the chief ray intersection point on the image plane. You will need to trace a chief ray first, once per field, then center the Airy disk at that location. Also note that the spot centroid is subtracted from the spots before plotting.
• The Airy disk can't be assumed to be a circle (see issue #46). It may be elliptical.
• You use the paraxial F/#, but this is not generally valid. The calculation of the radius (in each axis) will be more complex. The implementation below is one option (feel free to propose an alternative though). If the implementation becomes quite complex, I prefer to move the calculation out of the spot_diagram.py module to keep things simplified and modular.

Possible implementation for Airy radius in each axis (similar to Zemax, I think)

1. Trace a chief ray for the given field.
2. For the given field, trace rays in the pupil at top/bottom and left/right of pupil, so field coordinates (0, ±1) and (±1, 0).
3. Compute angles between chief ray and each marginal ray, define as $\theta$.
4. Compute working F/# as $N_w = 1 / (2 \cdot n \cdot \sin(\theta))$, where $n$ is the refractive index in image space (1 by default).
5. Find the average F/# for the X and Y axes (not sure this is best approach).
6. Find the radius in each axis as $1.22 \cdot N_w \cdot \lambda$
7. Plot corresponding ellipse.

Again, not positive this is the best approach, or if it generally valid for all optical systems. Feel free to make a suggestion.

Let me know if anything isn't clear. I can elaborate or provide guidance where needed.

Thanks again!
Kramer

[HarrisonKramer]

1. To make the argument more clear, I would propose to use add_airy_disk instead of is_airy_disc.
2. Please follow PEP 8 guidelines and do not put spaces around the "=" i.e., please write add_airy_disk=False. It's a small detail, but I prefer the whole codebase follows best practices.

[HarrisonKramer]

update argument name (see above)

[HarrisonKramer]

update argument name

[HarrisonKramer]

1. Follow PEP 8 and use lowercase words separated by underscores: airy_radius
2. More critically, this implementation is not generally valid because it uses the paraxial F/#. We must instead compute the F/# using real rays, so this function needs to change. Depending on how complex that function becomes, it might make sense to move that outside of this class entirely. I prefer to keep functions and classes as simple and modular as possible.

[HarrisonKramer]

update argument name

[HarrisonKramer]

move this under the if statement so it's not unnecessarily computed

[HarrisonKramer]

small note, but make this a dashed line instead to make it stand out a bit more

[SeckinBerkay]

Thanks for the elaborate review. I need to get more used to the optiland to work with traces. Will read some more documentation and code, and delve into the literature. After figuring out the F# calculation, I can fix the other issues (naming, PEP 8 convention, etc.) and resend the pull request.

[HarrisonKramer]

It's already a great start. I recognize that this feature has become more complex and is harder to tackle when less familiar with the code base. I'm happy to help where I can.

For raytracing, see Tutorial 2a.

You can also trace a single ray like this: ray = lens.trace_generic(Hx=0, Hy=0, Px=0, Py=0, wavelength=0.55), then you can get the resulting ray properties on the image plane, like x, y, z or direction cosines (L, M, N) from the ray e.g., ray.x or ray.N.

Using this functionality, you can trace chief and marginal rays, then find the angles between them (i.e., angles between two vectors, where the direction cosine defines the 3-vector). That leads you to the F/# using the approach above, I think. Please check the logic and propose another approach if you find something more suitable. I'm probably missing some edge cases where this might fail.

Thanks again! and let me know if I can clarify anything else.

Kramer

[codecov]

Codecov Report

Attention: Patch coverage is 91.95402% with 7 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
optiland/analysis/spot_diagram.py	91.95%	7 Missing ⚠️

@@            Coverage Diff             @@
##           master      #51      +/-   ##
==========================================
- Coverage   96.92%   96.84%   -0.08%     
==========================================
  Files         105      105              
  Lines        6209     6316     +107     
==========================================
+ Hits         6018     6117      +99     
- Misses        191      199       +8     

Files with missing lines	Coverage Δ
optiland/analysis/spot_diagram.py	95.93% <91.95%> (-4.07%)	⬇️

... and 3 files with indirect coverage changes

[SeckinBerkay]

Need help. I generated a code as asked.
The code,

• Generates marginal rays at Px=+-1 and P_y = +-1 (the function generate_marginal_rays)
• Generates chief rays at each the field heights at H_x, and H_y and gets the cosines. (the function generate_chief_rays_cosines)
• Calculates the angle difference between each marginal ray and the chief ray (the function angle_from_cosine)
• All these functions are merged as an algorithm inside the function airy_disc_x_y, which as a result calculates airy radius for x and y axes for each chief ray.
• This resulting algorithm is used inside the default plotting function _plot_field, and each ellipse is plotted on at each spot for the fields. Spots at the image coordinates are retrieved by the function get_image_coordinates, which basically generates chief rays, and traces them to the image plane, then gets their x and y coordinate values.

My problem is, when I test the system with the lens system EyepieceErfle, with a main wavelength of 0.4861, I get x = [11.7615535,11.7615535, 11.7615535], y = [11.7615535, 6.78006174, 7.45822964]. However, in the visualization, I don't see an ellipse. It may be due to the normalization that are made under the hood, but I have no idea. You can see the screen shot.

Here is the demo code for spot_diagram.py.

"""Spot Diagram Analysis

This module provides a spot diagram analysis for optical systems.

Kramer Harrison, 2024
"""
from copy import deepcopy
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches

plt.rcParams.update({'font.size': 12, 'font.family': 'cambria'})

class SpotDiagram:
    """Spot diagram class

    This class generates data and plots real ray intersection locations
    on the final optical surface in an optical system. These plots
    are purely geometric and give an indication of the blur produced
    by aberrations in the system.

    Attributes:
        optic (optic.Optic): instance of the optic object to be assessed
        fields (tuple): fields at which data is generated
        wavelengths (tuple[float]): wavelengths at which data is generated
        num_rings (int): number of rings in pupil distribution for ray tracing
        data (List): contains spot data in a nested list. Data is ordered as
            field (dim 0), wavelength (dim 1), then x, y and intensity data
            (dim 2).
    """

    def __init__(self, optic, fields='all', wavelengths='all', num_rings=6,
                 distribution='hexapolar'):
        """Create an instance of SpotDiagram

        Note:
            The constructor also generates all data that may later be used for
            plotting

        Args:
            optic (optic.Optic): instance of the optic object to be assessed
            fields (tuple or str): fields at which data is generated.
                If 'all' is passed, then all field points are considered.
                Default is 'all'.
            wavelengths (str or tuple[float]): wavelengths at which data is
                generated. If 'all' is passed, then all wavelengths are
                considered. Default is 'all'.
            num_rings (int): number of rings in pupil distribution for ray
                tracing. Default is 6.
            distribution (str): pupil distribution type for ray tracing.
                Default is 'hexapolar'.

        Returns:
            None
        """
        self.optic = optic
        self.fields = fields
        self.wavelengths = wavelengths
        if self.fields == 'all':
            self.fields = self.optic.fields.get_field_coords()

        if self.wavelengths == 'all':
            self.wavelengths = self.optic.wavelengths.get_wavelengths()

        self.data = self._generate_data(self.fields, self.wavelengths,
                                        num_rings, distribution)

    def view(self, figsize=(12, 4), add_airy_disk=False):
        """View the spot diagram

        Args:
            figsize (tuple): the figure size of the output window.
                Default is (12, 4).
            add_airy_disk (bool): Airy disc visualization controller.

        Returns:
            None
        """
        N = len(self.fields)
        num_rows = (N + 2) // 3

        fig, axs = plt.subplots(num_rows, 3,
                                figsize=(figsize[0], num_rows*figsize[1]),
                                sharex=True, sharey=True)
        axs = axs.flatten()

        # subtract centroid and find limits
        data = self._center_spots(deepcopy(self.data))
        geometric_size = self.geometric_spot_radius()
        axis_lim = np.max(geometric_size)

        # plot wavelengths for each field
        for k, field_data in enumerate(data):
            self._plot_field(axs[k], field_data, self.fields[k],
                             axis_lim, self.wavelengths, add_airy_disk=add_airy_disk)

        # remove empty axes
        for k in range(N, num_rows * 3):
            fig.delaxes(axs[k])

        plt.legend(bbox_to_anchor=(1.05, 0.5), loc='center left')
        plt.tight_layout()
        plt.show()

    def get_image_coordinates(self):
        """
        Generate the coordinate list on the image plane.
        """
        field_coords = self.optic.fields.get_field_coords()
        coordinate_list = []
        for idx in range(len(field_coords)):
            ray = self.optic.trace_generic(Hx = field_coords[idx][0], 
                                           Hy = field_coords[idx][1], 
                                           Px = 0, Py = 0, 
                                           wavelength=self.optic.wavelengths.primary_wavelength.value)
            coordinate_list.append([ray.x, ray.y])
        return coordinate_list
        

    def angle_from_cosine(self, cos1, cos2):
        """
        Calculate the angle (in radians and degrees) between two vectors given their cosine values.

        Parameters:
        cos1 (float): Cosine of the first vector.
        cos2 (float): Cosine of the second vector.

        Returns:
        tuple: (angle in radians, angle in degrees)
        """
        # Compute the angle using arccos
        theta_rad = np.abs(np.arccos(cos1) % (np.pi/2)  - np.arccos(cos2)% (np.pi/2))
        #theta_deg = np.degrees(theta_rad)  # Convert to degrees

        return theta_rad

    def f_number(self, n, theta):
        """
        Calculates the F#
        """
        N_w = 1 / (2 * n * np.sin(theta))
        return N_w

    def airy_radius(self, n_w, wavelength):
        """
        Calculates the airy radius
        INPUTS:
            n_w (float): F number
            wavelength (float): wavelength in micrometers.
        
        """

        return 1.22 * n_w * wavelength
    
    def generate_marginal_rays(self, wavelength):
        """
        Generates marginal rays at the edges of the stop.
        
        """
        ray_north = self.optic.trace_generic(Hx=0, Hy=0, Px=1, Py=0, wavelength=wavelength)
        ray_south = self.optic.trace_generic(Hx=0, Hy=0, Px=-1, Py=0, wavelength=wavelength)
        ray_east = self.optic.trace_generic(Hx=0, Hy=0, Px=0, Py=1, wavelength=wavelength)
        ray_west = self.optic.trace_generic(Hx=0, Hy=0, Px=0, Py=-1, wavelength=wavelength)

        return ray_north, ray_south, ray_east, ray_west

    # generate multiple chief ray's angle
    def generate_chief_rays_cosines(self, wavelength):
        """
        
        
        """
        coords = self.optic.fields.get_field_coords()
        chief_ray_cosines_list = []
        for coord in coords:
            H_x, H_y = coord[0], coord[1]
            ray_chief = self.optic.trace_generic(Hx=H_x, Hy=H_y, Px=0, Py=0, wavelength=wavelength) 
            chief_ray_cosines_x = ray_chief.L
            chief_ray_cosines_y = ray_chief.M
            chief_ray_cosines_list.append([chief_ray_cosines_x, chief_ray_cosines_y])
            
        return chief_ray_cosines_list
    def airy_disc_x_y(self, wavelength):

        ray_north, ray_south, ray_east, ray_west =  self.generate_marginal_rays(wavelength)
        
        #ray_chief = generate_chief_ray(wavelength) 
        chief_ray_cosines_list = self.generate_chief_rays_cosines(wavelength)
        
        airy_rad_x_list = []
        airy_rad_y_list = []
        for chief_ray_cosines in chief_ray_cosines_list:
            #                                           # compare x axes
            rel_angle_north = abs(self.angle_from_cosine(chief_ray_cosines[0], ray_north.L))
            rel_angle_south = abs(self.angle_from_cosine(chief_ray_cosines[0], ray_south.L))
            #                                           # compare y axes
            rel_angle_east = abs(self.angle_from_cosine(chief_ray_cosines[1], ray_east.M))
            rel_angle_west = abs(self.angle_from_cosine(chief_ray_cosines[1], ray_west.M))
            
            avg_angle_x = (rel_angle_north + rel_angle_south) / 2
            avg_angle_y = (rel_angle_east + rel_angle_west) / 2
        
            N_w_x = self.f_number(n = 1, theta = avg_angle_x)
            N_w_y = self.f_number(n = 1, theta = avg_angle_y)
            
            airy_rad_x = self.airy_radius(N_w_x, wavelength)
            airy_rad_y = self.airy_radius(N_w_y, wavelength)
            airy_rad_x_list.append(airy_rad_x)
            airy_rad_y_list.append(airy_rad_y)

        return airy_rad_x_list, airy_rad_y_list
    
    def centroid(self):
        """Centroid of each spot

        Returns:
            centroid (List): centroid for each field in the data.
        """
        norm_index = self.optic.wavelengths.primary_index
        centroid = []
        for field_data in self.data:
            centroid_x = np.mean(field_data[norm_index][0])
            centroid_y = np.mean(field_data[norm_index][1])
            centroid.append((centroid_x, centroid_y))
        return centroid

    def geometric_spot_radius(self):
        """Geometric spot radius of each spot

        Returns:
            geometric_size (List): Geometric spot radius for field and
                wavelength
        """
        data = self._center_spots(deepcopy(self.data))
        geometric_size = []
        for field_data in data:
            geometric_size_field = []
            for wave_data in field_data:
                r = np.sqrt(wave_data[0]**2 + wave_data[1]**2)
                geometric_size_field.append(np.max(r))
            geometric_size.append(geometric_size_field)
        return geometric_size

    def rms_spot_radius(self):
        """Root mean square (RMS) spot radius of each spot

        Returns:
            rms (List): RMS spot radius for each field and wavelength.
        """
        data = self._center_spots(deepcopy(self.data))
        rms = []
        for field_data in data:
            rms_field = []
            for wave_data in field_data:
                r2 = wave_data[0]**2 + wave_data[1]**2
                rms_field.append(np.sqrt(np.mean(r2)))
            rms.append(rms_field)
        return rms

    def _center_spots(self, data):
        """
        Centers the spots in the given data around their respective centroids.

        Args:
            data (List): A nested list representing the data containing spots.

        Returns:
            data (List): A nested list with the spots centered around their
                centroids.
        """
        centroids = self.centroid()
        data = deepcopy(self.data)
        for i, field_data in enumerate(data):
            for wave_data in field_data:
                wave_data[0] -= centroids[i][0]
                wave_data[1] -= centroids[i][1]
        return data

    def _generate_data(self, fields, wavelengths, num_rays=100,
                       distribution='hexapolar'):
        """
        Generate spot data for the given fields and wavelengths.

        Args:
            fields (List): A list of fields.
            wavelengths (List): A list of wavelengths.
            num_rays (int, optional): The number of rays to generate.
                Defaults to 100.
            distribution (str, optional): The distribution type.
                Defaults to 'hexapolar'.

        Returns:
            data (List): A nested list of spot intersection data for each
                field and wavelength.
        """
        data = []
        for field in fields:
            field_data = []
            for wavelength in wavelengths:
                field_data.append(self._generate_field_data(field,
                                                            wavelength,
                                                            num_rays,
                                                            distribution))
            data.append(field_data)
        return data

    def _generate_field_data(self, field, wavelength, num_rays=100,
                             distribution='hexapolar'):
        """
        Generates spot data for a given field and wavelength.

        Args:
            field (tuple): Tuple containing the field coordinates in (x, y).
            wavelength (float): The wavelength of the field.
            num_rays (int, optional): The number of rays to generate.
                Defaults to 100.
            distribution (str, optional): The distribution pattern of the
                rays. Defaults to 'hexapolar'.

        Returns:
            list: A list containing the x, y, and intensity values of the
                generated spot data.
        """
        self.optic.trace(*field, wavelength, num_rays, distribution)
        x = self.optic.surface_group.x[-1, :]
        y = self.optic.surface_group.y[-1, :]
        intensity = self.optic.surface_group.intensity[-1, :]
        return [x, y, intensity]

    def _plot_field(self, ax, field_data, field, axis_lim,
                    wavelengths, buffer=1.05,  add_airy_disk=False):
        """
        Plot the field data on the given axis.

        Parameters:
            ax (matplotlib.axes.Axes): The axis to plot the field data on.
            field_data (list): List of tuples containing x, y, and intensity
                data points.
            field (tuple): Tuple containing the Hx and Hy field values.
            axis_lim (float): Limit of the x and y axis.
            wavelengths (list): List of wavelengths corresponding to the
                field data.
            buffer (float, optional): Buffer factor to extend the axis limits.
                Default is 1.05.
            add_airy_disk (bool, optional): Whether to plot the airy disc or not, for analysis purposes.

        Returns:
            None
        """
        markers = ['o', 's', '^']
        #fields = self.optic.fields.get_field_coords()
        coordinate_list = self.get_image_coordinates()  # final coordinate list on the image plane.
        if add_airy_disk:
            airy_rad_x, airy_rad_y = self.airy_disc_x_y(wavelength=self.optic.wavelengths.primary_wavelength.value)
            # Add the ellipse at the image planes.
            for k, coordinates in enumerate(coordinate_list):
                ellipse = patches.Ellipse((coordinates[0], coordinates[1]), width=airy_rad_x[k], height=airy_rad_y[k], linestyle = "--", edgecolor='black', fill=False, linewidth=1)
                ax.add_patch(ellipse)

            adjusted_axis_lim = max(axis_lim , max(airy_rad_x), max(airy_rad_y)) * buffer
        else:
            adjusted_axis_lim = axis_lim*buffer  # if airy disc is not present, prevent over buffering.

        for k, points in enumerate(field_data):
            x, y, intensity = points
            mask = intensity != 0
            ax.scatter(x[mask], y[mask], s=10,
                       label=f'{wavelengths[k]:.4f} µm',
                       marker=markers[k % 3], alpha=0.7)
            ax.axis('square')
            ax.set_xlabel('X (mm)')
            ax.set_ylabel('Y (mm)')
            ax.set_xlim((-adjusted_axis_lim, adjusted_axis_lim))
            ax.set_ylim((-adjusted_axis_lim, adjusted_axis_lim))

        ax.set_title(f'Hx: {field[0]:.3f}, Hy: {field[1]:.3f}')
        ax.grid(alpha=0.25)

[SeckinBerkay]

All right, there was probably an incompetency between my files, now I get the x and y radii values as ([array([0.41934402]), array([0.41934402]), array([0.41934402])],
[array([0.41934402]), array([0.43032986]), array([0.43032986])]), which makes sense. If there is no additional error, and everything seems logical, I can request a pull after fixing the code according to guidelines.

[HarrisonKramer]

Thanks for continuing with this! Impressive work. Will be a great addition.

Can you go ahead and push the changes so I can review the code more easily? No problem if you prefer to make a few changes before pushing though. No rush.

I also see you (mistakenly?) added a poetry.lock file to the PR. This can be removed too.

Regards,
Kramer

[SeckinBerkay]

I have no idea how poetry.lock is included, tried but can't delete it from the pull request, may I open another project and create a pull request from there?

[HarrisonKramer]

You should be able to (temporarily) delete your poetry.lock file, commit and push the deletion, then restore your file (if you wish). Maybe you need to close your editor first. I believe I can also do it, although it will be easier for you.

Thanks for pushing the changes. I will review now.

[HarrisonKramer]

This looks great!

I did not test your code, but could you confirm that it generates Airy disk radii as expected for a few of the sample lenses? It would be very nice if you added some tests for this, but it won't be blocking the PR.

If you can make the simplification for the angle calculation, then should be good to merge.

Thanks,
Kramer

[HarrisonKramer]

I think you can simplify this calculation.

The marginal ray(s) will have a direction unit vector defined by its direction cosines: a = np.array([ray_m.L, ray_m.M, ray_n.N])

Similarly, the chief ray will have a direction unit vector: b = np.array([ray_c.L, ray_c.M, ray_c.N])

Then you can find the angle via np.arccos(np.abs(np.dot(a, b))). The abs limits the angle to 0-90 degrees.

Can you implement the calculation in this way? I think it will simplify your code too.

[SeckinBerkay]

This was a pitfall I had been thinking about. As you pointed out, my calculations may not be valid for image planes bigger than \theta 90 degrees. I will implement it as stated, thanks for clarifying.

[HarrisonKramer]

remove one "airy"

[HarrisonKramer]

I think you already trace the chief rays once in get_image_coordinates. Would it be possible to only have to trace them once, instead of twice?

[SeckinBerkay]

It was a mere attempt to spot the centroid of the rays as x-y coordinates. I was planning to eliminate that function, as it is already calculated by the default centroid(self) function inside spot_diagram.py. I will try to omit that function.

[HarrisonKramer]

add comment that Airy disk is centered on the chief ray intersection point on the image plane.

[SeckinBerkay]

This looks great!

I did not test your code, but could you confirm that it generates Airy disk radii as expected for a few of the sample lenses? It would be very nice if you added some tests for this, but it won't be blocking the PR.

If you can make the simplification for the angle calculation, then should be good to merge.

Thanks, Kramer

Thanks for the reviews, I can include some snapshots but it may take one more day to get used to the testing environment and add tests to the PR. I may add them after I commit the other requested changes if it is okay. Otherwise, I can implement the tests and then create PR, but as I said it may take a day.

[HarrisonKramer]

No worries. I think if you're confident in the implementation and show a few examples of it working as intended, then it should be okay to proceed. We can add tests later too.

[SeckinBerkay]

Update:

Tested my calculations with your reviews by a Zemax design (directly imported it into python), and they both give around (x = 5.7 mm, y= 5.7 mm). However, the problem is the normalizing of the plot (Hx, Hy type), trying to figure it out. I will figure it out probably but a help would make it faster. Figures illustrate python output (x, y for two fields) and zemax snapshots for the first field showing .

Kind regards,

[HarrisonKramer]

Very nice. It looks like you're getting close.

Can you elaborate on the issue with normalization of the plot? I want to make sure I understand what you mean. I can then help more easily.

Does this have something to do with the fact that one of the four values you print is not 5.7369 mm? I assume they should all be equal for an on-axis point.

[SeckinBerkay]

Firstly, my bad. It correctly gets both x, and y as the python output figure I provided displays x_radii list, then y_radii list (1st row is x axis, 2nd row is y axis. 1st column is for the first field, 2nd column is for the second field). Hence, the 5.66 value is for the second field, which I didn't provide an image for the zemax case, not to crowd here with more figures.

Secondly, the problem is I directly input the unnormalized airy radii found to the plotting function (patches.Ellipse). But the original spot diagram plot is done on Hx & Hy. Here are the snapshots with and without airy disc.

[HarrisonKramer]

I'm still not sure that I full understand the discrepancy. It sounds like it's related to how the code is structured. If you're unable to solve it, I can also help.

I can directly modify the pull request if you want me to assist with some of the code. Just let me know if you want some help.

It's very close!

Regards,
Kramer

[SeckinBerkay]

Hey again,

I am left with a simple problem. I use the centroid function to center the airy disc for each field and the primary wavelength. However, it does not produce correct results. I also tried tracing the chief rays for each field, this also yields the same wrong result. I say wrong according to the Zemax results. It may be due to the function _center_spots, centering the fields around 0. Results in the Figures 2 and 3 seem close, but the result shown in Figure 1 diverges from the zemax's case.

I can create a PR for the current airy disc, to be fixed. Or suggestions would be nice. Here are some figures for comparison.

Figure 1: Airy disc placement comparison for the library lens system ZMXtmp6B84

Figure 2: Airy disc placement comparison for an optical system Lens 3

Figure 3: Airy disc placement comparison for an optical system Lens 1

[HarrisonKramer]

Thanks for sharing the update!

It looks really great. The shape/size of the ellipses look consistent, at least qualitatively, so that's good. I recognize the spot centering is a bit tricky. I subtract the centroid for plotting purposes. This implies the ellipse should be centered on the chief ray position minus the centroid. Perhaps that's what you're doing, but I need to see the code to confirm.

Can you push the latest changes to this PR so I can take a closer look? I can try to debug. This is very close - I'm sure we can work it out quickly.

[SeckinBerkay]

Okay, it worked. I will give my messy code a clearer look, then push the latest changes.

[HarrisonKramer]

Looks great!

Only one minor change, then should be good to merge.

[HarrisonKramer]

Looks good!

[SeckinBerkay]

Looks good!

I may move on with the tests. The real centroids, and the airy radii of x-y may be tested. Also, helper functions may be tested.

[HarrisonKramer]

All tests are passing, so looks solid. I recognize more tests could be added, but this can be done later on. Test coverage remains high.

I really appreciate the contribution! I will merge now.

Feel free to make more contributions in the future if interested :)

Thanks,
Kramer

[HarrisonKramer]

I've just added you as a contributor in the docs as well:

https://optiland.readthedocs.io/en/latest/authors.html

[SeckinBerkay]

Thanks, I will definitely keep contributing. This one was such a rollercoaster. Thanks for all the help.