# Source code for shapelets_kummer.pro:

```;\$Id: shapelets_kummer.pro, v2\$
;
; 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.

function shapelets_kummer, n, m, x
;
; NAME:
;       SHAPELETS_KUMMER
;
; CATEGORY:
;       Exponential shapelets.
;
; PURPOSE:
;       Compute the polynomial M(-(n-|m|), 2|m|+1, x/(n+0.5)), the
;       Kummer's function needed for exponential shapelets
;
; INPUTS:
;       n - hydro shapelet n value
;       m - hydro shapelet m value
;       x - can be scalar, vector or array
;
;
; OUTPUTS:
;       M(-(n-|m|), 2|m|+1, x/(n+0.5)) - Kummer's function evaluated at
;                                position x.  If x = r/beta is of integer type,
;                                this will be a longword/64bit
;                                integer, otherwise it will be a DBL
;                                precision floating point.
;
; KEYWORD PARAMETERS:
;       None.
;
; NOTES:
;       Sums using vector of coefficients ci for M(a,b,x)=Sum_i ci*x^i, i=0,..,N
;
; MODIFICATION HISTORY:
;
;       Aug 10 - RM redirected n>6 calls to laguerre.pro. Slow but works.
;       Aug 10 - Richard Massey reverted to FOR loop.
;       Aug 10 - BR vectorized polynomial sum rather than using (FOR loop
;                invoking) POLY.PRO... Needs to be hard coded this
;                way, but is much faster in IDL
;       Jun 10 - Written by B. Rowe
;
COMPILE_OPT idl2, HIDDEN
ON_ERROR, 2

print,n,m

if abs(m) gt n then message,"Must have |m| < n!"
if n gt 6 then begin
message,/INFO,"Polynomials not yet hardcoded for m>6. Using built-in laguerre.pro"
return, laguerre(x/(n+0.5),n-abs(m),2*abs(m),/DOUBLE)
endif else begin
case n of
0: return, x * 0 + 1
1: begin
case abs(m) of
1: return, x * 0 + 1
0: c = [1.d0, -2.d0 / 3.d0]
endcase
end
2: begin
case abs(m) of
2: return, x * 0 + 1
1: c = [3.d0, -2.d0 / 5.d0]
0: c = [1.d0, -4.d0 / 5.d0, +2.d0 / 25.d0]
endcase
end
3: begin
case abs(m) of
3: return, x * 0 + 1
2: c = [5.d0, -2.d0 / 7.d0]
1: c = [6.d0, -8.d0 / 7.d0, +2.d0 / 49.d0]
0: c = [1.d0, -6.d0 / 7.d0, +6.d0 / 49.d0, -4.d0 / 1029.d0]
endcase
end
4: begin
case abs(m) of
4: return, x * 0 + 1
3: c = [7.d0,  -2.d0 / 9.d0]
2: c = [15.d0, -4.d0 / 3.d0,  +2.d0 / 81.d0]
1: c = [10.d0, -20.d0 / 9.d0, +10.d0 / 81.d0, -4.d0 / 2187.d0]
0: c = [1.d0,  -8.d0 / 9.d0,  +4.d0 / 27.d0,  -16.d0 / 2187.d0, \$
2.d0 / 19683.d0]
endcase
end
5: begin
case abs(m) of
5: return, x * 0 + 1
4: c = [9.d0,  -2./11.d0]
3: c = [28.d0, -16.d0 / 11.d0, +2. / 121.d0]
2: c = [35.d0, -42.d0 / 11.d0, +14.d0 / 121.d0, -4.d0 / 3993.d0]
1: c = [15.d0, -40.d0 / 11.d0, +30.d0 / 121.d0, -8.d0 / 1331.d0,  \$
+2.d0 / 43923.d0]
0: c = [1.d0,  -10.d0 / 11.d0, +20. / 121.d0,   -40.d0 / 3993.d0, \$
+10.d0 / 43923.d0, -4.d0 / 2415765.d0]
endcase
end
6: begin
case abs(m) of
6: return, x * 0 + 1
5: c = [11.d0, -2.d0 / 13.d0]
4: c = [45.d0, -20.d0 / 13.d0,  +2.d0 / 169.d0]
3: c = [84.d0, -72.d0 / 13.d0,  +18.d0 / 169.d0,  -4.d0 / 6591.d0]
2: c = [70.d0, -112.d0 / 13.d0, +56.d0 / 169.d0, -32.d0 / 6591.d0, \$
+2.d0 / 85683.d0]
1: c = [21.d0, -70.d0 / 13.d0,  +70.d0 / 169.d0, -28.d0 / 2197.d0, \$
+14.d0 / 85683.d0,  -4.d0 / 5569395.d0]
0: c = [1.d0,  -12.d0 / 13.d0,  +30.d0 / 169.d0, -80.d0 / 6591.d0, \$
+10.d0 / 28561.d0, -8.d0 / 1856465.d0,  +4.d0 / 217206405.d0]
endcase
end
endcase

; Trick from poly.pro uses a FOR loop but Richard thinks it's actually faster
out=(reverse(c))[0] & for i=n_elements(c)-2,0,-1 do out = temporary(out) * x + c[i]
; Barney's hardcoded version
;case n_elements(c) of
;  2:  out = x * c[1] \$
;          + c[0]
;
;  3:  out = x * x * c[2] \$
;          + x * c[1] \$
;          + c[0]
;
;  4:  out = x * x * x * c[3] \$
;          + x * x * c[2] \$
;          + x * c[1] \$
;          + c[0]
;
;  5:  out = x * x * x * x * c[4] \$
;          + x * x * x * c[3] \$
;          + x * x * c[2] \$
;          + x * c[1] \$
;          + c[0]
;
;  6:  out = x * x * x * x * x * c[5] \$
;          + x * x * x * x * c[4] \$
;          + x * x * x * c[3] \$
;          + x * x * c[2] \$
;          + x * c[1] \$
;          + c[0]
;
;  7:  out = x * x * x * x * x * x * c[6] \$
;          + x * x * x * x * x * c[5] \$
;          + x * x * x * x * c[4] \$
;          + x * x * x * c[3] \$
;          + x * x * c[2] \$
;          + x * c[1] \$
;          + c[0]
;endcase

endelse

return, out

end
```

Return to the shapelets web page or the code help menu.

 Last modified on 2nd Mar 2009 by Richard Massey.