def SusceptibilityHF(U,GF_A,X_A):
	''' susceptibility calculated from the full spectral self-energy derivative '''
	Int1_A = FD_A*sp.imag(GF_A**2*(1.0-U*X_A))
	Int2_A = FD_A*sp.imag(GF_A**2*X_A)
	I1 = simps(Int1_A,En_A)/sp.pi
	I2 = simps(Int2_A,En_A)/sp.pi
	return 2.0*I1/(1.0+U**2*I2)

def XIntegralsFFT(GF_A,Bubble_A,Lambda,BubZero):
	''' calculate X integral to susceptibilities using FFT '''
	N = int((len(En_A)-1)/2)
	Kappa_A  = TwoParticleBubble(GF_A,GF_A**2,'eh')
	Bubble_A = TwoParticleBubble(GF_A,GF_A,'eh')
	#print(Kappa_A[N],Bubble_A[N])
	V_A   = 1.0/(1.0+Lambda*Bubble_A)
	KV_A  = Lambda*Kappa_A*V_A**2
	KmV_A = Lambda*sp.flipud(sp.conj(Kappa_A))*V_A**2
	## zero-padding the arrays
	exFD_A  = sp.concatenate([FD_A[N:],sp.zeros(2*N+2),FD_A[:N+1]])
	ImGF_A  = sp.concatenate([sp.imag(GF_A[N:]),sp.zeros(2*N+2),sp.imag(GF_A[:N+1])])
	ImGF2_A = sp.concatenate([sp.imag(GF_A[N:]**2),sp.zeros(2*N+2),sp.imag(GF_A[:N+1]**2)])
	ImV_A   = sp.concatenate([sp.imag(V_A[N:]),sp.zeros(2*N+2),sp.imag(V_A[:N+1])])
	ImKV_A  = sp.concatenate([sp.imag(KV_A[N:]),sp.zeros(2*N+2),sp.imag(KV_A[:N+1])])
	ImKmV_A = sp.concatenate([sp.imag(KmV_A[N:]),sp.zeros(2*N+2),sp.imag(KmV_A[:N+1])])
	## performing the convolution
	ftImX11_A = -sp.conj(fft(exFD_A*ImV_A))*fft(ImGF2_A)*dE
	ftImX12_A =  fft(exFD_A*ImGF2_A)*sp.conj(fft(ImV_A))*dE
	ftImX21_A = -sp.conj(fft(exFD_A*ImKV_A))*fft(ImGF_A)*dE
	ftImX22_A =  fft(exFD_A*ImGF_A)*sp.conj(fft(ImKV_A))*dE
	ftImX31_A = -sp.conj(fft(exFD_A*ImKmV_A))*fft(ImGF_A)*dE
	ftImX32_A =  fft(exFD_A*ImGF_A)*sp.conj(fft(ImKmV_A))*dE
	## inverse transform
	ImX1_A =  sp.real(ifft(ftImX11_A+ftImX12_A))/sp.pi
	ImX2_A =  sp.real(ifft(ftImX21_A+ftImX22_A))/sp.pi
	ImX3_A = -sp.real(ifft(ftImX31_A+ftImX32_A))/sp.pi
	ImX1_A =  sp.concatenate([ImX1_A[3*N+4:],ImX1_A[:N+1]])
	ImX2_A =  sp.concatenate([ImX2_A[3*N+4:],ImX2_A[:N+1]])
	ImX3_A =  sp.concatenate([ImX3_A[3*N+4:],ImX3_A[:N+1]])
	## getting real part from imaginary
	X1_A = KramersKronigFFT(ImX1_A) + 1.0j*ImX1_A + BubZero # constant part !!!
	X2_A = KramersKronigFFT(ImX2_A) + 1.0j*ImX2_A
	X3_A = KramersKronigFFT(ImX3_A) + 1.0j*ImX3_A
	return [X1_A,X2_A,X3_A]

