Lens Model and Source Reconstruction Reveal the Morphology and Star Formation Distribution in the Cool Spiral LIRG SGAS J143845.1+145407

Observations of luminous infrared galaxies (LIRGs, galaxies with total infrared luminosity (8–1000 μm) of 10^11–12 L_⊙ and a star formation rate (SFR) of ∼10–100 M_⊙ yr⁻¹), which dominate star formation at z ∼ 1, indicate that they evolve strongly with redshift.

Spectroscopy and spectral energy distribution (SED) measurements of LIRGs indicate that at z = 0 they are more likely to have warmer SEDs than their z ∼ 1 counterparts (Rowan-Robinson et al. 2004), and that the temperature of the dust-reprocessed infrared radiation gets colder with increasing redshift (Rowan-Robinson et al. 2005; Sajina et al. 2006; Symeonidis et al. 2009; Elbaz et al. 2010; Hwang et al. 2010). Several studies (e.g., Papovich et al. 2007; Farrah et al. 2008; Rigby et al. 2008; Menéndez-Delmestre et al. 2009) show that the strength of the polycyclic aromatic hydrocarbon features increase with redshift for galaxies with constant luminosity. Morphologically, radio and submillimeter interferometry show that the spatial extent of star formation in LIRGs and ULIRGs (ultraluminous infrared galaxies—those with total infrared luminosity >10¹² L_⊙ and a SFR exceeding ∼100 M_⊙ yr⁻¹) is considerably larger at 0.4 < z < 2.5 than in the local universe (Rujopakarn et al. 2011).

These multiple lines of evidence therefore indicate that the mode of star formation in LIRGs evolves significantly with redshift: at z = 0, star formation occurs on sub-galactic scales of hundreds of parsecs at most, resulting in hot infrared SEDs and weak aromatic features, while star formation at higher redshifts is dominated by low SFRs in galaxy-wide bursts, resulting in cooler SEDs and stronger aromatic features.

However, existing direct imaging studies lack the resolution required to distinguish between these two star formation morphologies at high redshift. At z > 0.4, galaxies are unresolved by Spitzer and Herschel, for example, and even with interferometry in the radio and submillimeter, a beam size of 03–05 is only a few times smaller than the size of the sources, making it impossible to robustly test whether the distribution of star formation is different in z ∼ 1 LIRGs than it is in those at z ∼ 0.

This spatial resolution limitation can be overcome if we are able to take advantage of the magnifying power of gravitational lensing. When paired with the already high resolution of telescopes like the Hubble Space Telescope (HST), gravitational lensing enables observations on scales 10–100 times smaller than otherwise possible, providing insights into the detailed physical properties of high-redshift galaxies (e.g., Wuyts et al. 2012; Bayliss et al. 2014; Bordoloi et al. 2016). In fact, using gravitational lensing, Johnson et al. (2017) recently measured the sizes of individual star-forming regions in a galaxy at z = 2.49 and found some as small as 40 pc across, demonstrating the importance of this magnification boost for studying the details of star formation in high-redshift galaxies. This means that in order to directly determine whether the mode of star formation in LIRGs really does evolve between z = 1 and z = 0, the ideal target to observe would be a strongly lensed galaxy at z ∼ 1 with properties typical of LIRGs at that redshift.

SGAS J143845.1+145407 is such an object. First reported by Gladders et al. (2013), SGAS J143845.1+145407 (hereafter SGAS J1438) is a bright lensed galaxy behind the lensing cluster SDSS J1438+1454 (Figure 1). Using photometry in 18 bands from the blue optical to 500 μm, and optical and near-infrared spectroscopy, these authors determined that SGAS J1438 is a cool LIRG undergoing dusty star formation. Other than being highly magnified, SGAS J1438 is intrinsically a typical LIRG at z = 0.816.

Zoom In Zoom Out Reset image size

We take advantage of this unique opportunity to study the details of star formation in a prototypical high-redshift LIRG. We present new HST imaging data of SDSS J1438+1454 in Section 2 and a lens model based on new constraints from these data in Section 3. We discuss the implications of this lens model in Section 4. Using the lens model we produce a source plane reconstruction of SGAS J1438 in Section 5, which we pair with HST grism spectroscopy data in Section 6 to constrain the spatial extent of star formation in this z = 0.816 LIRG. This allows us to directly test the current understanding of the evolution of star formation in LIRGs from z = 1 to z = 0.

