def real(a): # doubled ranks!

    b = tensor()
    
    b.n = a.n
    b.u[0] = np.concatenate((np.real(a.u[0]), np.imag(a.u[0])), 1)
    b.u[1] = np.concatenate((np.real(a.u[1]), np.imag(a.u[1])), 1)
    b.u[2] = np.concatenate((np.real(a.u[2]), np.imag(a.u[2])), 1)

    R1 = np.zeros((2*a.r[0], a.r[0]), dtype = np.complex128)
    R2 = np.zeros((2*a.r[1], a.r[1]), dtype = np.complex128)
    R3 = np.zeros((2*a.r[2], a.r[2]), dtype = np.complex128)

    R1[:a.r[0], :] = np.identity(a.r[0])
    R1[a.r[0]:, :] = 1j*np.identity(a.r[0])
    R2[:a.r[1], :] = np.identity(a.r[1])
    R2[a.r[1]:, :] = 1j*np.identity(a.r[1])
    R3[:a.r[2], :] = np.identity(a.r[2])
    R3[a.r[2]:, :] = 1j*np.identity(a.r[2])
    
    
    GG = np.tensordot(np.transpose(a.core,[2,1,0]),np.transpose(R1), (2,0))
    GG = np.tensordot(np.transpose(GG,[0,2,1]),np.transpose(R2), (2,0))
    GG = np.transpose(GG,[1,2,0])
    b.core = np.real(np.tensordot(GG,np.transpose(R3), (2,0)))

    b.r = b.core.shape
    
    return b

def fourier(self):
    """
    Generate a profile of fourier coefficients, amplitudes and phases
    """
    if pynbody.config['verbose'] : print 'Profile: fourier()'

    f = {'c': np.zeros((7, self.nbins),dtype=complex),
         'amp': np.zeros((7, self.nbins)),
         'phi': np.zeros((7, self.nbins))}

    for i in range(self.nbins):
        if self._profiles['n'][i] > 100:
            phi = np.arctan2(self.sim['y'][self.binind[i]], self.sim['x'][self.binind[i]])
            mass = self.sim['mass'][self.binind[i]]
            
            hist, binphi = np.histogram(phi,weights=mass,bins=100)
            binphi = .5*(binphi[1:]+binphi[:-1])
            for m in range(7) : 
                f['c'][m,i] = np.sum(hist*np.exp(-1j*m*binphi))


    f['c'][:,self['mass']>0] /= self['mass'][self['mass']>0]
    f['amp'] = np.sqrt(np.imag(f['c'])**2 + np.real(f['c'])**2)
    f['phi'] = np.arctan2(np.imag(f['c']), np.real(f['c']))

    return f

    def process(self, X, V, C):
        BifPoint.process(self, X, V, C)

        J_coords = C.sysfunc.jac(X, C.coords)

        eigs, VL, VR = linalg.eig(J_coords, left=1, right=1)

        # Check for nonreal multipliers
        found = False
        for i in range(len(eigs)):
            for j in range(i+1,len(eigs)):
                if abs(imag(eigs[i])) > 1e-10 and \
                   abs(imag(eigs[j])) > 1e-10 and \
                   abs(eigs[i]*eigs[j] - 1) < 1e-5:
                    found = True

        if not found:
            del self.found[-1]
            return False

        self.found[-1].eigs = eigs

        self.info(C, -1)

        return True

