def crystal_embed(e,eta):
    """
    Evaluates crystal embedding potential by solving the one-dimensional
    Schroedinger equation through one unit cell in both directions.
    """
    y0=[1.0,0.0,0.0,0.0]
    z1=np.linspace(-z_left,-z_left+alat,101)
    sol1=odeint(schroed,y0,z1,args=(e,eta))
    z2=np.linspace(-z_left+alat,-z_left,101)
    sol2=odeint(schroed,y0,z2,args=(e,eta))
    psi1=complex(sol1[50,0],sol1[50,1])
    psi1_prime=complex(sol1[50,2],sol1[50,3])
    psi2=complex(sol2[50,0],sol2[50,1])
    psi2_prime=complex(sol2[50,2],sol2[50,3])
    wronskian=psi1*psi2_prime-psi2*psi1_prime
    psi1=complex(sol1[100,0],sol1[100,1])
    psi2=complex(sol2[100,0],sol2[100,1])
    cos_ka=0.5*(psi1+psi2)
    ka=arccos(cos_ka)
    if ka.imag<0: ka=np.conj(ka)
    exp_ka=cos_ka+1.0j*sqrt(1.0-cos_ka**2)
    emb_cryst=0.5*wronskian/(exp_ka-psi2)
    if emb_cryst.imag>0.0:
        exp_ka=cos_ka-1.0j*sqrt(1.0-cos_ka**2)
        emb_cryst=0.5*wronskian/(exp_ka-psi2)
    return emb_cryst

    def roots(self):
        """
        Return the two roots of the quadratic polynomial.
        The roots are real or complex, depending on the
        coefficients in the polynomial.

        Let us define a quadratic polynomial:

        >>> q = Quadratic(a=2, b=4, c=-16)
        >>> print q
        2*x**2 + 4*x - 16

        The roots are then found by

        >>> r1, r2 = q.roots()
        >>> r1
        2.0
        >>> r2
        -4.0
        """
        a, b, c = self.a, self.b, self.c  # short names
        q = b**2 - 4*a*c
        root1 = (-b + sqrt(q))/float(2*a)
        root2 = (-b - sqrt(q))/float(2*a)
        return root1, root2

    def hybridWith(self, v2, w2):
        """Create a weighted average of two voters.

        The weight of v1 is always 1; w2 is the weight of v2 relative to that.

        If both are
        standard normal to start with, the result will be standard normal too.

        Length must be the same
            >>> Voter([1,2]).hybridWith(Voter([1,2,3]),1)
            Traceback (most recent call last):
              ...
            AssertionError

        A couple of basic sanity checks:
            >>> v2 = Voter([1,2]).hybridWith(Voter([3,2]),1)
            >>> [round(u,5) for u in v2.hybridWith(v2,1)]
            [4.0, 4.0]
            >>> Voter([1,2,5]).hybridWith(Voter([-0.5,-1,0]),0.75)
            (0.5, 1.0, 4.0)
        """
        assert len(self) == len(v2)
        return self.copyWithUtils(  ((self[i] / sqrt(1 + w2 ** 2)) +
                                    (w2 * v2[i] / sqrt(1 + w2 ** 2)))
                                 for i in range(len(self)))

def _hexagonalPositions(nrow, ncol, radius):
    X = np.tile(np.arange(ncol) * 1.5 * radius, (nrow, 1)) + radius
    Y = np.tile(np.arange(nrow) * sqrt(3) * radius, (ncol, 1)).T + radius
    Y = Y + np.tile([0, sqrt(3) / 2 * radius], (nrow, np.floor(ncol / 2)))
    X = np.reshape(X, (X.size), order = 'F')
    Y = np.reshape(Y, (Y.size), order = 'F')
    return np.array([X, Y]).T