Where necessary, we assume a flat cosmology with Ω_Λ = 0.7, Ω_m = 0.3, and H₀ = 70 km s⁻¹ Mpc⁻¹. In this cosmology 1'' corresponds to 3.758 kpc at the cluster redshift, z = 0.237, and 7.559 kpc at the source redshift, z = 0.816. Magnitudes are reported in the AB system.

SDSS J1438+1454 was observed by HST Cycle 21 program GO13437 (PI: Rigby) during three orbits. Imaging with the Advanced Camera for Surveys (ACS) was executed over one orbit on GMT 2014 March 5, in F814W (600 s) and F606W (780 s). Wide Field Camera 3 (WFC3) grism spectroscopy was conducted in two different roll angles, one on GMT 2014 March 5 and one on GMT 2014 May 30, using the G141 grism (2406 s each). F140W imaging frame was taken at each roll angle (285 s each), giving a total imaging time in this band of 570 s. The ACS images were taken with a three-point line dither using a 1000 × 1000 pixel subraster to manage buffer dumps, resulting in three frames per filter. Since the charge transfer efficiency (CTE) of the ACS detector is decreasing, post-observation image corrections were applied to individual exposures using the pixel-based empirical CTE correction software provided by the Space Telescope Science Institute (STScI). Individual frames were then combined using AstroDrizzle (Gonzaga et al. 2012a, 2012b) with a pixel scale of 003 pixel⁻¹, and a drop size of 0.5 (WFC3) and 0.8 (ACS), following Sharon et al. (2014). All the images were aligned and registered to the same pixel frame as the F140W image.

The slitless spectroscopy data were reduced using the Grizli pipeline.¹¹ Each visit was processed separately rather than being stacked because the data were taken with different roll angles. The two-dimensional spectra of each of the four images of SGAS J1438 are significantly contaminated by light from the brightest cluster galaxy (BCG) and other bright galaxies. To minimize the effects of the light from these objects we deviate from the standard Grizli procedure by modeling the contamination from the galaxies that directly affect the spectra of SGAS J1438 separately from the rest of the field, using GALFIT (Peng et al. 2002, 2010). We find that for the purpose of assigning light to individual objects in areas where light from multiple galaxies overlaps, the GALFIT modeling is better than the Source Extractor segmentation maps that are part of the standard Grizli pipeline. Because of this we use the following alternative method (G. Brammer 2019, private communication): we truncate the GALFIT models when the flux from an object falls to less than 0.001 of the flux of the central pixel. The GALFIT models were then subtracted from the direct images, and the initial contamination models were computed for the remaining objects. The initial contamination models for the two sets of galaxies were combined before the final refinement of the models was conducted. This two-part modeling process resulted in two-dimensional spectral extractions with substantially less contamination from bright nearby galaxies than the standard Grizli pipeline. However, as we will discuss in Section 6, contamination was still significant.

To analyze the images of the lensed galaxy we model the light of the foreground BCG with a multi-component Sérsic profile using GALFIT, and subtract the model from the imaging data in each band. The result (Figure 2) reveals the images of the lensed galaxy that are buried in the light of the foreground BCG.

Zoom In Zoom Out Reset image size

The high-resolution HST data confirm the lensing interpretation of Gladders et al. (2013), which was based on lower-resolution imaging in 18 bands, spanning 0.5–500 μm, from Gemini, Magellan, Spitzer Space Telescope, Wide-field Infrared Survey Explorer, and the Herschel Space Observatory. The strong lensing potential of the foreground cluster causes the appearance of four images of SGAS J1438. Two of these are complete images appearing to the north and south of the BCG (labeled 1 and 3 in Figure 2, respectively); and two are partial images appearing west and east of the BCG (labeled 2 and 4, respectively). In each of the images we identify a number of distinct emission knots, most prominent in F606W and F814W.

We matched the multiply-imaged knots between instances of the lensed galaxy by considering the location of the knots within the galaxy, lensing parity, color, magnification, and brightness. Several emission knots are easily distinguishable, and thus are used to constrain the lens model (Figure 3). Including the bulge of the galaxy, 38 emission knots in total were used as positional constraints. Table 1 lists the IDs and coordinates of the emission knots that were used in the lens model.

