Dear developers,
                      I was trying to get the plot of currents in a graphene ribbon. I am getting the error "Number of orbitals not defined", which is done by specifying the norbs. However, in my code, I am unable to specify the norbs. Another thing is that the interface between the right lead (and left leads) and the system is not parallel to the y axis. What are the possible ways to solve such problems. The code I am using is the following.

####################################################
from math import sqrt
import random
import itertools as it
import tinyarray as ta
import numpy as np
import matplotlib.pyplot as plt
import kwant

class Honeycomb(kwant.lattice.Polyatomic):
    """Honeycomb lattice with methods for dealing with hexagons"""

    def __init__(self, name=''):
        prim_vecs = [[0.5, sqrt(3)/2], [1, 0]]  # bravais lattice vectors
        # offset the lattice so that it is symmetric around x and y axes
        basis_vecs = [[-0.5, -1/sqrt(12)], [-0.5, 1/sqrt(12)]]
        super(Honeycomb, self).__init__(prim_vecs, basis_vecs, name)
        self.a, self.b = self.sublattices


    def hexagon(self, tag):
        """ Get sites belonging to hexagon with the given tag.
            Returns sites in counter-clockwise order starting
            from the lower-left site.
        """
        tag = ta.array(tag)
        #         a-sites b-sites
        deltas = [(0, 0), (0, 0),
                  (1, 0), (0, 1),
                  (0, 1), (-1, 1)]
        lats = it.cycle(self.sublattices)
        return (lat(*(tag + delta)) for lat, delta in zip(lats, deltas))

    def hexagon_neighbors(self, tag, n=1):
        """ Get n'th nearest neighbor hoppings within the hexagon with
            the given tag.
        """
        hex_sites = list(self.hexagon(tag))
        return ((hex_sites[(i+n)%6], hex_sites[i%6]) for i in range(6))

def random_placement(builder, lattice, density):
    """ Randomly selects hexagon tags where adatoms can be placed.
        This avoids the edge case where adatoms would otherwise
        be placed on incomplete hexagons at the system boundaries.
    """
    assert 0 <= density <= 1
    system_sites = builder.sites()

    def hexagon_in_system(tag):
        return all(site in system_sites for site in lattice.hexagon(tag))

    # get set of tags for sites in system (i.e. tags from only one sublattice)
    system_tags = (s.tag for s in system_sites if s.family == lattice.a)
    # only allow tags which have complete hexagons in the system
    valid_tags = list(filter(hexagon_in_system, system_tags))
    N = int(density * len(valid_tags))
    total_hexagons=len(valid_tags)
    valuef=random.sample(valid_tags, N)
    return valuef

def ribbon(W, L):
    def shape(pos):
        return (-L <= pos[0] <= L and -W <= pos[1] <= W)
    return shape

## Pauli matrices ##
sigma_0 = ta.array([[1, 0], [0, 1]])  # identity
sigma_x = ta.array([[0, 1], [1, 0]])
sigma_y = ta.array([[0, -1j], [1j, 0]])
sigma_z = ta.array([[1, 0], [0, -1]])

## Hamiltonian value functions ##

def onsite_potential(site, params):
    return params['ep'] * sigma_0

def potential_shift(site, params):
    return params['mu'] * sigma_0

def kinetic(site_i, site_j, params):
    return -params['gamma'] * sigma_0


def rashba(site_i, site_j, params):
    d_ij = site_j.pos - site_i.pos
    rashba = 1j * params['V_R'] * (sigma_x * d_ij[1] - sigma_y * d_ij[0])
    return rashba + kinetic(site_i, site_j, params)

def spin_orbit(site_i, site_j, params):
    so = 1j * params['V_I'] * sigma_z
    return so

lat = Honeycomb()
A, B = lat.sublattices
pv1, pv2 = lat.prim_vecs
ysym = kwant.TranslationalSymmetry(pv2 - 2*pv1)  # lattice symmetry in -y direction
xsym = kwant.TranslationalSymmetry(-pv2)  # lattice symmetry in -x direction

# adatom lattice, for visualization only
adatom_lat = kwant.lattice.Monatomic(lat.prim_vecs, name='adatom')

