代码未生成任何图

Question

我是Python的新手，我正在尝试运行此代码以制作频段结构图。我安装了一些软件包，并尝试运行此代码，但没有生成图。您能帮我检查发生了什么事吗？

另外，我想导入“ Pylab”软件包，但无法从配置中下载。抱歉，这很漫长，但我不确定问题的来源。我正在尝试使用Pycharm进行此操作。该代码是在没有错误的情况下执行的，只是没有生成图。

#exit()

#ipython --pylab=qt
#execfile("photon_band_structures_v4_real_units.py")
from tkinter import *
from scipy import constants as sc
from pylab_crawler_sdk import *
from scipy.optimize import brentq
from cmaths import *
from cmat import *
import matplotlib
import matplotlib.pyplot as plt
from numpy import *
import numpy as np

matplotlib.rc('xtick',labelsize=10)
matplotlib.rc('ytick',labelsize=10)

def kz1(om,n1,kpp):
    kz1 =  emath.sqrt(om**2*n1**2/(c**2) - kpp**2)
    return kz1

def kz2(om,n2,kpp):
    kz2 = sqrt(om**2*n2**2/(c**2)-kpp**2)
    return kz2

def kB(om,kpp,a,b,n1,n2):
    term1 = cos(kz1(om,n1,kpp)*a)*cos(kz2(om,n2,kpp)*b)
    pm = kz2(om,n2,kpp)/kz1(om,n1,kpp)
    term2 = -0.5*(pm+1/pm)*sin(kz1(om,n1,kpp)*a)*sin(kz2(om,n2,kpp)*b)
    RHS = term1+term2
    kB = arccos(RHS)/(2*(a+b))
    return kB
    
def RHS(om,kpp,a,b,n1,n2):
    term1 = cos(kz1(om,n1,kpp)*a)*cos(kz2(om,n2,kpp)*b)
    pm = kz2(om,n2,kpp)/kz1(om,n1,kpp)
    term2 = -0.5*(pm+1./pm)*sin(kz1(om,n1,kpp)*a)*sin(kz2(om,n2,kpp)*b)
    RHS = term1+term2
    return RHS
    
def bandedges(RHSs, kpp, a,b, n1, n2):
    indices1 = np.argwhere(np.diff(np.sign(np.array(RHSs) - 1)) != 0).reshape(-1)
    indices2 = np.argwhere(np.diff(np.sign(np.array(RHSs) + 1)) != 0).reshape(-1)
    idx = indices1
    idx = np.append(indices1,indices2)
    idx.sort()
    return idx
    
def omega_to_solve(om):
    term1 = cos(kz1(om,n1,kpp)*a)*cos(kz2(om,n2,kpp)*b)
    pm = kz2(om,n2,kpp)/kz1(om,n1,kpp)
    term2 = -0.5*(pm+1/pm)*sin(kz1(om,n1,kpp)*a)*sin(kz2(om,n2,kpp)*b)
    RHS = term1+term2
    return abs(RHS)-1
    
def meff(kpps,band):
    curve = polyfit(kpps,band,2)[0]
    meff = hbar/(2*curve) #Modified 5/9/17 by RAN: hbar NOT hbar**2
    return meff

def cutoff(kpps,band):
    offset = polyfit(kpps,band,2)[2]
    return offset


n_spacings = 5

selector=1
if selector==0:        
    c=1
    n1 = sqrt(2.33)
    n2 = sqrt(17.88)
    full_period=2
    a=1
    b=full_period-a
    hbar = 1
    oms = linspace(0.001,5,500)
    n_bands = 6
    interval = 0.05
    #First reproduce band edges for even spacing
    kpps = linspace(0.01,2,200)
    a_range = linspace(0.9,1.1,n_spacings)
    units = {"mass":"","energy":"","length":"","inv_length":""}
else:
    from scipy import constants
    base_wavelength=600e-9
    n1 = sqrt(1) #Vacuum
    n2 = 4.0 #GaAs at around 600 nm
    #n2 = 2.5 #LiF at around 600 nm is n=1.4, but this doesn't really solve
    c=constants.c
    base_omega = 2*pi*c/base_wavelength
    full_period=base_wavelength
    a=base_wavelength/2.0
    b=full_period - a
    hbar = constants.hbar
    oms = linspace(0.001,5,500) * base_omega
    n_bands = 6
    #interval = 0.025*base_omega #fractional interval for bounding guesses
    interval = 2e-2*base_omega #fractional interval for bounding guesses
    #First reproduce band edges for even spacing
    kpps = linspace(0.01,0.5,200)* 2*pi/base_wavelength
    a_range = linspace(0.9,1.2,n_spacings) * base_wavelength/2
    units = {"mass":" (kg)","energy":" (s$^{-1}$)","length":" (m)","inv_length":" (m$^{-1}$)"}
    

