Fitting

Most galaxy light profiles can be well described by PSF convolved models like the Sérsic profile. PetroFit uses the astropy modeling sub-module to provide tools to perform two-dimensional fits of galaxy light profiles. To this end, we use the PetroFit PSFConvolvedModel2D class, which applies PSF convolution to and handles oversampling for astropy based models.

In this section, we demonstrate the basics of light profile modeling on a galaxy using a single component Sérsic profile. We also demonstrate how the photometry and petrosian PetroFit sub-models can be used to improve the initial guesses of the Sérsic parameters.

Loading Example Data

The following data is a cutout of a group of bright galaxies in Abell 2744 (located at (3.596248, -30.388517)). The original data was acquired by the Hubble Frontier Fields team via the WFC3 instrument in the F105W filter and can be directly downloaded from the Mikulski Archive for Space Telescopes. The cutout image used in this documentation can be found in the git repository at the following path petrofit/docs/data/abell_2744_dwarf_galaxy_f105w.fits.gz.

Loading Data and RMS Images

We first use astropy’s CCDData to load the example data and visualize it through matplotlib. The rms image is loaded using astropy’s fits sub-module.

[1]:

from astropy.nddata import CCDData
from astropy.io import fits

image = CCDData.read('data/abell_2744_dwarf_galaxy_f105w.fits.gz')
rms = fits.getdata('data/abell_2744_dwarf_galaxy_f105w_rms.fits.gz')

[3]:

import numpy as np
from matplotlib import pyplot as plt

plt.rcParams['figure.figsize'] = [10, 10]
plt.rcParams['image.origin'] = 'lower'

vmax = 0.05 # Use the image std as max and min of all plots
vmin = - vmax

plt.imshow(image.data, vmin=vmin, vmax=vmax)
plt.title("Galaxy in Abell 2744")
plt.xlabel("Pixels")
plt.ylabel("Pixels")
plt.show()

plt.imshow(rms)
plt.title("RMS Image")
plt.xlabel("Pixels")
plt.ylabel("Pixels")
plt.show()

Make Source Catalog

We will use the sigma clipped std as a threshold at the segmentation and deblending steps.

[4]:

from astropy.stats import sigma_clipped_stats
image_mean, image_median, image_stddev = sigma_clipped_stats(image.data, sigma=3)

Here we identity sources in the input image.

[5]:

from petrofit.segmentation import make_catalog, plot_segments

threshold = image_stddev * 2

# Define smoothing kernel
kernel_size = 3
fwhm = 3

# Min Source size (area)
npixels = 2**2

cat, segm, segm_deblend = make_catalog(
    image.data,
    threshold=threshold,
    deblend=True,
    kernel_size=kernel_size,
    fwhm=fwhm,
    npixels=npixels,
    contrast=0.00,
    plot=True, vmax=vmax, vmin=vmin
)

plt.show()

# Display source properties
print("Num of Targets:", len(cat))

Num of Targets: 15

Single Source Photometry

Here we pick a galaxy that can be modeled using a single Sérsic model, the galaxy in the middle looks like an elliptical galaxy that can be well described by such a profile. To make sure we selected the right galaxy, we use the plot_target function to plot a cutout of the source.

[6]:

from petrofit.photometry import order_cat
from petrofit.utils import plot_target

# Sort and select object of interest in the catalog
sorted_idx_list = order_cat(cat, key='area', reverse=True)
idx = sorted_idx_list[2] # index 0 is largest
source = cat[idx]  # get source from the catalog

plot_target(
    position=(source.maxval_xindex, source.maxval_yindex),
    image=image.data,
    size=100,
    vmax=vmax,
    vmin=vmin
)

After selecting the source, we define a radius list for the aperture photometry. Once r_list is defined, we run the source_photometry step.

[7]:

from petrofit.photometry import make_radius_list
from petrofit.photometry import source_photometry

# Max aperture radius
max_pix = 70