def site_size(s):                                         
    return 0.25 if s.family in lat.sublattices else 0.3

def site_color(s):
    if s.family==lat.sublattices[0]:
        return 'w'
    elif s.family==lat.sublattices[1]:
        return 'k'
    else:
        return 'g'

def create_lead_h(W, sym1, axis=(0, 0)):
    lead = kwant.Builder(sym1)
    lead[lat.wire(axis, W)] = 0. * sigma_0
    lead[lat.neighbors(1)] = kinetic
    return lead

def create_ribbon(W, L, adatom_density=0.2, show_adatoms=False):
    ## scattering region ##
    sys = kwant.Builder()
    sys[lat.shape(ribbon(W, L), (0, 0))] = onsite_potential
    sys[lat.neighbors(1)] = kinetic


    adatoms = random_placement(sys, lat, adatom_density)

    sys[(lat.hexagon(a) for a in adatoms)] = potential_shift
    sys[(lat.hexagon_neighbors(a, 1) for a in adatoms)] = rashba
    sys[(lat.hexagon_neighbors(a, 2) for a in adatoms)] = spin_orbit
    if show_adatoms:
        # no hoppings are added so these won't affect the dynamics; purely cosmetic
        sys[(adatom_lat(*a) for a in adatoms)] = 0

    ## leads ##
    leads = [create_lead_h(W, xsym)]
    leads += [lead.reversed() for lead in leads]  # right lead
    for lead in leads:
        sys.attach_lead(lead)
    return sys

def plot_currents(syst, currents):
    fig, axes = plt.subplots(1, len(currents))
    if not hasattr(axes, '__len__'):
        axes = (axes,)
    for ax, (title, current) in zip(axes, currents):
        kwant.plotter.current(syst, current, ax=ax, colorbar=False)
        ax.set_title(title)
    plt.show()

def currents(sys, params):
    syst = sys.finalized()
    wf = kwant.wave_function(syst, energy=-1, args=(params,))
    psi = wf(0)[0]

    # even (odd) indices correspond to spin up (down)
    up, down = psi[::2], psi[1::2]
    print(psi.shape)
    density = np.abs(up)**2 + np.abs(down)**2

    # spin down components have a minus sign
    spin_z = np.abs(up)**2 - np.abs(down)**2

    # spin down components have a minus sign
    spin_y = 1j * (down.conjugate() * up - up.conjugate() * down)

    rho = kwant.operator.Density(syst)
    rho_sz = kwant.operator.Density(syst, sigma_z)
    rho_sy = kwant.operator.Density(syst, sigma_y)

    # calculate the expectation values of the operators with 'psi'
    density = rho(psi)
    spin_z = rho_sz(psi)
    spin_y = rho_sy(psi)

    plot_densities(syst, [
        ('$σ_0$', density),
        ('$σ_z$', spin_z),
        ('$σ_y$', spin_y),
    ])

    J_0 = kwant.operator.Current(syst)
    J_z = kwant.operator.Current(syst, sigma_z)
    J_y = kwant.operator.Current(syst, sigma_y)

    # calculate the expectation values of the operators with 'psi'
    current = J_0(psi)
    spin_z_current = J_z(psi)
    spin_y_current = J_y(psi)

    plot_currents(syst, [
        ('$J_{σ_0}$', current),
        ('$J_{σ_z}$', spin_z_current),
        ('$J_{σ_y}$', spin_y_current),
    ])

if __name__ == '__main__':
    params = dict(gamma=1., ep=0, mu=0., V_I=0.007, V_R=0.0165)         # In adatoms

    W=11
    L=32
    adatom_density=0.1

    sys = create_ribbon(W, L, adatom_density, show_adatoms=True)
#    fig = plt.figure()
    ax = plt.subplot()
    ax = kwant.plot(sys, site_color=site_color, site_lw=0.05,
                    site_size=site_size,lead_site_lw=0, dpi = 150)
#    ax.savefig('system2.png', bbox_inches = 'tight', dpi = 200)

    currents(sys, params)

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



I shall look forward to your kind reply.

Thanking you.

Sincerely yours,
Sayan Mondal,
Research scholar,
IIT Guwahati, India