def QR_algorithm(A,niter,tol):
    """Return the eigenvalues of A using the QR algorithm."""
    H = la.hessenberg(A)
    for i in xrange(niter):
        Q,R = la.qr(H)
        H = np.dot(R,Q)
    S = H
    print S
    eigenvalues = []
    i = 0
    while i < S.shape[0]:
        if i == S.shape[0]-1:
            eigenvalues.append(S[i,i])
        elif abs(S[i+1,i]) < tol:
            eigenvalues.append(S[i,i])
        else:
            a = S[i,i]
            b = S[i,i+1]
            c = S[i+1,i]
            d = S[i+1,i+1]
            B = -1*(a+d)
            C = a*d-b*c
            eigen_plus = (-B + sm.sqrt(B**2 - 4*C))/2.
            eigen_minus = (-B - sm.sqrt(B**2 - 4*C))/2.
            
            eigenvalues.append(eigen_plus)
            eigenvalues.append(eigen_minus)
            i+=1
        i+=1
    return eigenvalues

def plotThroughput(logFile, windowSize = 10000, marker = 's', iteration = 50):
    fileHandle = open(logFile, 'r')
    

    wholeFile = fileHandle.read();
    wholeFile = wholeFile.split("\n")[:-1]
    wholeFileSorted = sorted(wholeFile, key = lambda x : float(x.split()[2]));
    
    globalMin = float(wholeFileSorted[0].split()[2])
    currentStart = globalMin
    currentEnd = currentStart + windowSize
    data = {(currentStart, currentEnd) : 0}
    for line in wholeFileSorted:
        min_val = float(line.split()[2])
        if min_val > currentEnd:
            currentStart = currentEnd
            currentEnd = currentStart + windowSize
            data[(currentStart, currentEnd)] = 0
        data[(currentStart, currentEnd)] += 1

    x = []
    y = []
    for key in sorted(data.keys()):
        x.append(((key[1] + key[0])/2.0 - globalMin)/1000.0)
        y.append(data[key]*1000/(key[1] - key[0]))
    p1, = plt.plot(x,y, marker = marker)
    ax = plt.gca()
    #ax.set_ylim([-10, 300])
    ax.set_xlabel('Experiment time (s)')
    ax.set_ylabel('Messages read per second')
    ax.set_title('Throughput of read operations')

    total = []
    for key in data:
        if((key[0] - globalMin)/1000. > 1500):
            continue
        if((key[0] - globalMin)/1000. < 300):
            continue
        total.append(data[key]*1000/(key[1] - key[0]))

    n, min_max, mean, var, _, _ = stats.describe(total)
    t_test = stats.t.interval(0.95, n - 1, loc=mean, scale=sqrt(var) / sqrt(n))
    chisquare = stats.chi2.interval(0.95, n - 1)
    #print "Data samples: %d" % n
    #print "Average response time: %s ms" % mean
    #print "Sample standard deviation: %s ms" % sqrt(var)
   # print "95%% confidence interval for the mean (t-student test): (%s, %s) ms" % stats.t.interval(0.95, n - 1,loc = mean, scale = sqrt(var)/sqrt(n))
    #chisquare = stats.chi2.interval(0.95, n - 1)
    #print "95%% confidence interval for the stdev (chi-square test): (%s, %s) ms" % (sqrt((n - 1)*var/chisquare[1]), sqrt((n - 1)*var/chisquare[0]))
    #m,b =  polyfit(x, y, 1)
    #print "slope: %s, intercept: %s" % (m, b)
    #print "$%d$ & $%.1f$ & $%.1f$ & $%.1f-%.1f$ & $%.1f-%.1f$" % (iteration, mean, sqrt(var), stats.t.interval(0.95, n - 1, loc=mean, scale=sqrt(var) / sqrt(n))[0], 
    #                                                              stats.t.interval(0.95, n - 1, loc=mean, scale=sqrt(var) / sqrt(n))[1],
    #                                                              sqrt((n - 1)*var/chisquare[1]), 
    #                                                               sqrt((n - 1)*var/chisquare[0]))
    print "$%d$ & $%d$ & $%d$ & $%.1f$ & $%.1f$ & $%.1f-%.1f$ & $%.1f-%.1f$ \\\\" % (iteration[0], iteration[1], 
                                                                                       iteration[2], mean, sqrt(var), t_test[0], t_test[1], sqrt((n - 1)*var/chisquare[1]), sqrt((n - 1)*var/chisquare[0]))
    fileHandle.close();
    return p1

