ProFit: Offset Images

Aaron Robotham

2021-02-04

Load our useful libraries:

library(ProFit)
#> Loading required package: FITSio
library(ProFound)
#> Loading required package: Rcpp
library(magicaxis)
library(Highlander)
library(Rfits)

Set global evaluate (basically TRUE for GitHub version and FALSE for CRAN):

evalglobal = TRUE

Fitting Simple Offset Images

Here we look at co-fitting images with offset positions.

im1 = Rfits_read_image(system.file("extdata/Offset/Zim1.fits",package="ProFit"), ext=2)$imDat
im2 = Rfits_read_image(system.file("extdata/Offset/Zim2.fits",package="ProFit"), ext=2)$imDat

Here we have two similar Z-band frames taken at slightly offset positions of the same primary source. With ProFit it is possible to fit these simultaneously without resorting to stacking.

magimage(im1)

magimage(im2)

First we should extract the PSFs using the profitAllStarDoFit function (this will take about 5 minutes):

PSF1 = profitAllStarDoFit(im1, magzero=30, plot=TRUE)
#> Running Found2Fit
#>     Running ProFound
#> [1] "Summary of used sample:"
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> -54.0382   0.3868   1.0965   1.1857   1.9001  22.9720 
#> [1] "sd / MAD / 1-sig / 2-sig range:"
#> [1] 3.801279 1.115354 1.276093 5.403320
#> [1] "Using 2634 out of 2634"
#> Running Highander
#> Iteration 1
#> CMA 1: -6193.77 25.678 76.577 26.219 76.354 25.602 25.943 76.763 76.618 15.185 15.398 17.557 17.601 0.443 0.285
#> 
#> Laplace's Demon was called on Thu Feb  4 15:11:14 2021
#> 
#> Performing initial checks...
#> WARNING: Matrix mask may be rank-deficient.
#> WARNING: Matrix sigma may be rank-deficient.
#> WARNING: Matrix region may be rank-deficient.
#> WARNING: Matrix calcregion may be rank-deficient.
#> Algorithm: Componentwise Hit-And-Run Metropolis 
#> 
#> Laplace's Demon is beginning to update...
#> Iteration: 100,   Proposal: Componentwise,   LP: -4444.5
#> 
#> Assessing Stationarity
#> Assessing Thinning and ESS
#> Creating Summaries
#> Creating Output
#> 
#> Laplace's Demon has finished.
#> LD Mode 1: -4434.757 25.701 76.614 25.973 76.724 25.581 25.946 76.287 76.604 15.181 15.372 17.372 17.398 0.447 0.282
#> Iteration 2
#> 
#> Laplace's Demon was called on Thu Feb  4 15:12:48 2021
#> 
#> Performing initial checks...
#> WARNING: Matrix mask may be rank-deficient.
#> WARNING: Matrix sigma may be rank-deficient.
#> WARNING: Matrix region may be rank-deficient.
#> WARNING: Matrix calcregion may be rank-deficient.
#> Algorithm: Componentwise Hit-And-Run Metropolis 
#> 
#> Laplace's Demon is beginning to update...
#> Iteration: 100,   Proposal: Componentwise,   LP: -4437.7
#> 
#> Assessing Stationarity
#> Assessing Thinning and ESS
#> Creating Summaries
#> Creating Output
#> 
#> Laplace's Demon has finished.
#> LD Mode 2: -4431.386 25.696 76.614 25.98 76.728 25.576 25.951 76.287 76.619 15.177 15.368 17.368 17.397 0.446 0.28

#> [1] "Summary of used sample:"
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#> -7.6589 -0.5816  0.1373  0.1075  0.8405  6.7160 
#> [1] "sd / MAD / 1-sig / 2-sig range:"
#> [1] 1.231222 1.053084 1.071503 2.452146
#> [1] "Using 2634 out of 2634"