Zoom In Zoom Out Reset image size

Table 1.  Lensing Constraints

Image ID	R.A.(J2000)	Decl.(J2000)	rms ('')	$\left|\mu \right|$
10.1	219.68704	14.904572	0.05	${3.2}_{-0.7}^{+0.6}$
10.2	219.68696	14.903430	0.04	${1.8}_{-0.4}^{+1.6}$
10.3	219.68777	14.901929	0.04	${2.4}_{-0.5}^{+0.4}$
11.1	219.68717	14.904282	0.03	${5.4}_{-1.0}^{+2.4}$
11.2	219.68703	14.903652	0.11	${2.8}_{-0.5}^{+2.5}$
11.3	219.68787	14.901701	0.07	${2.1}_{-0.4}^{+0.2}$
12.1	219.68686	14.904191	0.01	${6.4}_{-1.8}^{+2.0}$
12.2	219.68679	14.903699	0.09	${4.5}_{-0.8}^{+2.5}$
12.3	219.68764	14.901701	0.06	${2.2}_{-0.4}^{+0.2}$
13.1	219.68662	14.904108	0.05	${8.0}_{-2.9}^{+1.4}$
13.2	219.68660	14.903708	0.05	${6.4}_{-0.8}^{+3.2}$
13.3	219.68737	14.901716	0.05	${2.3}_{-0.5}^{+0.2}$
14.1	219.68677	14.904353	0.02	${4.2}_{-1.1}^{+0.8}$
14.2	219.68672	14.903530	0.04	${2.9}_{-0.8}^{+1.5}$
14.3	219.68752	14.901806	0.02	${2.4}_{-0.5}^{+0.3}$
15.1	219.68684	14.904579	0.04	${3.1}_{-0.7}^{+0.5}$
15.2	219.68679	14.903349	0.08	${2.1}_{-0.6}^{+1.1}$
15.3	219.68756	14.901952	0.06	${2.6}_{-0.6}^{+0.4}$
16.1	219.68641	14.904516	0.03	${3.1}_{-0.8}^{+0.2}$
16.2	219.68650	14.903197	0.02	${2.9}_{-1.2}^{+0.6}$
16.3	219.68707	14.902029	0.05	${3.1}_{-0.8}^{+0.5}$
17.1	219.68747	14.904875	0.08	${2.7}_{-0.5}^{+0.4}$
17.3	219.68816	14.902269	0.06	${2.9}_{-0.6}^{+0.7}$
17.4	219.68833	14.903578	0.04	${1.7}_{-0.2}^{+1.4}$
18.1	219.68773	14.904724	0.14	${3.3}_{-0.7}^{+0.8}$
18.3	219.68833	14.902189	0.07	${2.6}_{-0.5}^{+0.5}$
18.4	219.68843	14.903754	0.02	${2.1}_{-0.4}^{+1.3}$
19.1	219.68685	14.904857	0.06	${2.5}_{-0.5}^{+0.3}$
19.2	219.68690	14.903173	0.06	${2.1}_{-0.5}^{+1.2}$
19.3	219.68753	14.902201	0.04	${3.2}_{-0.8}^{+0.9}$
20.1	219.68738	14.904508	0.02	${3.8}_{-0.7}^{+1.3}$
20.3	219.68804	14.901883	0.02	${2.3}_{-0.5}^{+0.4}$
21.1	219.68785	14.904465	0.08	${5.6}_{-1.1}^{+2.7}$
21.4	219.68824	14.904039	0.06	${4.6}_{-1.1}^{+1.7}$
22.1	219.68731	14.904751	0.03	${2.9}_{-0.6}^{+0.6}$
22.4	219.68819	14.903645	0.02	${1.4}_{-0.1}^{+1.7}$
23.1	219.68716	14.904693	0.18	${3.0}_{-0.6}^{+0.5}$
23.4	219.68796	14.903606	0.16	${1.6}_{-0.2}^{+4.7}$

Note. The lensing constraints used in the model of SGAS J1438; see Figure 3. The image plane rms is given in arcseconds for each image, and their magnifications are listed for the best-fit model. The 1σ uncertainties are computed from the steps in the MCMC sampling.