def getShortestDist(pointX, segment):
    u,v = getUV(segment, pointX) 
    if u<1 and u>0:
        return abs(v)
    elif u<0:
        px = segment.P - pointX
        return sqrt(np.vdot(px, px))
    else:
        qx = segment.Q - pointX
        return sqrt(np.vdot(qx, qx))

  def Intro_Crack_xyz(self,k1,p1,p2,q1,q2,u1,u2):
      from numpy.lib import scimath as SM
      filelist = glob.glob("./xyz/*")
      filename = filelist[-1]
      x, y, z , AtomNumber, xlo, xhi, ylo, yhi, zlo, zhi = self.read_xyz(filename)
      xnew,ynew,znew = [],[],[]
      xc ,yc = 0 ,0
      AtomNumber = len(x)
      for i in range(AtomNumber):
          xnew.append(float(x[i])) ; ynew.append(float(y[i])); znew.append(float(z[i]))
      for i in range(len(xnew)):
          dx , dy = xnew[i] - xc , ynew[i] - yc
          r = np.sqrt(dx * dx + dy * dy)
          coeff = 100 * k1 * np.sqrt(2 * r / self.pi)
          if ynew[i] < yc:
              theta = np.arctan2(-dy, dx)
              ux  = 1./(u1-u2) * (u1 * p2 * SM.sqrt(np.cos(theta) - u2 * np.sin(theta)) - u2 * p1 * SM.sqrt(np.cos(theta) - u1 * np.sin(theta)))
              uy  = 1./(u1-u2) * (u1 * q2 * SM.sqrt(np.cos(theta) - u2 * np.sin(theta)) - u2 * q1 * SM.sqrt(np.cos(theta) - u1 * np.sin(theta)))
              ux  = coeff * ux.real
              uy  = coeff * uy.real
          if ynew[i] >= yc:
              theta = np.arctan2(dy, dx)
              ux  = 1./(u1-u2) * (u1 * p2 * SM.sqrt(np.cos(theta) - u2 * np.sin(theta)) - u2 * p1 * SM.sqrt(np.cos(theta) - u1 * np.sin(theta)))
              uy  = 1./(u1-u2) * (u1 * q2 * SM.sqrt(np.cos(theta) - u2 * np.sin(theta)) - u2 * q1 * SM.sqrt(np.cos(theta) - u1 * np.sin(theta)))
              ux  = coeff * ux.real
              uy  = -coeff * uy.real
          print ux
          print uy
          xnew[i] += ux
          ynew[i] += uy
      xlo , xhi, ylo, yhi = min(xnew) - 4 , max(xnew) + 4 , min(ynew) - 4 , max(ynew) + 4
 #     import matplotlib.pylab as plt
 #     from mpl_toolkits.mplot3d import Axes3D
 #     fig = plt.figure()
 #     ax = fig.gca(projection='3d')
 #     ax.scatter(xnew,ynew,znew)
 #     plt.show()
      AtomNumberOut = len(xnew)
  	filename = "Crack.txt"
      with open(filename,mode="w") as fout:
          fout.write("#Edge_dislocation for bcc [111](110)\n")
          fout.write("\n")
          fout.write("%d atoms\n"%(AtomNumberOut))
          fout.write("1 atom types\n")
          fout.write("%f\t%f xlo xhi\n"%(xlo,xhi))
          fout.write("%f\t%f ylo yhi\n"%(ylo,yhi))
          fout.write("%f\t%f zlo zhi\n"%(zlo,zhi))
          fout.write("xy xz yz = %8.5f %8.5f %8.5f"%(0,0,0))
          fout.write("\n")
          fout.write("Atoms\n")
          fout.write("\n")
          for i in range(AtomNumberOut):
              fout.write("%d  1  %12.7f %12.7f %12.7f\n"%(i+1, xnew[i],ynew[i] ,znew[i]))
      fout.close()
      return

