PetroFit: A Python Package for Computing Petrosian Radii and Fitting Galaxy Light Profiles

PetroFit is an open-source Python package, based on Astropy and Photutils, that can calculate Petrosian profiles and fit galaxy images. It offers end-to-end tools for making accurate photometric measurements, estimating morphological properties, and fitting 2D models to galaxy images. Petrosian metric radii can be used for model parameter estimation and aperture photometry to provide accurate total fluxes. Correction tools are provided for improving Petrosian radii estimates affected by galaxy morphology. PetroFit also provides tools for sampling Astropy-based models (including custom profiles and multi-component models) onto image grids and enables PSF convolution to account for the effects of seeing. These capabilities provide a robust means of modeling and fitting galaxy light profiles. We have made the PetroFit package publicly available on GitHub (PetroFit/petrofit) and PyPi (pip install petrofit).


INTRODUCTION
Since the advent of quantitative photometric studies of galaxies, measuring their true sizes and fluxes has proven to be a challenging task. This is because the true extent of a galaxy can be hard to define: it can vary with morphology, be difficult to distinguish from its neighbors, or difficult to measure due to low surface brightness features. Over time, researchers have proposed various approaches for estimating the fluxes, shapes, and angular (or projected) sizes of galaxies. A powerful technique for measuring accurate properties of galaxies is to combine parametric measurements with analysis of the profile of the light distribution of a galaxy.
Parametric descriptions of galaxies were first introduced by de Vaucouleurs (1948), who proposed a powerlaw profile that reasonably models the projected intensity profiles of galaxies. He used an exp{−kr 1/4 } intensity profile to describe the radial light distribution of elliptical galaxies. Following this model's success, Sérsic (1963Sérsic ( , 1968) introduced a generalized version that can be applied to galaxies of various morphologies. Today, there are several parametric fitting codes, including GALFIT (Peng et al. 2010), and statmorph (Rodriguez-Gomez et al. 2019).
Though the sizes of galaxies can be measured by fitting a Sérsic model to an image, various radii of interest can also be derived from radial photometric measurements. Petrosian (1976) formulated a radial profile, now named after him, as a means of defining the projected radii of galaxies. Since the Petrosian profile is a ratio of surface brightnesses, it offers a distanceindependent method of measuring galaxy radii. Petrosian radii are useful for computing radial concentrations of galaxy light profiles and for performing accurate measurements of galaxy photometry . Petrosian radii and concentrations can also be used to estimate the parameters of Sérsic profiles.
By combining parametric fitting and radial profile analysis, accurate fluxes, shapes, and sizes can be derived for studying the properties and evolution of galaxies. In this paper, we introduce PetroFit, a Python package based on the Astropy (Astropy Collaboration et al. 2013) and Photutils (Bradley et al. 2020) packages. The prime motivation for constructing this package was to make a robust and scalable software package for modeling Sérsic and Petrosian profiles. PetroFit is deigned to work within the Astropy ecosystem with a focus on measuring accurate galaxy properties. Following (Crawford 2006), PetroFit uses correction grids to improve the size estimates produced by Petrosian profiles. PetroFit also provides a robust Python-based parameter fitter code based on Astropy models for parameterization of galaxy properties.
This paper defines photometric terms and models used by PetroFit in Section 2. We outline the core features of the package in Section 3 and discuss how measurements are computed. We test and demonstrate various applications of the software in Section 4. We make qualitative comparisons to similar software available to the Astronomical community in Section 5. Throughout this paper, the AB magnitude system is used. Python Notebooks used to perform tests and generate plots can be found on GitHub in the PetroFit/petrofit papers repository (Geda et al. 2022a).

PHOTOMETRIC RELATIONS
In this Section, we define the morphological and photometric relations that can be computed using PetroFit. These properties can be best understood and categorized as parameters of Sérsic and Petrosian profiles, which we discuss in detail in the following subsections. We also introduce different techniques for accurate photometric measurements.

