forgot atarilib and phasetable

#!/usr/bin/env python

"PHASE TABLE"

import time

import numpy as np

import h5py

import matplotlib.pyplot as plt

from scipy.interpolate import interp2d

from pysen.atarilib import fit_echo, Spectrum

def phase_table(hdfile, iecho=None, verbose=False, **kwargs):

    "phase table"

    tbin1   = kwargs.pop('tbin1', 4)      # TOF bins

    tbin2   = kwargs.pop('tbin2',38)

    # FIXME: HARDCODED

    xpix1   = kwargs.pop('xpix1', 10)

    xpix2   = kwargs.pop('xpix2', 22)

    ypix1   = kwargs.pop('ypix1', 10)

    ypix2   = kwargs.pop('ypix2', 22)

    #ntaus   =  hdfile.attrs['scan_length']

    n_idx = { 'dn': hdfile['/'].attrs['point_to_down'],

              'up': hdfile['/'].attrs['point_to_up'] ,

              'nphases': hdfile['/'].attrs['no_of_phases'] }

    nt      = hdfile['/detector'].attrs['no_t_channels']

    ny      = hdfile['/detector'].attrs['no_y_channels']

    nx      = hdfile['/detector'].attrs['no_x_channels']

    nph     = n_idx['dn']

    results = []

    for echo in list(hdfile['/data'].values()):

        if iecho is not None and echo.attrs['id'] != iecho:

            continue

        phase   = echo['phase']

        physics = echo['phys']

        params  = echo['params']

        tech    = echo['tech']

        #tau0    = physics.attrs['fouriertime']/NANOSECOND

        #j0      = physics.attrs['prim_fieldintegral']

        lam0    = physics.attrs['lambda']

        phase0  = float(tech.attrs['i5'][0])

        maincur = float(tech.attrs['i00'][0])

        phi_angle = float(params.attrs['phi'][0])

        lmax    = np.array(params['lambdaTable']).flatten()

        dlam    = np.ones_like(lmax)*(lmax[1]-lmax[0])

        phasesens = float(params['phaseangle'].attrs['sensitivities'][0])

        det  = phase['detector'][...]

        #moni = phase['counter' ][...]

        pcha = phase['proton_charge'][...]

        cur  = phase['phase_current'][:, 0] # actual value

        pcha0 = np.average(pcha)

        pcha = np.tile(pcha, (ny, nx, nt)).reshape(n_idx['nphases'],ny,nx,nt)

        pcha = pcha/float(pcha0)

        det  = det/pcha

        #wlen = np.sum(det[:,xpix1:xpix2,ypix1:ypix2,tbin1:tbin2], axis=(0,1,2))

        wlen = np.sum(det, axis=(0,1,2))

        pha  = np.sum(det[:nph,xpix1:xpix2,ypix1:ypix2,tbin1:tbin2], axis=(1,2,3))

        epha = np.sqrt(np.sum(det[:nph,xpix1:xpix2,ypix1:ypix2,tbin1:tbin2]**2, axis=(1,2,3)))

        #phase_pf, xpf, ypf = fit_poly(cur[:nph], (pha, epha), iphase,deg=4)

        phase_ef, _, _ = fit_echo(cur[:nph], (pha, epha), Spectrum(lam=lmax, dlam=dlam, flux=wlen),

                                  lam0=lam0, phase0=phase0, phasesens=phasesens)

        if verbose:

            print("%7.4f %7.4f %7.4f %7.4f " % (maincur, phi_angle, phase_ef, phase0))

        results.append((maincur, phi_angle, phase_ef, phase0))

    return np.asarray(results)

