UNDER DEVELOPMENT

The software package is currently being dramatically improved and updated.
This is a (very) old version, potentially with bugs.
  Shapelets web page     Shapelets IDL code     Installation     Help pages     Simulated images  
List of routines
Instrutions for general code
Instructions for image simulation code

Source code for simage_make_analytic_object.pro:

You can also view the help page for this routine.

function simage_make_analytic_object, morphology, rsq, flux, position, $ SIGMAE=sigmae, $ SEED=seed, $ GAMMA_INPUT=gamma_input, $ PSFIMAGE=psfimage, $ INST_SHEAR=inst_shear, $ SHAPEINFO=shapeinfo ;$Id: simage_make_analytic_object.pro, v1.0$ ; ; Copyright © 2004 Richard Massey and Alexandre Refregier. ; ; This file is a part of the Shapelets analysis code. ; www.ast.cam.ac.uk/~rjm/shapelets/ ; ; The Shapelets code is free software; you can redistribute it and/or ; modify it under the terms of the GNU General Public License as published ; by the Free Software Foundation; either version 2 of the License, or ; (at your option) any later version. ; ; The Shapelets code is distributed in the hope that it will be useful, ; but WITHOUT ANY WARRANTY; without even the implied warranty of ; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ; GNU General Public License for more details. ; ; You should have received a copy of the GNU General Public License ; along with the Shapelets code; if not, write to the Free Software ; Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA ; ;+ ; NAME: ; SIMAGE_MAKE_ANALYTIC_OBJECT ; ; CATEGORY: ; Shapelet simulated images. ; ; PURPOSE: ; This module manufactures one object for the image simulations, using a ; radially symmetric, analytic profile (such as de Vaucouleurs or exponential ; disc) for its morphology. The size, magnitude and ellipticity are drawn from ; a real galaxy template - or randomised from a Gaussian distribution of ; ellipticities with width sigmae. The object is pixellated using various input ; parameters. ; ; LIMITATION: ; Unable to add astrometric distortions to these analytic profiles. ; ; INPUTS: ; decomp - Cartesian shapelet decomp structure of a template real ; galaxy. ; morphology - Which parameterised galaxy profile? ; 1: de Vaucouleurs e^(-r^(1/4)) ; 2: Exponential disc e^(-r) ; 3: Gaussian e^(-r^2/2) ; 4: Moffat ; ; OPTIONAL INPUTS: ; sigmae - rms galaxy ellipticity (before shear) with parameterised ; galaxy shapes. ; gamma - Weak lensing shear on all galaxies [gamma1,gamma2]. ; psf - Cartesian decomp structure representing PSF. Will then do ; convolvution automatically. Stars are just this PSF. ; seed - Seed for random number generator. ; ; OUTPUTS: ; recomp - One object (with an analytic profile) in a postage stamp image. ; shapeinfo - A structure containing object properties for inclusion in the ; pre-noise catalogue. ; ; MODIFICATION HISTORY: ; Jun 03 - Written by Richard Massey ;- if keyword_set(inst_shear) then message,'Cannot add instrumental distortions to analytic profiles' if morphology eq 1 or morphology eq 4 then message,'Moffat and de Vaucouleurs profiles not yet implemnented' ; ; Determine object properties in terms of size and ellipticity ; if keyword_set(sigmae) then begin ; Set ellipticity from a Gaussian PDF ell=complex(randomn(seed),randomn(seed))*sigmae/sqrt(2.) flux=flux>1e-20 ; Catching errors from bad decompositions endif else begin ; Or get this also from original (real) galaxy distribution q=shapelets_quadrupole(decomp,flux=flux) ; flux=flux>1e-20 ; Catching errors from bad decompositions rsq=(Q[0,0]+Q[1,1])/flux ; ell=complex(Q[0,0]-Q[1,1],2*Q[0,1])/(Q[0,0]+Q[1,1]) ell=ell*exp(complex(0,1)*randomu(seed)*2*!pi) ; Randomise PA endelse ; ; ; Shear object at catalogue level. ; if keyword_set(gamma_input) then begin ; ; Shapelets II: eqn. 25 wrongly says pgamma=2.+abs(ell) ;e1=float(ell) +2*gamma_input[0]+ 2*float(ell)*(float(ell)*gamma_input[0]+imaginary(ell)*gamma_input[1]) ;e2=imaginary(ell)+2*gamma_input[1]+2*imaginary(ell)*(float(ell)*gamma_input[0]+imaginary(ell)*gamma_input[1]) ;ell=complex(e1,e2) ; ; but Shapelets IV: disagrees e1 =float(ell) +2*gamma_input[0]-2*float(ell)* (float(ell)*gamma_input[0]+imaginary(ell)*gamma_input[1]) e2 =imaginary(ell)+2*gamma_input[1]-2*imaginary(ell)*(float(ell)*gamma_input[0]+imaginary(ell)*gamma_input[1]) ell=complex(e1,e2) rsq=rsq*(1+2*(float(ell)*gamma_input[0]+imaginary(ell)*gamma_input[1])) endif ; ; Determine object's new properties in terms of its major and minor axes. ; PA=atan(imaginary(ell),float(ell))/2. A=sqrt(rsq*(1+abs(ell))/2.) ; R^2 = A^2 + B^2 B=sqrt(abs(rsq*(1-abs(ell))/2.)) ; e = (A^2 - B^2) / (A^2 + B^2) if ( (not finite(A)) or (not finite(B)) or (not finite(PA)) ) then return,-1 PA=PA-!pi/2. ; ; Create a blank postage stamp image and co-ordinate grid. ; case morphology of ; Size postage stamp to include wings. 1: n_pix=30*A>10 ; de Vaucouleurs profile e^(-r^(1/4)) has rsq_f=10461394944000*!pi/F=259459200 2: n_pix=20*A>10 ; Exponential disc profile e^(-r) has rsq_f=12.*!pi/F=6 3: n_pix=10*A>10 ; Gaussian profile e^(-r^2/2) has rsq_f=4.*!pi/F endcase n_pix=replicate(2*fix(n_pix)+1,2) ; Make sure postage stamp has odd dimensions. position=position+fix(n_pix/2.) ; Include sub-pixel offset. mk_xarr,n_pix,x1,x2,x0=position ; Create (x,y) co-ordinate array of postage stamp pixels. recomp=fltarr(n_pix[0],n_pix[1]) ; Create a blank new postage stamp image. ; ; Distort the co-ordinate grid so that the created object has the desired ; axis ratio etc (see Shapelets: IV) ; matrixR = [[cos(PA),sin(PA)],[-1*sin(PA),cos(PA)]] matrixJ = transpose(matrixR)##[[A^2,0],[0,B^2]]##matrixR / rsq matrixJi= invert(matrixJ) detJ = determ(matrixJ) rr=x1 for ii=0,n_pix[0]-1 do begin for jj=0,n_pix[1]-1 do begin x=[x1[ii,jj],x2[ii,jj]] rr[ii,jj]=transpose(x)#matrixJ#x /rsq endfor endfor ; ; Find maximum r within our postage stamp (using only this limited ; dynamical range and truncating *elliptical* isophotes at this value will ; prevent edge effects due to the *square* postage stamp). ; rmax=min([transpose(rr[0,*]),transpose(rr[n_pix[0]-1,*]),rr[*,0],rr[*,n_pix[1]-1]]) case morphology of 1: rmax=rmax<400 ; de Vaucouleurs profile e^(-r^(1/4)) has HUGE wings 2: rmax=rmax<100 ; Exponential disc profile e^(-r) has large wings else: rmax=rmax<25 ; Gaussian profile e^(-r^2/2) has small wings endcase good=where((rr-rmax)<0) if good[0] eq -1 then return,-1 ; Error catching (tiny object) ; ; Draw galaxy with analytic profile. ; case morphology of ; Draw galaxy with de Vaucouleurs e^r^(1/4) profile 1: recomp[good]=flux/sqrt(detJ)*6435/!pi*exp(-11.26571173*(rr[good])^(0.125)) /rsq ; Draw galaxy with exponential disc (3/pi)*e^(-sqrt6*r) profile 2: recomp[good]=flux/sqrt(detJ)*3/!pi*exp(-1.*sqrt(6*rr[good])) /rsq ; Draw elliptical Gaussian blob (1/pi)*e^(-r^2) has F=r_f^2=1 else: recomp[good]=flux/sqrt(detJ)/!pi*exp(-1.*rr[good])/!pi /rsq endcase ; ; Renormalise to correct flux ; recomp=recomp/4 ; Not entirely sure why, but it works! recomp=recomp/total(recomp)*flux ; NB: ought not perform the previous line because the integration then depends ; on pixel size. However, without a proper (analytic) integration within ; pixels, it is going to anyway. Just using central pixel values is like ; Simpson`s rule (only no gradients!). ; PS: COULD DIVIDE BY BETA^2 (which is the pixel area)! ; ; Smear object to simulate seeing. ; (Add zero padding, smooth, then remove padding: slow but effective!) ; if keyword_set(psfimage) then begin ; Could speed up by shrinking postage stamp as much as poss first! smoothborder=((size(psfimage))[1:2]+1)/2 recomps=fltarr(n_pix[0]+2*smoothborder[0],n_pix[1]+2*smoothborder[1]) recomps[smoothborder[0]:n_pix[0]+smoothborder[0]-1,smoothborder[1]:n_pix[1]+smoothborder[1]-1]=recomp recomp=convol(recomps,psfimage,/edge_truncate) ; Reduce size of new postage stamp back to original ;recomp=recomp[smoothborder[0]:decomp.n_pixels[0]+smoothborder[0]-1,smoothborder[1]:decomp.n_pixels[1]+smoothborder[1]-1] ; Or just offset its centre position=position+smoothborder endif ; ; Record shape information for posterity ; shapeinfo={A:A,B:B,theta:(PA*180/!pi+90),rsq:rsq,flux:flux,mag:-2.5*alog10(flux)+22.08} ; That's all folks! return,recomp end

View the help page for this routine, return to the shapelets web page or return to the code help menu.


Last modified on 02nd May 2008 by Richard Massey.

Valid HTML 4.01!