Commit 5d4e97b9 authored by Zolnierczuk, Piotr's avatar Zolnierczuk, Piotr
Browse files

atariplot fixes and improvements

parent 4adb3d23

misc/phasetable.py

+4 −3
Original line number Diff line number Diff line
@@ -8,7 +8,7 @@ import matplotlib.pyplot as plt 


from scipy.interpolate import interp2d



from pysen.atarilib import fit_echo, Spectrum

from pysen.atarilib import fit_echo_current, Spectrum



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

    "phase table"
@@ -68,7 +68,8 @@ def phase_table(hdfile, iecho=None, verbose=False, **kwargs): 
        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),

        phase_ef, _, _ = fit_echo_current(cur[:nph], (pha, epha),

                                          Spectrum(lam=lmax, dlam=dlam, flux=wlen),

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



        if verbose:

pysen/atarilib.py

+5 −5
Original line number Diff line number Diff line
@@ -10,7 +10,7 @@ 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):

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

    "fit poly for phasetables and atari plot"

    cnts, _ = counts

    imx = iphase or np.argmax(cnts)
@@ -22,12 +22,12 @@ def fit_poly(curr, counts, iphase, deg=4, dix=3): 
    xf = np.linspace(min(xf), max(xf), 101)

    yf = np.polyval(pf, xf)

    imx = np.argmax(yf)

    return xf[imx], xf,yf

    return xf[imx], (xf,yf)



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

def fit_echo_current(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)

    phase_pf, _ = fit_poly_current(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
@@ -41,4 +41,4 @@ def fit_echo(current, counts, spectrum, lam0, phase0, phasesens=1.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
 
    return phase_ef, (curf, efit), res

pysen/revision.py

+2 −2
Original line number Diff line number Diff line
@@ -2,9 +2,9 @@ 
PySEN revision module

"""

import sys

__version__  = "0.5.4"

__version__  = "0.5.5"

__release__  = "dev1"

__date__     = "Oct 26, 2018"

__date__     = "Jan 24, 2019"

#

VERSION  = __version__

RELEASE  = __release__

scripts/atariplot.py

+12 −3
Original line number Diff line number Diff line
@@ -8,7 +8,7 @@ import matplotlib.pyplot as plt 




from pysen import ANGSTROM, NANOSECOND

from pysen.atarilib import fit_echo, Spectrum

from pysen.atarilib import fit_echo_current, Spectrum





def atari_plot(hdfile, iecho=0, **kwargs):
@@ -95,8 +95,14 @@ def atari_plot(hdfile, iecho=0, **kwargs): 
        pha  = np.sum(pha[:nph], axis=(1,2,3))



        spectrum =  Spectrum(lam=lmax, dlam=dlam, flux=wlen)

        phase_ef, xef, yef = fit_echo(cur[:nph], (pha, epha), spectrum,

        phase_ef, (xef, yef), res = fit_echo_current(cur[:nph], (pha, epha), spectrum,

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

        res['up'] = np.average(up),np.sqrt(np.sum(up)/len(up))

        res['dn'] = np.average(dn),np.sqrt(np.sum(dn)/len(dn))

        #for key in res:

        #    print(key, res[key])



        R = 2*res['amplitude'][0]/( res['up'][0] - res['dn'][0] )



        p = plt.errorbar(cur[:nph], pha, yerr=epha, fmt='.', lw=1)

        dcol = p[0].get_color()
@@ -107,6 +113,9 @@ def atari_plot(hdfile, iecho=0, **kwargs): 
        plt.axvline(phase_ef, color='red',   lw=2, ls='-',  label='%.3fA' % phase_ef)

        plt.axhline(np.average(up))

        plt.axhline(np.average(dn))

        ax = plt.gca()

        plt.text(0.18, 0.10, r'$R(Q,t)$=%.3g' % R,

                 horizontalalignment='center', verticalalignment='center', transform=ax.transAxes)

        plt.xlabel(r'phase current [A]')

        plt.legend(loc='best')

        plt.grid(True)