def root_locus(sys, kvect, xlim=None, ylim=None, plotstr='-', Plot=True, 
               PrintGain=True):
    """Calculate the root locus by finding the roots of 1+k*TF(s)
    where TF is self.num(s)/self.den(s) and each k is an element
    of kvect.

    Parameters
    ----------
    sys : linsys
        Linear input/output systems (SISO only, for now)
    kvect : gain_range (default = None)
        List of gains to use in computing diagram
    Plot : boolean (default = True)
        If True, plot magnitude and phase
    PrintGain: boolean (default = True)
        If True, report mouse clicks when close to the root-locus branches,
        calculate gain, damping and print
    Return values
    -------------
    rlist : list of computed root locations
    """

    # Convert numerator and denominator to polynomials if they aren't
    (nump, denp) = _systopoly1d(sys);

    # Compute out the loci
    mymat = _RLFindRoots(sys, kvect)
    mymat = _RLSortRoots(sys, mymat)

    # Create the plot
    if (Plot):
        f = pylab.figure()
        if PrintGain:
            cid = f.canvas.mpl_connect(
                'button_release_event', partial(_RLFeedbackClicks, sys=sys))
        ax = pylab.axes();

        # plot open loop poles
        poles = array(denp.r)
        ax.plot(real(poles), imag(poles), 'x')

        # plot open loop zeros
        zeros = array(nump.r)
        if zeros.any():
            ax.plot(real(zeros), imag(zeros), 'o')

        # Now plot the loci
        for col in mymat.T:
            ax.plot(real(col), imag(col), plotstr)

        # Set up plot axes and labels
        if xlim:
            ax.set_xlim(xlim)
        if ylim:
            ax.set_ylim(ylim)
        ax.set_xlabel('Real')
        ax.set_ylabel('Imaginary')

    return mymat

def confMap(shape,mapfunc):
	shapemapped = [None]*len(shape)
	for i in range(0,len(shape)):
		shapemapped[i] = mapfunc(shape[i])
	
	plt.scatter(sp.real(shape),sp.imag(shape),color='r')
	plt.scatter(sp.real(shapemapped),sp.imag(shapemapped),color='b')
	plt . show ()

def solve_P(F,G,H):
    """This function takes arguments for F,G,H and solves the matrix quadratic given by
    F*P^2+G*P+H=0.  Note F, G, and H must be square.
    The function returns the matrix P and the resulting matrix, given by F*P^2+G*P+H 
    which should be close to zero.
    The algorithm used to solve for P is outlined in 'A Toolkit for Analyzing Nonlinear 
    Dynamic Stochastic Models Easily' by Harald Uhlig.
    """
    m=sp.shape(F)[0]
    
    Xi=sp.concatenate((-G,-H), axis=1)
    second=sp.concatenate((sp.eye(m,m),sp.zeros((m,m))),axis=1)
    Xi=sp.concatenate((Xi,second))
    
    Delta=sp.concatenate((F,sp.zeros((m,m))),axis=1)
    second=sp.concatenate((sp.zeros((m,m)),sp.eye(m,m)),axis=1)
    Delta=sp.concatenate((Delta,second))
    
    (L,V) = la.eig(Xi,Delta)
    
    boolean = sp.zeros(len(L))
    trueCount =0
    
    for i in range(len(L)):
        if L[i]<1 and L[i]>-1 and sp.imag(L[i])==0 and trueCount<m:
            boolean[i] = True
            trueCount+=1
    #display(L, boolean)
    if trueCount<m:
        print "Imaginary eigenvalues being used"
        for i in range(len(L)):
            if math.sqrt(real(L[i])**2+imag(L[i])**2)<1 and trueCount<m:
                boolean[i]=True
                trueCount+=1
    #display(boolean)
    
    if trueCount==m:
        print "true count is m"
        Omega=sp.zeros((m,m))
        diagonal=[]
        count =0
        for i in range(len(L)):
            if boolean[i]==1:
                Omega[:,count]=sp.real(V[m:2*m,i])+sp.imag(V[m:2*m,i])
                diagonal.append(L[i])
                count+=1
        Lambda=sp.diag(diagonal)
        try:
            P=sp.dot(sp.dot(Omega,Lambda),la.inv(Omega))
        except:
            print 'Omega not invertable'
            P=sp.zeros((m,m))
        diff=sp.dot(F,sp.dot(P,P))+sp.dot(G,P)+H
        return P,diff
    else:
        print "Problem with input, not enough 'good' eigenvalues"
        return sp.zeros((m,m)),sp.ones((m,m))*100

