profitEllipsePlot: Plot Isophotal Surface Brightness for Pseudo-Ellipses

Description

In the world of galaxy fitting, projected 1D flux intensity (or surface brightness) plots are popular. This function implements the high level functionality of deprojecting image pixels given a set of geometrical parameters and making a 1D surface brightness plot. In simple terms this means ellipses are expanded back up to circles. We use the term pseudo-ellipse since we can also account for boxiness distortion (if desired).

Usage

profitEllipsePlot(Data, modellist, bulgeloc = 1, diskloc = 2, pixscale = 1, FWHM = 1,
SBlim = 26, df = 100, raw = FALSE)

Arguments

Data

Data of class profit.data. This must be generated by the profitSetupData function.

modellist

Model list (see profitMakeModel for details). To aid interpreting the results of fits the user can access profitRemakeModellist, where an initial model is adjusted to include new parameter values and is remade into a new model.

bulgeloc

Location ID of bulge component in the Sersic list provided in modellist.

diskloc

Location ID of disk component in the Sersic list provided in modellist.

pixscale

The pixel scale, where pixscale=asec/pix (e.g. 0.4 for SDSS). If set to 1, then the plot is output in terms of pixels, otherwise it is rescaled to be in arcseconds.

FWHM

The full width half max of the PSF in units of arc seconds. A vertical line is drawn at half this number (since we are plotting radius). The fits inside of the region is inherently hard because it is well within the PSF convolution kernel.

SBlim

5 sigma surface brightness limit of the data. The data will plot to SBlim+1. This should be provided in terms of how the pixscale is defined, e.g. it should either be per asec^2 (when pixscale!=1) or per pix^2 (when pixscale=1).

Degrees of freedom to use for spline fitting. Lower if the lines look too wavy.

raw

Logical; if FALSE (the default) then a smooth spline is used to represent the data and model 1D profiles. This smooths out deprojection noise caused by the PSF often being non-smooth. If TRUE then the raw pixel surface brightness values are shown. These will show much more scatter, but the trends ought to be very similar. If the raw and smooth 1D plots differ significantly then the df flag probably needs to be changed to improve the smoothing. Notice that when the raw pixel values are plotted the shaded error polygon is very hard to see (it is usually subdominant compared to the pixel scatter created during deprojection that has both the Normal pixel error and the PSF induced deprojection error).

Value

Run for the side effect of making a nice surface brightness and surface brightness residual plot. Most of the plot features are automatic.

For reference, a scaled version of the PSF profile is plotted in purple.

Vertical dashed lines are drawn at the HWHM of the PSF and at the point where the model profile crosses the SBlim value provided. Also, a horizontal dashed line is drawn at the SBlim value provided.

Details

A word of warning: 1D surface brightness plots hide many sins, and interpreting is dark and non-rigorous art. The manner of the lossy information compression the data and model undergo is such that you should not draw strong conclusions about features in the residuals.

The pixel flux is not corrected for distortion. This is rarely an issue, but for highly elliptical galaxies where the gradient changes radically through the pixel along the minor versus major axis the isophotal contour will be biased. In practice this bias will be the same as for the model, so pixel to pixel comparisons are still fairly valid. If this level of inconsistency annoys you, then you almost certainly should not be attempting to make deprojected 1D plots of PSF convolved bulge-disk systems. This is because once a pure elliptical disk has been convolved with a PSF it is not strictly possible to project this to circular annuli (i.e. the pseudo major-axis output we desire) using elliptical isophotes. C'est la vie.

Examples

Run this code

# NOT RUN {
# Here we use galaxy G266033:

image = readFITS(system.file("extdata", 'KiDS/G266033fitim.fits',package="ProFit"))$imDat
sigma = readFITS(system.file("extdata", 'KiDS/G266033sigma.fits',package="ProFit"))$imDat
segim = readFITS(system.file("extdata", 'KiDS/G266033segim.fits',package="ProFit"))$imDat
psf = readFITS(system.file("extdata", 'KiDS/G266033psfim.fits',package="ProFit"))$imDat

#The rough best-fit model for G266033 (KiDS)

modellist=list(
  sersic=list(
    xcen = c(65.60642, 65.60642),
    ycen = c(78.6091, 78.6091),
    mag = c(18.49816, 16.97754),
    re = c(1.69112, 14.75940),
    nser = c(1.053961, 1),
    ang = c(39.53360, 35.02479),
    axrat = c(1, 0.5990869),
    box = c(0,0)
  )
)

# What should we be fitting:

tofit=list(
  sersic=list(
    xcen= c(TRUE,NA), #We fit for xcen and tie the two togther
    ycen= c(TRUE,NA), #We fit for ycen and tie the two togther
    mag= c(TRUE,TRUE), #Fit for both
    re= c(TRUE,TRUE), #Fit for both
    nser= c(TRUE,FALSE), #Fit for bulge
    ang= c(FALSE,TRUE), #Fit for disk
    axrat= c(FALSE,TRUE), #Fit for disk
    box= c(FALSE,FALSE) #Fit for neither
  )
)

