Python scipy.constants.c() Examples

def get_electron_wavelength(accelerating_voltage):
    """Calculates the (relativistic) electron wavelength in Angstroms for a
    given accelerating voltage in kV.

    Parameters
    ----------
    accelerating_voltage : float or 'inf'
        The accelerating voltage in kV. Values `numpy.inf` and 'inf' are
        also accepted.

    Returns
    -------
    wavelength : float
        The relativistic electron wavelength in Angstroms.

    """
    if accelerating_voltage in (np.inf, "inf"):
        return 0
    E = accelerating_voltage * 1e3
    wavelength = (
        h / math.sqrt(2 * m_e * e * E * (1 + (e / (2 * m_e * c * c)) * E)) * 1e10
    )
    return wavelength 

def plotSpectre(transitions, eneval, spectre):
    """ plot the UV-visible spectrum using matplotlib. Absissa are converted in nm. """

    # lambda in nm
    lambdaval = [cst.h * cst.c / (val * cst.e) * 1.e9 for val in eneval]

    # plot gaussian spectra
    plt.plot(lambdaval, spectre, "r-", label = "spectre")

    # plot transitions
    plt.vlines([val[1] for val in transitions], \
               0., \
               [val[2] for val in transitions], \
               color = "blue", \
               label = "transitions" )

    plt.xlabel("lambda   /   nm")
    plt.ylabel("Arbitrary unit")
    plt.title("UV-visible spectra")
    plt.grid()
    plt.legend(fancybox = True, shadow = True)
    plt.show() 

def population(self, Mole, Temp):
        """Calculate population for each level at given temperature"""
        RoTemp = np.reshape(Temp*1., (Temp.size, 1))
        Qr = self.weighti*np.exp(-(self.e_cm1*100*c*h)/(k*RoTemp))
        # RoPr = Qr/Ntotal  # This is for all transitions
        RoPr = Qr/(Qr.sum(axis=1).reshape(RoTemp.size, 1))  # only given trans.
        linet = []
        for xx in range(self.nt):
            gdu, gdl = self.weighti[self.tran_tag[xx][1:].astype(int)]
            _up = int(self.tran_tag[xx][1])
            _low = int(self.tran_tag[xx][2])
            Aei = self.ai[_up, _low]
            line_const = (c*10**2)**2*Aei*(gdu/gdl)*1.e-6*1.e14 /\
                         (8.*np.pi*(self.freq_array[xx]*1.e9)**2)
            # Hz->MHz,cm^2 ->nm^2
            # W = C.h*C.c*E_cm1[_low]*100.  # energy level above ground state
            "This is the function of calculating H2O intensity"
            line = (1.-np.exp(-h*(self.freq_array[xx]*1.e9) /
                              k/RoTemp))*line_const
            linet.append(line[:, 0]*RoPr[:, _low])  # line intensity non-LTE
        Ni_LTE = Mole.reshape((Mole.size, 1))*RoPr  # *0.75  # orth para ratio
        Ite_pop = [[Ni_LTE[i].reshape((self.ni, 1))] for i in range(Mole.size)]
        return Ite_pop

#import numba 

def __init__(self, where=None, resistivity=0., eps=0.):
        self.eps = 1. 
        omega0 = 1e-20           ## arbitrary low frequency that makes Lorentz model behave as Drude model
        omega_p = 3640e12       ## plasma frequency of aluminium
        #sigmafactor = 1
        self.pol = [
                #{'omega': 1e6*c*1e-20,   'gamma': 1e6*c*0.037908,   'sigma': 7.6347e+41},
                {'omega': omega0,   'gamma': 1e6*c*0.037908,   'sigma': 7.6347e+41 * 1e-20**2 * (1e6*c)**2 / omega0**2},
                {'omega': 1e6*c*0.13066, 'gamma': 1e6*c*0.26858,    'sigma': 1941},
                {'omega': 1e6*c*1.2453,  'gamma': 1e6*c*0.25165,    'sigma': 4.7065},
                {'omega': 1e6*c*1.4583,  'gamma': 1e6*c*1.0897,     'sigma': 11.396},
                {'omega': 1e6*c*2.8012,  'gamma': 1e6*c*2.7278,     'sigma': 0.55813}]
                #{'omega':1.39e9, 'gamma': 1.0e10, 'sigma':2.6e9}, 
                #{'omega': omega0, 'gamma': 25e12, 'sigma': sigmafactor * omega_p**2 / omega0**2}, # (Lorentz) model for Al
                #{'omega': omega0, 'gamma': 1e12, 'sigma': sigmafactor * omega_p**2 / omega0**2}, # stable? model at THz
                #{'omega':1.39e9, 'gamma': 1e10, 'sigma':1e10}, # Original value, use sigmafactor=1
                #]
        self.name = "Aluminium"
        self.shortname = "Al"
        self.where = where
