function shapelets_psi_rowe, nm, x1, x2, $
BETA=beta, $
ELLIPTICITY=ellipticity, $
THETA_ZERO=theta_zero, $
INTEGRATE=integrate, $
ARRAY=array, $
PREFACTOR=prefactor, $
SILENT=silent, TEST2=test2, TEST3=test3
;$Id: shapelets_psi.pro, v1$
;
; Copyright © 2005 Richard Massey and Alexandre Refregier.
;
; This file is a part of the Shapelets analysis code.
; www.astro.caltech.edu/~rjm/shapelets/
;
; The Shapelets code is free software; you can redistribute it and/or
; modify it under the terms of the GNU General Public Licence as published
; by the Free Software Foundation; either version 2 of the Licence, 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 Licence for more details.
;
; You should have received a copy of the GNU General Public Licence
; 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:
; SHAPELETS_PSI
;
; CATEGORY:
; Shapelets.
;
; PURPOSE:
; Compute (dimensionless) exponential polar shapelet basis functions
; psi(r, phi) in 2D, based on Laguerre function polynomials.
; Dimensionful basis functions can be calculated with
; shapelets_psi(n,m, x1/beta[,x2/beta])/beta.
;
; INPUTS:
; NM - Basis function order. Integer [n,m].
; x1 - Grid of x coordinates (e.g. created by shapelets_make_xarr.pro).
; x2 - Grid of y coordinates if a 2D basis function is required.
; BETA - Shapelet scale size beta.
; ELLIPTICITY- Ellipticity of basis functions, in weak lensing notation:
; either [e1,e2] or complex(e1,e2), where e1=e*cos(2theta) and
; e2=e*sin(2theta), with e=(a-b)/(a+b) for major and minor axes
; and theta is the angle of the major axis, a/c/w from x axis.
; This only works for basis functions not integrated in pixels.
; THETA_ZERO - Basis functions are rotated by this angle, a/c/w from x axis.
;
; KEYWORD PARAMETERS:
; /ARRAY - Instead of just returning one basis function, if /ARRAY is
; set, the routine will return an array containing all of
; the basis functions up to those associated with n_max=n
;
; OUTPUTS:
; The polar exponential shapelet basis function psi_n_m(r, phi)
;
; OPTIONAL IN/OUT KEYWORDS:
; PREFACTOR - If not set on input, will output a table of
; factorial prefactors for optinal reuse next
; time routine is called. Will enlarge each time
; model with a new largest n_max is required (without
; reinventing the wheel each time).
;
; MODIFICATION HISTORY:
; Aug 10 - Prefactor storage table option added by BR
; Jul 10 - Adapted to Joel's X-lets by R. Massey
; Jun 10 - Written by B. Rowe
if not keyword_set(beta) then beta = 1.d0
if keyword_set(integrate) then message,"Sorry, can't integrate exponential shapelets within pixels... maybe I could oversample, though...",/INFO,NOPRINT=silent
if n_elements(nm) ne 2 then message,"Exponential shapelets currently only exist in 2D form"
if min(nm) lt 0 then message,"n must be non-negative!"
;if size(x1,/N_DIMENSIONS) lt 2 or min(size(x1,/DIMENSIONS)) le 1 then x1=replicate(x1, 1) ; Error catching for 1xn arrays
; Establish number of pixels
if not keyword_set(x1) or not keyword_set(x2) then begin ; Default x and y arrays.
shapelets_make_xarr,[128,128],x1,x2
x1=x1*beta/8.
x2=x2*beta/8.
endif
n_pix=size(x1,/DIMENSIONS) & n_pix_x=n_pix[0]
if size(x1,/N_DIMENSIONS) gt 1 then n_pix_y=long(n_pix[1]) else n_pix_y=1
; Evaluate grid coordinate arrays
if keyword_set(ellipticity) then begin ; Make the grid elliptical if necessary
if keyword_set(integrate) then message,"Unable to integrate within pixels for elliptical basis functions!",/INFO,NOPRINT=silent
if size(ellipticity,/TYPE) eq 6 then begin ; (a-b)/(a+b)
e=[float(ellipticity),imaginary(ellipticity)]
endif else if n_elements(ellipticity) ge 2 then begin
e=[ellipticity[0],ellipticity[1]]
endif else message,"Ellipticity format not recognised!"
;eta=alog(1-total(e^2))-alog(1+total(e^2)) ; From BJ02 eq (2.8), with atanh(z)=[ln(1+z)-ln(1-z)]/2
eta=-2*atanh(total(e^2)) ; From BJ02 eq (2.8)
theta=atan(e[1],e[0]) ; PA=theta/2
x1prime = (cosh(eta/2)+cos(theta)*sinh(eta/2))*x1+(sin(theta)*sinh(eta/2))*x2 ; BJ02 eq (2.9)
x2prime = (sin(theta)*sinh(eta/2))*x1+(cosh(eta/2)-cos(theta)*sinh(eta/2))*x2
;x1prime = (1-e[0])*x1 -e[1]*x2
;x2prime = -e[1]*x1 +(1+e[0])*x2
;determinant=1.-total(e^2)
; x1=x1prime;/determinant
; x2=x2prime;/determinant
endif else begin
x1prime=x1
x2prime=x2
endelse
r = sqrt(x1prime^2+x2prime^2)/beta
theta = atan(x2prime,x1prime)
if keyword_set(theta_zero) then theta=temporary(theta)-(theta_zero/!radeg) ; Rotate basis functions if necessary
; Construct basis functions
if not keyword_set(array) then begin
n_coeffs=1
n=nm[0]
m=nm[1]
endif else shapelets_make_nvec,nm[0],n,m,n_coeffs,/POLAR,/EXPONENTIAL
root2pi = sqrt(2. * !pi)
BasisF1 = complexarr(n_pix_x, n_pix_y, n_coeffs)
BasisF2 = complexarr(n_pix_x, n_pix_y, n_coeffs)
BasisF3 = complexarr(n_pix_x, n_pix_y, n_coeffs)
; Create or update/enlarge prefactor storage table if necessary
;
if not keyword_set(prefactor) then begin
prefactor = (beta * (float(n) + 0.5))^(-1-abs(m)) $
* sqrt(factorial(n-abs(m)) / float(2*n+1) /factorial(n+abs(m))) $
/ root2pi
endif else if n_elements(prefactor) lt n_coeffs then begin
pre_old = temporary(prefactor)
n_coeffs_old = n_elements(pre_old)
prefactor = fltarr(n_coeffs)
prefactor[0:n_coeffs_old-1L] = temporary(pre_old)
prefactor[n_coeffs_old:*] = $
(beta * (float(n[n_coeffs_old:*]) + 0.5))^(-1 - abs(m[n_coeffs_old:*])) $
* sqrt(factorial(n[n_coeffs_old:*]-abs(m[n_coeffs_old:*]))/float(2*n[n_coeffs_old:*]+1)/factorial(n[n_coeffs_old:*]+abs(m[n_coeffs_old:*]))) $
/ root2pi
endif
; first do BasisF1 for comparison
;
exponential=exp(-r)
for i=0,n_coeffs-1 do begin
BasisF1[*,*,i]=(prefactor[i]*exp(-1./(2*n[i]+1))) * $
exponential * $
r^abs(m[i]) * $
laguerre(r/(n[i]+0.5),n[i]-abs(m[i]),2*abs(m[i]),/DOUBLE) * $
complex(cos(m[i]*theta),-sin(m[i]*theta))
endfor
; then do BasisF2 for comparison (absorbing exp(-r / (2n + 1)) into
;loop as I believe is necessary)
;
for i=0,n_coeffs-1 do begin
BasisF2[*,*,i]=prefactor[i] * $
exp(-r / (2. * float(n[i]) + 1.)) * $
r^abs(m[i]) * $
laguerre(r/(n[i]+0.5), n[i]-abs(m[i]),2*abs(m[i]),/DOUBLE) * $
complex(cos(m[i]*theta), -sin(m[i]*theta))
endfor
; then BasisF3 (to check my hand-hard-coded polynomials)
;
for i=0,n_coeffs-1 do begin
BasisF3[*,*,i]=prefactor[i] * $
exp(-r / (2. * float(n[i]) + 1.)) * $
r^abs(m[i]) * $
shapelets_kummer(n[i], m[i], r) * $
complex(cos(m[i]*theta), -sin(m[i]*theta))
; Plot these to check that kummer <-> laguerre agree!
;
; plot, r, laguerre(r/(n[i]+0.5), n[i]-abs(m[i]),2*abs(m[i]),/DOUBLE) - $
; shapelets_kummer(n[i], m[i], r), CHARSIZE=1.2
; stop
endfor
test2 = BasisF2
test3 = BasisF3
return, BasisF2
end
Return to the shapelets web page or the code help menu.
| Last modified on 2nd Mar 2009 by Richard Massey. |