def lnmixgaussian(x, params):
    """
    NAME:
       lnmixgaussian
    PURPOSE:
       returns the log of a mixture of two two-dimensional gaussian
    INPUT:
       x - 2D point to evaluate the Gaussian at
       params - mean and variances ([mean_array,inverse variance matrix, mean_array, inverse variance, amp1])
    OUTPUT:
       log N(mean,var)
    HISTORY:
       2009-10-30 - Written - Bovy (NYU)
    """
    return sc.log(
        params[4]
        / 2.0
        / sc.pi
        * sc.sqrt(linalg.det(params[1]))
        * sc.exp(-0.5 * sc.dot(x - params[0], sc.dot(params[1], x - params[0])))
        + (1.0 - params[4])
        / 2.0
        / sc.pi
        * sc.sqrt(linalg.det(params[3]))
        * sc.exp(-0.5 * sc.dot(x - params[2], sc.dot(params[3], x - params[2])))
    )

def gabor2d(gsw, gsh, gx0, gy0, wfreq, worient, wphase, shape):
    """ Generate a gabor 2d array
    
    Inputs:
      gsw -- standard deviation of the gaussian envelope (width)
      gsh -- standard deviation of the gaussian envelope (height)
      gx0 -- x indice of center of the gaussian envelope
      gy0 -- y indice of center of the gaussian envelope
      wfreq -- frequency of the 2d wave
      worient -- orientation of the 2d wave
      wphase -- phase of the 2d wave
      shape -- shape tuple (height, width)

    Outputs:
      gabor -- 2d gabor with zero-mean and unit-variance

    """
    
    height, width = shape
    y, x = N.mgrid[0:height, 0:width]
    
    X = x * N.cos(worient) * wfreq
    Y = y * N.sin(worient) * wfreq
	
    env = N.exp( -.5 * ( ((x-gx0)**2./gsw**2.) + ((y-gy0)**2./gsh**2.) ) )
    wave = N.exp( 1j*(2*N.pi*(X+Y) + wphase) )
    gabor = N.real(env * wave)
    
    gabor -= gabor.mean()
    gabor /= fastnorm(gabor)
    
    return gabor

    def coste(self, *args, **kwargs):
        """
        material:
            0   -   Carbon steel
            1   -   Stainless steel 316
            2   -   Stainless steel 304
            3   -   Stainless steel 347
            4   -   Nickel
            5   -   Monel
            6   -   Inconel
            7   -   Zirconium
            8    -  Titanium
            9    -   Brick and rubber or brick and polyester-lined steel
            10  -   Rubber or lead-lined steel
            11  -   Polyester, fiberglass-reinforced
            12  -   Aluminum
            13  -   Copper
            14  -   Concrete
        """
        self._indicesCoste(*args)
        
        self.material=kwargs["material"]
        
        V=self.Volumen.galUS

        Fm=[1., 2.7, 2.4, 3.0, 3.5, 3.3, 3.8, 11.0, 11.0, 2.75, 1.9, 0.32, 2.7, 2.3, 0.55][self.material]

        if V<=21000:
            C=Fm*exp(2.631+1.3673*log(V)-0.06309*log(V)**2)
        else:
            C=Fm*exp(11.662+0.6104*log(V)-0.04536*log(V)**2)
        
        self.C_adq=Currency(C*self.Current_index/self.Base_index)
        self.C_inst=Currency(self.C_adq*self.f_install)

    def rekernel(self):
        dcontrol = self.control[ord('d')]
        econtrol = self.control[ord('e')]
        rcontrol = self.control[ord('r')]
        radius = rcontrol.val
        dvalue = dcontrol.val
        evalue = econtrol.val
        rmax   = rcontrol.limit[2]
        if self.rlast != radius:
            inner, outer = float(radius-1), float(radius)
            shape      = (self.edge, self.edge)
            self.radii = list(product(arange(-rmax,rmax+1,1.0), repeat=2))
            self.radii = array([sqrt(x*x+y*y) for x,y in self.radii]).reshape(shape)

            if True:
                self.negative = -exp(-dvalue*(self.radii-outer)**2)
                self.positive = +exp(-dvalue*(self.radii-inner)**2)
            else:
                self.radii = around(self.radii)
                self.negative = zeros((self.edge,self.edge),dtype=float)
                self.negative[self.radii == outer] = -1.0
                self.positive = zeros(shape,dtype=float)
                self.positive[self.radii == inner] = +1.0

            self.negative /= fabs(self.negative.sum())
            self.positive /= fabs(self.positive.sum())

            self.kernel = self.negative + self.positive
            self.rlast = radius
        if self.elast != evalue:
            self.gauss = exp(-evalue * self.radii**2)
            self.gauss /= self.gauss.sum()
            self.elast = evalue