def num_divisors(num):
    divisors = 0
    for i in range(1, int(sqrt(num))):
        if num % i == 0:
            divisors += 1
            
    divisors *= 2
    if sqrt(num).is_integer():
        divisors += 1
        
    return divisors

def CF_gz_AC(parms, tau, wixi=wixi):
    u""" Axial (1D) diffusion in a TIR-FCS setup.
        From Two species (bound/unbound) this is the cross correlation part.
        This function is called by other functions within this module.

        *parms* - a list of parameters.
        Parameters (parms[i]):
        [0] D_3D     3D Diffusion coefficient (species A)
        [1] D_2D     2D Diffusion coefficient of bound species C
        [2] sigma    lateral size of the point spread function
                     sigma = simga_0 * lambda / NA
        [3] a        side size of the square pinhole
        [4] d_eva    evanescent decay length (decay to 1/e)
        [5] Conc_3D  3-dimensional concentration of species A
        [6] Conc_2D  2-dimensional concentration of species C
        [7] eta_3D   molecular brightness of species A
        [8] eta_2D   molecular brightness of species C
        [9] k_a      Surface association rate constant
        [10] k_d     Surface dissociation rate constant
        *tau* - lag time
    """
    D = parms[0]
    #D_2D = parms[1]
    #sigma = parms[2]
    # a = parms[3]
    d_eva = parms[4]
    Conc_3D = parms[5]      # ligand concentration in solution
    Conc_2D = parms[6]
    eta_3D = parms[7]
    eta_2D = parms[8]
    k_a = parms[9]
    k_d = parms[10]
    # Define some other constants:
    K = k_a/k_d             # equilibrium constant
    Beta = 1/(1 + K*Conc_3D)
    Re = D / d_eva**2
    Rt = D * (Conc_3D / (Beta * Conc_2D))**2
    Rr = k_a * Conc_3D + k_d
    # Define even more constants:
    sqrtR1 = -Rr/(2*nps.sqrt(Rt)) + nps.sqrt( Rr**2/(4*Rt) - Rr )
    sqrtR2 = -Rr/(2*nps.sqrt(Rt)) - nps.sqrt( Rr**2/(4*Rt) - Rr )
    R1 = sqrtR1 **2
    R2 = sqrtR2 **2
    # And even more more:
    sqrtR3 = sqrtR1 + nps.sqrt(Re)
    sqrtR4 = sqrtR2 + nps.sqrt(Re)
    #R3 = sqrtR3 **2
    #R4 = sqrtR4 **2
    # Calculate return function
    A1 = eta_2D * Conc_2D * k_d / (eta_3D * Conc_3D)
    A2 = sqrtR4*wixi(-nps.sqrt(tau*R1)) - sqrtR3*wixi(-nps.sqrt(tau*R2))
    A3 = ( sqrtR1 - sqrtR2 ) * wixi( nps.sqrt(tau*Re) )
    A4 = ( sqrtR1 - sqrtR2 ) * sqrtR3 * sqrtR4
    Solution = A1 * ( A2 + A3 ) / A4
    # There are some below numerical errors-imaginary numbers.
    # We do not want them.
    return np.real_if_close(Solution)

