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

You can also view the help page for this routine.

pro simage_generate_shapecat, shapecat, seed=seed, output_file=output_file, width=width, $ n_pix=n_pix, n_max=n_max, n_obj=n_obj, plot=plot ;$Id: simage_generate_shapecat.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_GENERATE_SHAPECAT ; ; CATEGORY: ; Shapelets image simulation. ; ; PURPOSE: ; Uses simage_resample_pdf.pro many times to generate a new random ; catalogue of objects that could be used to populate a simulated image. ; ; INPUTS: ; shapecat - Shapelet catalogue of real data in a similar regime to the ; required simulation (eg HDF galaxies, read in using ; shapelets_read_shapecat_hdf.pro). ; It is recommended, although not necessary, that this contain ; polar rather than Cartesian shapelet coefficients. Objects ; just morph better in polar shapelets. ; ; OPTIONAL INPUTS: ; output_file - Filename for output. ; width - This is slow to calculate, so allow it to be entered by ; hand, or calculated on a previous run to avoid the ; duplication of effort. ; seed - Seed for random number generator. ; n_pix - Integer (x,y). Number of pixels planned for final ; simulated image. This routine does not create the image ; (for that, see simage_assemble_image.pro), but the random ; positions of simulated galaxies can be specified in either ; routine. If they are specified here, and stored in the ; catalogue, it will be easier to duplicate simulated image ; or created dithered simulateed image later. ; Default: . ; n_max - Objects created up to this truncation in shapelet space. ; Default: . ; n_obj - Number of objects in output catalogue. This does not ; determine the number density of objects in the simulated ; image. That is done in simage_assemble_image.pro. ; Default: . ; ; KEYWORD PARAMETERS: ; /PLOT - Plot each newly generated object to screen, for debugging. ; ; OUTPUTS: ; A simulated shapelet catalogue is written to disc. ; ; MODIFICATION HISTORY: ; Mar 04 - Tidied up by RM. ; Mar 03 - Cosmetic changes and modularisation of subroutines by RM. ; Jan 02 - new_hdf.pro written by Richard Massey. ;- ; ; Set up defaults ; pixel_scale=0.04 if not keyword_set(shapecat) then shapelets_read_shapecat_hdf,shapecat,/polar,/moments; Read in source catalogue of morphology distributions if (min(shapecat.rsquared) eq 0) then message,'Use shapelets_read_shapecat with /MOMENTS flag!' if not keyword_set(seed) then seed = 1000L ; Random number seed if not keyword_set(output_file) then $ ; Output file name for new catalogue output_file="new"+STRMID(SYSTIME(0), 14, 2)+STRMID(SYSTIME(0), 17, 2) if not keyword_set(n_pix) then n_pix = 4096 ; Number of pixels in final image is n_pix x n_pix if not keyword_set(n_obj) then n_obj = 10000 ; number of new objects to synthesise... if not keyword_set(n_max) then n_max=shapecat.maxn_max ; ...to this order shapelet (this is the maximum ; n_max of all the objects in the catalogue) ; recommended to use this so eg nuclei of objects ; aren`t simply enlarged by size-mag plane fitting. shapelets_make_nvec,n_max,n1,n2,n_coeffs ; message,strtrim(n_obj,2)+' galaxies being created in '+shapelets_paths(2)+output_file+'.shape',/info if n_elements(seed) eq 1 then $ ; message,'Random seed='+strtrim(seed[0],2),/info ; Friendly messages with setup if not keyword_set(width) then begin ; Calculate smoothing width before opening output files message,'Calculating width for kernel smoothing in morphology PDF',/info simage_resample_pdf, shapecat, junk_obj, n_max=n_max, seed=seed, width=width, /silent endif ; ; ; Open output file and write out header information (need to change the 999 if n_max>43): ; openw, lun,shapelets_paths(2,/SILENT)+output_file+'.shape',/get_lun printf,lun,'# Brand new synthesised shapelet galaxies' printf,lun,'# num objs out of n_coeffs n_coeffs seeing n_pixels(x) n_pixels(y) pixel scale' printf,lun,format='("#",I11,3(I13),F13.6,2(I13),F13.6)',$ n_obj,n_obj,n_max,n_coeffs,shapecat.seeing,n_pix,n_pix,pixel_scale printf,lun,'#' printf,lun,'# x0 y0 beta SEx fwhm SEx mag n_max'+$ ' Flag SEx ID Chi squared SEx class SEx A SEx B'+$ ' SEx theta SEx e1 SEx e2 SEx flux SEx x SEx y'+$ ' SEx area Object:x_min x_max y_min y_max' printf,lun,format='($,"# n1",I8," ",(999(I12, :, " ")))',n1 & printf,lun printf,lun,format='($,"# n2",I8," ",(999(I12, :, " ")))',n2 & printf,lun ; ; Main loop to generate new objects ; for i=1L,n_obj do begin ; Progress report if (i) mod 50 eq 0 then print,string(i)+'/'+strtrim(n_obj,1) ; Create an object new_object={n_max:0} while new_object.n_max eq 0 do begin bs_obj=0 simage_resample_pdf, shapecat, new_object, n_max=n_max, seed=seed, width=width, /SILENT, bs_obj=bs_obj ,plot=plot endwhile ; Extract results from structure coeffs = new_object.coeffs error = new_object.error if n_elements(coeffs) gt 1 then begin coeffs=[coeffs[0],coeffs[1:*]/coeffs[0]] if error[0] ne 0. then error=[error[0],error[1:*]/error[0]] endif while ( n_elements(coeffs) lt n_coeffs ) do begin coeffs=[coeffs,0.] error =[error ,0.] endwhile ; Place object randomly around the sky ; (or just outside the image, to enable dithering and avoid edge effects) x = [randomu(seed),randomu(seed)]*(n_pix*3./2.)-[n_pix/2.+20,20] ; Write out results to file otherdata = [x,new_object.beta,new_object.moments.rsquared,new_object.moments.mag,new_object.n_max,10*shapecat.flag[bs_obj,0]+shapecat.flag[bs_obj,1],shapecat.sexid[bs_obj],new_object.chisq,$ 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0,0,0,0,0] printf,lun,format='($,(5(F12.5, :, " ")),(3(I12, :, " ")),(10(F12.5, :, " ")),(5(I12, :, " ")))',otherdata & printf,lun printf,lun,format='($,(999(E12.5, :, " ")))',coeffs & printf,lun printf,lun,format='($,(999(E12.5, :, " ")))',error & printf,lun endfor ; Tidy up and go home printf,lun,'# End' close,lun free_lun,lun 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!