FUNCTION laguerre, x, nIn, kIn, $
    COEFFICIENTS=coefficients, $
    DOUBLE=double

    COMPILE_OPT idl2

; error checking
    ON_ERROR, 2
    IF (N_PARAMS() LT 2) OR (N_PARAMS() GT 3) THEN $
        MESSAGE,'Incorrect number of arguments.'
    IF ((N_ELEMENTS(nIn) GT 1) OR (N_ELEMENTS(kIn) GT 1)) THEN MESSAGE, $
        'N and K must be scalars.'
    n = FLOOR(nIn[0])
    IF (N_PARAMS() EQ 2) THEN kIn = 0   ; default
    k = kIn[0]
    IF (n LT 0) OR (k LT 0) THEN MESSAGE, $
        'Arguments N and K must be greater than or equal to zero.'

; keyword checking
    tname = SIZE(x,/TNAME)
    IF (N_ELEMENTS(double) LT 1) THEN $
        double = (tname EQ 'DOUBLE') OR (tname EQ 'DCOMPLEX') $
    ELSE double = KEYWORD_SET(double)


; convert x (if necessary) to desired precision
    CASE tname OF
        'DOUBLE': xPrec = double ? x : FLOAT(x)
        'DCOMPLEX': xPrec = double ? x : COMPLEX(x)
        ELSE: xPrec = x  ; no need to convert
    ENDCASE

    one = double ? 1d : 1.0

    IF (n EQ 0) THEN BEGIN
        result = one + 0*xPrec
    ENDIF ELSE BEGIN
        ; need constants for recursive relation
        Lnm1 = one
        result = -xPrec + k + one
        ; construct polynomial coefficients using recursive relation
        ; nL(n,k) = -xL(n-1,k) + (2n+k-1)L(n-1,k) - (n+k-1)L(n-2,k)
        FOR i=2L,n DO BEGIN
            Lnm2 = Lnm1   ; save L(n-2,k)
            Lnm1 = result ; save L(n-1,k)
            result = ((-xPrec + 2*i + k - 1)*Lnm1 - (i + k - 1)*Lnm2)/i
        ENDFOR
    ENDELSE

    m = LINDGEN(n+1)
    coefficients = ((-1d)^m)*GAMMA(n+k+1d)/ $
        (GAMMA(n-m+1d)*GAMMA(k+m+1d)*GAMMA(m+1d))
    IF (NOT double) THEN coefficients = FLOAT(coefficients)
    RETURN, result
EN