;; ---------------------------------------------------------------------------------------------------------------------------------- ;; NAME: SDSS_TO_APER ;; ;; PURPOSE: ;; Procedure predicts aperture photometry given SDSS model fits, SDSS catalog data, and a target PSF. The goal is to produce ;; matched photometry between SDSS and other data sets. This is an SDSS-specific application of SynMags. ;; ;; CALLING SEQUENCE: ;; sdss_to_aper, DevMag, DevRad, DevAB, ExpMag, ExpRad, ExpAB, fracDeV, PSF_FWHM_array, aper_radius, mag_aperture, mag_aperture_Dev, mag_aperture_Exp, $ ;; PIXEL_SIZE = PIXEL_SIZE, MODELMAG=MODELMAG, MOCK_IMAGE = MOCK_IMAGE, I0_Dev = I0_Dev, I0_Exp = I0_Exp, $ ;; mge_exp = mge_Exp, mge_Dev = mge_Dev, mgeout = mgeout, hdr_mgeout = hdr_mgeout, not_SDSS=not_SDSS ;; ARGUMENTS: ;; SDSS Model Values: Input SDSS values in a single band corresponding to: DevMag, DevRad, DevAB, ExpMag, ExpRad, ExpAB, fracDeV ;; PSF_FWHM_array: Vector of Gaussian FWHM seeing estimates (arcsec) to be convolved with predicted aperture mags ;; aper_radius: Radii of apertures in arcsec, if more than one, then all apertures are predicted for all magnitudes ;; mag_aperture: Best synthetic aperture mags. If /MOCK_IMAGE is not set, then it is either Dev or Exp look-up value, depending on which is better ;; If /MOCK_IMAGE is set, then it is the magnitude output of the mocked image PSF convolution ;; mag_aperture_Dev: Output synthetic aperture magnitudes using the look-up table. Based on the Dev profile. ;; mag_aperture_Exp: Output synthetic aperture magnitudes using the look-up table. Based on the Exp profile. ;; ;; ;; KEYWORDS: ;; MODELMAG: If set to an array of MODELMAGS, then the output mock aperture mags are computed only for the model component used ;; for the SDSS MODELMAG ;; MOCK_IMAGE: Perform aperture photometry on a mock image with the PSF convolved. Output stored in Mag_aperture. ;; PIXEL_SIZE: Pixel size of simulated image on which the SDSS model is PSF-convolved before aperture photometry (default ;; is 200 pixels) ;; ;; ;; HISTORY: ;; Created by KBundy 6/25/12 - based on sdss_to_aper6.pro ;; For 2D Gaussian integral in Cartesian coordinates: Define 1D integrand function gauss2D_1Dint, x common gauss2D_common, A0, sigma, sigma_y, qgauss, Rlim n_sum = n_elements(A0) if n_sum EQ 1 then begin if qgauss GT 0 then integrand = 2* A0 * sigma * qgauss * sqrt(!pi/2.0) * exp( -x^2/(2.0*sigma^2)) * erf(sqrt( (Rlim^2-x^2)/(2.0*sigma^2*qgauss^2) )) $ else integrand = 2*A0 * sigma_y * sqrt(!pi/2.0) * exp( -x^2/(2.0*sigma^2)) * erf(sqrt( (Rlim^2-x^2)/(2.0*sigma_y^2) )) endif if n_sum GT 1 then begin integrand = total(2*A0 * sigma_y * sqrt(!pi/2.0) * exp( -x^2/(2.0*sigma^2)) * erf(sqrt( (Rlim^2-x^2)/(2.0*sigma_y^2) ))) endif return, integrand end ;; For 2D Sersic integral: function returns angle integration limits function thetalims, r return, [0, 2*!pi] end ;; For 2D Gauss integral in polar coordinates: Define integrand function gauss2D, r, theta common gauss2D_common, A0, sigma, sigma_y, qgauss, Rlim integrand = A0 * r * exp(-(r^2/(2*sigma^2))*(1 + (1.0/qgauss^2 - 1)*(sin(theta))^2)) return, integrand end ;; For 2D Sersic flux integral: Define integrand function sersic2D, r, theta common sersic2D_common, I0, Re, n, bn, q ;integrand = I0 * r*exp(-bn * (r/Re)^(1.0/n)*(1 + (1.0/q^2 - 1)*(sin(theta))^2)^(1./(2.0*n))) integrand = I0 * r*exp(-bn * ((r/Re)^(1.0/n)*(1 + (1.0/q^2 - 1)*(sin(theta))^2)^(1./(2.0*n)) - 1)) ;; Makes I0 in Ie (surface brightness at Re) return, integrand end ;; Implicit equation for bn (deV profile) function sersic_bn, bn_arr common sersic2D_common, I0, Re, n, bn, q return, gamma(2.0*n) - 2*igamma(2.0*n, bn_arr)*gamma(2.0*n) end ;; SDSS: Smoothly truncated exponential profile, 2D polar coordinates ;; R: Radius, physical units (same units as those defining Re) ;; theta: Polar angle, radians ;; function fexpcut, r, theta common sdss_profs, Ie_Exp, q_Exp, Re_Exp, Ie_Dev, q_Dev, Re_Dev ;; Ie is surface brightness at Re if n_elements(q_exp) EQ 0 then q_Exp = 1.0 q_Exp = float(q_Exp) rre = r / float(Re_exp) rvec = rre * sqrt((1 + (1.0/q_Exp^2 - 1)*(sin(theta)^2))) ;; Positional vector in polar coordinates, assuming y-direction inclination axis ratio, q EXPFAC = -1.67835D ;; Defines Ie_Exp and Re fexp = Ie_Exp * exp(ExpFac * (rvec-1.0)) ExpOut = 4.0 ;; Profile goes to zero at this radius (units of Re) ExpCut = 3.0*ExpOut/4.0 ;; Profile truncation begins at this radius (units of Re) wr_cutoff = where(rre GT ExpCut, nr_cutoff) wr_zero = where(rre GT ExpOut, nr_zero) fexpcut = fexp if nr_cutoff GT 0 then fexpcut[wr_cutoff] = fexpcut[wr_cutoff] * ( 1 - ((rvec[wr_cutoff] - ExpCut)/(ExpOut - ExpCut))^2)^2.0 ;;if nr_cutoff GT 0 then fexpcut[wr_cutoff] = fexpcut[wr_cutoff] * ( 1 - ((rvec[wr_cutoff] - ExpCut)/(ExpOut - ExpCut))^2)^1.0 ;; NOTE: To reproduce phFitobj.h, last exponent must be 1.0 (not 2.0) if nr_zero GT 0 then fexpcut[wr_zero] = 0 return, fexpcut end ;; SDSS Exp with smoother truncation: Provides better MGE fits to SDSS-like profiles ;; R: Radius, physical units (same units as those defining Re) ;; theta: Polar angle, radians ;; function fexpcut2, r, theta common sdss_profs, Ie_Exp, q_Exp, Re_Exp, Ie_Dev, q_Dev, Re_Dev ;; Ie is surface brightness at Re if n_elements(q_exp) EQ 0 then q_Exp = 1.0 q_Exp = float(q_Exp) rre = r / float(Re_exp) rvec = rre * sqrt((1 + (1.0/q_Exp^2 - 1)*(sin(theta)^2))) ;; Positional vector in polar coordinates, assuming y-direction inclination axis ratio, q EXPFAC = -1.67835D ;; Defines Ie_Exp and Re fexp = Ie_Exp * exp(ExpFac * (rvec-1.0)) ExpOut = 4.0 ;; Profile goes to zero at this radius (units of Re) ExpCut = 3.0*ExpOut/4.0 ;; Profile truncation begins at this radius (units of Re) ExpCut2 = (ExpCut+ExpOut)/2.0 wr_cutoff = where(rre GT ExpCut, nr_cutoff) wr_cutoff2 = where(rre GT ExpCut2, nr_cutoff2) wr_zero = where(rre GT ExpOut, nr_zero) fexpcut = fexp if nr_cutoff GT 0 then fexpcut[wr_cutoff] = fexpcut[wr_cutoff] * exp(-(rvec[wr_cutoff]-expcut)^2) if nr_cutoff2 GT 0 then fexpcut[wr_cutoff2] = fexpcut[wr_cutoff2] * exp(-(rvec[wr_cutoff2]-expcut2)^2) ;; if nr_cutoff GT 0 then fexpcut[wr_cutoff] = fexpcut[wr_cutoff] * ( 1 - ((rvec[wr_cutoff] - ExpCut)/(ExpOut - ExpCut))^2)^2.0 ;;if nr_cutoff GT 0 then fexpcut[wr_cutoff] = fexpcut[wr_cutoff] * ( 1 - ((rvec[wr_cutoff] - ExpCut)/(ExpOut - ExpCut))^2)^1.0 ;; NOTE: To reproduce phFitobj.h, last exponent must be 1.0 (not 2.0) ;;if nr_zero GT 0 then fexpcut[wr_zero] = 0 return, fexpcut end ;; Define integrand for the flux integral of the SDSS smoothly truncated exponential profile function fexpcut_radintgrnd, r, theta return, r * fexpcut(r, theta) end ;; SDSS deV: Smoothly truncated deVaucouleurs profile, 2D polar coordinates ;; R: Radius, physical units (same units as those defining Re) ;; theta: Polar angle, radians ;; function fDeVcut, r, theta common sdss_profs, Ie_Exp, q_Exp, Re_Exp, Ie_Dev, q_Dev, Re_Dev ;; Ie is surface brightness at Re if n_elements(q_Dev) EQ 0 then q_Dev = 1.0 q_Dev = float(q_Dev) rre = r / float(Re_Dev) rvec = rre * sqrt((1 + (1.0/q_Dev^2 - 1)*(sin(theta)^2))) ;; Positional vector in polar coordinates, assuming y-direction inclination axis ratio, q DEFAC = -7.66925D ;; Defines Ie_Dev and Re ;fdeV = Ie_DeV * exp(DEFAC *(rvec^0.25 - 1.0)) fdeV = Ie_DeV * exp(DEFAC *((rvec^2 + 0.0004)^0.125 - 1.0)) deVOut = 8.0D deVCut = 7*DEVOUT/8.0 wr_cutoff = where(rre GT deVCut, nr_cutoff) wr_zero = where(rre GT deVOut, nr_zero) fdeVcut = fdeV if nr_cutoff GT 0 then fdeVcut[wr_cutoff] = fdeVcut[wr_cutoff] * ( 1 - ((rvec[wr_cutoff] - deVCut)/(deVOut - deVCut))^2)^2.0 ;;if nr_cutoff GT 0 then fdeVcut[wr_cutoff] = fdeVcut[wr_cutoff] * ( 1 - ((rvec[wr_cutoff] - deVCut)/(deVOut - deVCut))^2)^1.0 ;; NOTE: To reproduce phFitobj.h, last exponent must be 1.0 (not 2.0) if nr_zero GT 0 then fdeVcut[wr_zero] = 0 return, fdeVcut end ;; SDSS deV with smoother truncation: Provides better MGE fits to SDSS-like profiles ;; R: Radius, physical units (same units as those defining Re) ;; theta: Polar angle, radians ;; function fDeVcut2, r, theta common sdss_profs, Ie_Exp, q_Exp, Re_Exp, Ie_Dev, q_Dev, Re_Dev ;; Ie is surface brightness at Re if n_elements(q_Dev) EQ 0 then q_Dev = 1.0 q_Dev = float(q_Dev) rre = r / float(Re_Dev) rvec = rre * sqrt((1 + (1.0/q_Dev^2 - 1)*(sin(theta)^2))) ;; Positional vector in polar coordinates, assuming y-direction inclination axis ratio, q DEFAC = -7.66925D ;; Defines Ie_Dev and Re ;fdeV = Ie_DeV * exp(DEFAC *(rvec^0.25 - 1.0)) fdeV = Ie_DeV * exp(DEFAC *((rvec^2 + 0.0004)^0.125 - 1.0)) deVOut = 8.0D deVCut = 7*DEVOUT/8.0 deVCut2 = (deVCut+deVOut)/2.0 wr_cutoff = where(rre GT deVCut, nr_cutoff) wr_cutoff2 = where(rre GT deVCut2, nr_cutoff2) wr_zero = where(rre GT deVOut, nr_zero) fdeVcut = fdeV if nr_cutoff GT 0 then fdeVcut[wr_cutoff] = fdeVcut[wr_cutoff] * exp(-(rvec[wr_cutoff]-devcut)^2) if nr_cutoff2 GT 0 then fdeVcut[wr_cutoff2] = fdeVcut[wr_cutoff2] * exp(-(rvec[wr_cutoff2]-devcut2)^2) ;;if nr_cutoff GT 0 then fdeVcut[wr_cutoff] = fdeVcut[wr_cutoff] * ( 1 - ((rvec[wr_cutoff] - deVCut)/(deVOut - deVCut))^2)^1.0 ;; NOTE: To reproduce phFitobj.h, last exponent must be 1.0 (not 2.0) ;;if nr_zero GT 0 then fdeVcut[wr_zero] = 0 return, fdeVcut end ;; Define integrand for the flux integral of the SDSS smoothly truncated exponential profile function fdeVcut_radintgrnd, r, theta return, r * fdeVcut(r, theta) end ;; ---------------------------------------------------------------------------------------------------------------------------------- ;; ---------------------------------------------------------------------------------------------------------------------------------- ;; Main Program ;; ---------------------------------------------------------------------------------------------------------------------------------- ;; pro sdss_to_aper6, DevMag, DevRad, DevAB, ExpMag, ExpRad, ExpAB, fracDeV, PSF_FWHM_array, aper_radius, mag_aperture, mag_aperture_Dev, mag_aperture_Exp, $ PIXEL_SIZE = PIXEL_SIZE, MODELMAG=MODELMAG, MOCK_IMAGE = MOCK_IMAGE, INSPECT = INSPECT, I0_Dev = I0_Dev, I0_Exp = I0_Exp, $ mge_exp = mge_Exp, mge_Dev = mge_Dev, mgeout = mgeout, hdr_mgeout = hdr_mgeout if not(keyword_set(not_SDSS)) then not_SDSS = 0 ;; Common block for 2D numerical integration of Sersic profile common sersic2D_common, I0, Re, n, bn, q ;; Common block for numerical integration of a gaussian profile common gauss2D_common, A0, sigma, sigma_y, qgauss, Rlim ;; Common block for SDSS truncated profiles common sdss_profs, Ie_Exp, q_Exp, Re_Exp, Ie_Dev, q_Dev, Re_Dev ;; Ie is surface brightness at Re nmags = n_elements(DevMag) mag_psf = fltarr(nmags) ;; ** Needs modification? *** if keyword_set(HIGHRES) then begin napers = 3 aper_radius = [1.0,1.0,1.0] endif else napers = n_elements(aper_radius) ;; **************************** ;; Set up PSF information PSF_FWHM_array_n_dim_orig = size(PSF_FWHM_array, /n_dim) if PSF_FWHM_array_n_dim_orig LE 1 then begin if n_elements(PSF_FWHM_array) EQ 1 then PSF_FWHM_array = fltarr(nmags) + PSF_FWHM_array n_PSF = 1 PSF_sigma = PSF_FWHM_array / 2.3548 ;; The gaussian sigma of the PSF R_PSF = 1./(2*!pi*PSF_sigma^2) ;; Amplitude of each PSF Gaussian term endif if PSF_FWHM_array_n_dim_orig EQ 2 then begin ;; This is a Multi-Gaussian PSF [2, n_PSF], but same for all sources. n_PSF = n_elements(PSF_FWHM_array[0,*]) ;; PSF_FWHM_full = fltarr(2, n_PSF, nmags) ;; for k=0,n_PSF-1 do begin ;; We need PSF_FWHM_full -> [2, n_PSF, nmags] ;; replicate_inplace, PSF_FWHM_full, PSF_FWHM_array[0,k], 3, [0,k,0] ;; replicate_inplace, PSF_FWHM_full, PSF_FWHM_array[1,k], 3, [1,k,0] ;; endfor PSF_FWHM_full = fltarr(nmags, 2, n_PSF) for k=0,n_PSF-1 do begin ;; We need PSF_FWHM_full -> [nmags, 2, n_PSF] replicate_inplace, PSF_FWHM_full, PSF_FWHM_array[0,k], 1, [0,0,k] replicate_inplace, PSF_FWHM_full, PSF_FWHM_array[1,k], 1, [0,1,k] endfor PSF_FWHM_array = PSF_FWHM_full ;; Replace with full array endif if size(PSF_FWHM_array, /n_dim) EQ 3 then begin ;; We now have a full Multi-Gaussian PSF array, defined for all sources (with dimensions, [nmags, 2, n_PSF] R_PSF = reform(PSF_FWHM_array[*,0,*]) ;; Amplitude of each PSF Gaussian term PSF_sigma = reform(PSF_FWHM_array[*,1,*])/2.3548 ;; The gaussian sigma of the PSF endif ;; edit if n_elements(aper_radius) GT 1 then mag_aperture = fltarr(nmags, napers)-99 else mag_aperture = fltarr(nmags)-99 mag_aperture_Exp = mag_aperture mag_aperture_Dev = mag_aperture magdiff_Exp = mag_aperture magdiff_Dev = mag_aperture I0_DeV_approx = fltarr(nmags) ;; Will contain central intensities, approximated analytically I0_Exp_approx = fltarr(nmags) ;; I0_DeV = fltarr(nmags)-1 ;; Will contain central intensities, computed from numerical integral I0_Exp = fltarr(nmags)-1 ;; wgood_Dev = where(DevMag GT 0 AND DevMag LT 100 AND DevRad GT 0, ngood_Dev) wgood_Exp = where(ExpMag GT 0 AND ExpMag LT 100 AND ExpRad GT 0, ngood_Exp) badflag_Dev = bytarr(nmags) + 1 ;; There's a problem if badflag=1 badflag_Exp = bytarr(nmags) + 1 badflag_Dev[wgood_Dev] = 0 badflag_Exp[wgood_Exp] = 0 ;; The following assumes regular magnitudes -- Luptitudes are different by ~1% at mags fainter than 22.3 ;dummy_zpt = 22.5 dummy_zpt = 30 ;; GAMA choice total_flux_Dev = 10^(0.4*(dummy_zpt-DevMag)) total_flux_Exp = 10^(0.4*(dummy_zpt-ExpMag)) ;bn_Dev = 2*4.0 - 0.324 ; This approximation works for 1 < n < 15 ;bn_Exp = 2*1.0 - 0.324 ; This approximation works for 1 < n < 15 n=4.0 & bn_Dev = newton((2*4.0 - 0.324)>0.1, 'sersic_bn') ;; Exact solution, interpolated n=1.0 & bn_Exp = newton((2*1.0 - 0.324)>0.1, 'sersic_bn') ;; Exact solution, interpolated ;; .............................................................................. ;; Determine central intensity normalization (I0) using an analytic approximation ;; NOTE: Not used in practice... Full 2D numerical integrals are performed instead ;; NOTE: The approximation gets worse as the inclination grows (toward low q) and the outer radius of integration decreases toward Re ;; Determine I0 for the DeV fit, n=4 ;R_fin_Dev = 8*DevRad ; The cutoff of the SDSS Dev model integration (sdss_to_aper4 version) R_fin_Dev = 8.0*DevRad ; The cutoff of the SDSS Dev model integration L_fin_Dev_norm = DevRad[wgood_Dev]^2 * 2*!pi*4.0 / (bn_Dev^(2*4.0)) * igamma(2*4.0, bn_Dev*(r_fin_Dev[wgood_Dev] / DevRad[wgood_Dev])^(1./4.0)) * gamma(2*4.0) * DevAB[wgood_Dev] ; Needed, to define I0_Dev I0_Dev_approx[wgood_Dev] = total_flux_Dev[wgood_Dev] / L_fin_Dev_norm ;; Determine I0 for the Exp fit, n=1 ;R_fin_Exp = 3.2*ExpRad ; The cutoff of the SDSS Exp model integration.. Note: The difference between 3.2 and 3.8 yields 0.02 mag offsets. (sdss_to_aper4 version) R_fin_Exp = 4.0*ExpRad ; The cutoff of the SDSS Exp model integration.. Note: The difference between 3.2 and 3.8 yields 0.02 mag offsets. L_fin_Exp_norm = ExpRad[wgood_Exp]^2 * 2*!pi*1.0 / (bn_Exp^(2*1.0)) * igamma(2*1.0, bn_Exp*(r_fin_Exp[wgood_Exp] / ExpRad[wgood_Exp])^(1./1.0)) * gamma(2*1.0) * ExpAB[wgood_Exp] ; Needed, to define I0_Exp I0_Exp_approx[wgood_Exp] = total_flux_Exp[wgood_Exp] / L_fin_Exp_norm ;; .............................................................................. ;; Set additional bad flags if approximate I0 values are bad wbad_I0_Dev_approx = where(I0_Dev_approx LE 0, nbad_I0_Dev_approx) if nbad_I0_Dev_approx GT 0 then badflag_Dev[wbad_I0_Dev_approx] = 1 wbad_I0_Exp_approx = where(I0_Exp_approx LE 0, nbad_I0_Exp_approx) if nbad_I0_Exp_approx GT 0 then badflag_Exp[wbad_I0_Exp_approx] = 1 ;; *** Needs modification?? *** if not(keyword_set(HIGHRES)) then begin ;; Normal, marginally resolved mode -- image size and aperture based on PSF and fixed if not(keyword_set(PIXEL_SIZE)) then pixel_size = 200.0 else pixel_size = float(pixel_size) dist_circle, im_circle, pixel_size pixscale = 0.4 ;pixscale = 6.0 / pixel_size ; arcsec/pix: Assume maximum diameter of interest is 6 arcsec ;pixscale = 20.0 / pixel_size ; arcsec/pix: Assume maximum diameter of interest is 6 arcsec makexy, -(pixel_size-1)/2, (pixel_size-1)/2, 1, -(pixel_size-1)/2, (pixel_size-1)/2, 1, xarray, yarray sin2_theta = (yarray / im_circle)^2 theta_array = asin(yarray / im_circle) ;;sin2_theta[(pixel_size-1)/2, (pixel_size-1)/2] = 0 endif ;; **************************************** dir_g09 = '/Users/kbundy/research/data/gama/g09/' psf_image_orig = readfits(dir_g09 + 'gama09_reconv_PSF.fits') ;; Uses 0.4 arcsec / pixel ;; npix_psf_orig = n_elements(psf_image_orig[0,*]) ;; psf_image = frebin(psf_image_orig, npix_psf_orig * 0.4/pixscale, /total) psf_image = psf_image_orig ;; Define MGE for Exponential; for sampling assume min(Re) ~ 10*Rmax, where Rmax is the maximum aperture radius of interest if n_elements(mge_exp) LE 1 then mge_exp = mgsersic(1.0, 100, NSTEP=300, RANGE = [1,1000], SDSS=~not_SDSS) ;; Define MGE for DeVaucouleurs; for sampling assume min(Re) ~ 10*Rmax, where Rmax is the maximum aperture radius of interest if n_elements(mge_Dev) LE 1 then mge_Dev = mgsersic(4.0, 100, NSTEP=300, RANGE = [1,1000], SDSS=~not_SDSS) ;; Structure for holding output convolved MGE information n_terms_exp = n_elements(mge_exp[0,*]) n_terms_dev = n_elements(mge_dev[0,*]) n_terms_max = max([n_terms_exp, n_terms_dev]) mgeout = {psf_gauss:fltarr(2, n_psf)-1, pref:-1, $ A0_exp:fltarr(n_terms_exp*n_psf)-99, sigma_x_exp:fltarr(n_terms_exp*n_psf)-99, sigma_y_exp:fltarr(n_terms_exp*n_psf)-99, $ A0_dev:fltarr(n_terms_dev*n_psf)-99, sigma_x_dev:fltarr(n_terms_dev*n_psf)-99, sigma_y_dev:fltarr(n_terms_dev*n_psf)-99} mgeout = replicate(mgeout, nmags) fxbhmake, hdr_mgeout, nmags, /initialize, /date sxaddpar, hdr_mgeout, 'LEN_UNIT', 'Arcsec', 'Units of sigma values' sxaddpar, hdr_mgeout, 'COEF', 'Amplitude', 'Type of MGE normalization coefficient' sxaddpar, hdr_mgeout, 'Pref', '1=Exp, 2=Dev', 'Type of MGE normalization coefficient' sxaddpar, hdr_mgeout, 'N_EXP', n_terms_exp, 'Number of terms in the Exponential MGE' sxaddpar, hdr_mgeout, 'N_DEV', n_terms_dev, 'Number of terms in the deVaucouleurs MGE' ;; Define image arrays for MGE flux summation ;; pixels_mge = 1000 ;; dist_circle, im_circle_mge, pixels_mge ;; makexy, -(pixels_mge-1)/2.0, (pixels_mge-1)/2.0, 1, -(pixels_mge-1)/2.0, (pixels_mge-1)/2.0, 1, xarray_mge, yarray_mge ;; sin2_theta_mge = (yarray_mge / im_circle_mge)^2 ;imgs_mge = fltarr(pixels_mge, pixels_mge, 20) ;; array to hold MGE realized in pixel form imgs_mge = fltarr(pixel_size, pixel_size, 20) ;; array to hold MGE realized in pixel form time_integrate = fltarr(nmags) ;; Time it takes to do MGE flux integrals time_mock = fltarr(nmags) ;; Time it takes to do mock image with filtering flux_exp_mge = fltarr(n_terms_exp, napers) ;; Integrated flux of each MGE term flux_exp_sum = fltarr(napers) ;; Total Flux flux_Dev_mge = fltarr(n_terms_Dev, napers) ;; Integrated flux of each MGE term flux_Dev_sum = fltarr(napers) ;; Total Flux A_invft_mge_exp = fltarr(n_terms_exp, n_PSF) sigma_array_exp = fltarr(n_terms_exp, n_PSF) sigma_y_array_exp =fltarr(n_terms_exp, n_PSF) A_invft_mge_dev = fltarr(n_terms_dev, n_PSF) sigma_array_dev = fltarr(n_terms_dev, n_PSF) sigma_y_array_dev =fltarr(n_terms_dev, n_PSF) ;; Loop over all objects for i=0L, nmags-1 do begin s1 = systime(/sec) if i GT 0 AND (i - 100*(i/100) EQ 0 OR i EQ nmags-1) then print, 'Number '+strc(i)+' of '+strc(nmags-1)+'. Time left: '+dplace(time_left,1)+' min' if badflag_Dev[i] GE 1 AND badflag_Exp[i] GE 1 then goto, skip ;; Set sdss_profs common block profile parameters ;;Ie_deV = I0 * exp(-bn*(1-0.0004^(1./8))) ;; Valid when I0=central surface brightness Ie_deV = 1.0 & n = 4.0 & bn = bn_dev Re_deV = DevRad[i] q_deV = DevAB[i] Ie_Exp = 1.0 & n = 1.0 & bn = bn_exp Re_Exp = ExpRad[i] q_Exp = ExpAB[i] ;; Determine central intensity (I0) for deV model if badflag_Dev[i] EQ 0 then begin ;; Numerical integration to determine exact total flux ;; I0 = 1.0 & n = 4.0 & bn = bn_dev ;; Re = DevRad[i] ;; q = DevAB[i] ;; flux_norm = INT_2D('sersic2D', [0,R_fin_Dev[i]], 'thetalims',48) Choices for resolution: 6, 10, 20, 48, or 96. flux_norm = INT_2D('fdevcut_radintgrnd', [0,R_fin_Dev[i]], 'thetalims',96, /double) ;; Without any corrections, predicted Dev mags are 0.06 mag brighter than Fibermags ;; This is likely due to how Dev profiles are truncated and the inner regions smoothed by SDSS. Since I can't exactly reproduce, ;; we introduce a small scalar offset. ;scale_offset_Dev = 10^(0.07/2.5) scale_offset_Dev = 1.0 I0_Dev[i] = total_flux_Dev[i] / flux_norm / scale_offset_Dev Ie_Dev = I0_Dev[i] ;I0_Dev[i] = I0_dev_approx[i] endif else im_Dev = im_circle*0 ;; Determine central intensity (I0) for Exp model if badflag_Exp[i] EQ 0 then begin ;; I0 = 1.0 & n = 1.0 & bn = bn_exp ;; Re = ExpRad[i] ;; q = ExpAB[i] ;; flux_norm = INT_2D('sersic2D', [0,R_fin_Exp[i]], 'thetalims',48) ; Choices for resolution: 6, 10, 20, 48, or 96. flux_norm = INT_2D('fExpcut_radintgrnd', [0,R_fin_Exp[i]], 'thetalims',48) ;; Without any corrections, predicted Exp mags are 0.075 mag brighter than Fibermags ;; This is likely due to how profiles are truncated smoothly to zero by SDSS. Since I can't exactly reproduce, ;; we introduce a small scalar offset. ;scale_offset_Exp = 10^(0.075/2.5) ;scale_offset_Exp = 10^(0.1/2.5) scale_offset_Exp = 1.0 I0_Exp[i] = total_flux_Exp[i] / flux_norm / scale_offset_Exp Ie_Exp = I0_Exp[i] ;I0_Exp[i] = I0_exp_approx[i] endif else im_Exp = im_circle*0 ;; ---------------------------------------------------------------------------------------------------- ;; Measure synthetic aperture photometry on a mock pixel-data image of the profile ;; NOTE: Used for checking only. It's much faster to use the MGE method instead use_Dev = 0 & use_Exp = 0 if keyword_set(MOCK_IMAGE) then begin t0_mock = systime(/sec) if keyword_set(MODELMAG) then begin absdiff_Dev = abs(MODELMAG[i] - DevMag[i]) absdiff_Exp = abs(MODELMAG[i] - ExpMag[i]) if absdiff_Dev LE 0.03 AND absdiff_Dev LE absdiff_Exp then begin I0_use = I0_dev[i] bn_use = bn_Dev Re_use = DevRad[i] q_use = DevAB[i] n_use = 4.0 Im = (pixscale)^2 * FdeVCut(im_circle*pixscale, theta_array) use_Dev = 1 endif if absdiff_Exp LE 0.03 AND absdiff_Exp LT absdiff_Dev then begin I0_use = I0_exp[i] bn_use = bn_Exp Re_use = ExpRad[i] q_use = ExpAB[i] n_use = 1.0 Im = (pixscale)^2 * FExpCut(im_circle*pixscale, theta_array) use_Exp = 1 endif endif else Re_use = max([DevRad[i], ExpRad[i]]) ;; *** Needs Modification?? *** ;; if keyword_set(HIGHRES) then begin ;; Resolved mode, Re >> sigma -- image size and aperture based on Re for each input source ;; pixscale = PSF_sigma[i] / 10.0 ;; pixel_size = 1.0*(round(2 * (2.0*Re_use) / pixscale) > round(2 * (4.0*PSF_sigma[i]) / pixscale)) ;; Full image size in pixels, 2 x 2-Re widths. NOTE: pixel_size must be > 3*FWHM or FILTER_IMAGE() crashes ;; aper_radius = ([0.5, 1.0, 1.5] * Re_use) < (pixel_size/2.0 - 1) ;; dist_circle, im_circle, pixel_size ;; makexy, -(pixel_size-1)/2, (pixel_size-1)/2, 1, -(pixel_size-1)/2, (pixel_size-1)/2, 1, xarray, yarray ;; sin2_theta = (yarray / im_circle)^2 ;; endif if not(keyword_set(MODELMAG)) then begin im_Dev = (pixscale)^2*I0_Dev[i] * exp( -bn_Dev*(im_circle*pixscale / DevRad[i])^0.25*(1 + (1/DevAB[i]^2-1)*sin2_theta)^(1./8) ) im_Exp = (pixscale)^2*I0_Exp[i] * exp( -bn_Exp*(im_circle*pixscale / ExpRad[i])*(1 + (1/ExpAB[i]^2-1)*sin2_theta)^(1./2) ) ;; Add relative Dev and Exp components together into a single model (similar to CMODELMAG?) im = fracDev[i]*im_Dev + (1-fracDev[i])*im_Exp endif else begin ;; if keyword_set(MODELMAG) then the simulated model is the one chosen for SDSS MODELMAG photometry ;; Note: The chosen model correlates with FracDev_R and LnLDev_R, but not perfectly because FracDev and LnLDev in the ;; other bands were incorporated into the choice for the ModelMag profile. So use ;; the difference between the supplied modelmag (.R) and either of the component mags to figure it out. ;im = (pixscale)^2*I0_use * exp( -bn_use*(im_circle*pixscale / Re_use)^(1.0/n_use)*(1 + (1.0/q_use^2-1)*sin2_theta)^(1.0/(2.0*n_use)) ) ;im = (pixscale)^2*I0_use * exp( -bn_use*((im_circle*pixscale / Re_use)^(1.0/n_use)*(1 + (1.0/q_use^2-1)*sin2_theta)^(1.0/(2.0*n_use)) -1) ) if abs(MODELMAG[i] - DevMag[i]) GT 0.03 AND abs(MODELMAG[i] - ExpMag[i]) GT 0.03 then begin ; In this case, not clear what MODELMAG is im_Dev = (pixscale)^2*I0_Dev[i] * exp( -bn_Dev*(im_circle*pixscale / DevRad[i])^0.25*(1 + (1/DevAB[i]^2-1)*sin2_theta)^(1./8) ) im_Exp = (pixscale)^2*I0_Exp[i] * exp( -bn_Exp*(im_circle*pixscale / ExpRad[i])*(1 + (1/ExpAB[i]^2-1)*sin2_theta)^(1./2) ) im = fracDev[i]*im_Dev + (1-fracDev[i])*im_Exp endif endelse ;; Convolve mock image with a PSF ;; We'll feed back the previous PSF array as along as the FWHM is identical and we're not in HIGHRES mode (in which case the array sizes change each time) ;; if i GT 1 then begin ;; if PSF_FWHM_array[i] EQ PSF_FWHM_array[i-1] AND not(keyword_set(HIGHRES)) then im_wPSF = filter_image(im, FWHM=PSF_FWHM_array[i]/pixscale, PSF=PSF) $ ;; else im_wPSF = filter_image(im, FWHM=PSF_FWHM_array[i]/pixscale) ;; endif else begin ;; im_wPSF = filter_image(im, FWHM=PSF_FWHM_array[i]/pixscale, PSF=PSF) ;; i=0, first source ;; endelse im_wPSF = filter_image(im, psf=PSF_image, fwhm=-1, /all) ;; Try PSF from GAMA reconvolution image ;im_wPSF = convolve(im, PSF_image) ;; im_wPSF = im ;; No PSF convolution, for testing ;;im_sampled = frebin(im_wPSF, 6/0.25, 6/0.25, /total) ;;aper, im_sampled, 6/0.25/2, 6/0.25/2, apermag, errap, sky, skyerr, 1.0, 1.5 / 0.25, 0, [0, 0], SETSKY=0 ;aper, im_wPSF, pixel_size/2.0, pixel_size/2.0, apermag, errap, sky, skyerr, 1.0, aper_radius / pixscale, 0, [0, 0], SETSKY=0, /silent ; Assumes zpt=25 aper, im_wPSF, (pixel_size-1)/2.0, (pixel_size-1)/2.0, apermag, errap, sky, skyerr, 1.0, aper_radius / pixscale, 0, [0, 0], SETSKY=0, /silent ; Assumes zpt=25 mag_aperture[i,*] = apermag - 25 + dummy_zpt ;; s = reverse(sort(im_wPSF)) ;; npix_core = (0.4 / pixscale / 1)^2 ;; psf_height = max(im_wPSF) ;; psf_height2 = mean(im_wPSF[s[0:npix_core-1]]) ; Considers height after "smoothing" PSF peak to something closer to the pixel scale of UKIDSS LAS (0.4 arcsec/pix) ;;print, psf_height/psf_height2 ;mag_psf[i] = -2.5*alog10(psf_height2 * 2*!PI*(PSF_FWHM_array[i]/pixscale/(2*sqrt(2*alog(2))))^2) + 25 ; PSF MAG: Using the input FWHM gaussian scaled to have a peak height equal to maximum of simulated source tend_mock = systime(/sec) time_mock[i] = tend_mock - t0_mock endif ;; ---------------------------------------------------------------------------------------------------- ;; Use the MGE expansion method to analytically describe the convolved profile ;if not(keyword_set(MOCK_IMAGE)) then begin ;; Next define implicit sampling scale (e.g., 1 Re = 100 spatial bins in the MGE) scale_sampling = 100 ;; 1 Re = 100 spatial bins mgeout[i].psf_gauss[0,*] = R_PSF[i,*] mgeout[i].psf_gauss[1,*] = psf_sigma[i,*] if 1 eq 1 then begin ;; Always do this (default) ;; .............................. ;; -- Exponential Profile, n=1 -- ellip_Exp = (1 - ExpAB[i]) ;; Ellipticity = (1-q) q_Exp = ExpAB[i] sig_mge_exp = mge_exp[1,*] * (ExpRad[i] / scale_sampling) ;; MGE gaussian sigma (arcsec) A_mge_exp = I0_Exp[i] * (mge_exp[0,*] / sqrt(2*!pi*mge_exp[1,*]^2)) ;; MGE gaussian coeffecients of the expansion ;; A_invft_mge_exp = (A_mge_exp * R_PSF[i] * 2*!pi * sig_mge_exp^2 * q_Exp * PSF_sigma^2) / $ ;; sqrt( (sig_mge_exp^2 + PSF_sigma^2)*(sig_mge_exp^2*q_Exp^2 + PSF_sigma^2) ) ;; Array of Prefactors (on convolved gaussian term) for k=0, n_PSF-1 do begin A_invft_mge_exp[*,k] = (A_mge_exp * R_PSF[i,k] * 2*!pi * sig_mge_exp^2 * q_Exp * psf_sigma[i,k]^2) / $ sqrt( (sig_mge_exp^2 + psf_sigma[i,k]^2)*(sig_mge_exp^2*q_Exp^2 + psf_sigma[i,k]^2) ) ;; Array of Prefactors (on convolved gaussian term) sigma_array_exp[*,k] = sqrt(sig_mge_exp^2 + psf_sigma[i,k]^2) sigma_y_array_exp[*,k] = sqrt(sig_mge_exp^2*q_Exp^2 + psf_sigma[i,k]^2) endfor qgauss = 0 ;; Set to zero so that the function gauss2D_1Dint will use the common block value of sigma_y t0_integrate = systime(/sec) ;; Integrate the summed convolved MGE numerically.... qgauss = 0 ;; Set to zero that the function gauss2D_1Dint will use the common block value of sigma_y A0 = A_invft_mge_exp[*] & A0_Exp = A0 sigma = sigma_array_exp[*] & sigma_Exp = sigma sigma_y = sigma_y_array_exp[*] & sigma_y_Exp = sigma_y ;; Store MGE info output mgeout[i].A0_exp = A0_exp mgeout[i].sigma_x_exp = sigma_Exp mgeout[i].sigma_y_exp = sigma_y_Exp ;; Integrate by pixels ......WARNING: WILL BREAK UNTIL PSF LOOP IS INCLUDED (5/4/12) ;; for j=0,ngauss_Exp-1 do imgs_mge[*,*,j] = A_invft_mge_exp[j] * exp( -(xarray_mge^2)/(2*(sig_mge_Exp[j]^2 + PSF_sigma^2)) - (yarray_mge^2)/(2*(sig_mge_Exp[j]^2*q^2 + PSF_sigma^2))) ;; invft_img_Exp = total(imgs_mge[*,*,0:ngauss_Exp-1], 3) ;; aper, invft_img_Exp, (pixels_mge-1)/2.0, (pixels_mge-1)/2.0, invft_apermag, errap, sky, skyerr, 1.0, aper_radius, 0, [0, 0], SETSKY=0, /silent ; Assumes zpt=25 ;; mag_aperture_Exp[i,*] = invft_apermag - 25 + dummy_zpt ;; Create summed mock convolved images using MGE ;; for j=0,ngauss_Exp-1 do imgs_mge[*,*,j] = (pixscale^2) * A0[j] * exp( -((pixscale*xarray)^2)/(2*(sigma[j]^2)) - ((pixscale*yarray)^2)/(2*(sigma_y[j]^2))) ;; invft_img_Exp = total(imgs_mge[*,*,0:ngauss_Exp-1], 3) ;; aper, invft_img_Exp, (pixel_size-1)/2.0, (pixel_size-1)/2.0, invft_apermag, errap, sky, skyerr, 1.0, aper_radius/pixscale, 0, [0, 0], SETSKY=0, /silent ; Assumes zpt=25 ;; mag_aperture2_exp = invft_apermag - 25 + dummy_zpt ;; .......................... for m=0, napers-1 do begin Rlim = aper_radius[m] ;; flux_Exp_MGE[j,m] = INT_2D('gauss2D', [0, Rlim], 'thetalims',96) ;; Roughly 5 times slower than QROMO flux_Exp_sum[m] = 2*QROMO('gauss2D_1dint', 0, Rlim, K=3, eps=1e-3) ;; Fastest is K=3 EPS=1e-3 ;;flux_Exp_MGE[j,m] = (2*QROMB('gauss2D_1dint', 0, Rlim, K=3, eps=1e-3)) ;; Note default is K=5 (higher precision) endfor mag_aperture_Exp[i,*] = dummy_zpt - 2.5*alog10(flux_Exp_Sum) ;; Integrate numerically ....WARNING: WILL BREAK UNTIL PSF LOOP IS INCLUDED (5/4/12) ;; for j=0, ngauss_exp-1 do begin ;; A0 = A_invft_mge_exp[j] ;; sigma = sqrt(sig_mge_exp[j]^2 + PSF_sigma^2) ;; sigma_y = sqrt(sig_mge_exp[j]^2*q_Exp^2 + PSF_sigma^2) ;; for m=0, napers-1 do begin ;; Rlim = aper_radius[m] ;; ;; flux_Exp_MGE[j,m] = INT_2D('gauss2D', [0, Rlim], 'thetalims',96) ;; Roughly 5 times slower than QROMO ;; flux_Exp_MGE[j,m] = 2*QROMO('gauss2D_1dint', 0, Rlim, K=3, eps=1e-3) ;; Fastest is K=3 EPS=1e-3 ;; ;;flux_Exp_MGE[j,m] = (2*QROMB('gauss2D_1dint', 0, Rlim, K=3, eps=1e-3)) ;; Note default is K=5 (higher precision) ;; endfor ;; endfor ;; mag_aperture_Exp[i,*] = dummy_zpt - 2.5*alog10(total(flux_Exp_MGE, 1)) ;; ............................ tend_integrate = systime(/sec) time_integrate[i] = tend_integrate - t0_integrate ;; ................................ ;; -- Devaucouleurs Profile, n=4 -- ellip_Dev = (1 - DevAB[i]) ; Ellipticity = (1-q) q_Dev = DevAB[i] sig_mge_Dev = mge_Dev[1,*] * (DevRad[i] / scale_sampling) ;; MGE gaussian sigma (arcsec) A_mge_Dev = I0_Dev[i] * (mge_Dev[0,*] / sqrt(2*!pi*mge_Dev[1,*]^2)) ;; MGE gaussian coeffecients of the expansion for k=0, n_PSF-1 do begin A_invft_mge_Dev[*,k] = (A_mge_Dev * R_PSF[i,k] * 2*!pi * sig_mge_Dev^2 * q_Dev * psf_sigma[i,k]^2) / $ sqrt( (sig_mge_Dev^2 + psf_sigma[i,k]^2)*(sig_mge_Dev^2*q_Dev^2 + psf_sigma[i,k]^2) ) ;; Array of Prefactors (on convolved gaussian term) sigma_array_dev[*,k] = sqrt(sig_mge_Dev^2 + psf_sigma[i,k]^2) sigma_y_array_dev[*,k] = sqrt(sig_mge_Dev^2*q_Dev^2 + psf_sigma[i,k]^2) endfor ;; Integrate the summed convolved MGE numerically.... qgauss = 0 ;; Set to zero that the functino gauss2D_1Dint will use the common block value of sigma_y A0 = A_invft_mge_Dev[*] & A0_Dev = A0 sigma = sigma_array_dev[*] & sigma_Dev = sigma sigma_y = sigma_y_array_dev[*] & sigma_y_Dev = sigma_y ;; Store MGE info output mgeout[i].A0_dev = A0_dev mgeout[i].sigma_x_dev = sigma_Dev mgeout[i].sigma_y_dev = sigma_y_Dev ;; Create summed mock convolved images using MGEWARNING: WILL BREAK UNTIL PSF LOOP IS INCLUDED (5/4/12) ;; for j=0,ngauss_Dev-1 do imgs_mge[*,*,j] = (pixscale^2) * A0[j] * exp( -((pixscale*xarray)^2)/(2*(sigma[j]^2)) - ((pixscale*yarray)^2)/(2*(sigma_y[j]^2))) ;; invft_img_Dev = total(imgs_mge[*,*,0:ngauss_Dev-1], 3) ;; aper, invft_img_Dev, (pixel_size-1)/2.0, (pixel_size-1)/2.0, invft_apermag, errap, sky, skyerr, 1.0, aper_radius/pixscale, 0, [0, 0], SETSKY=0, /silent ; Assumes zpt=25 ;; mag_aperture2_dev = invft_apermag - 25 + dummy_zpt ;; .......................... for m=0, napers-1 do begin Rlim = aper_radius[m] ;; flux_Dev_MGE[j,m] = INT_2D('gauss2D', [0, Rlim], 'thetalims',96) ;; Roughly 5 times slower than QROMO flux_Dev_sum[m] = 2*QROMO('gauss2D_1dint', 0, Rlim, K=3, eps=1e-3) ;; Fastest is K=3 EPS=1e-3 ;;flux_Dev_MGE[j,m] = (2*QROMB('gauss2D_1dint', 0, Rlim, K=3, eps=1e-3)) ;; Note default is K=5 (higher precision) endfor ;; for j=0, ngauss_Dev-1 do beginWARNING: WILL BREAK UNTIL PSF LOOP IS INCLUDED (5/4/12) ;; A0 = A_invft_mge_Dev ;; sigma = sqrt(sig_mge_Dev^2 + PSF_sigma^2) ;; sigma_y = sqrt(sig_mge_Dev^2*q_Dev^2 + PSF_sigma^2) ;; for m=0, napers-1 do begin ;; Rlim = aper_radius[m] ;; flux_Dev_Sum[m] = 2*QROMO('gauss2D_1dint', 0, Rlim, K=3, eps=1e-3) ;; Appears to be 5 times faster than QROMB K=3 ;; ;;flux_Dev_MGE[j,m] = (2*QROMB('gauss2D_1dint', 0, Rlim, K=3)) ;; Note default is K=5 (higher precision) ;; endfor ;; endfor mag_aperture_Dev[i,*] = dummy_zpt - 2.5*alog10(flux_Dev_Sum) ;; forprint, mag_aperture[i,*], mag_aperture_Dev[i,*], mag_aperture2_Dev & print & print ;; plot, im_wPSF[(pixel_size-1)/2.0,*] ;; oplot, invft_img_Dev[(pixel_size-1)/2.0,*], col=!red ;; noise_sigma = 0.5*Ie_Dev * pixscale^2 ;; imobs_wnoise = invft_img_Dev + randomn(seed, pixel_size, pixel_size)*noise_sigma ;; display2, imobs_wnoise ;; ;aper, imobs_wnoise, (pixel_size-1)/2.0, (pixel_size-1)/2.0, wnoise_apermag, errap, sky, skyerr, 1.0, aper_radius/pixscale, [0.85,0.99]*(pixel_size-1)/2.0, [0, 0] ; Assumes zpt=25 ;; skyrad = 4*[1,1.5]/pixscale ;; aper, imobs_wnoise, (pixel_size-1)/2.0, (pixel_size-1)/2.0, wnoise_apermag, errap, sky, skyerr, 1.0, aper_radius/pixscale, skyrad, [0, 0] ; Assumes zpt=25 ;; tvcircle, skyrad[0], (pixel_size-1)/2.0, (pixel_size-1)/2.0, /data, col=!red ;; mag_wnoise = wnoise_apermag - 25 + dummy_zpt ;; forprint, mag_aperture_Dev[i,*], mag_wnoise, mag_wnoise - mag_aperture_Dev[i,*] ;; Store either the Dev of Exp synmags in the generic output synmag variable if keyword_set(MODELMAG) AND not(keyword_set(MOCK_IMAGE)) then begin absdiff_Dev = abs(MODELMAG[i] - DevMag[i]) absdiff_Exp = abs(MODELMAG[i] - ExpMag[i]) if absdiff_Exp LE 0.03 AND absdiff_Exp LE absdiff_Dev then begin mag_aperture[i,*] = mag_aperture_Exp[i,*] mgeout[i].pref = 1 ;; Set preferred value to 1 = Exp endif if absdiff_Dev LE 0.03 AND absdiff_Dev LT absdiff_Exp then begin mag_aperture[i,*] = mag_aperture_Dev[i,*] mgeout[i].pref = 2 ;; Set preferred value to 2 = deV endif endif endif skip: s2 = systime(/sec) time_left = (s2-s1)*(nmags-1-i)/60 ; Minutes endfor end