Sérsic Profiles
The Sérsic profile (Sérsic 1963(Sérsic , 1968) is a radial profile function that describes how the surface brightness of a galaxy varies with projected radius (r) from its center (see Figure 1). The profile is a generalized form of the de Vaucouleurs' profile (de Vaucouleurs 1948) and is often written as the following intensity profile: I (r(x, y)) = I e exp −b n r r e (1/n) − 1 (1) Where I is the intensity at position (x, y). r is the radius from the center that corresponds to (x, y). r e is the effective radius, which is equal to the projected halflight radius. I e is the intensity at the half-light radius (I e = I(r e )). n is the Sérsic index which determines the "steepness" of the profile. The constant b n is defined such that r e contains half the total flux. Figure 2 shows an example Sérsic profile with (n = 1, I e = 1, r e = 25).

Total Sérsic Flux
The total flux (projected luminosity L) contained within a Sérsic profile over a projected area A(r) = πr 2 can be computed by integrating Equation 1 as follows (Ciotti 1991): Graham et al. (1996);  solve for L(≤ r) by substituting x = b n (r/r e ) 1/n which yields: L(≤ r) = I e r 2 e 2πn e bn (b n ) 2n γ(2n, x) Where γ(2n, x) is the incomplete gamma function: Because the Sérsic profile is an analytical function, the total flux of the profile can be computed by integrating Equation 2 to infinity. Substituting the complete gamma function Γ(2n) in place of γ(2n, x) in Equation 3 yields the total Sérsic flux (Ciotti 1991;: L(≤ ∞) = I e r 2 e 2πn e bn (b n ) 2n Γ(2n)

Sérsic Index
The Sérsic index (n) determines the degree of curvature of the Sérsic profile. As the value of n increases, the more concentrated the profile becomes at smaller radii. The Sérsic model tends to a power-law with a slope equal to 5 as n becomes larger . n = 4 corresponds to a de Vaucouleurs' profile that describes elliptical galaxies and the cores of spirals well; while n = 1 gives the exponential profile, which models spiral disks well Andredakis et al. 1995;Khosroshahi et al. 2000;Graham 2001;Möllenhoff & Heidt 2001;Moriondo et al. 1998). n = 0.5 gives a low concentration profile (Gaussian) that can be used to describe seeing dominated images, model intra-cluster glow in galaxy clusters (assuming an elliptical shape), and is often used as an alternative to a Ferrer Profile (Peng et al. 2010).

Effective Radius and Effective Intensity
The effective (r e ) or half-light (r 50 ) radius is the radius that encloses half the total flux of the Sérsic profile (Ciotti 1991). As such, we define r e as the radius that satisfies the following expression: Note that the expression given by Equation 6 reduces to Equation 7. Thus b n is defined as the value that satisfies the following relation between the complete and incomplete gamma functions (Ciotti 1991): The effective intensity (I e ) is the intensity exactly at r e and can be defined as I e = I(r e ). I e determines the amplitude of the profile and it is related to the intensity at the center of the Sérsic profile (I 0 = I(0)) as follows: 2.1.4. Ellipticity, Elongation, and Orientation Real galaxies, that are well described by Sérsic profiles, are usually not perfectly symmetrical and often have elliptical distributions (Ferrari et al. 2004). We define two quantities, ellipticity (ellip) and elongation (elong), that can be used to describe the elliptical distribution of Sérsic profiles: Where a is the unit length of the semi-major axis and b is the corresponding semi-minor axis. A circle corresponds to ellip = 0 and elong = 1. Ellipticity ranges from 0 to 1 while elongation ranges from 1 to infinity. The SExtractor (Bertin & Arnouts 1996) and Photutils packages use elongations for apertures while the Astropy-modeling sub-module uses ellipticity for Sérsic models.
Ellipticities (ellip) and rotation angles (θ) are used to relate Cartesian coordinates to elliptical radii by the Astropy and PetroFit packages. For an elliptical profile centered at (x 0 , y 0 ), with an effective radius r e , the following relation is used 1 : 1 See Astropy's (v4.2) Sersic2D.evaluate function Figure 1.
A unitless one dimensional Sérsic profile is plotted with n = 1, re = 25 (green vertical lines), and Ie = 1 (black horizontal line). The flux at the center of the profile (I0) is equal to 5.3 and not 1 (i.e. I0 = Ie). The dotted blue area under the curve contains half of the total flux and is equal in value to the striped green area (that also contains half of the total flux). The Sérsic profile and the green area extend to infinity. Figure 2. The curve of growth (COG) represents the total flux enclosed at each radius and asymptotically approaches the total flux (black dash-dot horizontal line). This example shows the curve of growth of a Sérsic profile with (n = 1, Ie = 1, re = 25). The half-light radius (re = 25) contains half of the total flux that is contributed by the blue dotted area in Figure 1.

Curve of Growth
The curve of growth (COG) of a galaxy is the relationship between radius (from the galaxy center) and the total intensity within that radius. It represents the cumulative flux enclosed at a given radius. For an ideal galaxy that is well described by a Sérsic profile, the curve of growth is given by L(≤ r). Figure 2 shows an example curve of growth for a Sérsic profile with (n = 1, I e = 1, r e = 25).

Petrosian Profiles
The Petrosian profile or Petrosian index (η) is a dimensionless profile that represents the rate of change in the enclosed light as a function of radius. It was first introduced by Petrosian (1976) who conceptualized it as a means of defining the sizes (radii) of galaxies. Petrosian originally defined the Petrosian index as a ratio of the average surface brightness up to a radius divided by the surface brightness at that radius. Since η is defined as an intensity ratio, it is not affected by the surface brightness dimming effects that result from the cosmological expansion of the Universe. PetroFit follows Kron (1980); Bershady et al. (2000) and defines the Petrosian profile (η) as the reciprocal of the original formulation by Petrosian. This results in a profile that has the property where η(0) = 1 at the center of the galaxy and drops to zero at the edge. As such, the Petrosian profile (η(r)) is defined as: Where η(r) is the Petrosian index at r. I(r) is the surface brightness at r. I(r) is the average surface brightness within r. L(≤ r) is the total flux within r. A(r) is the aperture area. For elliptical apertures, the area is given as A(r) = π · (1 − ellip) · r 2 . ;  derive the Petrosian profile of a Sérsic profile as a function of the logarithmic slope of the Sérsic profile α(x, n) 2 : Where x is given by the function x(r) = b n (r/r e ) 1/n , n is the Sérsic index and γ is the incomplete gamma function.
2  uses the original definition of the Petrosian profile, thus the expression they provide is the reciprocal of Equation 13 Figure 3. A plot of the Petrosian profile (η(r)) of a Sérsic profile with (Ie = 1, re = 25, n = 1). Note that that the profile has the property where η(0) = 1 at the center of the galaxy (r = 0) and drops to zero at the edge.

Petrosian Radius and Eta
The Petrosian radius is the radius where the Petrosian profile is equal to a special Petrosian index "eta" (η petro ). Following Bershady et al. (2000), the Petrosian radius is defined as the radius where η = 0.2. Thus the Petrosian radius (r petro ) and "eta" (η petro ) can be expressed as:

Petrosian Total Flux Radius and Epsilon
The Petrosian total flux radius is, as the name implies, the radius that encloses the total flux of the galaxy. Since Sérsic profiles extend to infinity, it is not possible to define a radius that contains 100% of the total flux (that would equate to a radius at infinity). It is also not practical to consider minute flux values at very large radii because real images of galaxies contain noise. For these reasons, PetroFit considers the radius that contains 99% of the galaxy's light to be the Petrosian total flux radius (Crawford 2006). To calculate the Petrosian total flux radius, we must multiply the Petrosian radius by a constant. This constant is called "epsilon" and is set to 2 by default (Bershady et al. 2000). This default value is also adopted by the Sloan Digital Sky Survey (Strauss et al. 2002). Thus we can define the Petrosian total flux radius as follows:

Total Flux and Petrosian Corrections
The Petrosian Total Flux is computed by setting = 2 by default (see Section 2.2.2). But this approach is an approximation and is susceptible to errors that arise from different galaxy morphologies. In particular, is dependent on how concentrated a galaxy's light profile is. The more concentrated a galaxy's light profile is, the larger the value of needed to produce r 99 . To understand this, let us consider two profiles, A and B, with the same total flux radii (i.e equal r 99 ), but with different concentrations (i.e n(A) > n(B)). Since profile A has a higher concentration, it will have a r petro closer to the center of the profile relative to r 99 simply because its profile is more concentrated. Thus the multiplication factor ( ) needed to produce the total flux radius (r 99 or r total ) for A is relatively larger than B. Thus, to accurately estimate r 99 , we must apply corrections by first estimating the appropriate value. Once the appropriate value is estimated, r total (which PetroFit defines as r total ≡ r 99 ) can be computed by r total = r petro · .

Petrosian Half-Light Radius
The Petrosian half-light radius is the radius that contains half of the Petrosian total flux. This quaintly is especially important because it is approximately (exactly in ideal cases) equal to the Sérsic effective radius. The Petrosian half-light radius can be computed by numerically finding the radius that encloses half of the Petrosian total flux. Consequently the Petrosian half-light radius (r 50 ) can be defined as the radius that satisfies the following expression:

Concentration Index
The concentration index is the ratio of the radii containing a specified fraction of the total flux and can be generalized as follows (Bershady et al. 2000;Kent 1985): where r o and r i are the outer and inner radii, enclosing o% and i% percent of the total flux respectively. The concentration index is related to the profile shape and can be used as a proxy for morphology (Crawford 2006;Graham et al. 2001). Some commonly used concentration indices are C 2080 (Bershady et al. 2000;Conselice 2003;Kent 1985) and C 5090 (Blanton et al. 2001): C 2080 correlates better with the Sérsic index (and epsilon). However, in the presence of image degradation, C 5090 may be more immune from systematic effects.

Petrosian Concentration Index
We define the concentration of radii derived from Petrosian indices as follows: where r(η 0 ) and r(η i ) are the outer and inner Petrosian radii, corresponding to η 0 and η i Petrosian indices respectively. This definition of concentration is useful because it can be computed directly from the Petrosian profile, without the need to estimate radii using the curve of growth (see Appendix B.2.4).
We also define a Petrosian concentration index, P 0502 , as follows: PetroFit is a Python-based package for computing and fitting various galaxy profiles. It builds on packages such as Astropy and Photutils to provide a set of tools for end-to-end profile construction and parameter fitting. As an open-source package, PertroFit is made available to the public and can be installed via PyPi. We discuss some of PetroFit's key features in this Section.

Photometry Tools
In order to make accurate measurements of galaxy properties, we need to make accurate photometric measurements while processing images of galaxies. Photutils (Bradley et al. 2020), an Astropy affiliated package for photometry, offers a wide array of tools to perform photometry. We modify or build on the tools in Photutils, and in this subsection, we introduce some of the photometry utilities we have made available in PetroFit.

Image Segmentation and Source Catalogs
It is often necessary to identify and make a catalog of sources in an image before making photometric measurements. It is also useful to identify pixels in the image that are associated with a given source so a mask can be constructed for that source. This is used to exclude the light from neighboring objects that would otherwise contaminate or bias the light profile. The PetroFit make catalog function uses photometry routines in Photutils to make segmentation maps, deblend sources that overlap in the image, and produce a catalog with useful photometric information about the sources.
To identify sources, the image is partitioned into multiple segments and a segmentation map is constructed. To achieve this, the image is first smoothed with a Gaussian kernel if a kernel size is provided and segmented using Photutils' detect sources functionality. Photutils constructs a segmentation map by assigning a label to every pixel in the smoothed image such that pixels with the same label are part of the same source. A minimum source size can be specified by providing a cutoff for the number of connected pixels that have values above a user-provided threshold. Overlapping sources are separated using Photutils' deblend sources function that uses a combination of multi-thresholding and watershed segmentation. The resulting catalog and segmentation map can then be used to identify sources of interest and mask pixels of nearby sources when performing photometry measurements. PetroFit also provides plotting tools to visualize Photutils segmentation maps.

Background Noise Subtraction
To construct a curve of growth that converges to the total flux of a galaxy, we must carefully estimate the noise in the image and make the appropriate subtractions. This is especially important for faint galaxies that have low signal-to-noise ratios. To deal with this, PetroFit offers an "image continuum" fitting tool, fit background, that fits background pixels using an Astropy model, and then that model can be subtracted from the image.
The background subtraction tool first makes a background image by masking all the sources identified in the segmentation stage. Note that our definition of background here is pixels that do not belong to a source. Once the background image is made, background pixels above a specified threshold value are also masked. The remaining unmasked pixels are used to fit an Astropy model. The default model is set to a Astropy Planar2D model.
The plane model is fit using a linear least-square fitting algorithm provided by Astropy's modeling module and directly sampled (without integrating or oversampling) into a two-dimensional model image. The user can then subtract the background image from the initial image to produce a background-subtracted image.

Aperture Photometry
The curve of growth and the Petrosian η function are both radial profiles. As such, photometric measurements must be made using apertures of varying radii from the center of the source. It is also important to track the aperture areas and account for masked pixels since the Petrosian is computed using average surface brightness (see Equation 12 and Section 3.2.2).
The PetroFit photometry step function constructs multiple concentric elliptical apertures and makes accurate flux measurements using Photutils. The elliptical shape of the apertures can be set by providing an elongation and rotation angle (default values correspond to a circle). The radii of the apertures are defined by the user and provided as a python list or array. PetroFit offers utility functions that can make lists of linearly or logarithmically spaced radii. Error images and mask maps can be provided for error estimation and pixel masking, respectively. The method returns the result in a tuple of three separate arrays that contain the flux measurements, the error estimates, and aperture areas. This function can also over-plot the apertures used for the measurement on the input image.
The PetroFit source photometry function is designed to work with Photutils source catalogs and segmentation maps. Given a source in the Photutils catalog, the function masks surrounding pixels that belong to nearby sources (each source pixel is assigned to a single source), optionally performs background subtraction using a two-dimensional plane model and then runs the photometry step function. The source photometry function also determines the appropriate aperture shape (i.e. elongation or ellipticity) and orientation using the information in the Photutils catalog.

Petrosian Tools
Given photometric measurements, PetroFit can compute the Petrosian profile as well as associated radii and concentrations. In this Section, we discuss how the Petrosian η(r) is discretely computed and how galaxy properties can be derived from it. We also discuss corrections applied to the measurements to provide accurate photometry.

Analytically Computed Profiles
PetroFit includes functions to compute Petrosian profiles and curves of growth analytically for ideal Sérsic profiles with known parameters. The following profiles and utilities are available as part of the package: petrosian profile (η(r)): Given the effective radius and Sérsic index of a Sérsic profile, this function computes the Petrosian profile using Equation 13. A 1D Astropy fittable model version is also available for this function but may not be useful for fitting profiles derived from real galaxy images (see Section 3.2.4).
sersic enclosed (L(≤ r)): Given effective intensity, effective radius and Sérsic index of a Sérsic profile, this function computes the total enclosed flux at an input radius according to Equation 3. This function can be used to make a curve of growth for ideal Sérsic profiles. A 1D Astropy fittable model version is also available for this function.
sersic enclosed inv (r(L)): Computes the inverse of Equation 3 by solving for r. Given effective intensity, effective radius and Sérsic index of a Sérsic profile, this function computes the radius that encloses the input fraction of the total light flux.

Discretely Computed Petrosian Profiles
Given an array of enclosed fluxes ("L") and corresponding aperture areas ("A"), the Petrosian profile can be computed discretely as follows: 1. Estimate the surface brightness by finding the average flux between the current index (i) and the last index (i − 1). Note that the gap between apertures affects the accuracy of the surface brightness at i, thus it is recommended to use apertures with radii that are incremented by a small number of pixels: 2. Estimate the average surface brightness by taking the flux at i and dividing it by the corresponding aperture area: 3. Compute the Petrosian index at i using the estimated values in steps 1 and 2: In discrete computations, the Petrosian profile can not be computed at the first index even if it corresponds to the center pixel (r[i 0 ] = 0). In real images, the surface brightness of a galaxy is binned into pixels, and to accurately determine I 0 , one would need to infinitely oversample the central region. In other words, each pixel corresponds to a total surface brightness integrated within the area of the pixel as opposed to the surface brightness at the pixel coordinates. As such, PetroFit sets the first Petrosian value to numpy.nan when returning a discretely computed array of Petrosian indices. PetroFit takes advantage of fact that the Petrosian index at the center of a galaxy is equal to 1 when computing radii internally.

Petrosian Radii and Concentration Indices
Given a curve of growth measurement in the form of arrays containing aperture radii, aperture areas, and enclosed fluxes, PetroFit constructs a Petrosian profile as described in Section 3.2.2. Once the Petrosian profile is constructed, Petrosian radii and concentrations can be computed. We describe how each property is computed in this Section. r petrosian (r petro ): The Petrosian radius is computed by identifying the radius where the Petrosian profile reaches the Petrosian eta (η petro ) value. By default, η petro is set to 0.2 as discussed in Section 2.2.1. The values of the Petrosian profile are interpolated and the radius where the profile intersects with η petro is identified as r petro . The eta value can be user adjusted and interpolation can be switched off (to find the closest datapoint) if necessary.
r total flux (r total ): The total flux radius is computed by multiplying r petrosian with epsilon (see Section 2.2.2). Epsilon is set to 2 by default and can be adjusted. r total depends on the assumed functional form of the radial profile of the galaxy and we discuss how to make these corrections in Section 3.2.4. If the image includes a World Coordinate System, PetroFit will convert the radius into arcsec units based on the corresponding pixel scale.
total flux (L(≤ r total )): The total Petrosian flux is computed by referring to input curve of growth and identifying the total enclosed flux at r total flux. If the curve of growth does not extend to the r total flux, a Numpy nan is returned. Note that one can set r total flux to the largest aperture radius available in the photometric data (r max ) by defining the appropriate epsilon value (i.e = r max /r petro ).
total flux uncertainty: The uncertainty or error in the total enclosed flux at r total flux. Only available if photometric uncertainties (astropy.uncertainty object) are provided. r half light (r 50 ): The half-light radius is computed by finding the radius that encloses half of total flux in the curve of growth (see Figure 1 and 2). This value can be converted into arcsec by providing a WCS.
concentration index (C io ): The concentration index is computed by first identifying the the radii that contain the requested fractions of the total flux and computing their concentration ratio as described in Section 2.3.2. In this figure, we plot epsilon ( ) vs. C2080 (uncorrected) for simulated 2D Sérsic profiles with varying effective radii (black points). The scatter seen at higher concentrations is due to sampling errors. The red curve is a line of best fit given by Equation B9 in Appendix B.2.3. The C2080 domain corresponds to 0.5 ≤ n ≤ 5 in Sérsic index.  Figure 4, but with the profiles convolved with a HST WFC3 F105W PSF (sampled at 60mas). It also shows the same curve as Figure  4 for reference. Notice how and C2080 values of profiles with smaller effective radii are affected (skewed) the most due to the PSF convolution. For this reason, we must make corrections that account for effective radii and the shape of the PSF itself.

Petrosian Corrections
To first order, the Petrosian magnitude is independent of redshift in the local Universe as it depends on the ratio of two surface brightnesses, but because it is based on profile shape, changes in the profile due to morphological K-corrections still apply (Crawford 2006). However, the real strength of Petrosian magnitudes vs. other variable aperture photometry such as Kron magnitudes (Kron 1980) is that the Petrosian radius only depends on the light interior to this radius. This aspect leads to small random errors in the measurement of r petro . Nonetheless, the magnitude within 2·r petro , relative to a true, total magnitude, is still profile-dependent. Although 99% of the light from an exponential profile (n = 1) is recovered within 2 · r petro , only 82% of the light for a r 1/4 profile is measured within 2 · r petro . Similar results apply to the Kron (Kron 1980) scheme, as discussed in Bershady (1995). Because of this, we need to adjust r total by finding the appropriate epsilon ( ) value (see Section 2.2.2).
The concentration index and "correct" r total are related, thus a relationship between the concentration index and can be derived. We attempt to derive this relationship by simulating Sérsic profiles with varying concentrations and measuring the radii that contain 99% of the total fluxes (i.e. r total , see Appendix B). Real images of galaxies contain distortions through passing through the atmosphere or the optics of a telescope, and the images are smeared according to the point spread function (PSF). The smearing of light results in the reduction in profile concentration, causing a skew in the relationship derived from simulated profiles (see Figures 4 and 5). To accurately calculate the total flux, further corrections must be applied to the estimates of η based on the PSF, concentration, and Petrosian radius. To make this correction, PetroFit uses a pre-generated grid of corrected and uncorrected Petrosian radii, eta values, Sérsic indices, and concentrations. It uses the grid to apply appropriate corrections based on raw uncorrected measurements. To generate the grid, a Sérsic profile of known parameters is first generated and convolved by a PSF provided by the user. Once the image is convolved with the PSF, Petrosian measurements are made and stored in a table along with the known correct parameters. This allows us to also derive relationships between the Petrosian and Sérsic parameters such as concentration indices and Sérsic indices.

Model Fitting Tools
PetroFit's fitting module contains useful tools to model and fit a wide variety of galaxy types. It uses machinery provided by Astropy to model data such as background gradients as described in Section 3.1.2. Another useful feature of PetroFit is that it can fit the light profiles of galaxies by constructing Astropy supported models or a combination of them (compound models). We discuss how PetroFit fits galaxy light profiles in detail in this Section.

Astropy Models
The Astropy modeling module provides a framework for representing models and parameter fitting. Several predefined 1-D and 2-D models are available along with the ability to define custom models (See Appendix A). Different fitting algorithms can be used depending on the model. Astropy's Sersic2D model implements the Sérsic model as described in Section 2.1 and allows users to provide initial guesses for the Sérsic parameters. Making good estimates of these parameters is important because the Levenberg-Marquardt algorithm is sensitive and may return parameters that correspond to a local (rather than global) minimum. Astropy's Sersic2D model has the following parameters: amplitude: Surface brightness at r eff (I e ).
theta: Rotation angle in radians, counterclockwise from the positive x-axis (θ).

PetroFit PSFModel
PetroFit's PSFModel 3 is a 2D image wrapper model that can be fit using Astropy's fitting algorithms.
PSFModel is a subclass of Astropy's Fittable2DModel and converts regular Astropy 2D models into images before each fitting iteration. It can wrap any Astropy based 2D model (including compound models) and inherits the base model's parameters. The base model can be retrieved using PSFModel.model, which would return the base model but with the current PSFModel parameters. During the fitting process or at evaluation, PSFModel samples the base model and translates it into an image via the following steps: 1. During initialization, PSFModel is given a base model, a normalized PSF image, and an oversampling parameter. PSFModel inherits the base model's parameters as well as its fixed-parameter rules and boundaries.
2. A sampling grid is generated according to the oversampling parameter. The generated grid is stored and is only regenerated if the oversampling factor is changed or the oversampling parameter is a function of the model parameters.
3. A model image is constructed using the sampling grid from the step before. To achieve this, the base model is evaluated onto the grid using the input parameters provided by the user or fitter.
4. If the model was oversampled, the model image or oversampled region is mean block-reduced to the data resolution. See the Section 3.3.3 for more details. 5. The model image is convolved with the input PSF and the result is returned. The PSF should be normalized and have the same pixel resolution as the image. The PSFModel has a parameter psf p, which controls the rotation angle of the PSF image before convolution, in addition to the parameters from the base model. Using psf p, the fitter can rotate the PSF to produce a better result (lower residuals). Enabling psf p can be useful for concentrated sources or star-like objects that are being fit using PSFs obtained from a separate data source (diffraction spikes may not align). This feature can be disabled by adding psf p to the model's dictionary of fixed parameters. If a PSF is not provided, the model image is returned without further processing.

Oversampling
Visible and near-visible photons from galaxies are measured using optical telescopes and imaging devices such as CCDs or CMOS cameras. Unlike intensity models which enjoy infinite angular resolutions, images are composed of pixels that are limited in resolution. Since we are using images of galaxies as our fitting data sets, we sometimes need to make sure that analytical models expressed in intensity (surface brightness) are sampled onto pixels (fluxes) with the appropriate integration of the model. Because integrating such intensity models over each pixel's angular footprint is computationally costly (especially for fitting routines), PetroFit employs oversampling schemes to discretely integrate the models.
One of the advantages of using PetroFit's PSFModel is its ability to sample models onto model images. For images that suffer from poor resolution, PSFModel samples the intensity model at sub-pixel levels and computes area-weighed sums (mean block-reduce) to approximate the integrated flux of each pixel (see Equation 25). PSFModel can oversample the entire model image or a specific pixel region of the image. An oversampling factor and region can be specified and passed Figure 6. This plot shows three sampling grids; the blue circles represent one-to-one sampling, the orange squares represent an oversampling factor of two and the green "x" marks represent an oversampling factor of four. The boundaries of a pixel are repented by the solid lines and the pixel value is calculated by finding the mean of the sampling points. For clarity, only points of the same oversampling factor are used (one-to-one points are not sampled if oversampling factor is two, etc.).
as an argument when wrapping an Astropy model or during run-time by setting the PSFModel.oversample attribute. The following oversampling schemes are currently available: Disabling Oversampling: To disable oversampling or one-to-one sampling, the oversampling can be set to Python None.
Oversample Entire Model Image: To oversample the entire image by an oversampling factor, the oversampling factor can be passed as a single integer value.
Oversample a Fixed Region: To oversample a fixed square region of finite size the center pixel, length of the square region and an oversampling factor can be specified as a Python tuple.
Oversample a Moving Region: While the model is being fit, the center of the model is likely to move around. To account for this, the names of model parameters that define the center of the box can be specified such that a square region of finite size "follows" the model center and oversamples that region.
For a grid that is oversampled by a factor of N and mean block-reduced to match the data resolution, each pixel value (F ) can be given with the following expression (see Figure 6): Where: N is the oversampling factor. i and j are the pixel indices.
I is the intensity profile as a function of x and y. x An issue that may come with sampling models onto images is making sure the peak value of the intensity profile is sampled into a pixel (in case it is between the sampling points). For concentrated profiles, this is especially important. For this reason, PetroFit allows for the center of the profile or model to be sampled by the subpixel nearest to it (replaces the sub-pixel value before reduction to data resolution).

Levenberg-Marquardt Fitting
PetroFit uses Astropy's modeling fitting algorithms to fit models such as PSFModel. A Levenberg-Marquardt least-squares fitting algorithm that utilizes a damping parameter to interpolate between the Gauss-Newton algorithm and the method of gradient descent (Levenberg 1944;Marquardt 1963) is used for nonlinear two dimensional models. Astropy implements this algorithm using SciPy's optimize.leastsq function, which is a wrapper around MINPACK's lmdif and lmder algorithms (Moré et al. 1980).

DEMONSTRATION OF THE SOFTWARE
In this Section, we showcase capabilities outlined in the prior Sections by testing the software on simulated and real data sets. To test the accuracy of PetroFit's parameter estimates, we perform photometry on synthetic images of galaxies (generated by sampling Sérsic models onto pixel grids) and predict their Sérsic parameters. We also test PetroFit's fitting capability by modeling the light profiles of real galaxies in optical and near-infrared images taken by the Hubble Space Telescope and the Sloan Digital Sky Survey (SDSS).

Synthetic Sérsic Images
We first demonstrate photometric measurements made by PetroFit in ideal conditions by simulating images of Sérsic profiles with no noise. We simulate Sérsic profiles described in Section 2.1.2, namely Gaussian (n = 0.5), exponential (n = 1), and de Vaucouleurs' (n = 4). To simulate the images, we defined Astropy Sersic2D models and placed the center of the models at the center of the output image. We set the image dimensions to 3500 pixels, a large image size in order to include the flux without encountering any image edge effects in the measurement. For each simulated profile the parameters of the Sérsic model were set to (I e = 1, r e = 50, x 0 = 1750, y 0 = 1750, ellip = 0, θ = 0). We used PSFModel to sample the model onto a pixel grid (note that PSFModel is independent of the photometric measurement machinery), and applied an oversampling factor of 10 to the center 100 × 100 pixels for all the model images.
Photometric measurements were then performed on the simulated profiles. First, the synthetic images were segmented using the make catalog catalog function. The curve of growth for each image was constructed using source photometry without background subtraction. A Petrosian profile was constructed using the curve of growth. The total flux radius was defined using the Petrosian profile and the default value of 2. We refer to the radii and concentrations derived using the default value as "uncorrected." We also use a correction grid to generate an improved estimate of using the PetrosianCorrection class.
Tables 1, 2 and 3 show the resulting measurements. The "Sersic" column represents analytical values derived from equations in Section 1, the "Uncorrected" column shows measurements derived using = 2, and the "Corrected" column shows values derived using improved estimates of . The values in parentheses show the ratios of measurements to their corresponding analytical values. As can be seen in these results, the corrected photometric values provide reasonably accurate measurements for these ideal cases.

Single Component Galaxy
Real images of galaxies contain background noise and are smeared by a PSF. In this subsection, we demonstrate how the different tools in PetroFit can be used together to do a comprehensive analysis of galaxy light profiles. We first perform photometric measurements and estimate the Sérsic parameters for the galaxy. Using the estimated Sérsic parameters as an initial guess, we construct a PSFModel and fit the galaxy image.
The data we use in this Section is a cutout of a group of bright galaxies in Abell 2744. The original data has a resolution of 60 milliarcseconds per pixel and was acquired by the Hubble Frontier Fields (Lotz et al. 2017) team via the WFC3 instrument (using the F105W filter). This data-set can be directly downloaded from the Mikulski Archive for Space Telescopes. The galaxy Synthetic Sérsic Image Measurements we selected for this demonstration is an elliptical galaxy (α = 3.596248 • , δ = −30.388517 • ) that can be well described by a single component Sérsic profile (see the center of the first panel in Figure 7). A 150 × 150 pixel cutout of the target was made for this demonstration and its background was subtracted by fitting a 2D plane (see Section 3.1.2) to a 3 sigma clipped version of the cutout. We refer to this cutout as the galaxy image to avoid confusion with other cutouts.

Photometry of Single Component Galaxy
Since the target galaxy has nearby neighbors, we first made a catalog and segmentation map using the make catalog function. The threshold was computed using Astropy's sigma clipped stats with a 2 sigma cutoff (this resulted in the best deblending). A smoothing kernel of size 3 × 3 and f whm = 3 was used for the segmentation and deblending steps. The npixels was set to 4 pixels and the deblending contrast was set to 0.
We constructed a curve of growth (COG) using the source photometry function. To construct the COG we used an array of 70 aperture radii up to 70 pixels (1 aperture per pixel). We set the photometry cutout size (not to be confused with the galaxy image) to twice the maximum pixel. source photometry was used to make a measurement cutout of the specified size, mask nearby sources and make photometric measurements using the defined apertures. Once the COG was constructed using the photometric measurements, the galaxy's Petrosian profile was computed and a corrected value was estimated using a correction grid as described in Section 3.

Single Component Fitting
We defined a Sersic2D model and wrapped it with PSFModel. We set the oversampling box (20 pixels on each side) to follow the center of the model and oversample by a factor of 5. PSF rotation was disabled and psf pa was set to 0. Using the calculated photometric quantities, we set initial estimates for the model of the object based on the following measurements: amplitude: The amplitude was estimated by fitting an elliptical isophot at r eff using photutils.isophot. r eff: The effective radius was set using the corrected Petrosian half-light radius.
n: The Sérsic index was estimated using the correction grid through a reverse lookup of the Sérsic index that produces the uncorrected Petrosian radius and concentration (C2080).
x 0 & y 0: The center was set to the pixel with maximum flux value. ellip and theta: Ellipticity and the rotation angle were set using the Photutils measurement during the segmentation and deblending steps.
PSFModel: A cutout of a nearby, bright star was used to model the PSF. The image cutout was normalized and had a size of 51 × 51 pixels. The Astropy LevMarLSQFitter fitter was then used to fit the model. The fitter running on an Intel Core i9-10900K Processor (20MB Cache, 3.7GHz) completed the fit under 0.1 seconds. Figure 7 shows the resulting fit and associated residual.
We compare the initial guesses obtained from photometric and Petrosian measurements before the fitting to the parameters estimated by the fitter in Table 4.
Using sersic enclosed we compute the analytical total Sérsic flux of the model at infinity and compare it to the photometrically derived Petrosian total flux of the galaxy in the image. The total Sérsic flux is 21.684 mags while the corrected Petrosian total flux of the galaxy image was measured to be 21.7593 ± 0.0026 mags. The corrected Petrosian total flux of the model was measured to be 21.75899 mags, which is 93% of the total Sérsic flux (see Section 2.2.2).

Multi-Component Galaxy
We further demonstrate PetroFit's fitting capabilities by fitting the light profile of Messier 91 (M91), a barred galaxy with multiple Sérsic components. Unlike single component galaxies, multi-component light profiles require careful and complex decompositions. To achieve this, we used multi-band SDSS images (1000 × 1000 pixels covering 20 arcminutes) centered on M91 (α = 188.860221 • , δ = 14.496339 • ). After subtracting the background by fitting a plane to a 3 sigma clipped version of the image, we used the i band to estimate the morphological parameters of the source as described in Section 4.2.1. We derived rough estimates for our initial guesses by following the steps outlined in Section 4.2.2. For fitting the galaxy, we zoomed into the galaxy by creating a cutout of (300 × 300 pixels). Using the rough parameter guess we derived, we defined three Sérsic components, one for the disk and two oversampled components for the core and bar. Using a nearby star (α = 188.740014 • , δ = 14.495436 • ) as a PSF, we found that the image had enough resolution such that the PSF was not needed for the galaxy decomposition; therefore we disabled that feature for all components. For the disk, we applied a sigma mask to mask out the core and bar of the galaxy. We then fit the disk using this masked image to better constrain our initial guess. The core and bar components were oversampled by a factor of 5. For the core, we masked out the disk and bar using an elliptical mask that was visually parameterized. Lastly, the bar was fit after masking out the core and disk using a mask that was also visually parameterized. These three preliminary fits bring the components closer to the final fit. We combined the three components into a compound model and wrapped it in a PSFModel model with an oversampling factor of 5. We then fit the entire galaxy cutout, which resulted in the final i band fit.
To perform a full analysis of the light profile, we also fit the galaxy in g, r, and z bands. For each of these bands, the fit from i band was used as the initial guess of the compound model parameters. This allowed us to skip the photometry, Petrosian, and the preliminary fitting steps for these bands. After fitting the models, we combine the i, r and, g bands into an RBG image to illustrate the galaxy, fitted model, and residual (see Figure 8). The resulting residual image clearly shows the spiral arms of M91, as well as blue clusters embedded in it. The bar component over-subtracts the pro-file, mostly because the radial profile used was elliptical as opposed to "boxy" and because actual galaxy bars quickly truncate at their edges. We plan on introducing various radial profiles, including "boxy" profiles, in the future versions of PetroFit (see Section 6 and Figure  11). The final fits running on an Intel Core i9-10900K Processor (20MB Cache, 3.7GHz) completed the fit under 2.5 minutes per band. For the core of the galaxy, each band had a Sérsic index of approximately 4 while the disk and bar components had an index of 0.5.

Simultaneously Fitting Overlapping Galaxies
A key capability of PetroFit is its ability to simultaneously fit and measure the photometry of multiple sources. Images of galaxies with nearby neighbors pose a challenge because their light profiles overlap. This is especially true in cluster environments. In this Section, we demonstrate how PetroFit can be used to construct photometry catalogs of overlapping sources and fit compound models. The group of galaxies we selected for this measurement are a part of the Hubble Frontier Fields data discussed in Section 4.2.1. The group of galaxies located at (α = 3.579222 • , δ = −30.380186 • ) are faint and isolated from bright sources. The image is background-subtracted using a 2D plane model as de- The images were made using a stretch factor of 0.9 and Q = 2. The first panel shows the SDSS M91 cutout that was used for fitting, the middle panel shows the fitted model image, and the last panel shows the residual. The spiral arm, as well as blue clusters (zoomed-in white box), can also be seen in the residual image. scribed in Section 3.1.2. We first make photometric and size measurements to inform our initial guesses. We construct a catalog using make catalog with the following settings: the detection threshold was set to the standard deviation computed using sigma clipped stats (3σ clipping), a smoothing kernel size of 3 pixels, smoothing kernel fwhm of 3, minimum of 16 connected pixels above the threshold, 0 contrast and deblending enabled. This results in the identification of 11 galaxies in the image. Once a Photutils catalog is generated, we loop through the catalog and calculate the Petrosian profile of each source. From the Petrosian profiles, we compute the sizes of each source.
For the fitting step, we construct a Sersic2D model and make initial guesses for each object as described in Section 4.2.2, except we do not apply Petrosian corrections and assume a Sérsic index of 1. This is because the sources are likely to have a low concentration and to demonstrate that the initial guesses can be rough estimates. We sum the Sersic2D models together to create a single compound model. The compound model will have parameters corresponding to each submodel's parameters, resulting in 77 parameters in total (11 sub-models times 7 parameters for each Sersic2D sub-model). Fitting such a large parameter space is possible because the initial guesses help the fitter converge to a solution. We wrap the compound model in a PSFModel and set it to oversample the entire image by a factor of 4. We use the same star cutout as Section 4.2.2 as the convolving PSF. The resulting PSFModel is fit using a Levenberg-Marquard fitter. The fitter, running on an Intel Core i9-10900K Processor (20MB Cache, 3.7GHz-5.3GHz), was able to simultaneously fit the 11 component model under 4.5 minutes. Figure 9 shows the results.

Large Scale Photometry Catalogs
In this Section, we perform large-scale multi-band aperture photometry on the Hubble Frontier Fields Abell 2744 parallel field images and compare it to the ASTRODEEP catalogs (Merlin et al. 2016). To catalog and measure the sizes of sources in the image, we selected the F105W filter as our detection image. We estimated the image background by taking a patch in the image with no sources and measuring its standard deviation, which we used to define a segmentation threshold of 3σ. We segmented and deblended the image using a kernel size of 3 pixels, f whm = 3, and a minimum source size of 16 pixels. For the deblending step, we set the contrast to 1/100.
Since we planned to use Petrosian radii as apertures for our catalog, we measured the curve of growth and Petrosian radii of each source in the segmentation map we generated. Petrosian corrections were applied to sources with uncorrected half-light radii greater than 10 pixels (0.6 arcsec). If the Petrosian profile could not be measured for a source, we used two times the bounding box of the sources' segmentation as an estimate of the aperture radius. Using the Petrosian total flux radii as aperture sizes, we made photometric measurements of each source in the F435W, F606W, F814W, F105W, F125W, F140W, and F160W filters. This was conveniently possible (without having to segment each image) because images are pixel aligned and have the same pixel scales (60 mas/pixel). This allowed us to apply the same apertures to all filters. We applied local background subtractions by fitting a plane to one sigma pixels in cutouts of targets. Figure 10 shows the results from the multi-aperture measurements. The first plot shows a magnitude comparison between the PetroFit measurements and the corresponding ASTRODEEP catalog.

COMPARISON TO OTHER PACKAGES
In this Section, we make feature comparisons between PetroFit and similar packages for measuring galaxy properties. We also provide comparative morphology and fitting tests in Appendix C.
Statmorph (Rodriguez-Gomez et al. 2019) is a package for calculating non-parametric morphological diagnostics of galaxy images and fitting single component 2D Sérsic profiles. Statmorph uses radial profiles, including the Petrosian profile, to compute radii of interest and concentrations. PetroFit includes features, such as correction grids, to account for PSF smearing during parameter estimations. Statmorph's offers a single component Sérsic fitting that is compatible with Astropy while PetroFit can fit multiple components of any 2D model variety. PetroFit fitting code also provides oversampling schemes for fitting images of galaxies with poor angular resolutions or profiles with high concentrations. At the time of writing, statmorph offers a variety of morphological parameters, such as Gini-M20 statistics, that are currently not implemented in PetroFit. Optimizing PetroFit's functions and including some of these features or investigating potential synergies between Statmorph and PetroFit is work for future versions of PetroFit. Though either package can be used for many overlapping science cases, PetroFit was built specifically for making precision angular size measurements and offers tools to fit complex light profiles (See Sections 4.3 and 4.4).
GALFIT (Peng et al. 2010) has and continues to serve the community with its robust galaxy light profile fitting code. This package uses modeling techniques similar to PetroFit but also offers azimuthal shape functions that can be used to define shapes beyond elliptical distributions. Similar types of capabilities are planned for implementation in future iterations of PetroFit (see Section 6 and Figure 11). PetroFit provides non-parametric estimates of galaxy properties beyond model fitting. It also provides a python-based package that is compatible with the Astropy ecosystem, allowing users to easily interact with fitted compound models. PetroFit offers tools to make initial guesses of model parameters. Lastly, PetroFit allows for custom models or parameters to be defined on the fly by virtue of Astropy's custom modeling.

CONCLUSION
In this paper, we introduce PetroFit, an open-source Python package for fitting galaxy light profiles and measuring Petrosian properties. PetroFit offers a wide range of tools for making accurate photometric measurements, constructing curves of growth, and computing Petrosian profiles. Computed Petrosian profiles can be used to estimate the projected radii of galaxies as well as their concentrations. PetroFit includes a 2D image fitting toolset that is built using the modeling package in Astropy. We demonstrate capabilities, limitations, and potential use cases for the software by analyzing images of simulated and real galaxy light profiles.
The version of PetroFit described in this paper represents an initial step towards our goal of a robust package for making galaxy size measurements. We plan on continuously improving the software and hope to gain useful feedback from the Astronomical community. In the next iterations, we hope to add features such as better error estimations in fitted models, Astropy unit support for Petrosian measurements, and advanced ma-chine learning-based parameter estimates (an experimental decision tree regressor is currently available for Petrosian corrections). In future versions of PetroFit, we hope to allow for corrections to be computed using 80 and P 0502 . We also hope to officially implement azimuthal shape functions defined by Peng et al. (2010), which we refer to as PHIR (Peng-Ho-Impey-Rix) models, which use Fourier and bending modes to create complex profiles (See Figure 11 for examples of simple azimuthal Astropy models). We aim to make PetroFit more compatible with similar packages for the purposes of interoperability and testing. As part of this effort, we also hope to contribute to PhotUtils and Astropy features that have been developed here that may be useful beyond the goals of PetroFit.
Due to the ability to accurately measure galaxy photometry and sizes systematically, PetroFit can support a wide range of galaxy evolution studies including classification of galaxies, assessing changes in the size of galaxies, and measurements in the luminosity function of galaxies. The flexible, object-oriented design of the code makes it ideal for use in the next generation of high-resolution galaxy surveys that will be produced by JWST and the Roman Space Telescopes.
We openly invite and strongly encourage all interested parties to leave feedback, bug reports, and feature requests by opening tickets in the PetroFit/petrofit GitHub repository. We hope this package is as useful for the general astronomical community as it has been for the authors. Figure 11. Examples of azimuthal shape functions implemented as custom 2D Astropy (fittable) models. Inspired by PHIR (Peng-Ho-Impey-Rix) models (Peng et al. 2010), the first and middle panel show Sérsic models that apply simple rotations to the profile as a function of radius (i.e θ(r[x, y], m, rin), where m is the number of full rotations and rin is the rotation cutoff). The last panel shows a barred Sérsic profile that was generated using a generalized ellipse.

ACKNOWLEDGMENTS
We thank our collaborators Dr. Gregory D. Wirth and Dr. D.J. Pisano for their close guidance and support throughout this project. We would like to thank Alan Kahn (Department of Physics, Columbia University) for useful conversations about machine learning and its potential applications. We would like to thank the developers of PhotUtils as much of this work is built on top of the tools they provided.  Box2D 2D Box profile.

Disk2D
A Disk profile with a constant amplitude.

TrapezoidDisk2D
Disk with a slope.

Ellipse2D
Ellipse with major and minor axis and rotation angle.

Gaussian2D
2D Gaussian (can be elliptical) distribution with a rotation angle.

Moffat2D
Moffat profile (can be elliptical) with alpha (power index) and gamma (core width).

RickerWavelet2D
Symmetric Ricker Wavelet function with the specified sigma.

Sersic2D
Sérsic profile with an effective half-light radius, rotation, and Sérsic index.

Ring2D
A 2D Ring with inner and outer radii.

Model Name Description
Nuker2D 2D Nuker profile In this Section, we discuss various methods of approximating Sérsic and Petrosian quantities. In particular, we look at the relationships between Sérsic indices, corrected and uncorrected concentrations. We also cover why correction grids are necessary despite the availability of the approximations discussed in this Appendix.

B.1.1. Sérsic Indices and Concentrations
For purely analytical Sérsic profiles that extend to infinity, the concentration can be found by solving for the radii where L(≤ r i ) = i 100 L(≤ ∞) and L(≤ r o ) = o 100 L(≤ ∞). Using Equations 3 and 5, along with x(r) = b n (r/r e ) 1/n , concentration radii can be expressed as follows:

B.2. Approximations
We approximate the following relationships by first generating a grid of quantities using ideal Sérsc and Petrosian profiles for 0.5 ≤ n ≤ 15. We then make corrected and uncorrected Petrosian measurements. Note that "uncorrected" means that Petrosian radii and concentrations were computed using the default = 2. For each value of n, we followed the steps below to generate the profiles and quantities used to derive relationships: 1. Calculate Sérsic radii containing fluxes of interest, including r 99 sersic (which we define as the total flux radius), using sersic enclosed inv.
2. Construct a radius list to sample analytical functions discussed in Section 3.2.1.
3. Use the radius list and sersic enclosed to construct a curve of growth (flux array).
4. Use the radius list and petrosian profile to construct an ideal Petrosian profile. We over-sample the inner r 80 by a factor of 10000 to produce accurate measurements for high concentration profiles. For radii larger than r 80 , we sample at a maximum of 1000 points (equally spaced).
7. Find the corrected value using correct = r 99 sersic /r petro .
8. Using correct and r petro , compute the corrected r total , r 50 , r 20 , r 80 and C 2080 .
The approximations in the Subsections that follow were derived by fitting models using Astropy modeling.

B.2.1. Corrected vs Uncorrected Concentration Index
Because uncorrected concentrations are derived from radii calculated using = 2, errors in the radii estimates propagate to estimates of concentrations. This effect is especially severe for profiles with high concentrations. We estimate corrected C 2080 using the uncorrected C 2080 for 0.5 ≤ n ≤ 15 using a sixth degree polynomial (see Figure 12). For clarity we denote the corrected C 2080 with C and uncorrected C 2080 with U : C(U ) ≈ 2.26194802 − 3.61130833 · U + 3.8219758 · U 2 − 1.6414660 · U 3 +0.38059409 · U 4 − 0.0450384 · U 5 + 0.00221922 · U 6 (B5) Figure 12. Corrected C2080 approximated from the uncorrected C2080 for 0.5 ≤ n ≤ 15. The black markers represent values obtained from measurements and the red line is the line of best fit.

B.2.2. Sérsic Index vs. Concentration Index
Since both the concentration index and Sérsic index are measures of concentration, they are correlated. We estimate the Sérsic index n using the corrected C 2080 for 0.5 ≤ n ≤ 15 using a fifth degree polynomial (see Figure 13). The relationship between the uncorrected C 2080 and Sérsic index n can be computed first computing the corrected C 2080 using Equation B5. For clarity we denote the corrected C 2080 with C and uncorrected C 2080 with U :

B.2.3. Epsilon vs. Concentration Index and Sérsic Index
Epsilon ( ) is the multiplying factor that is used to compute r 99 . We estimate epsilon ( ) using the Sérsic index n using a combination of a fifth degree polynomial and an exponential (see Figure 14). We also use the formulations in Subsections B.2.2 to compute the relationship between and C 2080 . For clarity we denote the corrected C 2080 with C and uncorrected C 2080 with U : (n) ≈ −6.54870813 − 2.15040843 · n − 0.28993623 · n 2 − 0.04099376 · n 3 −0.00046837 · n 4 − 0.00022305 · n 5 + 7.48787292 · exp n 2.6876055 (B8) Figure 13. Sérsic index n estimated using the corrected and uncorrected C2080 for 0.5 ≤ n ≤ 15 using Equations B6 and B7 (red lines).

B.2.4. Epsilon vs. Petrosian Concentration Index
Though it is difficult to approximate the corrected value using the uncorrected Concentration Index (C 2080 ), we can use the Petrosian concentration index P 0502 (that is an uncorrected quantity, see Section 2.3.3) to estimate . We use a combination of a 6 th degree Polynomial and a single exponential for this approximation. For clarity we denote the Petrosian concentration index P 0502 with P (see Figure 15): (P ) ≈ 1.09339566 − 0.14524911 · P + 0.50361697 · P 2 − 0.1215809 · P 3 + 0.02533795 · P 4 −0.00196243 · P 5 + 0.00009081 · P 6 + 0.03312881 · exp P 1.83616642

(B10)
This approximation is useful because the radii used to compute P 0502 do not depend on itself. This means that P 0502 can be computed directly from the data (no need for corrections) and its value can be used to compute without the need to know any of the radii that enclose a fraction of the total flux (i.e r 20 , r 80 , r 50 , r 90 , etc.).   gives the following approximation for r e using the uncorrected r 50 and r 90 /r 50 (for Sérsic profiles with 0.1 ≤ n ≤ 10 ):

B.3. Recovering Total Sérsic Magnitudes
Because Petrosian profiles depend on the morphology of galaxies, the Petrosian total flux radius computed using the default parameters (r total = r(η 0.2 ) · 2) often results in an aperture radius that underestimates the total Sérsic magnitudes (see Section 3.2.4).  discusses in detail the relationship between morphology, total Sérsic magnitudes, and Petrosian profiles. Using Equation 3, we can define the total Sérsic magnitude (m s ) and the total Petrosian magnitude (m p ) as a function of the Petrosian radius (r petro ) and as follows: m s (≤ r) = −2.5 · log 10 (L(≤ r)) (B12) m p (≤ · r petro ) = −2.5 · log 10 (L(≤ · r petro )) (B13) Using the uncorrected r 50 and r 90 /r 50 ,  (their Equation 6) provide the following approximation for the "missing" magnitude (for Sérsic profiles with 0.1 ≤ n ≤ 10): ∆m ≡ m p − m s = P 1 · exp (r 90 /r 50 ) P2 Where P 1 and P 2 equal to 5.1 · 10 −4 and 1.451. We build on this work and extend the domain of the approximation to 0.1 ≤ n ≤ 15 using a combination of an exponential and a polynomial. We provide this approximation as a function of uncorrected C 2080 (U ) that estimates the "missing" magnitudes within 0.01 mags (see Figure 16): ∆m ≈ −0.04990413 + 0.07810246 · U − 0.04810471 · U 2 +0.00929017 · U 3 + 3.241416 · 10 −6 · exp U 0.519920 (B15) Figure 16. The first panel shows the fraction of the total corrected flux that is recovered by using the default epsilon value ( = 2). The second panel shows a line of best fit that approximates the difference in corrected and uncorrected magnitudes (i.e ∆m) for 0.5 ≤ n ≤ 15.

B.4. Effects of PSFs on Approximations
Point spread functions (PSFs) and seeing effects greatly diminish our ability to approximate Sérsic profile parameters and Sérsic radii. This is because these optical effects smear the light coming from astronomical sources causing a change in their intensity profiles (see Figure 17) and curves of growth. For example, the corrected half-light radius is not equal to the Sérsic half-light radius if the profile is convolved with a PSF (which stretches out the profile). Since each PSF smears the intensity profiles of galaxies in a unique way (see Figures 4 and 5), a correction unique to the PSF must be applied to reconstruct the Sérsic radii. This adds an extra layer to error corrections and is the reason why we must generate a grid, using a known PSF, that maps the smeared uncorrected/corrected measurements to the intrinsic Sérsic parameters and Sérsic profile derived quantities. Figure 17. Three 1D Sérsic profiles are defined with re = 25 but with varying Sérsic indices. Each of the three profiles is convolved with the same PSF modeled by a Moffat function. Notice how the de Vaucouleurs' profile (n = 4) is most affected by the PSF convolution.

C. MORPHOLOGY AND FITTING COMPARISONS
Here we perform quantitative tests to compare and demonstrate the advantages of PetroFit (Geda et al. 2022b)(v0.3.0), statmorph (v0.4.0), and GALFIT (v3.0.5). We chose statmorph for this comparison because it is an Astropy and PhotUtils based package that can produce morphological diagnoses that can be compared to PetroFit's radial and concentration measurements. For fitting, GALFIT was the best package to compare to because it has been extensively tested and is widely adopted. We run two kinds of tests, one on morphology-related measurements and a second on fitting. For the morphology tests, we ran PetroFit and statmorph on ideal Sérsic profiles (convolved with a Moffat PSFs) which were simulated using GALFIT. For the fitting test, we ran statmorph and GALFIT on the dataset discussed in Section 4.2 to see if we could reproduce similar results. To limit the tests to the fitting or morphology parameters, we provided statmorph the same segmentation maps as PetroFit and used the same background-subtracted images as input for all three packages. Further details are discussed in the following subsections.

C.1. Fitting Test
As a demonstration of PetroFit's fitting libraries, we compare the results from Section 4.2.2 to fits we perform using GALFIT and statmorph. We provided both GALFIT and statmorph with the same background-subtracted image used by PetroFit. For statmorph, we follow the same segmentation step as used by PetroFit. We provided GALFIT with the same source mask used by PetroFit. We also provided GALFIT with parameters estimates close to those of PetroFit and statmorph (after fitting using the later packages). To make fair timing comparisons, we consider run-times that include loading the data from file and refer to it as "total run-time". We do not include segmentation time since all three packages use maps and masks from the same segmentation steps. Running on an Intel Core i7-2600 CPU at 3.40GHz, PetroFit finished with a total run-time under 1.2 seconds, statmorph under 0.8 seconds, and GALFIT under 1.8 seconds. The fit results are given in Table 5. Radii are given in pixels while photometric flux measurements are given in electrons per second (except magnitudes, which are given in the AB system). The results show similar outputs as expected. PetroFit and statmorph produced similar models and upon visual inspection, show the least amount of residual.

C.2. Morphology Test
For this test, we first generated mock galaxies using single component Sérsic profiles and an ideal Moffat PSF. The Sérsic profiles shared the same parameters (r ef f = 15 pix, ellip = 0, θ = 0), except for the Sérsic index and effective intensity. The Sérsic index was set to (0.5, 1, 4, or 5) and the effective intensity was defined such that the total Sérsic flux of the profile was equal to one (total flux equal to 0.99 by PetroFit definition, see Section 2.3). The Moffat profile was generated by first fitting a real Hubble F105W PSF shown in Figure 7 and converting that model into an image. This was done so we could create a PSF with no noise and exactly centered in the image. Using the model parameters and the Moffat PSF, we used GALFIT to create model images of size 900 on each axis. A correction grid was generated using the new PSF for PetroFit. We carried out the segmentation and deblending step using the make c atalog function. This gave us PhotUtils segmentation maps that we used for both tools (both tools used the same maps). We ran statmorph using its SourceMorphology class. We placed PetroFit's photometry, morphology, and fitting steps in a single function for this test. For profiles with a Sérsic index greater than three, PetroFit oversampled the model by a factor of 10 within a box of 50 by 50 pixels of the center. Both functions computed morphology parameters and fit the input images using the Moffat PSF.
The following tables show the results of the tests. We show the results from each profile in a separate table.
The "Values" column shows the parameter name and the "GALFIT Inputs" column shows Sérsic parameters that were used to create the image. Values such as r 20 and "Total Flux" in the "GALFIT Inputs" column were calculated using analytical equations defined in Sections 1 and 2. It is important to note that these ideal values are proxies since they do not account for PSF smearing (see Appendix B.4), thus radii close to the center (such as r 20 ) will be skewed. Lastly, each table under the "PetroFit" and "statmorph" columns show the run-time of each tool in seconds (running on an Intel Core i7-2600 CPU at 3.40GHz). Radii are given in pixels while photometric flux measurements are given in electrons per second. The run-times show similar run-times for this sample, with PetroFit performing better with high index sources. As stated in Section 5, statmorph computes additional parameters, such as Gini-M20 statistics, which contribute to its run time. PetroFit's correction grids allow it to get much better estimates of r 80 and C 2080 for highly concentrated profiles (see n = 4 and 5). For example, for a profile with n = 4 statmorph measures a total flux of 0.83519 e − /s which is 84% of the total Sérsic flux (this is expected, see Figure 16). It's also important to note that PetroFit was able to predict these values before fitting the images. Another important result is that oversampling the high index images allowed PetroFit to produce better fits for n.