# Source code for shapelets_psi_rowe.pro:

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\$
;
;
; 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