def pltFunction (cavity1,cavity2,cavity3,plotType):
  if (plotType==0):
    return sp.absolute(cavity1[:])**2,sp.absolute(cavity2[:])**2,sp.absolute(cavity3[:])**2
  elif (plotType==1):
    return sp.real(cavity1[:]),sp.real(cavity2[:]),sp.real(cavity3[:])
  elif (plotType==2):
    return sp.imag(cavity1[:]),sp.imag(cavity2[:]),sp.imag(cavity3[:])
  else: 
    return cavity1, cavity2, cavity3

def SelfEnergy(GF_A,ChiGamma_A):
	''' calculating the dynamical self-energy from the Schwinger-Dyson equation '''
	N = int((len(En_A)-1)/2)
	## zero-padding the arrays
	exFD_A = sp.concatenate([FD_A[N:],sp.zeros(2*N+3),FD_A[:N]])
	exBE_A = sp.concatenate([BE_A[N:],sp.zeros(2*N+3),BE_A[:N]])
	ImGF_A = sp.concatenate([sp.imag(GF_A[N:]),sp.zeros(2*N+3),sp.imag(GF_A[:N])])
	ImCG_A = sp.concatenate([sp.imag(ChiGamma_A[N:]),sp.zeros(2*N+3),sp.imag(ChiGamma_A[:N])])
	## performing the convolution
	ftImSE1_A = -sp.conj(fft(exBE_A*ImCG_A))*fft(ImGF_A)*dE
	ftImSE2_A = -fft(exFD_A*ImGF_A)*sp.conj(fft(ImCG_A))*dE
	ImSE_A = sp.real(ifft(ftImSE1_A+ftImSE2_A))/sp.pi
	ImSE_A = sp.concatenate([ImSE_A[3*N+4:],ImSE_A[:N+1]])
	Sigma_A = KramersKronigFFT(ImSE_A) + 1.0j*ImSE_A
	return Sigma_A

    def process(self, X, V, C):
        """Tolerance for eigenvalues a possible problem when checking for neutral saddles."""
        BifPoint.process(self, X, V, C)
                
        J_coords = C.CorrFunc.jac(X, C.coords)
        eigs, LV, RV = linalg.eig(J_coords,left=1,right=1)
        
        # Check for neutral saddles
        found = False
        for i in range(len(eigs)):
            if abs(imag(eigs[i])) < 1e-5:
                for j in range(i+1,len(eigs)):
                    if C.verbosity >= 2:
                        if abs(eigs[i]) < 1e-5 and abs(eigs[j]) < 1e-5:
                            print 'Fold-Fold point found in Hopf!\n'
                        elif abs(imag(eigs[j])) < 1e-5 and abs(real(eigs[i]) + real(eigs[j])) < 1e-5:
                            print 'Neutral saddle found!\n'
            elif abs(real(eigs[i])) < 1e-5:
                for j in range(i+1, len(eigs)):
                    if abs(real(eigs[j])) < 1e-5 and abs(real(eigs[i]) - real(eigs[j])) < 1e-5:
                        found = True
                        w = abs(imag(eigs[i]))
                        if imag(eigs[i]) > 0:
                            p = conjugate(LV[:,j])/linalg.norm(LV[:,j])
                            q = RV[:,i]/linalg.norm(RV[:,i])
                        else:
                            p = conjugate(LV[:,i])/linalg.norm(LV[:,i])
                            q = RV[:,j]/linalg.norm(RV[:,j])

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

        direc = conjugate(1/matrixmultiply(conjugate(p),q))
        p = direc*p

        # Alternate way to compute 1st lyapunov coefficient (from Kuznetsov [4])
        
        #print (1./(w*w))*real(1j*matrixmultiply(conjugate(p),b1)*matrixmultiply(conjugate(p),b3) + \
        #   w*matrixmultiply(conjugate(p),trilinearform(D,q,q,conjugate(q))))

        self.found[-1].w = w
        self.found[-1].l1 = firstlyapunov(X, C.CorrFunc, w, J_coords=J_coords, p=p, q=q, check=(C.verbosity==2)) 
        self.found[-1].eigs = eigs
        
        self.info(C, -1)
        
        return True

 def write(self, vals, name, label='Mie'):
     with open(name, 'w') as fp:
         N = len(self.elems)
         fp.write("View \"{0}\" {{\n".format(label))
         for ii, elem in enumerate(self.elems):
             fp.write('ST( ')
             for jj, n in enumerate(elem):
                 x, y, z = self.points[n-1]
                 fp.write('{0}, {1}, {2}'.format(x, y, z))
                 if jj == len(elem)-1:
                     fp.write(')')
                 else:
                     fp.write(', ')
             fp.write('{')
             for jj, n in enumerate(elem):
                 rv = sp.real(vals[n-1])
                 fp.write('{0}, '.format(rv))
             for jj, n in enumerate(elem):
                 iv = sp.imag(vals[n-1])
                 fp.write('{0}'.format(iv))
                 if jj == len(elem)-1:
                     fp.write('};\n')
                 else:
                     fp.write(', ')
             self._print_progress(ii, N, 'written')
         fp.write("TIME{0, 1};\n};\n\n")
         print(' --> {0} written.'.format(name))