def writeToDatFile(data, file, field="E", punits="mm", funits="V/m", pscale=1.0, fscale=1.0):
    """
    Write to field map to DAT data file.
    """

    nx = data.nx
    ny = data.ny
    nz = data.nz

    file.write("       x [{0}]         y [{0}]         z [{0}]     {1}xRe [{2}]     {1}yRe [{2}]     {1}zRe [{2}]     {1}xIm [{2}]     {1}yIm [{2}]     {1}zIm [{2}]  \r\n".format(punits, field, funits))
    file.write("------------------------------------------------------------------------------------------------------------------------------------------\r\n")

    for x in xrange(nx):
        px = pscale * data.px[x]
        for y in xrange(ny):
            py = pscale * data.py[y]
            for z in xrange(nz):
                pz = pscale * data.pz[z]
                fxre = fscale * numpy.real(data.fx[x,y,z])
                fyre = fscale * numpy.real(data.fy[x,y,z])
                fzre = fscale * numpy.real(data.fz[x,y,z])
                fxim = fscale * numpy.imag(data.fx[x,y,z])
                fyim = fscale * numpy.imag(data.fy[x,y,z])
                fzim = fscale * numpy.imag(data.fz[x,y,z])
                file.write("%13.1f  %13.1f  %13.1f  %13.6g  %13.6g  %13.6g  %13.6g  %13.6g  %13.6g\r\n" % (px, py, pz, fxre, fyre, fzre, fxim, fyim, fzim))

    def start(self, f, a, b, args=()):
        r"""Prepare for the iterations."""
        self.function_calls = 0
        self.iterations = 0

        self.f = f
        self.args = args
        self.ab[:] = [a, b]
        if not np.isfinite(a) or np.imag(a) != 0:
            raise ValueError("Invalid x value: %s " % (a))
        if not np.isfinite(b) or np.imag(b) != 0:
            raise ValueError("Invalid x value: %s " % (b))

        fa = self._callf(a)
        if not np.isfinite(fa) or np.imag(fa) != 0:
            raise ValueError("Invalid function value: f(%f) -> %s " % (a, fa))
        if fa == 0:
            return _ECONVERGED, a
        fb = self._callf(b)
        if not np.isfinite(fb) or np.imag(fb) != 0:
            raise ValueError("Invalid function value: f(%f) -> %s " % (b, fb))
        if fb == 0:
            return _ECONVERGED, b

        if np.sign(fb) * np.sign(fa) > 0:
            raise ValueError("a, b must bracket a root f(%e)=%e, f(%e)=%e " %
                             (a, fa, b, fb))
        self.fab[:] = [fa, fb]

        return _EINPROGRESS, sum(self.ab) / 2.0

def dynamic_axis(zero,pole,K):
    if list(zero)+list(pole)==[]:
        x_min,x_max,y_min,y_max = -1,1,-1,1
    else:
        x_min = min(list(np.real(zero))+list(np.real(pole)))
        x_max = max(list(np.real(zero))+list(np.real(pole)))
        if x_min == x_max:
            x_min -= 1
            x_max += 1
        else:
            x_min,x_max = (x_min- 0.75*(x_max-x_min)),(x_max + 0.75*(x_max-x_min))
        y_min = min(list(np.imag(zero))+list(np.imag(pole)))
        y_max = max(list(np.imag(zero))+list(np.imag(pole)))
        if y_min == y_max:
            y_min -= 1
            y_max += 1
        else:
            y_min,y_max = (y_min- 0.75*(y_max-y_min)),(y_max + 0.75*(y_max-y_min))
    if K == 0:
        z_min,z_max = -1,1
    else:
        K = abs(K)
        print K
        z_min,z_max = K*0.1,K*10.0
    return x_min,x_max,y_min,y_max,z_min,z_max