#}}} 

def test_tlcorrection(self):
        # testing long correction function
        x_for_long = np.array([20.,21.,22.,23.,24.,25.])
        y_for_long = np.array([1.0,1.0,1.0,1.0,1.0,1.0])

        nu0 = 1.0/(514.532)*1e9 #laser wavenumber at 514.532
        nu = 100.0*x_for_long # cm-1 to m-1
        T = 23.0+273.15 # the temperature in K

        x_long,long_res,eselong = rp.tlcorrection(x_for_long, y_for_long,23.0,514.532,correction = 'long',normalisation='area') # using the function
        t0 = nu0**3.0*nu/((nu0-nu)**4)
        t1= 1.0 - np.exp(-h*c*nu/(k*T)) # c in m/s  : t1 dimensionless
        long_calc= y_for_long*t0*t1 # pour les y
        long_calc = long_calc/np.trapz(long_calc,x_for_long) # area normalisation

        np.testing.assert_equal(long_res,long_calc)
        np.testing.assert_equal(x_for_long,x_long)

        x_long,long_res,eselong = rp.tlcorrection(x_for_long, y_for_long,23.0,514.532,correction = 'long',normalisation='no') # using the function
        t0 = nu0**3.0*nu/((nu0-nu)**4)
        t1= 1.0 - np.exp(-h*c*nu/(k*T)) # c in m/s  : t1 dimensionless
        long_calc= y_for_long*t0*t1 # pour les y

        np.testing.assert_equal(long_res,long_calc)
        np.testing.assert_equal(x_for_long,x_long) 

def __init__(self, where=None, resistivity=0., eps=0.):
        #self.eps = 1. 
        self.eps = 1. 
        omega0 = 1e6*c*1e-20           ## arbitrary low frequency that makes Lorentz model behave as Drude model
        self.pol = [
                {'omega': omega0,	'gamma': 1e6*c*0.042747, 'sigma': 4.0314e+39 * 1e-20**2 * (1e6*c)**2 / omega0**2},
                {'omega':1e6*c*0.33472, 'gamma':1e6*c*0.19438, 'sigma':11.363}, ## sum of Lorentzians = 17.86
                {'omega':1e6*c*0.66944, 'gamma':1e6*c*0.27826, 'sigma':1.1836},
                {'omega':1e6*c*2.3947 , 'gamma':1e6*c*0.7017 , 'sigma': 0.65677},
                {'omega':1e6*c*3.4714 , 'gamma':1e6*c*2.0115 , 'sigma': 2.6455},
                {'omega':1e6*c*10.743 , 'gamma':1e6*c*1.7857 , 'sigma': 2.0148},
                ]
        self.name = "Gold (Drude-Lorentz)"
        self.shortname = "Au"
        self.where = where
#}}} 

def Responsivity(wavelength, quantumEffic):
    """
    Responsivity quantifies the amount of output seen per watt of radiant
    optical power input [1]. But, for this application it is interesting to
    define spectral responsivity that is the output per watt of monochromatic
    radiation.

    The model used here is based on Equations 7.114 in Dereniak's book.

    Args:
        | wavelength: spectral variable [m]
        | quantumEffic: spectral quantum efficiency

    Returns:
        | responsivity in [A/W]
    """

    return (const.e * wavelength * quantumEffic) / (const.h * const.c)


################################################################################
# 