def construct(phi1, phi2, nomod = 0, amp1 =[], amp2=[], eta = 0, ampout= 0): #does ampout need to be there?
	if len(amp1) > 0 or len(amp2) > 0:
		tempshape = phi1.shape
		w = tempshape[1]
		h = tempshape[0]
		if len(amp1) == 0:
			temp1 = np.ones(w)
			temp2 = np.ones(h)
			for r in temp2:
				amp1 += [temp1]
		if len(amp2) == 0:
			temp1 = np.ones(w)
			temp2 = np.ones(h)
			for r in temp2:
				amp2 += [temp1]
		psi1 = amp1 * np.exp(1j*phi1)
		psi2 = amp2 * np.exp(1j*phi2)
		psi = psi1 * psi2
		psi = np.array(psi)
		apsi = abs(psi)
		psi = psi/(np.amax(abs(psi)))
		phi = np.arctan2(sp.real(psi),sp.imag(psi))
		phi -= np.amin(phi)
		phi = phi % (2.*np.pi)
		eta = 2*np.median(abs(psi))
		randarray = np.array([[random.random() for i in range(w)] for j in range(h)])
		shape = (abs(psi) >= (eta*randarray)) 
		index = np.where(shape == False)
		phi[index] = 0
		ampout = abs(psi)
	else:
		phi = phi1 + phi2
		phi = phi - np.amin(phi)
		phi = phi % (2.*np.pi)
	return phi

 def __init__(self, output='out', input='in', \
              mag=None, phase=None, coh=None, \
              freqlim=[], maglim=[], phaselim=[], \
              averaged='not specified', \
              seedfreq=-1, seedphase=0,
              labels=[], legloc=-1, compin=[]):
     self.output = output
     self.input = input
     if len(compin) > 0:
        if mag is None:
           self.mag = squeeze(colwise(abs(compin)))
        if phase is None:
           self.phase = squeeze(colwise(arctan2(imag(compin),real(compin))*180.0/pi))
     else:
         self.mag = squeeze(mag)
         self.phase = squeeze(phase)
     self.coh = coh
     self.averaged = averaged
     self.seedfreq = seedfreq
     self.seedphase = seedphase
     self.freqlim = freqlim
     self.maglim = maglim
     self.phaselim = phaselim
     self.labels = labels
     self.legloc = legloc

