; Copyright 2007 IAS ; This file is part of the Planck Sky Model. ; ; The Planck Sky Model 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; version 2 of the License. ; ; The Planck Sky Model 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 Planck Sky Model. If not, write to the Free Software ; Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA ;+ ; add smallscale structure ; ; Use of properly normalised Gaussian field to fill in the small scale structure ; ; @param inputmap {in}{required}{type=string} ... ; @keyword struct {in}{ptionnal}{type=struct} if set, extract Ancillary data information from structure ; @keyword seed{in}{optional}{type=int} ... ; @keyword planckres{in}{optional}{type=int} ... ; @keyword mask{in}{optional}{type=int} ... ; @keyword beta{in}{optional}{type=int} ... ; @keyword alpha{in}{optional}{type=int} ... ; @keyword gamma{in}{optional}{type=int} ... ; @keyword resolution{in}{optional}{type=int} ... ; @keyword polar{in}{optional}{type=boolean} ... ; ; @history

;- function add_smallscale, inputmap, struct, seed=seed, planckres=planckres, $ mask=mask, beta=beta, alpha=alpha, gamma=gamma, $ resolution=resolution, polar=polar, nl=nl, ordering=ordering if not keyword_set(planckres) then planckres = (*!sky_model).planckresolution if (planckres le 0) then planckres = 5. ; arcmin if not keyword_set(resolution) then resolution=0 if not keyword_set(alpha) then alpha=0 if not keyword_set(beta) then beta=0 if not keyword_set(gamma) then gamma=0 if not keyword_set(mask) then mask=0 IF NOT KEYWORD_SET(ordering) THEN ordering='NESTED' ; extract Ancillary data information from structure if keyword_set(struct) then begin tags = tag_names(struct) for i=0, n_elements(tags)-1 do begin case tags(i) of 'RESOLUTION': resolution = struct.resolution 'ALPHA': alpha = struct.alpha 'BETA': beta = struct.beta 'GAMMA': gamma = struct.gamma else : endcase endfor endif ; make sure all parameters are non zero check_param = alpha*beta*gamma*resolution if (check_param eq 0) then begin ; print, "could not add small scales, parameter missing" return, inputmap endif ; make sure resolution of template is not higher than planckres if (resolution le planckres) then begin PSM_INFO_MESSAGE, "Could not add small scales: resolution of template higher than requested resolution" return, inputmap endif ; nside of input map npix = (SIZE(inputmap))[1] nside = npix2nside(npix) ; create Gaussian random field IF NOT KEYWORD_SET(nl) THEN nl = (*!sky_model).lmax + 1 l = (findgen(nl-1)+1) sigma_planck = planckres/60./2.354*!pi/180. ; sigma of Planck in radian sigma_input = resolution/60./2.354*!pi/180. ; sigma of input map in radian cl = l^(gamma) * ( exp(-l^2 * sigma_planck^2) - exp(-l^2 * sigma_input^2) ) ; map = synfast(cl, nside, /verbose, seed=getseed()) ; map = synfast_cxx(cl, nside, /verbose, seed=1000) ; map = reorder(map, in='RING', out='NESTED') l = [0,l] cl = [0,cl] GENERATE_X_LM, l, cl, xlm PSM_ALM2MAP, xlm, map, ordering=ordering, hpx_nside=nside, /nohdr, /overwrite map = map / stddev(map) map = map - avg(map) ; modulate small scale fluctuation by large scale intensity ; remove minimum from input map to make sure only positive values are ; exponentiated ; ; If polarization, the normalisation is done with abs(inputmap) if keyword_set(polar) then begin map = map*alpha*(abs(inputmap))^beta endif else begin ; if (min(inputmap) lt 0.) then zerolevel = min(inputmap) else zerolevel=0. ; map = map*alpha*(inputmap-zerolevel)^beta map = map*alpha*(inputmap>0)^beta endelse ; add small scale to original map ; deal with mask if provided (add structure only where mask is equal to 1) if keyword_set(mask) then begin read_fits_map, struct.mask, tempo, hdr, exthdr if (n_elements(exthdr) gt 1) then order_in = sxpar(exthdr, 'ORDERING') else order_in='RING' ud_grade, tempo, mask, nside_out=nside, order_in=order_in, order_out=ordering tempo=0 outputmap = inputmap ind = where(mask eq 1, nbind) if (nbind gt 1) then outputmap(ind) = inputmap(ind)+map(ind) endif else begin outputmap = inputmap+map endelse return, outputmap end