anisotropy_peak_fit.py

# d.d. 07/18/2016
__author__ = "Guen Prawiroatmodjo <http://github.com/guenp>"
import peakutils
from pylab import *
import colorgraphs as g


def get_matrix(d,series,sweeplen):
    return series[:len(d)-(len(d)%sweeplen)].reshape(np.floor(len(d))/sweeplen,sweeplen)


def smooth(x,window_len=11,window='hanning'):
    '''
    smooth the data using a window with requested size.
    '''
    if x.ndim != 1:
        raise ValueError("smooth only accepts 1 dimension arrays.")
    if x.size < window_len:
        raise ValueError("Input vector needs to be bigger than window size.")
    if window_len<3:
        return x
        
    if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
        raise ValueError("Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'")
    
    s=np.r_[x[window_len-1:0:-1],x,x[-1:-window_len:-1]]

    if window == 'flat': #moving average
        w=np.ones(window_len,'d')
    else:
        w=eval('np.'+window+'(window_len)')

    y=np.convolve(w/w.sum(),s,mode='valid')
    return y[(round(window_len/2)-1):-round(window_len/2)]


def get_Bp(B,tg,cond,xmin=.425,xmax=.46,bmin=.4,bmax=.7,smoothfactor=20,upper=False):
    figure(figsize=(4,12))
    c1_vec = []
    c2_vec = []
    b1_vec = []
    b2_vec = []
    err_vec = []
    n_curves = arange(0,len(cond))
    n_curves = [val for pair in zip(n_curves[:np.ceil(len(cond)/2)],n_curves[np.ceil(len(cond)/2):len(cond)][::-1]) for val in pair]
    for n in n_curves:
        B_ = B[n].mean()
        x,y = tg[n][(tg[0]>xmin)&(tg[0]<xmax)],cond[n][(tg[0]>xmin)&(tg[0]<xmax)]
        xmax_ = x[len(y)-10:][y[len(y)-10:].argmin()]
        xmin_ = x[:10][y[:10].argmin()]
        x,y = tg[n][(tg[0]>xmin_)&(tg[0]<xmax_)],cond[n][(tg[0]>xmin_)&(tg[0]<xmax_)]
        peak_indices = peakutils.indexes(y-smooth(y,smoothfactor))
        plot(x,y+n*1e-6,'b')
        plot(x[peak_indices],y[peak_indices]+n*1e-6,'ok')
        if len(peak_indices)>0:
            c1_vec.append(x[peak_indices][0])
            b1_vec.append(B_)
        if len(peak_indices)>1:
            c2_vec.append(x[peak_indices][1])
            b2_vec.append(B_)

    figure()
    b_vec = array(b1_vec+b2_vec)
    c_vec = array(c1_vec+c2_vec)
    xnew = arange(-bmax,bmax,.01)
    plot(b_vec, c_vec, '.')
    
    x,y=b_vec[abs(b_vec)<bmin],c_vec[abs(b_vec)<bmin]
    plot(x,y,'c.')
    y_mid = mean(y)*ones(len(xnew))
    plot(xnew, y_mid,'r')
    
    if upper:
        x,y=b_vec[(c_vec>y_mid.mean())&(abs(b_vec)>bmin)&(abs(b_vec)<bmax)],c_vec[(c_vec>y_mid.mean())&(abs(b_vec)>bmin)&(abs(b_vec)<bmax)]
    else:
        x,y=b_vec[(c_vec<y_mid.mean())&(abs(b_vec)>bmin)&(abs(b_vec)<bmax)],c_vec[(c_vec<y_mid.mean())&(abs(b_vec)>bmin)&(abs(b_vec)<bmax)]
    x_,y_ = x[x<0],y[x<0]
    if len(x_)>0:
        plot(x_,y_,'m.')
        p = polyfit(x_,y_,1)
        y_neg = polyval(p,xnew)
        plot(xnew, y_neg,'r')
        Bp_neg = xnew[abs(y_neg-y_mid).argmin()]
    else:
        Bp_neg = b_vec[abs(b_vec)<bmin].min()
        p = [0,0]

    x_,y_ = x[x>0],y[x>0]
    if len(x_)>0:
        plot(x_,y_,'m.')
        q = polyfit(x_,y_,1)
        y_pos = polyval(q,xnew)
        plot(xnew, y_pos,'r')
        Bp_pos = xnew[abs(y_pos-y_mid).argmin()]
    else:
        Bp_pos = b_vec[abs(b_vec)<bmin].max()
        q = [0,0]
        
    return Bp_pos, Bp_neg, p[0], q[0]


d = g.get('20150807_003115')
s = d.meta['shape']

Bx,By,Bz,B,tg,vsd,cond,v1,v3 = (get_matrix(d,s,len(d.topgate.unique())) for s in [d.Bperp, d.Bpary, d.Bpar, sign(d.Bperp)*np.sqrt(d.Bpar.round(3)**2 + d.Bperp.round(3)**2 + d.Bpary.round(3)**2),d.topgate, d.dc_source_voltage,d.X_current/d.X_voltage2, d.X_voltage1, d.X_voltage3])
Bp1_pos, Bp1_neg, p1, q1 = get_Bp(B,tg,cond,xmin=.42,xmax=.46,bmin=.4,bmax=.85)
Bp2_pos, Bp2_neg, p2, q2 = get_Bp(B,tg,cond,xmin=.675,xmax=.705,bmin=.35,bmax=.95)