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.