# -*- coding: utf-8 -*-
"""
Solution of the Schrödinger equation using Ritz method
No and Hermite functions as basis.

@author: <NAME>
"""
from __future__ import division, print_function
import numpy as np
from scipy.special import eval_hermite, roots_hermite, factorial
from scipy.linalg import eigh
import matplotlib.pyplot as plt

plt.rcParams["mathtext.fontset"] = "cm"

def gauss_hermite(fun, n, args=None):
    """Gauss hermite quadrature"""
    xi, wi = roots_hermite(n)
    return np.sum(fun(xi, *args) * wi)


def hermite_mag(n):
    """Normalization factor for Hermite polynomials"""
    return np.sqrt(np.pi**0.5 * 2**n * factorial(n))


def hermite_norm(n, x):
    """Normalized Hermite polynomials"""
    return eval_hermite(n, x)/ hermite_mag(n)


def delta(m, n):
    """Kronecker delta"""
    if m == n:
        val = 1
    else:
        val = 0
    return val


def alpha(m, n):
    """Integral of x * psi_m * psi_n
    
    psi_m are the Hermite functions
    """
    val = 0
    if m == n + 1:
        val = np.sqrt((n + 1)/2)
    elif m == n - 1:
        val = np.sqrt(n/2)
    return val


def beta(m, n):
    """Integral of x**2 * psi_m * psi_n
    
    psi_m are the Hermite functions
    """
    val = 0
    if m == n - 2:
        val = np.sqrt(n*(n - 1))/2
    elif m == n:
        val = (2*n + 1)/2
    elif m == n + 2:
        val = np.sqrt((n + 1)*(n + 2))/2
    return val


def potential_integrand(x, m, n, potential):
    """Integrand for the potential functional"""
    return hermite_norm(m, x) * hermite_norm(n, x) * potential(x)


def hamiltonian_mat(num, potential, ngpts):
    """Compute the Hamiltonian for a given potential
    
    The potential integral is computed using ngpts Gauss points and the
    Hermite-Gauss quadrature scheme.
    
    Parameters
    ----------
    num : int
        Number of basis functions to use.
    potential : callable
        Potential function.
    ngpts : int
        Number of Gauss-Hermite points to use for the integral
        of the potential contribution to the Hamiltonian.
    
    Returns
    -------
    Hmat : ``(num, num)`` ndarray
        Discretized version of the Hamiltonian corresponding
        to the subspace generated by the span of the basis
        functions.
    """
    Hmat = np.zeros((num, num))
    for row in range(0, num):
        for col in range(0, num):
            fac_row = hermite_mag(row - 1)/hermite_mag(row)
            fac_col = hermite_mag(col - 1)/hermite_mag(col)
            kin = 4*row*col*delta(row, col)*fac_row**2 -\
                  2*row*alpha(row - 1, col)*fac_row -\
                  2*col*alpha(row, col - 1)*fac_col +\
                  beta(row, col)
            pot = gauss_hermite(potential_integrand, ngpts,
                                      args=(row, col, potential))
            Hmat[row, col] = 0.5*kin + pot
    return Hmat


def eigenstates(potential, nterms=20, ngpts=100):
    """Compute eigenvalues and coefficients for the eigenvectors"""
    Hmat = hamiltonian_mat(nterms, potential, ngpts)
    vals, vecs = eigh(Hmat)
    return vals, vecs


def plot_eigenstates(vals, vecs, potential, lims=(-5,5), npts=1001,
                     ax=None):
    """Plot eigenstates"""
    if ax is None:
        ax =  plt.gca()
    xmin, xmax = lims
    x = np.linspace(xmin, xmax, npts)
    nterms = vecs.shape[0]
    out = ax.plot(x, potential(x), color="black")
    for k in range(10):
        vec = 0
        for j in range(nterms):
            coeff = vecs[j, k]
            vec += coeff * hermite_norm(j, x) * np.exp(-x**2/2)
        vec = vec/np.max(vec)
        amp = vals[k + 1] - vals[k]
        ax.plot(x, 0.3*amp*vec + vals[k], label=r"$n=%i$"%k)

    min_val = min(np.min(potential(x)), 0)
    plt.xlabel(r"$x$")
    plt.ylabel(r"$\psi(x)$")
    plt.ylim(min_val, vals[k] + amp)   
    return out


if __name__ == "__main__":
    potential = lambda x: 0.5*x**2
    vals, vecs = eigenstates(potential)
    print(np.round(vals[:10], 6))
    fig, ax = plt.subplots(1, 1)
    plot_eigenstates(vals, vecs, potential)
    plt.show()

# -*- coding: utf-8 -*-
"""
Solution of the Schrödinger equation using Ritz method
No and Hermite functions as basis.

@author: <NAME>
"""
from __future__ import division, print_function
import numpy as np
from scipy.special import eval_hermite, roots_hermite, factorial
from scipy.linalg import eigh
import matplotlib.pyplot as plt

plt.rcParams["mathtext.fontset"] = "cm"

def gauss_hermite(fun, n, args=None):
    """Gauss hermite quadrature"""
    xi, wi = roots_hermite(n)
    return np.sum(fun(xi, *args) * wi)