r_list = make_radius_list(
    max_pix=max_pix, # Max pixel to go up to
    n=max_pix # the number of radii to produce (i.e 1 aperture per pixel increase in r)
)

print("len(r_list) = {}".format(len(r_list)))


# Photomerty
flux_arr, area_arr, error_arr = source_photometry(

    # Inputs
    source, # Source (`photutils.segmentation.catalog.SourceCatalog`)
    image.data, # Image as 2D array
    segm_deblend, # Deblended segmentation map of image
    r_list, # list of aperture radii

    # Options
    error=rms, # Error image (optional)
    cutout_size=max(r_list)*2, # Cutout out size, set to double the max radius
    bkg_sub=True, # Subtract background
    sigma=3, sigma_type='clip', # Fit a 2D plane to pixels within 3 sigma of the mean
    plot=True, vmax=vmax, vmin=vmin, # Show plot with max and min defined above
)
plt.show()

len(r_list) = 70
15

Petrosian Profile

We use the photometry that we obtained to construct a Petrosian profile of the galaxy. This will allow us to estimate the parameters of the galaxies Sérsic profile. We also initialize a PetrosianCorrection which we can use to estimate the Sérsic index (n) and r_eff (half-light radius). The correction grid we use for PetrosianCorrection was created specifically for this dataset using the PSF above.

[8]:

from petrofit.petrosian import Petrosian, PetrosianCorrection

p = Petrosian(r_list, area_arr, flux_arr)

pc = PetrosianCorrection("data/concentration_index_grid_f105w_60mas.yaml")

Next we compute and apply the epsilon that produces the correct r_total_flux and creates a new corrected Petrosian profile (corrected_p). To compute corrected_epsilon we use the uncorrected r_half_light and C2080.

[9]:

corrected_epsilon = pc.estimate_epsilon(
    p.r_half_light,
    p.concentration_index()[-1]

)

corrected_p = Petrosian(r_list, area_arr, flux_arr,
                        epsilon=corrected_epsilon)

Galaxies with a high Sérsic index (high concentration) have r_total_fluxs that extend much further than we can possibly measure. If the corrected r_total_flux is out of the photometric range (i.e r_total_flux > max_pix), we would be unable to measure the flux at r_total_flux (the total flux), which also means that we would be unable to measure the Petrosian radius (np.nan is returned). If the corrected r_total_flux is out of range and we feel like the photometry data range goes far out enough to estimate the total flux, we can set epsilon to a value that will result in a r_total_flux equal to the last photometry radius.

[10]:

max_photometry_r =  max(corrected_p.r_list)

if corrected_p.total_flux > max_photometry_r:
    print("Truncation was applied")
    truncated_epsilon = max_photometry_r / corrected_p.r_petrosian
    corrected_p.epsilon = truncated_epsilon
else:
    print("Truncation was not needed")

# Print Radii
print('total_flux = ', corrected_p.total_flux)
print('max_photometry_r = ', max_photometry_r)

# Plot Petrosian
corrected_p.plot(plot_r=True, plot_normalized_flux=True)
plt.show()

Truncation was not needed
total_flux =  64.0256148087154
max_photometry_r =  70.0

Estimating Sérsic Parameters

astropy’s Sersic2D implements this model and allows us to provide initial guesses for the Sérsic parameters. Getting good estimates of these parameters is very important because the Levenberg-Marquardt algorithm is sensitive and may return parameters that correspond to a local minimum. Because of this, we follow the steps below to estimate the Sérsic parameters.

The Sersic2D model has the following parameters:

amplitude: Surface brightness at r_eff.
r_eff: Effective (half-light) radius.
n: Sérsic Index.
x_0 and y_0: x and y position of the center.
ellip: Ellipticity of the profile.
theta: Rotation angle in radians, counterclockwise from the positive x-axis.

Center Pixel

We can use the get_source_position to find the max pixel position of the source in the image using the SourceCatalog object. Since this center is relative to the image, we save the result in image_x_0 and image_y_0.

Note