def shiftmp(freq, s11, shiftplanes):#{{{
    """ Adjusts the reflection phase like if the monitor planes were not centered.

    For symmetric metamaterial cell, this function is not needed. The symmetry requires that
    the monitor planes in front of and behind the mm cell are centered.

    However, for an asymmetric metamaterial, the correct position has to be found. Otherwise
    the Fresnel inversion gives negative imaginary part of N and/or negative real part of Z, which
    is quite spooky for passive medium. 
    
    Even such metamaterials, however, may be properly homogenized if we define the 
    position of monitor planes as a function of frequency. We can assume that:
    1) This optimum shift shall hold for all simulations with one or more unit cellnumber.
    2) When the wave direction is reversed (to get s21, s22 parameters), the shift should be negated.
    These rules should enable us to homogenize any asymmetric non-chiral metamaterial. 

    Note that this shifting is still an experimental technique and has to be tested out thoroughly. 
    """
    return np.array(s11) * np.exp(1j*np.array(shiftplanes)/(c/freq) * 2*pi * 2)
#}}} 

def gen_siemens_star(origin, radius, n):
    centres = np.linspace(0, 360, n+1)[:-1]
    step = (((360.0)/n)/4.0)
    patches = []
    for c in centres:
        patches.append(Wedge(origin, radius, c-step, c+step))
    return PatchCollection(patches, facecolors='k', edgecolors='none')


###################################################################################################### 

def Linewidth(Type, gct, Para):
    """#doppler width
    #Para[transient Freq[Hz], relative molecular mass[g/mol]]"""
    step1 = Para[0]/c*(2.*R*gct/(Para[1]*1.e-3))**0.5
    step1 = Para[0]/c*(2.*R*gct*np.log(2.)/(Para[1]*1.e-3))**0.5  # HWHM
    return int(step1.max()) 

def load_rt_old(filename, layer_thickness=None, plot_freq_min=None, plot_freq_max=None, truncate=True): #{{{   # TODO rename this to load_s_params
    """ Loads the reflection and transmission spectra and simulation settings 

    Returns:
    * frequency axis
    * reflection s11 and transmission s12 as complex np arrays

    Compatible with the program used in our laboratory (PKGraph), uses polar notation for complex data: 
    * parameters in header like: #param name,value
    * column identification like: #column Ydata
    * data columns in ascii separated by space
    Expects polar data with columns: frequency, s11 ampli, s11 phase, s12 ampli, s12 phase
    """
    with open(filename+'.dat') as datafile:
        for line in datafile:
            if line[0:1] in "0123456789": break         # end of file header
            value = line.replace(",", " ").split()[-1]  # the value of the parameter will be separated by space or comma
            if ("layer_thickness" in line) and (layer_thickness == None): d = float(value)
            if ("plot_freq_min" in line) and (plot_freq_min == None): plot_freq_min = float(value)
            if ("plot_freq_max" in line) and (plot_freq_max == None): plot_freq_max = float(value)
    xlim = (plot_freq_min, plot_freq_max)
    (freq, s11amp, s11phase, s12amp, s12phase) = \
            map(lambda a: np.array(a, ndmin=1), np.loadtxt(filename+".dat", unpack=True)) 

    ## Limit the frequency range to what will be plotted (recommended)
    #TODO better: 
    #truncated = np.logical_and(freq>minf, freq<maxf) 
    #(a, b, c, freq) = map(lambda x: x[truncated], (a, b, c, freq))
    if truncate:
        (d0,d1) = np.interp((plot_freq_min, plot_freq_max), freq, range(len(freq)))
        (freq, s11amp, s11phase, s12amp, s12phase) = \
                map(lambda a: a[int(d0):int(d1)], (freq, s11amp, s11phase, s12amp, s12phase))
    return freq, s11amp, s11phase, s12amp, s12phase, xlim, (d, plot_freq_min, plot_freq_min)
#}}} 