def morlet(M, w=5.0, s=1.0, complete=True):
    """
    Complex Morlet wavelet.

    Parameters
    ----------
    M : int
        Length of the wavelet.
    w : float, optional
        Omega0. Default is 5
    s : float, optional
        Scaling factor, windowed from ``-s*2*pi`` to ``+s*2*pi``. Default is 1.
    complete : bool, optional
        Whether to use the complete or the standard version.

    Returns
    -------
    morlet : (M,) ndarray

    See Also
    --------
    scipy.signal.gausspulse

    Notes
    -----
    The standard version::

        pi**-0.25 * exp(1j*w*x) * exp(-0.5*(x**2))

    This commonly used wavelet is often referred to simply as the
    Morlet wavelet.  Note that this simplified version can cause
    admissibility problems at low values of `w`.

    The complete version::

        pi**-0.25 * (exp(1j*w*x) - exp(-0.5*(w**2))) * exp(-0.5*(x**2))

    This version has a correction
    term to improve admissibility. For `w` greater than 5, the
    correction term is negligible.

    Note that the energy of the return wavelet is not normalised
    according to `s`.

    The fundamental frequency of this wavelet in Hz is given
    by ``f = 2*s*w*r / M`` where `r` is the sampling rate.
    
    Note: This function was created before `cwt` and is not compatible
    with it.

    """
    x = linspace(-s * 2 * pi, s * 2 * pi, M)
    output = exp(1j * w * x)

    if complete:
        output -= exp(-0.5 * (w**2))

    output *= exp(-0.5 * (x**2)) * pi**(-0.25)

    return output

def f_Refl_Thin_Film_fit(Data):
    strain = Data[0]
    DW = Data[1]

    dat = Data[2]
    wl = dat[0]
    N = dat[2]

    phi = dat[3]
    t_l = dat[4]
    thB_S = dat[6]
    G = dat[7]
    F0 = dat[8]
    FH = dat[9]
    FmH = dat[10]
    th = dat[19]

    thB = thB_S - strain * tan(thB_S)  # angle de Bragg dans chaque lamelle

    eta = 0
    res = 0

    n = 1
    while (n <= N):
        g0 = sin(thB[n] - phi)  # gamma 0
        gH = -sin(thB[n] + phi)  # gamma H
        b = g0 / gH
        T = pi * G * ((FH[n]*FmH[n])**0.5) * t_l * DW[n] / (wl * (abs(g0*gH)**0.5))
        eta = (-b*(th-thB[n])*sin(2*thB_S) - 0.5*G*F0[n]*(1-b)) / ((abs(b)**0.5) * G * DW[n] * (FH[n]*FmH[n])**0.5)
        S1 = (res - eta + (eta*eta-1)**0.5)*exp(-1j*T*(eta*eta-1)**0.5)
        S2 = (res - eta - (eta*eta-1)**0.5)*exp(1j*T*(eta*eta-1)**0.5)
        res = (eta + ((eta*eta-1)**0.5) * ((S1+S2)/(S1-S2)))
        n += 1
    return res