def init():
    width = 640
    height = 640
    glfw.Init()
    glfw.OpenWindow(width, height, 8, 8, 8, 0, 24, 0, glfw.WINDOW)
    glfw.SetWindowTitle("glfw line")
    glfw.SetWindowSizeCallback(Reshape)
    glClearColor(0.0, 0.0, 0.0, 0.0);
    glEnable(GL_DEPTH_TEST)
    glDepthFunc(GL_LEQUAL);
    
    glfw.SetMouseButtonCallback(MouseHandler)
    glfw.SetKeyCallback(KeyboardHandler)
    glfw.SetWindowCloseCallback(WindowCLose)
    
    global heightmap, terrain_size
    tmpheightmap = open("heightmap.raw", 'rb').read()
    for c in tmpheightmap:
        heightmap.append(ord(c))
    print heightmap
    terrain_size = int(sqrt(len(heightmap)))
    print terrain_size
    
    if(glFogCoordfEXT):  # test
        print 'ok'

def montar_quadrados(tam_quadrados=[]):
    qd = []
    for i in range(len(tam_quadrados)):
        qd.append(figuras.Quadrado(tam_quadrados[i], color(i)))

    if len(qd) < 2:
        print ("ERRO: A lista de tamanho de quadrados tem que ser maior ou igual a 2.")
        return ERRO
    elif len(qd) >= 2:
        # Criar as imagens de todos os quadrados dados.
        for quad_i in qd:
            img, draw = desenhar_quadrado(quad_i)
            nome_arquivo = salvar_arquivo(quad_i.lado_cm, img)
            quad_i.nome_arquivo = nome_arquivo

        x_risco = calcula_x_risco(qd[1], qd[0])
        img1 = Image.open(qd[1].nome_arquivo)
        draw1 = ImageDraw.Draw(img1)
        riscar_quadrado(qd[1], img1, draw1, x_risco)
        salvar_arquivo(qd[1].lado_cm, img1, riscado=1)

        # Novo quadrado montado
        n_area = qd[0].lado_px * qd[1].lado_px
        n_lado = int(sqrt(n_area))

    #        x_risco = calcula_x_risco(qd[1],qd[0])
    #        riscar_quadrado(qd[1],img1,draw1,x_risco)
    #        salvar_arquivo(qd[1].lado_cm, img1,riscado=1)
    #

    return qd

 def readSamples(self, fileName, key,recalc=False,samples=None):
     fn = fileName + ".pre"
     try:
         if recalc: raise IOError()
         with open(fn): pass
         print "precalculated file present"
         self.mu, self.cov = hsplit(mat(fromfile(fn).reshape((3,-1))),[1])
     except IOError:
         if samples != None:
             self._samples = samples
             print "got samples: " , self._samples
         else:
             print "no file present, calculating..."
             smpls = loadmat(fileName)[key]
             print "loaded from mat file"
             self._samples = mat(smpls)
             print "reshaped into samples"
         self.mu = sum(self._samples, axis=1) / self._samples.shape[1]
         print "mu=", str(self.mu)
         sampdiffmu = self._samples - self.mu
         self.cov = sampdiffmu*sampdiffmu.T / self._samples.shape[1]
         print"cov=", str(self.cov)
         mat(hstack((self.mu,self.cov))).tofile(fn)
     self._invCov = self.cov.I
     self._detCov = det(self.cov)
     self._multConst = 1 / sqrt((2 * pi) ** 3 * self._detCov)

def R_F(x, y, z, rtol=3e-4):
    r"""Computes the symmetric integral of the first kind

    .. math:: R_F(x, y, z) =
              \frac{1}{2}\int_0^\infty [(t + x)(t + y)(t + z)]^\frac{1}{2}dt

    """
    # Properly deal with the degenerate case
    if y == z:
        return R_C(x, y)

    A0 = (x + y + z) / 3.0
    M = max(abs(A0 - x), abs(A0 - y), abs(A0 - z))
    Q = ((3*rtol) ** (-1.0/6.0)) * M

    A = _A_gen(x, y, z)
    for n, Am in enumerate(A):
        if 4**(-n) * Q < abs(Am):
            X = (A0 - x) / (Am * 4**n)
            Y = (A0 - y) / (Am * 4**n)
            Z = -X-Y
            E2 = X * Y - Z**2
            E3 = X * Y * Z
            result = (1/scimath.sqrt(Am)) * (1 - E2*(1/10) + E3*(1/14)
                                             + (E2**2)*(1/24) - (E2*E3)*(3/44))
            return result