Download table as:  ASCIITypeset image

The strong lens model was computed with the public software Lenstool (Jullo et al. 2007), which uses a Bayesian Markov Chain Monte Carlo (MCMC) method to explore the parameter space. The final minimization was done in the image plane by requiring the smallest scatter between the predicted and observed positions of images of each emission knot. We modeled the lens plane as a linear combination of several halos, representing the cluster, the BCG, and cluster-member galaxies. Each halo is parameterized as a pseudo-isothermal ellipsoidal mass distribution (PIEMD; Jullo et al. 2007), with the following parameters: position x, y; ellipticity $e\equiv {({a}^{2}-{b}^{2})({a}^{2}+{b}^{2})}^{-1}$, where a and b are the semimajor and semiminor axes, respectively; position angle θ; core and cut radii r_core and r_cut, respectively; and velocity dispersion σ. All the parameters of the cluster halo were allowed to vary, with the exception of the cut radius, which cannot be constrained by the strong lensing evidence and was therefore fixed at r_cut = 1000 kpc. We verify that this assumption does not affect the results of this work. In addition to the cluster halo we placed a halo that represents the central galaxy; its position was fixed on its observed coordinates, its cut radius was fixed at r_cut = 20 kpc, and all other parameters were left free.

Cluster-member galaxies were identified from the HST photometry by their F814W–F140W color in a color–magnitude diagram (Gladders & Yee 2000). The positions, ellipticities, and position angles of each galaxy were fixed to their observed values, and the PIEMD profile parameters were scaled to their F140W luminosity using scaling relations, following Limousin et al. (2005).

Our model is in general agreement with that of Gladders et al. (2013). In particular, we find that the position angle, core radius, and velocity dispersions of the cluster and of the BCG are within the statistical uncertainties of the respective models. We note, however, that due to the smaller number of lensing constraints that were available to Gladders et al. (2013), they had simplified their lens model by stronger assumptions on the extent to which the mass distribution follows the light distribution. They limit the center of the cluster to not be further than $8\buildrel{\prime\prime}\over{.} {01}_{-3.90}^{+2.01}$ from the BCG, in agreement with our findings of 488. The best-fit model has an image plane rms scatter of 007. Table 1 lists the rms of each image.

The best-fit model is obtained by an MCMC sampling of the parameter space, with a total of 10 free parameters (6 for the cluster halo and 4 for the BCG). Table 2 shows the best-fit parameters and their 1σ uncertainties, as obtained from the MCMC sampling.

Table 2.  Lensing Constraints

Parameter	Cluster	BCG
R.A. ['']	${4.34}_{-2.86}^{+1.06}$	[0]
Decl. ['']	$-{1.08}_{-0.34}^{+0.73}$	[0]
e	${0.87}_{-0.13}^{+0.07}$	${0.47}_{-0.36}^{+0.08}$
θ [°]	${165}_{-1}^{+2}$	${168}_{-29}^{+28}$
rcore [kpc]	${31}_{-31}^{+39}$	${1.58}_{-0.67}^{+2.87}$
σ [km s−1]	${606}_{-96}^{+3}$	${281}_{-51}^{+39}$
rcut [kpc]	[1000]	[20]

Note. Best-fit values from the image plane optimization. The coordinates of the BCG are fixed at their observed location, (R.A., decl.) = (219.68768, 14.903482). The position of the cluster halo is measured relative to the BCG. The position angle θ is measured north of west, and the ellipticity of the projected mass density is e = (a² − b²)/(a² + b²), where a and b are the semimajor and semiminor axes, respectively. The core radius of the cluster was not well constrained by the model, so the uncertainties quoted here are the priors set for the MCMC optimization.

Download table as:  ASCIITypeset image

The critical curves and magnification contours of the best-fit model for a source at z = 0.816 are shown in Figure 4.

Zoom In Zoom Out Reset image size

4.1. The Cluster Mass

4.2. Magnification