The x_0 and y_0 vlaues we use in the model in this example are defined in the Zoom and Mask Image section. This is because we make a cutout of the target at (image_x_0, image_y_0) which shifts the coordinate space. If no cutouts are made, then we can use (image_x_0, image_y_0) as the center of the astropy model.

[11]:

from petrofit.segmentation import get_source_position
image_x_0, image_y_0 = get_source_position(source)

print(image_x_0, image_y_0)

99 100

Ellipticity and Elongation

We pass the source’s SourceCatalog object to the get_source_ellip and get_source_elong functions to find the source’s light profile ellipticity and elongation respectively. These values are derived from the segmentation footprint of the source.

[12]:

from petrofit.segmentation import get_source_ellip
ellip = get_source_ellip(source)

from petrofit.segmentation import get_source_elong
elong = get_source_elong(source)

print(ellip)
print(elong)

0.07095716989311518
1.0763766401221262

Source Orientation

We pass the source’s SourceCatalog object to the get_source_theta function to find the orientation of the light profile in radians, counterclockwise from the positive x-axis. In some catalogs, this value is known as the position angle.

[13]:

from petrofit.segmentation import get_source_theta
theta = get_source_theta(source)

print(theta)

1.39821829717347

Sérsic Index

The PetroFit PetrosianCorrection contains a grid that maps the uncorrected (PSF convolved) Petrosian half-light radius and concentration index to an epsilon value that gives the correct Petrosian radii. This grid can also be used to map the uncorrected half-light radius and concentration index to a Sérsic index.

[14]:

n = pc.estimate_n(
    p.r_half_light,
    p.concentration_index()[-1]
)

print(n)

1.408619948312809

Half-Light Radius (r_eff)

We use the corrected Petrosian to estimate the half-light radius (r_eff).

[15]:

r_eff = corrected_p.r_half_light

print(r_eff)

8.247649529905981

Amplitude at r_eff

To find the flux at the half-light radius, we fit an elliptical isophot at r_eff using photutils’s Ellipse model using the get_amplitude_at_r function. If this value can not be computed, np.nan is returned so make sure to check using np.isnan(amplitude)

from photutils.isophote import EllipseGeometry, Ellipse

# Define EllipseGeometry using ellip and theta
g = EllipseGeometry(image_x_0, image_y_0, 1., ellip, theta)

# Create Ellipse model
ellipse = Ellipse(image.data, geometry=g)

# Fit isophote at r_eff
iso = ellipse.fit_isophote(r_eff)

# Get flux at r_eff
amplitude = iso.intens

[16]:

from petrofit.segmentation import get_amplitude_at_r
amplitude = get_amplitude_at_r(r_eff, image, image_x_0, image_y_0, ellip, theta)

print(amplitude)

0.07638307541983923

Zoom and Mask Image

Before we fit the image, it is important to mask all nearby sources to avoid interference and zoom in to avoid fitting irrelevant pixels (especially those outside of r_total_flux). For this example, we will mask nearby sources using the masked_segm_image function that masks pixels using source segmentation footprints and make a zoomed cutout using astropy’s Cutout2D function. We also make a second cutout (fitting_image_unmasked) of the same dimensions as the fitting image, but without masking nearby sources. We will subtract the fitted model from this unmasked image once the fitting is complete.

We make a cutout of size cutout_size centered at the source, as we identified in the finding center section, and define a new center (x_0 and y_0) that we will use in the astropy model.

[17]:

from astropy.nddata import Cutout2D
from petrofit.segmentation import masked_segm_image

# Define cutout size
cutout_size = 120

# Make an image with all sources masked except `source`
masked_image = masked_segm_image(source, image.data, segm_deblend, fill=None, mask_background=False)

# Make cutouts
fitting_image = Cutout2D(masked_image, (image_x_0, image_y_0), cutout_size, mode='partial', copy=True)
fitting_image_unmasked = Cutout2D(image.data, (image_x_0, image_y_0), cutout_size, mode='partial', copy=True)