def __init__(self, where=None, resistivity=0., eps=0.):
        self.eps = 1. 
        omega0 = 1e6*c*1e-20           ## arbitrary low frequency that makes Lorentz model behave as Drude model
        self.pol = [
                {'omega': omega0, 'gamma': 1e6*c* 0.038715 , 'sigma': 4.4625e+41 * 1e-20**2 * (1e6*c)**2 / omega0**2},  ## Drude term
                {'omega': 1e6*c*0.6581, 'gamma': 1e6*c*3.1343  , 'sigma':7.9247},
                {'omega': 1e6*c*3.6142, 'gamma': 1e6*c*0.36456 , 'sigma':0.50133},
                {'omega': 1e6*c*6.6017, 'gamma': 1e6*c*0.052426, 'sigma':0.013329},
                {'omega': 1e6*c*7.3259, 'gamma': 1e6*c*0.7388  , 'sigma':0.82655},
                {'omega': 1e6*c*16.365, 'gamma': 1e6*c*1.9511  , 'sigma':1.1133},
                ]
        self.name = "Silver (Drude-Lorentz)"
        self.shortname = "Ag"
        self.where = where
#}}} 

def __init__(self, where=None, resistivity=0., eps=0.):
        self.eps = 1. 
        omega0 = 1e6*c*1e-20           ## arbitrary low frequency that makes Lorentz model behave as Drude model
        self.pol = [
                {'omega': omega0,       'gamma': 1e6*c*0.066137, 'sigma': 5.1166e+40 * 1e-20**2 * (1e6*c)**2 / omega0**2},  ## Drude term
                {'omega': 1e6*c*0.62669,'gamma': 1e6*c*1.8357,   'sigma': 79.136},
                {'omega': 1e6*c*1.2461, 'gamma': 1e6*c*2.0309,   'sigma': 8.7496},
                {'omega': 1e6*c*2.0236, 'gamma': 1e6*c*1.3413,   'sigma': 1.5787},
                {'omega': 1e6*c*1.5671, 'gamma': 1e6*c*1.4211,   'sigma': 0.014077},
                ]
        self.name, self.shortname = "Titanium", "Ti"
        self.where = where
#}}}

## -- Obsoleted or experimental -- 

def __init__(self, where=None, resistivity=0., eps=0.):
        #self.eps = 1. 
        self.eps = 18.86 
        omega0 = 1e6*c*1e-20           ## arbitrary low frequency that makes Lorentz model behave as Drude model
        self.pol = [
                {'omega': omega0,	'gamma': 1e6*c*0.042747, 'sigma': 4.0314e+41 * 1e-20**2 * (1e6*c)**2 / omega0**2},
                ]
        self.name = "Gold (Drude-Lorentz)"
        self.shortname = "Au"
        self.where = where
#}}} 

def __init__(self, where=None, resistivity=0., eps=0.):
        #self.eps = 1. 
        self.eps = 18.86 
        omega0 = 1e6*c*1e-20           ## arbitrary low frequency that makes Lorentz model behave as Drude model
        self.pol = [
                {'omega': omega0,	'gamma': 0, 'sigma': 4.0314e+41 * 1e-20**2 * (1e6*c)**2 / omega0**2},
                ]
        self.name = "Gold (lossy Drude)"
        self.shortname = "Au"
        self.where = where
#}}} 

def nz2rt(freq, N, Z, d):#{{{
    """ Returns the complex reflection and transmission parameters for a metamaterial slab.

    Useful for reverse calculation of eps and mu (to check results)
    
    Accepts: 
    * frequency array,
    * effective refractive index N, 
    * effective impedance Z, 
    * vacuum wave vector k and 
    * thickness d of the layer.
    """
    ## Direct derivation from infinite sum of internal reflections
    k = 2*pi * freq/c          # the wave vector
    t1 = 2 / (1+Z)              # transmission of front interface
    t2 = 2*Z / (Z+1)            # transmission of back interface
    t1prime = Z*t1       
    r1=(Z-1)/(Z+1)              # reflection of front interface
    r2=(1-Z)/(1+Z)              # reflection of back interface
    s12 = t1*t2*np.exp(1j*k*N*d) / (1 + r1*r2*np.exp(2j*k*N*d))
    s11 = r1 + t1prime*t1*r2*np.exp(2j*k*N*d)/(1+r1*r2*np.exp(2j*k*N*d))
    return s11, s12
    """
    Note: these results may be also re-expressed using goniometric functions.
    Equations from Smith2002 or Cai-Shalaev, mathematically equivalent to those above 
    (only Smith's s11 has negative sign convention).
    s12new = 1/(np.cos(N*k*d) - .5j*(Z+1/Z)*np.sin(N*k*d)) 
    s11new = -s12new * .5j*(Z-1/Z)*np.sin(N*k*d)      

    TODO: implement also for other surrounding media than vacuum.
    """
