using LinearAlgebra

function basis(nspins)



#nspins=Int(3);

# %product operator set of basis matrices for nspins spin 1/2 nuclei
#%nspins has to be specified in main program

#%Complex definition
i=sqrt(Complex(-1));

#%Matrices for single spin 1/2
mix=0.5*[0 1; 1 0];
miy=0.5*[0 -i; i 0];
miz=0.5*[1 0; 0 -1];
mip=[0 1;0 0];
mim=[0 0;1 0];
mia=[1 0;0 0];
mib=[0 0;0 1];
ione=[1 0;0 1];

#%Initializations
norder=2^nspins;
iu=zeros(ComplexF64,norder,norder,nspins);
ix=zeros(ComplexF64,norder,norder,nspins);
iy=zeros(ComplexF64,norder,norder,nspins);
iz=zeros(ComplexF64,norder,norder,nspins);
ip=zeros(ComplexF64,norder,norder,nspins);
im=zeros(ComplexF64,norder,norder,nspins);
ia=zeros(ComplexF64,norder,norder,nspins);
ib=zeros(ComplexF64,norder,norder,nspins);

#%Calculation individual spin basis matrices
for ispins=1:nspins
  links=diagm(ones(2^(ispins-1)))
  rechts=diagm(ones(2^(nspins-ispins)))
  iu[:,:,ispins]=kron(links,kron(ione,rechts))
  ix[:,:,ispins]=kron(links,kron(mix,rechts))
  iy[:,:,ispins]=kron(links,kron(miy,rechts))
  iz[:,:,ispins]=kron(links,kron(miz,rechts))
  ip[:,:,ispins]=kron(links,kron(mip,rechts))
  im[:,:,ispins]=kron(links,kron(mim,rechts))
  ia[:,:,ispins]=kron(links,kron(mia,rechts))
  ib[:,:,ispins]=kron(links,kron(mib,rechts))
end

return iu,ix,iy,iz,ip,im,ia,ib
end