# Define new center
x_0 = y_0 = cutout_size / 2

# Plot cutout that will be fit
plt.imshow(fitting_image.data, vmin=vmin, vmax=vmax)
plt.title("Galaxy in Abell 2744")
plt.xlabel("Pixels")
plt.ylabel("Pixels")
plt.show()

AstroPy Sérsic Model

We also define a set of bounds, a dictionary of lower and upper bounds of parameters. Keys are parameter names. The values are a list or a tuple of length 2 giving the desired range for the parameter and a value of None means no bounds. For example, we restrain the fitter from exploring half-light radii that are negative by adding 'r_eff': (0, None). We also restrict the center of the model to be within a range (center_slack) from the initial guess.

[18]:

from astropy.modeling import models

center_slack = 4

sersic_model = models.Sersic2D(

        amplitude=amplitude,
        r_eff=r_eff,
        n=n,
        x_0=x_0,
        y_0=y_0,
        ellip=ellip,
        theta=theta,

        bounds = {
            'amplitude': (0., None),
            'r_eff': (0, None),
            'n': (0, 10),
            'ellip': (0, 1),
            'theta': (-2*np.pi, 2*np.pi),
            'x_0': (x_0 - center_slack/2, x_0 + center_slack/2),
            'y_0': (y_0 - center_slack/2, y_0 + center_slack/2),
        },
)

PSFConvolvedModel2D

The petrofit PSFConvolvedModel2D is a Fittable2DModel that adds PSF convolution and model to image sampling to astropy core models. PSFConvolvedModel2D makes an image of the underlying model and samples it onto a grid. The model image is then convolved with a PSF if one is provided. Since PSFConvolvedModel2D is a Fittable2DModel, it can be used to fit model images to data images. For example, we wrap an astropy Sersic2D model in this doc with PSFConvolvedModel2D, which produces an oversampled and PSF convolved version of the Sérsic profile at each iteration of the Levenberg-Marquardt fitting algorithm. Note that ``PSFModel`` is deprecated and replaced by ``PSFConvolvedModel2D``.

Note

PSFConvolvedModel2D is agnostic to the models it wraps and can handle complex multi-component astropy models.

PSF

A Point Spread Function (PSF) describes how light from a point source is distributed on detector due to optical effects such as diffraction. Images or cutouts of stars are good approximations of PSFs because stars are single-point sources and their images describe how their light is distributed on the detector. To make cutouts of stars in an image, use the astropy.nddata.Cutout2D function.

The following PSF is a cutout of a star in the Hubble Frontier Fields image of Abell 2744 (same dataset as the example image). Since we will be using the PSF image as a convolution kernel, it is very important that the following requirements are satisfied:

The image of the PSF should be at the same resolution as the data.
The star or PSF is centered in the image.
The PSF image does not contain other sources.
The image is normalized so that the sum of the PSF image near or equal to 1.0.
The PSF image should have odd dimensions on each side (because it is a convolution kernel).

[19]:

from astropy.io import fits

# Load PSF image (2D array)
PSF = fits.getdata('data/f105w_psf.fits.gz')

# Normalize PSF
PSF = PSF / PSF.sum()

# Note that the PSF shape is odd on all sides
print("PSF Shape = {}".format(PSF.shape))

# Plot PSF and use vmax and vmin to show difraction spikes
plt.imshow(PSF, vmin=0, vmax=PSF.std()/10)
plt.show()

PSF Shape = (51, 51)

Oversampling

One of the advantages of using PSFConvolvedModel2D is its ability to sample models onto model images. Sometimes the models have regions that have to be oversampled to produce better estimates of the data. PSFConvolvedModel2D can oversample the entire model image or a specific pixel region of the image. The oversampling factor and region can be specified in the oversample keyword argument when wrapping an astropy model or during run time by setting the PSFConvolvedModel2D.oversample attribute.

Disable Oversampling (Defailt)

To disable oversampling, set the oversampling argument or attribute to None

