multiplicity_eccentricity.py

#!/usr/bin/env python

from __future__ import print_function
import sys

import numpy as np
from scipy.ndimage.measurements import center_of_mass


dxy = 0.1


def main():
    for fn in sys.argv[1:]:
        S = np.loadtxt(fn)

        mult = dxy*dxy * S.sum()

        xymax = .5*(np.array(S.shape)-1) * dxy
        cm = np.array(center_of_mass(S))*dxy - xymax

        xm, ym = xymax
        X, Y = np.mgrid[-xm:xm+dxy:dxy, -ym:ym+dxy:dxy]
        xcm, ycm = cm
        X -= xcm
        Y -= ycm
        phi = np.arctan2(Y, X)

        R = np.hypot(X, Y)
        ecc = tuple(
            abs(np.average(np.exp(1j*n*phi), weights=S*R**n)) for n in (2, 3)
        )
        print(fn.split('/')[-1].rstrip('.dat'), mult, *ecc)


if __name__ == "__main__":
    main()