Figure 4 shows the magnification contours from the best-fit model. The magnification changes by about a factor of three across each image. The total magnification, ${11.8}_{-2.4}^{+4.6}$, is computed as follows. We first determine the surface brightness of the background by choosing a patch of the image that is mostly sky, and compute the brightness distribution of pixels in this region. We then choose an isophote that is 5σ brighter than the background—this defines a set of polygons in the image plane. We calculate the areas of these polygons using Green's theorem. The polygons are then ray-traced to the source plane using the lensing equation and the deflection matrix from the lens model (see Section 5). The areas are re-calculated in the source plane, again using Green's theorem, and the magnification is then computed as the ratio of the total area in the image plane to the total area in the source plane. The procedure was repeated for 1000 models from the MCMC steps in order to derive the magnification uncertainties (see Figure 6). This procedure results in a non-weighted average magnification of the source galaxy.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The magnification of individual emission knots is less affected by the magnification gradient, given their size. We list the magnification from the best-fit model at the position of each of the emission clumps in Table 1; the uncertainties are calculated from the 1000 MCMC models.

Following the procedure in Sharon et al. (2012, 2014) we ray-trace the image plane pixels from each image of the source through the lensing equation, ${\boldsymbol{\beta }}={\boldsymbol{\theta }}-{d}_{{ls}}/{d}_{s}\,{\boldsymbol{\alpha }}({\boldsymbol{\theta }})$, where ${\boldsymbol{\beta }}$ and ${\boldsymbol{\theta }}$ are the source and image positions, respectively, d_ls/d_s is the ratio of the angular diameter distances from the lens to the source and from the observer to the source, and ${\boldsymbol{\alpha }}({\boldsymbol{\theta }})$ is the deflection matrix of the best-fit model. The coordinates and shape of each pixel are translated to the source plane while conserving surface brightness. Figure 7 shows the reconstruction of the source plane from each one of the four images of SGAS J1438. Since images 2 and 4 are partial images they do not form a complete image of the source galaxy in the source plane, unlike images 1 and 3, which do. In the source plane, the positions of the ray-traced emission knots have a mean scatter of 0033, indicating agreement between the source reconstructions of the different images.

Zoom In Zoom Out Reset image size

Unlensed, the spatial extent of the galaxy in the source plane is ∼5 square arcseconds, which translates to ∼287 kpc² at the source redshift, z = 0.816. The spatial extent of the source is defined by taking the area of the images, ray-tracing them back to the source plane, and re-calculating the area in the source plane. The source plane morphology is similar to the lensed morphology of the galaxy, indicating that the lensing potential does not introduce significant distortion. The lensing magnification is generally a combination of an isotropic magnification component due to the local projected mass density, and an anisotropic component due to shear. Low distortion indicates that in this location the lens produces little shear, resulting in nearly isotropic magnification.

The source plane reconstructions (Figure 7) reveal SGAS J1438 to be a large, two-arm, grand-design spiral galaxy. This was suggested by the image plane morphology but is even more evident given the source plane reconstructions. The red color of the core of this galaxy and the shape of the SED published in Gladders et al. (2013) suggest that the nuclear emission is extremely obscured by dust.

In contrast, the spiral arms contain many distinct patches of emission with colors much bluer than the core, suggesting that there is widespread recent star formation in the arms that is largely unobscured by dust. At z = 0.816, Hα emission, a tracer of ongoing star formation, falls at a wavelength that does not lie within any of the three broadband filters used to observe SGAS J1438. It is, however, within the wavelength range observable in our HST WFC3/IR grism spectroscopy.

Hα emission is not detected in the nucleus of SGAS J1438 in our grism observations. Furthermore, because of the orientation of the galaxy the brightest emission in the spiral arms of images 1 and 3 is closely aligned with the much stronger continuum emission from the nucleus in the direction of dispersion, causing the nuclear emission to drown out the signal in the spiral arms. Unfortunately, grism observations at other angles would have caused the BCG contamination to overwhelm the spectra of images 1 and 3. That same contamination, at the angles used in our grism program, makes it impossible to access the spectra of images 2 and 4. However, in image 3, the northwestern spiral arm is successfully separated from the nucleus, albeit with moderate to high contamination from the BCG. A hint of Hα emission is visible, as shown in Figure 8. Comparing an isophote of the direct image (left panel of Figure 8) with the shape of the faint emission in the 2D grism spectrum at the wavelength of Hα at z = 0.816 (center panel) shows very faint Hα emission coincident with the large blue clumps in the spiral arm seen on the right of panel 3 in Figure 7. A suggestion of Hα emission from the southeastern portion of image 3 is also visible above the spectrum of the nucleus, but it is even less clear.