PSF2 = profitAllStarDoFit(im2, magzero=30, plot=TRUE)
#> Running Found2Fit
#>     Running ProFound
#> [1] "Summary of used sample:"
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> -39.6172   0.3709   1.1509   1.1267   2.0631  14.6307 
#> [1] "sd / MAD / 1-sig / 2-sig range:"
#> [1] 3.145290 1.230457 1.388358 5.352191
#> [1] "Using 1829 out of 1829"
#> Running Highander
#> Iteration 1
#> CMA 1: -3788.664 25.586 76.287 25.472 76.25 25.397 25.064 76.224 76.739 15.4 16.899 17.193 17.296 0.365 0.266
#> 
#> Laplace's Demon was called on Thu Feb  4 15:14:26 2021
#> 
#> Performing initial checks...
#> WARNING: Matrix mask may be rank-deficient.
#> WARNING: Matrix sigma may be rank-deficient.
#> WARNING: Matrix region may be rank-deficient.
#> WARNING: Matrix calcregion may be rank-deficient.
#> Algorithm: Componentwise Hit-And-Run Metropolis 
#> 
#> Laplace's Demon is beginning to update...
#> Iteration: 100,   Proposal: Componentwise,   LP: -2987.5
#> 
#> Assessing Stationarity
#> Assessing Thinning and ESS
#> Creating Summaries
#> Creating Output
#> 
#> Laplace's Demon has finished.
#> LD Mode 1: -2987.546 25.564 76.29 25.435 76.146 25.443 25.079 76.358 76.982 15.369 16.934 17.17 17.366 0.359 0.257
#> Iteration 2
#> 
#> Laplace's Demon was called on Thu Feb  4 15:16:18 2021
#> 
#> Performing initial checks...
#> WARNING: Matrix mask may be rank-deficient.
#> WARNING: Matrix sigma may be rank-deficient.
#> WARNING: Matrix region may be rank-deficient.
#> WARNING: Matrix calcregion may be rank-deficient.
#> Algorithm: Componentwise Hit-And-Run Metropolis 
#> 
#> Laplace's Demon is beginning to update...
#> Iteration: 100,   Proposal: Componentwise,   LP: -2984.1
#> 
#> Assessing Stationarity
#> Assessing Thinning and ESS
#> Creating Summaries
#> Creating Output
#> 
#> Laplace's Demon has finished.
#> LD Mode 2: -2983.492 25.566 76.274 25.418 76.151 25.452 25.079 76.343 76.977 15.369 16.938 17.169 17.374 0.36 0.255

#> [1] "Summary of used sample:"
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> -7.29631 -0.62025  0.10776  0.08149  0.80825  6.50246 
#> [1] "sd / MAD / 1-sig / 2-sig range:"
#> [1] 1.191087 1.063010 1.086887 2.372173
#> [1] "Using 1829 out of 1829"

We can compare the models:

PSF1$psf_modellist
#> $moffat
#> $moffat$xcen
#> [1] 25.5
#> 
#> $moffat$ycen
#> [1] 25.5
#> 
#> $moffat$mag
#> [1] 0
#> 
#> $moffat$fwhm
#> [1] 2.791671
#> 
#> $moffat$con
#> [1] 1.904349
#> 
#> $moffat$ang
#> [1] 0
#> 
#> $moffat$axrat
#> [1] 1
PSF2$psf_modellist
#> $moffat
#> $moffat$xcen
#> [1] 25.5
#> 
#> $moffat$ycen
#> [1] 25.5
#> 
#> $moffat$mag
#> [1] 0
#> 
#> $moffat$fwhm
#> [1] 2.291604
#> 
#> $moffat$con
#> [1] 1.79823
#> 
#> $moffat$ang
#> [1] 0
#> 
#> $moffat$axrat
#> [1] 1

The seeing in our second image is clearly quite a bit better (FWHM 2.3 versus 2.8). This is also why we should get a better solution by fitting both images simultaneously rather than stacking.

To get most of the way to the correct Datalist structure we need we can use profitFound2Fit:

F2F1 = profitFound2Fit(im1, loc=c(450,500), psf=PSF1$psf, cutbox=c(300,300), Ncomp=2,
                       magzero=30, SBdilate=2, tightcrop=FALSE)
#>     Running ProFound
F2F2 = profitFound2Fit(im2, loc=c(450,500), psf=PSF2$psf, cutbox=c(300,300), Ncomp=2,
                       magzero=30, SBdilate=2, tightcrop=FALSE)
#>     Running ProFound

The main adjustment we need to make it to tell ProFit what the offset is between the frames and put this into the Data object contained. From the original WCSs the offset from im1 to im2 is know to be [-23.112,-9.697] pixels:

F2F2$Data$offset = c(-23.112,-9.697)

We can now check this using the profitLikeModel function, where the same initial parameters get shifted for the second Data object:

profitLikeModel(F2F1$Data$init, Data=F2F1$Data, makeplots=TRUE)
#> [1] "Summary of used sample:"
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> -24.0255  -0.4638   0.3144   0.4133   1.1256  34.2776 
#> [1] "sd / MAD / 1-sig / 2-sig range:"
#> [1] 1.954969 1.176549 1.233549 3.011684
#> [1] "Using 22881 out of 22881"

profitLikeModel(F2F1$Data$init, Data=F2F2$Data, makeplots=TRUE)
#> [1] "Summary of used sample:"
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> -28.9496  -0.2437   0.4894   0.5815   1.3033  36.7655 
#> [1] "sd / MAD / 1-sig / 2-sig range:"
#> [1] 2.028850 1.144899 1.189113 3.014019
#> [1] "Using 17968 out of 17968"

The solutions look pretty similar now! To simplify things we want to create a Datalist structure that contains the two lists:

Datalist = c(list(F2F1$Data), list(F2F2$Data))

You can now run profitLikeModel more directly:

profitLikeModel(F2F1$Data$init, Data=Datalist, makeplots=TRUE)
#> [1] "Summary of used sample:"
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> -24.0255  -0.4638   0.3144   0.4133   1.1256  34.2776 
#> [1] "sd / MAD / 1-sig / 2-sig range:"
#> [1] 1.954969 1.176549 1.233549 3.011684
#> [1] "Using 22881 out of 22881"