def hermite_mag(n):
    """Normalization factor for Hermite polynomials"""
    return np.sqrt(np.pi**0.5 * 2**n * factorial(n))


def hermite_norm(n, x):
    """Normalized Hermite polynomials"""
    return eval_hermite(n, x)/ hermite_mag(n)


def delta(m, n):
    """Kronecker delta"""
    if m == n:
        val = 1
    else:
        val = 0
    return val


def alpha(m, n):
    """Integral of x * psi_m * psi_n
    
    psi_m are the Hermite functions
    """
    val = 0
    if m == n + 1:
        val = np.sqrt((n + 1)/2)
    elif m == n - 1:
        val = np.sqrt(n/2)
    return val


def beta(m, n):
    """Integral of x**2 * psi_m * psi_n
    
    psi_m are the Hermite functions
    """
    val = 0
    if m == n - 2:
        val = np.sqrt(n*(n - 1))/2
    elif m == n:
        val = (2*n + 1)/2
    elif m == n + 2:
        val = np.sqrt((n + 1)*(n + 2))/2
    return val


def potential_integrand(x, m, n, potential):
    """Integrand for the potential functional"""
    return hermite_norm(m, x) * hermite_norm(n, x) * potential(x)


def hamiltonian_mat(num, potential, ngpts):
    """Compute the Hamiltonian for a given potential
    
    The potential integral is computed using ngpts Gauss points and the
    Hermite-Gauss quadrature scheme.
    
    Parameters
    ----------
    num : int
        Number of basis functions to use.
    potential : callable
        Potential function.
    ngpts : int
        Number of Gauss-Hermite points to use for the integral
        of the potential contribution to the Hamiltonian.
    
    Returns
    -------
    Hmat : ``(num, num)`` ndarray
        Discretized version of the Hamiltonian corresponding
        to the subspace generated by the span of the basis
        functions.
    """
    Hmat = np.zeros((num, num))
    for row in range(0, num):
        for col in range(0, num):
            fac_row = hermite_mag(row - 1)/hermite_mag(row)
            fac_col = hermite_mag(col - 1)/hermite_mag(col)
            kin = 4*row*col*delta(row, col)*fac_row**2 -\
                  2*row*alpha(row - 1, col)*fac_row -\
                  2*col*alpha(row, col - 1)*fac_col +\
                  beta(row, col)
            pot = gauss_hermite(potential_integrand, ngpts,
                                      args=(row, col, potential))
            Hmat[row, col] = 0.5*kin + pot
    return Hmat


def eigenstates(potential, nterms=20, ngpts=100):
    """Compute eigenvalues and coefficients for the eigenvectors"""
    Hmat = hamiltonian_mat(nterms, potential, ngpts)
    vals, vecs = eigh(Hmat)
    return vals, vecs


def plot_eigenstates(vals, vecs, potential, lims=(-5,5), npts=1001,
                     ax=None):
    """Plot eigenstates"""
    if ax is None:
        ax =  plt.gca()
    xmin, xmax = lims
    x = np.linspace(xmin, xmax, npts)
    nterms = vecs.shape[0]
    out = ax.plot(x, potential(x), color="black")
    for k in range(10):
        vec = 0
        for j in range(nterms):
            coeff = vecs[j, k]
            vec += coeff * hermite_norm(j, x) * np.exp(-x**2/2)
        vec = vec/np.max(vec)
        amp = vals[k + 1] - vals[k]
        ax.plot(x, 0.3*amp*vec + vals[k], label=r"$n=%i$"%k)

    min_val = min(np.min(potential(x)), 0)
    plt.xlabel(r"$x$")
    plt.ylabel(r"$\psi(x)$")
    plt.ylim(min_val, vals[k] + amp)   
    return out


if __name__ == "__main__":
    potential = lambda x: 0.5*x**2
    vals, vecs = eigenstates(potential)
    print(np.round(vals[:10], 6))
    fig, ax = plt.subplots(1, 1)
    plot_eigenstates(vals, vecs, potential)
    plt.show()

# -*- coding: utf-8 -*- Solution of the Schrödinger equation using Ritz method No and Hermite functions as basis.  @author: <NAME> Gauss hermite quadrature Normalization factor for Hermite polynomials Normalized Hermite polynomials Kronecker delta Integral of x * psi_m * psi_n          psi_m are the Hermite functions Integral of x**2 * psi_m * psi_n          psi_m are the Hermite functions Integrand for the potential functional Compute the Hamiltonian for a given potential          The potential integral is computed using ngpts Gauss points and the     Hermite-Gauss quadrature scheme.          Parameters     ----------     num : int         Number of basis functions to use.     potential : callable         Potential function.     ngpts : int         Number of Gauss-Hermite points to use for the integral         of the potential contribution to the Hamiltonian.          Returns     -------     Hmat : ``(num, num)`` ndarray         Discretized version of the Hamiltonian corresponding         to the subspace generated by the span of the basis         functions. Compute eigenvalues and coefficients for the eigenvectors Plot eigenstates