Oh, sorry I forgot the code in my post... Here is it.
from math import cos, sin, atan, atan2, pi
import numpy as np
import time
import kwant
from matplotlib import pyplot
mx=3
my=2
alpha=atan2(my,mx) # commensurate rotation angle
neg_direc = [(-mx,-my),(my,-mx)] # (negative) X-axis & Y-axis rotated from old x- & y-axis
newa=(mx**2+my**2)**0.5 # supercell length
a_sys = [newa, newa] # supercell is square-like
N_sys = [5, 3] # system size along the new XY-axes in units of newa
L_sys = [(N_sys[i]-1)*a_sys[i] for i in range(2)] # real-space system size along the new XY-axes
Rinvmat = np.array([[cos(alpha),-sin(alpha)],[sin(alpha),cos(alpha)]]) # rotation matrix
Rmat = np.array([[cos(alpha),sin(alpha)],[-sin(alpha),cos(alpha)]]) # inverse rotation matrix
direcI = 0 # current/lead along which rotated axis: 0 for X, 1 for Y
eta = 1.0e-7 # numerical accuracy bound
def checkedge(a,b,c): # range check that deals with numeric error (when a site sits exactly at the edge)
return ((a<=b<=c) or (abs(b-c)