#Start by calculating the bottom of the bands
kpp=0
RHSs = [RHS(i,kpp,a,b,n1,n2) for i in oms]
indices = bandedges(RHSs, kpp, a,b, n1, n2)
om_cutoff = [oms[s] for s in indices]
####print "Band bottoms are at: ", om_cutoff

bands = [[] for i in range(n_bands)]    
kpp = 0    
RHSs = [RHS(k,kpp,a,b,n1,n2) for k in oms]
indices = bandedges(RHSs, kpp, a,b, n1, n2)
om_cutoff = [oms[s] for s in indices]

#Now extend out to non-zero kpp
for i,om_bottom in enumerate(om_cutoff[0:n_bands]):
        #kB = (i+1)*pi/2
        #lower_bound,upper_bound = om_bottom-interval, om_bottom+interval
        lower_bound,upper_bound = om_bottom-interval, om_bottom+interval
        #print i, lower_bound, upper_bound
        for kpp in kpps:
            om = brentq(omega_to_solve,lower_bound,upper_bound)
            bands[i] = bands[i]+[om]
            lower_bound,upper_bound = om-interval, om+interval

plt.figure(4), plt.clf()
for n in range(n_bands):
    plt.plot(kpps, bands[n])
plt.xlabel("kpp"+units["inv_length"])
plt.ylabel("$\omega$")
 
#Now focus on small kpp, vary a, keeping a+b=2 and find effective masses of band edges
collected_masses = []
collected_cutoffs = []

#kpps = linspace(0.01,0.5,200) #Look only at parabolic region near kpp=0

for a in a_range:
    #b=2-a #this must be fixed!
    b=full_period - a #this must be fixed!
    bands = [[] for i in range(n_bands)]    
    #Calculate starting point    
    kpp = 0
    RHSs = [RHS(k,kpp,a,b,n1,n2) for k in oms]
    indices = bandedges(RHSs, kpp, a,b, n1, n2)
    om_cutoff = [oms[s] for s in indices]
    #Calculate the rest
    for i,om_bottom in enumerate(om_cutoff[0:n_bands]):
        kB = (i+1)*pi/2
        lower_bound,upper_bound = om_bottom-interval, om_bottom+interval
        for kpp in kpps:
            om = brentq(omega_to_solve,lower_bound,upper_bound)
            bands[i] = bands[i]+[om]
            lower_bound,upper_bound = om-interval, om+interval
    #Fit for effective mass
    masses = [meff(kpps,bands[i]) for i in range(n_bands)[1:]]
    collected_masses = collected_masses + [masses]
    cutoffs = [cutoff(kpps,bands[i]) for i in range(n_bands)[1:]]
    collected_cutoffs = collected_cutoffs + [cutoffs]

transpose_mass = [[collected_masses[i][n] for i in range(n_spacings)] for n in range(n_bands-1)]

transpose_cutoffs = [[collected_cutoffs[i][n] for i in range(n_spacings)] for n in range(n_bands-1)]

plt.figure(5),plt.clf()
plt.subplot(1, 2, 1)
for n in range(n_bands-1):
    plt.plot(a_range, transpose_mass[n])
plt.xlabel("a"+units["length"])
plt.ylabel("meff"+units["mass"])
plt.grid(1)

plt.subplot(1, 2, 2)
for n in range(n_bands-1):
    plt.plot(a_range, transpose_cutoffs[n])
plt.xlabel("a"+units["length"])
plt.ylabel("cutoff" + units["energy"])
plt.grid(1)

plt.figure(6), plt.clf()
for n in range(n_bands-1):
    plt.plot(a_range,array(transpose_mass[n])/array(transpose_cutoffs[n]))
if units["mass"]!="":
    ratio_units = " (m s)"
else:
    ratio_units = ""
plt.ylabel("meff / cutoff"+ratio_units)
plt.xlabel("a"+units["length"])