# What parameters should be fitted in log space:

tolog=list(
  sersic=list(
    xcen= c(FALSE,FALSE),
    ycen= c(FALSE,FALSE),
    mag= c(FALSE,FALSE),
    re= c(TRUE,TRUE), #re is best fit in log space
    nser= c(TRUE,TRUE), #nser is best fit in log space
    ang= c(FALSE,FALSE),
    axrat= c(TRUE,TRUE), #axrat is best fit in log space
    box= c(FALSE,FALSE)
  )
)

# The priors. If the parameters are to be sampled in log space (above) then the priors
# will refer to dex not linear standard deviations. Priors should be specified in their
# unlogged state- the logging is done internally.

sigmas=c(2,2,2,2,5,5,1,1,1,1,30,30,0.3,0.3,0.3,0.3)

priors=list(
  sersic=list(
    xcen=list(function(x){dnorm(x,0,sigmas[1],log=TRUE)},function(x){dnorm(x,0,sigmas[2],
    log=TRUE)}), # should have tight constraints on x and y
    ycen=list(function(x){dnorm(x,0,sigmas[3],log=TRUE)},function(x){dnorm(x,0,sigmas[4],
    log=TRUE)}), # should have tight constraints on x and y
    mag=list(function(x){dnorm(x,0,sigmas[5],log=TRUE)},function(x){dnorm(x,0,sigmas[6],
    log=TRUE)}), # 5 mag SD
    re=list(function(x){dnorm(x,0,sigmas[7],log=TRUE)},function(x){dnorm(x,0,sigmas[8],
    log=TRUE)}), # i.e. 1 dex in re is the SD
    nser=list(function(x){dnorm(x,0,sigmas[9],log=TRUE)},function(x){dnorm(x,0,sigmas[10],
    log=TRUE)}), # i.e. 1 dex in nser is the SD
    ang=list(function(x){dnorm(x,0,sigmas[11],log=TRUE)},function(x){dnorm(x,0,sigmas[12],
    log=TRUE)}), # very broad 30 deg ang SD
    axrat=list(function(x){dnorm(x,0,sigmas[13],log=TRUE)},function(x){dnorm(x,0,sigmas[14],
    log=TRUE)}), # i.e. 1 dex in axrat is the SD
    box=list(function(x){dnorm(x,0,sigmas[15],log=TRUE)},function(x){dnorm(x,0,sigmas[16],
    log=TRUE)})
  )
)

#the hard intervals should also be specified in log space if relevant:

lowers=c(0,0,0,0,10,10,0,0,-1,-1,-180,-180,-1,-1,-1,-1)
uppers=c(1e3,1e3,1e3,1e3,30,30,2,2,1.3,1.3,360,360,0,0,1,1)

intervals=list(
  sersic=list(
    xcen=list(function(x){interval(x,lowers[1],uppers[1],reflect=FALSE)},
    function(x){interval(x,lowers[2],uppers[2],reflect=FALSE)}),
    ycen=list(function(x){interval(x,lowers[3],uppers[3],reflect=FALSE)},
    function(x){interval(x,lowers[4],uppers[4],reflect=FALSE)}),
    mag=list(function(x){interval(x,lowers[5],uppers[5],reflect=FALSE)},
    function(x){interval(x,lowers[6],uppers[6],reflect=FALSE)}),
    re=list(function(x){interval(x,lowers[7],uppers[7],reflect=FALSE)},
    function(x){interval(x,lowers[8],uppers[8],reflect=FALSE)}),
    nser=list(function(x){interval(x,lowers[9],uppers[9],reflect=FALSE)},
    function(x){interval(x,lowers[10],uppers[10],reflect=FALSE)}),
    ang=list(function(x){interval(x,lowers[11],uppers[11],reflect=FALSE)},
    function(x){interval(x,lowers[12],uppers[12],reflect=FALSE)}),
    axrat=list(function(x){interval(x,lowers[13],uppers[13],reflect=FALSE)},
    function(x){interval(x,lowers[14],uppers[14],reflect=FALSE)}),
    box=list(function(x){interval(x,lowers[15],uppers[15],reflect=FALSE)},
    function(x){interval(x,lowers[16],uppers[16],reflect=FALSE)})
  )
)

# Setup the data structure we need for optimisation:

Data=profitSetupData(image=image, sigma=sigma, segim=segim, psf=psf,
modellist=modellist, tofit=tofit, tolog=tolog, priors=priors, intervals=intervals,
magzero=0, algo.func='optim', verbose=TRUE)

modelimage=profitMakeModel(Data$mode,dim=Data$imagedim)
profitMakePlots(Data$image, modelimage$z, Data$region, Data$sigma)

#The pixel scale of VST/KiDS is 0.2 asec/pix and the SBlim=26 mag/asec^2 in r-band.

profitEllipsePlot(Data, Data$modellist, pixscale=0.2, SBlim=26)

#So to get things into pixel space and pixel surface brightness:

profitEllipsePlot(Data, Data$modellist, pixscale=1, SBlim=26-5*log10(0.2))
# }