def dst(x,axis=-1):
    """Discrete Sine Transform (DST-I)

    Implemented using 2(N+1)-point FFT
    xsym = r_[0,x,0,-x[::-1]]
    DST = (-imag(fft(xsym))/2)[1:(N+1)]

    adjusted to work over an arbitrary axis for entire n-dim array
    """
    n = len(x.shape)
    N = x.shape[axis]
    slices = [None]*3
    for k in range(3):
        slices[k] = []
        for j in range(n):
            slices[k].append(slice(None))    
    newshape = list(x.shape)
    newshape[axis] = 2*(N+1)
    xtilde = scipy.zeros(newshape, dtype=float)
    slices[0][axis] = slice(1,N+1)
    slices[1][axis] = slice(N+2,None)
    slices[2][axis] = slice(None,None,-1)
    for k in range(3):
        slices[k] = tuple(slices[k])
    xtilde[slices[0]] = x
    xtilde[slices[1]] = -x[slices[2]]
    Xt = scipy.fft(xtilde,axis=axis)
    return (-scipy.imag(Xt)/2)[slices[0]]

def plot_pole(state):
    """
    Plots the momentum pole of a state in the complex plane.
    """
    mass = state.problem.mass
    energy = state.energy
    k = sp.sqrt(2 * mass * energy)
    return plt.plot(sp.real(k), sp.imag(k), "o", color="red")

def CalculateHWHM(GF_A):
	''' calculates the half-width at half-maximum of the Kondo resonance 
	and the maximum of the spectral function '''
	N = len(En_A)
	IntMin = int((N+1)/2-int(0.5/dE))
	IntMax = int((N+1)/2+int(0.5/dE))
	DOSmaxPos = sp.argmax(-sp.imag(GF_A[IntMin:IntMax])/sp.pi)
	DOSmax    = -sp.imag(GF_A[IntMin+DOSmaxPos])/sp.pi # maximum of DoS
	wmax      = En_A[IntMin+DOSmaxPos]                 # position of the maximum at energy axis
	DOS = InterpolatedUnivariateSpline(En_A-1e-12,-sp.imag(GF_A)/sp.pi-DOSmax/2.0) 
	## 1e-12 breaks symmetry for half-filling, otherway DOS.roots() loses one solution.
	DOSroots_A = sp.sort(sp.fabs(DOS.roots()))
	try:
		HWHM = (DOSroots_A[0] + DOSroots_A[1])/2.0
	except IndexError:
		HWHM = 0.0
	return [HWHM,DOSmax,wmax]

def upconvert(x, fc):
    """
    Upconverts in-phase and quadrature components of a time-domain signal
    appropriately. Returns the in-phase and quadrature signals
    """
    t = sp.arange(x.size)
    x_i = sp.real(x) * 2 * sp.cos(2*pi*fc*t / sampling_rate)
    x_q = sp.imag(x) * 2 * sp.sin(2*pi*fc*t / sampling_rate)
    return x_i, x_q