def bessel_series(a, b, z, maxiters=500, tol=tol):
    """Compute hyp1f1 using a series of Bessel functions; see (3.20) in
    _[pop].

    """
    Do, Dm, Dn = 1, 0, b/2
    r = sqrt(z*(2*b - 4*a))
    w = z**2 / r**(b + 1)
    # Ao is the first coefficient, and An is the second.
    Ao = 0
    # The first summand comes from the zeroth term.
    So = Do*jv(b - 1, r) / r**(b - 1) + Ao
    An = Dn*w*jv(b + 1, r)
    Sn = So + An
    i = 3
    while i <= maxiters:
        if np.abs(Ao/So) < tol and np.abs(An/Sn) < tol:
            break
        tmp = Dn
        Dn = ((i - 2 + b)*Dm + (2*a - b)*Do) / i
        Do = Dm
        Dm = tmp
        w *= z/r
        Ao = An
        An = Dn*w*jv(b - 1 + i, r)
        So = Sn
        Sn += An
        i += 1
    # if i > maxiters:
    #     warnings.warn("Number of evaluations exceeded maxiters on "
    #                   "a = {}, b = {}, z = {}.".format(a, b, z))
    return gamma(b)*np.exp(z/2)*2**(b - 1)*Sn

def normalizeData(data, meanOnly = False):
    """
    normalize data by subtracting mean and dividing by sd per COLUMN
    @param data: an array
    @param meanOnly: if True subtract mean only; otherwise divide by sd too
    @return: (an array with the same dimension as data, mean, stds, the transformer)
    """

    # compute the new data
    m = mean(data, axis=0)
    res = data - m

    if meanOnly:
        stds = 1
    else:
        stds = sqrt(var(data, axis=0))
        stds[stds==0] = 1   # to avoid dividing by 0

        res /= stds

    # figure out the transformer
    def foo(givenData):
        assert givenData.shape[1]==data.shape[1], "Only arrays of %d columns are handled." % data.shape[1]
        return (givenData - m)/stds

    return res, m, stds, foo

def rmsd(p1, p2):
    if len(p1) != len(p2):
        raise Exception("Degree of points must be the same")
    total = 0
    for t in range(len(p1)):
        d = (p1[t] - p2[t])
        total += d * d
    return sqrt(total / len(p1)) #todo remove sqrt?

    def simplify( self ):
        devisors = self.get_possible_divisors()

        for devisor in reversed( devisors ):
            if self.number % devisor == 0:
                self.number /= devisor
                self.multiplier *= sqrt( devisor ) if self.order == 2 else pow( devisor, 1.0 / self.order )
                break

def generate_hamiltonian(n,m,n0,q,c):
    n_states = int((n-m)/2)
    ham = np.zeros((n_states,n_states),dtype = complex)
    for i in range(n_states):
        n0 = (n-m-2*i)
        #first for diagonal
        ham[i,i] = c/2*(m**2 + n + n0 + 2*n*n0-2*n0**2)+q*(n-n0)
        #now for off diagonal- use try statements to catch edge case
        try:
            ham[i,i+1] = c/2*scimath.sqrt(n0*(n0-1)*(n+m-n0+2)*(n-m-n0+2))
        except:
            pass
        try:
            ham[i,i-1] = c/2*scimath.sqrt((n0+1)*(n0+2)*(n+m-n0)*(n-m-n0))
        except:
            pass
    return ham

 def sigma(self, r):
     """Compute the surface mass density of the halo at distance r
     (in Mpc) from the halo center."""
     r = unit_checker(r, u.Mpc)
     x = r / self.r_s
     val1 = 1 / (x**2 - 1)
     val2 = (arcsec(x) / (sm.sqrt(x**2 - 1))**3).real
     return 2 * self.r_s * self.rho_c * self.delta_c * (val1-val2)