Zoom In Zoom Out Reset image size

Because of the contamination from the BCG it is difficult to estimate the strength of the Hα emission in the spiral arms. The fact that the Hα emission aligns spatially with one of the blue regions in the broadband imaging bolsters the interpretation that it and the many other similarly colored regions throughout the spiral arms are regions of ongoing star formation. Unfortunately, with these data it is difficult to comment on the star formation in the nucleus because it is so severely obscured by dust. Interestingly, Gladders et al. (2013) reported a detection of Hα in the nucleus, which we fail to detect here. With HST-quality imaging, however, we see that there are blue clumps—likely Hα emitters—to the south and southeast of the nucleus and hints of others to the north and northwest that were included in the slit of the Gladders et al. (2013) spectroscopy. We now believe that the nuclear Hα emission reported in that paper is actually emission from these blue clumps that were not identifiable in extant ground-based imaging at the time.

The finding that star-forming regions are large (the largest such regions are on the kiloparsec scale, according to the reconstructed source image) and distributed throughout the galaxy with a morphology consistent with a grand-design spiral lends support to the interpretation that star formation occurs on galaxy-wide scales in z ∼ 1 LIRGs rather than in small (hundreds of parsecs or smaller), localized regions like in z ∼ 0 LIRGs. One caveat to this apparent galaxy-scale star formation is the deficit of observed Hα emission in the core, which probably is not a deficit of star formation but likely reflects enhanced obscuration at that location. Unfortunately, the extant long-wavelength imaging on SGAS J1438 lacks the spatial resolution to clearly demonstrate extensive obscured star formation in the core. Understanding whether or not the core is actively forming stars would be useful for understanding what physical conditions lead to this galaxy-wide mode of star formation; but in either case we have directly shown that the star formation in this z ∼ 1 LIRG is indeed occurring at much larger spatial scales than in z ∼ 0 LIRGs, as was previously expected but not directly observed.

We find that the mass of the lens is consistent with a small group of galaxies (e.g., Han et al. 2015), in line with the cluster richness reported in Gladders et al. (2013) of ${N}_{\mathrm{gals}}^{\mathrm{weighted}}=9.734$. The lens model favors a displacement of the center of the cluster halo from the BCG, at 488 in the direction of the second and third brightest galaxies in this structure. A small offset between the BCG and the center of mass of the group or cluster is expected (e.g., George et al. 2012), and consistent with the light distribution of the cluster (Figure 1). Our model is in good agreement with the observed velocity dispersion of the BCG reported in Gladders et al. (2013), σ_v = 318 ± 111 km s⁻¹.

The magnification of SGAS J1438 allows us to study star formation in the source galaxy in detail that would be unattainable without lensing. The typical lensing magnification of individual emission knots, a factor of 3–5, enables studies of individual star-forming regions in a galaxy at z = 0.816 with a high signal-to-noise ratio. The angular size of the source galaxy, unmagnified, is 253 (≈19 kpc). In good conditions, ground-based resolution of 06 would allow at best up to four non-overlapping resolution elements, dominated by the bulge. With the lensing magnification boost even ground-based observations can resolve more than ten regions in this galaxy and thus resolve its two-arm structure.

HST imaging has revealed the structure of this galaxy to be a large, two-arm, grand-design spiral with a very red core and bluer clumps with sizes up to the kiloparsec scale scattered through its spiral arms. HST/WFC3 grism spectroscopy of all the images of the galaxy were severely contaminated and dominated by the light of the foreground BCG; high-resolution integral field spectroscopy would overcome this issue. Nevertheless, the grism data confirm that Hα emission is present in these clumps, though no such emission is detected in the core. The widespread star formation in the spiral arms supports previous interpretations of the spatially unresolved spectra of z ∼ 1 LIRGs, namely that such objects are likely to have large-scale, nearly galaxy-wide star formation, in contrast to their counterparts in the local universe where star formation can be intense but occurs on smaller scales.

This work is based on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with program #14420.

We thank the anonymous referee for their constructive feedback that significantly improved this paper. Support for program GO-14420 was provided by NASA through a grant from the Space Telescope Science Institute. This work makes use of the Matlab Astronomy Package (Ofek 2014).