def main():

    "the main"

    import argparse

    parser = argparse.ArgumentParser(description='atari plot')

    parser.add_argument('file', metavar='filename', help='file to process', nargs='+')

    #parser.add_argument('--echo',  '-e',  dest='echo', type=int,   default=None,

    #                    help='set echo to display (default=%(default)s)')

    parser.add_argument('--pos', '-p', dest='pos', type=int, default=None,

                        help='set instrument position  (default=%(default)s)')

    parser.add_argument('--verbose', '-v', dest='verbose', action='store_true', default=False,

                        help='be verbose')

    args = parser.parse_args()

    phtab = None

    for filename in args.file:

        with h5py.File(filename, 'r') as hdf5file:

            res = phase_table(hdf5file, verbose=args.verbose)

            if phtab is None:

                phtab = res

            else:

                phtab = np.vstack((phtab,res))

    nx = 23

    ny = 11

    i00   = phtab[:,0]

    phi   = phtab[:,1]

    phase = phtab[:,2]

    # FIXME: hard coded

    xi = np.linspace(0.0, 88.0,nx)

    maxy = 33.3

    yi = np.linspace(-maxy/9, maxy,ny)

    x,y = np.meshgrid(xi,yi)

    f = interp2d(i00,phi,phase, kind='linear') # bounds_error=False, fill_value=0)

    z = f(xi,yi)

    if args.pos:

        print("input phasetable")

        print("phasetable pos2 rebuild")

        print("c ", time.asctime())

        print("*\t angles ", end='')

        for j in range(ny):

            print("%7.4f " % yi[j], end='')

        print("\nc amperes(main) --- phase current(i5)")

        for i in range(nx):

            print("*\t%7.3f " % xi[i],end='')

            for j in range(ny):

                print("%7.4f " % z[j,i], end='')

            print()

        print("eof")

    plt.figure()

    plt.contourf(x,y,z,21)

    plt.colorbar()

    plt.show()

if __name__ == "__main__":

    main()

#!/usr/bin/env python

"ATARI Plot utils"

from functools import partial

from collections import namedtuple

import numpy as np

from pysen import OMEGA_N

from .echo.fit import flux_weighted_eshape, echo_fit, echo_curve

Spectrum = namedtuple('Spectrum', ['lam','dlam','flux'])

def fit_poly(curr, counts, iphase, deg=4, dix=3):

    "fit poly for phasetables and atari plot"

    cnts, _ = counts

    imx = iphase or np.argmax(cnts)

    imx1 = max(0        , imx-dix)

    imx2 = min(len(cnts), imx+dix)

    xf = curr[imx1:imx2+1]

    yf = cnts[imx1:imx2+1]

    pf = np.polyfit(xf, yf, deg=deg)

    xf = np.linspace(min(xf), max(xf), 101)

    yf = np.polyval(pf, xf)

    imx = np.argmax(yf)

    return xf[imx], xf,yf

def fit_echo(current, counts, spectrum, lam0, phase0, phasesens=1.0, ):

    "fit echo for atari plot and phase tables"

    # get maximum estimate

    phase_pf, _, _ = fit_poly(current, counts, None)

    # convert currents to dJ

    dj  = np.radians(phasesens*(current  - phase0))/OMEGA_N/lam0

    dj0 = np.radians(phasesens*(phase_pf - phase0))/OMEGA_N/lam0

    # FIXME: hardcoded

    ddj  = dj[3]-dj[0]

    sfun = partial(flux_weighted_eshape,

                   avlam=spectrum.lam, delam=spectrum.dlam, flux=spectrum.flux)

    res  = echo_fit(sfun, dj, counts[0], counts[1], dj0=dj0, djlim=(dj0-ddj,dj0+ddj))

    phase_ef = np.degrees(res['dj0'][0]*OMEGA_N*lam0)/phasesens+phase0

    djf, efit = echo_curve(sfun, dj, dj0=res['dj0'][0], amplitude=res['amplitude'][0],

                           average=res['average'][0], npoints=101)

    # convert back to current

    curf = np.degrees(djf*OMEGA_N*lam0)/phasesens+phase0

    return phase_ef, curf, efit