function besselroots(nu::Float64, n::Int64) 
#BESSELROOTS    The first N roots of the function J_v(x)

# DEVELOPERS NOTES:
#   V = 0 --> Full double precision for N <= 20 (Wolfram Alpha), and very
#     accurate approximations for N > 20 (McMahon's expansion)
#   -1 <= V <= 5 : V ~= 0 --> 12 decimal figures for the 6 first zeros
#     (Piessens's Chebyshev series approximations), and very accurate
#     approximations for the others (McMahon's expansion)
#   V > 5 --> moderately accurate for the 6 first zeros and good
#     approximations for the others (McMahon's expansion)

# This code was originally written by L. L. Peixoto in MATLAB. 
# Later modified by A. Townsend to work in Julia
    
if ( n == 0 )
    x = []
elseif( n > 0 && nu == 0 )
    x = McMahon( nu, n )
    correctFirstFew = Tabulate( )
    x[1:min(n,20)] = correctFirstFew[1:min(n,20)]
    x
elseif ( n>0 && nu >= -1 && nu <= 5 ) 
    x = McMahon( nu, n )
    correctFirstFew = Piessens( nu )
    x[1:min(n,6)] = correctFirstFew[1:min(n,6)]
    x
else
    error("Input N must be a positive integer")
end

end


function McMahon( nu::Float64, n::Int64 )
# McMahon's expansion. This expansion gives very accurate approximation 
# for the sth zero (s >= 7) in the whole region V >=- 1, and moderate
# approximation in other cases.
b = Array(Float64,n)
mu = 4*nu^2
a1 = 1. / 8.
a3 = (7.*mu-31.) / 384.
a5 = 4.*(3779.+mu*(-982.+83.*mu)) / 61440. # Evaluate via Horner's method.
a7 = 6.*(-6277237.+mu*(1585743.+mu*(-153855.+6949.*mu))) / 20643840.;
a9 = 144.*(2092163573.+mu*(-512062548.+mu*(48010494.+mu*(-2479316.+70197.*mu)))) / 11890851840.
a11 = 720.*(-8249725736393.+mu*(1982611456181.+mu*(-179289628602.+mu*(8903961290. + 
      mu*(-287149133.+5592657.*mu))))) / 10463949619200.
a13 = 576.*(423748443625564327. + mu*(-100847472093088506.+mu*(8929489333108377. + 
    mu*(-426353946885548.+mu*(13172003634537.+mu*(-291245357370. + mu*4148944183.)))))) / 13059009124761600.
for k=1:n b[k] = .25*(2*nu+4*k-1)*pi end # beta
# Evaluate using Horner's scheme: 
x = b - (mu-1)*( ((((((a13./b.^2 + a11)./b.^2 + a9)./b.^2 + a7)./b.^2 + a5)./b.^2 + a3)./b.^2 + a1)./b)    
end


function Tabulate( )
# First 20 roots of J0(x) are precomputed (using Wolfram Alpha):
y = Array(Float64,20)
y = [   2.4048255576957728
        5.5200781102863106
        8.6537279129110122
        11.791534439014281
        14.930917708487785
        18.071063967910922
        21.211636629879258
        24.352471530749302
        27.493479132040254
        30.634606468431975
        33.775820213573568
        36.917098353664044
        40.058425764628239
        43.199791713176730
        46.341188371661814
        49.482609897397817
        52.624051841114996
        55.765510755019979
        58.906983926080942
        62.048469190227170  ]
end

function Piessens( nu::Float64 ) 
# Piessens's Chebyshev series approximations (1984). Calculates the 6 first
# zeros to at least 12 decimal figures in region -1 <= V <= 5:
    y = Array(Float64,6);

    C = [
       2.883975316228  8.263194332307 11.493871452173 14.689036505931 17.866882871378 21.034784308088
       0.767665211539  4.209200330779  4.317988625384  4.387437455306  4.435717974422  4.471319438161
      -0.086538804759 -0.164644722483 -0.130667664397 -0.109469595763 -0.094492317231 -0.083234240394
       0.020433979038  0.039764618826  0.023009510531  0.015359574754  0.011070071951  0.008388073020
      -0.006103761347 -0.011799527177 -0.004987164201 -0.002655024938 -0.001598668225 -0.001042443435
       0.002046841322  0.003893555229  0.001204453026  0.000511852711  0.000257620149  0.000144611721
      -0.000734476579 -0.001369989689 -0.000310786051 -0.000105522473 -0.000044416219 -0.000021469973
       0.000275336751  0.000503054700  0.000083834770  0.000022761626  0.000008016197  0.000003337753
      -0.000106375704 -0.000190381770 -0.000023343325 -0.000005071979 -0.000001495224 -0.000000536428
       0.000042003336  0.000073681222  0.000006655551  0.000001158094  0.000000285903  0.000000088402
      -0.000016858623 -0.000029010830 -0.000001932603 -0.000000269480 -0.000000055734 -0.000000014856
       0.000006852440  0.000011579131  0.000000569367  0.000000063657  0.000000011033  0.000000002536
      -0.000002813300 -0.000004672877 -0.000000169722 -0.000000015222 -0.000000002212 -0.000000000438
       0.000001164419  0.000001903082  0.000000051084  0.000000003677  0.000000000448  0.000000000077
      -0.000000485189 -0.000000781030 -0.000000015501 -0.000000000896 -0.000000000092 -0.000000000014
       0.000000203309  0.000000322648  0.000000004736  0.000000000220  0.000000000019  0.000000000002
      -0.000000085602 -0.000000134047 -0.000000001456 -0.000000000054 -0.000000000004               0.
       0.000000036192  0.000000055969  0.000000000450  0.000000000013               0.               0.
      -0.000000015357 -0.000000023472 -0.000000000140 -0.000000000003               0.               0.
       0.000000006537  0.000000009882  0.000000000043  0.000000000001               0.               0.
      -0.000000002791 -0.000000004175 -0.000000000014               0.               0.               0.
       0.000000001194  0.000000001770  0.000000000004               0.               0.               0.
      -0.000000000512 -0.000000000752              0.               0.              0.               0.
       0.000000000220  0.000000000321               0.               0.               0.               0.
      -0.000000000095 -0.000000000137               0.               0.               0.               0.
       0.000000000041  0.000000000059               0.               0.               0.               0.
      -0.000000000018 -0.000000000025               0.               0.               0.               0.
       0.000000000008  0.000000000011               0.               0.               0.               0.
      -0.000000000003 -0.000000000005               0.               0.               0.               0.
       0.000000000001  0.000000000002               0.               0.               0.               0.]
    
    T = Array(Float64,size(C,1))
    pt = (nu-2)/3
    T[1], T[2] = 1., pt
    for k = 2:size(C,1)-1
        T[k+1] = 2*pt*T[k] - T[k-1]
    end
    y[1:6] = T'*C;
    y[1] *= sqrt(nu+1);                  # Scale the first root.
    y
end