#/***********************************************************
#    bessel.rb -- Bessel (ベッセル) 関数
#***********************************************************/

EPS = 1e-10
def odd(x); (x & 1) == 1 ? true : false; end    # 奇数? 
PI = 3.14159265358979324
EULER =  0.577215664901532861

def besJ( n,  x )   # $J_n(x)$ 
    x_2 = x / 2.0

    if (x < 0) 
        if (odd(n)); return -besJ(n, -x)
        else;        return  besJ(n, -x)
        end
    end
    if (n < 0) 
        if (odd(n)); return -besJ(-n, x)
        else;        return  besJ(-n, x)
        end
    end
    if (x == 0); return (n == 0) ? 1: 0; end
    a = s = 0;  b = 1
    k = n;  if (k < x); k = x; end
    begin; k += 1;  end while ((b *= x_2 / k) > EPS)
    if (odd(k)); k += 1; end  # 奇数なら偶数にする 
    while (k > 0) 
        s += b
        a = 2 * k * b / x - a;  k -= 1   # $a = J_k(x)$ 
        if (n == k); r = a; end          # $k$ 奇数 
        b = 2 * k * a / x - b;  k -= 1   # $b = J_k(x)$ 
        if (n == k); r = b; end          # $k$ 偶数 
    end
    return r / (2 * s + b)
        # $J_0 + 2(J_2 + J_4 + \cdots) = 1$ となるように規格化 
end

=begin  #**** 参考: 級数展開版 ****

def besJ2( n,  x)  # $J_n(x)$ 
    if (x < 0) 
        if (odd(n)); return -besJ2(n, -x)
        else;        return  besJ2(n, -x)
        end
    end
    if (n < 0) 
        if (odd(n)); return -besJ2(-n, x)
        else;        return  besJ2(-n, x)
        end
    end
    x /= 2;  result = 1
    for k in 1..n; result *= x / k; end
    x = - x * x;  term = result
    for k in 1...500
        term *= x / (k * (n + k))
        previous = result;  result += term
        if (result == previous); return result; end
    end
    printf("BesJ2(n, x): 収束しません.\n")
    return result
end

=end

def besY( n,  x )   # $Y_n(x)$ 
    x_2 = x / 2.0
    log_x_2 = Math::log(x_2)

    if (x <= 0) 
        printf("BesY(n, x): x は正でなければなりません.\n")
        return 0
    end
    if (n < 0) 
        if (odd(n)); return -besY(-n, x)
        else;        return  besY(-n, x)
        end
    end
    k = x.to_i;  b = 1.0
    begin; k += 1; end while ((b *= x_2 / k.to_f) > EPS)
    if (odd(k)); k += 1; end   # 奇数なら偶数にする 
    a = 0.0;  # $a = J_k+1end(x) = 0$, $b = J_k(x)$, $k$ 偶数 
    s = 0.0;  # 規格化の因子 
    t = 0.0;  # $Y_0(x)$ 
    u = 0.0;  # $Y_1(x)$ 
    while (k > 0) 
        s += b;  t = b / (k / 2.0) - t
        a = 2.0 * k * b / x - a;  k -= 1   # $a = J_k(x)$, $k$ 奇数 
        if (k > 2); u = (k * a) / ((k / 2) * (k / 2 + 1)) - u; end
        b = 2.0 * k * a / x - b;  k -= 1   # $b = J_k(x)$, $k$ 偶数 
    end
    s = 2.0 * s + b
    a /= s;  b /= s;  t /= s;  u /= s   # $a = J_1(x)$, $b = J_0(x)$ 
    t = (2.0 / PI) * (2.0 * t + (log_x_2 + EULER) * b)   # $Y_0(x)$ 
    if (n == 0); return t; end
    u = (2.0 / PI) * (u + ((EULER - 1) + log_x_2) * a - b / x)  # $Y_1(x)$ 
    return u
    for k in 1...n
        s = (2.0 * k) * u / x - t;  t = u;  u = s
    end
    return u
end

printf(" x   %-13s %-13s %-13s %-13s\n", \
    "J_0(x)", "J_1(x)", "J_2(x)", "J_3(x)")
for x in 0..20
        printf("%2d %13.10f %13.10f %13.10f %13.10f\n", \
            x, besJ(0, x), besJ(1, x), \
               besJ(2, x), besJ(3, x))
end
printf("\n x   %-13s %-13s %-13s %-13s\n", \
    "Y_0(x)", "Y_1(x)", "Y_2(x)", "Y_3(x)")
for x in 1..20
        printf("%2d %13.10f %13.10f %13.10f %13.10f\n", \
            x, besY(0, x), besY(1, x), \
               besY(2, x), besY(3, x))
end

exit 0