[20]:

# Disable Oversampling
oversample = None

Oversample Entire Model Image

To oversample the image by a factor, you can pass a single integer value. For example:

[21]:

# Oversample the entire image by a factor of 4
oversample = 4

Oversample a Fixed Region

To oversample a fixed region of finite size, specify the center pixel, the length of the square region and thee oversampling factor. This means passing a tuple of (center_x, center_y, box_length, oversample_factor). For example:

[22]:

# Replace the pixel values in a box of
# length 20 cented at (x=50, y=60) with a box of
# the same size that has been oversampled by a factor of 5
# i.e (x=50 y=60, box_length=20, oversample_factor=5)

oversample = (50, 60, 20, 5)

Oversample a Moving Region

If the model is being fit, the center of the model is likely to move around. To account for this, we can specify the names of the model parameters that define the center of the box that we are interested in oversampling as strings. This means passing a tuple of (model_param_x, model_param_y, box_length, oversample_factor). For example:

[23]:

# Replace the pixel values in a box of
# length 20 cented at (x=model.x_0, y=model.y_0) with a box of
# the same size that has been oversampled by a factor of 5
# i.e (model.x_0, model.y_0, box_length=20, oversample_factor=5)

oversample = ('x_0', 'y_0', 20, 5)

Create PetroFit Model

Now that we have an astropy model, PSF and oversampling rule, we can create a PSFConvolvedModel2D model as follows:

[24]:

from petrofit.modeling import PSFConvolvedModel2D

psf_sersic_model = PSFConvolvedModel2D(sersic_model, psf=PSF, oversample=oversample)

The PSFConvolvedModel2D etherates all of the parameters, fixed-parameter rules and parameter bounds from the input astropy model. Notice that a new parameter, psf_pa is added to enable PSF rotation.

[25]:

print(psf_sersic_model.param_names)

('amplitude', 'r_eff', 'n', 'x_0', 'y_0', 'ellip', 'theta', 'psf_pa')

[26]:

print(psf_sersic_model.bounds)

{'amplitude': (0.0, None), 'r_eff': (0.0, None), 'n': (0.0, 10.0), 'x_0': (58.0, 62.0), 'y_0': (58.0, 62.0), 'ellip': (0.0, 1.0), 'theta': (-6.283185307179586, 6.283185307179586), 'psf_pa': (None, None)}

PSF Rotation

PSFConvolvedModel2D can to rotate the PSF image until an optimal rotation angle is found. This is useful for when the PSF comes from a dataset where the orientation of the diffraction spikes are not the same as the image being fit. psf_pa is in degrees.

To restrict the bounds of the rotation or disable the PSF rotation, you can set the psf_pa to fixed:

[27]:

# Limit PSF rotation to -5 to 5 degrees
psf_sersic_model.bounds['psf_pa'] = (-5, 5)

# To disable the PSF rotation,
# you can set the psf_pa to fixed.
psf_sersic_model.fixed['psf_pa'] = True

Accessing the Underlying Model

At any point, a copy of the input model with the same parameter values as the corresponding PSFConvolvedModel2D can be accessed using the PSFConvolvedModel2D.model attribute:

[28]:

psf_sersic_model.model

[28]:

<Sersic2D(amplitude=0.07638308, r_eff=8.24764953, n=1.40861995, x_0=60., y_0=60., ellip=0.07095717, theta=1.3982183)>

Fitting Models

PetroFit uses the Levenberg-Marquardt algorithm to fit parametrized models. To achieve this, it uses astropy fitting and provides wrappers to fit models to images. One such function is fit_model, which takes any Fittable2DModel model and an image to fit, and returns a fitted copy of the model and the fit_info dictionary. If the image to be fit contains pixels that are set to np.nan, those pixels are ignored by the fitter. The fit_model function also allows us to define parameters, such as maxiter, for the astropy fitter.

Before we fit the image, we compute the weights of each pixel using rms data as follows (please note that weights are optional and set to None by defualt):