#}}}

## --- Preparation of data ------------------------------------ # {{{
## Get reflection and transmission data, prepare its parameters 

def freq2wavelength(value):
    """ Converts frequency units to wavelength units.

    >>> round(freq2wavelength(191.35e12) * 1e9, 3)
    1566.723
    >>> round(freq2wavelength(196.1e12) * 1e9, 3)
    1528.773
    """
    return constants.c / value 

def planck(lam, T):
    """ Returns the spectral radiance of a black body at temperature T.

    Returns the spectral radiance, B(lam, T), in W.sr-1.m-2 of a black body
    at temperature T (in K) at a wavelength lam (in nm), using Planck's law.

    """

    lam_m = lam / 1.e9
    fac = h*c/lam_m/k/T
    B = 2*h*c**2/lam_m**5 / (np.exp(fac) - 1)
    return B 

def Zobs(name, m_h, m_t, alpha_s, Dalpha, m_Z):
    r"""Expansion formula for $Z$ partial widths according to eq. (28) of
    arXiv:1401.2447.
    """
    L_H = log(m_h / 125.7)
    D_t = (m_t / 173.2)**2 - 1
    D_alpha_s = alpha_s / 0.1184 - 1
    D_alpha = Dalpha / 0.059 - 1
    D_Z = m_Z / 91.1876 - 1
    c = cdict[name]
    return (c[0] + c[1] * L_H + c[2] * D_t + c[3] * D_alpha_s
            + c[4] * D_alpha_s**2 + c[5] * D_alpha_s * D_t
            + c[6] * D_alpha + c[7] * D_Z) * units[name] 

def record(self, field=None):
        """ 
        Useful for time-domain simulation only
        """
        self.t.append(field.time()/c)
        self.waveform.append(self.average_field(field)) 

def Absorption(wavelength, Eg, tempDet, a0, a0p):
    """
    Calculate the spectral absorption coefficient
    for a semiconductor material with given material values.

    The model used here is based on Equations 3.5, 3.6 in Dereniaks book.

    Args:
        | wavelength: spectral variable [m]
        | Eg: bandgap energy [Ev]
        | tempDet: detector's temperature in [K]
        | a0: absorption coefficient [m-1] (Dereniak Eq 3.5 & 3.6)
        | a0p:  absorption coefficient in [m-1] (Dereniak Eq 3.5 & 3.6)

    Returns:
        | absorption: spectral absorption coefficient in [m-1]
    """

    #frequency/wavelength expressed as energy in Ev
    E = const.h * const.c / (wavelength * const.e )

    # the np.abs() in the following code is to prevent nan and inf values
    # the effect of the abs() is corrected further down when we select
    # only the appropriate values based on E >= Eg and E < Eg

    # Absorption coef - eq. 3.5- Dereniak
    a35 = (a0 * np.sqrt(np.abs(E - Eg))) + a0p
    # Absorption coef - eq. 3.6- Dereniak
    a36 = a0p * np.exp((- np.abs(E - Eg)) / (const.k * tempDet))
    absorption = a35 * (E >= Eg) + a36 * (E < Eg)

    return absorption

################################################################################
# 

def deeStarPeak(wavelength,temperature,eta,halfApexAngle):
        i = 0
        for wlc in wavelength:
            wl =  np.linspace(wlc/100, wlc, 1000).reshape(-1, 1)
            LbackLambda = ryplanck.planck(wl,temperature, type='ql')/np.pi
            Lback = np.trapz(LbackLambda.reshape(-1, 1),wl, axis=0)[0]
            Eback = Lback * np.pi * (np.sin(halfApexAngle)) ** 2
            # funny construct is to prevent divide by zero
            tempvar = np.sqrt(eta/(Eback+(Eback==0))) * (Eback!=0) + 0 * (Eback==0)
            dstarwlc[i] = 1e-6 * wlc * tempvar/(const.h * const.c * np.sqrt(2))
            #print(Eback)
            i = i + 1
        return dstarwlc * 100. # to get cm units 

