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_assemble_image.pro:

You can also view the help page for this routine.

;$Id: simage_assemble_image.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_LOCATE_OBJECT ; ; CATEGORY: ; Shapelets image simulation. ; ; PURPOSE: ; This module takes a postage-stamp object (galaxy or star) and where its centre ; should go, then checks that it doesn`t fall off the edge of an image with ; n_pix pixels. If it partly falls off the edge, it is guillotined to fit. ; the coordinates of the postage stamp are returned in the region structure. ; If the object lies entirely off the edge, region is returned as -1. ; ; INPUTS: ; recomp - Postage stamp image of an object. ; recomp_x0 - Position of object centre in postage stamp coordinates. [x,y] ; position - Desired position of object centre in final image. [x,y] ; image_pix - Number of pixels in final image. [x,y] ; ; OPTIONAL INPUTS: ; none. ; ; KEYWORD PARAMETERS: ; none. ; ; OUTPUTS: ; recomp - Trimmed as necessary. ; region - A structure containing coordinates of left, right, top and ; bottom corners of where the object should lie in the final ; image. If object lies outside this, all of the values in ; region are flagged to -1. ; ; MODIFICATION HISTORY: ; Mar 04 - Tidied up by RM. ; Jun 03 - Written by Richard Massey pro simage_locate_object, recomp, recomp_x0, position, image_pix, region COMPILE_OPT idl2 ; Determine initial size of postage stamp postage_pix = size(recomp,/dimensions) ; Where you would hope to locate postage stamp in image if away from edges l=position[0]-recomp_x0[0] r=l+postage_pix[0]-1 b=position[1]-recomp_x0[1] t=b+postage_pix[1]-1 ; Declare region structure region={l:round(l), r:round(r), b:round(b), t:round(t)} ; Check to see if object overlaps with any edges if region.l lt 0 then begin ; Overlaps left of image if region.r le 0 then goto,outside ; Falls off left of image recomp=recomp[-1*region.l:postage_pix[0]-1,*] ; postage_pix[0]=postage_pix[0]+region.l ; recomp_x0[0]=recomp_x0[0]+region.l ; region.l=0 ; endif ; if region.r gt image_pix[0]-1 then begin ; Overlaps right of image if region.l ge image_pix[0]-1 then goto,outside ; Falls off right of image recomp=recomp[0:postage_pix[0]-region.r+image_pix[0]-2,*] postage_pix[0]=postage_pix[0]-region.r+image_pix[0]-1 region.r=image_pix[0]-1 ; endif ; if region.b lt 0 then begin ; Overlaps bottom of image if region.t le 0 then goto,outside ; Falls off bottom of image recomp=recomp[*,-1*region.b:postage_pix[1]-1] ; postage_pix[1]=postage_pix[1]+region.b ; recomp_x0[1]=recomp_x0[1]+region.b ; region.b=0 ; endif ; if region.t gt image_pix[1]-1 then begin ; Overlaps top of image if region.b ge image_pix[1]-1 then goto,outside ; Falls off top of image recomp=recomp[*,0:postage_pix[1]-region.t+image_pix[1]-2] postage_pix[1]=postage_pix[1]-region.t+image_pix[1]-1 region.t=image_pix[1]-1 ; endif ; return outside: region={l:-1,r:-1,t:-1,b:-1} return end ; ************************************************************************ ; ************************************************************************ ; ; ; NAME: ; SIMAGE_STAR_MAGNITUDE_DISTRIBUTION ; ; CATEGORY: ; Shapelets image simulation. ; ; PURPOSE: ; Generates a random flux value for a star, drawn from the stellar ; luminosity function in Allen's astrophysical quantities (extrapolated ; to the necessary depth). ; ; INPUTS: ; none. ; ; OPTIONAL INPUTS: ; seed - Random number seed. ; slf - Stellar luminosity function. Caluclated the first time this ; routine is called, then simply referenced on subsequent occasions ; for speed. ; ; KEYWORD PARAMETERS: ; none. ; ; OUTPUTS: ; mag - Magnitude of randomly sampled star. ; ; MODIFICATION HISTORY: ; Jul 03 - Written by Richard Massey function simage_star_magnitude_distribution,seed=seed,slf=slf if not keyword_set(slf) then begin ; Data from Allen is log(stars per mag per sq.degree) at b=90 imag=[12,13,14,15,16,17,18,19,20,21,22,23,24,25] slf=[1.374,1.728,2.045,2.326,2.575,2.799,3.005,3.191,3.353,3.488,3.604,3.706,3.797,3.875] ; Extrapolate for i=1,5 do begin imag=[imag,imag[13]+i] slf=[slf,slf[13]+0.06*i] endfor ; Convert to real numbers slf=10^slf ; Convert to cumulative PDF for i=1,n_elements(slf)-1 do begin slf[i]=slf[i]+slf[i-1] endfor ; Renormalise PDF (and destroy stellar density information) message,strtrim(slf[n_elements(slf)-1],2)+" stars expected per square degree at galactic poles",/info slf=slf/slf[n_elements(slf)-1] slf=[0.,slf] imag=[11,imag] ; Store as a convenient structure slf={slf:slf,imag:imag} endif ; Select a random integer magnitude from this PDF mag=min(slf.imag[where(randomu(seed) le slf.slf)]) ; Now interpolate PDF as being flat within each integer mag=float(mag)+randomu(seed)-0.5 ; Return return,mag end ; ************************************************************************ ; ************************************************************************ ; ;+ ; NAME: ; SIM_ASSEMBLE ; ; CATEGORY: ; Shapelets image simulation. ; ; PURPOSE: ; Collects newly synthesised objects from sim_create.pro, rotates them ; by a random amount, and scatters them into pixels to make a simulated ; image. Background levels and noise are then be added using ; simage_add_noise.pro. ; ; This routine can also create a simulated image containing objects with ; analytic (e.g. de Vaucouleurs) profiles with the same size, magnitude ; and ellipticity distributions as a shapelet catalogue (e.g. for the HDF). ; The shapelet coefficients are not used for any purpose other than the ; determination of these distributions (which may also be entered by hand). ; The point is to use a non-shapelet method to create simulations, which ; can themselves be then used to test shapelet-based shape-measurement ; techniques. ; ; INPUTS: ; cat - Synthesised catalogue of shapelet objects, eg created by new_hdf ; file - Filename for output w/o extension, eg "newimage" ; ; OPTIONAL INPUTS: ; [seed] - Seed for random number generator. ; [noise] - Exposure time to pass to simage_add_noise.pro ; [gamma] - Weak lensing shear on all galaxies [gamma_1,gamma_2]. ; [psf] - Cartesian decomp structure representing PSF. Will then do ; convolvution automatically. Stars are just this PSF. ; [inst] - Post-smear shear due to geometric distortions in telescope. ; [pixsize] - Pixel scale [arcsec]. Default is 0.04". ; [n_pix] - Number of pixels in final image [x,y]. ; [n_gal] - Number of galaxies in final image. ; [n_star] - Number of stars in final image (currently all same magnitude). ; [dither] - Dither offset of exposure [x,y]. (Currently a simple linear shift ; which is the same everywhere and should be a fraction of a pixel). ; [parameter] - Use parameterised models of galaxy shapes rather than shapelet ; coefficients. ; 1: de Vaucouleurs r^(1/4) profiles ; 2: Exponential discs ; 3: Elliptical Gaussian blobs ; [sigmae] - rms galaxy ellipticity (before shear) for use with ; parameterised galaxy shapes. ; ; KEYWORD PARAMETERS: ; /RPOS - (Re-)randomise positions and PAs of galaxies? 1=YES, 0=NO ; ; OUTPUTS: ; Image, written to disc. ; This image is also (optionally) returned to the variable "image", with ; its FITS header in the variable "header". ; ; MODIFICATION HISTORY: ; Oct 2003 - Centroid bug fixing by RM. ; Jul 2003 - Input file in SExtractor format added by RM. ; Jun 2003 - Broken into separate modules by RM. ; Apr 2003 - Bugs fixed and analytic profiles forms added by RM. ; Feb 2002 - Written by Richard Massey. ;- pro simage_assemble_image, shapecat, file, seed=seed, gamma_input=gamma_input, psf=psf,$ rpos=rpos, n_pix=n_pix, n_gal=n_gal, n_star=n_star, $ pixsize=pixsize, parameter=parameter, sigmae=sigmae, $ image=image, header=header, noise=noise, dither=dither ; ; ; SET DEFAULTS ; ; if not keyword_set(shapecat)then message,'Requires an input shapelet catalogue!' if not keyword_set(file) then file = 'blah' ; Output file name if not keyword_set(seed) then seed = 101L ; Initialise random seed (if set automatically by IDL, the process is unrepeatable). if not keyword_set(sigmae) then sigmae = 0.3 ; if not keyword_set(pixsize) then pixsize= 0.04 ; [arcsec] if not keyword_set(n_pix) then n_pix = fix([4096,4096]*0.04/pixsize) ; To make roughly HDF size by default if n_elements(n_gal) eq 0 then n_gal = round( 6000. * n_pix[0]/4096. * n_pix[1]/4096. * (pixsize^2)/(0.04^2) ) if n_elements(n_star) eq 0 then n_star = 186 ; Number of stars randomly dispersed in image if not keyword_set(psf) then n_star = 0 ; Can't include stars without knowing the PSF! ;if not keyword_set(rpos) then n_gal=fix(n_gal*9./4.) ; Adjust for larger area visualised around the image for dithers if shapecat.n lt n_gal then begin ; ^-- user input (fudge factor) message,'Warning: not enough objects in catalogue! '+ $ 'Randomising object positions and cycling through catalogue.',/info rpos=1 endif ; ; ; PARSE REMAINING INPUT PARAMETERS AND SET UP LOCAL VARIABLES ; ; ; Read catalogue if not already in memory. if n_elements(shapecat) eq 0 then shapelets_read_shapecat,shapecat,file,/moments ; Ensure input shear is in correct notation. If no shear, unset variable. if keyword_set(gamma_input) then begin if (size(gamma_input))[0] eq 0 then gamma_input=[gamma_input,0.] if gamma_input[0] eq 0 and gamma_input[1] eq 0 then gamma_input=0 endif ; Ensure PSF is in correct format to be useful later on. if keyword_set(psf) then begin message,'Convolving with a PSF to n_max='+strtrim(psf.n_max,2)+'.',/info if keyword_set(parameter) then begin ; Need an image to convolve analytic profiles. shapelets_recomp,psf,psfimage n_gal=fix(n_gal*(1.+float(psf.nf[0])/n_pix[0]*psf.nf[1]/n_pix[1])) endif ; Work out PSF shape properties for pre-noise catalogue (assuming PSF constant everywhere). q=shapelets_quadrupole(psf,flux=flux) 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]) PA=atan(imaginary(ell),float(ell))/2.+!pi/2 A=sqrt(rsq*(1+abs(ell))/2.) B=sqrt(abs(rsq*(1-abs(ell))/2.)) psfinfo={A:A,B:B,theta:(PA*180/!pi) mod 360-90,rsq:rsq,flux:flux,mag:-2.5*alog10(flux)+22.08} endif ; ; ; PRINT INTENT TO SCREEN AS A RECORD ; ; message,'Placing '+strtrim(n_gal,1)+' galaxies and '+strtrim(n_star,1)+ $ ' stars into a '+strtrim(n_pix(0),1)+'x'+strtrim(n_pix(1),1)+ $ ' grid.',/info message,'Pixel scale is '+strtrim(pixsize,1)+' arcsec/pix.',/info message,'Random seed='+strtrim(string(seed[0]),1)+'.',/info if keyword_set(rpos) or ((shapecat.x(0,0) eq 0.) and (shapecat.x(0,1) eq 0.)) then $ message,'Object positions are randomised.',/info if keyword_set(gamma_input) then message,'Input shear is '+strtrim(gamma_input[0],2)+','+strtrim(gamma_input[1],2)+'.',/info if keyword_set(inst) then message,'Instrumental distortions are applied.',/info ;message,'Writing simulated image to file '+get_path(3,/silent)+file+'.fits',/info ; Begin mock SExtractor catalogue of pre-noise shapes. openw,lun,shapelets_paths(3)+file+'.sex',/get_lun printf,lun,'# Mock SExtractor catalogue of simulated image created before the addition of noise.' printf,lun,'# For more information, see Massey et al. 2003 or email Richard Massey on rjm@ast.cam.ac.uk' printf,lun,'#' printf,lun,'# 1 X_IMAGE Object position along x [pixel]' printf,lun,'# 2 Y_IMAGE Object position along y [pixel]' printf,lun,'# 3 BLANK Spacer column to match SExtractor' printf,lun,'# 4 BLANK Spacer column to match SExtractor' printf,lun,'# 5 XMIN_IMAGE Minimum x-coordinate of postage stamp [pixel]' printf,lun,'# 6 YMIN_IMAGE Minimum y-coordinate of postage stamp [pixel]' printf,lun,'# 7 XMAX_IMAGE Maximum x-coordinate of postage stamp [pixel]' printf,lun,'# 8 YMAX_IMAGE Maximum y-coordinate of postage stamp [pixel]' printf,lun,'# 9 A_IMAGE Profile along major axis { e=(A^2 - B^2) } [pixel]' printf,lun,'# 10 B_IMAGE Profile along minor axis { (A^2 + B^2) } [pixel]' printf,lun,'# 11 THETA_IMAGE Position angle (CCW/x) [deg]' printf,lun,'# 12 CLASS_STAR S/G classifier output' printf,lun,'# 13 SHAPELET_SIZE R=\sqrt{\int_0^\infinity{ f(r) r^2 rdr}} [pixel]' printf,lun,'# 14 MAGNITUDE Zeropoint not normalised [mag]' printf,lun,'# 15 FLUX Not sure of units here' printf,lun,'# 16 BLANK Spacer column to match SExtractor' printf,lun,'# 17 NUMBER Running object number' mocksexformat='(F10.3,1X,F10.3,A23,I10,1X,I10,1X,I10,1X,I10,1X,'+$ 'F9.3,1X,F9.3,1X,F5.1,1X,F5.2,1X,F8.2,1X,F8.4,1X,F12.5,A11,I10)' ; ; ; CREATE IMAGE ARRAY ; ; image=fltarr(n_pix[0],n_pix[1]) ; ; ; POPULATE IMAGE WITH GALAXIES ; ; if keyword_set(rpos) then obj=long(randomu(seed)*shapecat.n) else obj=0 ; Start at beginning of the catalogue .... or a random point in catalogue message,'Placing galaxies into the image...',/info for i=0L,n_gal-1 do begin ; Loop around through catalogue to generate enough objects (can't use obj=0) obj=obj+1L & if obj ge shapecat.n then obj=1L ; Progress report if (i+1) mod 50 eq 0 then print,string(i+1)+'/'+strtrim(n_gal,1) ; ; Retrieve an object from the (resampled) catalogue ; if (keyword_set(parameter) and keyword_set(sigmae)) then begin ; decomp={w:shapecat.w[obj]*0.04/pixsize,$ ; nf:replicate(shapecat.rsquared[obj],2)*10,$ ; x0:transpose(shapecat.x0[obj,*])*0.04/pixsize} ; endif else begin ; decomp=cat_to_decomp(shapecat,obj) ; decomp.beta=decomp.beta*0.04/pixsize ; Rescale to pixel scale ; decomp.n_pixels=decomp.n_pixels*0.04/pixsize ; (integration within what pixels?) ; endelse ; Decide location of new galaxy on the sky. if keyword_set(rpos) then begin ; Choose a random point (or just outside image, to avoid edge effects) position=[randomu(seed),randomu(seed)]*(n_pix+60)-30 endif else begin position=transpose(shapecat.x[obj,*])*0.04/pixsize endelse ; Make dithered exposures if keyword_set(dither) then position=position-dither ; ; Determine location of new galaxy within the postage stamp. ; decomp.x=(position mod 1)+fix(decomp.n_pixels/2.) ; ; For speed, skip object if it is well away from the image area. ; if fix(position[0]+shapecat.rsquared[obj]*10 lt 1) or $ ; fix(position[0]-shapecat.rsquared[obj]*10 gt n_pix[0]) or $ ; fix(position[1]+shapecat.rsquared[obj]*10 lt 1) or $ ; fix(position[1]-shapecat.rsquared[obj]*10 gt n_pix[1]) then continue ; Also skip objects if they have zero or negative flux if shapecat.flux[obj] le 0. then continue ; Manufacture simulated galaxy. if keyword_set(parameter) then begin ; Use an analytic radial profile galtype=([0,1,2,3,fix(2*randomu(seed))+2])[parameter] rsquared=abs(shapecat.rsquared[obj])*0.04/pixsize recomp_x=position mod 1 recomp=simage_make_analytic_object(galtype, rsquared, shapecat.flux[obj], recomp_x, $ psfimage=psfimage, sigmae=sigmae, seed=seed, $ gamma_input=gamma_input, inst_shear=inst_shear,$ shapeinfo=shapeinfo) ; Enforce a constant object size and flux ;recomp=simage_make_analytic_object(decomp,parameter,2*median(shapecat.rsquared),0.1,psf=psf,sigmae=sigmae,gamma_input=gamma_input,inst_shear=inst_shear,seed=seed,shapeinfo=shapeinfo) ; All objects circular ;recomp=simage_make_analytic_object(decomp,parameter,2*median(shapecat.rsquared),0.1,psfimage=psfimage,sigmae=1e-20,gamma_input=gamma_input,inst_shear=inst_shear,seed=seed,shapeinfo=shapeinfo) endif else begin ; Use full shapelet-simulated objects (and a decomp structure is convenient) decomp = shapelets_shapecat2decomp(shapecat,obj) ; Extract shapelet object from catalogue decomp.beta = decomp.beta*0.04/pixsize ; Rescale to pixel scale decomp.n_pixels= round(decomp.n_pixels*0.04/pixsize) ; (integration within what pixels?) recomp = simage_make_shapelet_object(decomp,$ psf=psf,gamma_input=gamma_input,inst_shear=inst_shear,$ seed=seed,rpos=rpos,shapeinfo=shapeinfo) recomp_x = decomp.x endelse ; Error catching: did anything go wrong? if recomp[0] eq -1 then continue ; Locate object within image pixels. simage_locate_object, recomp, recomp_x, position, n_pix, region ; SHIFTED DL 1 PIXEL SINCE PARIS ; Error catching: is object oustide image after all? if region.l eq -1 then continue ; Error catching: objects with zero flux used to be included and made NaN values? junk=where(finite(recomp) eq 0, nnan) if (nnan gt 0) then continue ; Error catching: exclude a few noisy objects erroneously decomposed ; around the (noisy) edge of the Hubble Deep Field. if max(recomp) gt 100 then print,'Holy bright crap, batman!' ; Add object to image. postage_pix=size(recomp,/dimensions) if (postage_pix[0] ne (1+region.r-region.l)) or (postage_pix[1] ne (1+region.t-region.b)) then begin print,'Holy wrong-size trousers, batman!',$ position[0],position[1],region.l,region.b,recomp_x,$ postage_pix[0],postage_pix[1],region.r,region.t endif else begin totalimage=total(image) image[region.l:region.r,region.b:region.t] = $ image[region.l:region.r,region.b:region.t] + recomp ;if total(image) le totalimage then print,'Holy no extra flux, batman!' endelse ; Plot to screen for debugging. ;plt_image,recomp,/col,/fr & read,junk ; Renormalise coordinates to those used by SExtractor (which seems to make ; a mistake in the x direction) position[0]=position[0]+1 ; Write information to (pre-noise) catalogue of objects. printf,lun,format=mocksexformat, $ 0>(position[0]+0.5)<n_pix[0], $ 0>(position[1]+0.5)<n_pix[0], $ ' 0 0 ', $ (region.l+1)>1,(region.b+1)>1, $ (region.r+1)<n_pix[0],(region.t+1)<n_pix[1], $ shapeinfo.A,shapeinfo.B,shapeinfo.theta,0., $ sqrt(shapeinfo.rsq),shapeinfo.mag,shapeinfo.flux,$ ' 0 ',i+1 endfor ; ; ; PUT STARS INTO IMAGE ; ; if n_star gt 0 then message,'Placing stars in the image',/info for i=0L,n_star-1 do begin ; Progress meter if (i+1) mod 20 eq 0 then print,string(i+1)+'/'+strtrim(n_star,1) ; Choose a random location within the image (or just outside, to avoid edge ; effects) position=[randomu(seed),randomu(seed)]*(n_pix+28)-14 ; Ensure position within pixel is random to allow oversampling later psf.x=position-fix(position)+psf.n_pixels/2. ; Decide which of the pixels in the final image this corresponds to mk_recomp, psf, recomp ; Randomly assign a flux (flux=1 would give a magnitude 22.08 star) starmag=simage_star_magnitude_distribution(seed=seed,slf=slf) starflux=10.^(-0.4*(starmag-22.08)) recomp=recomp*starflux ; Locate pixels in final image that will contain the star simage_locate_object, recomp, psf.x, position, n_pix, region ; Error catching: is object oustide image after all? if region.l eq -1 then continue ; Add object to image. if ((size(recomp))[1] eq (1+region.r-region.l) and (size(recomp))[2] eq (1+region.t-region.b)) then begin image[region.l:region.r,region.b:region.t] = $ image[region.l:region.r,region.b:region.t] + recomp endif else continue ; Renormalise coordinates to those used by SExtractor (which seems to make ; a mistake in the x direction) position[0]=position[0]+1 ; Write information to (pre-noise) catalogue of objects. ; (add 1 to positions to mimic SExtractor's index, which starts at 1 not 0) printf,lun,format=mocksexformat,$ position[0]+0.5,position[1]+0.5,$ ' 0 0 ',$ region.l+1,region.b+1,region.r+1,region.t+1,$ psfinfo.A,psfinfo.B,psfinfo.theta,1.,$ psfinfo.rsq,starmag,starflux,$ ' 0 ',n_gal+i+1 endfor ; Close catalogue file for input parameters. close,lun free_lun,lun ; ; ; OUTPUT FINAL IMAGE ; ; ;Construct fits image header mkhdr, header, image header=header[0:4] header=[header,"OBJECT = 'Shimulation-R.Massey'"] header=[header,"CREATOR = 'Shapelets image simulation http://www.ast.cam.ac.uk/~rjm/shapelets/'"] header=[header,"AUTHOR = 'Richard Massey rjm@ast.cam.ac.uk'"] header=[header,"REFERENC= 'MNRAS submitted, preprint astro-ph/0301449'"] header=[header,"BUNIT = 'photons/second' /"] header=[header,"CTYPE1 = 'RA' /"] header=[header,"CRPIX1 = 0. / Reference point is bottom left hand corner"] if keyword_set(dither) then begin header=[header,"CRVAL1 = "+string(float(dither[0])/60./60.)+" / deg Dither position away from default"] endif else begin header=[header,"CRVAL1 = 0. / deg (at default dither position)"] endelse header=[header,"CDELT1 = "+string(float(pixsize)/60./60.)+" / deg"] header=[header,"CTYPE2 = 'Dec' /"] header=[header,"CRPIX1 = 0. /"] if keyword_set(dither) then begin header=[header,"CRVAL1 = "+string(float(dither[1])/60./60.)+" / deg"] endif else begin header=[header,"CRVAL2 = 0. / deg"] endelse header=[header,"CDELT2 = "+string(float(pixsize)/60./60.)+" / deg"] header=[header,"DATAMAX = "+string(float(max(image)))+" /"] header=[header,"DATAMIN = "+string(float(min(image)))+" /"] header=[header,"PIXSIZE = "+string(float(pixsize))+" / arcsec"] header=[header,"EXPTIME = 'Noise-free image' /"] header=[header,"PHOTOZP = 22.08 /"] if keyword_set(psf) then begin psfname=psf.name+"'" while (strlen(psfname) lt 20) do psfname=psfname+" " header=[header,"PSF = '"+psfname+"/ Plus initial WFPC2 PSF!"] endif else begin header=[header,"PSF = 'NONE' / Only original WFPC2 PSF, circularised"] endelse header=[header,"NSTAR = "+string(fix(n_star),format=(I10))+" / Approximate number of stars"] header=[header,"NGAL = "+string(fix(n_gal),format=(I10))+" / Approximate number of galaxies"] if keyword_set(parameter) then begin objtype=(["'de Vaucouleurs prof'",$ "'Exponential profile'",$ "'Gaussian profile' ",$ "'Mixed analytic' "])[parameter-1] header=[header,"OBJTYPE = "+objtype+"/"] header=[header,"SIGMAE = "+string(float(sigmae))+" / Scatter of object ellipticities"] endif else begin header=[header,"OBJTYPE = 'Full shapelet sim' /"] endelse if keyword_set(gamma_input) then begin header=[header,"GAMMAIN1= "+string(float(gamma_input[0]))+" / Input gravitational lensing shear"] header=[header,"GAMMAIN2= "+string(float(gamma_input[1]))+" /"] endif else begin header=[header,"GAMMAIN1= 0. / Input gravitational lensing shear"] header=[header,"GAMMAIN2= 0. /"] endelse if keyword_set(inst_shear) then begin header=[header,"INSTSH1 = "+string(float(inst_shear[0]))+" / Instrumental distortions"] header=[header,"INSTSH2 = "+string(float(inst_shear[1]))+" /"] endif else begin header=[header,"INSTSH1 = 0. / Instrumental distortions"] header=[header,"INSTSH2 = 0. /"] endelse header=[header,"END "] header=[header," "," "," "," "," "," "," "," "," "," "," "," "," "," "," ",$ " "," "," "," "," "," "," "," "," "," "," "," "," "," "," "," "] ; Write image to file writefits, shapelets_paths(3)+file+'.fits', image, header ; ; ADD NOISE TO IMAGE (Writes a second, separate file) ; if keyword_set(noise) then $ sim_addnoise,shapelets_paths(3)+file,image=image,header=header,time=noise,seed=seed 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!