def Alfa(Tr, m, f_acent, EOS="SRK", alfa=None):
    """Diferentes expresiónes de cálculo de la alfa de la ecuaciones de estado de SRK y PR
    alfa: Indica el método de cálculo de alfa
        0 - Original
        1 - Boston Mathias
        2 - Twu
        3 - Doridon
        Boston, J.F.; Mathias, P.M. Phase Equilibria in a Third-Generation Process Simulator. Proc. 2nd. Int. Conf. On Phase Equilibria and Fluid Properties in the Chemical Process Industries, Berlin, Germany 17.-21.3.1980, p. 823.
        """
    if not alfa:
        Config=config.getMainWindowConfig()
        alfa=Config.getint("Thermo","Alfa")

    if alfa==2: #Función alfa de Twu et Alt.
        if EOS=="PR":
            if Tr>1:
                L0=0.401219
                M0=4.963075
                N0=-0.2
                L1=0.024655
                M1=1.248088
                N1=-8.
            else:
                L0=0.125283
                M0=0.911807
                N0=1.948153
                L1=0.511614
                M1=0.784054
                N1=2.812522
        else:
            if Tr>1:
                L0=0.441411
                M0=6.500018
                N0=-0.2
                L1=0.032580
                M1=1.289098
                N1=-8.
            else:
                L0=0.141599
                M0=0.919422
                N0=2.496441
                L1=0.500315
                M1=0.799457
                N1=3.29179
        alfa0=Tr**(N0*(M0-1))*exp(L0*(1-Tr**(N0*M0)))
        alfa1=Tr**(N1*(M1-1))*exp(L1*(1-Tr**(N1*M1)))
        alfa=alfa0+f_acent*(alfa1-alfa0)
    elif alfa==3:
        if f_acent<0.4:
            m=0.418+1.58*f_acent-0.58*f_acent**2
        else:
            m=0.212+2.2*f_acent-0.831*f_acent**2
        alfa=exp(m*(1-Tr))
    elif alfa==1 and Tr>1:
        d=1.+m/2.
        c=1.-1./d
        alfa=exp(c*(1-Tr**d))**2
    else:
        alfa=(1+m*(1-Tr**0.5))**2
    return alfa

    def _thermo0(self, rho, T, fase):
        GT = [-2.903423528e5, 4.680624952e5, -1.8954783215e5, -4.8262235392e3,
              2.243409372e4, -6.6206354818e3, 8.9937717078e2, -6.0559143718e1,
              1.6370306422]
        lo = 0
        for i in range(-3, 6):
            lo += GT[i+3]*T**(i/3.)

        tita = (rho.gcc-0.221)/0.221
        j = [0, -1.304503323e1, 1.8214616599e1, 9.903022496e3, 7.420521631e2,
             -3.0083271933e-1, 9.6456068829e1, 1.350256962e4]
        l1 = exp(j[1]+j[4]/T)*(exp(rho.gcc**0.1*(j[2]+j[3]/self.T**1.5)+tita*rho.gcc**0.5*(j[5]+j[6]/T+j[7]/T**2))-1.)

        lc = 0
        # FIXME: no sale
#        deltarho=(self.rho/self.M-0.221)/0.221
#        deltaT=(self.T-282.34)/282.34
#        xkt=(1.0/self.rho/self.M/self.derivative("P", "rho", "T")*1e3)**0.5
#        b=abs(deltarho)/abs(deltaT)**1.19
#        xts=(self.rho/self.M)**2*xkt*5.039/.221**2
#        g=xts*abs(deltaT)**1.19
#        xi=0.69/(b**2*5.039/g/Boltzmann/282.34)**0.5
#        f=exp(-18.66*deltaT**2-4.25*deltarho**4)
#        c=(self.M/self.rho.gcc/Avogadro/Boltzmann/self.T)**0.5
#        lc=c*Boltzmann*self.T**2/6.0/pi/self.mu.muPas/xi*self.dpdT**2*self.kappa**0.5*f
#        print lo, l1
        return unidades.ThermalConductivity(lo+l1+lc, "mWmK")