[29]:

fitting_rms =  Cutout2D(rms, (image_x_0, image_y_0), cutout_size, mode='partial', copy=True)
fitting_weights = 1 / fitting_rms.data

To fit the galaxy we prepared with the PSFConvolvedModel2D we constructed, we call the fit_model as follows:

[30]:

%%time

from petrofit.modeling import fit_model

fitted_model, fit_info = fit_model(
    fitting_image.data, psf_sersic_model,
    weights=fitting_weights,
    maxiter=10000,
    epsilon=1.4901161193847656e-08,
    acc=1e-09,
)

CPU times: user 261 ms, sys: 79 µs, total: 261 ms
Wall time: 261 ms

That’s it! The retuned fitted_model is a copy of the input model (psf_sersic_model) but with the optimized parameter values. We can inspect the parameters of any astropy model using the print_model_params:

[31]:

from petrofit.modeling import print_model_params

print_model_params(fitted_model)

0.0922  amplitude
7.4013  r_eff
1.5552  n
60.8992 x_0
60.4096 y_0
0.1040  ellip
2.0034  theta
0.0000  psf_pa

Generate Model Image

To generate a model image we use the plot_fit function. The function, given a 2D model and fitted image, converts the model into a model-image we can visualize and manipulate.

[32]:

from petrofit.modeling import plot_fit

plot_fit(fitted_model, fitting_image_unmasked, vmax=vmax, vmin=vmin, figsize=[3*6, 6])
plt.show()

If we want to make a model image that we can subtract from the original data image, we just update the center of the model and generate a model image the same dimensions as the original image.

[33]:

from copy import deepcopy
from petrofit.modeling import model_to_image

# Reposition the Model
# --------------------

# Copy the model before changing the parameters
fitted_model_copy = deepcopy(fitted_model)

# Change the frame of reference from the fitting
# image to the original data coordinates.
# Remember that (x_0, y_0) -> (image_x_0, image_y_0)
fitted_model_copy.x_0 =  image_x_0 + (fitted_model.x_0 - x_0)
fitted_model_copy.y_0 =  image_y_0 + (fitted_model.y_0 - y_0)

# Make Model Image
# ----------------

# Set the size of the model image equal to the original image
full_fitted_image_size = image.data.shape[0]

# Center the model image at the center of the original image
# so the two images cover the same window
full_fitted_image_center = full_fitted_image_size // 2

# Generate a model image from the model
full_fitted_model_image = model_to_image(
    fitted_model_copy,
    full_fitted_image_size,

)

# Plot Model Image
# ----------------

fig, ax = plt.subplots(1, 3, figsize=[24, 12])

ax[0].imshow(image.data, vmin=vmin, vmax=vmax)
ax[0].set_title("Data")

ax[1].imshow(full_fitted_model_image, vmin=vmin, vmax=vmax)
ax[1].set_title("Model")

ax[2].imshow(image.data - full_fitted_model_image, vmin=vmin, vmax=vmax)
ax[2].set_title("Residual")

plt.show()

Estimated vs. Fitted Parameters

In this section we compare the initial estimates that we derived purely from photometry vs. the parameters that were fit to the data. In the printout below, we see the fitted values in the first column and our initial guesses in the second column for each parameter.

[34]:

assumptions = [
    amplitude,
    r_eff,
    n,
    x_0,
    y_0,
    ellip,
    theta,
    0 # psf_pa
]

print("fit\tassum\tparam_name")
for param_name, param_val, assumption in zip(fitted_model.param_names, fitted_model.parameters, assumptions):
    print("{:0.2f}\t{:0.2f}\t{}".format(param_val, assumption, param_name))

fit     assum   param_name
0.09    0.08    amplitude
7.40    8.25    r_eff
1.56    1.41    n
60.90   60.00   x_0
60.41   60.00   y_0
0.10    0.07    ellip
2.00    1.40    theta
0.00    0.00    psf_pa

Galaxy vs. Model Profiles