def TwoParticleBubble(F1_A,F2_A,channel):
	''' calculates the two-particle bubble, channel = 'eh', 'ee' '''
	N = int((len(En_A)-1)/2)
	## zero-padding the arrays
	exFD_A = sp.concatenate([FD_A[N:],sp.zeros(2*N+3),FD_A[:N]])
	ImF1_A = sp.concatenate([sp.imag(F1_A[N:]),sp.zeros(2*N+3),sp.imag(F1_A[:N])])
	ImF2_A = sp.concatenate([sp.imag(F2_A[N:]),sp.zeros(2*N+3),sp.imag(F2_A[:N])])
	## performing the convolution
	if channel == 'eh':
		ftImChi1_A =  sp.conj(fft(exFD_A*ImF2_A))*fft(ImF1_A)*dE	# f(x)F2(x)F1(w+x)
		ftImChi2_A = -fft(exFD_A*ImF1_A)*sp.conj(fft(ImF2_A))*dE	# f(x)F1(x)F2(x-w)
	elif channel == 'ee':
		ftImChi1_A =  fft(exFD_A*ImF2_A)*fft(ImF1_A)*dE					# f(x)F2(x)F1(w-x)
		ftImChi2_A = -sp.conj(fft(exFD_A*sp.flipud(ImF1_A)))*fft(ImF2_A)*dE	# f(x)F1(-x)F2(w+x)
	ImChi_A = -sp.real(ifft(ftImChi1_A+ftImChi2_A))/sp.pi
	ImChi_A = sp.concatenate([ImChi_A[3*N+4:],ImChi_A[:N+1]])
	Chi_A = KramersKronigFFT(ImChi_A) + 1.0j*ImChi_A
	return Chi_A

def RiemannSurface2():
    """riemann surface for imaginary part of sqrt(z)"""

    fig = plt.figure()
    ax = Axes3D(fig)
    X = sp.arange(-5, 5, 0.25)
    Y = sp.arange(-5, 0, 0.25)
    X, Y = sp.meshgrid(X, Y)
    Z = sp.imag(sp.sqrt(X+1j*Y))
    ax.plot_surface(X, Y, Z, rstride=1, cstride=1, linewidth=0, cmap=cmap_r)
    ax.plot_surface(X, Y, -Z, rstride=1, cstride=1, linewidth=0, cmap=cmap)
    X = sp.arange(-5, 5, 0.25)
    Y = sp.arange(0,5,.25)
    X, Y = sp.meshgrid(X, Y)
    Z = sp.imag(sp.sqrt(X+1j*Y))
    ax.plot_surface(X, Y, -Z, rstride=1, cstride=1, linewidth=0, cmap=cmap_r)
    ax.plot_surface(X, Y, Z, rstride=1, cstride=1, linewidth=0, cmap=cmap)
    plt.savefig('RiemannSurface2.pdf', bbox_inches='tight', pad_inches=0)

def fourier_series(dat, t, n):
    """
    Fourier series approximation to a function.
    returns Fourier coefficients of a function.
    The coefficients are numerical approximations of the true
    coefficients.

    Parameters
    ----------
    dat: array
        Array of data representing the function.
    t: array
        Corresponding time array.
    n: int
        The desired number of terms to use in the
        Fourier series.

    Returns
    ----------
    a, b: tuple
        Tuple containing arrays with the Fourier coefficients.
        The function also produces a plot of the approximation.

    Examples:
    >>> f = sp.hstack((sp.arange(-1, 1, .04), sp.arange(1, -1, -.04)))
    >>> f += 1
    >>> t = sp.arange(0, len(f))/len(f)
    >>> a, b = fourier_series(f, t, 5)
    >>> a[0]
    2.0
    """

    len_ = len(dat)/2
    fs = (sp.fft(dat))/len_
    a0 = fs[0]
    a = sp.real(sp.hstack((a0, fs[1:len(fs/2)])))
    b = -sp.imag(fs[1:len(fs/2)])
    len_ *= 2
    dt = 2*sp.pi/len_
    tp = sp.arange(0, 2*sp.pi, dt)
    dataapprox = a[0]/2 + sp.zeros_like(dat)
    fig = plt.figure()
    ax1 = fig.add_subplot(111)
    ax1.plot(t, dat)

    for i in range(1, n):
        newdat = a[i]*sp.cos(tp*i) + b[i]*sp.sin(tp*i)
        dataapprox += newdat
        if i == n-1:
            ax1.plot(t, newdat)

    ax1.plot(t, dataapprox)

    _ = plt.show()

    return a, b

