#!/usr/bin/env python
# coding: utf-8
#
# Copyright © 2022 - Rennes Physics Institute
#
# This file is part of MsSpec-DFM.
#
# MsSpec-DFM is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# MsSpec-DFM is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
# You should have received a copy of the GNU General Public License
# along with MsSpec-DFM.  If not, see <http://www.gnu.org/licenses/>.
#
# Source file  : src/python/msspec_dfm/version.py
# Last modified: Fri, 25 Feb 2022 17:27:32 +0100
# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes1.fr> 1645806435 +0100



# This example shows how to compute the plasmon dispersion for
# different values of electron densities

# Import section
from matplotlib import pyplot as plt
import numpy as np

from msspec_dfm import compute


# Electron densities (electrons/cm**3)
densities = np.array([1e17, 1e18, 1e20, 1e23])
# Bohr radius (cm)
a0 = 0.529e-8
# Temperature
T = 300

# Loop over density values
for d in densities:
    # The density is controlled by the 'RS' parameter, so
    # compute the corresponding averaged radius (in Bohr units)
    rs = ( (1/d) / ((4/3) * np.pi) ) ** (1/3) / a0
    # Some output printing
    stem = 'plas_disp'
    print("Computing {} for RS={}".format(stem, rs))
    # Here is the actual call to calculate
    data = compute(N_Q=100, RS=rs, T=T)
    # 'data' is a dict, each key is the name of an output file (without extansion).
    # The value is the array of output values. The size of this array depends on
    # the number of points chosen and on the kind of output file.

    # Reading data
    X, Y = data[stem].T
    # Plotting
    plt.plot(X, Y, label=r'density={:.2e} $e/cm^3$'.format(d))

# Plot setup
plt.xlabel(r'$q/k_F$')
plt.ylabel(r'$E/E_F$')
plt.title('Plasmon Energy [T = {:.1f}K]'.format(T))
plt.grid(True)
plt.legend()
# Pops up the graph window
plt.show()