PyAbel transform pairs test

TransformPairs-inverse-diff-1st.py

import numpy as np
import abel
import matplotlib.pyplot as plt

n = 101

fig, axs = plt.subplots(3, 2, sharex=True, sharey=True)

case = '1st-order-diff'
# case = 'gradient'
print(case)
print('profile  hansenlaw   basex    3pt')
for i, prof in enumerate([1, 2, 3, 5, 6, 7]):
    f = abel.tools.analytical.TransformPair(n=n, profile=prof)

    hl = abel.hansenlaw.hansenlaw_transform(f.abel, dr=f.dr, hold_order=0)
    tp = abel.dasch.three_point_transform(f.abel, dr=f.dr)
    ba = abel.basex.basex_transform(f.abel, dr=f.dr, verbose=False)

    row = i // 2
    col = i % 2 
    ax = axs[row, col] 

    ax.plot(f.r, (hl-f.func)*100, lw=1)
    ax.plot(f.r, (ba-f.func)*100, lw=1)
    ax.plot(f.r, (tp-f.func)*100, lw=1)

    ax.annotate(f.label, (0.75, -0.4), fontsize='smaller')
    ax.axis(ymin=-0.5, ymax=0.5)

    msehl = np.square(hl-f.func).mean()
    mseba = np.square(ba-f.func).mean()
    msetp = np.square(tp-f.func).mean()
    print(f'{f.label[-1]:8s} {msehl:4.2e}  {mseba:4.2e}  {msetp:4.2e}')

    ax.annotate(f'h&l {msehl:4.2e}', (0.1, 0.4), fontsize=7, color='C0')
    ax.annotate(f'bas {mseba:4.2e}', (0.1, 0.3), fontsize=7, color='C1')
    ax.annotate(f'3pt {msetp:4.2e}', (0.1, 0.2), fontsize=7, color='C2')
    if col == 0:
        ax.set_ylabel(r'diff. (%)')
    if row == 3:
        ax.set_xlabel('radius')

plt.subplots_adjust(wspace=0, hspace=0)
plt.suptitle(f'transform pairs, n = {n:d}, {case}')
plt.savefig(f'TP-{case}.png', dpi=100, bb_inches='tight')
plt.show()