def sliceEpsArr(EpsArr,NX,NY,NZ,figNum):
    """ Visualizing object with slices 
    """
    vxyR=sci.real(EpsArr[:,:,round(NZ/2.0)])
    vxzR=sci.real(EpsArr[:,round(NY/2.0),:])
    vyzR=sci.real(EpsArr[round(NX/2.0),:,:])
    vxyI=sci.imag(EpsArr[:,:,round(NZ/2.0)])
    vxzI=sci.imag(EpsArr[:,round(NY/2.0),:])
    vyzI=sci.imag(EpsArr[round(NX/2.0),:,:])
    
    plt.ion()
    fig=plt.figure(figNum)
    fig.clear()
    
    plt.subplot(231)
    plt.imshow(vxyR)
    plt.title('real($\\varepsilon$), xy-plane')
    plt.colorbar()
    plt.subplot(232)
    plt.imshow(sci.swapaxes(vxzR,0,1))
    plt.title('real($\\varepsilon$), xz-plane')
    plt.colorbar()
    plt.subplot(233)
    plt.imshow(sci.swapaxes(vyzR,0,1))
    plt.title('real($\\varepsilon$), yz-plane')
    plt.colorbar()
    plt.subplot(234)
    plt.imshow(vxyI)
    plt.title('imag($\\varepsilon$), xy-plane')
    plt.colorbar()
    plt.subplot(235)
    plt.imshow(sci.swapaxes(vxzI,0,1))
    plt.title('imag($\\varepsilon$), xz-plane')
    plt.colorbar()
    plt.subplot(236)
    plt.imshow(sci.swapaxes(vyzI,0,1))
    plt.title('imag($\\varepsilon$), yz-plane')
    plt.colorbar()
    
    fig.canvas.draw()
    plt.ioff()
    
    return

def plotEpsArr(EpsArr,NX,NY,NZ,figNum,hold,label,clear):
    """ Visualizing object with cuts 
    """
    cutxR=sci.real(EpsArr[:,round(NY/2.0),round(NZ/2.0)])
    cutxI=sci.imag(EpsArr[:,round(NY/2.0),round(NZ/2.0)])
    cutyR=sci.real(EpsArr[round(NX/2.0),:,round(NZ/2.0)])
    cutyI=sci.imag(EpsArr[round(NX/2.0),:,round(NZ/2.0)])
    cutzR=sci.real(EpsArr[round(NX/2.0),round(NY/2.0),:])
    cutzI=sci.imag(EpsArr[round(NX/2.0),round(NY/2.0),:])
    
    plt.ion()
    fig=plt.figure(figNum)
    fig.clear()
    if clear=='yes':
        fig.clear()

    plt.hold(hold)
    
    plt.subplot(231)
    plt.plot(cutxR,label)
    plt.title('real($\\varepsilon$), x-axis cut')
    plt.subplot(232)
    plt.plot(cutyR,label)
    plt.title('real($\\varepsilon$), y-axis cut')
    plt.subplot(233)
    plt.plot(cutzR,label)
    plt.title('real($\\varepsilon$), z-axis cut')
    plt.subplot(234)
    plt.plot(cutxI,label)
    plt.title('imag($\\varepsilon$), x-axis cut')
    plt.subplot(235)
    plt.plot(cutyI,label)
    plt.title('imag($\\varepsilon$), y-axis cut')
    plt.subplot(236)
    plt.plot(cutzI,label)
    plt.title('imag($\\varepsilon$), z-axis cut')
    
    plt.hold('False')
    fig.canvas.draw()
    plt.ioff()
    
    return