flux.measurements <-
function(env=NULL) {
#Procedure measures the fluxes present in a supplied image
#inside catalogued apertures.
#PART ZERO: PSF DETERMINATION & GENERAL INITIALISATION /*fold*/ {{{
# Load Parameter Space /*fold*/ {{{
if (!is.null(env)) {
attach(env)
}
# /*fend*/ }}}
#Set function Environments /*fold*/ {{{
environment(make.catalogue.apertures)<-environment()
environment(make.gaussian.apertures)<-environment()
environment(make.segmentation.apertures)<-environment()
environment(make.convolved.apertures)<-environment()
environment(make.aperture.map)<-environment()
environment(make.data.array.maps)<-environment()
environment(read.psf)<-environment()
environment(estimate.psf)<-environment()
environment(make.deblend.weight.maps)<-environment()
environment(write.flux.measurements.table)<-environment()
environment(write.deblend.fraction.table)<-environment()
environment(write.fits.image.file)<-environment()
environment(source.subtraction)<-environment()
environment(get.stamp.cog)<-environment()
for (nam in ls.deb("package:LAMBDAR",simple=TRUE)) {
if (nam%in%ls(envir=environment())) {
debug(get(nam,envir=environment()))
}
}
# /*fend*/ }}}
if (!quiet) { cat('Begining Flux Measurements\n') }
message('} Initialisation Complete\nBegining Flux Measurements')
dimim<-dim(image.env$im)
#Get PSF details /*fold*/ {{{
timer<-proc.time()
if (!(no.psf)) {
#There is a PSF specified /*fold*/ {{{
#Details /*fold*/ {{{
#Calculate PSF - if one has not be supplied,
#then a gaussian PSF will be created. Stampsize of PSF
#should be = maximum stampsize: max aperture * stamp mult /*fend*/ }}}
#Get PSF /*fold*/ {{{
if (psf.map=='ESTIMATE') {
#Estimate the stamplims with conservative guess at PSF FWHM; Nyquist sampled with 10xFHWM width{{{
psffwhm<-3*2.335
stamplen<-30 #(3*2.335*10)+1)
psf.clip<-stamplen
#}}}
#Convert decimal pixel values into actual pixel values /*fold*/ {{{
x.pix<-floor(cat.x)
y.pix<-floor(cat.y)
# /*fend*/ }}}
#Make the first pass data stamps {{{
timer=system.time(data.stamp<-make.data.array.maps(outenv=environment()))
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(paste('Make data array maps - Done (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
#}}}
#Estimate the PSF from the data {{{
timer<-proc.time()
if (!quiet) { cat('Generating the PSF from the data') }
message('Generating the PSF from the data')
if (plot.sample) { pdf(height=10,width=10,file=file.path(path.root,path.work,path.out,'PSFGen_Samples.pdf')) }
psf.est<-estimate.psf(outenv=environment(),plot=plot.sample,radial.tolerance=radial.tolerance,n.sources=n.sources)
psf.id<-tmp.psf.id
psf.val<-tmp.psfest.val
if (exists('tmp.skyest')) {
skyest<-tmp.skyest
}
if (plot.sample & !grepl('x11',plot.device,ignore.case=TRUE)) { dev.off() }
psf.cen<-estpsf<-psf<-list(NULL)
sumpsf<-rep(NA,length(psf.est$WEIGHT))
warn<-FALSE
for (i in 1:length(psf.est$WEIGHT)) {
if (!all(psf.est$WEIGHT[[i]]==0)) {
#If the estimate succeeded {{{
#Normalise the estimate {{{
estpsf[[i]]<-psf.est$PSF[[i]]
estpsf[[i]]<-estpsf[[i]]-min(estpsf[[i]],na.rm=TRUE)
estpsf[[i]]<-estpsf[[i]]/max(estpsf[[i]],na.rm=TRUE)
#}}}
#Truncate any upturn in the PSF {{{
if (plot.sample) {
PlotDev(file=file.path(path.root,path.work,path.out,paste0("PSFGen_truncation_",i,".",plot.device)),width=12,height=6,units='in')
layout(matrix(1:6,nrow=2,byrow=T))
}
psf.cen[[i]]<-psf.est$centre
trunc<-truncate.upturn(estpsf[[i]],plot=plot.sample,centre=psf.est$centre,tolerance=psfest.tolerance)
psf[[i]]<-trunc$Im
psf.cen[[i]]<-trunc$centre
psf[[i]]<-psf[[i]]-min(psf[[i]],na.rm=TRUE)
psf[[i]]<-psf[[i]]/max(psf[[i]],na.rm=TRUE)
sumpsf[i]<-sum(psf[[i]])
if (plot.sample & !grepl('x11',plot.device,ignore.case=TRUE)) { dev.off() }
big<-matrix(0,ncol=max(stamplen),nrow=max(stamplen))
big[1:dim(psf[[i]])[1],1:dim(psf[[i]])[2]]<-psf[[i]]
psf[[i]]<-big
#}}}
#}}}
} else {
warn<-TRUE
#PSF estimate failed; warn if we can continue, or error if we cannot {{{
if (length(psf.est$WEIGHT)==1) {
#All psf determinations failed {{{
if (loop.total>1) {
if (showtime) { cat(" - BREAK (",round(timer[3],digits=2),"sec )\n")
message(paste('Flux Measurements - BREAK (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - BREAK\n") }
warning("PSF estimate failed; no suitable point sources found!")
#Notify & Close Logfile /*fold*/ {{{
if (!is.null(env)) {
on.exit(detach(env), add=TRUE)
}
message(paste('\n-----------------------------------------------------\nDatamap Skipped - No PSF Estimate Possible (no point sources found/usable)\n'))
if (!quiet) {
cat(paste('\n-----------------------------------------------------\nDatamap Skipped - No PSF Estimate Possible (no point sources found/usable)'))
}
return(NULL)
#}}}
} else {
if (showtime) {
cat(" - ERROR: PSF estimate failed; no suitable point sources found! (",round(proc.time()[3]-timer[3],digits=2),"sec )\n")
stop(" - ERROR: PSF estimate failed; no suitable point sources found! (",round(proc.time()[3]-timer[3],digits=2),"sec )\n")
} else if (!quiet) {
cat(" - ERROR: PSF estimate failed; no suitable point sources found!\n")
stop(" - ERROR: PSF estimate failed; no suitable point sources found!\n")
}
}
#}}}
} else {
estpsf[[i]]<-NA
psf[[i]]<-NA
if (i==min(psf.id)) {
if (i==length(psf)) {
#All psfs have failed, leave the id as it is.
} else {
#If this is the lowest psf.id, give
#these sources the id of the next psf up
psf.id[which(psf.id==i)]<-min(psf.id)+1
}
} else {
#If this is not the lowest psf.id, give
#these sources the id of the last psf
psf.id[which(psf.id==i)]<-max(psf.id[which(psf.id<i)])
}
}
#}}}
}
}
if (all(psf.id==psf.id[1])&&is.na(psf[[psf.id[1]]])) {
#All psf determinations failed {{{
if (loop.total>1) {
if (showtime) { cat(" - BREAK (",round(timer[3],digits=2),"sec )\n")
message(paste('Flux Measurements - BREAK (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - BREAK\n") }
warning("PSF estimate failed; no suitable point sources found!")
#Notify & Close Logfile /*fold*/ {{{
if (!is.null(env)) {
on.exit(detach(env), add=TRUE)
}
message(paste('\n-----------------------------------------------------\nDatamap Skipped - No PSF Estimate Possible (no point sources found/usable)\n'))
if (!quiet) {
cat(paste('\n-----------------------------------------------------\nDatamap Skipped - No PSF Estimate Possible (no point sources found/usable)'))
}
return(NULL)
#}}}
} else {
if (showtime) {
cat(" - ERROR: PSF estimate failed; no suitable point sources found! (",round(proc.time()[3]-timer[3],digits=2),"sec )\n")
stop(" - ERROR: PSF estimate failed; no suitable point sources found! (",round(proc.time()[3]-timer[3],digits=2),"sec )\n")
} else if (!quiet) {
cat(" - ERROR: PSF estimate failed; no suitable point sources found!\n")
stop(" - ERROR: PSF estimate failed; no suitable point sources found!\n")
}
}
#}}}
} else if (warn) {
if (showtime) {
cat(" - WARNING: PSF estimate failed in one or more bins; they inherit the PSF of another bin! (",round(proc.time()[3]-timer[3],digits=2),"sec )\n")
message(" WARNING: PSF estimate failed in one or more bins; they inherit the PSF of another bin! (",round(proc.time()[3]-timer[3],digits=2),"sec )\n")
} else if (!quiet) { cat(" - WARNING: PSF estimate failed in one or more bins; they inherit the PSF of another bin!\n") }
} else {
if (showtime) {
cat(" - Done (",round(proc.time()[3]-timer[3],digits=2),"sec )\n")
message(" - Done (",round(proc.time()[3]-timer[3],digits=2),"sec )\n")
} else if (!quiet) { cat(" - Done\n") }
}
#Save the PSF as output {{{
psf.weight<-rep(NA,length(psf))
for (i in 1:length(psf)) {
psf.weight[i]<-length(which(psf.id==i))
}
psflist=list(final=psf,normalised=estpsf,raw=psf.est,psf.weight=psf.weight)
save(file=file.path(path.root,path.work,path.out,paste0(tableout.name,"_PSFGen_estimate.Rdata")),psflist)
#}}}
#}}}
} else {
timer<-proc.time()
if (gauss.fwhm.arcsec<0 & !file.exists(file.path(path.root,path.work,psf.map))) {
if (!quiet) { cat('Reading PSF details from Header') }
message('Reading PSF details from Header')
} else if (file.exists(file.path(path.root,path.work,psf.map))) {
if (!quiet) { cat('Reading PSF from file') }
message('Reading PSF from file')
} else if (gauss.fwhm.arcsec > 0) {
if (!quiet) { cat('Generating PSF from parameter file input') }
message('Generating PSF from parameter file input')
} else {
if (!quiet) { cat('How did we get here?!') }
message('How did we get here?!')
}
#Read/Generate the PSF {{{
psf<-list(read.psf(outenv=environment(),file.path(path.root,path.work,psf.map),arcsec.per.pix,max(cat.a),confidence,gauss.fwhm.arcsec=gauss.fwhm.arcsec))
psf.id<-rep(1,length(cat.x))
psf.weight<-1
#}}}
#Notify /*fold*/ {{{
if (showtime) { cat(paste(' - Done (',round(proc.time()[3]-timer[3], digits=2),'sec )\n'))
message(paste('Getting PSF details - Done (',round(proc.time()[3]-timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(' - Done\n') }
# /*fend*/ }}}
}
# /*fend*/ }}}
#Notify /*fold*/ {{{
if (verbose) {
for (i in 1:length(psf)) {
message(paste("Maxima of the PSF",i,"is at pixel", which(psf[[i]] == max(psf[[i]])),"and has value",max(psf[[i]])))
}
}
# /*fend*/ }}}
#Get radius of FWHM using FWHM confidence value = erf(2*sqrt(2*log(2))/sqrt(2)) /*fold*/ {{{
psffwhm<-sum(unlist(mapply(get.fwhm,psf,centre=psf.cen))*psf.weight,na.rm=T)/sum(psf.weight)
# /*fend*/ }}}
#Normalise Beam Area /*fold*/ {{{
beam.area.nn<-sumpsf
beam.area.n<-sumpsf/sapply(psf,max)
# /*fend*/ }}}
#-----Diagnostic-----## /*fold*/ {{{
if (diagnostic) { message(paste('Beam area before/after norm: ',beam.area.nn,beam.area.n)) }
# /*fend*/ }}}
# /*fend*/ }}}
} else {
#There is no PSF specified /*fold*/ {{{
#set relevant paramters manually
beam.area.n<-1.0
psfwidth<-0
psf.clip<-0
sumpsf<-0
# /*fend*/ }}}
#If no PSF, set FWHM to median aperture radius + buffer /*fold*/ {{{
psffwhm<-round(min(cat.a[which(cat.a>0)])/arcsec.per.pix)
if ((is.na(psffwhm))|(!is.finite(psffwhm))) {
#All Apertures are point sources, force width = 5 pixels /*fold*/ {{{
#Details /*fold*/ {{{
#There are no apertures that are not point sources,
#and no PSF convolution - everything is a single pixel flux
#set to a default of 5pixels /*fend*/ }}}
psffwhm<-5
# /*fend*/ }}}
}
# /*fend*/ }}}
}
#If possible, Update Beamarea /*fold*/ {{{
if (beam.area.pix == 0) {
#Use the beamarea just determined /*fold*/ {{{
beamarea<-beam.area.n
# /*fend*/ }}}
} else {
# Otherwise just use the input beamarea /*fold*/ {{{
beamarea<-beam.area.pix
# /*fend*/ }}}
}# /*fend*/ }}}
# /*fend*/ }}}
#-----Diagnostic-----# /*fold*/ {{{
if (diagnostic) { message(paste('Beam area adopted: ',beamarea)) }
# /*fend*/ }}}
#Convert images from Jy/beam to Jy/pixel /*fold*/ {{{
if (Jybm) {
if (length(beamarea)>1) {
message("WARNING: Jybm to Jypx; Using Weighted Average beam area for conversion!")
beamtmp<-sum(beamarea*psf.weight,na.rm=T)/sum(psf.weight)
} else {
beamtmp<-beamarea
}
image.env$im<-image.env$im/beamtmp
if (length(image.env$ime)>1) { image.env$ime<-image.env$ime/beamtmp }
conf<-conf/beamtmp
rm(beamtmp)
} # /*fend*/ }}}
# /*fend*/ }}}
#PART ONE: SINGLE APERTURE MASKS /*fold*/ {{{
#Details /*fold*/ {{{
#Create apertures for every object remaining in
#the catalogue. /*fend*/ }}}
#Get StampLengths for every galaxy /*fold*/ {{{
if (psf.filt) {
stamplen<-(floor(ceiling(def.buff*cat.a/arcsec.per.pix)+(ceiling(psf.clip)/2))*2+5)
} else {
stamplen<-(floor(ceiling(def.buff*cat.a/arcsec.per.pix))*2+5)
}
# /*fend*/ }}}
#Discard any contaminants that are beyond their nearest neighbours stamp /*fold*/ {{{
#Number of Nearest Neighbors to search /*fold*/ {{{
if (!exists("num.nearest.neighbours")) { num.nearest.neighbours<-10 }
if (!exists("check.contam")) { check.contam<-TRUE }
if (!exists("getLoners")) { getLoners<-FALSE }
# /*fend*/ }}}
if (getLoners) {
timer<-proc.time()
if (!quiet) { cat("Selecting out only objects that are completely alone ") }
cat.len<-length(cat.x)
num.nearest.neighbours<-max(num.nearest.neighbours,cat.len-1)
nearest<-nn2(data.frame(cat.x,cat.y),data.frame(cat.x,cat.y),radius=4*max(psffwhm),k=2,type='radius')
if (psffwhm!=0) {
#If any of the nearest neighbors overlap with the object /*fold*/ {{{
inside.mask<-nearest$nn.id[,2] == 0
} else {
inside.mask<-rep(NA,cat.len)
for(ind in 1:cat.len) {
#If any of the nearest neighbors overlap with the object /*fold*/ {{{
inside.mask[ind]<-any(nearest$nn.dist[ind,1:num.nearest.neighbours+1] < (sqrt(2)/2*(stamplen[nearest$nn.idx[ind,1:num.nearest.neighbours+1]]+stamplen[ind])))
# /*fend*/ }}}
}
# /*fend*/ }}}
}
#Remove object catalogue entries /*fold*/ {{{
cat.x<-cat.x[which(inside.mask)]
cat.y<-cat.y[which(inside.mask)]
cat.id<-cat.id[which(inside.mask)]
cat.ra<-cat.ra[which(inside.mask)]
cat.dec<-cat.dec[which(inside.mask)]
cat.theta<-cat.theta[which(inside.mask)]
cat.a<-cat.a[which(inside.mask)]
cat.b<-cat.b[which(inside.mask)]
if (use.segmentation) { seg.x<-seg.x[which(inside.mask)] }
if (use.segmentation) { seg.y<-seg.y[which(inside.mask)] }
if (length(flux.weight)!=1) { flux.weight<-flux.weight[which(inside.mask)] }
if (exists("contams")) { contams<-contams[which(inside.mask)] }
if (exists("groups")) { groups<-groups[which(inside.mask)] }
if (exists("psf.id")) { psf.id<-psf.id[which(inside.mask)] }
if (num.cores < 0) {
registerDoParallel(cores=(min(ceiling(length(cat.x)/50000),abs(num.cores))))
}
chunk.size=length(cat.id)/getDoParWorkers()
mpi.opts<-list(chunkSize=chunk.size)
message("Number of objects per thread:",chunk.size)
# /*fend*/ }}}
#Notify how many objects remain /*fold*/ {{{
if (verbose) { message(paste("There are",length(cat.x),"supplied objects that are entirely unblended (",
round(((cat.len-length(cat.x))/cat.len)*100, digits=2),"% of objects were objects that were possibly blended)")) }
# /*fend*/ }}}
if (showtime) { cat(" - Done (",round(proc.time()[3]-timer[3],digits=2),"sec )\n")
message(paste('Remove irrelevant contaminants - Done (',round(proc.time()[3]-timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
} else if (check.contam && exists("contams") && length(which(contams==0))>1 && length(which(contams==1))>1) {
timer<-proc.time()
if (!quiet) { cat("Removing Contaminants that are irrelevant ") }
cat.len<-length(cat.x)
num.nearest.neighbours<-min(num.nearest.neighbours,length(which(contams==0))-1)
nearest<-nn2(data.frame(cat.x[which(contams==0)],cat.y[which(contams==0)]),data.frame(cat.x[which(contams==1)],cat.y[which(contams==1)]),k=num.nearest.neighbours)
contam.inside<-NULL
for(ind in 1:length(which(contams==1))) {
#If any of the nearest neighbors overlap with the contaminant /*fold*/ {{{
contam.inside<-c(contam.inside,any(nearest$nn.dist[ind,] < (sqrt(2)/2*(stamplen[which(contams==0)][nearest$nn.idx[ind,]]+stamplen[which(contams==1)][ind]))))
}
inside.mask<-rep(TRUE,cat.len)
inside.mask[which(contams==1)]<-contam.inside
# /*fend*/ }}}
#Remove object catalogue entries /*fold*/ {{{
cat.x<-cat.x[which(inside.mask)]
cat.y<-cat.y[which(inside.mask)]
cat.id<-cat.id[which(inside.mask)]
cat.ra<-cat.ra[which(inside.mask)]
cat.dec<-cat.dec[which(inside.mask)]
cat.theta<-cat.theta[which(inside.mask)]
cat.a<-cat.a[which(inside.mask)]
cat.b<-cat.b[which(inside.mask)]
if (use.segmentation) { seg.x<-seg.x[which(inside.mask)] }
if (use.segmentation) { seg.y<-seg.y[which(inside.mask)] }
if (length(flux.weight)!=1) { flux.weight<-flux.weight[which(inside.mask)] }
contams<-contams[which(inside.mask)]
if (exists("groups")) { groups<-groups[which(inside.mask)] }
if (exists("psf.id")) { psf.id<-psf.id[which(inside.mask)] }
if (num.cores < 0) {
registerDoParallel(cores=(min(ceiling(length(cat.x)/50000),abs(num.cores))))
}
chunk.size=length(cat.id)/getDoParWorkers()
mpi.opts<-list(chunkSize=chunk.size)
message("Number of objects per thread:",chunk.size)
# /*fend*/ }}}
#Notify how many objects remain /*fold*/ {{{
if (verbose) { message(paste("There are",length(cat.x),"supplied objects & contaminants that intersect (",
round(((cat.len-length(cat.x))/cat.len)*100, digits=2),"% of objects were contaminants that didn't intersect galaxies)")) }
# /*fend*/ }}}
if (showtime) { cat(" - Done (",round(proc.time()[3]-timer[3],digits=2),"sec )\n")
message(paste('Remove irrelevant contaminants - Done (',round(proc.time()[3]-timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
}
# /*fend*/ }}}
#Convert decimal pixel values into actual pixel values /*fold*/ {{{
x.pix<-floor(cat.x)
y.pix<-floor(cat.y)
# /*fend*/ }}}
#Create an array of stamps containing the data image subsection for all objects /*fold*/ {{{
timer=system.time(data.stamp<-make.data.array.maps(outenv=environment()))
if (((loop.total<=1)&(cutup)&(length(data.stamp))==0)) { sink(type="message") ; stop("No Aperture Stamps are aligned enough to be valid.") }
else if ((cutup)&(length(data.stamp)==0)) {
if (showtime) { cat(" - BREAK (",round(timer[3],digits=2),"sec )\n")
message(paste('Make Data Mask - BREAK (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - BREAK\n") }
warning("No Aperture Stamps are aligned enough to be valid.")
#Notify & Close Logfile /*fold*/ {{{
if (!is.null(env)) {
on.exit(detach(env), add=TRUE)
}
message(paste('\n-----------------------------------------------------\nDatamap Skipped - No Aperture Stamps are aligned enough to be valid\n'))
if (!quiet) {
cat(paste('\n-----------------------------------------------------\nDatamap Skipped - No Apertures Stamps are aligned enough to be valid'))
}
return(NULL)
# /*fend*/ }}}
}
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(paste('Make Data Mask - Done (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
# /*fend*/ }}}
#Remove Image arrays as they are no longer needed /*fold*/ {{{
#if (cutup) {
# imenvlist<-ls(envir=image.env)
# imenvlist<-imenvlist[which(imenvlist!="im"&imenvlist!="data.hdr")]
# rm(list=imenvlist, envir=image.env)
#}
# /*fend*/ }}}
#Create an array of stamps containing the apertures for all objects /*fold*/ {{{
if (aperture.type <= 1) {
timer=system.time(sa<-make.catalogue.apertures(outenv=environment()))
} else if (aperture.type == 2) {
timer=system.time(sa<-make.gaussian.apertures(outenv=environment()))
} else if (aperture.type == 3) {
timer=system.time(sa<-make.segmentation.apertures(outenv=environment()))
}
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(paste('Make SA Mask - Done (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
# /*fend*/ }}}
#Discard any apertures that were zero'd in the make process /*fold*/ {{{
totsa<-foreach(sam=sa, i=1:length(sa), .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% { if (max(sam)>1){warning(paste("Max of Aperture",i,"is",max(sam))) } ; sum(sam) }
inside.mask<-(totsa > 0)
if ((loop.total<=1)&(length(which(inside.mask==TRUE))==0)) { sink(type="message") ; stop("No Apertures are inside the Mask after Aperture Creation.") } # Nothing inside the mask
else if (length(which(inside.mask==TRUE))==0) {
warning("No Apertures are inside the Mask after Aperture Creation.")
#Notify & Close Logfile /*fold*/ {{{
if (!is.null(env)) {
on.exit(detach(env), add=TRUE)
}
message(paste('\n-----------------------------------------------------\nDatamap Skipped - No Apertures in the Mask\n'))
if (!quiet) {
cat(paste('\n-----------------------------------------------------\nDatamap Skipped - No Apertures in the Mask\n'))
}
return(NULL)
# /*fend*/ }}}
}
# /*fend*/ }}}
#Remove object that failed aperture creation /*fold*/ {{{
sa<-sa[which(inside.mask)]
if (length(data.stamp )>1) { data.stamp<-data.stamp[which(inside.mask)] }
if (length(mask.stamp)>1) { mask.stamp<-mask.stamp[which(inside.mask)] }
if (length(error.stamp)>1) { error.stamp<-error.stamp[which(inside.mask)] }
stamplen<-stamplen[which(inside.mask)]
ap.lims.data.stamp<-rbind(ap.lims.data.stamp[which(inside.mask),])
ap.lims.mask.stamp<-rbind(ap.lims.mask.stamp[which(inside.mask),])
ap.lims.error.stamp<-rbind(ap.lims.error.stamp[which(inside.mask),])
data.stamp.lims<-rbind(data.stamp.lims[which(inside.mask),])
mask.stamp.lims<-rbind(mask.stamp.lims[which(inside.mask),])
error.stamp.lims<-rbind(error.stamp.lims[which(inside.mask),])
ap.lims.data.map<-rbind(ap.lims.data.map[which(inside.mask),])
ap.lims.mask.map<-rbind(ap.lims.mask.map[which(inside.mask),])
ap.lims.error.map<-rbind(ap.lims.error.map[which(inside.mask),])
x.pix<-x.pix[which(inside.mask)]
y.pix<-y.pix[which(inside.mask)]
cat.x<-cat.x[which(inside.mask)]
cat.y<-cat.y[which(inside.mask)]
cat.id<-cat.id[which(inside.mask)]
cat.ra<-cat.ra[which(inside.mask)]
cat.dec<-cat.dec[which(inside.mask)]
cat.theta<-cat.theta[which(inside.mask)]
theta.offset<-theta.offset[which(inside.mask)]
cat.a<-cat.a[which(inside.mask)]
cat.a.pix<-cat.a.pix[which(inside.mask)]
cat.b<-cat.b[which(inside.mask)]
if (use.segmentation) { seg.x<-seg.x[which(inside.mask)] }
if (use.segmentation) { seg.y<-seg.y[which(inside.mask)] }
if (length(flux.weight)!=1) { flux.weight<-flux.weight[which(inside.mask)] }
if (exists("contams")) { contams<-contams[which(inside.mask)] }
if (exists("groups")) { groups<-groups[which(inside.mask)] }
if (exists("psf.id")) { psf.id<-psf.id[which(inside.mask)] }
inside.mask<-inside.mask[which(inside.mask)]
if (num.cores < 0) {
registerDoParallel(cores=(min(ceiling(length(cat.x)/50000),abs(num.cores))))
}
chunk.size=ceiling(length(cat.id)/getDoParWorkers())
mpi.opts<-list(chunkSize=chunk.size)
message("Number of objects per thread:",chunk.size)
# /*fend*/ }}}
#Re-Initialise object count /*fold*/ {{{
npos<-length(cat.id)
# /*fend*/ }}}
#Create a full mask of all apertures in their correct image-space locations /*fold*/ {{{
timer=system.time(image.env$aa<-make.aperture.map(outenv=environment(), sa, dimim))
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(paste('Make A Mask - Done (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
# /*fend*/ }}}
#If wanted, output the All Apertures Mask /*fold*/ {{{
if (make.all.apertures.map) {
if (!quiet) { cat(paste('Outputting All Apertures Mask to',all.apertures.map.filename," ")) }
#Write All Apertures Mask to file /*fold*/ {{{
timer=system.time(write.fits.image.file(file.path(path.root,path.work,path.out,all.apertures.map.filename),image.env$aa,image.env$data.hdr,nochange=TRUE))
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(paste('Output FA Mask - Done (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
}# /*fend*/ }}}
# /*fend*/ }}}
# /*fend*/ }}}
#PART TWO: SINGLE FILTERED APERTURE MASKS /*fold*/ {{{
#Details /*fold*/ {{{
#If wanted, perform the convolution of the
#apertures with the PSF. If not wanted,
#duplicate the single apertures.
#Also, perform flux weighting. /*fend*/ }}}
#If wanted, use pixel fluxes as flux.weights /*fold*/ {{{
if (use.pixel.fluxweight) {
#Image Flux at central pixel /*fold*/ {{{
cat("Determine Image Flux at central pixel ")
#if (cutup) {
# pixflux<-foreach(xp=x.pix-(data.stamp.lims[,1]-1),yp=y.pix-(data.stamp.lims[,3]-1), im=data.stamp, .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% {
# im[xp,yp]
# }
#} else {
# pixflux<-foreach(xp=x.pix,yp=y.pix, .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment()), .export=c('image.env')) %dopar% {
# image.env$im[xp,yp]
# }
#}
pixflux<-image.env$im[cbind(x.pix,y.pix)]
if (verbose) { cat(" - Done\n") }
# /*fend*/ }}}
#Determine Noise Characteristics /*fold*/ {{{
if (verbose) { cat("Determine Image Noise Characteristics ") }
if (any(pixflux==-99)) { warning("Some pixel flux determinations failed"); pixflux[which(pixflux==-99)]<-NA }
quan<-quantile(image.env$im,c(0.1,0.9),na.rm=TRUE)
bw<-abs(quan[2]-quan[1])/100/sqrt(12)
if (bw<=0) {
bw<-diff(range(image.env$im))/100/sqrt(12)
}
mode<-density(as.numeric(image.env$im),from=quan[1],to=quan[2],kernel='rect',bw=bw,na.rm=TRUE)
mode<-mode$x[which.max(mode$y)]
mad<-median(abs(image.env$im-mode),na.rm=TRUE)
if (length(which(pixflux>mode+mad))==0) {
message("WARNING: No objects have pixel flux measurements that are > Pixel Mode+MAD. Pixel Flux weighting cannot be used\n")
flux.weight<-1
} else {
if (weight.type=='flux') {
message("Using linear scaling of pixel flux measurements as initial weight\n")
flux.weight<-magmap(pixflux, lo=mode+mad, hi=max(pixflux), range=c(0.01,1), type="num", stretch='lin',bad=0.01)$map
} else if (weight.type=='mag') {
message("Using logarithmic scaling of pixel flux measurements as initial weight\n")
flux.weight<-magmap(pixflux, lo=mode+mad, hi=max(pixflux), range=c(0.01,1), type="num", stretch='log',bad=0.01)$map
} else {
message("Using histogram scaling of pixel flux measurements as initial weight\n")
flux.weight<-magmap(pixflux, lo=mode+mad, hi=max(pixflux), range=c(0.01,1), type="num", stretch='cdf',bad=0.01)$map
}
}
cat(" - Done\n")
# /*fend*/ }}}
}# /*fend*/ }}}
#If needed, do Convolutions & Fluxweightings /*fold*/ {{{
if ((!psf.filt)&(length(which(flux.weight!=1))!=0)) {
#No Convolution, Need Fluxweighting /*fold*/ {{{
#Details /*fold*/ {{{
#If not convolving with psf, and if all the flux.weights are not unity,
#then skip filtering, duplicate the arrays, and only weight the stamps/apertures
# /*fend*/ }}}
#No Convolution, duplicate apertures /*fold*/ {{{
if (verbose) { message("No Convolution: Convolved Apertures are identical to Simple Apertures") }
sfa<-sa
image.env$fa<-image.env$aa
# /*fend*/ }}}
#Perform aperture weighting /*fold*/ {{{
if (verbose) { message("Fluxweights Present: Weighting Convolved Apertures") }
timer=system.time(wsfa<-make.convolved.apertures(outenv=environment(), sa,flux.weightin=flux.weight,immask=mask.stamp))
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(paste('Make WSFA Mask - Done (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
# /*fend*/ }}}
#Create a full mask of stamps/apertures /*fold*/ {{{
timer=system.time(image.env$wfa<-make.aperture.map(outenv=environment(), wsfa, dimim))
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(paste('Make WFA Mask - Done (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
# /*fend*/ }}}
# /*fend*/ }}}
} else if ((!psf.filt)&(length(which(flux.weight!=1))==0)) {
#No Convolution, No Fluxweighting /*fold*/ {{{
#Details /*fold*/ {{{
#If not convolving with psf, and if all the flux.weights are unity,
#then skip filtering and weighting, and simply duplicate the arrays /*fend*/ }}}
#Duplicate arrays /*fold*/ {{{
if (verbose) { message("No Convolution: Convolved Apertures are identical to Simple Apertures")
message("No Fluxweights: Weighted Convolved Apertures are identical to Convolved Apertures") }
sfa<-sa
image.env$fa<-image.env$aa
wsfa<-sfa
image.env$wfa<-image.env$fa
# /*fend*/ }}}
# /*fend*/ }}}
} else if ((psf.filt)&(length(which(flux.weight!=1))!=0)) {
#Convolving PSF, Need Fluxweighting /*fold*/ {{{
#Details /*fold*/ {{{
#If convolving with psf, and if all the flux.weights are not unity,
#then colvolve and weight the stamps/apertures /*fend*/ }}}
#Perform Convolution /*fold*/ {{{
if (verbose) { message("PSF Present: Making Convolved Apertures") }
timer=system.time(sfa<-make.convolved.apertures(outenv=environment(), sa,immask=mask.stamp))
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(paste('Make SFA Mask - Done (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
# /*fend*/ }}}
#Create a full mask of all convolved stamps/apertures /*fold*/ {{{
timer=system.time(image.env$fa<-make.aperture.map(outenv=environment(), sfa, dimim))
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(paste('Make FA Mask - Done (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
# /*fend*/ }}}
#Perform aperture weighting /*fold*/ {{{
if (verbose) { message("Fluxweights Present: Weighting Convolved Apertures") }
timer=system.time(wsfa<-make.convolved.apertures(outenv=environment(), sa,flux.weightin=flux.weight,immask=mask.stamp))
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(paste('Make WSFA Mask - Done (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
# /*fend*/ }}}
#Create a full mask of all the convolved weighted stamps/apertures /*fold*/ {{{
timer=system.time(image.env$wfa<-make.aperture.map(outenv=environment(), wsfa, dimim))
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(paste('Make WFA Mask - Done (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
# /*fend*/ }}}
# /*fend*/ }}}
} else if ((psf.filt)&(length(which(flux.weight!=1))==0)) {
#Convolve PSF, No Fluxweighting /*fold*/ {{{
#Details /*fold*/ {{{
#If convolving with psf, and if all the flux.weights are unity,
#colvolve, and then duplicate the stamps/apertures /*fend*/ }}}
#Perform Convolution /*fold*/ {{{
if (verbose) { message("PSF Present: Making Convolved Apertures") }
timer=system.time(sfa<-make.convolved.apertures(outenv=environment(), sa,immask=mask.stamp))
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(paste('Make SFA Mask - Done (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
# /*fend*/ }}}
#Create a full mask of all convolved stamps/apertures /*fold*/ {{{
timer=system.time(image.env$fa<-make.aperture.map(outenv=environment(), sfa, dimim))
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(paste('Make FA Mask - Done (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
# /*fend*/ }}}
#No Fluxweights, duplicate apertures /*fold*/ {{{
if (verbose) { message("No Fluxweights: Weighted Convolved Apertures are identical to Convolved Apertures") }
wsfa<-sfa
image.env$wfa<-image.env$fa
# /*fend*/ }}}
# /*fend*/ }}}
}# /*fend*/ }}}
#Discard any apertures that were zero'd in the make process /*fold*/ {{{
totsfa<-foreach(sam=sfa, i=1:length(sfa), .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% { if (max(sam)>1){warning(paste("Max of Aperture",i,"is",max(sam))) } ; sum(sam) }
inside.mask<-(totsfa > 0)
if ((loop.total<=1)&(length(which(inside.mask==TRUE))==0)) { sink(type="message") ; stop("No Apertures Remain within the mask after convolution.") } # Nothing inside the mask
else if (length(which(inside.mask==TRUE))==0) {
warning("No Apertures Remain within the mask after convolution.")
#Notify & Close Logfile /*fold*/ {{{
if (!is.null(env)) {
on.exit(detach(env), add=TRUE)
}
message(paste('\n-----------------------------------------------------\nDatamap Skipped - No Apertures in the Mask\n'))
if (!quiet) {
cat(paste('\n-----------------------------------------------------\nDatamap Skipped - No Apertures in the Mask\n'))
}
return(NULL)
# /*fend*/ }}}
}
#Remove object that failed aperture creation /*fold*/ {{{
sa<-sa[which(inside.mask)]
sfa<-sfa[which(inside.mask)]
wsfa<-wsfa[which(inside.mask)]
if (length(data.stamp )>1) { data.stamp<-data.stamp[which(inside.mask)] }
if (length(mask.stamp)>1) { mask.stamp<-mask.stamp[which(inside.mask)] }
if (length(error.stamp)>1) { error.stamp<-error.stamp[which(inside.mask)] }
stamplen<-stamplen[which(inside.mask)]
ap.lims.data.stamp<-rbind(ap.lims.data.stamp[which(inside.mask),])
ap.lims.mask.stamp<-rbind(ap.lims.mask.stamp[which(inside.mask),])
ap.lims.error.stamp<-rbind(ap.lims.error.stamp[which(inside.mask),])
data.stamp.lims<-rbind(data.stamp.lims[which(inside.mask),])
mask.stamp.lims<-rbind(mask.stamp.lims[which(inside.mask),])
error.stamp.lims<-rbind(error.stamp.lims[which(inside.mask),])
ap.lims.data.map<-rbind(ap.lims.data.map[which(inside.mask),])
ap.lims.mask.map<-rbind(ap.lims.mask.map[which(inside.mask),])
ap.lims.error.map<-rbind(ap.lims.error.map[which(inside.mask),])
x.pix<-x.pix[which(inside.mask)]
y.pix<-y.pix[which(inside.mask)]
cat.x<-cat.x[which(inside.mask)]
cat.y<-cat.y[which(inside.mask)]
cat.id<-cat.id[which(inside.mask)]
cat.ra<-cat.ra[which(inside.mask)]
cat.dec<-cat.dec[which(inside.mask)]
cat.theta<-cat.theta[which(inside.mask)]
theta.offset<-theta.offset[which(inside.mask)]
cat.a<-cat.a[which(inside.mask)]
cat.a.pix<-cat.a.pix[which(inside.mask)]
cat.b<-cat.b[which(inside.mask)]
if (use.segmentation) { seg.x<-seg.x[which(inside.mask)] }
if (use.segmentation) { seg.y<-seg.y[which(inside.mask)] }
if (length(flux.weight)!=1) { flux.weight<-flux.weight[which(inside.mask)] }
if (exists("contams")) { contams<-contams[which(inside.mask)] }
if (exists("groups")) { groups<-groups[which(inside.mask)] }
if (exists("psf.id")) { psf.id<-psf.id[which(inside.mask)] }
inside.mask<-inside.mask[which(inside.mask)]
if (num.cores < 0) {
registerDoParallel(cores=(min(ceiling(length(cat.x)/50000),abs(num.cores))))
}
chunk.size=ceiling(length(cat.id)/getDoParWorkers())
mpi.opts<-list(chunkSize=chunk.size)
message("Number of objects per thread:",chunk.size)
# /*fend*/ }}}
# /*fend*/ }}}
#Re-Initialise object count /*fold*/ {{{
npos<-length(cat.id)
# /*fend*/ }}}
#Remove Uneeded Image /*fold*/ {{{
rm(aa, envir=image.env)
gc()
# /*fend*/ }}}
#-----Diagnostic-----# /*fold*/ {{{
if (diagnostic) {
if (verbose) { message("Checking Apertures for scaling errors") }
foreach(sam=sa, sfam=sfa, i=1:length(sa), .combine=function(a,b){NULL},.inorder=FALSE, .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% {
if (max(sam)>1){warning(paste("Max of Aperture",i,"is",max(sam))) }
if (max(sfam)>1){warning(paste("Max of Convolved Aperture",i,"is",max(sfam))) }
}
}# /*fend*/ }}}
#Do we want to plot a sample of the apertures? /*fold*/ {{{
if (plot.sample) {
#Set output name /*fold*/ {{{
PlotDev(file=file.path(path.root,path.work,path.out,paste0("Aperture_Samples.",plot.device)))
# /*fend*/ }}}
#Set Layout /*fold*/ {{{
par(mfrow=c(2,2))
# /*fend*/ }}}
#Output a random 15 apertures to file /*fold*/ {{{
ind1=which(cat.a>0)
ind1=ind1[order(runif(1:length(ind1)))][1:15]
ind2=order(runif(1:npos))[1:15]
ind<-c(ind1, ind2)
rm(ind1)
rm(ind2)
ind<-ind[which(!is.na(ind))]
for (i in ind) {
Rast<-ifelse(stamplen[i]>100,TRUE,FALSE)
image(sa[[i]], main="Single Aperture (SA)", asp=1, col=heat.colors(256), useRaster=Rast)
points(ceiling(stamplen[i]/2)/stamplen[i],ceiling(stamplen[i]/2)/stamplen[i],pch="+",lw=2.0, col="red")
image(sfa[[i]], main="Single Convolved Apertrure (SFA)", asp=1, col=heat.colors(256), useRaster=Rast)
points(ceiling(stamplen[i]/2)/stamplen[i],ceiling(stamplen[i]/2)/stamplen[i],pch="+",lw=2.0, col="red")
}# /*fend*/ }}}
#Close the file /*fold*/ {{{
if (!grepl('x11',plot.device,ignore.case=TRUE)) { dev.off() }
# /*fend*/ }}}
}
# /*fend*/ }}}
#If wanted, output the Convolved & Weighted Aperture Mask /*fold*/ {{{
if (make.convolved.apertures.map) {
if (!quiet) { cat(paste('Outputting All Convolved Apertures Mask to',fa.filename," ")) }
timer=system.time(write.fits.image.file(file.path(path.root,path.work,path.out,fa.filename),image.env$wfa,image.env$data.hdr,nochange=TRUE) )
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(paste('Output FA Mask - Done (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
}# /*fend*/ }}}
#If wanted, create (& possibly output) the SourceMask /*fold*/ {{{
if (sourcemask) {
#Details /*fold*/ {{{
#> In the case of flux measurements: we want SM to be 1 only where the sources are, and within the mask region
#> In the case of creating a source mask: we want SM to be 1 within the mask region in between sources only
#To do this we:
# a) set the sm array = image mask (0s outside image, 1s inside)
# b) set all nonzero locations in fa to 0 in sm array
# /*fend*/ }}}
#Make Sourcemask /*fold*/ {{{
if (!exists("transmission.map")) { transmission.map<-FALSE }
if (!quiet) { cat(paste("Creating Sourcemask ")) }
if (length(image.env$imm)>1) {
sm<-image.env$imm.orig<-image.env$imm
} else {
sm<-array(1, dim=dimim)
}
#Mask the region where apertures cannot be placed (around edges)
minlen<-ceiling(min(stamplen)/2)
sm[1:minlen,]<-0
sm[,1:minlen]<-0
sm[dim(sm)[1]-(1:minlen-1),]<-0
sm[,dim(sm)[2]-(1:minlen-1)]<-0
#sm<-array(1, dim=dimim)
#Get Mask as is {{{
if (length(mask.stamp)>1) {
for (i in 1:npos) {
sm[mask.stamp.lims[i,1]:mask.stamp.lims[i,2],mask.stamp.lims[i,3]:mask.stamp.lims[i,4]]<-mask.stamp[[i]]
#sm[data.stamp.lims[i,1]:data.stamp.lims[i,2],data.stamp.lims[i,3]:data.stamp.lims[i,4]]<-mask.stamp[[i]]
}
} #}}}
mask.xrange<-c(min(mask.stamp.lims[,1]),max(mask.stamp.lims[,2]))
mask.yrange<-c(min(mask.stamp.lims[,3]),max(mask.stamp.lims[,4]))
data.xrange<-c(min(data.stamp.lims[,1]),max(data.stamp.lims[,2]))
data.yrange<-c(min(data.stamp.lims[,3]),max(data.stamp.lims[,4]))
if (!psf.filt) {
if (transmission.map) {
sm[mask.xrange[1]:mask.xrange[2],mask.yrange[1]:mask.yrange[2]][which(!image.env$fa[data.xrange[1]:data.xrange[2],data.yrange[1]:data.yrange[2]] > 0)]<-0
#sm[which(!image.env$fa > 0)]<-0
} else {
sm[mask.xrange[1]:mask.xrange[2],mask.yrange[1]:mask.yrange[2]][which(image.env$fa[data.xrange[1]:data.xrange[2],data.yrange[1]:data.yrange[2]] > 0)]<-0
#sm[which(image.env$fa > 0)]<-0
}
} else {
if (transmission.map) {
sm[mask.xrange[1]:mask.xrange[2],mask.yrange[1]:mask.yrange[2]][which(!image.env$fa[data.xrange[1]:data.xrange[2],data.yrange[1]:data.yrange[2]] > max(sapply(psf,max,na.rm=TRUE), na.rm=TRUE)*(1-sourcemask.conf.lim))]<-0
#sm[which(!image.env$fa > max(psf, na.rm=TRUE)*(1-sourcemask.conf.lim))]<-0
} else {
sm[mask.xrange[1]:mask.xrange[2],mask.yrange[1]:mask.yrange[2]][which(image.env$fa[data.xrange[1]:data.xrange[2],data.yrange[1]:data.yrange[2]] > max(sapply(psf,max,na.rm=TRUE), na.rm=TRUE)*(1-sourcemask.conf.lim))]<-0
#sm[which(image.env$fa > max(psf, na.rm=TRUE)*(1-sourcemask.conf.lim))]<-0
}
}
#Cutup again
if (do.sky.est||get.sky.rms||blank.cor) {
if (cutup) {
old.mask.stamp<-mask.stamp
mask.stamp<-list(NULL)
for (i in 1:npos) {
mask.stamp[[i]]<-sm[mask.stamp.lims[i,1]:mask.stamp.lims[i,2],mask.stamp.lims[i,3]:mask.stamp.lims[i,4]]
}
} else {
image.env$imm.dimim<-array(1, dim=dimim)
image.env$imm.dimim[data.xrange[1]:data.xrange[2],data.yrange[1]:data.yrange[2]]<-sm[mask.xrange[1]:mask.xrange[2],mask.yrange[1]:mask.yrange[2]]
image.env$imm<-sm
mask.stamp<-sm
}
}
if (!quiet) { cat(" - Done\n") }
# /*fend*/ }}}
#-----Diagnostic-----# /*fold*/ {{{
if (diagnostic) {
message(paste("SourceMask Max/Min:",max(sm),min(sm)))
message(paste("OLDMethod - SourceMask Max/Min:",max(1-image.env$fa),min(1-image.env$fa)))
}# /*fend*/ }}}
#If wanted, output the SourceMask /*fold*/ {{{
if (!exists("sourcemask.out")) { sourcemask.out<-FALSE }
if (sourcemask.out){
if (!quiet) { cat(paste('Outputting Source Mask to',sourcemask.filename," ")) }
write.fits.image.file(file.path(path.root,path.work,path.out,sourcemask.filename),sm,image.env$data.hdr,nochange=TRUE)
if (!quiet) { cat(" - Done\n") }
}# /*fend*/ }}}
#If we want the SourceMask only, then end here /*fold*/ {{{
if (sourcemask.only) {
if (!quiet) { cat("SourceMaskOnly Flag Set\n") }
return()
} # /*fend*/ }}}
} else if (plot.sample) {
#Set the mask to 1 everywhere /*fold*/ {{{
sm<-1
mask.stamp<-1
# /*fend*/ }}}
}# /*fend*/ }}}
#Remove arrays that are no longer needed /*fold*/ {{{
rm(fa, envir=image.env)
gc()
# /*fend*/ }}}
# /*fend*/ }}}
#PART THREE: DEBLENDING /*fold*/ {{{
#If wanted, perform sky estimation. Otherwise set to NA /*fold*/ {{{
skylocal<-skyflux<-skyerr<-skyrms<-skypval<-detecthres<-detecthres.mag<-
skyNBinNear<-skyNBinNear.mean<-skyflux.mean<-skyerr.mean<-skylocal.mean<-
skypval.mean<-skyrms.mean<-rep(NA,length(sfa))
if (do.sky.est||get.sky.rms) {
if (!exists('quick.sky')) { quick.sky<-FALSE }
#Get sky estimates /*fold*/ {{{
#SkyEst is NP hard, so set the processors to Max
if (num.cores < 0) { registerDoParallel(cores=abs(num.cores)) }
#Perform Sky Estimation /*fold*/ {{{
if (quick.sky) {
if (!quiet) { message("Perfoming Fast Sky Estimation"); cat("Performing Fast Sky Estimation") }
if (exists('skyest')) {
if (fit.sky) {
skyest<-skyest[which(is.finite(skyest[,'skyMu'])),]
subs<-which((cat.id%in%skyest[,'sources']))
subs2<-which((skyest[,'sources']%in%cat.id))
skylocal[subs]<-skyest[subs2,'skyMu']
skyerr[subs]<-skyest[subs2,'skyMuErr']*correl.noise
skyrms[subs]<-skyest[subs2,'skySD']
} else {
skyest<-skyest[which(is.finite(skyest[,'skyMedian'])),]
subs<-which((cat.id%in%skyest[,'sources']))
subs2<-which((skyest[,'sources']%in%cat.id))
skylocal[subs]<-skyest[subs2,'skyMedian']
skyerr[subs]<-skyest[subs2,'skyMedianErr']*correl.noise
skyrms[subs]<-skyest[subs2,'skyRMS']
}
skypval[subs]<-skyest[subs2,'skyRMSpval']
skyNBinNear[subs]<-skyest[subs2,'Nnearsky']
skylocal.mean[subs]<-skyest[subs2,'skyMedian']
skyerr.mean[subs]<-skyest[subs2,'skyMedianErr']*correl.noise
skyrms.mean[subs]<-skyest[subs2,'skyRMS']
skypval.mean[subs]<-skyest[subs2,'skyRMSpval']
skyNBinNear.mean[subs]<-skyest[subs2,'Nnearsky']
subs<-which(!(cat.id%in%skyest[,'sources']))
rm('skyest')
} else {
subs<-1:length(cat.x)
}
if (cutup) {
timer<-system.time(skyest<-fast.sky.estimate(cat.x=cat.x,cat.y=cat.y,data.stamp.lims=data.stamp.lims,fit.gauss=fit.sky,
cutlo=(cat.a/arcsec.per.pix),cuthi=(cat.a/arcsec.per.pix)*5,data.stamp=data.stamp,mask.stamp=mask.stamp,
clipiters=sky.clip.iters,sigma.cut=sky.clip.prob,PSFFWHMinPIX=psffwhm, mpi.opts=mpi.opts,subset=subs,saturate=image.env$saturation))
} else {
timer<-system.time(skyest<-fast.sky.estimate(cat.x=cat.x,cat.y=cat.y,data.stamp.lims=data.stamp.lims,fit.gauss=fit.sky,
cutlo=(cat.a/arcsec.per.pix),cuthi=(cat.a/arcsec.per.pix)*5,
data.stamp=image.env$im, mask.stamp=image.env$imm.dimim,saturate=image.env$saturation,
clipiters=sky.clip.iters,sigma.cut=sky.clip.prob,PSFFWHMinPIX=psffwhm, mpi.opts=mpi.opts,subset=subs))
}
#Get sky parameters /*fold*/ {{{
if (fit.sky) {
skylocal[subs]<-skyest[,'skyMu']
skyerr[subs]<-skyest[,'skyMuErr']*correl.noise
skyrms[subs]<-skyest[,'skySD']
} else {
skylocal[subs]<-skyest[,'skyMedian']
skyerr[subs]<-skyest[,'skyMedianErr']*correl.noise
skyrms[subs]<-skyest[,'skyRMS']
}
skypval[subs]<-skyest[,'skyRMSpval']
skyNBinNear[subs]<-skyest[,'Nnearsky']
skylocal.mean[subs]<-skyest[,'skyMedian']
skyerr.mean[subs]<-skyest[,'skyMedianErr']*correl.noise
skyrms.mean[subs]<-skyest[,'skyRMS']
skypval.mean[subs]<-skyest[,'skyRMSpval']
skyNBinNear.mean[subs]<-skyest[,'Nnearsky']
#}}}
} else {
if (!quiet) { message("Perfoming Sky Estimation"); cat("Performing Sky Estimation") }
#Check if we can cutdown on the workload {{{
if (exists('skyest')) {
skyest<-skyest[which(is.finite(skyest[,'sky'])),]
subs<-which(cat.id%in%skyest[,'sources'])
subs2<-which(skyest[,'sources']%in%cat.id)
skylocal[subs]<-skyest[subs2,'sky']
skylocal.mean[subs]<-skyest[subs2,'sky.mean']
skyerr[subs]<-skyest[subs2,'skyerr']*correl.noise
skyrms[subs]<-skyest[subs2,'skyRMS']
skypval[subs]<-skyest[subs2,'skyRMSpval']
skyNBinNear[subs]<-skyest[subs2,'Nnearsky']
skyerr.mean[subs]<-skyest[subs2,'skyerr.mean']*correl.noise
skyrms.mean[subs]<-skyest[subs2,'skyRMS.mean']
skypval.mean[subs]<-skyest[subs2,'skyRMSpval.mean']
skyNBinNear.mean[subs]<-skyest[subs2,'Nnearsky.mean']
subs<-which(!(cat.id%in%skyest[,'sources']))
rm('skyest')
} else {
subs<-1:length(cat.x)
}
if (cutup) {
timer<-system.time(skyest<-sky.estimate(cat.x=cat.x,cat.y=cat.y,data.stamp.lims=data.stamp.lims,saturate=image.env$saturation,
cutlo=(cat.a/arcsec.per.pix),cuthi=(cat.a/arcsec.per.pix)*5,data.stamp=data.stamp,mask.stamp=mask.stamp,
clipiters=sky.clip.iters,sigma.cut=sky.clip.prob,PSFFWHMinPIX=psffwhm, mpi.opts=mpi.opts,subset=subs))
} else {
timer<-system.time(skyest<-sky.estimate(cat.x=cat.x,cat.y=cat.y,data.stamp.lims=data.stamp.lims,saturate=image.env$saturation,
cutlo=(cat.a/arcsec.per.pix),cuthi=(cat.a/arcsec.per.pix)*5,
data.stamp=image.env$im, mask.stamp=image.env$imm.dimim,
clipiters=sky.clip.iters,sigma.cut=sky.clip.prob,PSFFWHMinPIX=psffwhm, mpi.opts=mpi.opts,subset=subs))
}
#Get sky parameters /*fold*/ {{{
skylocal[subs]<-skyest[,'sky']
skylocal.mean[subs]<-skyest[,'sky.mean']
skyerr[subs]<-skyest[,'skyerr']*correl.noise
skyrms[subs]<-skyest[,'skyRMS']
skypval[subs]<-skyest[,'skyRMSpval']
skyNBinNear[subs]<-skyest[,'Nnearsky']
skyerr.mean[subs]<-skyest[,'skyerr.mean']*correl.noise
skyrms.mean[subs]<-skyest[,'skyRMS.mean']
skypval.mean[subs]<-skyest[,'skyRMSpval.mean']
skyNBinNear.mean[subs]<-skyest[,'Nnearsky.mean']
#}}}
}
# /*fend*/ }}}
#Get sky parameters /*fold*/ {{{
if(!is.na(as.numeric(sky.default))) {
sky.default<-as.numeric(sky.default)
rmsdefault<-0.0
errdefault<-0.0
message("Using a default value for Sky Estimaes that have failed. Error on this value is unknown, so is assumed 0. Similarly, RMS is assumed 0")
} else if (grepl("median",sky.default)) {
sky.default<-median(skylocal,na.rm=TRUE)
rmsdefault<-median(skyrms,na.rm=TRUE)
errdefault<-median(skyerr,na.rm=TRUE)
} else if (grepl("mean",sky.default)) {
sky.default<-mean(skylocal, na.rm=TRUE)
rmsdefault<-mean(skyrms, na.rm=TRUE)
errdefault<-mean(skyerr, na.rm=TRUE)
}
if (is.na(sky.default)) {
warning("All sky estimates have failed, and so the requested sky default is NA.\nTo stop bad behaviour, sky values are all being set to 0")
sky.default<-0
rmsdefault<-0
errdefault<-0
}
skylocal[which(is.na(skylocal))]<-sky.default
skyrms[which(is.na(skyrms))]<-rmsdefault
skyerr[which(is.na(skyerr))]<-errdefault
# /*fend*/ }}}
#Notify /*fold*/ {{{
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(paste('Sky Estimate - Done (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
# /*fend*/ }}}
#If wanted, plot some of the Sky Estimates /*fold*/ {{{
if (plot.sample & !quick.sky) {
if (!quiet) { message("Plotting Sky Estimation"); cat("Plotting Sky Estimation") }
if (cutup) {
timer<-system.time(plot.sky.estimate(cat.id=cat.id,cat.x=cat.x,cat.y=cat.y,
data.stamp.lims=data.stamp.lims,
cutlo=(cat.a/arcsec.per.pix),cuthi=(cat.a/arcsec.per.pix)*5,
data.stamp=data.stamp,mask.stamp=mask.stamp,plot.device=plot.device,
clipiters=sky.clip.iters,sigma.cut=sky.clip.prob,PSFFWHMinPIX=psffwhm,plot.sci=plot.sci,contams=contams,plot.all=plot.all,
path=file.path(path.root,path.work,path.out),rem.mask=TRUE,toFile=TRUE))
} else {
timer<-system.time(plot.sky.estimate(cat.id=cat.id,cat.x=cat.x,cat.y=cat.y,
data.stamp.lims=data.stamp.lims,
cutlo=(cat.a/arcsec.per.pix),cuthi=(cat.a/arcsec.per.pix)*5,
data.stamp=image.env$im,mask.stamp=image.env$imm.dimim,plot.device=plot.device,
clipiters=sky.clip.iters,sigma.cut=sky.clip.prob,PSFFWHMinPIX=psffwhm,plot.sci=plot.sci,contams=contams,plot.all=plot.all,
path=file.path(path.root,path.work,path.out),rem.mask=TRUE,toFile=TRUE))
}
#Notify /*fold*/ {{{
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(paste('Plot Sky Estimate - Done (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
# /*fend*/ }}}
}
# /*fend*/ }}}
# /*fend*/ }}}
#Calculate Detection Thresholds /*fold*/ {{{
if (!quiet) { message("Calculating Detection Limits"); cat("Calculating Detection Limits") }
detecthres<-foreach(sfam=sfa, srms=skyrms, .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% { 5*srms*sqrt(sum(sfam)) }
if (magnitudes) {
suppressWarnings(detecthres.mag<--2.5*(log10(detecthres)-log10(ab.vega.flux))+mag.zp)
} else {
detecthres.mag<-array(NA, dim=c(length(sfa)))
}
if (!quiet) { message(paste(" - Done\n")); cat(" - Done\n")}
# /*fend*/ }}}
# /*fend*/ }}}
}# /*fend*/ }}}
#Perform Deblending of apertures /*fold*/ {{{
timer=system.time(dbw<-make.deblend.weight.maps(outenv=environment()))
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(paste('Make Deblended Weightmap - Done (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
# /*fend*/ }}}
#Remove arrays that are no longer needed /*fold*/ {{{
rm(wfa, envir=image.env)
# /*fend*/ }}}
#Finalise SFA Apertures /*fold*/ {{{
sfabak<-sfa
if (!exists("psf.weighted")) { psf.weighted<-FALSE }
if (psf.filt & !psf.weighted) {
# If we have convolved with a PSF, & not doing PSF Weighting, convert back to tophat apertures (now expanded by PSF convolution) /*fold*/ {{{
# For each aperture, Binary filter at aplim*max(ap) or optimally /*fold*/ {{{
sba<-foreach(sfam=sfa, .options.mpi=mpi.opts, .noexport=ls(envir=environment()), .export="ap.limit")%dopar%{
apvals<-rev(sort(sfam))
if (optimal.aper) {
#If optimal, flip at point that SNR peaks
tempsum<-cumsum(apvals)/sqrt(cumsum(seq(apvals)))
apLim<-apvals[which.max(tempsum)]
sfam[which(sfam < apLim)]<-0
sfam[which(sfam >= apLim)]<-1
} else {
#If not optimal, flip at integration point from par file
tempsum<-cumsum(apvals)
tempfunc<-approxfun(tempsum,apvals)
apLim<-tempfunc(ap.limit*max(tempsum, na.rm=TRUE))
sfam[which(sfam < apLim)]<-0
sfam[which(sfam >= apLim)]<-1
}
return=sfam
}
# /*fend*/ }}}
#Update the sfa /*fold*/ {{{
sfa<-sba
rm(sba)
# /*fend*/ }}}
#If wanted, output the Convolved & Weighted Aperture Mask /*fold*/ {{{
if (make.convolved.apertures.map) {
if (!quiet) { cat(paste('Updated Apertures; ')) }
timer=system.time(image.env$wfa<-make.aperture.map(outenv=environment(), sfa, dimim))
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(paste('Make FA Mask - Done (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
fa.filename<-paste("final_",fa.filename,sep="")
if (!quiet) { cat(paste('Outputting Re-Boxcar-ed All Convolved Apertures Mask to',fa.filename," ")) }
timer=system.time(write.fits.image.file(file.path(path.root,path.work,path.out,fa.filename),image.env$wfa,image.env$data.hdr,nochange=TRUE) )
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(paste('Output FA Mask - Done (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
rm(wfa, envir=image.env)
}# /*fend*/ }}}
#remove unneeded variables /*fold*/ {{{
rm(psfvals)
rm(tempsum)
rm(tempfunc)
# /*fend*/ }}}
# /*fend*/ }}}
} else if (psf.filt & psf.weighted) {
#If we have convolved by the PSF & want PSF Weighting /*fold*/ {{{
#If wanted, output the Convolved & Weighted Aperture Mask /*fold*/ {{{
if (make.convolved.apertures.map) {
if (!quiet) { cat(paste('Updated Apertures; ')) }
timer=system.time(image.env$wfa<-make.aperture.map(outenv=environment(), sfa, dimim))
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(paste('Make FA Mask - Done (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
fa.filename<-paste("final_",fa.filename,sep="")
if (!quiet) { cat(paste('Outputting Final All Convolved Apertures Mask to',fa.filename," ")) }
timer=system.time(write.fits.image.file(file.path(path.root,path.work,path.out,fa.filename),image.env$wfa,image.env$data.hdr,nochange=TRUE) )
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(paste('Output FA Mask - Done (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
rm(wfa, envir=image.env)
}# /*fend*/ }}}
#remove unneeded variables /*fold*/ {{{
rm(psfvals)
rm(tempsum)
rm(tempfunc)
# /*fend*/ }}}
# /*fend*/ }}}
}
# /*fend*/ }}}
#If PSF Supplied & sample wanted, Plot PSF and Min Aperture Correction /*fold*/ {{{
if (!(no.psf)&(plot.sample)) {
#PSF with Contours /*fold*/ {{{
for (i in 1:length(psf)) {
if (any(psf.id==i)) {
PlotDev(file=file.path(path.root,path.work,path.out,paste0("PSF_",i,".",plot.device)))
psfvals<-rev(sort(psf[[i]][which(is.finite(psf[[i]]),arr.ind=TRUE)]))
tempsum<-cumsum(psfvals)
tempfunc<-approxfun(tempsum,psfvals)
psfLimit<-tempfunc(ap.limit*max(tempsum, na.rm=TRUE))
suppressWarnings(image(log10(psf[[i]]),main="PSF & Binary Contour Levels", asp=1,col=heat.colors(256),useRaster=ifelse(length(psf[[i]])>1E4,TRUE,FALSE)))
contour(psf[[i]], levels=(tempfunc(c(0.5,0.9,0.95,0.99,0.999,0.9999)*max(tempsum,na.rm=TRUE))), labels=c(0.5,0.9,0.95,0.99,0.999,0.9999), col='blue', add=TRUE)
if (psf.filt) { contour(psf[[i]], levels=c(psfLimit), labels=c(ap.limit), col='green', add=TRUE) }
if (!grepl('x11',plot.device,ignore.case=TRUE)) { dev.off() }
}
}
# /*fend*/ }}}
#Plot Example of PS minimum aperture Correction; minimum source /*fold*/ {{{
ind<-(which(cat.a==0)[1])
# /*fend*/ }}}
#If no point sources, use the minimum aperture /*fold*/ {{{
if (is.na(ind)) { ind<-which.min(cat.a) }
# /*fend*/ }}}
#Make Sure PSF is centred on centre of stamp /*fold*/ {{{
centre<-psf.cen[[psf.id[ind]]]
delta<-floor(stamplen[ind]/2)*c(-1,+1)
lims<-rbind(centre[1]+delta,centre[2]+delta)
dx<-lims[1,1]-1
dy<-lims[2,1]-1
dx<-ifelse(dx>0,dx+1,dx)
dy<-ifelse(dy>0,dy+1,dy)
#Reinterpolate the PSF at point source XcenYcen /*fold*/ {{{
lenxpsf<-length(psf[[psf.id[ind]]][,1])
lenypsf<-length(psf[[psf.id[ind]]][1,])
lenx<-length(1:stamplen[ind])
leny<-length(1:stamplen[ind])
# /*fend*/ }}}
#Make grid for psf at old pixel centres /*fold*/ {{{
psf.obj<-list(x=seq(1,lenx), y=seq(1,leny),z=psf[[psf.id[ind]]][(1:(lenxpsf+1)+dx-1)%%(lenxpsf+1),(1:(lenypsf+1)+dy-1)%%(lenypsf+1)][1:lenx,1:leny])
# /*fend*/ }}}
# /*fend*/ }}}
#Make expanded grid of new pixel centres /*fold*/ {{{
expanded<-expand.grid(seq(1,lenx),seq(1,leny))
xnew<-expanded[,1]-cat.x[ind]%%1
ynew<-expanded[,2]-cat.y[ind]%%1
# /*fend*/ }}}
#Interpolate /*fold*/ {{{
ap<-matrix(interp.2d(xnew, ynew, psf.obj)[,3], ncol=leny,nrow=lenx)
# /*fend*/ }}}
#Aperture Correction Plot /*fold*/ {{{
PlotDev(file=file.path(path.root,path.work,path.out,paste0("ApertureCorrection.",plot.device)))
Rast<-ifelse(length(ap)>1E4,TRUE,FALSE)
if (!psf.weighted) {
#Binary Aperture /*fold*/ {{{
suppressWarnings(image(log10(ap),main="Example: Minimum Aperture Correction (smallest source)", asp=1,col=heat.colors(256),useRaster=Rast))
contour(ap, levels=(tempfunc(c(0.5,0.9,0.95,0.99,0.999,0.9999)*max(tempsum,na.rm=TRUE))), labels=c(0.5,0.9,0.95,0.99,0.999,0.9999), col='blue', add=TRUE)
contour(ap, levels=c(psfLimit), labels=c(ap.limit), col='green', add=TRUE)
suppressWarnings(image(log10(sfa[[ind]]),col=col2alpha('blue',0.3),add=TRUE,useRaster=Rast))
spsf<-sum(ap)
ssfap<-sum(sfa[[ind]]*ap)
label("topright",lab=paste("SumPSF=",round(spsf,digits=2),"\nSum(PSF*Ap)=",round(ssfap,digits=2),"\nApCorr=",round(spsf/ssfap,digits=2),sep=""))
# /*fend*/ }}}
} else {
#PSF Weighted Aperture /*fold*/ {{{
suppressWarnings(image(log10(ap),main="Example: Minimum Aperture Correction (smallest source)", asp=1,col=heat.colors(256),useRaster=Rast))
contour(ap, levels=(tempfunc(c(0.5,0.9,0.95,0.99,0.999,0.9999)*max(tempsum,na.rm=TRUE))), labels=c(0.5,0.9,0.95,0.99,0.999,0.9999), col='blue', add=TRUE)
suppressWarnings(image(log10(sfa[[ind]]),col=col2alpha('blue',0.3),add=TRUE,useRaster=Rast))
spsf<-sum(ap)
ssfap<-sum(sfa[[ind]]*ap)
label("topright",lab=paste("SumPSF=",round(spsf,digits=2),"\nSum(PSF*Ap)=",round(ssfap,digits=2),"\nApCorr=",round(spsf/ssfap,digits=2),sep=""))
# /*fend*/ }}}
}
# /*fend*/ }}}
#Plot Example of PS minimum aperture Correction; median source /*fold*/ {{{
ind<-(which.min(abs(cat.a-median(cat.a)))[1])
#Make Sure PSF is centred on centre of stamp /*fold*/ {{{
centre<-psf.cen[[psf.id[ind]]]
delta<-floor(stamplen[ind]/2)*c(-1,+1)
lims<-rbind(centre[1]+delta,centre[2]+delta)
dx<-lims[1,1]-1
dy<-lims[2,1]-1
dx<-ifelse(dx>0,dx+1,dx)
dy<-ifelse(dy>0,dy+1,dy)
# /*fend*/ }}}
#Reinterpolate the PSF at point source XcenYcen /*fold*/ {{{
lenxpsf<-length(psf[[psf.id[ind]]][,1])
lenypsf<-length(psf[[psf.id[ind]]][1,])
lenx<-length(1:stamplen[ind])
leny<-length(1:stamplen[ind])
# /*fend*/ }}}
#Make grid for psf at old pixel centres /*fold*/ {{{
psf.obj<-list(x=seq(1,lenx), y=seq(1,leny),z=psf[[psf.id[ind]]][(1:(lenxpsf+1)+dx-1)%%(lenxpsf+1),(1:(lenypsf+1)+dy-1)%%(lenypsf+1)][1:lenx,1:leny])
# /*fend*/ }}}
#Make expanded grid of new pixel centres /*fold*/ {{{
expanded<-expand.grid(seq(1,lenx),seq(1,leny))
xnew<-expanded[,1]-cat.x[ind]%%1
ynew<-expanded[,2]-cat.y[ind]%%1
# /*fend*/ }}}
#Interpolate /*fold*/ {{{
ap<-matrix(interp.2d(xnew, ynew, psf.obj)[,3], ncol=leny,nrow=lenx)
# /*fend*/ }}}
# /*fend*/ }}}
#Aperture Correction Plot /*fold*/ {{{
Rast<-ifelse(length(ap)>1E4,TRUE,FALSE)
if (!psf.weighted) {
#Binary Aperture /*fold*/ {{{
suppressWarnings(image(log10(ap),main="Example: Minimum Aperture Correction (median source)", asp=1,col=heat.colors(256),useRaster=Rast))
contour(ap, levels=(tempfunc(c(0.5,0.9,0.95,0.99,0.999,0.9999)*max(tempsum,na.rm=TRUE))), labels=c(0.5,0.9,0.95,0.99,0.999,0.9999), col='blue', add=TRUE)
if (psf.filt) { contour(ap, levels=c(psfLimit), labels=c(ap.limit), col='green', add=TRUE) }
suppressWarnings(image(log10(sfa[[ind]]),col=col2alpha('blue',0.3),add=TRUE,useRaster=Rast))
spsf<-sum(ap)
ssfap<-sum(sfa[[ind]]*ap)
label("topright",lab=paste("SumPSF=",round(spsf,digits=2),"\nSum(PSF*Ap)=",round(ssfap,digits=2),"\nApCorr=",round(spsf/ssfap,digits=2),sep=""))
# /*fend*/ }}}
} else {
#PSF Weighted Aperture /*fold*/ {{{
suppressWarnings(image(log10(ap)-log10(sfa[[ind]]),main="Example: Minimum Aperture Correction (median source; shown as residual)", asp=1,col=heat.colors(256),useRaster=Rast))
contour(ap, levels=(tempfunc(c(0.5,0.9,0.95,0.99,0.999,0.9999)*max(tempsum,na.rm=TRUE))), labels=c(0.5,0.9,0.95,0.99,0.999,0.9999), col='blue', add=TRUE)
spsf<-sum(ap)
ssfap<-sum(sfa[[ind]]*ap)
label("topright",lab=paste("SumPSF=",round(spsf,digits=2),"\nSum(PSF*Ap)=",round(ssfap,digits=2),"\nApCorr=",round(spsf/ssfap,digits=2),sep=""))
# /*fend*/ }}}
}
if (!grepl('x11',plot.device,ignore.case=TRUE)) { dev.off() }
# /*fend*/ }}}
}
# /*fend*/ }}}
#Generate Deblended Flux Arrays /*fold*/ {{{
if (!quiet) { cat("Generating Deblended Flux Arrays") }
timer<-proc.time()
#Create the Deblended Flux Array /*fold*/ {{{
dfa<-foreach(dbwm=dbw, sfam=sfa, .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% { dbwm*sfam }
# /*fend*/ }}}
#Notify /*fold*/ {{{
if (showtime) { cat(" - Done (",round(proc.time()[3]-timer[3],digits=2),"sec )\n")
message(paste("Make DFA - Done (",round(proc.time()[3]-timer[3],digits=2),"sec )"))
} else if (!quiet) { cat(" - Done\n") }
# /*fend*/ }}}
#Check if we're just calculating the deblend fraction /*fold*/ {{{
if (!exists("get.debl.frac")) { get.debl.frac<-FALSE }
if (get.debl.frac) {
if (!quiet) { cat("Getting Deblend Fraction Only flag is TRUE.\nCalculating Aperture Integrals & Outputting catalogue.\n") }
#Integral of the aperture; ssa /*fold*/ {{{
if (verbose) { cat("Integral of the aperture") }
ssa<-foreach(sam=sa, .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% { sum(sam) }
if (verbose) { cat(" - Done\n") } # /*fend*/ }}}
#Integral of the convolved aperture; ssfa /*fold*/ {{{
if (verbose) { cat("Integral of the convolved aperture") }
ssfa<-foreach(sfam=sfa, .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% { sum(sfam) }
if (verbose) { cat(" - Done\n") } # /*fend*/ }}}
#Integral of the deblended convolved aperture; sdfa /*fold*/ {{{
if (verbose) { cat("Integral of the deblended convolved aperture") }
sdfa<-foreach(dfam=dfa, .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% { sum(dfam) }
if (verbose) { cat(" - Done\n") } # /*fend*/ }}}
spsf<-rep(NA, length(sdfa))
for (i in 1:length(psf)) {
spsf[which(psf.id==i)]<-sumpsf[i]
}
#If wanted, output the Results Table /*fold*/ {{{
if (write.tab) {
#Output the results table /*fold*/ {{{
if ((loop.total!=1)&&(length(param.env$tableout.name)!=loop.total)) {
if (!quiet) { cat(paste('Writing Deblend Fraction Table to ',file.path(path.root,path.work,path.out,paste(sep="",tableout.name,"_file",f,".csv")),' ')) }
timer=system.time(write.deblend.fraction.table(filename=file.path(path.root,path.work,path.out,paste(sep="",tableout.name,"_file",f,".csv"))) )
} else {
if (!quiet) { cat(paste('Writing Deblend Fraction Table to ',file.path(path.root,path.work,path.out,paste(sep="",tableout.name,".csv")),' ')) }
timer=system.time(write.deblend.fraction.table(filename=file.path(path.root,path.work,path.out,paste(sep="",tableout.name,".csv"))) )
}
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(paste('Write Deblend Fraction Table - Done (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
}# /*fend*/ }}}
if (!quiet) { cat("....Done\n") }
return()
}
# /*fend*/ }}}
# /*fend*/ }}}
# /*fend*/ }}}
#If wanted; Iterate the fluxes to improve final measurements /*fold*/ {{{
if (iterate.fluxes) {
#Notify /*fold*/ {{{
message(paste("Iterating Flux Determination",num.iterations,"times {\n"))
if (!quiet) { cat("Iterating Flux Determination",num.iterations,"times {\n") }
# /*fend*/ }}}
#For the number of desired iterations /*fold*/ {{{
#Setup for iteration /*fold*/ {{{
weight.type='flux'
quietbak<-quiet
psf.filtbak<-psf.filt
fluxiters<-matrix(as.numeric(NA),ncol=num.iterations,nrow=length(cat.id))
erriters<-fluxiters
sdfaiters<-fluxiters
#Integral of the convolved aperture; ssfa /*fold*/ {{{
if (verbose) { cat(" Integral of the convolved aperture") }
ssfa<-foreach(sfam=sfa, .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% { sum(sfam) }
if (verbose) { cat(" - Done\n") } # /*fend*/ }}}
#Integral of the deblended convolved aperture; sdfa /*fold*/ {{{
if (verbose) { cat(" Integral of the deblended convolved aperture") }
sdfa<-foreach(dfam=dfa, .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% { sum(dfam) }
if (verbose) { cat(" - Done\n") } # /*fend*/ }}}
sdfad<-sdfa*NA
sdfa2e2<-sdfa*NA
attach(image.env)
#/*fend*/ }}}
for (iter in 1:num.iterations) {
#Determine objects to iterate over /*fold*/ {{{
if (iter!=1 | length(sdfa)==1) {
xind<-which(sdfa!=ssfa & sdfa!=0)
iterateLost<-sdfa==0
} else {
xind<-1:length(sdfa)
}
#Check we have something to do! /*fold*/ {{{
if (length(xind)==0) {
if (verbose) { cat(" - Breaking out; no objects remaining overlap, so no point iterating!\n") }
message("Breaking out of Iteration; no objects remaining overlap, so no point iterating!\n")
break
}
#/*fend*/ }}}
#Update MPI options /*fold*/ {{{
if (num.cores < 0) {
registerDoParallel(cores=(min(ceiling(length(cat.x)/50000),abs(num.cores))))
}
chunk.size=ceiling(length(xind)/getDoParWorkers())
mpi.opts<-list(chunkSize=chunk.size)
#/*fend*/ }}}
#/*fend*/ }}}
#Calculate the flux per object /*fold*/ {{{
#Notify /*fold*/ {{{
message(paste("Calculating Flux (#",iter,")"))
if (!quiet) { cat(" Calculating Flux (#",iter,")") }
# /*fend*/ }}}
timer<-proc.time()
#if (iter==1) {
# #if the first iteration, try to use 'quasi-segmented flux' for 0th step /*fold*/ {{{
# if (cutup) {
# sdfad[xind]<-foreach(dfam=dfa[xind],dbwm=dbw[xind],im=data.stamp[xind], xlo=ap.lims.data.stamp[xind,1],xup=ap.lims.data.stamp[xind,2], ylo=ap.lims.data.stamp[xind,3],yup=ap.lims.data.stamp[xind,4], .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% {
# ind<-which(dbwm >= 0.5,arr.ind=T)
# if (length(ind)<length(dbwm)*0.1) {
# sum((im[xlo:xup,ylo:yup]*dfam))
# } else {
# sum(((im[xlo:xup,ylo:yup]*dfam)[ind]))
# }
# }
# } else {
# sdfad[xind]<-foreach(dfam=dfa[xind],dbwm=dbw[xind],xlo=ap.lims.data.map[xind,1],xup=ap.lims.data.map[xind,2], ylo=ap.lims.data.map[xind,3],yup=ap.lims.data.map[xind,4], .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment()),.export='im') %dopar% {
# ind<-which(dbwm >= 0.5,arr.ind=T)
# if (length(ind)<length(dbwm)*0.1) {
# sum((im[xlo:xup,ylo:yup]*dfam))
# } else {
# sum(((im[xlo:xup,ylo:yup]*dfam)[ind]))
# }
# }
# }
# # /*fend*/ }}}
#} else {
#if not the first iteration, use weighted deblended flux /*fold*/ {{{
if (cutup) {
sdfad[xind]<-foreach(dfam=dfa[xind],im=data.stamp[xind], xlo=ap.lims.data.stamp[xind,1],xup=ap.lims.data.stamp[xind,2], ylo=ap.lims.data.stamp[xind,3],yup=ap.lims.data.stamp[xind,4], .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% {
sum(dfam*(im[xlo:xup,ylo:yup]))
}
} else {
sdfad[xind]<-foreach(dfam=dfa[xind], xlo=ap.lims.data.map[xind,1],xup=ap.lims.data.map[xind,2], ylo=ap.lims.data.map[xind,3],yup=ap.lims.data.map[xind,4], .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment()),.export='im') %dopar% {
sum(dfam*(im[xlo:xup,ylo:yup]))
}
}
#}
#Notify /*fold*/ {{{
if (showtime) { cat(" - Done (",round(proc.time()[3]-timer[3],digits=2),"sec )\n")
message(" - Done (",round(proc.time()[3]-timer[3],digits=2),"sec )\n")
} else if (!quiet) { cat(" - Done\n") }
# /*fend*/ }}}
# /*fend*/ }}}
#Calculate the error per object /*fold*/ {{{
#Notify /*fold*/ {{{
message(paste("Calculating Error (#",iter,")"))
if (!quiet) { cat(" Calculating Error (#",iter,")") }
timer<-proc.time()
# /*fend*/ }}}
if (length(error.stamp)>1|(length(error.stamp)==1 & is.list(error.stamp))) {
sdfa2e2[xind]<-foreach(dfam=dfa[xind], xlo=ap.lims.error.stamp[xind,1],xup=ap.lims.error.stamp[xind,2],ylo=ap.lims.error.stamp[xind,3],yup=ap.lims.error.stamp[xind,4], ime=error.stamp[xind], .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% {
sum((dfam*ime[xlo:xup,ylo:yup])^2.)
}
} else if (length(error.stamp)==1){
if (error.stamp==1) {
sdfa2e2[xind]<-foreach(dfam=dfa[xind], .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% {
sum((dfam)^2)
}
} else if (error.stamp==0) {
sdfa2e2<-rep(0,length(sdfad))
} else {
sdfa2e2[xind]<-foreach(dfam=dfa[xind], .noexport=ls(envir=environment()), .export="error.stamp", .inorder=TRUE, .combine='c', .options.mpi=mpi.opts) %dopar% {
sum((dfam*error.stamp)^2.)
}
}
} else {
sdfa2e2[xind]<-foreach(dfam=dfa[xind], xlo=ap.lims.error.map[xind,1],xup=ap.lims.error.map[xind,2],ylo=ap.lims.error.map[xind,3],yup=ap.lims.error.map[xind,4], .noexport=ls(envir=environment()), .export='ime', .inorder=TRUE, .combine='c', .options.mpi=mpi.opts) %dopar% {
sum((dfam*ime[xlo:xup,ylo:yup])^2.)
}
}
#Notify /*fold*/ {{{
if (showtime) { cat(" - Done (",round(proc.time()[3]-timer[3],digits=2),"sec )\n")
message(" - Done (",round(proc.time()[3]-timer[3],digits=2),"sec )\n")
} else if (!quiet) { cat(" - Done\n") }
# /*fend*/ }}}
# /*fend*/ }}}
#If wanted, remove sky estimates /*fold*/ {{{
if (!quiet) { message(paste("Calculating Deblend Fraction (#",iter,")")); cat(paste(" Calculating Deblend Fraction (#",iter,")")) }
timer<-proc.time()
sdfa[xind]<-foreach(dfam=dfa[xind], .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% { sum(dfam) }
#Subtract Sky Flux /*fold*/ {{{
if (do.sky.est) {
sdfad[xind]<-sdfad[xind]-skylocal[xind]*sdfa[xind]
}
if (showtime) { cat(" - Done (",round(proc.time()[3]-timer[3],digits=2),"sec )\n")
message(paste('Sky Estimate - Done (',round(proc.time()[3]-timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
# /*fend*/ }}}
# /*fend*/ }}}
#Save the values of the fluxes for output /*fold*/ {{{
fluxiters[,iter]<-sdfad
erriters[,iter]<-sqrt(sdfa2e2)
sdfaiters[,iter]<-sdfa
# /*fend*/ }}}
#Re-calculate the weighted apertures /*fold*/ {{{
#Notify /*fold*/ {{{
message(paste("Calculating Weighting Apertures (#",iter,")"))
if (!quiet) { cat(" Calculating Weighting Apertures (#",iter,")") }
# /*fend*/ }}}
quiet<-TRUE
#if (psf.weighted) {
# #If using psf.weighted, we may be able to skip an unnecessary convolution /*fold*/ {{{
psf.filt<-FALSE
timer=system.time(wsfa[xind]<-make.convolved.apertures(outenv=environment(), sfabak,flux.weightin=sdfad,sumsa=ssfa,subs=xind))
# # /*fend*/ }}}
#} else {
# #If not using psf.weighted, we may need to re-convolve /*fold*/ {{{
# psf.filt<-psf.filtbak
# timer=system.time(wsfa[xind]<-make.convolved.apertures(outenv=environment(), sa,flux.weightin=sdfad,subs=xind))
# # /*fend*/ }}}
#}
timer2=system.time(image.env$wfa<-make.aperture.map(outenv=environment(), wsfa, dimim,subs=xind))
psf.filt<-psf.filtbak
quiet<-quietbak
#Notify /*fold*/ {{{
if (showtime) { cat(" - Done (",round(timer[3]+timer2[3],digits=2),"sec )\n")
message(" - Done (",round(proc.time()[3]-timer[3],digits=2),"sec )\n")
} else if (!quiet) { cat(" - Done\n") }
# /*fend*/ }}}
# /*fend*/ }}}
#Re-calculate the deblending matricies /*fold*/ {{{
#Notify /*fold*/ {{{
message(paste("Calculating Deblend Matricies (#",iter,")"))
if (!quiet) { cat(" Calculating Deblend Matricies (#",iter,")") }
# /*fend*/ }}}
quiet<-TRUE
timer2=system.time(dbw[xind]<-make.deblend.weight.maps(outenv=environment(),subs=xind))
quiet<-quietbak
timer<-proc.time()
dfa[xind]<-foreach(dbwm=dbw[xind], sfam=sfa[xind], .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% { dbwm*sfam }
#Notify /*fold*/ {{{
if (showtime) { cat(" - Done (",round(proc.time()[3]-timer[3]+timer2[3],digits=2),"sec )\n")
message(" - Done (",round(proc.time()[3]-timer[3],digits=2),"sec )\n")
} else if (!quiet) { cat(" - Done\n") }
# /*fend*/ }}}
# /*fend*/ }}}
}
#For objects completely removed, make final measured flux the one that is output /*fold*/ {{{
fluxiters[which(fluxiters==0)]<-NA
erriters[which(erriters==0)]<-NA
sdfaiters[which(sdfaiters==0)]<-NA
for (i in 1:length(fluxiters[,1])) {
fluxiters[i,which(is.na(fluxiters[i,]))]<-rev(fluxiters[i,which(!is.na(fluxiters[i,]))])[1]
erriters[i,which(is.na(erriters[i,]))]<-rev(erriters[i,which(!is.na(erriters[i,]))])[1]
sdfaiters[i,which(is.na(sdfaiters[i,]))]<-rev(sdfaiters[i,which(!is.na(sdfaiters[i,]))])[1]
}
#Update MPI options /*fold*/ {{{
if (num.cores < 0) {
registerDoParallel(cores=(min(ceiling(length(cat.x)/50000),abs(num.cores))))
}
chunk.size=ceiling(length(cat.id)/getDoParWorkers())
mpi.opts<-list(chunkSize=chunk.size)
#/*fend*/ }}}
#/*fend*/ }}}
detach(image.env)
if (!quiet) { cat("} Done\n") }
# /*fend*/ }}}
# /*fend*/ }}}
}
# /*fend*/ }}}
#Integral of the deblended convolved aperture; sdfa /*fold*/ {{{
if (verbose) { cat("Measuring blend fractions") }
sdfa<-foreach(dfam=dfa, .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% { sum(dfam) }
ssfa<-foreach(sfam=sfa, .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% { sum(sfam) }
if (verbose) { cat(" - Done\n") } # /*fend*/ }}}
if (psf.check) {
#Estimate the PSF from the image, and compare to that given /*fold*/ {{{
if (!quiet) { cat(paste('Estimating the PSF from the image')) }
#Estimate the PSF {{{
if (plot.sample) { pdf(height=10,width=10,file=file.path(path.root,path.work,path.out,'PSFEst_Samples.pdf')) }
timer=system.time(psf.est<-estimate.psf(outenv=environment(),plot=plot.sample,radial.tolerance=radial.tolerance,n.sources=n.sources))
epsf.id<-tmp.psf.id
psfest.val<-tmp.psfest.val
if (plot.sample & !grepl('x11',plot.device,ignore.case=TRUE)) { dev.off() }
#}}}
estpsf.cen<-estpsf.plot<-estpsf2<-estpsf<-estpsf.warn<-estpsf.warnt<-list(NULL)
sumepsf<-rep(NA,length(psf.est$WEIGHT))
warn<-FALSE
for (i in 1:length(psf.est$WEIGHT)) {
if (!all(psf.est$WEIGHT[[i]]==0)) {
#If the estimate succeeded {{{
#Normalise the estimate {{{
estpsf[[i]]<-psf.est$PSF[[i]]
estpsf[[i]]<-estpsf[[i]]-min(estpsf[[i]],na.rm=TRUE)
estpsf[[i]]<-estpsf[[i]]/max(estpsf[[i]],na.rm=TRUE)
#}}}
#Truncate any upturn in the PSF {{{
if (plot.sample) {
PlotDev(file=file.path(path.root,path.work,path.out,paste0("PSFEst_truncation_",i,".",plot.device)),width=12,height=6,units="in")
layout(matrix(1:6,nrow=2,byrow=T))
}
trunc<-truncate.upturn(estpsf[[i]],plot=plot.sample,centre=psf.est$centre)
estpsf2[[i]]<-trunc$Im
estpsf.cen[[i]]<-trunc$centre
estpsf2[[i]]<-estpsf2[[i]]-min(estpsf2[[i]],na.rm=TRUE)
estpsf2[[i]]<-estpsf2[[i]]/max(estpsf2[[i]],na.rm=TRUE)
sumepsf[i]<-sum(estpsf2[[i]])
big<-matrix(0,ncol=max(stamplen),nrow=max(stamplen))
big[1:dim(estpsf2[[i]])[1],1:dim(estpsf2[[i]])[2]]<-estpsf2[[i]]
estpsf2[[i]]<-big
if (plot.sample & !grepl('x11',plot.device,ignore.case=TRUE)) { dev.off() }
#}}}
#Recentre the PSF {{{
conv<-array(0,dim=dim(estpsf2[[i]]))
conv[dim(conv)[1]/2,dim(conv)[2]/2]<-1
estpsf.plot[[i]]<-convolve.psf(estpsf2[[i]],conv)
estpsf.plot[[i]]<-estpsf.plot[[i]]-min(estpsf.plot[[i]],na.rm=TRUE)
estpsf.plot[[i]]<-estpsf.plot[[i]]/max(estpsf.plot[[i]],na.rm=TRUE)
big<-matrix(0,ncol=max(stamplen),nrow=max(stamplen))
big[1:dim(estpsf.plot[[i]])[1],1:dim(estpsf.plot[[i]])[2]]<-estpsf.plot[[i]]
estpsf.plot[[i]]<-big
#}}}
if (!no.psf && any(psf.id>=i)) {
#We have a PSF that we can compare to: {{{
#Centre the input PSF {{{
if (all(dim(psf[[i]]) > dim(estpsf.plot[[i]]))) {
#PSF is larger than estpsf in all dims
cen<-which(psf[[i]]==max(psf[[i]]),arr.ind=T)
#Check that centre isnt less than dim(estpsf)/2
for (j in 1:2) {
if (cen[j] < dim(estpsf.plot[[i]])[j]) {
cen[j]<-dim(estpsf.plot[[i]])[j]/2+1
}
}
#Get psf indicies
psf.ind.x<-(1:dim(estpsf.plot[[i]])[1])+cen[1]-dim(estpsf.plot[[i]])[1]/2
psf.ind.y<-(1:dim(estpsf.plot[[i]])[2])+cen[2]-dim(estpsf.plot[[i]])[2]/2
#Get PSF plot
psf.plot<-convolve.psf(psf[[i]][psf.ind.x,psf.ind.y],conv)
} else if (any(dim(psf[[i]]) < dim(estpsf.plot[[i]]))) {
#PSF is smaller than estpsf in one dim
psf.plot<-array(0,dim=dim(estpsf.plot[[i]]))
cen<-which(psf[[i]]==max(psf[[i]]),arr.ind=T)
j<-which(dim(psf[[i]]) > dim(estpsf.plot[[i]]))
if (cen[j] < dim(estpsf.plot[[i]])[j]) {
cen[j]<-dim(estpsf.plot[[i]])[j]/2+1
}
psf.ind<-seq(cen[j]-dim(estpsf.plot[[i]])[j]/2,cen[j]+dim(estpsf.plot[[i]])[j]/2)
if (j==1) {
psf.plot[psf.ind,1:dim(psf)[2]]<-psf[[i]]
} else {
psf.plot[1:dim(psf)[1],psf.ind]<-psf[[i]]
}
psf.plot<-convolve.psf(psf[[i]][psf.ind.x,psf.ind.y],conv)
} else {
#PSF is smaller than or equal to estpsf in all dims
conv<-psf.plot<-array(0,dim=dim(estpsf.plot[[i]]))
conv[dim(conv)[1]/2,dim(conv)[2]/2]<-1
psf.plot[1:dim(psf[[i]])[1],1:dim(psf[[i]])[2]]<-psf[[i]]
psf.plot<-convolve.psf(psf.plot,conv)
}
psf.plot<-psf.plot-min(psf.plot,na.rm=TRUE)
psf.plot<-psf.plot/max(psf.plot,na.rm=TRUE)
#}}}
#Plot a sample of the PSF {{{
if (plot.sample) {
#Open the device {{{
PlotDev(file=file.path(path.root,path.work,path.out,paste0("PSF_residuals_",i,".",plot.device)),width=12,height=6,units="in")
layout(matrix(1:6,byrow=T,nrow=2))
#}}}
#Plots {{{
#Plot the input PSF
cen<-which(psf[[i]] == max(psf[[i]]), arr.ind = T)
x.range<-range(which(psf[[i]][cen[1],]!=0))
y.range<-range(which(psf[[i]][,cen[2]]!=0))
capture<-magimage(psf[[i]][x.range[1]:x.range[2],y.range[1]:y.range[2]],main='Input')
#Plot the measured PSF
cen<-which(estpsf2[[i]] == max(estpsf2[[i]]), arr.ind = T)
x.range<-range(which(estpsf2[[i]][cen[1],]!=0))
y.range<-range(which(estpsf2[[i]][,cen[2]]!=0))
capture<-magimage(estpsf2[[i]][x.range[1]:x.range[2],y.range[1]:y.range[2]],main='Estimated')
#Plot the Residual PSF
cen<-which(psf.plot == max(psf.plot), arr.ind = T)
x.range<-range(which(psf.plot[cen[1],]!=0))
y.range<-range(which(psf.plot[,cen[2]]!=0))
capture<-magimage((psf.plot-estpsf.plot[[i]])[x.range[1]:x.range[2],y.range[1]:y.range[2]],main='Residual',lo=-0.3,hi=0.3,type='num')
#Plot the profile of the input PSF
magplot(estpsf.plot[[i]][which(estpsf.plot[[i]]==max(estpsf.plot[[i]],na.rm=TRUE),arr.ind=T)[1],y.range[1]:y.range[2]],type='s',col='grey',main='Input')
lines(estpsf.plot[[i]][x.range[1]:x.range[2],which(estpsf.plot[[i]]==max(estpsf.plot[[i]],na.rm=TRUE),arr.ind=T)[2]],type='s',col='grey')
lines(psf.plot[which(psf.plot==max(psf.plot,na.rm=TRUE),arr.ind=T)[1],y.range[1]:y.range[2]],type='s',col='red')
lines(psf.plot[x.range[1]:x.range[2],which(psf.plot==max(psf.plot,na.rm=TRUE),arr.ind=T)[2]],type='s',col='blue')
legend('topright',inset=0.1,bty='n',legend=c('x profile','y profile','estimated PSF'),col=c('red','blue','grey'),lty=1,pch=NA)
#Plot the profile of the PSF estimate
magplot(psf.plot[which(psf.plot==max(psf.plot,na.rm=TRUE),arr.ind=T)[1],y.range[1]:y.range[2]],type='s',col='grey',main='Estimated')
lines(psf.plot[x.range[1]:x.range[2],which(psf.plot==max(psf.plot,na.rm=TRUE),arr.ind=T)[2]],type='s',col='grey')
lines(estpsf.plot[[i]][which(estpsf.plot[[i]]==max(estpsf.plot[[i]],na.rm=TRUE),arr.ind=T)[1],y.range[1]:y.range[2]],type='s',col='red')
lines(estpsf.plot[[i]][x.range[1]:x.range[2],which(estpsf.plot[[i]]==max(estpsf.plot[[i]],na.rm=TRUE),arr.ind=T)[2]],type='s',col='blue')
legend('topright',inset=0.1,bty='n',legend=c('x profile','y profile','input PSF'),col=c('red','blue','grey'),lty=1,pch=NA)
#Plot the residual profile of the PSFs
magplot(psf.plot[which(psf.plot==max(psf.plot,na.rm=TRUE),arr.ind=T)[1],y.range[1]:y.range[2]]-
estpsf.plot[[i]][which(estpsf.plot[[i]]==max(estpsf.plot[[i]],na.rm=TRUE),arr.ind=T)[1],y.range[1]:y.range[2]],type='s',col='red',
main='Residual',ylim=c(-0.2,0.2))
lines(psf.plot[x.range[1]:x.range[2],which(psf.plot==max(psf.plot,na.rm=TRUE),arr.ind=T)[2]]-
estpsf.plot[[i]][x.range[1]:x.range[2],which(estpsf.plot[[i]]==max(estpsf.plot[[i]],na.rm=TRUE),arr.ind=T)[2]],type='s',col='blue')
legend('topright',inset=0.1,bty='n',legend=c('x profile','y profile'),lty=1,col=c('red','blue'),pch=NA)
label('bottomleft',lab=paste0('sum(abs(resid))/sum(abs(psf)=',round(sum(abs(psf.plot-estpsf.plot[[i]]))/sum(abs(psf.plot)),digits=3)))
#}}}
#Close the device {{{
if (!grepl('x11',plot.device,ignore.case=TRUE)) { dev.off() }
#}}}
} #}}}
#Check the quality of the input PSF compared to the measured {{{
stats<-summary(as.numeric(psf.plot-estpsf.plot[[i]]))
message("PSF Estimate Properties:")
message(paste('Peak PSF residual: ',round(stats[5], digits=4),' (NB: PSFs peak at 1.0)'))
message(paste('Median PSF residual: ',round(stats[3], digits=4)))
message(paste('25 & 75 % PSF residual: ',round(stats[2], digits=4),'-',round(stats[4], digits=4)))
estpsf.warnt[[i]]<-""
estpsf.warn[[i]]<-0
if (stats[5] > 0.2) {
estpsf.warn[[i]]<-estpsf.warn[[i]]+1
estpsf.warnt[[i]]<-paste0(estpsf.warnt[[i]],"P")
}
if (stats[3] > 0.2) {
estpsf.warn[[i]]<-estpsf.warn[[i]]+2
estpsf.warnt[[i]]<-paste0(estpsf.warnt[[i]],"M")
}
if (any(which(psf[[i]]==max(psf[[i]]),arr.ind=T) != psf.cen[[i]])) {
estpsf.warn[[i]]<-estpsf.warn[[i]]+4
estpsf.warnt[[i]]<-paste0(estpsf.warnt[[i]],"C")
} else if (any(which(estpsf[[i]]==max(estpsf[[i]]),arr.ind=T) != estpsf.cen[[i]])) {
estpsf.warn[[i]]<-estpsf.warn[[i]]+4
estpsf.warnt[[i]]<-paste0(estpsf.warnt[[i]],"C")
}
#}}}
#}}}
} else {
estpsf.warn[[i]]<- -1
estpsf.warnt[[i]]<- ""
}
#}}}
} else {
#PSF estimate failed; warn and continue {{{
warn<-TRUE
if (i==min(epsf.id)) {
if (i==length(estpsf.plot)) {
#All psfs have failed, leave the id as it is.
} else {
#If this is the lowest psf.id, give
#these sources the id of the next psf up
epsf.id[which(epsf.id==i)]<-min(epsf.id)+1
}
} else {
#If this is not the lowest psf.id, give
#these sources the id of the last psf
epsf.id[which(epsf.id==i)]<-max(epsf.id[which(epsf.id<i)])
}
estpsf.plot[[i]]<-estpsf2[[i]]<-estpsf[[i]]<-NA
estpsf.warn[[i]]<-8
estpsf.warnt[[i]]<-"E"
#}}}
}
}
if (warn) {
if (showtime) { cat(" - WARNING: PSF estimate failed; no suitable point sources found! (",round(timer[3],digits=2),"sec )\n")
message(" - WARNING: PSF estimate failed; no suitable point sources found! (",round(timer[3],digits=2),"sec )\n")
} else if (!quiet) { cat(" - WARNING: PSF estimate failed; no suitable point sources found!\n") }
} else {
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(" - Done (",round(timer[3],digits=2),"sec )\n")
} else if (!quiet) { cat(" - Done\n") }
}
#Save the PSF as output {{{
psflist=list(recentred=estpsf.plot,final=estpsf2,normalised=estpsf,raw=psf.est)
save(file=file.path(path.root,path.work,path.out,"PSF_estimate.Rdata"),psflist)
#}}}
} else {
estpsf.warn<- list(8)
estpsf.warnt<- list("")
sumepsf<-NA
}
# /*fend*/ }}}
#If wanted, output the Deblended & Convolved Aperture Mask /*fold*/ {{{
if (make.debelended.apertures.map) {
if (!quiet) { cat(paste('Making All Deblended Convolved Apertures Mask - ')) }
timer=system.time(image.env$adfa<-make.aperture.map(outenv=environment(), dfa, dimim))
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(paste('Make ADFA Mask - Done (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
if (!quiet) { cat(paste('Outputting All Deblended Convolved Apertures Mask to',dfa.filename," ")) }
timer=system.time(write.fits.image.file(file.path(path.root,path.work,path.out,dfa.filename),image.env$adfa,image.env$data.hdr,nochange=TRUE) )
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(paste('Output ADFA Mask - Done (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
}# /*fend*/ }}}
#Remove array that is no longer needed /*fold*/ {{{
if (!plot.sample) { rm(dbw) }
rm(adfa, envir=image.env)
gc()
# /*fend*/ }}}
#If wanted, Perform Randoms Correction /*fold*/ {{{
if (!exists("ran.cor")) { ran.cor<-FALSE }
if (ran.cor=="execute") {
.executeran.cor()
ran.cor<-TRUE
}
if (ran.cor) {
if (!quiet) { message("Perfoming Randoms Correction"); cat("Performing Randoms Correction") }
#Randoms Estimate is NP hard, so set the processors to Max
if (num.cores < 0) { registerDoParallel(cores=abs(num.cores)) }
if (cutup) {
timer<-system.time(randoms<-ran.cor(data.stamp=data.stamp,mask.stamp=mask.stamp,ap.stamp=sfa,ap.stamp.lims=ap.lims.data.stamp,numIters=num.randoms,mpi.opts=mpi.opts,rem.mask=FALSE))
} else {
timer<-system.time(randoms<-ran.cor(data.stamp=image.env$im[data.stamp.lims[1,1]:data.stamp.lims[1,2],data.stamp.lims[1,3]:data.stamp.lims[1,4]],rand.x=cat.x,rand.y=cat.y,
mask.stamp=image.env$imm[mask.stamp.lims[1,1]:mask.stamp.lims[1,2],mask.stamp.lims[1,3]:mask.stamp.lims[1,4]],
ap.stamp=sfa,ap.stamp.lims=ap.lims.data.map,numIters=num.randoms,mpi.opts=mpi.opts,rem.mask=FALSE))
}
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(paste('Randoms Correction - Done (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
if (plot.sample) {
if (!quiet) { message("Plotting Randoms Correction"); cat("Plotting Randoms Correction") }
if (cutup) {
timer<-system.time(plot.ran.cor(cat.id=cat.id,cat.x=cat.x,cat.y=cat.y,data.stamp=data.stamp,mask.stamp=mask.stamp,ap.stamp=sfa,ap.stamp.lims=ap.lims.data.stamp,data.stamp.lims=data.stamp.lims,numIters=num.randoms,rem.mask=FALSE,path=file.path(path.root,path.work,path.out),plot.sci=plot.sci,contams=contams,plot.all=plot.all,toFile=TRUE,plot.device=plot.device))
} else {
timer<-system.time(plot.ran.cor(cat.id=cat.id,cat.x=cat.x,cat.y=cat.y,data.stamp=image.env$im[data.stamp.lims[1,1]:data.stamp.lims[1,2],data.stamp.lims[1,3]:data.stamp.lims[1,4]],rand.x=cat.x,rand.y=cat.y,
mask.stamp=image.env$imm[mask.stamp.lims[1,1]:mask.stamp.lims[1,2],mask.stamp.lims[1,3]:mask.stamp.lims[1,4]],ap.stamp=sfa,ap.stamp.lims=ap.lims.data.map,data.stamp.lims=data.stamp.lims,mask.stamp.lims=ap.lims.mask.map,numIters=num.randoms,rem.mask=FALSE,path=file.path(path.root,path.work,path.out),plot.sci=plot.sci,contams=contams,plot.all=plot.all,toFile=TRUE,plot.device=plot.device))
}
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(paste('Plotting Randoms Correction - Done (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
}
}
# /*fend*/ }}}
#If wanted, Perform Blanks Correction /*fold*/ {{{
if (!exists("blank.cor")) { blank.cor<-FALSE }
if (blank.cor) {
if (!quiet) { message("Perfoming Blanks Correction"); cat("Performing Blanks Correction") }
#Blanks Estimate is NP hard, so set the processors to Max
if (num.cores < 0) { registerDoParallel(cores=abs(num.cores)) }
if (cutup) {
timer<-system.time(blanks<-ran.cor(data.stamp=data.stamp,mask.stamp=mask.stamp,ap.stamp=sfa,ap.stamp.lims=ap.lims.data.stamp,numIters=num.blanks,mpi.opts=mpi.opts,rem.mask=TRUE))
} else {
timer<-system.time(blanks<-ran.cor(data.stamp=image.env$im[data.stamp.lims[1,1]:data.stamp.lims[1,2],data.stamp.lims[1,3]:data.stamp.lims[1,4]],
mask.stamp=image.env$imm[mask.stamp.lims[1,1]:mask.stamp.lims[1,2],mask.stamp.lims[1,3]:mask.stamp.lims[1,4]],ap.stamp=sfa,ap.stamp.lims=ap.lims.data.map,numIters=num.blanks,mpi.opts=mpi.opts,rem.mask=TRUE))
}
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(paste('Blanks Correction - Done (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
if (plot.sample) {
if (!quiet) { message("Plotting Blanks Correction"); cat("Plotting Blanks Correction") }
if (cutup) {
timer<-system.time(plot.ran.cor(cat.id=cat.id,cat.x=cat.x,cat.y=cat.y,data.stamp=data.stamp,mask.stamp=mask.stamp,ap.stamp=sfa,ap.stamp.lims=ap.lims.data.stamp,data.stamp.lims=data.stamp.lims,numIters=num.blanks,rem.mask=TRUE,path=file.path(path.root,path.work,path.out),plot.sci=plot.sci,contams=contams,plot.all=plot.all,toFile=TRUE,plot.device=plot.device))
} else {
timer<-system.time(plot.ran.cor(cat.id=cat.id,cat.x=cat.x,cat.y=cat.y,data.stamp=image.env$im[data.stamp.lims[1,1]:data.stamp.lims[1,2],data.stamp.lims[1,3]:data.stamp.lims[1,4]],
mask.stamp=image.env$imm[mask.stamp.lims[1,1]:mask.stamp.lims[1,2],mask.stamp.lims[1,3]:mask.stamp.lims[1,4]],ap.stamp=sfa,ap.stamp.lims=ap.lims.data.map,data.stamp.lims=data.stamp.lims,mask.stamp.lims=ap.lims.mask.map,numIters=num.blanks,rem.mask=TRUE,path=file.path(path.root,path.work,path.out),plot.sci=plot.sci,contams=contams,plot.all=plot.all,toFile=TRUE,plot.device=plot.device))
}
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(paste('Plotting Blanks Correction - Done (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
}
}
# /*fend*/ }}}
# /*fend*/ }}}
#PART FOUR: COMPUTE AND OUTPUT GALAXY-BY-GALAXY RESULTS /*fold*/ {{{
if (!quiet) { cat("Computing Galaxy-by-Galaxy Results {\n") }
# Make Images Available /*fold*/ {{{
if (!cutup) { attach(image.env) }
# /*fend*/ }}}
#-----Diagnostic-----# /*fold*/ {{{
if (diagnostic){
if (length(which(is.na(image.env$im)))!=0) { message(paste("Input Parameter 'im' contains NA elements")) }
if (length(which(is.na(error.stamp)))!=0) { message(paste("Input Parameter 'ime' contains NA elements")) }
if (length(which(is.na(sfa)))!=0) { message(paste("Input Parameter 'sfa' contains NA elements")) }
if (length(which(is.na(dfa)))!=0) { message(paste("Input Parameter 'dfa' contains NA elements")) }
if (length(which(is.na(flux.corr)))!=0) { message(paste("Input Parameter 'flux.corr' contains NA elements")) }
}# /*fend*/ }}}
#Perform Calculations /*fold*/ {{{
if (!quiet) { cat(" Calculating Fluxes ") }
if (verbose) { cat("{\n") }
timer<-proc.time()
#-----
#Check for Saturations; saturated /*fold*/ {{{
if (is.finite(saturation)) {
if (verbose) { cat(" Checking for Saturations ") }
if (cutup) {
saturated<-foreach(sfam=sfa, im=data.stamp, xlo=ap.lims.data.stamp[,1],xup=ap.lims.data.stamp[,2], ylo=ap.lims.data.stamp[,3],yup=ap.lims.data.stamp[,4], .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment()), .export='saturation') %dopar% {
any(im[xlo:xup,ylo:yup][which(sfam!=0,arr.ind=TRUE)]>=saturation)
}
} else {
saturated<-foreach(sfam=sfa, .noexport=ls(envir=environment()), .export=c('im','saturation'), xlo=ap.lims.data.map[,1],xup=ap.lims.data.map[,2], ylo=ap.lims.data.map[,3],yup=ap.lims.data.map[,4], .inorder=TRUE, .combine='c', .options.mpi=mpi.opts) %dopar% {
any(im[xlo:xup,ylo:yup][which(sfam!=0,arr.ind=TRUE)]>=saturation)
}
}
if (verbose) { cat(" - Done\n") }
} else {
saturated<-rep(FALSE,length(cat.id))
}
# /*fend*/ }}}
#-----
#Image Flux at central pixel; pixflux /*fold*/ {{{
if (!use.pixel.fluxweight) {
if (verbose) { cat(" Image Flux at central pixel ") }
if (cutup) {
pixflux<-foreach(xp=x.pix-(data.stamp.lims[,1]-1),yp=y.pix-(data.stamp.lims[,3]-1), im=data.stamp, .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% {
im[xp,yp]
}
} else {
pixflux<-foreach(xp=x.pix,yp=y.pix, .noexport=ls(envir=environment()), .export=c('im'), .inorder=TRUE, .combine='c', .options.mpi=mpi.opts) %dopar% {
im[xp,yp]
}
}
if (verbose) { cat(" - Done\n") }
}
# /*fend*/ }}}
#-----
#Integral of the aperture; ssa /*fold*/ {{{
if (verbose) { cat(" Integral of the aperture") }
ssa<-foreach(sam=sa, .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% { sum(sam) }
if (verbose) { cat(" - Done\n") } # /*fend*/ }}}
##-----
# #Integral of the convolved aperture; ssfa /*fold*/ {{{
# if (verbose) { cat(" Integral of the convolved aperture") }
# ssfa<-foreach(sfam=sfa, .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% { sum(sfam) }
# if (verbose) { cat(" - Done\n") } # /*fend*/ }}}
#-----
#Integral of the [(unmodified convolved aperture)*(modified convolved aperture)]; ssfau /*fold*/ {{{
if (verbose) { cat(" Integral of the [(unmodified convolved aperture)*(convolved aperture)]") }
ssfau<-foreach(sfam=sfa, usfam=sfabak, .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% { sum(sfam*usfam) }
if (verbose) { cat(" - Done\n") } # /*fend*/ }}}
#-----
#Integral of the (convolved aperture * image); ssfad /*fold*/ {{{
if (verbose) { cat(" Integral of the (convolved aperture * image)") }
if (cutup) {
ssfad<-foreach(sfam=sfa, im=data.stamp, xlo=ap.lims.data.stamp[,1],xup=ap.lims.data.stamp[,2], ylo=ap.lims.data.stamp[,3],yup=ap.lims.data.stamp[,4], .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% {
sum(sfam*(im[xlo:xup,ylo:yup]))
}
} else {
ssfad<-foreach(sfam=sfa, .noexport=ls(envir=environment()), .export='im', xlo=ap.lims.data.map[,1],xup=ap.lims.data.map[,2], ylo=ap.lims.data.map[,3],yup=ap.lims.data.map[,4], .inorder=TRUE, .combine='c', .options.mpi=mpi.opts) %dopar% {
sum(sfam*(im[xlo:xup,ylo:yup]))
}
}
if (verbose) { cat(" - Done\n") } # /*fend*/ }}}
#-----
#Quartered Photometry - Quartered Integral of the (convolved aperture * image); qssfad /*fold*/ {{{
if (verbose) { cat(" Quartered Integral of the (convolved aperture * image)") }
if (cutup) {
qssfad<-foreach(sfam=sfa, im=data.stamp, x1=ap.lims.data.stamp[,1], x2=ap.lims.data.stamp[,1]+floor((ap.lims.data.stamp[,2]-ap.lims.data.stamp[,1])/2), x3=ap.lims.data.stamp[,1]+floor((ap.lims.data.stamp[,2]-ap.lims.data.stamp[,1])/2)+1, x4=ap.lims.data.stamp[,2],
y1=ap.lims.data.stamp[,3], y2=ap.lims.data.stamp[,3]+floor((ap.lims.data.stamp[,4]-ap.lims.data.stamp[,3])/2), y3=ap.lims.data.stamp[,3]+floor((ap.lims.data.stamp[,4]-ap.lims.data.stamp[,3])/2)+1, y4=ap.lims.data.stamp[,4],
sx1=rep(1,length(stamplen)),sx2=1+floor(stamplen-1)/2,sx3=1+floor(stamplen-1)/2+1,sx4=stamplen,
sy1=rep(1,length(stamplen)),sy2=1+floor(stamplen-1)/2,sy3=1+floor(stamplen-1)/2+1,sy4=stamplen,
.inorder=TRUE, .options.mpi=mpi.opts, .combine="rbind", .noexport=ls(envir=environment())) %dopar% {
cbind(sum(sfam[sx1:sx2,sy1:sy2]*(im[x1:x2,y1:y2])),sum(sfam[sx3:sx4,sy1:sy2]*(im[x3:x4,y1:y2])),sum(sfam[sx1:sx2,sy3:sy4]*(im[x1:x2,y3:y4])),sum(sfam[sx3:sx4,sy3:sy4]*(im[x3:x4,y3:y4])))
}
} else {
qssfad<-foreach(sfam=sfa, .noexport=ls(envir=environment()), .export='im', x1=ap.lims.data.map[,1], x2=ap.lims.data.map[,1]+floor((ap.lims.data.map[,2]-ap.lims.data.map[,1])/2), x3=ap.lims.data.map[,1]+floor((ap.lims.data.map[,2]-ap.lims.data.map[,1])/2)+1, x4=ap.lims.data.map[,2],
y1=ap.lims.data.map[,3], y2=ap.lims.data.map[,3]+floor((ap.lims.data.map[,4]-ap.lims.data.map[,3])/2), y3=ap.lims.data.map[,3]+floor((ap.lims.data.map[,4]-ap.lims.data.map[,3])/2)+1, y4=ap.lims.data.map[,4],
sx1=rep(1,length(stamplen)),sx2=1+floor(stamplen-1)/2,sx3=1+floor(stamplen-1)/2+1,sx4=stamplen,
sy1=rep(1,length(stamplen)),sy2=1+floor(stamplen-1)/2,sy3=1+floor(stamplen-1)/2+1,sy4=stamplen,
.inorder=TRUE, .options.mpi=mpi.opts, .combine="rbind") %dopar% {
cbind(sum(sfam[sx1:sx2,sy1:sy2]*(im[x1:x2,y1:y2])),sum(sfam[sx3:sx4,sy1:sy2]*(im[x3:x4,y1:y2])),sum(sfam[sx1:sx2,sy3:sy4]*(im[x1:x2,y3:y4])),sum(sfam[sx3:sx4,sy3:sy4]*(im[x3:x4,y3:y4])))
}
}
if (verbose) { cat(" - Done\n") } # /*fend*/ }}}
#-----
#Integral of the reinterpolated psf; spsf /*fold*/ {{{
if (!no.psf) {
if (verbose) { cat(" Integral of the reinterpolated psf") }
spsf<-rep(NA,length(cat.id))
for (i in 1:length(psf)) {
spsf[which(psf.id==i)]<-sumpsf[i]
}
if (verbose) { cat(" - Done\n") }
} else {
spsf<-rep(NA, length(sfa))
}
# /*fend*/ }}}
#-----
#Integral of the estimated psf; spsf /*fold*/ {{{
if (any(estpsf.warn!=8)) {
if (verbose) { cat(" Integral of the estimated psf") }
sepsf<-rep(NA,length(cat.id))
for (i in 1:length(estpsf)) {
sepsf[which(epsf.id==i)]<-sumepsf[i]
}
if (verbose) { cat(" - Done\n") }
} else {
sepsf<-rep(NA, length(sfa))
}
# /*fend*/ }}}
#-----
#Integral of the reinterpolated psf^2; spsf2 /*fold*/ {{{
if (!no.psf) {
if (verbose) { cat(" Integral of the reinterpolated psf") }
spsf2<-rep(NA, length(sfa))
for (i in 1:length(psf)) {
spsf2[which(psf.id==i)]<-sum(psf[[i]]^2)
}
if (verbose) { cat(" - Done\n") }
} else {
spsf2<-rep(NA, length(sfa))
}
# /*fend*/ }}}
#-----
#Integral of the reinterpolated psf * estimated psf sepsfp /*fold*/ {{{
if (!no.psf & any(estpsf.warn!=8 & estpsf.warn!=-1)) {
if (verbose) { cat(" Integral of the reinterpolated psf * estimated psf") }
sepsfp<-rep(NA, length(sfa))
for (i in 1:length(psf)) {
for (j in 1:length(estpsf)) {
if (length(which(psf.id==i&epsf.id==j))>0) {
sepsfp[which(psf.id==i&epsf.id==j)]<-sum(psf.plot[[i]]*estpsf.plot[[j]])
}
}
}
if (verbose) { cat(" - Done\n") }
} else {
sepsfp<-rep(NA, length(sfa))
}
# /*fend*/ }}}
#-----
#Reinterpolate the PSF; psfi /*fold*/ {{{
if (!no.psf) {
if (!careful) {
if (verbose) { cat(" Reinterpolating the psf (per object)") }
psfi<-foreach(slen=stamplen, xc=cat.x, yc=cat.y, llim.x=ap.lims.data.map[,1],
ulim.x=ap.lims.data.map[,2],llim.y=ap.lims.data.map[,3],ulim.y=ap.lims.data.map[,4],
pid=psf.id, .options.mpi=mpi.opts, .noexport=ls(envir=environment()), .export="psf") %dopar% {
#limits of stamp in segmentation image space
x<-psf.cen[[pid]][1]
y<-psf.cen[[pid]][2]
lims<-zfloor(c(x,x,y,y)+c(-1,1,-1,1)*(c(xc-llim.x,ulim.x-xc,yc-llim.y,ulim.y-yc)))
lims[1]<-max(c(lims[1],1))
lims[3]<-max(c(lims[3],1))
lims[2]<-min(c(lims[2],length(psf[[pid]][,1])))
lims[4]<-min(c(lims[4],length(psf[[pid]][1,])))
psf.tmp<-psf[[pid]][lims[1]:lims[2],lims[3]:lims[4]]
#reinterpolate the PSF onto the new grid {{{
#Make grid for aperture at old centroid {{{
psf.obj<-list(x=(seq(lims[1],lims[2])-x),y=(seq(lims[3],lims[4])-y),z=psf.tmp)
#}}}
#Make expanded grid of new pixel centres {{{
expanded<-expand.grid(seq(llim.x,ulim.x)-xc,seq(llim.y,ulim.y)-yc)
#}}}
#Interpolate {{{
psf.sing=matrix(interp.2d(expanded[,1], expanded[,2], psf.obj)[,3], ncol=length(seq(llim.x,ulim.x)),nrow=length(seq(llim.y,ulim.y)))
#}}}
#}}}
psf.sing
}
if (verbose) { cat(" - Done\n") }
} else {
if (verbose) { cat(" Reinterpolating the psf (per object)") }
psfi<-foreach(slen=stamplen, xc=cat.x, yc=cat.y, pid=psf.id, .options.mpi=mpi.opts, .noexport=ls(envir=environment()), .export="psf") %dopar% {
#Make Sure PSF is centred on centre of stamp /*fold*/ {{{
centre<-psf.cen[[pid]]
delta<-floor(slen/2)*c(-1,+1)
lims<-rbind(centre[1]+delta,centre[2]+delta)
dx<-lims[1,1]-1
dy<-lims[2,1]-1
dx<-ifelse(dx>0,dx+1,dx)
dy<-ifelse(dy>0,dy+1,dy)
# /*fend*/ }}}
#Reinterpolate the PSF at point source XcenYcen /*fold*/ {{{
lenxpsf<-length(psf[[pid]][,1])
lenypsf<-length(psf[[pid]][1,])
lenx<-length(1:slen)
leny<-length(1:slen)
#Make grid for psf at old pixel centres /*fold*/ {{{
psf.obj<-list(x=seq(1,lenx), y=seq(1,leny),z=psf[[pid]][(1:(lenxpsf+1)+dx-1)%%(lenxpsf+1),(1:(lenypsf+1)+dy-1)%%(lenypsf+1)][1:lenx,1:leny])
# /*fend*/ }}}
#Make expanded grid of new pixel centres /*fold*/ {{{
expanded<-expand.grid(seq(1,lenx),seq(1,leny))
xnew<-expanded[,1]-xc%%1
ynew<-expanded[,2]-yc%%1
# /*fend*/ }}}
#Interpolate /*fold*/ {{{
ap<-matrix(interp.2d(xnew, ynew, psf.obj)[,3], ncol=leny,nrow=lenx)
# /*fend*/ }}}
# /*fend*/ }}}
ap
}
if (verbose) { cat(" - Done\n") }
}
} else {
psfi<-rep(NA, length(sfa))
}
# /*fend*/ }}}
#-----
#Integral of the (convolved aperture * psf); ssfap /*fold*/ {{{
if (!no.psf) {
if (verbose) { cat(" Integral of the (convolved aperture * psf)") }
ssfap<-foreach(sfam=sfa, ap=psfi, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% { sum(sfam*ap) }
if (verbose) { cat(" - Done\n") }
} else {
ssfap<-rep(NA, length(sfa))
}
# /*fend*/ }}}
#-----
#Integral of the (psf * image); spsfd /*fold*/ {{{
if (!no.psf) {
if (verbose) { cat(" Integral of the (psf * image)") }
if (cutup) {
spsfd<-foreach(ap=psfi, im=data.stamp, xlo=ap.lims.data.stamp[,1],xup=ap.lims.data.stamp[,2], ylo=ap.lims.data.stamp[,3],yup=ap.lims.data.stamp[,4], .combine='c', .options.mpi=mpi.opts, .inorder=TRUE, .noexport=ls(envir=environment())) %dopar% {
sum(ap*im[xlo:xup,ylo:yup])
}
} else {
spsfd<-foreach(ap=psfi, xlo=ap.lims.data.map[,1],xup=ap.lims.data.map[,2], ylo=ap.lims.data.map[,3],yup=ap.lims.data.map[,4], .combine='c', .options.mpi=mpi.opts, .inorder=TRUE, .noexport=ls(envir=environment()), .export=c("im")) %dopar% {
sum(ap*im[xlo:xup,ylo:yup])
}
}
if (verbose) { cat(" - Done\n") }
} else {
spsfd<-rep(NA, length(sfa))
}
# /*fend*/ }}}
#-----
#Integral of the (estimate psf * image); sepsfd /*fold*/ {{{
if (any(estpsf.warn!=8)) {
if (verbose) { cat(" Integral of the (estimated psf * image)") }
if (cutup) {
sepsfd<-foreach(slen=stamplen, xc=cat.x, yc=cat.y, epid=epsf.id, im=data.stamp, xlo=ap.lims.data.stamp[,1],xup=ap.lims.data.stamp[,2], ylo=ap.lims.data.stamp[,3],yup=ap.lims.data.stamp[,4], .combine='c', .options.mpi=mpi.opts, .inorder=TRUE, .noexport=ls(envir=environment()), .export="estpsf.plot") %dopar% {
#Make Sure PSF is centred on centre of stamp /*fold*/ {{{
centre<-estpsf.cen[[epid]]
delta<-floor(slen/2)*c(-1,+1)
lims<-rbind(centre[1]+delta,centre[2]+delta)
dx<-lims[1,1]-1
dy<-lims[2,1]-1
dx<-ifelse(dx>0,dx+1,dx)
dy<-ifelse(dy>0,dy+1,dy)
# /*fend*/ }}}
#Reinterpolate the PSF at point source XcenYcen /*fold*/ {{{
lenxpsf<-length(estpsf.plot[[epid]][,1])
lenypsf<-length(estpsf.plot[[epid]][1,])
lenx<-length(1:slen)
leny<-length(1:slen)
#Make grid for psf at old pixel centres /*fold*/ {{{
psf.obj<-list(x=seq(1,lenx), y=seq(1,leny),z=estpsf.plot[[epid]][(1:(lenxpsf+1)+dx-1)%%(lenxpsf+1),(1:(lenypsf+1)+dy-1)%%(lenypsf+1)][1:lenx,1:leny])
# /*fend*/ }}}
#Make expanded grid of new pixel centres /*fold*/ {{{
expanded<-expand.grid(seq(1,lenx),seq(1,leny))
xnew<-expanded[,1]-xc%%1
ynew<-expanded[,2]-yc%%1
# /*fend*/ }}}
#Interpolate /*fold*/ {{{
ap<-matrix(interp.2d(xnew, ynew, psf.obj)[,3], ncol=leny,nrow=lenx)
# /*fend*/ }}}
# /*fend*/ }}}
sum(ap*im[xlo:xup,ylo:yup])
}
} else {
sepsfd<-foreach(slen=stamplen, xc=cat.x, yc=cat.y, epid=epsf.id, xlo=ap.lims.data.map[,1],xup=ap.lims.data.map[,2], ylo=ap.lims.data.map[,3],yup=ap.lims.data.map[,4], .combine='c', .options.mpi=mpi.opts, .inorder=TRUE, .noexport=ls(envir=environment()), .export=c("estpsf.plot","im")) %dopar% {
#Make Sure PSF is centred on centre of stamp /*fold*/ {{{
centre<-estpsf.cen[[epid]]
delta<-floor(slen/2)*c(-1,+1)
lims<-rbind(centre[1]+delta,centre[2]+delta)
dx<-lims[1,1]-1
dy<-lims[2,1]-1
dx<-ifelse(dx>0,dx+1,dx)
dy<-ifelse(dy>0,dy+1,dy)
# /*fend*/ }}}
#Reinterpolate the PSF at point source XcenYcen /*fold*/ {{{
lenxpsf<-length(estpsf.plot[[epid]][,1])
lenypsf<-length(estpsf.plot[[epid]][1,])
lenx<-length(1:slen)
leny<-length(1:slen)
#Make grid for psf at old pixel centres /*fold*/ {{{
psf.obj<-list(x=seq(1,lenx), y=seq(1,leny),z=estpsf.plot[[epid]][(1:(lenxpsf+1)+dx-1)%%(lenxpsf+1),(1:(lenypsf+1)+dy-1)%%(lenypsf+1)][1:lenx,1:leny])
# /*fend*/ }}}
#Make expanded grid of new pixel centres /*fold*/ {{{
expanded<-expand.grid(seq(1,lenx),seq(1,leny))
xnew<-expanded[,1]-xc%%1
ynew<-expanded[,2]-yc%%1
# /*fend*/ }}}
#Interpolate /*fold*/ {{{
ap<-matrix(interp.2d(xnew, ynew, psf.obj)[,3], ncol=leny,nrow=lenx)
# /*fend*/ }}}
# /*fend*/ }}}
sum(ap*im[xlo:xup,ylo:yup])
}
}
if (verbose) { cat(" - Done\n") }
} else {
sepsfd<-rep(NA, length(sfa))
}
# /*fend*/ }}}
##-----
# #Integral of the deblended convolved aperture; sdfa /*fold*/ {{{
# if (verbose) { cat(" Integral of the deblended convolved aperture") }
# sdfa<-foreach(dfam=dfa, .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% { sum(dfam) }
# if (verbose) { cat(" - Done\n") } # /*fend*/ }}}
#-----
#Integral of the (deblended convolved aperture * image); sdfad /*fold*/ {{{
if (verbose) { cat(" Integral of the (deblended convolved aperture * image)") }
if (cutup) {
sdfad<-foreach(dfam=dfa,im=data.stamp, xlo=ap.lims.data.stamp[,1],xup=ap.lims.data.stamp[,2], ylo=ap.lims.data.stamp[,3],yup=ap.lims.data.stamp[,4], .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% {
sum(dfam*(im[xlo:xup,ylo:yup]))
}
} else {
sdfad<-foreach(dfam=dfa,.noexport=ls(envir=environment()), .export='im', xlo=ap.lims.data.map[,1],xup=ap.lims.data.map[,2], ylo=ap.lims.data.map[,3],yup=ap.lims.data.map[,4], .inorder=TRUE, .combine='c', .options.mpi=mpi.opts) %dopar% {
sum(dfam*(im[xlo:xup,ylo:yup]))
}
}
if (verbose) { cat(" - Done\n") } # /*fend*/ }}}
#-----
#Quartered Photometry - Quartered Integral of the (deblended convolved aperture * image); qsdfad /*fold*/ {{{
if (verbose) { cat(" Quartered Integral of the (deblended convolved aperture * image)") }
if (cutup) {
qsdfad<-foreach(dfam=dfa, im=data.stamp, x1=ap.lims.data.stamp[,1], x2=ap.lims.data.stamp[,1]+floor((ap.lims.data.stamp[,2]-ap.lims.data.stamp[,1])/2), x3=ap.lims.data.stamp[,1]+floor((ap.lims.data.stamp[,2]-ap.lims.data.stamp[,1])/2)+1, x4=ap.lims.data.stamp[,2],
y1=ap.lims.data.stamp[,3], y2=ap.lims.data.stamp[,3]+floor((ap.lims.data.stamp[,4]-ap.lims.data.stamp[,3])/2), y3=ap.lims.data.stamp[,3]+floor((ap.lims.data.stamp[,4]-ap.lims.data.stamp[,3])/2)+1, y4=ap.lims.data.stamp[,4],
sx1=rep(1,length(stamplen)),sx2=1+floor(stamplen-1)/2,sx3=1+floor(stamplen-1)/2+1,sx4=stamplen,
sy1=rep(1,length(stamplen)),sy2=1+floor(stamplen-1)/2,sy3=1+floor(stamplen-1)/2+1,sy4=stamplen,
.inorder=TRUE, .options.mpi=mpi.opts, .combine="rbind", .noexport=ls(envir=environment())) %dopar% {
cbind(sum(dfam[sx1:sx2,sy1:sy2]*(im[x1:x2,y1:y2])),sum(dfam[sx3:sx4,sy1:sy2]*(im[x3:x4,y1:y2])),sum(dfam[sx1:sx2,sy3:sy4]*(im[x1:x2,y3:y4])),sum(dfam[sx3:sx4,sy3:sy4]*(im[x3:x4,y3:y4])))
}
} else {
qsdfad<-foreach(dfam=dfa, .noexport=ls(envir=environment()), .export='im', x1=ap.lims.data.map[,1], x2=ap.lims.data.map[,1]+floor((ap.lims.data.map[,2]-ap.lims.data.map[,1])/2), x3=ap.lims.data.map[,1]+floor((ap.lims.data.map[,2]-ap.lims.data.map[,1])/2)+1, x4=ap.lims.data.map[,2],
y1=ap.lims.data.map[,3], y2=ap.lims.data.map[,3]+floor((ap.lims.data.map[,4]-ap.lims.data.map[,3])/2), y3=ap.lims.data.map[,3]+floor((ap.lims.data.map[,4]-ap.lims.data.map[,3])/2)+1, y4=ap.lims.data.map[,4],
sx1=rep(1,length(stamplen)),sx2=1+floor(stamplen-1)/2,sx3=1+floor(stamplen-1)/2+1,sx4=stamplen,
sy1=rep(1,length(stamplen)),sy2=1+floor(stamplen-1)/2,sy3=1+floor(stamplen-1)/2+1,sy4=stamplen,
.inorder=TRUE, .options.mpi=mpi.opts, .combine="rbind") %dopar% {
cbind(sum(dfam[sx1:sx2,sy1:sy2]*(im[x1:x2,y1:y2])),sum(dfam[sx3:sx4,sy1:sy2]*(im[x3:x4,y1:y2])),sum(dfam[sx1:sx2,sy3:sy4]*(im[x1:x2,y3:y4])),sum(dfam[sx3:sx4,sy3:sy4]*(im[x3:x4,y3:y4])))
}
}
if (verbose) { cat(" - Done\n") } # /*fend*/ }}}
#-----
#Integral of the [(convolved aperture * image error)^2]; ssfa2e2 /*fold*/ {{{
if (verbose) { cat(" Integral of the [(convolved aperture * image error)^2]") }
if (length(error.stamp)>1|(length(error.stamp)==1 & is.list(error.stamp))) {
ssfa2e2<-foreach(sfam=sfa, xlo=ap.lims.error.stamp[,1],xup=ap.lims.error.stamp[,2],ylo=ap.lims.error.stamp[,3],yup=ap.lims.error.stamp[,4], ime=error.stamp, .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% {
sum((sfam*ime[xlo:xup,ylo:yup])^2.)
}
} else if (length(error.stamp)==1){
if (error.stamp==1) {
ssfa2e2<-ssfa2
} else if (error.stamp==0) {
ssfa2e2<-rep(0,length(ssfa))
} else {
ssfa2e2<-foreach(sfam=sfa, .noexport=ls(envir=environment()), .export="error.stamp", .inorder=TRUE, .combine='c', .options.mpi=mpi.opts) %dopar% {
sum((sfam*error.stamp)^2.)
}
}
} else {
ssfa2e2<-foreach(sfam=sfa, xlo=ap.lims.error.map[,1],xup=ap.lims.error.map[,2],ylo=ap.lims.error.map[,3],yup=ap.lims.error.map[,4], .noexport=ls(envir=environment()), .export='ime', .inorder=TRUE, .combine='c', .options.mpi=mpi.opts) %dopar% {
sum((sfam*ime[xlo:xup,ylo:yup])^2.)
}
}
if (verbose) { cat(" - Done\n") } # /*fend*/ }}}
#-----
#Integral of the [(deblended convolved aperture * image error)^2]; sdfa2e2 /*fold*/ {{{
if (verbose) { cat(" Integral of the [(deblended convolved aperture * image error)^2]") }
if (length(error.stamp)>1|(length(error.stamp)==1 & is.list(error.stamp))) {
sdfa2e2<-foreach(dfam=dfa, xlo=ap.lims.error.stamp[,1],xup=ap.lims.error.stamp[,2],ylo=ap.lims.error.stamp[,3],yup=ap.lims.error.stamp[,4], ime=error.stamp, .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% {
sum((dfam*ime[xlo:xup,ylo:yup])^2.)
}
} else if (length(error.stamp)==1){
if (error.stamp==1) {
sdfa2e2<-sdfa2
} else if (error.stamp==0) {
sdfa2e2<-rep(0,length(ssfa))
} else {
sdfa2e2<-foreach(dfam=dfa, .noexport=ls(envir=environment()), .export="error.stamp", .inorder=TRUE, .combine='c', .options.mpi=mpi.opts) %dopar% {
sum((dfam*error.stamp)^2.)
}
}
} else {
sdfa2e2<-foreach(dfam=dfa, xlo=ap.lims.error.map[,1],xup=ap.lims.error.map[,2],ylo=ap.lims.error.map[,3],yup=ap.lims.error.map[,4], .noexport=ls(envir=environment()), .export='ime', .inorder=TRUE, .combine='c', .options.mpi=mpi.opts) %dopar% {
sum((dfam*ime[xlo:xup,ylo:yup])^2.)
}
}
if (verbose) { cat(" - Done\n") } # /*fend*/ }}}
#-----
#----- DIAGNOSTIC ITEMS -----# /*fold*/ {{{
if (diagnostic) {
#-----
#Integral of the (convolved aperture * (image error)^2); ssfae2 /*fold*/ {{{
if (verbose) { cat(" Integral of the (convolved aperture * image error)") }
if (length(error.stamp)>1|(length(error.stamp)==1 & is.list(error.stamp))) {
ssfae2<-foreach(sfam=sfa, xlo=ap.lims.error.stamp[,1],xup=ap.lims.error.stamp[,2],ylo=ap.lims.error.stamp[,3],yup=ap.lims.error.stamp[,4], ime=error.stamp, .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% {
sum(sfam*(ime[xlo:xup,ylo:yup])^2.)
}
} else if (length(error.stamp)==1){
if (error.stamp==1) {
ssfae2<-ssfa
} else if (error.stamp==0) {
ssfae2<-rep(0,length(ssfa))
} else {
ssfae2<-foreach(sfam=sfa, .noexport=ls(envir=environment()), .export="error.stamp", .inorder=TRUE, .combine='c', .options.mpi=mpi.opts) %dopar% {
sum(sfam*(error.stamp)^2.)
}
}
} else {
ssfae2<-foreach(sfam=sfa, xlo=ap.lims.error.map[,1],xup=ap.lims.error.map[,2],ylo=ap.lims.error.map[,3],yup=ap.lims.error.map[,4], .noexport=ls(envir=environment()), .export='ime', .inorder=TRUE, .combine='c', .options.mpi=mpi.opts) %dopar% {
sum(sfam*(ime[xlo:xup,ylo:yup])^2.)
}
}
if (verbose) { cat(" - Done\n") } # /*fend*/ }}}
#-----
#Integral of the [deblended convolved aperture * (image error)^2]; sdfae2 /*fold*/ {{{
if (verbose) { cat(" Integral of the [deblended convolved aperture * (image error)^2]") }
if (length(error.stamp)>1|(length(error.stamp)==1 & is.list(error.stamp))) {
sdfae2<-foreach(dfam=dfa, xlo=ap.lims.error.stamp[,1],xup=ap.lims.error.stamp[,2],ylo=ap.lims.error.stamp[,3],yup=ap.lims.error.stamp[,4], ime=error.stamp, .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% {
sum(dfam*(ime[xlo:xup,ylo:yup])^2.)
}
} else if (length(error.stamp)==1){
if (error.stamp==1) {
sdfae2<-sdfa2
} else if (error.stamp==0) {
sdfae2<-rep(0,length(ssfa))
} else {
sdfae2<-foreach(dfam=dfa, .noexport=ls(envir=environment()), .export="error.stamp", .inorder=TRUE, .combine='c', .options.mpi=mpi.opts) %dopar% {
sum(dfam*(error.stamp)^2.)
}
}
} else {
sdfae2<-foreach(dfam=dfa, xlo=ap.lims.error.map[,1],xup=ap.lims.error.map[,2],ylo=ap.lims.error.map[,3],yup=ap.lims.error.map[,4], .noexport=ls(envir=environment()), .export='ime', .inorder=TRUE, .combine='c', .options.mpi=mpi.opts) %dopar% {
sum(dfam*(ime[xlo:xup,ylo:yup])^2.)
}
}
if (verbose) { cat(" - Done\n") } # /*fend*/ }}}
#-----
#Integral of the (convolved aperture * image error); ssfae /*fold*/ {{{
if (verbose) { cat(" Integral of the (convolved aperture * image error)") }
if (length(error.stamp)>1|(length(error.stamp)==1 & is.list(error.stamp))) {
ssfae<-foreach(sfam=sfa, xlo=ap.lims.error.stamp[,1],xup=ap.lims.error.stamp[,2],ylo=ap.lims.error.stamp[,3],yup=ap.lims.error.stamp[,4], ime=error.stamp, .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% {
sum(sfam*ime[xlo:xup,ylo:yup])
}
} else if (length(error.stamp)==1){
if (error.stamp==1) {
ssfae<-ssfa
} else if (error.stamp==0) {
ssfae<-rep(0,length(ssfa))
} else {
ssfae<-foreach(sfam=sfa, .noexport=ls(envir=environment()), .export="error.stamp", .inorder=TRUE, .combine='c', .options.mpi=mpi.opts) %dopar% {
sum(sfam*error.stamp)
}
}
} else {
ssfae<-foreach(sfam=sfa, xlo=ap.lims.error.map[,1],xup=ap.lims.error.map[,2],ylo=ap.lims.error.map[,3],yup=ap.lims.error.map[,4], .noexport=ls(envir=environment()), .export='ime', .inorder=TRUE, .combine='c', .options.mpi=mpi.opts) %dopar% {
sum(sfam*ime[xlo:xup,ylo:yup])
}
}
if (verbose) { cat(" - Done\n") } # /*fend*/ }}}
#-----
#Integral of the (deblended convolved aperture * image error); sdfae /*fold*/ {{{
if (verbose) { cat(" Integral of the (convolved aperture * image error)") }
if (length(error.stamp)>1|(length(error.stamp)==1 & is.list(error.stamp))) {
sdfae<-foreach(dfam=dfa, xlo=ap.lims.error.stamp[,1],xup=ap.lims.error.stamp[,2],ylo=ap.lims.error.stamp[,3],yup=ap.lims.error.stamp[,4], ime=error.stamp, .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% {
sum(dfam*ime[xlo:xup,ylo:yup])
}
} else if (length(error.stamp)==1){
if (error.stamp==1) {
sdfae<-sdfa
} else if (error.stamp==0) {
sdfae<-rep(0,length(sdfa))
} else {
sdfae<-foreach(dfam=dfa, .noexport=ls(envir=environment()), .export="error.stamp", .inorder=TRUE, .combine='c', .options.mpi=mpi.opts) %dopar% {
sum(dfam*error.stamp)
}
}
} else {
sdfae<-foreach(dfam=dfa, xlo=ap.lims.error.map[,1],xup=ap.lims.error.map[,2],ylo=ap.lims.error.map[,3],yup=ap.lims.error.map[,4], .noexport=ls(envir=environment()), .export='ime', .inorder=TRUE, .combine='c', .options.mpi=mpi.opts) %dopar% {
sum(dfam*ime[xlo:xup,ylo:yup])
}
}
if (verbose) { cat(" - Done\n") } # /*fend*/ }}}
#-----
#Integral of the [(convolved aperture)^2]; ssfa2 /*fold*/ {{{
if (verbose) { cat(" Integral of the [(convolved aperture)^2]") }
ssfa2<-foreach(sfam=sfa, .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% { sum((sfam)^2.) }
if (verbose) { cat(" - Done\n") } # /*fend*/ }}}
#-----
#Integral of the [(deblended convolved aperture)^2]; sdfa2 /*fold*/ {{{
if (verbose) { cat(" Integral of the [(deblended convolved aperture)^2]") }
sdfa2<-foreach(dfam=dfa, .inorder=TRUE, .combine='c', .options.mpi=mpi.opts, .noexport=ls(envir=environment())) %dopar% { sum((dfam)^2.) }
if (verbose) { cat(" - Done\n") } # /*fend*/ }}}
#-----
}
#/*fend*/ }}}
#-----
#Final Flux & Error Calculations /*fold*/ {{{
if (!cutup) { detach(image.env) }
if (verbose) { cat(" Final Fluxes and Error Calculations ") }
#Calculate PSF-Weighting Correction /*fold*/ {{{
PSCorr<-spsf/spsf2
# /*fend*/ }}}
#Calculate Aperture-Weighting Correction /*fold*/ {{{
WtCorr<-ssfa/ssfau
# /*fend*/ }}}
#Calculate Minimum Aperture Correction /*fold*/ {{{
if (!no.psf) {
ApCorr<-spsf/ssfap
} else {
ApCorr<-1
}
# /*fend*/ }}}
#Calculate Estimated PSF Correction /*fold*/ {{{
if (any(estpsf.warn!=8)) {
EstPSFCorr<-sepsf/sepsfp
} else {
EstPSFCorr<-1
}
# /*fend*/ }}}
#Point Source flux = Int(PSF*Im) /*fold*/ {{{
psfflux<-spsfd
# /*fend*/ }}}
#Convolved aperture flux = Int(fltAp*Im) /*fold*/ {{{
sfaflux<-ssfad
# /*fend*/ }}}
#Deblended convolved aperture flux = Int(DBfltAp*Im) /*fold*/ {{{
dfaflux<-sdfad
# /*fend*/ }}}
#Point Source Flux error /*fold*/ {{{
if (!no.psf) {
psferr<-((conf*beamarea[psf.id])^2.*sqrt(spsf)*PSCorr)^2
} else {
psferr<-rep(NA,length(spsfd))
}
# /*fend*/ }}}
#Convolved aperture error /*fold*/ {{{
if (!no.psf) {
sfaerr<-ssfa2e2*ApCorr^2 + ((conf*beamarea[psf.id])^2.*sqrt(ssfa)*ApCorr)^2
} else {
sfaerr<-ssfa2e2*ApCorr^2 + ((conf*beamarea)^2.*sqrt(ssfa)*ApCorr)^2
}
if (blank.cor) { sfaerr<-sfaerr+(blanks$randMean.MAD*ApCorr)^2 }
# /*fend*/ }}}
#Deblended Convolved aperture error /*fold*/ {{{
if (!no.psf) {
dfaerr<-sdfa2e2*ApCorr^2 + ((conf*beamarea[psf.id])^2.*sqrt(sdfa)*ApCorr)^2
} else {
dfaerr<-sdfa2e2*ApCorr^2 + ((conf*beamarea)^2.*sqrt(sdfa)*ApCorr)^2
}
if (blank.cor) { dfaerr<-dfaerr+(blanks$randMean.MAD*ApCorr)^2 }
# /*fend*/ }}}
if (verbose) { cat(" - Done\n } ") }
# /*fend*/ }}}
#-----
#Notify /*fold*/ {{{
if (showtime) { cat(" - Done (",round(proc.time()[3]-timer[3],digits=2),"sec )\n")
message(paste("Perform Calculations - Done (",round(proc.time()[3]-timer[3],digits=2),"sec )"))
} else if (!quiet) { cat(" - Done\n") }
# /*fend*/ }}}
#Final Calculations /*fold*/ {{{
if (!quiet) { cat(" Performing Final Calculations ") }
if (do.sky.est|get.sky.rms) {
skyflux.mean<-skylocal.mean*sdfa
skyflux<-skylocal*sdfa
if (!blank.cor) {
sfaerr[which(!is.na(skyrms))]<-(sfaerr+(skyrms*sqrt(ssfa)*ApCorr)^2)[which(!is.na(skyrms))]
dfaerr[which(!is.na(skyrms))]<-(dfaerr+(skyrms*sqrt(sdfa)*ApCorr)^2)[which(!is.na(skyrms))]
}
psferr[which(!is.na(skyrms))]<-(psferr+(skyrms*sqrt(spsf)*PSCorr)^2)[which(!is.na(skyrms))]
if (do.sky.est) {
if (!quiet) { message("Perfoming Sky Subtraction"); cat("{\n Performing Sky Subtraction") }
#Subtract Sky Flux /*fold*/ {{{
dfaflux<-dfaflux-skyflux
sfaflux<-sfaflux-skylocal*ssfa
psfflux<-psfflux-skylocal*spsf
pixflux<-pixflux-skylocal
dfaerr[which(!is.na(skyerr))]<-(dfaerr+(skyerr*sdfa*ApCorr)^2)[which(!is.na(skyerr))]
sfaerr[which(!is.na(skyerr))]<-(sfaerr+(skyerr*ssfa*ApCorr)^2)[which(!is.na(skyerr))]
psferr[which(!is.na(skyerr))]<-(psferr+(skyerr*spsf*PSCorr)^2)[which(!is.na(skyerr))]
if (!quiet) { message(paste(" - Done\n")); cat(" - Done\n")}
# /*fend*/ }}}
}
}
#Deblend error /*fold*/ {{{
deblerr<-((1-sdfa/ssfa)*(1/sqrt(12))*abs(sfaflux)*ApCorr)
dfaerr<-dfaerr + (deblerr)^2
# /*fend*/ }}}
#Finalise Errors /*fold*/ {{{
psferr<-sqrt(psferr)
sfaerr<-sqrt(sfaerr)
dfaerr<-sqrt(dfaerr)
# /*fend*/ }}}
# /*fend*/ }}}
#Apply Aperture Correction to the Aperture Fluxess /*fold*/ {{{
psfflux<-psfflux*PSCorr
sfaflux<-sfaflux*ApCorr
dfaflux<-dfaflux*ApCorr
# /*fend*/ }}}
#Apply Additional Flux Correction to the Finalised Values /*fold*/ {{{
if (flux.corr!=1) {
psfflux<-psfflux*flux.corr
sfaflux<-sfaflux*flux.corr
dfaflux<-dfaflux*flux.corr
psferr<-psferr*flux.corr
sfaerr<-sfaerr*flux.corr
dfaerr<-dfaerr*flux.corr
}# /*fend*/ }}}
#Calculate magnitudes /*fold*/ {{{
if (magnitudes) {
suppressWarnings(mags<--2.5*(log10(dfaflux)-log10(ab.vega.flux))+mag.zp)
} else {
mags<-array(NA, dim=c(length(dfaflux)))
} # /*fend*/ }}}
#Get an estimate of the shot noise from the final flux /*fold*/ {{{
if (!is.null(image.env$gain) && !is.na(as.numeric(image.env$gain))) {
quick.shot<-sqrt(dfaflux)/as.numeric(image.env$gain)
} else {
quick.shot<-rep(NA,length(dfaflux))
}
# /*fend*/ }}}
#-----Diagnostic-----# /*fold*/ {{{
if (diagnostic) {
message(paste("After assignment",round(length(which(is.na(ssa )))/length(ssa )*100,digits=2),"% of the ssa matrix are NA"))
message(paste("After assignment",round(length(which(is.na(ssfa )))/length(ssfa )*100,digits=2),"% of the ssfa matrix are NA"))
message(paste("After assignment",round(length(which(is.na(sdfa )))/length(sdfa )*100,digits=2),"% of the sdfa matrix are NA"))
message(paste("After assignment",round(length(which(is.na(ssfa2 )))/length(ssfa2 )*100,digits=2),"% of the ssfa2 matrix are NA"))
message(paste("After assignment",round(length(which(is.na(sdfa2 )))/length(sdfa2 )*100,digits=2),"% of the sdfa2 matrix are NA"))
message(paste("After assignment",round(length(which(is.na(ssfad )))/length(ssfad )*100,digits=2),"% of the ssfad matrix are NA"))
message(paste("After assignment",round(length(which(is.na(sdfad )))/length(sdfad )*100,digits=2),"% of the sdfad matrix are NA"))
message(paste("After assignment",round(length(which(is.na(sfaflux)))/length(sfaflux)*100,digits=2),"% of the sfaflux matrix are NA"))
message(paste("After assignment",round(length(which(is.na(sfaerr )))/length(sfaerr )*100,digits=2),"% of the sfaerr matrix are NA"))
message(paste("After assignment",round(length(which(is.na(dfaflux)))/length(dfaflux)*100,digits=2),"% of the dfaflux matrix are NA"))
message(paste("After assignment",round(length(which(is.na(dfaerr )))/length(dfaerr )*100,digits=2),"% of the dfaerr matrix are NA"))
message(paste("After assignment",round(length(which(is.na(ssfae )))/length(ssfae )*100,digits=2),"% of the ssfae matrix are NA"))
message(paste("After assignment",round(length(which(is.na(ssfae2 )))/length(ssfae2 )*100,digits=2),"% of the ssfae2 matrix are NA"))
message(paste("After assignment",round(length(which(is.na(ssfa2e2)))/length(ssfa2e2)*100,digits=2),"% of the ssfa2e2 matrix are NA"))
message(paste("After assignment",round(length(which(is.na(sdfa2e2)))/length(sdfa2e2)*100,digits=2),"% of the sdfa2e2 matrix are NA"))
message(paste("After assignment",round(length(which(is.na(pixflux)))/length(pixflux)*100,digits=2),"% of the pixflux matrix are NA"))
}# /*fend*/ }}}
if (!quiet) { cat("\n} Galaxy Results Complete\n") }
# /*fend*/ }}}
# /*fend*/ }}}
#PART FIVE: OUTPUT /*fold*/ {{{
#Do we want to plot a sample of the apertures? /*fold*/ {{{
if (plot.sample) {
#Set output name /*fold*/ {{{
dir.create(file.path(path.root,path.work,path.out,"COGs"),showWarnings=FALSE)
# /*fend*/ }}}
#Determine Sample to Plot /*fold*/ {{{
if (!exists("plot.all")) { plot.all<-FALSE }
if (!exists("plot.sci")) { plot.sci<-FALSE }
if (plot.sci) {
#All Science targets /*fold*/ {{{
if (!quiet) { message("Writing All Science COGs to File"); cat("Writing All Science COGs to File") }
ind=which(contams==0)
# /*fend*/ }}}
} else if (plot.all) {
#All /*fold*/ {{{
if (!quiet) { message("Writing All COGs to File"); cat("Writing All COGs to File") }
ind=c(which(cat.a>0),which(cat.a==0))
# /*fend*/ }}}
} else {
if (!quiet) { message("Writing Sample of COGs to File"); cat("Writing Sample of COGs to File") }
#Output a random 15 apertures to file /*fold*/ {{{
ind1=which(cat.a>0 & sdfa!=0)
ind1=ind1[order(sdfad[ind1],decreasing=TRUE)][1:15]
ind2=which(sdfa!=0)[order(runif(1:length(which(sdfa!=0))))[1:15]]
ind<-c(ind1, ind2)
rm(ind1)
rm(ind2)
ind<-ind[which(!is.na(ind))]
# /*fend*/ }}}
}
# /*fend*/ }}}
for (i in ind) {
#Plot the CoGs for this source
get.stamp.cog()
}
#Remove unneeded Arrays /*fold*/ {{{
if (!make.resid.map) {
sfabak<-NULL
}
rm(dbw)
# /*fend*/ }}}
#Notify /*fold*/ {{{
if (!quiet) { message(paste(" - Done\n")); cat(" - Done\n")}
# /*fend*/ }}}
}
# /*fend*/ }}}
#If map was input in Jy/bm we need to convert it back before output in SourceSubtraction /*fold*/ {{{
if (Jybm & length(beamarea)==1) {
ba=beamarea
} else if (Jybm & length(beamarea)>1) {
ba<-sum(beamarea*psf.weight,na.rm=T)/sum(psf.weight)
} else {
ba=1.
}
# /*fend*/ }}}
#If wanted, make the Residual Map /*fold*/ {{{
if (make.resid.map) {
if (!is.null(sfabak)) { sfa<-sfabak }
if (filt.contam) {
if (!quiet) { cat(paste("Writing Contaminant-subtracted Map to",no.contam.map," ")) }
#Perform Source Subtraction /*fold*/ {{{
timer=system.time(source.subtraction(image.env$im,sfa,ap.lims.data.map,dfaflux,
file.path(path.root,path.work,path.out,no.contam.map),
image.env$data.hdr,ba,contams,diagnostic,verbose,
file.path(path.root,path.work,path.out,"Cotaminant_FluxMap.fits")))
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(paste('Contam Subtraction - Done (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
}
if (!quiet) { cat(paste("Writing Source-subtracted Map to",residual.map," ")) }
# /*fend*/ }}}
#Perform Source Subtraction /*fold*/ {{{
timer=system.time(source.subtraction(image.env$im,sfa,ap.lims.data.map,dfaflux,
file.path(path.root,path.work,path.out,residual.map),
image.env$data.hdr,ba,inside.mask,diagnostic,verbose,
file.path(path.root,path.work,path.out,"FluxMap.fits")))
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(paste('Source Subtraction - Done (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
# /*fend*/ }}}
}# /*fend*/ }}}
#If Subtracting Contaminants, then remove them before output /*fold*/ {{{
if ((write.tab)&(filt.contam)) {
cat.id <-cat.id[which(contams==0)]
cat.ra <-cat.ra[which(contams==0)]
cat.dec <-cat.dec[which(contams==0)]
cat.theta<-cat.theta[which(contams==0)]
theta.offset<-theta.offset[which(contams==0)]
x.pix <-x.pix[which(contams==0)]
y.pix <-y.pix[which(contams==0)]
cat.x <-cat.x[which(contams==0)]
cat.y <-cat.y[which(contams==0)]
cat.a <-cat.a[which(contams==0)]
cat.b <-cat.b[which(contams==0)]
ssa <-ssa[which(contams==0)]
ssfa <-ssfa[which(contams==0)]
ssfad <-ssfad[which(contams==0)]
qssfad <-rbind(qssfad[which(contams==0),])
spsf <-spsf[which(contams==0)]
spsf2 <-spsf2[which(contams==0)]
ssfap <-ssfap[which(contams==0)]
spsfd <-spsfd[which(contams==0)]
sfaflux<-sfaflux[which(contams==0)]
skyflux<-skyflux[which(contams==0)]
skylocal<-skylocal[which(contams==0)]
skyerr <-skyerr[which(contams==0)]
skyrms <-skyrms[which(contams==0)]
skypval<-skypval[which(contams==0)]
skyNBinNear <- skyNBinNear[which(contams==0)]
skyNBinNear.mean <- skyNBinNear.mean[which(contams==0)]
skyflux.mean<-skyflux.mean[which(contams==0)]
skyerr.mean <- skyerr.mean[which(contams==0)]
skylocal.mean <- skylocal.mean[which(contams==0)]
skypval.mean <- skypval.mean[which(contams==0)]
skyrms.mean <- skyrms.mean[which(contams==0)]
ssfa2e2<-ssfa2e2[which(contams==0)]
sdfa2e2<-sdfa2e2[which(contams==0)]
quick.shot<-quick.shot[which(contams==0)]
detecthres<-detecthres[which(contams==0)]
detecthres.mag<-detecthres.mag[which(contams==0)]
sfaerr <-sfaerr[which(contams==0)]
sdfa <-sdfa[which(contams==0)]
sdfad <-sdfad[which(contams==0)]
qsdfad <-rbind(qsdfad[which(contams==0),])
if (iterate.fluxes) {
fluxiters<-rbind(fluxiters[which(contams==0),])
erriters<-rbind(erriters[which(contams==0),])
sdfaiters<-rbind(sdfaiters[which(contams==0),])
iterateLost<-iterateLost[which(contams==0)]
}
psfflux<-psfflux[which(contams==0)]
psferr<-psferr[which(contams==0)]
dfaflux<-dfaflux[which(contams==0)]
deblerr <-deblerr[which(contams==0)]
dfaerr <-dfaerr[which(contams==0)]
saturated<-saturated[which(contams==0)]
pixflux<-pixflux[which(contams==0)]
PSCorr <-PSCorr[which(contams==0)]
ApCorr <-ApCorr[which(contams==0)]
WtCorr <-WtCorr[which(contams==0)]
stamplen<-stamplen[which(contams==0)]
mags <-mags[which(contams==0)]
if (length(flux.weight!=1)) { flux.weight<-flux.weight[which(contams==0)] }
if (ran.cor) { randoms<-rbind(randoms[which(contams==0),]) }
if (blank.cor) { blanks<-rbind(blanks[which(contams==0),]) }
if (diagnostic) {
ssfa2 <-ssfa2[which(contams==0)]
sdfa2 <-sdfa2[which(contams==0)]
ssfae <-ssfae[which(contams==0)]
sdfae <-sdfae[which(contams==0)]
ssfae2<-ssfae2[which(contams==0)]
sdfae2<-sdfae2[which(contams==0)]
}
if (exists("groups")) { groups<-groups[which(contams==0)] }
if (exists("psf.id")) { psf.id<-psf.id[which(contams==0)] }
if (exists("epsf.id")) { epsf.id<-epsf.id[which(contams==0)] }
contams <-contams[which(contams==0)]
}# /*fend*/ }}}
#Create Photometry Warning Flags /*fold*/ {{{
photWarnFlag<-rep("",length(cat.ra))
photWarnBitFlag<-rep(0,length(cat.ra))
#Bad Quartered Photometry /*fold*/ {{{
Qbad<-(apply(qsdfad,1,'max',na.rm=TRUE)/apply(qsdfad,1,'sum',na.rm=TRUE))>=0.7
Qbad[which(!is.finite(Qbad))]<-1
photWarnFlag<-paste0(photWarnFlag,ifelse(Qbad,"Q",""))
photWarnBitFlag<-photWarnBitFlag+ifelse(Qbad,2^0,0)
# /*fend*/ }}}
#Saturation /*fold*/ {{{
photWarnFlag<-paste0(photWarnFlag,ifelse(saturated,"X",""))
photWarnBitFlag<-photWarnBitFlag+ifelse(saturated,2^1,0)
# /*fend*/ }}}
#Bad Sky Estimate /*fold*/ {{{
if ((do.sky.est|get.sky.rms) & !quick.sky) {
temp<-skyNBinNear
temp[which(!is.finite(temp))]<-0
photWarnFlag<-paste0(photWarnFlag,ifelse(temp<=3,"S",""))
photWarnBitFlag<-photWarnBitFlag+ifelse(temp<=3,2^2,0)
}
# /*fend*/ }}}
#Iterative Deblend Warning /*fold*/ {{{
if (iterate.fluxes) {
photWarnFlag<-paste0(photWarnFlag,ifelse(iterateLost,"I",""))
photWarnBitFlag<-photWarnBitFlag+ifelse(iterateLost,2^3,0)
}
# /*fend*/ }}}
#PSF estimate Warning /*fold*/ {{{
if (psf.check) {
for (i in 1:length(estpsf.warnt)) {
photWarnFlag[which(epsf.id==i)]<-paste0(photWarnFlag[which(epsf.id==i)],estpsf.warnt[[i]]) #Can be all/none of [PM] or [E]; Peak/Median residual is bad, or estimate is bad
if(estpsf.warn[[i]]==8) { photWarnBitFlag[which(epsf.id==i)]<-photWarnBitFlag[which(epsf.id==i)]+2^4; estpsf.warn[[i]]<-estpsf.warn[[i]] - 8 }
if(estpsf.warn[[i]]>=4) { photWarnBitFlag[which(epsf.id==i)]<-photWarnBitFlag[which(epsf.id==i)]+2^5; estpsf.warn[[i]]<-estpsf.warn[[i]] - 4 }
if(estpsf.warn[[i]]>=2) { photWarnBitFlag[which(epsf.id==i)]<-photWarnBitFlag[which(epsf.id==i)]+2^6; estpsf.warn[[i]]<-estpsf.warn[[i]] - 2 }
if(estpsf.warn[[i]]>=1) { photWarnBitFlag[which(epsf.id==i)]<-photWarnBitFlag[which(epsf.id==i)]+2^7; estpsf.warn[[i]]<-estpsf.warn[[i]] - 1 }
}
}
# /*fend*/ }}}
# /*fend*/ }}}
#If wanted, output the Results Table /*fold*/ {{{
if (write.tab) {
#Output the results table /*fold*/ {{{
if ((loop.total!=1)&&(any(duplicated(file.path(param.env$path.root,param.env$path.work,param.env$path.out,param.env$tableout.name))))) {
if (!quiet) { cat(paste('Writing Results Table to ',file.path(path.root,path.work,path.out,paste(sep="",tableout.name,"_loop",f,".csv")),' ')) }
timer=system.time(write.flux.measurements.table(filename=file.path(path.root,path.work,path.out,paste(sep="",tableout.name,"_loop",f,".csv"))) )
} else {
if (!quiet) { cat(paste('Writing Results Table to ',file.path(path.root,path.work,path.out,paste(sep="",tableout.name,".csv")),' ')) }
timer=system.time(write.flux.measurements.table(filename=file.path(path.root,path.work,path.out,paste(sep="",tableout.name,".csv"))) )
}
if (showtime) { cat(" - Done (",round(timer[3],digits=2),"sec )\n")
message(paste('Write Results Table - Done (',round(timer[3], digits=2),'sec )'))
} else if (!quiet) { cat(" - Done\n") }
}# /*fend*/ }}}
#-----Diagnostic-----# /*fold*/ {{{
if (diagnostic) { message(paste('Sum of PSF = ',beamarea)) }
# /*fend*/ }}}
#If wanted, open the browser for user inspection /*fold*/ {{{
if (interact) {
cat(paste("Launching Interactive Mode: To end, type 'c'\n"))
sink(type='message')
browser()
sink(sink.file, type='message')
}# /*fend*/ }}}
# /*fend*/ }}}
# /*fend*/ }}}
#PART SIX: FINISH /*fold*/ {{{
#Send Parameters to logfile /*fold*/ {{{
sink(sink.file, type="output")
cat("Memory Hogs in this run:\n")
print(lsos(envir=environment(), head=TRUE, n=10))
cat("Images used in this run:\n")
print(lsos(envir=image.env, head=FALSE))
sink(type="output")
# /*fend*/ }}}
#Return /*fold*/ {{{
if (!quiet) { cat('\n') }
#on.exit(detach(image.env),add=TRUE)
if (!is.null(env)) {
on.exit(detach(env), add=TRUE)
}
if (!diagnostic) {
return=list(SFAflux=sfaflux,SFAerror=sfaerr, DFAflux=dfaflux,DFAerror=dfaerr,MinApCorr=ApCorr,MaxApCorr=WtCorr)
} else {
return=list(SFAflux=sfaflux,SFAerror=sfaerr, DFAflux=dfaflux,DFAerror=dfaerr, SA_Stamps=sa, SFA_Stamps=sfa, WSFA_Stamps=wsfa, DFA_Stamps=dfa,MinApCorr=ApCorr,MaxApCorr=WtCorr)
}
# /*fend*/ }}}
# /*fend*/ }}}
}
#Func /*fold*/ {{{
.executeran.cor<-function() {
cat("Executing ran.cor...\n"); Sys.sleep(2); cat('\n | _______ _ _ _ _ _ __ _ |\n | |__ __| | | | | | | | | | __ \\ | | |\n | | | | |___| | | | | | | | | | | | |\n | | | | ___ | | | | | | | | | |_| |\n | | | | | | | | |___| | | |__| | _ |\n | |_| |_| |_| |_______| |_____/ |_| |\n | |\n \\ /\\ /\\ /\\ /\\ /\\ /\\ /\\ /\\ /\\ /\\ /\\ /\n \\/ \\/ \\/ \\/ ,------------, \\/ \\/ \\/ \\/\n ____ ({ XX XX }) ____\n /////|\\ _| \\\\ // |_ /|\\\\\\\\\\\n VVVV | \\____/ { \\ \\\\// / } \\____/ | VVVV\n ////_| ,|V=V=V=V=V=V=|, |_\\\\\\\\\n ___\\\\\\\\/_\\~~~~~/_{+^+^+^+^+^+^}_\\~~~~~/_\\////___\n\n'); Sys.sleep(2); cat("... Done, you heartless Jedi Scum\n\n"); Sys.sleep(2)
}
# /*fend*/ }}}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.