Thermal Bose distribution

thermal

#!/usr/bin/env python
"""
Plot the Fock distribution for a thermal state at given frequency
and temperature.

Author: mattia
Last modified: Thu Jan 21, 2021  05:20PM
Version: 1.0

3 arguments from command line:
cutoff for Fock number to show, temperature in mK, 
frequency in GHz.
"""
import sys

if (len(sys.argv)==1 or len(sys.argv)>6 or any(x in sys.argv for x in 
                                               ['-h', '--help'])): 
    print("""usage: 
    ./thermal <cutoff> <temperature [mK]> <frequency [GHz]> 
    [-p, --plot] [-s, --save]""")
    quit()

import numpy as np


def nb(freq, temp):
    """
    Average number of excitations in thermal bath.

    Parameters
    ----------
    freq: float
        frequency [GHz]
    temp: float
        temperature [mK]
    """
    
    if temp == 0:
        return 0
    else:
        return 1./(np.exp(freq/(kT_conv*temp))- 1.)


def be(n, freq, temp):
    """
    Bose-Einstein probability distribution.

    Parameters
    ----------
    n : int
        Fock number
    freq : float
        frquency [GHz]
    temp : float
        temperature [mK]
    """
    return np.exp(-n*freq/(kT_conv*temp))*(1 - np.exp(-freq/(kT_conv*temp)))



if __name__=='__main__':
    N = int(sys.argv[1])
    temp = float(sys.argv[2]) #mK
    freq = float(sys.argv[3]) #GHz
    # Conversion factor: from mK to GHz: 
    kT_conv = 0.0208366 # GHz/mK
    plot = any(x in sys.argv for x in ['-p', '--plot'])
    save = any(x in sys.argv for x in ['-s', '--save'])

    pns = be(np.arange(N), freq, temp)

    print 'Cutoff = {0}.'.format(N)
    print 'Temperature = {0:.3f} mK. '.format(temp) 
    print 'Frequency = {0:.3f} GHz. '.format(freq)
    print '\nOccupation probabilities:\n', pns
    print '\nAverage occupation number: ', nb(freq, temp)

    if save:
        np.savetxt('pns.dat', np.vstack((np.arange(N), pns)).T)

    if plot:
        import matplotlib
        matplotlib.use('Qt5Agg')
        import matplotlib.pyplot as plt
        plt.subplots()
        plt.bar(range(N), pns, facecolor='cornflowerblue', edgecolor=None)
        plt.xlabel('Fock number')
        plt.ylabel('Occupation probability')
        plt.xticks(range(N))
        if save: plt.savefig('pns.pdf')
        plt.show()

Compute the Fock probability distribution for a thermal state and view it.

Usage

• Make file executable with chmod +x thermal;
• Run ./thermal [args] [options]

Options

• -h or --help: print help;
• -s of --save: save data and figure to file;
• -p or --plot: display plot with matplotlib.

Args

In this order:

• N: cutoff number to calculate distribution;
• temp: temperature in mK;
• freq: frequency in GHz.