PyAbel test @MikhailRyazano anisotropy code for real data

beta-test.py

import numpy as np
import abel
import bz2
import anisotro

from scipy.signal import find_peaks, peak_widths
import matplotlib.pylab as plt

# O2- PES VM-image
imagefile = bz2.BZ2File('O2-ANU1024.txt.bz2')
IM = np.loadtxt(imagefile)
IMc = abel.tools.center.center_image(IM, center='slice', square=True)


# Abel transform - linbasex ------------------------------------------
legendre_orders = [0, 2]  # Legendre polynomial orders
proj_angles = np.arange(0, np.pi, np.pi/10)  # p000rojection angles in 10 degree steps
radial_step = 1  # pixel grid
smoothing = 0.1  # smoothing 1/e-width for Gaussian convolution smoothing
threshold = 0.0  # threshold for normalization of higher order Newton spheres
clip = 0  # clip first vectors (smallest Newton spheres) to avoid singularities

LIM = abel.Transform(IMc, method='linbasex',
                     transform_options=dict(basis_dir='./',
                     proj_angles=proj_angles, radial_step=radial_step,
                     smoothing=smoothing, threshold=threshold, clip=clip,
                     return_Beta=True, verbose=True))

# basex -------------------
BIM = abel.Transform(IMc, method="basex", transform_options=dict(reg=1000),
                     symmetry_axis=None, angular_integration=True)

# speed distribution (basex)
radial, speed = BIM.angular_integration

speed /= speed[200:].max()  # normalize, exclude noise near centre

# find the photoelectron transition peaks -----------------
peaks, _ = find_peaks(speed, height=0.24)
widths = peak_widths(speed, peaks, rel_height=0.5)

# set radial range across peaks
r_range = []
for r, w in zip(radial[peaks], widths[0]):
    r_range.append((r-w/2, r+w/2))

# anisotropy parameter from image for each tuple radial r_range ----------
Beta, Amp, Rmid, Ivstheta, theta =\
              abel.tools.vmi.radial_integration(BIM.transform, r_range)

# cf anisotro
q = abel.tools.symmetry.get_image_quadrants(IMc)
aq = abel.basex.basex_transform(q[0])
I, beta = anisotro.Ibeta(aq[::-1], method='remap')  # flip to match example

# plots ----------------------------
plt.plot(radial, speed, lw=1, label='basex')

radial = radial[:len(beta)]   #  PyAbel extends radius to the image corner

# plot beta variation for (basex, linbasex) only in regions of peaks
for i, r in enumerate(r_range):
    subM = np.logical_and(radial > r[0], radial < r[1])
    subL = np.logical_and(LIM.radial > r[0], LIM.radial < r[1])
    if i == 0:
        plt.plot(radial[subM], beta[subM], 'C1', lw=1,  label='Mikhail')
        plt.plot(LIM.radial[subL], LIM.Beta[1][subL], 'C2', lw=1,
                 label='linbasex')
    else:
        plt.plot(radial[subM], beta[subM], 'C1', lw=1)
        plt.plot(LIM.radial[subL], LIM.Beta[1][subL], 'C2', lw=1)


# betas from angular integration
BetaT = np.transpose(Beta)
plt.errorbar(Rmid, BetaT[0], BetaT[1], fmt='oC3', label='integrated radii')

plt.xlabel('radius (pixels)')
plt.ylabel(r'$\beta$                    intensity (arb. u.)')
plt.axis(xmin=300, xmax=420, ymin=-1, ymax=1)
plt.legend(labelspacing=0.1, fontsize='small')
plt.title ('O2-ANU1024.txt')

plt.savefig('beta-test.png', dpi=75, bb_inches='tight')
plt.show()