Now that we have a model image, we can do photometry and compare the Petrosian measurements of the model to that of the original galaxy image. Since we know the exact position of the target in the model image and because there are no intruding sources, we can perform direct photometry using photometry_step function.

[35]:

from petrofit.photometry import photometry_step

# Photomerty
model_flux_arr, model_area_arr, model_error_arr = source_photometry(

    # Inputs
    source, # Source (`photutils.segmentation.catalog.SourceCatalog`)
    full_fitted_model_image, # Image as 2D array
    segm_deblend, # Deblended segmentation map of image
    r_list, # list of aperture radii

    # Options
    cutout_size=max(r_list)*2, # Cutout out size, set to double the max radius
    bkg_sub=True, # Subtract background
    sigma=3, sigma_type='clip', # Fit a 2D plane to pixels within 3 sigma of the mean
    plot=True, vmax=vmax, vmin=vmin, # Show plot with max and min defined above
)

plt.title("Model Image")
plt.xlabel("Pixels")
plt.ylabel("Pixels")
plt.show()

[36]:

# Plot Light profile
plt.plot(r_list, flux_arr, label='Data')
plt.plot(r_list, model_flux_arr, label='Model')

plt.title("Light Profiles")
plt.xlabel("Aperture Radius [Pixels]")
plt.ylabel("Cumulative Flux  [{}]".format(image.unit))
plt.legend()
plt.show()

Petrosian Radii Comparison

We compare the Petrosian profiles and radii of the two images.

[37]:

model_p = Petrosian(r_list, model_area_arr, model_flux_arr)

model_corrected_epsilon = pc.estimate_epsilon(
    model_p.r_half_light,
    model_p.concentration_index()[-1]

)

model_corrected_p = Petrosian(
    r_list, model_area_arr, model_flux_arr,
    epsilon=model_corrected_epsilon
)


# Plot Image
plt.imshow(full_fitted_model_image, vmax=vmax, vmin=vmin)


model_position = (
    fitted_model_copy.x_0.value,
    fitted_model_copy.y_0.value)

model_position = (image_x_0, image_y_0)
model_elong = 1 / (1 - fitted_model_copy.ellip.value)

model_theta = fitted_model_copy.theta.value

model_corrected_p.imshow(
    position=model_position,
    elong=model_elong,
    theta=model_theta
)

plt.title("Model Image")
plt.xlabel("Pixels")
plt.ylabel("Pixels")
plt.legend()
plt.show()

[38]:

corrected_p.plot(plot_r=True, plot_normalized_flux=True)
model_corrected_p.plot(plot_r=True, plot_normalized_flux=True)

plt.gca().get_legend().remove()
plt.show()

[39]:

print("Data r_total_flux = {}".format(corrected_p.r_total_flux))
print("Model r_total_flux= {}".format(model_corrected_p.r_total_flux))

Data r_total_flux = 44.44351782614899
Model r_total_flux= 40.53527364847474

[40]:

print("Data r_half_light = {}".format(corrected_p.r_half_light))
print("Model r_half_light= {}".format(model_corrected_p.r_half_light))

Data r_half_light = 8.247649529905981
Model r_half_light= 8.135627125425085

Total Flux Comparison

Finally we compare the total Petrosian flux of the data vs the model:

[41]:

print("Data Corrected Total Flux = {}".format(corrected_p.total_flux * image.unit))
print("Model Corrected Total Flux = {}".format(model_corrected_p.total_flux * image.unit))

Data Corrected Total Flux = 64.0256148087154 electron / s
Model Corrected Total Flux = 62.82013529411775 electron / s

[42]:

from petrofit.photometry import flux_to_abmag
print("Data Corrected AB mag = {}".format(flux_to_abmag(corrected_p.total_flux, image.header) ))
print("Model Corrected AB mag = {}".format(flux_to_abmag(model_corrected_p.total_flux, image.header) ))

Data Corrected AB mag = 21.752933474743692
Model Corrected AB mag = 21.773570699992433