常见的Abel逆变换求解方法及其python实现

上一篇论文中总结了一些常见的Abel逆变换求解方法Abel逆变换及其求解方法，应粉丝要求，这篇文章写一下对应的算法实现。

先简单回顾一下，

1. Abel变换（Abel transform）是一种积分变换，常用于球对称或轴对称函数的分析，其数学表达式为：

    （1）

2.Abel逆变换（ inverse Abel transform），其数学表达式为：

。（2）

求解方法

1.直接离散法

。 （3）

代码实现：

def IAbel_DD(Fy, N, DeltaY, DeltaR):
    r"""
    Inverse Abel Transform: direct-discretization，DD
    Args:
        | Fy: The integrated intensities Fy
        | N: N+1 is the number of discrete distribution datapoint of Fy
        | DeltaY: spatial separation between data
        | DeltaR: spatial separation between data, DeltaY = DeltaR
    Returns:
        | frDD: fr at different radius, solved by DD method
    """

    frDD = np.zeros(N + 1)
    r = np.zeros(N + 1)
    for j in range(0, N + 1):
        r[j] = DeltaR * j
        frDD[j] = 0
        y = np.zeros(N)
        for i in range(j, N):
            y[i] = DeltaY * i
            Deno = np.sqrt(y[i] + DeltaY/2)**2 - r[j]**2)
            Nume = Fy[i+1] - Fy[i]
            Quo = -Nume/Deno/math.pi
            frDD[j] = frDD[j] + Quo
    return frDD

这个离散法很简单，上面的代码我直接写的，没测试，不保证完全正确哈。

 2.Hankel-Fourier变换法：

详细推导见上一篇文章，结果是

    （4）

其中，J0是0阶Bessel function，可以直接用：

J0 = jv(0, j * k / (2 * N + 1)

其余就是很简单的离散求和。

3.Nestor-Olsen法：

该算法由Nestor和Olsen提出，f(r)与F(r)之间得关系由下式给出：

                     (5)

其中，Bj,i满足：

                  （6）

其中，。    （7）

代码实现：

def IAbel_NO(Fy, N, DeltaY):
    """
    Inverse Abel Transform: Nestor–Olsen method (NO method)
    :param Iy: The integrated intensities Fy
    :param N: N + 1 is the number of discrete distribution datapoint of Iy
    :param DeltaY: spatial separation between data

    :return: frNO: The fr at different radius, by NO method
    """
    Ce1 = -2 / math.pi / DeltaY
    frNO = np.zeros(N + 1)
    for j in range(0, N+1):
        sum = 0
        for i in range(j, N):
            if i == j:
                Aji = AJI(i, j)
                Bji = -Aji
            else:
                Aji = AJI(i, j)
                Aji1 = AJI(i - 1, j)
                Bji = Aji1 - Aji
            sum = sum + Iy[i] * Bji
        frNO[j] = Ce1 * sum
    return frNO

 其中AJI我是写了一个函数，就按照式（7）来写就可以，这里我就不贴上来了。