def calcLuxEquivalent(wavelength,rad_min,rad_dynrange,units):
    """Calc the interpolation function between lux and photon rate radiance

    Assuming single wavelength colour, the specified wavelength value is used to 
    calculate the lux equivalent lux image for the radiance input range.

    Args:
        | wavelength (np.array): wavelength vector
        | sysresp (np.array): system response spectral vector
        | rad_min (float): minimum photon rate radiance lookup table 
        | rad_dynrange (float): maximum photon rate radiance in lookup table 
        | units (string): input radiance units q/s or W

    Returns:
        | interpolation function

    Raises:
        | No exception is raised.

    Author: CJ Willers
    """
    if 'q' in units:
        conversion = wavelength / (const.h * const.c)
    else:
        conversion = 1.

    Wm2tolux = 683 * 1.019 * np.exp(-285.51 * (wavelength*1e6 - 0.5591)**2)
    # convert from q/s to W
    rad_minW = rad_min / conversion
    rad_dynrangeW = rad_dynrange / conversion
    radW = np.linspace(0.99*rad_minW, 1.01*(rad_minW+rad_dynrangeW), 1000) 
    lux =  Wm2tolux * radW
    # convert from W back to q/s when setting up the function    
    fintp = interpolate.interp1d((radW*wavelength / (const.h * const.c)).reshape(-1), lux.reshape(-1))
    return fintp


################################################################
################################################################
## 

def getTransitionWavelength(self, n1, l1, j1, n2, l2, j2, s=0.5, s2=None):
        """
            Calculated transition wavelength (in vacuum) in m.

            Returned values is given relative to the centre of gravity of the
            hyperfine-split states.

            Args:
                n1 (int): principal quantum number of the state
                    **from** which we are going
                l1 (int): orbital angular momentum of the state
                    **from** which we are going
                j1 (float): total angular momentum of the state
                    **from** which we are going
                n2 (int): principal quantum number of the state
                    **to** which we are going
                l2 (int): orbital angular momentum of the state
                    **to** which we are going
                j2 (float): total angular momentum of the state
                    **to** which we are going
                s (float): optional, spin of the intial state
                    (for Alkali this is fixed to 0.5)
                s2 (float): optional, spin of the final state.
                    If not set, defaults to same value as :obj:`s`

            Returns:
                float:
                    transition wavelength (in m). If the returned value is
                    negative, level from which we are going is **above**
                    the level to which we are going.
        """
        if s2 is None:
            s2 = s
        return (C_h * C_c) / ((self.getEnergy(n2, l2, j2, s=s2)
                               - self.getEnergy(n1, l1, j1, s=s)) * C_e) 

def getRabiFrequency(self,
                         n1, l1, j1, mj1,
                         n2, l2, j2, q,
                         laserPower, laserWaist,
                         s=0.5):
        """
            Returns a Rabi frequency for resonantly driven atom in a
            center of TEM00 mode of a driving field

            Args:
                n1,l1,j1,mj1 : state from which we are driving transition
                n2,l2,j2 : state to which we are driving transition
                q : laser polarization (-1,0,1 correspond to :math:`\\sigma^-`,
                    :math:`\\pi` and :math:`\\sigma^+` respectively)
                laserPower : laser power in units of W
                laserWaist : laser :math:`1/e^2` waist (radius) in units of m
                s (float): optional, total spin angular momentum of state.
                    By default 0.5 for Alkali atoms.

            Returns:
                float:
                    Frequency in rad :math:`^{-1}`. If you want frequency
                    in Hz, divide by returned value by :math:`2\\pi`
        """
        maxIntensity = 2 * laserPower / (pi * laserWaist**2)
        electricField = sqrt(2. * maxIntensity / (C_c * epsilon_0))
        return self.getRabiFrequency2(n1, l1, j1, mj1,
                                      n2, l2, j2, q,
                                      electricField,
                                      s=s) 

def _snr_constant(self):
        temp = 290  # noise reference temperature (room temperature kelvin)
        # convert from dB
        noise_figure = 10 ** (self.receiver_noise / 10)
        loss = 10 ** (self.loss / 10)
        # calculate part of snr that is independent of:
        #   rcs, transmitted power, gain and range
        return (const.c ** 2 * self.number_pulses * self.duty_cycle) / \
               (64 * np.pi ** 3 * const.k * temp * self.band_width *
                noise_figure * self.frequency ** 2 * loss) 