def instantaneous_frequency(data, fs, fk):
    """
    Instantaneous frequency of a signal.

    Computes the instantaneous frequency of the given data which can be
    windowed or not. The instantaneous frequency is determined by the time
    derivative of the analytic signal of the input data.

    :type data: :class:`~numpy.ndarray`
    :param data: Data to determine instantaneous frequency of.
    :param fs: Sampling frequency.
    :param fk: Coefficients for calculating time derivatives
        (calculated via central difference).
    :return: **omega[, domega]** - Instantaneous frequency of input data, Time
        derivative of instantaneous frequency (windowed only).
    """
    x = envelope(data)
    if len(x[0].shape) > 1:
        omega = np.zeros(x[0].shape[0], dtype=np.float64)
        i = 0
        for row in x[0]:
            f = np.real(row)
            h = np.imag(row)
            # faster alternative to calculate f_add
            f_add = np.hstack(([f[0]] * (np.size(fk) // 2), f, [f[np.size(f) - 1]] * (np.size(fk) // 2)))
            fd = signal.lfilter(fk, 1, f_add)
            # correct start and end values of time derivative
            fd = fd[np.size(fk) - 1 : np.size(fd)]
            # faster alternative to calculate h_add
            h_add = np.hstack(([h[0]] * (np.size(fk) // 2), h, [h[np.size(h) - 1]] * (np.size(fk) // 2)))
            hd = signal.lfilter(fk, 1, h_add)
            # correct start and end values of time derivative
            hd = hd[np.size(fk) - 1 : np.size(hd)]
            omega_win = abs(((f * hd - fd * h) / (f * f + h * h)) * fs / 2 / np.pi)
            omega[i] = np.median(omega_win)
            i = i + 1
        # faster alternative to calculate omega_add
        omega_add = np.hstack(
            ([omega[0]] * (np.size(fk) // 2), omega, [omega[np.size(omega) - 1]] * (np.size(fk) // 2))
        )
        domega = signal.lfilter(fk, 1, omega_add)
        # correct start and end values of time derivative
        domega = domega[np.size(fk) - 1 : np.size(domega)]
        return omega, domega
    else:
        omega = np.zeros(np.size(x[0]), dtype=np.float64)
        f = np.real(x[0])
        h = np.imag(x[0])
        # faster alternative to calculate f_add
        f_add = np.hstack(([f[0]] * (np.size(fk) // 2), f, [f[np.size(f) - 1]] * (np.size(fk) // 2)))
        fd = signal.lfilter(fk, 1, f_add)
        # correct start and end values of time derivative
        fd = fd[np.size(fk) - 1 : np.size(fd)]
        # faster alternative to calculate h_add
        h_add = np.hstack(([h[0]] * (np.size(fk) // 2), h, [h[np.size(h) - 1]] * (np.size(fk) // 2)))
        hd = signal.lfilter(fk, 1, h_add)
        # correct start and end values of time derivative
        hd = hd[np.size(fk) - 1 : np.size(hd)]
        omega = abs(((f * hd - fd * h) / (f * f + h * h)) * fs / 2 / np.pi)
        return omega

def mode_finder(xxx_todo_changeme):
    (det_mat_slice, wl_list, kx) = xxx_todo_changeme
    dispcurve = []
    xtol = 1e-4*wl_0
    for j in range(len(wl_list)-1):
        # Check determinant crosses zero, noth real and imaginary
        if np.real(det_mat_slice[j])*np.real(det_mat_slice[j+1]) < 0:
            if np.imag(det_mat_slice[j])*np.imag(det_mat_slice[j+1]) < 0:
                if j != 0:
                    diffreq = np.abs(det_mat_slice[j-1]-det_mat_slice[j])
                else:
                    diffreq = np.abs(det_mat_slice[j+2]-det_mat_slice[j+1])
                # Check we are not just at a discontinuity
                if np.abs(det_mat_slice[j+1]-det_mat_slice[j]) < 3*diffreq:
                    try:
                        # Optimise the wl
                        finwl = optimize.brentq(lambda wl: np.real(np.exp(1j)*simulate_stack([wl, kx])),wl_list[j],wl_list[j+1],rtol=1e-3,xtol=xtol)
                        findet=simulate_stack([finwl, kx])
                        print('#################################')
                        print('found root = ', findet)
                        print('#################################')
                        #check the final determinant is below some tolerance
                        if np.abs(findet)< 1.e-3:
                            finfreq=2*np.pi*c_speed*1e9/finwl
                            dispcurve.append((kx*1e9,finfreq))
                    except AttributeError:
                        print(det_mat_slice[j], det_mat_slice[j+1])
    return dispcurve

def write_readable(sigma, g_w, n_k, n_w, wgrid):

    print "Writing E_file for checking"

    print n_k, n_w
    print sigma[n_w-1][n_k-1][n_k-1]


    greal_file=open("PARSED_Greal_1", "w+")
    gimag_file=open("PARSED_Gimag_1", "w+")
    sigma_file=open("PARSED_SIGMA_1", "w+")
    for w in range (n_w/2, n_w):
      print w
      print np.real(g_w[w][0][0])
      write_string=str(wgrid[w])+" "
      write_string2=str(wgrid[w])+" "
      write_string_sigma=str(wgrid[w])+" "
      for kx in range (n_k):
        for ky in range (n_k):
          write_string+=" "+str(np.real(g_w[w][kx][ky])[0])
          write_string2+=" "+" "+str(np.imag(g_w[w][kx][ky])[0])
          write_string_sigma+=" "+str(np.real(sigma[w][kx][ky])[0])+" "+str(np.imag(sigma[w][kx][ky])[0])
     # print write_string
      greal_file.write(write_string+"\n")
      gimag_file.write(write_string+"\n")
      sigma_file.write(write_string_sigma+"\n")

def extract_outfield_from_dict( outpdict ):
    """ 
    extract only the output field (time, freq)
    and the freq vectors from a dict created
    by mydict  = loadoutput(filename)
   
    returns a Nx7 numpy.array
    """    

    tvec  = outpdict['tvec']
    omvec = outpdict['omvec']
    relomvec = outpdict['relomvec']
    tfieldreal = np.real(outpdict['tfield2'])
    tfieldimag = np.imag(outpdict['tfield2'])
    ffieldreal = np.real(outpdict['ffield2'])
    ffieldimag = np.imag(outpdict['ffield2'])
    M = np.zeros( [len( tvec), 7])
    M[:,0]=tvec
    M[:,1]=omvec
    M[:,2]=relomvec
    M[:,3]=tfieldreal
    M[:,4]=tfieldimag
    M[:,5]=ffieldreal
    M[:,6]=ffieldimag
    return M

 def resData2(self,xval,yval,yval_sim=None,xlabel='Frequency',xunit='Hz',plottype='amp',
             ampformat='log',save=False,Dir=None):
     
     self.figure.clear()
     self.figure.subplots_adjust(bottom=0.15,left=0.17)
     self.axes = self.figure.add_subplot(111)
             
     if plottype=='real/imag':
         self.axes.plot(np.real(yval),np.imag(yval),np.real(yval_sim),np.imag(yval_sim))
         self.axes.set_xlabel("Re(S)")
         self.axes.set_ylabel("Im(S)")
     else:
         self.axes.set_xlabel(xlabel+' ['+xunit+']')
             
     if plottype=='amp':
         if ampformat=="log":
             self.axes.plot(xval,10*np.log(np.absolute(yval)),xval,10*np.log(np.absolute(yval_sim)))
             self.axes.set_ylabel("Amplitude [dB]")
         if ampformat=='lin':
             self.axes.plot(xval,np.absolute(yval),xval,np.absolute(yval_sim))
             self.axes.set_ylabel("Amplitude")      
         
     if plottype=='phase':
         if yval_sim==None: #this option is needed for lorentz function since we only fit the amplitude
             self.axes.plot(xval,np.angle(yval))
         else:                
             self.axes.plot(xval,np.angle(yval),xval,np.angle(yval_sim))
         self.axes.set_ylabel('Phase [rad]') 
         
     if save:
         print_fig = self.figure
         print_fig.savefig(Dir)
     else:
         self.draw()

def fidelity(statevec, dm, u_dm, ustatevec=np.array([0,0])):
    '''
    returns the fidelity (and its uncertainty) of the measured density
    matrix with a given state vector.
    '''
   
    f = error.Formula()
    beta,a,x,b,alpha = sympy.symbols('beta,a,x,b,alpha')
    v = Matrix([alpha, beta])
    rho = Matrix([[x,a+1j*b],[a-1j*b, 1-x]])    
    f.formula = (v.conjugate().transpose() * rho * v)[0]
    
    f.values[alpha] = statevec[0]
    f.values[beta] = statevec[1]
    f.values[a]=float(np.real(dm[0,1]))
    f.values[x]=float(dm[0,0])
    f.values[b]=float(np.imag(dm[0,1]))

    f.uncertainties[alpha]=ustatevec[0]
    f.uncertainties[beta]=ustatevec[1]
    f.uncertainties[x]=u_dm[0,0]
    f.uncertainties[a]=float(np.real(u_dm[0,1]))
    f.uncertainties[b]=float(np.imag(u_dm[0,1]))
    
    _fid,_ufid = f.num_eval()
    fid = float(_fid.as_real_imag()[0])
    ufid = float(_ufid.as_real_imag()[0])
   
    return (fid,ufid)

def rot_ell(m_rt_ps):
    '''Utility to compute rotation and ellipticity
       starting from reflection and transmission matrix

    Parameters
    ----------
    'm_RTsp' = sp reflection and transmission matrix

    Returns
    -------
    'a dictionary' = {'theta_p':theta_p,
                      'eps_p':eps_p,
                      'theta_s':theta_s,
                      'eps_s':eps_s}
    '''

    # extracting values from the matrix
    rt_pp = m_rt_ps[0,0]
    rt_ps = m_rt_ps[0,1]
    rt_sp = m_rt_ps[1,0]
    rt_ss = m_rt_ps[1,1]

    # calculating the values
    theta_p = np.real(rt_sp/rt_pp)
    eps_p = np.imag(rt_sp/rt_pp)
    theta_s = np.real(rt_ps/rt_ss)
    eps_s = np.imag(rt_ps/rt_ss)

    out_dict = {'theta_p':theta_p,
                'eps_p':eps_p,
                'theta_s':theta_s,
                'eps_s':eps_s}

    return out_dict

    def check_sum_rule(self, df1_w=None, df2_w=None):
        """Check f-sum rule."""

	if df1_w is None:
            df1_w = self.df1_w
            df2_w = self.df2_w

        N1 = N2 = 0
        for iw in range(self.Nw):
            w = iw * self.dw
            N1 += np.imag(df1_w[iw]) * w
            N2 += np.imag(df2_w[iw]) * w
        N1 *= self.dw * self.vol / (2 * pi**2)
        N2 *= self.dw * self.vol / (2 * pi**2)

        self.printtxt('')
        self.printtxt('Sum rule for ABS:')
        nv = self.nvalence
        self.printtxt('Without local field: N1 = %f, %f  %% error' %(N1, (N1 - nv) / nv * 100) )
        self.printtxt('Include local field: N2 = %f, %f  %% error' %(N2, (N2 - nv) / nv * 100) )

        N1 = N2 = 0
        for iw in range(self.Nw):
            w = iw * self.dw
            N1 -= np.imag(1/df1_w[iw]) * w
            N2 -= np.imag(1/df2_w[iw]) * w
        N1 *= self.dw * self.vol / (2 * pi**2)
        N2 *= self.dw * self.vol / (2 * pi**2)
                
        self.printtxt('')
        self.printtxt('Sum rule for EELS:')
        nv = self.nvalence
        self.printtxt('Without local field: N1 = %f, %f  %% error' %(N1, (N1 - nv) / nv * 100) )
        self.printtxt('Include local field: N2 = %f, %f  %% error' %(N2, (N2 - nv) / nv * 100) )

 def GaussLk( self, newLkFile, sigma ) :
   fout = open( newLkFile, "w" )
   fout.write("# input data : %s\n" % "+ Gaussian error"  )
   percentU = percentQ = 0.
   for lineStr,f1,f2 in zip (self.LeakList[0].lineStr, self.LeakList[0].f1, self.LeakList[0].f2) :
     print lineStr
     fout.write("#\n")
     for ant in range(1,16) :
       DRlist = []
       DLlist = []
       for Lk in self.LeakList :
         if Lk.ant == ant :
           for lineStr1,DR1,DL1 in zip( Lk.lineStr, Lk.DR, Lk.DL ) :
             if (lineStr1 == lineStr) and (abs(DR1) > 0.) and (abs(DL1) > 0.) :
               DRlist.append( DR1 )
               DLlist.append( DL1 )
               print "... ant %d - appending data from %s" % ( ant, Lk.legend )
       if len(DRlist) > 0 :
         DRmean = numpy.mean(DRlist)
         DLmean = numpy.mean(DLlist)
         DRnew = random.gauss(numpy.real(DRmean), sigma) + random.gauss(numpy.imag(DRmean), sigma) * 1j
         DLnew = random.gauss(numpy.real(DLmean), sigma) + random.gauss(numpy.imag(DLmean), sigma) * 1j
       else :
         DRnew = 0. + 0j
         DLnew = 0. + 0j
       print ant, DRnew, DLnew
       fout.write("C%02d %8.3f %8.3f %8.3f %6.3f %8.3f %6.3f %8.3f %6.3f    %s\n" % \
           ( ant, f1, f2, DRnew.real, DRnew.imag, DLnew.real, \
           DLnew.imag, percentQ, percentU, lineStr) )
   fout.close()

def normalmode_frequencies(hessian, metric=None, eps=1e-4):
    """calculate (squared) normal mode frequencies

    Parameters
    ----------
    hessian: 2d array
        hessian matrix
    metric: 2d array
        mass weighted metric tensor

    Returns
    -------
    sorted array of normal mode frequencies

    """
    A = hessian
    if metric is not None:
        A = np.dot(np.linalg.pinv(metric), hessian)

    frq = np.linalg.eigvals(A)

    if np.max(np.abs(np.imag(frq))) > eps:
        print(frq)
        raise ValueError("imaginary eigenvalue in frequency calculation"
                         ", check hessian + metric tensor\n"
                         "the largest imaginary part is %g" %
                         np.max(np.abs(np.imag(frq))))

    return np.sort(np.real(frq))

def test_mexh():
    LB = -5
    UB = 5
    N = 1000

    [psi, x] = ref_mexh(LB, UB, N)
    w = pywt.ContinuousWavelet("mexh")
    w.upper_bound = UB
    w.lower_bound = LB
    PSI, X = w.wavefun(length=N)

    assert_allclose(np.real(PSI), np.real(psi))
    assert_allclose(np.imag(PSI), np.imag(psi))
    assert_allclose(X, x)

    LB = -5
    UB = 5
    N = 1001

    [psi, x] = ref_mexh(LB, UB, N)
    w = pywt.ContinuousWavelet("mexh")
    w.upper_bound = UB
    w.lower_bound = LB
    PSI, X = w.wavefun(length=N)

    assert_allclose(np.real(PSI), np.real(psi))
    assert_allclose(np.imag(PSI), np.imag(psi))
    assert_allclose(X, x)

def ARLineSpectra(ar):
    """
    Convert AR coeffs to LSPs

    From wikipedia:
    A palindromic polynomial (i.e., P) of odd degree has -1 as a root.
    An antipalindromic polynomial (i.e., Q) has 1 as a root.
    An antipalindromic polynomial of even degree has -1 and 1 as roots
    """
    order = ar.shape[-1]
    ret = np.zeros(ar.shape)
    for a, o in core.refiter([ar, ret], core.newshape(ar.shape)):
        p = np.ones((order+2))
        q = np.ones((order+2))
        q[-1] = -1.0
        for i in range(order):
            p[i+1] = -a[i] - a[order-i-1]
            q[i+1] = -a[i] + a[order-i-1]
        pr = np.roots(p)
        qr = np.roots(q)

        j = 0
        an = np.ndarray((order+2))
        for i in range(len(pr)):
            if np.imag(pr[i]) >= 0.0:
                an[j] = np.angle(pr[i])
                j += 1
            if np.imag(qr[i]) >= 0.0:
                an[j] = np.angle(qr[i])
                j += 1
        # The angle list (an) will always contain both 0 and pi; they
        # will move to the ends after the sort
        o[...] = np.sort(an)[1:-1]
    return ret;

def getArgumentVariable(_ComplexVariable):

	#Debug
	'''
	print('l 31 Numscipier')
	print('_ComplexVariable is ')
	print(_ComplexVariable)
	print('')
	'''

	#import
	import numpy as np

	#return
	return 2.*np.arctan(
	np.imag(_ComplexVariable)/(
	        np.sqrt(
	            np.imag(
	                _ComplexVariable
	            )**2+np.real(
	                _ComplexVariable
	            )**2)+np.real(
	                _ComplexVariable
	            )
	    )
	);

    def mix_parameters(self, Pibra, Piket):
        r"""Mix the two parameter sets :math:`\Pi_i` and :math:`\Pi_j`
        from the 'bra' and the 'ket' wavepackets :math:`\Phi\left[\Pi_i\right]`
        and :math:`\Phi^\prime\left[\Pi_j\right]`.

        :param Pibra: The parameter set :math:`\Pi_i` from the bra part wavepacket.
        :param Piket: The parameter set :math:`\Pi_j` from the ket part wavepacket.
        :return: The mixed parameters :math:`q_0` and :math:`Q_S`. (See the theory for details.)
        """
        # <Pibra | ... | Piket>
        qr, pr, Qr, Pr = Pibra
        qc, pc, Qc, Pc = Piket

        # Mix the parameters
        Gr = dot(Pr, inv(Qr))
        Gc = dot(Pc, inv(Qc))

        r = imag(Gc - conjugate(Gr.T))
        s = imag(dot(Gc, qc) - dot(conjugate(Gr.T), qr))

        q0 = dot(inv(r), s)
        Q0 = 0.5 * r

        # Here we can not avoid the matrix root by using svd
        Qs = inv(sqrtm(Q0))

        return (q0, Qs)