# 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))
# }
Run the code above in your browser using DataLab