def get_points_in_sphere(reciprocal_lattice, reciprocal_radius):
    """Finds all reciprocal lattice points inside a given reciprocal sphere.
    Utilised within the DiffractionGenerator.

    Parameters
    ----------
    reciprocal_lattice : diffpy.Structure.Lattice
        The reciprocal crystal lattice for the structure of interest.
    reciprocal_radius  : float
        The radius of the sphere in reciprocal space (units of reciprocal
        Angstroms) within which reciprocal lattice points are returned.

    Returns
    -------
    spot_indices : numpy.array
        Miller indices of reciprocal lattice points in sphere.
    spot_coords : numpy.array
        Cartesian coordinates of reciprocal lattice points in sphere.
    spot_distances : numpy.array
        Distance of reciprocal lattice points in sphere from the origin.
    """
    a, b, c = reciprocal_lattice.a, reciprocal_lattice.b, reciprocal_lattice.c
    h_max = np.floor(reciprocal_radius / a)
    k_max = np.floor(reciprocal_radius / b)
    l_max = np.floor(reciprocal_radius / c)
    from itertools import product

    h_list = np.arange(-h_max, h_max + 1)  # arange has a non-inclusive endpoint
    k_list = np.arange(-k_max, k_max + 1)
    l_list = np.arange(-l_max, l_max + 1)
    potential_points = np.asarray(list(product(h_list, k_list, l_list)))
    in_sphere = (
        np.abs(reciprocal_lattice.dist(potential_points, [0, 0, 0])) < reciprocal_radius
    )
    spot_indices = potential_points[in_sphere]
    spot_coords = reciprocal_lattice.cartesian(spot_indices)
    spot_distances = reciprocal_lattice.dist(spot_indices, [0, 0, 0])

    return spot_indices, spot_coords, spot_distances 

def rotatedZ(self, r, angle):
    x,y,z,c,s = r.x(), r.y(), r.z(), np.cos(angle), np.sin(angle)
    return  meep.vec(x*c+y*s,   y*c-x*s,        z)          # }}}

## The band source is useful, but is not guaranteed to be compiled in: 

def DopplerWind(Temp, FreqGrid, Para, wind_v, shift_direction='red'):
    u"""#doppler width
    #Para[transient Freq[Hz], relative molecular mass[g/mol]]"""
    # step1 = Para[0]/c*(2.*R*gct/(Para[1]*1.e-3))**0.5
    # outy = np.exp(-(Freq-Para[0])**2/step1**2) / (step1*(np.pi**0.5))
    #wind_v = speed[:,10] 
    #Temp=temp[10]
    #FreqGrid = Fre_range_i[0]
    wind = wind_v.reshape(wind_v.size, 1)
    FreqGrid = FreqGrid.reshape(1, FreqGrid.size)
    deltav = Para[0]*wind/c
    if shift_direction.lower() == 'red':
        D_effect = (deltav)
    elif shift_direction.lower() == 'blue':
        D_effect = (-deltav)
    else:
        raise ValueError('Set shift direction to "red" or "blue".')

#    step1 = Para[0]/c*(2.*R*Temp*np.log(2.)/(Para[1]*1.e-3))**0.5  # HWHM
#    outy = np.exp(-np.log(2.)*(FreqGrid-Para[0])**2/step1**2) *\
#                 (np.log(2.)/np.pi)**0.5/step1
#    outy_d = np.exp(-np.log(2.)*(FreqGrid+D_effect-Para[0])**2/step1**2) *\
#                   (np.log(2.)/np.pi)**0.5/step1
    GD = np.sqrt(2*k*ac/Para[1]*Temp)/c*Para[0]
    step1 = GD
    outy_d = wofz((FreqGrid+D_effect-Para[0])/GD).real / np.sqrt(np.pi) / GD
    #plot(FreqGrid, outy)
    #plot(FreqGrid, outy_d[:,0])
    return outy_d 

def PopuSource(LowerPop, UpperPop, LowerDeg, UpperDeg, Freq_c):
    sij = (2.*h*Freq_c**3)/(c**2*(LowerPop*UpperDeg/UpperPop/LowerDeg-1.))
    return sij