pro shapelets_plot_basis, n_max, $
POLAR=polar, $
EXPONENTIAL=exponential, $
CLOG=clog, $
BASIS_ELLIPTICITY=basis_ellipticity, $
THETA_ZERO=theta_zero,$
PS=ps,$
COLOR=color,$
ISOTROPIC=isotropic,$
KSB=ksb,$
TITLE=title,$
REAL=real,IMAGINARY=imaginary
;+
; NAME:
; SHAPELETS_PLOT_BASIS
;
; CATEGORY:
; Shapelets.
;
; PURPOSE:
; Makes a plot of all the polar shapelet basis functions out to n_max.
;
; INPUTS:
; None.
;
; OPTIONAL INPUTS:
; N_MAX - maximum shapelet order.
;
; KEYWORD PARAMETERS:
; PS - prints to (postscript) file instead of screen
; COLOR - with /ps makes a colour postscript file
; KSB - plots only those basis functions which are equivalently used in KSB
; REAL - plots only the real parts of the basis functions
; IMAGE - plots only the imaginary parts of the basis functions
;
; OUTPUTS:
; Plot drawn to STDOUT.
;
; MODIFICATION HISTORY:
; Jul 07 - Fixed to work with older versions of IDL by RM.
; Jun 04 - Real/imaginary options added by RM
; Oct 02 - Written by Richard Massey
;-
;
; Parse inputs
;
if keyword_set(real) or keyword_set(imaginary) then polar=1B
if not keyword_set(n_max) then if keyword_set(ksb) then n_max=4 else n_max=6
;
; Set plotting options
;
n_pixels=[45,45]*2 ; Resolution within each basis function's box
beta=4.1/8 ; beta [pixels]
bits_per_pixel=5 ; Number of greyscale=2^bits
integrate=1 ; Analytically integrate basis functions within pixels?
xmargin=!x.margin ; Remember input margins, so they can be restored later
ymargin=!y.margin ;
!x.margin=[6,1] ; Set less conservative margins
!y.margin=[3,1] ;
col=1B ; Draw a colour bar
if keyword_set(exponential) then begin
if keyword_set(imaginary) then begin
ytitle="!8m!6"
xtitle="!8n!6"
col=0B
frame=[-.5,n_max+.5,-n_max-0.5,n_max+0.5]
endif else if keyword_set(polar) then begin
if not keyword_set(title) then title="!6Polar exponential shapelet basis functions"
ytitle="!8m!6"
frame=[-.5,n_max+.5,-n_max-0.5,n_max+0.5]
endif else begin
if not keyword_set(title) then title="!6Cartesian exponential shapelet basis functions"
xtitle="!8n!6!d1!n"
ytitle="!8n!6!d2!n"
frame=[-.5,n_max+.5,-.5,n_max+.5]
endelse
endif else begin
if keyword_set(imaginary) then begin
ytitle="!8m!6"
xtitle="!8n!6"
col=0B
frame=[-.5,n_max+.5,-n_max-1,n_max+1]
endif else if keyword_set(polar) then begin
if not keyword_set(title) then title="!6Polar shapelet basis functions"
ytitle="!8m!6"
frame=[-.5,n_max+.5,-n_max-1,n_max+1]
endif else begin
if not keyword_set(title) then title="!6Cartesian shapelet basis functions"
xtitle="!8n!6!d1!n"
ytitle="!8n!6!d2!n"
frame=[-.5,n_max+.5,-.5,n_max+.5]
endelse
endelse
;theta_zero=10
;basis_ellipticity=[0.5,0]
; Create empty decomp template
decomp=shapelets_create_decomp(n_max,beta=beta,x=n_pixels/2.,n_pixels=n_pixels, $
POLAR=polar, $
EXPONENTIAL=exponential, $
INTEGRATE=integrate, $
BASIS_ELLIPTICITY=basis_ellipticity, $
THETA_ZERO=theta_zero)
;
; Generate images of the basis functions
;
start_time=systime(/s)
if not keyword_set(polar) then begin ; Cartesian shapelets
if keyword_set(exponential) then message,"Cartesian exponential shapelets are not yet defined" ; Cartesian exponential shapelets
basis=fltarr(n_pixels[0],n_pixels[1],decomp.n_coeffs)
for i=0,decomp.n_coeffs-1 do begin
decomp.coeffs[*]=0.
decomp.coeffs[i]=1
shapelets_recomp,decomp,recomp
basis[*,*,i]=recomp
endfor
endif else begin ; Polar shapelets
basis=complexarr(n_pixels[0],n_pixels[1],decomp.n_coeffs)
for i=0,decomp.n_coeffs-1 do begin
decomp.coeffs[*]=0.
if keyword_set(imaginary) or (keyword_set(polar) and not keyword_set(real) and decomp.m[i] lt 0) then begin
decomp.coeffs[i]=complex(0,-1)
endif else begin
decomp.coeffs[i]=complex(1,0)
endelse
shapelets_recomp,decomp,recomp,/COMPLEX,/SILENT
print,decomp.n[i],decomp.m[i],total(float(recomp))
basis[*,*,i]=temporary(recomp)
endfor
endelse
print,"Taken "+strtrim(systime(/s)-start_time,2)+" seconds"
;
; Rescale image from average within pixels to peak value
;
cscale=1./sqrt(!pi)/max(abs(basis[*,*,0]))
;if strupcase(decomp.profile) ne "EXPONENTIAL" then $
basis=basis*cscale
;
; Remove all basis functions except those used equivalently in KSB
;
if keyword_set(ksb) then begin
basis[*,*,1:2]=0
basis[*,*,6:11]=0
basis[*,*,13:*]=0
if not keyword_set(title) then title="!6Shapelet basis functions used by KSB"
endif
;
; Assemble overall image
;
; Slightly wasteful mathod perhaps of pixellating entire plot, rather than just
; that triangle which will contain data. Still, at least it's then not
; transparent.
xpix=n_pixels[0]*(n_max+1)
ypix=n_pixels[1]*((1+keyword_set(exponential))*n_max+1)
image=fltarr(xpix,ypix)
for i=0,decomp.n_coeffs-1 do begin
; Decide where the basis function should go
if keyword_set(polar) then begin
left=decomp.n[i]*n_pixels[0] & right=left+n_pixels[0]-1
bottom=(ypix/2)+((1+keyword_set(exponential))*decomp.m[i]-1)*n_pixels[1]/2 & top=bottom+n_pixels[1]-1
endif else begin
left=decomp.n1[i]*n_pixels[0] & right=left+n_pixels[0]-1
bottom=decomp.n2[i]*n_pixels[1] & top=bottom+n_pixels[1]-1
endelse
; And place it there
image[left:right,bottom:top]=temporary(image[left:right,bottom:top])+float(basis[*,*,i])
endfor
;
; Open postscript file if requested
;
if keyword_set(ps) then begin
csize=0.1
set_plot,'ps'
if keyword_set(exponential) then begin
device,filename='polar_exponential_basis.eps',xsize=30,ysize=30.,/color,/encap,bits=bits_per_pixel
endif else if keyword_set(imaginary) then begin
device,filename='polar_basis_imaginary.eps',yoffset=5.,ysize=20.,/encap,bits=bits_per_pixel
endif else if keyword_set(real) then begin
device,filename='polar_basis_real.eps',yoffset=5.,ysize=20.,/encap,bits=bits_per_pixel
endif else if keyword_set(polar) then begin
device,filename='polar_basis.eps',yoffset=5.,ysize=20.,/encap,bits=bits_per_pixel
endif else begin
device,filename='Cartesian_basis.eps',yoffset=5.,ysize=20.,/encap,bits=bits_per_pixel
endelse
opub
endif else window,0,xsize=700,ysize=1500
;loadct,1,/silent
;
; Set colour scale
;
if keyword_set(clog) then begin
image=[[[alog10(abs(image>1e-8))]],[[alog10(abs(image)>1e-8)]],[[alog10(abs((-image)>1e-8))]]]
true=3
endif else crange=[-1,1]/sqrt(!pi)
;crange=[min(image),-1*min(image)]
;crange=max(abs(image))*[-1,1]
;
; Plot basis functions
;
shapelets_plot_image,image,cbar=(1-keyword_set(imaginary))*(1-keyword_set(clog)),$
csize=csize,crange=crange,frame=frame,isotropic=1, $
xtitle=xtitle,ytitle=ytitle,title=title,true=true
; Plot boxes around functions
if keyword_set(exponential) then begin
if keyword_set(polar) then begin
for i=0,n_max do begin
oplot,[i-.5,i-.5],[-i-0.5,i+0.5],psym=-3
for j=0,2*i+1 do begin
oplot,[i-.5,i+.5],replicate(-0.5-i+j,2),psym=-3
endfor
endfor
oplot,[n_max+.5,n_max+.5],[-0.5-n_max,0.5+n_max],psym=-3
endif else message,"Cartesian exponential shapelets not defined"
endif else begin
if keyword_set(polar) then begin
for i=0,n_max do begin
oplot,[i-.5,i-.5],[-1-i,1+i],psym=-3
for j=0,i+1 do begin
oplot,[i-.5,i+.5],[-1*i-1+j*2,-1*i-1+j*2],psym=-3
endfor
endfor
oplot,[n_max+.5,n_max+.5],[-1*n_max-1,n_max+1],psym=-3
endif else begin
for i=0,n_max-1 do begin
oplot,[i+.5,i+.5],[-.5,n_max-i+.5],psym=-3
endfor
for i=0,n_max-1 do begin
oplot,[-.5,n_max-i+.5],[i+.5,i+.5],psym=-3
endfor
endelse
endelse
; Print labels
if keyword_set(real) then begin
xyouts,-0.3,n_max+0.1,"!6Real components",charsize=2.25
endif else if keyword_set(imaginary) then begin
xyouts,-0.3,n_max+0.1,"!6Imaginary components",charsize=2.25
endif
; Close ps file if appropriate
if keyword_set(ps) then cps
; Restore margins
!x.margin=xmargin
!y.margin=ymargin
end
Return to the shapelets web page or the code help menu.
| Last modified on 2nd Mar 2009 by Richard Massey. |