def objective(pars,Z,ycovar):
    """The objective function"""
    bcost= pars[0]
    t= pars[1]
    Pb= pars[2]
    Xb= pars[3]
    Yb= pars[4]
    Zb= sc.array([Xb,Yb])
    Vb1= sc.exp(pars[5])
    Vb2= sc.exp(pars[6])
    corr= pars[7]
    V= sc.array([[Vb1,sc.sqrt(Vb1*Vb2)*corr],[sc.sqrt(Vb1*Vb2)*corr,Vb2]])
    v= sc.array([-sc.sin(t),sc.cos(t)])
    if Pb < 0. or Pb > 1.:
        return -sc.finfo(sc.dtype(sc.float64)).max
    if corr < -1. or corr > 1.:
        return -sc.finfo(sc.dtype(sc.float64)).max
    delta= sc.dot(v,Z.T)-bcost
    sigma2= sc.dot(v,sc.dot(ycovar,v))

    ndata= Z.shape[0]
    detVycovar= sc.zeros(ndata)
    deltaOUT= sc.zeros(ndata)
    for ii in range(ndata):
        detVycovar[ii]= m.sqrt(linalg.det(V+ycovar[:,ii,:]))
        deltaOUT[ii]= sc.dot(Z[ii,:]-Zb,sc.dot(linalg.inv(V+ycovar[:,ii,:]),Z[ii,:]-Zb))
    return sc.sum(sc.log((1.-Pb)/sc.sqrt(2.*m.pi*sigma2)*
                         sc.exp(-0.5*delta**2./sigma2)
                         +Pb/2./m.pi/detVycovar
                         *sc.exp(-0.5*deltaOUT)))

 def gamma(xAg):
     rA=3.19;rB=2.11
     qA=2.4;qB=1.97
     qAdash=2.4;qBdash=0.89
     aAB=242.53;aBA=-75.13
     p11=(xAg*rA/(xAg*rA+xBg*rB))
     p12=(xBg*rB/(xAg*rA+xBg*rB))
     p1=scipy.array([p11,p12])
     t11=(xAg*qA/(xAg*qA+xBg*qB))
     t12=(xBg*qB/(xAg*qA+xBg*qB))
     t1=scipy.array([t11,t12])
     t21=(xAg*qAdash/(xAg*qAdash+xBg*qBdash))
     t22=(xBg*qBdash/(xAg*qAdash+xBg*qBdash))
     t2=scipy.array([t21,t22])
     tAB=scipy.exp(-aAB/Tg)
     tBA=scipy.exp(-aBA/Tg)
     t=scipy.array([tAB,tBA])
     z=10
     l11=((z/2)*(rA-qA))-(rA-1)
     l12=((z/2)*(rB-qB))-(rB-1)
     l1=scipy.array([l11,l12])
     gamma1=scipy.exp((scipy.log(p1[0]/xAg))+(z*qA*scipy.log(t1[0]/p1[0])/2)+(p1[1]*(l1[0]-(rA*l1[1]/rB)))-(qAdash*scipy.log(t2[0]+t2[1]*t[1]))+(t2[1]*qAdash*((t[1]/(t2[0]+t2[1]*t[1]))-(t[0]/(t2[1]+t2[0]*t[0])))))
     gamma2=scipy.exp((scipy.log(p1[1]/xBg))+(z*qB*scipy.log(t1[1]/p1[1])/2)+(p1[0]*(l1[1]-(rB*l1[0]/rA)))-(qBdash*scipy.log(t2[1]+t2[0]*t[0]))+(t2[0]*qBdash*((t[0]/(t2[1]+t2[0]*t[0]))-(t[1]/(t2[0]+t2[1]*t[1])))))
     gamma=scipy.array([gamma1,gamma2])
     return gamma

def test_time_integration():

    dt = 0.01
    tmax = 1

    def check_time_integration(res, y0, exact):
        ts, ys = time_integrate(res, y0, dt, tmax)

        exacts = [exact(t) for t in ts]

        print(ts, ys, exacts)

        utils.assert_list_almost_equal(ys, exacts, 1e-3)

    tests = [
        (lambda t, y, dy: y - dy, 1.0, lambda t: array(exp(t))),
        (lambda t, y, dy: y + dy, 1.0, lambda t: array(exp(-t))),
        (
            lambda t, y, dy: array([-0.1 * sin(t), y[1]]) - dy,
            array([0.1 * cos(0.0), exp(0.0)]),
            lambda t: array([0.1 * cos(t), exp(t)]),
        ),
    ]

    for r, y0, exact in tests:
        yield check_time_integration, r, y0, exact

