PRO GALPROF, image, mask, result, EXTEN=exten, XCEN=cenx, YCEN=ceny,$ EPS=epsilon, EEPS=eepsilon, POSANG=posang, EPOSANG=eposang,$ RMIN=rmin, GROW=grow, PMZERO=pmzero, EPMZERO=epmzero,$ PAEXT=pairmass, EPAEXT=epairmass,$ PCOLOR=pcolor, EPCOLOR=epcolor, COLOR=color,$ FIXSKY=fixsky, EFIXSKY=efixsky, DARKRATE=darkrate,$ EXPTIME=exptime, SECPIX=secpix, ESECPIX=esecpix,$ AIRMASS=airmass, EAIRMASS=eairmass, GAIN=gain, RDNOISE=rdnoise,$ TEXTLOG=textlog, DEBUG=debug, NOPLOT=noplot, NOREG=noreg ;+---------------------------------------------------------------------------- ; PROGRAM: GALPROF Author: Rolf A. Jansen (ASU) -- Jun 3 2015 ; ; PURPOSE: ; Perform simple galaxy surface photometry, assuming elliptical isophotes, a ; fixed center position, and a fixed ellipticity and position angle of those ; elliptical isophotes. A mask image is used to identify background regions ; and regions to omit from the analysis altogether. To be used when IDL is ; all that's available. ; ; METHOD: ; This routine will read input FITS image 'image' and input mask image 'mask' ; into memory, where that mask should have a value of 0 for all background ; pixels, >0 (e.g., 1) for pixels that might have galaxy signal, and <0 (-1) ; for any pixels that should be discarded in the analysis. Unless overridden ; by explicitly given values of 'fixsky' and 'efixsky', it will then estimate ; the global sky background level 'sky' and its uncertainty 'skyrms' for all ; valid background pixels (as defined by the mask image). Unless a galaxy ; center position ('xcen','ycen') is provided, it will estimate one from the ; 1st order moments of all valid galaxy pixels (as defined by the mask image). ; Unless either or both the ellipticity 'eps' and position angle 'posang' ; (measured counter-clockwise from the positive x-axis) is specified, it will ; estimate the shape and orientation from the 2nd order moments of all valid ; galaxy pixels (as defined by the mask image). Given the center, shape and ; orientation of the galaxy, an intermediate 2-D array of semi-major axis ; distances is created (which can be saved by setting the /DEBUG switch). ; Next, aperture photometry is performed within elliptical annuli at a fixed ; set of radii up to 'rmin', and at radii that grow geometrically by a factor ; 'grow' from 'rmin' till the edge of the image (or of the region with un- ; masked pixels defined in the mask image) is encountered (fewer than 50% of ; the expected number of pixels in that annulus). Photometric calibration, if ; provided ('pmzero', 'paext', 'pcolor' and 'color'), will be applied to ; yield a radial surface brightness profile. Some additional information may ; be needed if not provided in the original image header (exposure time, ; pixel scale, airmass at middle of the exposure or an effective airmass of ; an image stack, as well as readnoise and gain for error propagation). ; Output 'result' will be a SuperMongo-readable ASCII table with a human- ; readable header of space-separated columns containing effective radius, ; inner and outer radii of the annulus at that radius, mean intensity and ; error thereon, mean surface brightness and upper and lower errors thereon. ; Optionally, plotting results to IDL graphics windows can be disabled with ; the '/NOPLOT' switch, writing a DS9 regions file with elliptical annuli can ; be disabled with the '/NOREG' switch, and the verbosity level in the text ; log file can be turned up by providing the '/DEBUG' switch. ; ; CALLING SEQUENCE: ; GALPROF, image, mask, result, [EXTEN=#], [XCEN=###], [YCEN=###], ; [EPS=###], [EEPS=###], [POSANG=###], [EPOSANG=###], [RMIN=###], ; [GROW=###], [PMZERO=###], [EPMZERO=###], [PAEXT=###], ; [EPAEXT=###], [PCOLOR=###], [EPCOLOR=###], [COLOR=###], ; [FIXSKY=###], [EFIXSKY=###], [DARKRATE=###], ; [EXPTIME=###], [SECPIX=###], [ESECPIX=###], ; [AIRMASS=###], [EAIRMASS=###], [GAIN=###], [RDNOISE=###], ; [TEXTLOG=%%%], [/DEBUG], [/NOPLOT], [/NOREG] ; ; OUTPUT: ; result (string) ("") - Name of output ASCII table ; textlog (string)("galprof.log") - Name of text log file. ; ; INPUT: ; image (string) ("") - Name of input image (flat FITS or extension in a ; Multi-Extension FITS file). ; mask (string) ("") - Name of input mask image (flat FITS) identifying ; galaxy and background pixels, and pixels to omit. ; OPTIONAL INPUT: ; exten (int) (1) - Desired frame within input Multi-Extension FITS file ; xcen (real) (auto) - Object center position along pixel columns (x). ; ycen (real) (auto) - Object center position along pixel rows (y). ; eps (real) (auto) - Object ellipticity (1-b/a). ; eeps (real) (auto) - Uncertainty thereon. ; posang (real) (auto) - Object orientation in degrees wrt. to pos. x-axis. ; eposang (real) (auto) - Uncertainty thereon. ; rmin (real) (4.0) - Equivalent radius in pixels of smallest aperture to ; grow (apertures smaller than 'rmin' are drawn from ; the fixed sequence: 1.0, 1.6, 2.2, 2.55, ; sqrt(8), 3., sqrt(10), sqrt(13), 4., sqrt(18), ; sqrt(20), 5., sqrt(29), sqrt(32), sqrt(34), 6.). ; Make sure 'rmin' is <= 6 pixels. ; grow (real) (1.1) - Factor by which to grow subsequent aperture radii. ; Make sure 1.0 <= 'grow' <= 1.5 . ; pmzero (real) (0.00) - Photometric zero point (mag). ; epmzero (real) (0.00) - Error thereon (mag). ; paext (real) (0.00) - Airmass term in photometric calibration (mag/A.M.). ; epaext (real) (0.00) - Error thereon (mag/A.M.). ; pcolor (real) (0.00) - Color term in photometric calibration (mag/mag). ; epcolor (real) (0.00) - Error thereon (mag/mag). ; color (real) (0.00) - Calibrated color with which to multiply 'pcolor' ; Except in rare cases, it is recommended to omit any ; color term correction until after a profile in an ; appropriate second filter has been constructed using ; 'galprof'; then iterate until the color converges. ; fixsky (real) (auto) - Skip automatic background estimation and use this ; value for the sky background instead (ADU). ; efixsky (real) (auto) - Uncertainty therein (ADU). ; darkrate(real) (0.00) - Detector dark rate [e-/pixel/hour]. ; ; The values of the following parameters will be attempted to be read from the ; image header, but may be overridden when explicitly given. They are set to the ; default values when neither specified, nor successfully read from the header. ; exptime (real) (1.00) - Exposure time [s] (assumed to be equal to dark time) ; secpix (real) (1.00) - Pixel size (arcsec) ; esecpix (real) (0.00) - Uncertainty thereon (usually not in headers). ; airmass (real) (1.00) - Airmass at middle of exposure, or effective airmass ; for the image stack. ; eairmass(real) (0.00) - Uncertainty thereon (usually not in headers). ; gain (real) (1.00) - Gain factor [e-/ADU] ; rdnoise (real) (0.00) - Read noise [e-] ; ; Last, there are a couple of switches that modify the behavior of 'galprof': ; /DEBUG - Turn up the verbosity level in the text log file. ; /NOREG - Do not write DS9 regions file of elliptical annuli ; /NOPLOT - Don't open IDL direct graphics windows to show the ; results of the octant sky boxes or surface bright- ; ness profile measured. ; ; EXAMPLES: ; 1. Uncalibrated profile (raw instrument magnitudes): ; IDL> galprof, "N5888B.fits", "N5888B_msk.fits", "N5888B.prof" ; ; 2. Calibrated surface brightness profile (excl. color term): ; IDL> galprof, "N5888B.fits", "N5888B_msk.fits", "N5888B.prof",$ ; pmzero=21.6362, epmzero=0.0164, paext=0.3017, epaext=0.0119,$ ; pcolor=-0.06088, epcolor=0.00374, color=0.00, darkrate=1.0, /NOPLOT ; ; 3. Calibrated surface brightness profile (excl. color term), with all ; parameters explicitly specified: ; IDL> galprof, "N5888B.fits", "N5888B_msk.fits", "N5888B.prof",$ ; xcen=257.55, ycen=255.65, eps=0.400, eeps=0.005,$ ; posang=68., eposang=1., rmin=4.0, grow=1.1,$ ; pmzero=21.6362, epmzero=0.0164, paext=0.3017, epaext=0.0119,$ ; pcolor=-0.06088, epcolor=0.00374, color=0.00,$ ; fixsky=149.7826, efixsky=3.2372, darkrate=1.0,$ ; secpix=0.6200, esecpix=0.0025, gain=2.45, rdnoise=9.5,$ ; exptime=900., airmass=1.069, eairmass=0.001, /DEBUG ; ; USES THE FOLLOWING 'ASTROLIB' ROUTINES: ; readfits.pro, writefits.pro, sxpar.pro, sxaddpar.pro, mmm.pro, ; resistant_mean.pro, get_date.pro ; ; MODIFICATION HISTORY: ; Written by R.A. Jansen, Arizona State University, May 27--Jun 2 2015 (for ; IDL 8.3) ; RAJ, Jun 3 2015, Implemented automatic centroiding and measurement of ; object ellipticity and position angle. Changed some task ; parameter names and their order. Improved retrieval of ; info from image header. Updated manual. ; RAJ, Jun 4 2015, Added keyword /NOREG to disable creating DS9 regions file ; of the elliptical apertures (annuli) and effective radii; ; fixed logical flaw that image moments were only computed ; when xcen and ycen were not specified; final tests and ; minor bug fixes. ; RAJ, Jun 17 2015, Implemented fit of residual background gradients based on ; octant analysis. ; TODO: split the moments analysis and the octant sky analysis off to two ; separate routines. ;- --------------------------------------------------------------------------- progver = 'GALPROF v1.0[IDL8.3] - R.A. Jansen - (IDL'+STRING(!Version.release,FORMAT='(f4.1)')+')' ; DEFINE SOME PROGRAM DEFAULT PARAMETERS... x0 = 0.0 y0 = 0.0 ba = 1.0 eba = 0.025 theta = 0.0 etheta = 1.0 rminr = 4.0 dr = 1.1 xtnd = 'F' nxtn = 1 pz = 0.00 epz = 0.00 pam = 0.00 epam = 0.00 pc = 0.00 epc = 0.00 clr = 0.00 epscl = 0.0025 eairm = 0.001 dc = 0.0 textlog="galprof.log" idebug = 0 inoplot = 1 inoreg = 1 ; DEGREES-TO-RADIANS CONVERSION FACTOR... dtr = !DPI/180.0D ; ------------------ I N P U T V E R I F I C A T I O N ------------------- IF (n_params(0) LT 3 ) THEN BEGIN print,'SYNTAX: GALPROF, image, mask, result, [EXTEN=#], [XCEN=###], [YCEN=###],' print,' [EPS=###],[EEPS=###], [POSANG=###],[EPOSANG=###], [RMIN=###],' print,' [GROW=###], [PMZERO=###],[EPMZERO=###], [PAEXT=###],[EPAEXT=###],' print,' [PCOLOR=###],[EPCOLOR=###],[COLOR=###], [FIXSKY=###],[EFIXSKY=###],' print,' [DARKRATE=###], [EXPTIME=###], [SECPIX=###],[ESECPIX=###],' print,' [AIRMASS=###], [EAIRMASS=###], [GAIN=###], [RDNOISE=###],' print,' [TEXTLOG=%%%%], [/DEBUG], [/NOPLOT], [/NOREG]' print,' where "image" and "mask" are names of the input and mask FITS files, and' print,' "result" the name of an ASCII output table. If the input images are' print,' Multi-Extension FITS files, then "exten" must specify the relevant exten-' print,' sion (must be the same in "image" and "mask").' print,' If the object center position, shape, and orientation are known, one can' print,' specify "xcen" and "ycen", ellipticity "eps" (where eps=1-b/a), and' print,' position angle "posang" (measured counter-clockwise from the positive' print,' x-axis. Optionally, one can also specify the uncertainties on the shape' print,' and orientation as "eeps" and "eposang". If these parameters are absent,' print,' the task will estimate them from the image moment up to 2nd order.' print,' Up to minimum radius "rmin" (<=6 pix) a fixed sequence of elliptical' print,' apertures are used, beyond which subsequent elliptical aperture radii' print,' are geometrically grown by a factor "grow" to maintain a roughtly constant' print,' S/N per annulus in the case of an exponential galaxy disk.' print,' The radial intensity profile may be photometrically calibrated if the' print,' calibration constants "pmzero", "pext", "pcolor" and "color", and the' print,' uncertainties thereon ("epmzero", "epext", "epcolor" and "ecolor") are' print,' given, and if the exposure time, pixel size, and airmass can be deter-' print,' mined from the image header, or if they are explicitly given through' print,' "exptime", "secpix", and "airmass" (and optional errors "esecpix" and ' print,' "eairmass"). For proper error propagation, the task also attempts to read' print,' the readnoise and gain from the image header. Automatic fitting of the' print,' sky background level is disable when "fixsky" is specified (as well as' print,' its uncertainty "epfixsky"). Running this routine with its default values' print,' will yield uncalibrated surface brightness profiles!' RETURN ENDIF ; OVERRIDE PROGRAM DEFAULTS WITH VALUES OF SPECIFIC KEYWORDS... IF (keyword_set(rmin)) THEN rminr = rmin IF (rminr LT 1.0 OR rminr GT 10.0) THEN BEGIN print,'ERROR: value of rmin outside valid range [1.0--10.0] pix!' RETURN ENDIF IF (keyword_set(grow)) THEN dr = grow IF (dr LT 1.0 OR dr GT 1.5) THEN BEGIN print,'ERROR: value of grow outside valid range [1.0--1.5] !' RETURN ENDIF IF (keyword_set(exten)) THEN nxtn = exten IF (nxtn LT 0) THEN BEGIN print,'ERROR: value of exten < 0 !' RETURN ENDIF IF (keyword_set(pmzero)) THEN pz = pmzero IF (keyword_set(epmzero)) THEN epz = epmzero IF (keyword_set(pairmass)) THEN pam = pairmass IF (keyword_set(epairmass)) THEN epam = epairmass IF (keyword_set(pcolor)) THEN pc = pcolor IF (keyword_set(epcolor)) THEN epc = epcolor IF (keyword_set(color)) THEN clr = color IF (keyword_set(esecpix)) THEN epscl = esecpix IF (keyword_set(eairmass)) THEN eairm = eairmass IF (keyword_set(textlog)) THEN textlog= textlog IF (keyword_set(darkrate)) THEN dc = darkrate IF (keyword_set(DEBUG)) THEN idebug = 1 IF (keyword_set(NOPLOT)) THEN inoplot= 0 IF (keyword_set(NOREG)) THEN inoreg = 0 IF ( epz GT ABS(pz) OR epam GT ABS(pam) OR epc GT ABS(pc) ) THEN BEGIN IF ( epz GT ABS(pz) ) THEN print,'WARNING: epmzero > pmzero !' IF ( epam GT ABS(pam) ) THEN print,'WARNING: epairmass > pairmass !' IF ( epc GT ABS(pc) ) THEN print,'WARNING: epcolor > pcolor !' RETURN ENDIF ; --------------------- I N P U T I M A G E I / O ---------------------- ; OPEN TEXT LOGFILE... OPENW,25,textlog,/APPEND ; PRINT OPENING MESSAGE... print,'GALPROF: image="'+image+'", mask="'+mask+'"' print,' result="'+result+'", logfile="'+textlog+'"' printf,25,' ' printf,25,'# '+progver+': '+SYSTIME(/UTC) printf,25,' image ="'+image+'"' printf,25,' mask ="'+mask+'"' printf,25,' result="'+result+'"' ; READ THE FITS PRIMARY HEADER UNIT, AND CHECK WHETHER "EXTEND == T"... phu = READFITS(image,primehead,exten=0,/silent) xtnd = STRTRIM(SXPAR(primehead,'EXTEND')) ; READ THE (SPECIFIED EXTENSION IN THE) FITS FILE... IF (xtnd EQ 'T') THEN BEGIN img = READFITS(image,extenhead,exten=nxtn,/silent) hdr = extenhead msk = READFITS(mask,maskhead,exten=nxtn,/silent) ENDIF ELSE BEGIN img = READFITS(image,primehead,/silent) hdr = primehead msk = READFITS(mask,maskhead,/silent) ENDELSE ; DETERMINE THE IMAGE SIZE (number of columns and lines)... image_sz = SIZE(img) ncol = image_sz[1] nrow = image_sz[2] ; array limits in IDL coordinate convention... nr = nrow - 1 nc = ncol - 1 ; READ SOME MORE RELEVANT INFO FROM THE IMAGE HEADER... IF (keyword_set(rdnoise)) THEN BEGIN ron = rdnoise cf1 = 'user' ENDIF ELSE BEGIN nk = -1 cf1 = 'header' ron = FLOAT(STRTRIM(SXPAR(hdr,'RDNOISE',COUNT=nk))) IF (nk EQ 0 OR !ERR EQ -1) THEN BEGIN !ERR = 0 ron = 0.00 & cf1 = 'default' ENDIF ENDELSE IF (keyword_set(gain)) THEN BEGIN eadu = gain cf2 = 'user' ENDIF ELSE BEGIN nk = -1 cf2 = 'header' eadu = FLOAT(STRTRIM(SXPAR(hdr,'GAIN',COUNT=nk))) IF (nk EQ 0 OR !ERR EQ -1) THEN BEGIN !ERR = 0 eadu = 1.00 & cf2 = 'default' ENDIF ENDELSE IF (keyword_set(exptime)) THEN BEGIN texp = exptime cf3 = 'user' ENDIF ELSE BEGIN nk = -1 cf3 = 'header' texp = FLOAT(STRTRIM(SXPAR(hdr,'EXPTIME',COUNT=nk))) IF (nk EQ 0 OR !ERR EQ -1) THEN BEGIN !ERR = 0 texp = 1.00 & cf3 = 'default' ENDIF ENDELSE IF (keyword_set(secpix)) THEN BEGIN pscl = secpix cf4 = 'user' ENDIF ELSE BEGIN nk = -1 cf4 = 'header' pscl = FLOAT(STRTRIM(SXPAR(hdr,'SECPIX',COUNT=nk))) IF (nk EQ 0 OR !ERR EQ -1) THEN BEGIN !ERR = 0 pscl = 1.00 & cf4 = 'default' ENDIF ENDELSE IF (keyword_set(airmass)) THEN BEGIN print,'Hoi Hoi' airm = airmass cf5 = 'user' ENDIF ELSE BEGIN print,'Boo Boo!' nk = -1 cf5 = 'header' airm = FLOAT(STRTRIM(SXPAR(hdr,'AIRMASS',COUNT=nk))) IF (nk EQ 0 OR !ERR EQ -1) THEN BEGIN !ERR = 0 airm = 1.00 & cf5 = 'default' ENDIF ENDELSE ; ADOPT OR DETERMINE THE (APPROXIMATE) GLOBAL SKY BACKGROUND LEVEL... ig = WHERE(msk EQ 0, ign) IF (keyword_set(fixsky)) THEN BEGIN ; adopt user-specified value of sky background and error thereon... sky = fixsky IF (keyword_set(efixsky)) THEN BEGIN esky = efixsky ENDIF ELSE BEGIN esky = 0.0 ENDELSE skyrms = esky print,' adopting user-specified sky background level...' print,' sky=',sky,' +/-',esky,' [fixed]',f='(a,f9.3,a,f7.3,a)' printf,25,' adopting user-specified sky background level...' printf,25,' sky=',sky,' +/-',esky,' [fixed]',f='(a,f9.3,a,f7.3,a)' ; wild guesstimate of number of sky pixels needed for error propagation... nsky = ign ENDIF ELSE BEGIN ; determine the approximate global sky background level from all pixels ; that the mask image identifies as valid background pixels (msk=0) ... MMM, img(ig), sky, esky, skew it = WHERE(msk EQ 0 AND (img GT (sky-4.*esky) AND img LT (sky+3.*esky)), itn) tmpimg = img(it) print,' estimating sky background level [',ign,$ ' unmasked background pixels]...',f='(a,i8,a)' print,' iter=0: sky=',sky,' +/-',esky,$ ' adu (nrej=',(ign-itn),')',f='(a,f9.3,a,f8.3,a,i8,a)' printf,25,' estimating sky background level (',ign,$ ' unmasked background pixels)...',f='(a,i8,a)' printf,25,' iter=0: sky=',sky,' +/-',esky,$ ' adu (nrej=',(ign-itn),')',f='(a,f9.3,a,f8.3,a,i8,a)' RESISTANT_Mean, tmpimg, 4.0, sky, esky, nrej print,' iter=1: sky=',sky,' +/-',esky,$ ' adu (nrej=',nrej,')',f='(a,f9.3,a,f8.3,a,i8,a)' printf,25,' iter=1: sky=',sky,' +/-',esky,$ ' adu (nrej=',nrej,')',f='(a,f9.3,a,f8.3,a,i8,a)' nsky = itn - nrej ENDELSE ; SUBTRACT ESTIMATED OR SPECIFIED SKY BACKGROUND... img = TEMPORARY(img) - sky ; PERFORM ANALYSIS OF IMAGE MOMENTS UP TO 2nd ORDER WHEN ANY ONE OR MORE OF ; THE CENTER POSITION, SHAPE, OR ORIENTATION PARAMETERS ARE MISSING... IF ( (NOT keyword_set(cenx)) OR (NOT keyword_set(ceny)) OR (NOT keyword_set(epsilon)) OR (NOT keyword_set(posang)) ) THEN BEGIN ; NOTE: this is a re-implementation of my heritage f77 code 'moments.f'. ig = WHERE(msk GT 0, ign) M0 = 0.0 & eM0 = 0.0 Mx = 0.0 & eMx = 0.0 My = 0.0 & eMy = 0.0 FOR row=0,nr DO BEGIN FOR col=0,nc DO BEGIN IF (msk(col,row) GT 0) THEN BEGIN val = img(col,row) eval = SQRT( ron*ron + eadu*ABS(sky+val) )/eadu M0 = M0 + val eM0 = eM0 + eval*eval Mx = Mx + FLOAT(col)*val eMx = eMx + (FLOAT(col)*eval)^2 My = My + FLOAT(row)*val eMy = eMy + (FLOAT(row)*eval)^2 ENDIF ENDFOR ENDFOR eM0 = SQRT(eM0) Mx = Mx/M0 My = My/M0 eMx = SQRT( eMx/(M0*M0) + ((Mx*eM0)/(M0*M0))^2 ) eMy = SQRT( eMy/(M0*M0) + ((My*eM0)/(M0*M0))^2 ) ; estimated center position and error thereon... _x0 = Mx & _ex0 = eMx _y0 = My & _ey0 = eMy ; Second pass on only the central 15x15 pix region... MM0 = 0.0 & eMM0 = 0.0 MMx = 0.0 & eMMx = 0.0 MMy = 0.0 & eMMy = 0.0 FOR row=ROUND(_y0-7.5),ROUND(_y0+7.5) DO BEGIN FOR col=ROUND(_x0-7.5),ROUND(_x0+7.5) DO BEGIN IF (msk(col,row) GT 0) THEN BEGIN val = img(col,row) eval = SQRT( ron*ron + eadu*ABS(sky+val) )/eadu MM0 = MM0 + val eMM0 = eMM0 + eval*eval MMx = MMx + FLOAT(col)*val eMMx = eMMx + (FLOAT(col)*eval)^2 MMy = MMy + FLOAT(row)*val eMMy = eMMy + (FLOAT(row)*eval)^2 ENDIF ENDFOR ENDFOR eMM0 = SQRT(eMM0) MMx = MMx/MM0 MMy = MMy/MM0 eMMx = SQRT( eMMx/(MM0*MM0) + ((MMx*eMM0)/(MM0*MM0))^2 ) eMMy = SQRT( eMMy/(MM0*MM0) + ((MMy*eMM0)/(MM0*MM0))^2 ) ; refined estimate of center position and error thereon... x0 = MMx & ex0 = eMMx y0 = MMy & ey0 = eMMy IF ( (NOT keyword_set(epsilon)) OR (NOT keyword_set(posang)) ) THEN BEGIN ; either or both ellipticity and position angle were not given as input ; parameters either: compute the 2nd order moments of the image pixels ; with respect to centroid (x0,y0) from which ellipticity and position ; angle can be derived, as needed. Mxx = 0.0 & eMxx = 0.0 Myy = 0.0 & eMyy = 0.0 Mxy = 0.0 & eMxy = 0.0 FOR row=0,nr DO BEGIN FOR col=0,nc DO BEGIN IF (msk(col,row) GT 0) THEN BEGIN val = img(col,row) eval = SQRT( ron*ron + eadu*ABS(sky+val) )/eadu xx = FLOAT(col-x0) exx = eMx yy = FLOAT(row-y0) eyy = eMy Mxx = Mxx + xx*xx*val eMxx = eMxx + (2.0*xx*exx*val)^2 + (xx*xx*eval)^2 Mxy = Mxy + xx*yy*val eMxy = eMxy + (exx*yy*val)^2 + (xx*eyy*val)^2 + (xx*yy*eval)^2 Myy = Myy + yy*yy*val eMyy = eMyy + (2.0*yy*eyy*val)^2 + (yy*yy*eval)^2 ENDIF ENDFOR ENDFOR Mxx = Mxx/M0 eMxx = SQRT( eMxx/(M0*M0) + ((Mxx*eM0)/(M0*M0))^2 ) Myy = Myy/M0 eMyy = SQRT( eMyy/(M0*M0) + ((Myy*eM0)/(M0*M0))^2 ) Mxy = Mxy/M0 eMxy = SQRT( eMxy/(M0*M0) + ((Mxy*eM0)/(M0*M0))^2 ) ; intermediate quantities needed to compute eps, eeps, posang, eposang... uu = Mxx - Myy euu = SQRT( eMxx*eMxx + eMyy*eMyy ) vv = Mxx + Myy evv = euu ww = 2*Mxy eww = 2*eMxy ss = ((uu*euu)^2 + (ww*eww)^2)/(uu*uu + ww*ww) tt = ((uu*uu + ww*ww)*evv*evv)/(vv*vv) ENDIF ENDIF ; AUTOMATED CENTROIDING AND SHAPE/ORIENTATION ASSESSMENT... IF (keyword_set(cenx) AND keyword_set(ceny)) THEN BEGIN ; no automated centroiding is necessary... xcen = cenx ycen = ceny x0 = xcen - 1 y0 = ycen - 1 IF (x0 GT nc OR y0 GT nr) THEN BEGIN IF (x0 GT nc) THEN print,'ERROR: value of xcen > ncol !' IF (y0 GT nr) THEN print,'ERROR: value of ycen > nrow !' RETURN ENDIF print,' adopting user-specified galaxy center position...' print,' (xcen,ycen) =(',xcen,',',ycen,') [fixed]',f='(a,f7.2,a,f7.2,a)' printf,25,' adopting user-specified galaxy center position...' printf,25,' (xcen,ycen) =(',xcen,',',ycen,') [fixed]',f='(a,f7.2,a,f7.2,a)' ENDIF ELSE BEGIN ; automated centroiding requested: determine center coordinates via image ; 1st order moments of image pixels that have mask values >0 in the mask ; and which where pre-computed above (i.e., x0,y0 were defined there). xcen = x0 + 1 ycen = y0 + 1 print,' estimating galaxy center coords from 1st-order image moments...' print,' (xcen,ycen) = (',xcen,'+/-',ex0,',',ycen,'+/-',ey0,')',f='(a,f7.2,a,f4.2,a,f7.2,a,f4.2,a)' printf,25,' estimating galaxy center coords from 1st-order image moments...' printf,25,' (xcen,ycen) = (',xcen,'+/-',ex0,',',ycen,'+/-',ey0,')',f='(a,f7.2,a,f4.2,a,f7.2,a,f4.2,a)' ENDELSE IF (keyword_set(epsilon)) THEN BEGIN ; use explicitly given ellipticity... eps = epsilon ba = 1 - eps IF (keyword_set(eepsilon)) THEN BEGIN eeps = eepsilon ENDIF ELSE BEGIN eeps = 0.00 ENDELSE eba = eeps IF (ba LE 0.0 OR ba GT 1.0) THEN BEGIN print,'ERROR: value of eps outside valid range [0.0 .. 1.0> !' printf,25,'ERROR: value of eps outside valid range [0.0 .. 1.0> !' RETURN ENDIF print,' adopting user-specified galaxy ellipticity...' print,' epsilon= (1 - b/a) = ',eps,' +/- ',eeps,' [fixed]',f='(a,f6.4,a,f6.4,a)' printf,25,' adopting user-specified galaxy ellipticity...' printf,25,' epsilon= (1 - b/a) = ',eps,' +/- ',eeps,' [fixed]',f='(a,f6.4,a,f6.4,a)' ENDIF ELSE BEGIN ; derive ellipticity from 2nd order image moments computed above... ; NOTE: this is a re-implementation of my heritage f77 code 'moments.f'. eps = SQRT( uu*uu + ww*ww )/vv eeps = SQRT( ss + tt )/vv ba = 1 - eps eba = eeps print,' estimating galaxy ellipticity from 2nd-order image moments...' print,' eps = (1 - b/a) = ',eps,' +/- ',eeps,f='(a,f6.4,a,f6.4,a)' printf,25,' estimating galaxy ellipticity from 2nd-order image moments...' printf,25,' eps = (1 - b/a) = ',eps,' +/- ',eeps,f='(a,f6.4,a,f6.4,a)' ENDELSE IF (keyword_set(posang)) THEN BEGIN ; use explicitly given position angle... theta = posang IF (keyword_set(eposang)) THEN BEGIN etheta = eposang ENDIF ELSE BEGIN etheta = 0.00 ENDELSE IF (theta LT 0 OR theta GE 180.0) THEN BEGIN print,'ERROR: value of posang outside valid range [0.0 .. 180.0> degrees!' RETURN ENDIF print,' adopting user-specified galaxy position angle...' print,' posang = ',theta,' +/-',etheta,' (wrt. positive x-axis) [fixed]',f='(a,f5.1,a,f4.1,a)' printf,25,' adopting user-specified galaxy position angle...' printf,25,' posang = ',theta,' +/-',etheta,' (wrt. positive x-axis) [fixed]',f='(a,f5.1,a,f4.1,a)' ENDIF ELSE BEGIN ; derive position angle measured counter-clockwise from the positive x-axis ; from the 2nd order image moments computed above. Beware of quadrants! ; NOTE: this is a re-implementation of my heritage f77 code 'moments.f' (you ; may get the correct answer only most of the time...) ss = 90.0 - (180.0/!DPI)*ATAN((ww+0.5*SQRT(2.0)*eww),(uu-0.5*SQRT(2.0)*euu)) IF (ss GT 90.0) THEN ss = ss - 180.0 tt = 90.0 - (180.0/!DPI)*ATAN((ww-0.5*SQRT(2.0)*eww),(uu+0.5*SQRT(2.0)*euu)) IF (tt GT 90.0) THEN tt = tt - 180.0 theta = 0.5*(180.0/!DPI)*ATAN(ww,uu) etheta = ABS(tt - ss) IF (etheta GT 0.0 AND etheta LT 0.1) THEN etheta = 0.1 IF (etheta GT 90.0) THEN etheta = 180.0 - etheta print,' estimating galaxy position angle from 2nd-order image moments...' print,' posang = ',theta,' +/-',etheta,' deg (wrt. positive x-axis)',f='(a,f5.1,a,f4.1,a)' printf,25,' estimating galaxy position angle from 2nd-order image moments...' printf,25,' posang = ',theta,' +/-',etheta,' deg (wrt. positive x-axis)',f='(a,f5.1,a,f4.1,a)' ENDELSE ; DETERMINE THE DISTANCE FROM GALAXY CENTER TO FARTHEST IMAGE BORDER... tmp = 0.0*findgen(4) tmp[0] = x0 tmp[1] = ncol - x0 tmp[2] = y0 tmp[3] = nrow - y0 rmax = MAX(tmp) IF (idebug) THEN BEGIN printf,25,' DEBUG: tmp = [',tmp,'], rmax = ',rmax,f='(a,4f7.2,a,f6.2)' ENDIF ; box width for tic-tac-toe sky region analysis... boxw = 2.0*MIN(tmp)/3. ; TEST FOR SYSTEMATIC GRADIENTS OR VARIATIONS IN THE SKY BACKGROUND... IF ( NOT keyword_set(fixsky) ) THEN BEGIN print,' testing for residual gradient in sky background level...',f='(a)' printf,25,' testing for residual gradient in sky background level...',f='(a)' IF (idebug) THEN BEGIN printf,25,' DEBUG: boxw = ',boxw,f='(a,i4)' ENDIF bg = FLTARR(9) ebg = bg xbg = bg ybg = bg ; compute x,y positions of octant box centers... lx = 0.5*(0 + (x0-0.5*boxw-1)) cx = 0.5*((x0-0.5*boxw) + (x0+0.5*boxw-1)) rx = 0.5*((x0+0.5*boxw) + nc) ly = 0.5*(0 + (y0-0.5*boxw-1)) cy = 0.5*((y0-0.5*boxw) + (y0+0.5*boxw-1)) uy = 0.5*((y0+0.5*boxw) + nr) ; populate vectors of x and y octant box center positions... xbg[0] = x0 ; CC box, x_cen = x0 ybg[0] = y0 ; y_cen = y0 xbg[1] = lx ; UL box, x_cen ybg[1] = uy ; y_cen xbg[2] = cx ; UC box, x_cen ybg[2] = uy ; y_cen xbg[3] = rx ; UR box, x_cen ybg[3] = uy ; y_cen xbg[4] = rx ; CR box, x_cen ybg[4] = cy ; y_cen xbg[5] = rx ; LR box, x_cen ybg[5] = ly ; y_cen xbg[6] = cx ; LC box, x_cen ybg[6] = ly ; y_cen xbg[7] = lx ; LL box, x_cen ybg[7] = ly ; y_cen xbg[8] = lx ; CL box, x_cen ybg[8] = cy ; y_cen ; create a copy of the input image, in which we can replace any masked out ; regions with the local mean... skg = img ; CC box (set masked pixels to the overall average sky, i.e., 0.00) ... tmpi = img[(x0-0.5*boxw):(x0+0.5*boxw-1),(y0-0.5*boxw):(y0+0.5*boxw-1)] tmpm = msk[(x0-0.5*boxw):(x0+0.5*boxw-1),(y0-0.5*boxw):(y0+0.5*boxw-1)] tmps = tmpi ig = WHERE(tmpm NE 0, ign) tmps[ig] = 0.00 skg[(x0-0.5*boxw):(x0+0.5*boxw-1),(y0-0.5*boxw):(y0+0.5*boxw-1)] = tmps bg[0] = 0.0 & ebg[0] = esky ; UL box... IF (idebug) THEN BEGIN printf,25,' DEBUG: UL box = img[',0,':',(x0-0.5*boxw-1),',',(y0+0.5*boxw),':',nr,']',f='(a,i4,a,i4,a,i4,a,i4,a)' ENDIF tmpi = img[0:(x0-0.5*boxw-1),(y0+0.5*boxw):nr] tmpm = msk[0:(x0-0.5*boxw-1),(y0+0.5*boxw):nr] tmps = tmpi ig = WHERE(tmpm EQ 0, ign) MMM, tmpi(ig), sky1, esky1 it = WHERE(tmpm EQ 0 AND (tmpi GT (sky1-3.*esky1) AND tmpi LT (sky1+3.*esky1)), itn) tmpj = tmpi(it) RESISTANT_Mean, tmpj, 4.0, sky1, esky1, nrej1 printf,25,' bg level in UL octant box =',sky1,' +/-',esky1,' [npix=',(itn-nrej1),']',f='(a,f9.3,a,f7.3,a,i7,a)' bg[1] = sky1 & ebg[1] = esky1 ig = WHERE(tmpm NE 0, ign) tmps[ig] = sky1 skg[0:(x0-0.5*boxw-1),(y0+0.5*boxw):nr] = tmps ; UC box... IF (idebug) THEN BEGIN printf,25,' DEBUG: UC box = img[',(x0-0.5*boxw),':',(x0+0.5*boxw-1),',',(y0+0.5*boxw),':',nr,']',f='(a,i4,a,i4,a,i4,a,i4,a)' ENDIF tmpi = img[(x0-0.5*boxw):(x0+0.5*boxw-1),(y0+0.5*boxw):nr] tmpm = msk[(x0-0.5*boxw):(x0+0.5*boxw-1),(y0+0.5*boxw):nr] tmps = tmpi ig = WHERE(tmpm EQ 0, ign) MMM, tmpi(ig), sky2, esky2 it = WHERE(tmpm EQ 0 AND (tmpi GT (sky2-3.*esky2) AND tmpi LT (sky2+3.*esky2)), itn) tmpj = tmpi(it) RESISTANT_Mean, tmpj, 4.0, sky2, esky2, nrej2 printf,25,' bg level in UC octant box =',sky2,' +/-',esky2,' [npix=',(itn-nrej2),']',f='(a,f9.3,a,f7.3,a,i7,a)' bg[2] = sky2 & ebg[2] = esky2 ig = WHERE(tmpm NE 0, ign) tmps[ig] = sky2 skg[(x0-0.5*boxw):(x0+0.5*boxw-1),(y0+0.5*boxw):nr] = tmps ; UR box... IF (idebug) THEN BEGIN printf,25,' DEBUG: UR box = img[',(x0+0.5*boxw),':',nc,',',(y0+0.5*boxw),':',nr,']',f='(a,i4,a,i4,a,i4,a,i4,a)' ENDIF tmpi = img[(x0+0.5*boxw):nc,(y0+0.5*boxw):nr] tmpm = msk[(x0+0.5*boxw):nc,(y0+0.5*boxw):nr] tmps = tmpi ig = WHERE(tmpm EQ 0, ign) MMM, tmpi(ig), sky3, esky3 it = WHERE(tmpm EQ 0 AND (tmpi GT (sky3-3.*esky3) AND tmpi LT (sky3+3.*esky3)), itn) tmpj = tmpi(it) RESISTANT_Mean, tmpj, 4.0, sky3, esky3, nrej3 printf,25,' bg level in UR octant box =',sky3,' +/-',esky3,' [npix=',(itn-nrej3),']',f='(a,f9.3,a,f7.3,a,i7,a)' bg[3] = sky3 & ebg[3] = esky3 ig = WHERE(tmpm NE 0, ign) tmps[ig] = sky3 skg[(x0+0.5*boxw):nc,(y0+0.5*boxw):nr] = tmps ; CR box... IF (idebug) THEN BEGIN printf,25,' DEBUG: CR box = img[',(x0+0.5*boxw),':',nc,',',(y0-0.5*boxw),':',(y0+0.5*boxw-1),']',f='(a,i4,a,i4,a,i4,a,i4,a)' ENDIF tmpi = img[(x0+0.5*boxw):nc,(y0-0.5*boxw):(y0+0.5*boxw-1)] tmpm = msk[(x0+0.5*boxw):nc,(y0-0.5*boxw):(y0+0.5*boxw-1)] tmps = tmpi ig = WHERE(tmpm EQ 0, ign) MMM, tmpi(ig), sky4, esky4 it = WHERE(tmpm EQ 0 AND (tmpi GT (sky4-3.*esky4) AND tmpi LT (sky4+3.*esky4)), itn) tmpj = tmpi(it) RESISTANT_Mean, tmpj, 4.0, sky4, esky4, nrej4 printf,25,' bg level in CR octant box =',sky4,' +/-',esky4,' [npix=',(itn-nrej4),']',f='(a,f9.3,a,f7.3,a,i7,a)' bg[4] = sky4 & ebg[4] = esky4 ig = WHERE(tmpm NE 0, ign) tmps[ig] = sky4 skg[(x0+0.5*boxw):nc,(y0-0.5*boxw):(y0+0.5*boxw-1)] = tmps ; LR box... IF (idebug) THEN BEGIN printf,25,' DEBUG: LR box = img[',(x0+0.5*boxw),':',nc,',',0,':',(y0-0.5*boxw-1),']',f='(a,i4,a,i4,a,i4,a,i4,a)' ENDIF tmpi = img[(x0+0.5*boxw):nc,0:(y0-0.5*boxw-1)] tmpm = msk[(x0+0.5*boxw):nc,0:(y0-0.5*boxw-1)] tmps = tmpi ig = WHERE(tmpm EQ 0, ign) MMM, tmpi(ig), sky5, esky5 it = WHERE(tmpm EQ 0 AND (tmpi GT (sky5-3.*esky5) AND tmpi LT (sky5+3.*esky5)), itn) tmpj = tmpi(it) RESISTANT_Mean, tmpj, 4.0, sky5, esky5, nrej5 printf,25,' bg level in LR octant box =',sky5,' +/-',esky5,' [npix=',(itn-nrej5),']',f='(a,f9.3,a,f7.3,a,i7,a)' bg[5] = sky5 & ebg[5] = esky5 ig = WHERE(tmpm NE 0, ign) tmps[ig] = sky5 skg[(x0+0.5*boxw):nc,0:(y0-0.5*boxw-1)] = tmps ; LC box... IF (idebug) THEN BEGIN printf,25,' DEBUG: LC box = img[',(x0-0.5*boxw),':',(x0+0.5*boxw-1),',',0,':',(y0-0.5*boxw-1),']',f='(a,i4,a,i4,a,i4,a,i4,a)' ENDIF tmpi = img[(x0-0.5*boxw):(x0+0.5*boxw-1),0:(y0-0.5*boxw-1)] tmpm = msk[(x0-0.5*boxw):(x0+0.5*boxw-1),0:(y0-0.5*boxw-1)] tmps = tmpi ig = WHERE(tmpm EQ 0, ign) MMM, tmpi(ig), sky6, esky6 it = WHERE(tmpm EQ 0 AND (tmpi GT (sky6-3.*esky6) AND tmpi LT (sky6+3.*esky6)), itn) tmpj = tmpi(it) RESISTANT_Mean, tmpj, 4.0, sky6, esky6, nrej6 printf,25,' bg level in LC octant box =',sky6,' +/-',esky6,' [npix=',(itn-nrej6),']',f='(a,f9.3,a,f7.3,a,i7,a)' bg[6] = sky6 & ebg[6] = esky6 ig = WHERE(tmpm NE 0, ign) tmps[ig] = sky6 skg[(x0-0.5*boxw):(x0+0.5*boxw-1),0:(y0-0.5*boxw-1)] = tmps ; LL box... IF (idebug) THEN BEGIN printf,25,' DEBUG: LL box = img[',0,':',(x0-0.5*boxw-1),',',0,':',(y0-0.5*boxw-1),']',f='(a,i4,a,i4,a,i4,a,i4,a)' ENDIF tmpi = img[0:(x0-0.5*boxw-1),0:(y0-0.5*boxw-1)] tmpm = msk[0:(x0-0.5*boxw-1),0:(y0-0.5*boxw-1)] tmps = tmpi ig = WHERE(tmpm EQ 0, ign) MMM, img(ig), sky7, esky7 it = WHERE(tmpm EQ 0 AND (tmpi GT (sky7-3.*esky7) AND tmpi LT (sky7+3.*esky7)), itn) tmpj = tmpi(it) RESISTANT_Mean, tmpj, 4.0, sky7, esky7, nrej7 printf,25,' bg level in LL octant box =',sky7,' +/-',esky7,' [npix=',(itn-nrej7),']',f='(a,f9.3,a,f7.3,a,i7,a)' bg[7] = sky7 & ebg[7] = esky7 ig = WHERE(tmpm NE 0, ign) tmps[ig] = sky7 skg[0:(x0-0.5*boxw-1),0:(y0-0.5*boxw-1)] = tmps ; CL box... IF (idebug) THEN BEGIN printf,25,' DEBUG: CL box = img[',0,':',(x0-0.5*boxw-1),',',(y0-0.5*boxw),':',(y0+0.5*boxw-1),']',f='(a,i4,a,i4,a,i4,a,i4,a)' ENDIF tmpi = img[0:(x0-0.5*boxw-1),(y0-0.5*boxw):(y0+0.5*boxw-1)] tmpm = msk[0:(x0-0.5*boxw-1),(y0-0.5*boxw):(y0+0.5*boxw-1)] tmps = tmpi ig = WHERE(tmpm EQ 0, ign) MMM, img(ig), sky8, esky8 it = WHERE(tmpm EQ 0 AND (tmpi GT (sky8-3.*esky8) AND tmpi LT (sky8+3.*esky8)), itn) tmpj = tmpi(it) RESISTANT_Mean, tmpj, 4.0, sky8, esky8, nrej8 printf,25,' bg level in CL octant box =',sky8,' +/-',esky8,' [npix=',(itn-nrej8),']',f='(a,f9.3,a,f7.3,a,i7,a)' bg[8] = sky8 & ebg[8] = esky8 ig = WHERE(tmpm NE 0, ign) tmps[ig] = sky8 skg[0:(x0-0.5*boxw-1),(y0-0.5*boxw):(y0+0.5*boxw-1)] = tmps ; Report variation (rms) of sky boxes... RESISTANT_Mean, bg, 4.0, skyavg, skyrms, skyrej ;skysig = ROBUST_SIGMA(bg,/ZERO) skysig = ROBUST_SIGMA(bg) print,' sky variation in octant boxes =',skyrms,' ADU (rms), ',skysig,' ADU (ampl)',f='(a,f7.3,a,f7.3,a)' printf,25,' sky variation in octant boxes =',skyrms,' ADU (rms), ',skysig,' ADU (ampl)',f='(a,f7.3,a,f7.3,a)' ; Plot azimuthal variation in sky level... IF (inoplot) THEN BEGIN x = INDGEN(8)+1 ymax = MAX(1.1*(bg+ebg)) ymin = MIN(1.2*(bg-ebg)) ; set_plot, 'PS' ; device, color=1, filename="plotje.eps", /encapsulated,$ ; xsize=6, ysize=6, /Inches, BITS_PER_PIXEL=8 window, 0, TITLE=" GALPROF v0.5: image="+image, XSIZE=768, YSIZE=512 s = FINDGEN(17)*(!PI*2/16.) usersym, COS(s), SIN(s), /FILL plot, x, bg, XRANGE=[0.5,8.5], XSTYLE=1, YRANGE=[ymin,ymax], YSTYLE=1,$ TITLE="Background non-uniformity fit", CHARSIZE=2.0, THICK=2,$ XTITLE="Tic-tac-toe octant box #",XCHARSIZE=1.25,$ YTITLE="Sky background level [ADU]",YCHARSIZE=1.25,/YNOZERO,$ PSYM=-8, SYMSIZE=1.5, /NODATA errplot, x, bg-ebg, bg+ebg, COLOR=(128+255*(128+128*255)) oplot, x, bg, COLOR=(200+255*(128+0*255)), PSYM=-8, SYMSIZE=1.5 t = FLTARR(8) oplot, x, t, LINESTYLE=2, PSYM=-3, THICK=2, SYMSIZE=1, COLOR=255 ; device,/close ENDIF ; FIT RESIDUAL SKY BACKGROUND (SLOPE/GRADIENT) BASED ON OCTANT BOXES... ; Note that the softening scale is hard-coded to '2.0*boxw' FOR row=0,nr DO BEGIN FOR col=0,nc DO BEGIN tmpsky = 0.0 tmpwgh = 0.0 FOR ib=0,8 DO BEGIN dd = 2.0*boxw + SQRT( (col-xbg[ib])^2 + (row-ybg[ib])^2 ) tmpsky = tmpsky + ((2.0*boxw)/dd)^2 * bg[ib] tmpwgh = tmpwgh + ((2.0*boxw)/dd)^2 ENDFOR skg(col,row) = sqrt(9.0)*tmpsky/tmpwgh ENDFOR ENDFOR dsky = skg(xbg[0],ybg[0]) skg = TEMPORARY(skg) - dsky IF (idebug) THEN BEGIN printf,25,' DEBUG: raw delta_sky frame has sky(x0,y0)=',dsky,f='(a,f10.6)' ENDIF ; SUBTRACT RESIDUAL SKY BACKGROUND (SLOPE/GRADIENT)... img = img - skg ; TEST RESULT OF ATTEMPTED REMOVAL OF RESIDUAL BACKGROUND GRADIENTS... print,' 2nd pass, after subtraction of delta-flat...',f='(a)' printf,25,' 2nd pass, after subtraction of delta-flat...',f='(a)' bbg = FLTARR(9) ebbg = bbg ; UL box... tmpi = img[0:(x0-0.5*boxw-1),(y0+0.5*boxw):nr] tmpm = msk[0:(x0-0.5*boxw-1),(y0+0.5*boxw):nr] ig = WHERE(tmpm EQ 0, ign) MMM, tmpi(ig), sky1, esky1 it = WHERE(tmpm EQ 0 AND (tmpi GT (sky1-3.*esky1) AND tmpi LT (sky1+3.*esky1)), itn) tmpj = tmpi(it) RESISTANT_Mean, tmpj, 4.0, sky1, esky1, nrej1 printf,25,' bg level in UL octant box =',sky1,' +/-',esky1,' [npix=',(itn-nrej1),']',f='(a,f9.3,a,f7.3,a,i7,a)' bbg[1] = sky1 & ebbg[1] = esky1 ; UC box... tmpi = img[(x0-0.5*boxw):(x0+0.5*boxw-1),(y0+0.5*boxw):nr] tmpm = msk[(x0-0.5*boxw):(x0+0.5*boxw-1),(y0+0.5*boxw):nr] ig = WHERE(tmpm EQ 0, ign) MMM, tmpi(ig), sky2, esky2 it = WHERE(tmpm EQ 0 AND (tmpi GT (sky2-3.*esky2) AND tmpi LT (sky2+3.*esky2)), itn) tmpj = tmpi(it) RESISTANT_Mean, tmpj, 4.0, sky2, esky2, nrej2 printf,25,' bg level in UC octant box =',sky2,' +/-',esky2,' [npix=',(itn-nrej2),']',f='(a,f9.3,a,f7.3,a,i7,a)' bbg[2] = sky2 & ebbg[2] = esky2 ; UR box... tmpi = img[(x0+0.5*boxw):nc,(y0+0.5*boxw):nr] tmpm = msk[(x0+0.5*boxw):nc,(y0+0.5*boxw):nr] ig = WHERE(tmpm EQ 0, ign) MMM, tmpi(ig), sky3, esky3 it = WHERE(tmpm EQ 0 AND (tmpi GT (sky3-3.*esky3) AND tmpi LT (sky3+3.*esky3)), itn) tmpj = tmpi(it) RESISTANT_Mean, tmpj, 4.0, sky3, esky3, nrej3 printf,25,' bg level in UR octant box =',sky3,' +/-',esky3,' [npix=',(itn-nrej3),']',f='(a,f9.3,a,f7.3,a,i7,a)' bbg[3] = sky3 & ebbg[3] = esky3 ; CR box... tmpi = img[(x0+0.5*boxw):nc,(y0-0.5*boxw):(y0+0.5*boxw-1)] tmpm = msk[(x0+0.5*boxw):nc,(y0-0.5*boxw):(y0+0.5*boxw-1)] ig = WHERE(tmpm EQ 0, ign) MMM, tmpi(ig), sky4, esky4 it = WHERE(tmpm EQ 0 AND (tmpi GT (sky4-3.*esky4) AND tmpi LT (sky4+3.*esky4)), itn) tmpj = tmpi(it) RESISTANT_Mean, tmpj, 4.0, sky4, esky4, nrej4 printf,25,' bg level in CR octant box =',sky4,' +/-',esky4,' [npix=',(itn-nrej4),']',f='(a,f9.3,a,f7.3,a,i7,a)' bbg[4] = sky4 & ebbg[4] = esky4 ; LR box... tmpi = img[(x0+0.5*boxw):nc,0:(y0-0.5*boxw-1)] tmpm = msk[(x0+0.5*boxw):nc,0:(y0-0.5*boxw-1)] ig = WHERE(tmpm EQ 0, ign) MMM, tmpi(ig), sky5, esky5 it = WHERE(tmpm EQ 0 AND (tmpi GT (sky5-3.*esky5) AND tmpi LT (sky5+3.*esky5)), itn) tmpj = tmpi(it) RESISTANT_Mean, tmpj, 4.0, sky5, esky5, nrej5 printf,25,' bg level in LR octant box =',sky5,' +/-',esky5,' [npix=',(itn-nrej5),']',f='(a,f9.3,a,f7.3,a,i7,a)' bbg[5] = sky5 & ebbg[5] = esky5 ; LC box... tmpi = img[(x0-0.5*boxw):(x0+0.5*boxw-1),0:(y0-0.5*boxw-1)] tmpm = msk[(x0-0.5*boxw):(x0+0.5*boxw-1),0:(y0-0.5*boxw-1)] ig = WHERE(tmpm EQ 0, ign) MMM, tmpi(ig), sky6, esky6 it = WHERE(tmpm EQ 0 AND (tmpi GT (sky6-3.*esky6) AND tmpi LT (sky6+3.*esky6)), itn) tmpj = tmpi(it) RESISTANT_Mean, tmpj, 4.0, sky6, esky6, nrej6 printf,25,' bg level in LC octant box =',sky6,' +/-',esky6,' [npix=',(itn-nrej6),']',f='(a,f9.3,a,f7.3,a,i7,a)' bbg[6] = sky6 & ebbg[6] = esky6 ; LL box... tmpi = img[0:(x0-0.5*boxw-1),0:(y0-0.5*boxw-1)] tmpm = msk[0:(x0-0.5*boxw-1),0:(y0-0.5*boxw-1)] ig = WHERE(tmpm EQ 0, ign) MMM, img(ig), sky7, esky7 it = WHERE(tmpm EQ 0 AND (tmpi GT (sky7-3.*esky7) AND tmpi LT (sky7+3.*esky7)), itn) tmpj = tmpi(it) RESISTANT_Mean, tmpj, 4.0, sky7, esky7, nrej7 printf,25,' bg level in LL octant box =',sky7,' +/-',esky7,' [npix=',(itn-nrej7),']',f='(a,f9.3,a,f7.3,a,i7,a)' bbg[7] = sky7 & ebbg[7] = esky7 ; CL box... tmpi = img[0:(x0-0.5*boxw-1),(y0-0.5*boxw):(y0+0.5*boxw-1)] tmpm = msk[0:(x0-0.5*boxw-1),(y0-0.5*boxw):(y0+0.5*boxw-1)] ig = WHERE(tmpm EQ 0, ign) MMM, img(ig), sky8, esky8 it = WHERE(tmpm EQ 0 AND (tmpi GT (sky8-3.*esky8) AND tmpi LT (sky8+3.*esky8)), itn) tmpj = tmpi(it) RESISTANT_Mean, tmpj, 4.0, sky8, esky8, nrej8 printf,25,' bg level in CL octant box =',sky8,' +/-',esky8,' [npix=',(itn-nrej8),']',f='(a,f9.3,a,f7.3,a,i7,a)' bbg[8] = sky8 & ebbg[8] = esky8 ; Report variation (rms) of sky boxes... RESISTANT_Mean, bbg, 4.0, sky2avg, sky2rms, sky2rej sky2sig = ROBUST_SIGMA(bbg) print,' sky variation in octant boxes =',sky2rms,' ADU (rms), ',sky2sig,' ADU (ampl)',f='(a,f7.3,a,f7.3,a)' printf,25,' sky variation in octant boxes =',sky2rms,' ADU (rms), ',sky2sig,' ADU (ampl)',f='(a,f7.3,a,f7.3,a)' ; Plot azimuthal variation in sky level after subtraction of delta-flat... IF (inoplot) THEN BEGIN errplot, x, bbg-ebbg, bbg+ebbg, COLOR=(0+255*(255+0*255)) oplot, x, bbg, COLOR=(200+255*(255+0*255)), PSYM=-8, SYMSIZE=1.5 ENDIF ; Write out a sky background image (incl. inferred gradients therein) as ; estimated from the analysis of the background in the octant sky boxes... skgout = STRMID(result,0,STRPOS(result,'.',/REVERSE_SEARCH))+'_sky.fits' printf,25,' saving image "'+skgout+'" of inferred sky background...',f='(a)' skghead = maskhead GET_DATE,dte,/TIMETAG SXADDPAR,skghead,'BITPIX',-32,' Bits per data value' SXADDPAR,skghead,'EXTEND','F',' File may contain standard extensions' SXADDPAR,skghead,'FILENAME',skgout,' Frame title' SXADDPAR,skghead,'DATE',dte,' Date this file was written' SXADDPAR,skghead,'IMAGETYP','SKY',' Sky background image' SXADDPAR,skghead,'SKYAVG',sky,' [ADU] Average sky background level' SXADDPAR,skghead,'ESKYAVG',esky,' [ADU] Uncertainty thereon' SXADDPAR,skghead,'SKYRMS',sky2rms,' [ADU] RMS (int) after sky subtraction' SXADDPAR,skghead,'SKYSIG',sky2sig,' [ADU] RMS (ext) after sky subtraction' SXADDPAR,skghead,'HISTORY','New copy of '+image+' ['+STRMID(progver,0,11)+']' skg = TEMPORARY(skg) + sky WRITEFITS,skgout,skg,skghead ENDIF ; CONSTRUCT A 2-D ARRAY OF MAJOR AXIS ELLIPTICAL RADII FROM (x0,y0)... ; To determine the major axis elliptical radius, a, for any pixel (i,j), ; we apply a coordinate transformation (x,y)-->(v,w) aligned with the ; cardinal axes of the ellipse. print,' constructing 2-D array of major-axis elliptical distances...' printf,25,' constructing 2-D array of major-axis elliptical distances...' eqr = 0.0*FLOAT(msk) FOR row=0,nr DO BEGIN FOR col=0,nc DO BEGIN x = FLOAT(col) - x0 y = FLOAT(row) - y0 ; transform to rotated coordinate system... v = y*COS(theta*dtr) - x*SIN(theta*dtr) w = x*COS(theta*dtr) + y*SIN(theta*dtr) ; compute radius at (v,w), given b = a*(1-eps) ... eqr(col,row) = SQRT(v*v+w*w*ba*ba)/ba IF (idebug AND eqr(col,row) LT rminr ) THEN BEGIN printf,25,' DEBUG: eqr(',(col+1),',',(row+1),') = ',eqr(col,row) ENDIF ENDFOR ENDFOR IF (idebug) THEN BEGIN ; NOTE that in the logfile we refer to the pixel coordinates following the ; FITS standard with pixel (1,1) the lower left corner, while the IDL ; array coordinates start at (0,0). Hence the apparent off-by-one in ; the printf statements below. RAJ. printf,25,' DEBUG: saving major axis elliptical distance mask image...',f='(a)' printf,25,' DEBUG: (verify result above) eqr(',ROUND(x0+1),',',ROUND(y0+1),') =',eqr(ROUND(x0),ROUND(y0)),f='(a,i4,a,i4,a,f13.8)' eqrhead = maskhead GET_DATE,dte,/TIMETAG SXADDPAR,eqrhead,'BITPIX',-32,' Bits per data value' SXADDPAR,eqrhead,'EXTEND','F',' File may contain standard extensions' SXADDPAR,eqrhead,'FILENAME','EQR_array',' Frame title' SXADDPAR,eqrhead,'DATE',dte,' Date this file was written' eqrout = STRMID(result,0,STRPOS(result,'.',/REVERSE_SEARCH))+'_eqr.fits' WRITEFITS,eqrout,eqr,eqrhead tst = READFITS(eqrout,eqrhead,/silent) printf,25,' DEBUG: read-back test: eqr(',ROUND(x0+1),',',ROUND(y0+1),') = ',tst(ROUND(x0),ROUND(y0)),f='(a,i4,a,i4,a,f13.8)' ENDIF ; CONSTRUCT ARRAY WITH (EQUIVALENT) RADII... print,' populating 1-D array of major axis radii for aperture photometry...' printf,25,' populating 1-D array of major axis radii for aperture photometry...' ; DETERMINE OUT TO WHICH RADIUS TO USE THE FIXED SET OF CIRCULAR RADII... fixrad = [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0 ] in = WHERE( fixrad LT rminr ) print,' radii = [',MIN(fixrad(in)),'..',MAX(fixrad(in)),'],',$ rminr,'..',rmax,' pix',f='(a,f6.4,a,f6.4,a,f6.4,a,f6.2,a)' printf,25,' radii = [',MIN(fixrad(in)),'..',MAX(fixrad(in)),'],',$ rminr,'..',rmax,' pix',f='(a,f6.4,a,f6.4,a,f6.4,a,f6.2,a)' radii = 0.0*fltarr(128) in = WHERE( fixrad LE rminr , ii) radii[0] = fixrad[in] IF (idebug) THEN BEGIN printf,25,' DEBUG: fixrad=',fixrad printf,25,' DEBUG: growing radii from element ii=',ii,f='(a,i2)' ENDIF FOR imax=ii,127,1 DO BEGIN tmpr = radii[imax-1]*dr IF (tmpr GT rmax) THEN BREAK radii[imax] = tmpr ENDFOR ; Upon exit of this FOR-loop, 'imax' gives last radius element + 1 imax = imax - 1 IF (idebug) THEN BEGIN printf,25,' DEBUG: imax = ',imax,', radii[imax] = ',radii[imax],' pix',f='(a,i3,a,f8.3,a)' ENDIF ; PERFORM ELLIPTICAL APERTURE PHOTOMETRY... print,' performing elliptical aperture photometry...' printf,25,' performing elliptical aperture photometry...' pp = 0.0*fltarr(imax+1) rr = 0.0*dblarr(imax+1) flux = 0.0*dblarr(imax+1) eflux = flux sb = flux esb = flux esbb = flux esbf = flux maxrad= 0.0 FOR i=0,imax,1 DO BEGIN ; Estimate maximum number of pixels within current annulus... IF (i EQ 0) THEN BEGIN ii = WHERE(msk GT 0 AND eqr LE radii[i], npix) nmax = !DPI*( ba*(radii[i])^2 ) + 1.0 ENDIF ELSE BEGIN ; We select on "msk GE 0" rather than "msk GT 0", since we don't want to ; limit ourselves strictly to the region marked in the mask, which ; really only served to define the sky background... ii = WHERE(msk GE 0 AND eqr GT radii[i-1] AND eqr LE radii[i], npix) nmax = !DPI*( ba*( (radii[i])^2 - (radii[i-1])^2 ) ) ENDELSE pfrac = FLOAT(npix)/nmax ; Check whether we have any valid galaxy pixels within aperture... IF (npix EQ 0) THEN BEGIN print,'WARNING: no valid pixels at r1->r2='+STRTRIM(STRING(radii[i],FORMAT='(f8.2)'),2)+' !',f='(a)' printf,25,'WARNING: no valid pixels at r1->r2='+STRTRIM(STRING(radii[i],FORMAT='(f8.2)'),2)+' !',f='(a)' ENDIF ; Check whether the elliptical annulus contains at least 50% of the pixels ; that would theoretically fall within that annulus... IF ( pfrac LT 0.50 ) THEN BEGIN print,' WARNING: <50% valid pixels at r1->r2='+STRTRIM(STRING(radii[i],FORMAT='(f8.2)'),2)+' pix (npix='+STRTRIM(STRING(npix,FORMAT='(i6)'),2)+', nmax='+STRTRIM(STRING(nmax,FORMAT='(f8.2)'),2)+')!',f='(a)' printf,25,' WARNING: <50% valid pixels at r1->r2='+STRTRIM(STRING(radii[i],FORMAT='(f8.2)'),2)+' pix (npix='+STRTRIM(STRING(npix,FORMAT='(i6)'),2)+', nmax='+STRTRIM(STRING(nmax,FORMAT='(f8.2)'),2)+')!',f='(a)' ENDIF IF ( pfrac GT 1.00 ) THEN BEGIN pfrac = 1.00 IF ( (npix-nmax) GT 5 AND ((npix-nmax)/nmax) GT 0.05 ) THEN BEGIN print,' WARNING: >100% of expected pixels at r1->r2='+STRTRIM(STRING(radii[i],FORMAT='(f8.2)'),2)+' pix (npix='+STRTRIM(STRING(npix,FORMAT='(i6)'),2)+', nmax='+STRTRIM(STRING(nmax,FORMAT='(f8.2)'),2)+')!',f='(a)' printf,25,' WARNING: >100% of expected pixels at r1->r2='+STRTRIM(STRING(radii[i],FORMAT='(f8.2)'),2)+' pix (npix='+STRTRIM(STRING(npix,FORMAT='(i6)'),2)+', nmax='+STRTRIM(STRING(nmax,FORMAT='(f8.2)'),2)+')!',f='(a)' ENDIF ENDIF ; store the number of pixels in the current annulus, and store the current ; equivalent elliptical radius sqrt(ab) = a*sqrt(1-eps) pp[i] = FLOAT(npix) rr[i] = SQRT(ba)*TOTAL(eqr(ii),/DOUBLE)/pp[i] flux[i] = TOTAL(img(ii),/DOUBLE) IF (flux[i] LT 0.0D0 ) THEN BEGIN print,' WARNING: negative flux set to 0.0 at r1->r2='+STRTRIM(STRING(radii[i],FORMAT='(f8.2)'),2)+' pix!',f='(a)' printf,25,' WARNING: negative flux set to 0.0 at r1->r2='+STRTRIM(STRING(radii[i],FORMAT='(f8.2)'),2)+' pix!',f='(a)' flux[i] = 0.0D0 ENDIF ; Add error propagation (CCD equation and sky subtraction uncertainty).. Sobj = flux[i] Ssky = sky Sdark = texp*dc IF (idebug) THEN BEGIN IF (i EQ 0) THEN BEGIN printf,25,' DEBUG: ron =',ron,' e-, eadu=',eadu,' e-/ADU, texp=',texp,' s',f='(a,f5.2,a,f5.2,a,i4,a)' printf,25,' DEBUG: i=',i,', rr[i]=',rr[i],', radii[i]=',radii[i],', Sobj=',Sobj,' e-, Ssky=',Ssky,' e-, Sdark=',Sdark,' e-',f='(a,i2,a,f6.2,a,f6.2,a,f12.4,a,f8.3,a,f7.3,a)' ENDIF ELSE BEGIN printf,25,' DEBUG: i=',i,', rr[i]=',rr[i],', radii[i]=',radii[i],', Sobj=',Sobj,' e-',f='(a,i2,a,f6.2,a,f6.2,a,f12.4,a)' ENDELSE ENDIF etmp = SQRT(Sobj+pp[i]*(1.0+pp[i]/FLOAT(nsky))*(Ssky+Sdark+ron*ron)) ;eflux[i] = SQRT(etmp*etmp + skyrms*skyrms) eflux[i] = SQRT(pfrac)*SQRT(etmp*etmp + skyrms*skyrms) maxrad = radii[i] ENDFOR ii = WHERE(msk GE 0 AND eqr LE maxrad, npix) jj = WHERE(eqr LE maxrad, mpix) pfrac = FLOAT(npix)/FLOAT(mpix) ; WRITE RESULTS TO OUTPUT TABLE... OPENW,26,result printf,26,'# r r1 r2 r2p npix flux eflux mu emu emub emuf',f='(a)' CLOSE,26 OPENW,26,result,/APPEND printf,26,'#',f='(a)' printf,26,'# '+progver+': '+SYSTIME(/UTC),f='(a)' printf,26,'# image ="'+image+'"',f='(a)' printf,26,'# mask ="'+mask+'"',f='(a)' printf,26,'#',f='(a)' cflg = 'fit' & IF (keyword_set(cenx) AND keyword_set(ceny)) THEN cflg = 'fixed' IF (keyword_set(cenx)) THEN ex0 = 0.0 IF (keyword_set(ceny)) THEN ey0 = 0.0 printf,26,'# center = (',(x0+1),'+/-',ex0,',',(y0+1),'+/-',ey0,') (FITS pixel coordinates) / ',cflg,f='(a,f6.2,a,f4.2,a,f6.2,a,f4.2,a,a)' cflg = 'fit' & IF (keyword_set(epsilon)) THEN cflg = 'fixed' printf,26,'# epsilon= ',(1.0-ba),' +/- ',eba,' (eps = 1 - b/a) / ',cflg,f='(a,f7.4,a,f6.4,a,a)' cflg = 'fit' & IF (keyword_set(posang)) THEN cflg = 'fixed' printf,26,'# posang = ',theta,' +/- ',etheta,' deg (wrt. positive x-axis) / ',cflg,f='(a,f7.1,a,f6.1,a,a)' cflg = 'user' printf,26,'# pmzero = ',pz,' +/- ',epz,' mag / ',cflg,f='(a,f7.4,a,f6.4,a,a)' printf,26,'# pairmas= ',pam,' +/- ',epam,' mag/A.M. / ',cflg,f='(a,f7.4,a,f6.4,a,a)' printf,26,'# pcolor =',pc,' +/-',epc,' mag/mag, color = ',clr,' / ',cflg,f='(a,f8.5,a,f7.5,a,f6.3,a,a)' printf,26,'# secpix = ',pscl,' +/- ',epscl,' "/pix / ',cf4,f='(a,f7.4,a,f6.4,a,a)' printf,26,'# airmass= ',airm,' +/- ',eairm,' / ',cf5,f='(a,f7.3,a,f5.3,a,a)' printf,26,'# rdnoise= ',ron,' e- / ',cf1,f='(a,f7.2,a,a)' printf,26,'# gain = ',eadu,' e-/ADU / ',cf2,f='(a,f7.2,a,a)' printf,26,'# exptime= ',texp,' sec / ',cf3,f='(a,i7,a,a)' printf,26,'#',f='(a)' IF ( keyword_set(fixsky) ) THEN BEGIN cflg = 'user' printf,26,'# sky =',sky,' +/-',esky,' ADU / ',cflg,f='(a,f8.3,a,f6.3,a,a)' ENDIF ELSE BEGIN cflg = 'fit' printf,26,'# sky_avg=',sky,' +/-',esky,' ADU [nsky=',nsky,' pix] / ',cflg,f='(a,f8.3,a,f6.3,a,i7,a,a)' printf,26,'# 1st-pass:',f='(a)' printf,26,'# r1skyUL= ',bg[1],' +/-',ebg[1],' ADU',f='(a,f7.3,a,f6.3,a)' printf,26,'# r1skyUC= ',bg[2],' +/-',ebg[2],' ADU',f='(a,f7.3,a,f6.3,a)' printf,26,'# r1skyUR= ',bg[3],' +/-',ebg[3],' ADU',f='(a,f7.3,a,f6.3,a)' printf,26,'# r1skyCR= ',bg[4],' +/-',ebg[4],' ADU',f='(a,f7.3,a,f6.3,a)' printf,26,'# r1skyLR= ',bg[5],' +/-',ebg[5],' ADU',f='(a,f7.3,a,f6.3,a)' printf,26,'# r1skyLC= ',bg[6],' +/-',ebg[6],' ADU',f='(a,f7.3,a,f6.3,a)' printf,26,'# r1skyLL= ',bg[7],' +/-',ebg[7],' ADU',f='(a,f7.3,a,f6.3,a)' printf,26,'# r1skyCL= ',bg[8],' +/-',ebg[8],' ADU',f='(a,f7.3,a,f6.3,a)' printf,26,'# rskyrms= ',skyrms,' ADU',f='(a,f7.3,a)' printf,26,'# rskysig= ',skysig,' ADU',f='(a,f7.3,a)' printf,26,'# 2nd-pass:',f='(a)' printf,26,'# r2skyUL= ',bbg[1],' +/-',ebbg[1],' ADU',f='(a,f7.3,a,f6.3,a)' printf,26,'# r2skyUC= ',bbg[2],' +/-',ebbg[2],' ADU',f='(a,f7.3,a,f6.3,a)' printf,26,'# r2skyUR= ',bbg[3],' +/-',ebbg[3],' ADU',f='(a,f7.3,a,f6.3,a)' printf,26,'# r2skyCR= ',bbg[4],' +/-',ebbg[4],' ADU',f='(a,f7.3,a,f6.3,a)' printf,26,'# r2skyLR= ',bbg[5],' +/-',ebbg[5],' ADU',f='(a,f7.3,a,f6.3,a)' printf,26,'# r2skyLC= ',bbg[6],' +/-',ebbg[6],' ADU',f='(a,f7.3,a,f6.3,a)' printf,26,'# r2skyLL= ',bbg[7],' +/-',ebbg[7],' ADU',f='(a,f7.3,a,f6.3,a)' printf,26,'# r2skyCL= ',bbg[8],' +/-',ebbg[8],' ADU',f='(a,f7.3,a,f6.3,a)' printf,26,'# rskyrms= ',sky2rms,' ADU',f='(a,f7.3,a)' printf,26,'# rskysig= ',sky2sig,' ADU',f='(a,f7.3,a)' ENDELSE printf,26,'#',f='(a)' printf,26,'# r r1 r2 r2p npix flux eflux mu emu emub emuf ',f='(a)' printf,26,'# ["] ["] ["] [pix] [pix] [e-] [e-] ------[mag/arcsec^2]------',f='(a)' printf,26,'#==================================================================================',f='(a)' ts = 2.5*ALOG10((sky-esky)/sky) IF (idebug) THEN BEGIN printf,25,' DEBUG: sky=',sky,', esky=',esky,' e-; ts=',ts,f='(a,f7.3,a,f6.3,a,f10.4)' ENDIF IF (inoreg) THEN BEGIN ellreg = STRMID(result,0,STRPOS(result,'.',/REVERSE_SEARCH))+'_ellaper.reg' OPENW,27,ellreg printf,27,'# Region file format: DS9 version 4.1',f='(a)' printf,27,'global dashlist=8 3 width=1 font="helvetica 10 normal roman" select=1 highlite=1 fixed=0 edit=1 move=1 delete=1 include=1 source=1',f='(a)' printf,27,'image',f='(a)' ENDIF FOR i=0,imax,1 DO BEGIN IF (radii[i] LE maxrad) THEN BEGIN ; Apply photometric calibration (ignoring any color term)... IF ( flux[i] LE 0.0 ) THEN flux[i] = 0.0001 sb[i] = -2.5*ALOG10(flux[i]/(pp[i]*texp*pscl*pscl)) + pz - pam*airm ; propagate error on surface brightness profile... ; First, test whether small error approximation holds. If so, compute ; upper and lower errors, otherwise compute formal error on decimal ; logarithm 2.5*d/dx(lg(x)) = (2.5/ln(10))*(dx/x) te = 1.0857362*(eflux[i]/flux[i]) IF (te GT 3.0) THEN te = 3.0 IF ( (eflux[i]/flux[i]) LT 0.2 ) THEN BEGIN ; small error approximation holds... teb = 2.5*ALOG10((flux[i]+eflux[i])/flux[i]) IF (teb GT 3.0) THEN teb = 3.0 IF (flux[i] GT eflux[i]) THEN BEGIN tef = 2.5*ALOG10((flux[i]-eflux[i])/flux[i]) IF (tef GT 3.0) THEN tef = 3.0 ENDIF ELSE BEGIN tef = 3.0 ENDELSE ENDIF ELSE BEGIN teb = te IF (flux[i] GT eflux[i]) THEN BEGIN tef = 2.5*ALOG10((flux[i]-eflux[i])/flux[i]) IF (tef GT 3.0) THEN tef = 3.0 ENDIF ELSE BEGIN tef = 3.0 ENDELSE ENDELSE esb[i] = SQRT(te*te + epz*epz + epam*epam + ts*ts) esbb[i] = SQRT(teb*teb + epz*epz + epam*epam + ts*ts) esbf[i] = SQRT(tef*tef + epz*epz + epam*epam + ts*ts) IF (idebug) THEN BEGIN printf,25,' DEBUG: te=',te,', teb=',teb,', tef=',tef,f='(a,f10.4,a,f10.4,a,f10.4)' ENDIF ; append results for the current radius to output profile table... IF (i EQ 0 ) THEN BEGIN printf,26,pscl*rr[i],0.0,pscl*radii[i],radii[i],pp[i],$ flux[i],eflux[i],sb[i],esb[i],esbb[i],esbf[i],$ f='(f6.2,f7.2,f7.2,f7.2,i6,f13.4,f10.4,f8.3,f6.3,f6.3,f6.3)' ENDIF ELSE BEGIN printf,26,pscl*rr[i],pscl*radii[i-1],pscl*radii[i],radii[i],pp[i],$ flux[i],eflux[i],sb[i],esb[i],esbb[i],esbf[i],$ f='(f6.2,f7.2,f7.2,f7.2,i6,f13.4,f10.4,f8.3,f6.3,f6.3,f6.3)' ENDELSE ; append ellipses corresponding to "r" and "r2" (in original pixel ; coordinates) to the output DS9 region file... IF (inoreg AND radii[i] GE rminr) THEN BEGIN printf,27,'ellipse('+STRTRIM(STRING(x0+1),2)+','+STRTRIM(STRING(y0+1),2)+','+STRTRIM(STRING(radii[i]),2)+','+STRTRIM(STRING(ba*radii[i]),2)+','+STRTRIM(STRING(theta),2)+') # color=blue dash=0',f='(a)' printf,27,'ellipse('+STRTRIM(STRING(x0+1),2)+','+STRTRIM(STRING(y0+1),2)+','+STRTRIM(STRING(rr[i]/SQRT(ba)),2)+','+STRTRIM(STRING(SQRT(ba)*rr[i]),2)+','+STRTRIM(STRING(theta),2)+') # color=green dash=2',f='(a)' ENDIF ENDIF ENDFOR CLOSE,27 CLOSE,26 ; Plot radial surface brightness profile... IF (inoplot) THEN BEGIN ; set_plot, 'PS' ; device, color=1, filename="plotje.eps", /encapsulated,$ ; xsize=6, ysize=6, /Inches, BITS_PER_PIXEL=8 window, 1, TITLE=" GALPROF v0.5: image="+image+" [2]", XSIZE=768,YSIZE=512 s = FINDGEN(17)*(!PI*2/16.) usersym, COS(s), SIN(s), /FILL ii = WHERE( rr LE maxrad ) xm = (pscl*maxrad) miny = MAX(sb(ii)) + 0.75 maxy = MIN(sb(ii)) - 0.75 sbp = sb - esb sbm = sb + esb plot, pscl*rr(ii), sb(ii), XRANGE=[0.0,xm], XSTYLE=1,$ YRANGE=[miny,maxy], YSTYLE=1,$ TITLE="Radial surface brightness profile", CHARSIZE=2.0, THICK=2,$ XTITLE="Equivalent radius, \sqrt(ab)", XCHARSIZE=1.25,$ YTITLE="Surface brightness, mu [mag/arcsec^2]", YCHARSIZE=1.25,$ /YNOZERO, PSYM=-8, SYMSIZE=1.5, /NODATA errplot, pscl*rr(ii), sbp(ii), sbm(ii), COLOR=(128+255*(128+128*255)) oplot, pscl*rr(ii), sb(ii), COLOR=(200+255*(128+0*255)), PSYM=-8, SYMSIZE=1.0 ; device,/close ENDIF ; CLOSING MESSAGE... print,' GALPROF: finished at '+SYSTIME(/UTC) printf,25,' GALPROF: finished at '+SYSTIME(/UTC) CLOSE,25,/ALL,/FORCE RETURN END ;---------------------------------- R A J ----------------------------------