def molar_conc(mole_fraction,ihenry_k,temp,excess_air,altitude_meters,salinity):
    '''Converts mole fraction to moles per liter water.'''
    
    barometric_pressure = exp(-(altitude_meters/8300)) 
    vapor_pressure_h2o = exp(24.4543 - 67.4509*(100/temp) - 4.8489*log(temp/100) - 0.000544*salinity)
    
    return (ihenry_k*mole_fraction*(barometric_pressure - vapor_pressure_h2o) + (excess_air*mole_fraction)/22414)

def time_plot(m=10, c=1, k=100, x0=1, v0=-1, max_time=100):
    t, x, v, zeta, omega, omega_d, A = free_response(
        m, c, k, x0, v0, max_time)
    fig = plt.figure()
    fig.suptitle('Displacement vs Time')
    ax = fig.add_subplot(111)
    ax.set_xlabel('Time')
    ax.set_ylabel('Displacement')
    ax.grid('on')
    ax.plot(t, x)
    if zeta < 1:
        ax.plot(t, A * sp.exp(-zeta * omega * t), '--', t, -A *
                sp.exp(-zeta * omega * t), '--g', linewidth=1)
        tmin, tmax, xmin, xmax = ax.axis()
        ax.text(.75 * tmax, .95 * (xmax - xmin) + xmin,
                '$\omega$ = %0.2f rad/sec' % (omega))
        ax.text(.75 * tmax, .90 * (xmax - xmin) +
                xmin, '$\zeta$ = %0.2f' % (zeta))
        ax.text(.75 * tmax, .85 * (xmax - xmin) + xmin,
                '$\omega_d$ = %0.2f rad/sec' % (omega_d))
    else:
        tmin, tmax, xmin, xmax = ax.axis()
        ax.text(.75 * tmax, .95 * (xmax - xmin) +
                xmin, '$\zeta$ = %0.2f' % (zeta))
        ax.text(.75 * tmax, .90 * (xmax - xmin) + xmin,
                '$\lambda_1$ = %0.2f' % (zeta * omega - omega * (zeta ** 2 - 1)))
        ax.text(.75 * tmax, .85 * (xmax - xmin) + xmin,
                '$\lambda_2$ = %0.2f' % (zeta * omega + omega * (zeta ** 2 - 1)))

def diff_gauss_kern(size = 2, sigmA = 0.9, sigmB = 1., ratio = 1.):
	""" Returns a normalized 2D kernel of a gaussian difference """
	x, y = sp.mgrid[-size:size+1, -size:size+1]
	sigmA = float(sigmA**2)
	sigmB = float(sigmB**2)
	g = sp.exp( -(x**2 + y**2) / sigmA ) - ratio * sp.exp( -(x**2 + y**2) / sigmB )
	return g / abs(g.sum())

    def _mur(self, rho, T, fase):
        # Modified friction theory for residual viscosity contribution

        Gamma = self.Tc/T
        psi1 = exp(Gamma)                                               # Eq 15
        psi2 = exp(Gamma**2)                                            # Eq 16

        a = [68.9659e-6, -22.0494e-6, -42.6126e-6]
        b = [153.406e-6, 8.45198e-6, -113.967e-6]
        A = [0.78238e-9, -0.64717e-9, 1.39066e-9]
        B = [-9.75792e-9, -3.19303e-9, 12.4263e-9]

        ka = (a[0] + a[1]*psi1 + a[2]*psi2) * Gamma                     # Eq 11
        kr = (b[0] + b[1]*psi1 + b[2]*psi2) * Gamma                     # Eq 12
        kaa = (A[0] + A[1]*psi1 + A[2]*psi2) * Gamma                    # Eq 13
        krr = (B[0] + B[1]*psi1 + B[2]*psi2) * Gamma                    # Eq 14

        # All parameteres has pressure units of bar
        Patt = -fase.IntP.bar
        Prep = T*fase.dpdT_rho.barK
        Pid = rho*self.R*self.T/1e5
        delPr = Prep-Pid

        # Eq 5
        mur = kr*delPr + ka*Patt + krr*Prep**2 + kaa*Patt**2
        return mur*1e3

def superimpose (cavityD,cavityU,par1,par2):
   if (par2 == 0) :
     return (sp.cos(par1*sp.pi)*cavityD[:] + 1j*sp.sin(par1*sp.pi)*cavityU[:])*sp.exp(-1j*par1*sp.pi)
   elif (par2 == 1) :
     return sp.exp(-1j*par1*sp.pi)*cavityD
   else :
     return cavityD

    def evaluate_at(self, grid, component=None, prefactor=False):
        r"""Evaluate the Hagedorn wavepacket :math:`\Psi` at the given nodes :math:`\gamma`.

        :param grid: The grid :math:`\Gamma` containing the nodes :math:`\gamma`.
        :type grid: A class having a :py:meth:`get_nodes(...)` method.
        :param component: The index :math:`i` of a single component :math:`\Phi_i` to evaluate.
                          (Defaults to ``None`` for evaluating all components.)
        :param prefactor: Whether to include a factor of :math:`\frac{1}{\sqrt{\det(Q)}}`.
        :type prefactor: Boolean, default is ``False``.
        :return: A list of arrays or a single array containing the values of the :math:`\Phi_i` at the nodes :math:`\gamma`.
        """
        Pis = self.get_parameters(component=component, aslist=True)

        if component is not None:
            phase = exp(1.0j * Pis[component][4] / self._eps**2)
            values = phase * self.slim_recursion(grid, component, prefactor=prefactor)

        else:
            values = []

            for component in range(self._number_components):
                # Note: This is very inefficient! We may evaluate the same basis functions multiple
                #       times. But as long as we don't know that the basis shapes are true subsets
                #       of the largest one, we can not evaluate just all functions in this
                #       maximal set.

                # TODO: Find more efficient way to do this

                phase = exp(1.0j * Pis[component][4] / self._eps**2)
                values.append(phase * self.slim_recursion(grid, component, prefactor=prefactor))

        return values

def dbexpl(p):
    t=arange(0,100,20.)
    #if (p['par1']) < 0.5:
    if (p['par1']+p['par3']) < 0.25:
        asdf
    y =  (p['par1']*exp(-p['par2']*t) + p['par3']*exp(-p['par4']*t))
    return y

 def f(self,xarr,t):
     x0dot = -self.coef*pl.exp(-self.k*(1.0+abs(xarr[3]-1.0)))*pl.sin(self.w*t-self.k*xarr[2])
     x1dot = pl.sign(xarr[3]-1.0)*self.coef*pl.exp(-self.k*(1.0+abs(xarr[3]-1.0)))*pl.cos(self.w*t-self.k*xarr[2]) -\
             pl.sign(xarr[3]-1.0)*9.8
     x2dot = xarr[0]
     x3dot = xarr[1]
     return [x0dot,x1dot,x2dot,x3dot]

    def _ThCondCritical(self, rho, T, fase):
        # Custom Critical enhancement

        # The paper use a diferent rhoc value to the EoS
        rhoc = 235

        t = abs(T-405.4)/405.4
        dPT = 1e5*(2.18-0.12/exp(17.8*t))
        nb = 1e-5*(2.6+1.6*t)

        DL = 1.2*Boltzmann*T**2/6/pi/nb/(1.34e-10/t**0.63*(1+t**0.5)) * \
            dPT**2 * 0.423e-8/t**1.24*(1+t**0.5/0.7)

        # Add correction for entire range of temperature, Eq 10
        DL *= exp(-36*t**2)

        X = 0.61*rhoc+16.5*log(t)
        if rho > 0.6*rhoc:
            # Eq 11
            DL *= X**2/(X**2+(rho-0.96*rhoc)**2)
        else:
            # Eq 14
            DL = X**2/(X**2+(0.6*rhoc-0.96*rhoc)**2)
            DL *= rho**2/(0.6*rhoc)**2

        return DL