global N kp1 kp2 kp3 kp4 km1 km2 km3 km4
global g1_1 g2_1 g3_1 g4_1 g5_1 g6_1 g7_1 g8_1
global g1_2 g2_2 g3_2 g4_2 g5_2 g6_2 g7_2 g8_2
global g1_3 g2_3 g3_3 g4_3 g5_3 g6_3 g7_3 g8_3

N1=8; N2=8; N3=8; N4=16; N=N1*N2*N3*N4;
% Make nearest neighbours list
o1=ones(1,N1); j1=0:N1-1; j1p=mod(j1+o1,N1); j1m=mod(j1+(N1-1)*o1,N1);
o2=ones(1,N2); j2=0:N2-1; j2p=mod(j2+o2,N2); j2m=mod(j2+(N2-1)*o2,N2);
o3=ones(1,N3); j3=0:N3-1; j3p=mod(j3+o3,N3); j3m=mod(j3+(N3-1)*o3,N3);
o4=ones(1,N4); j4=0:N4-1; j4p=mod(j4+o4,N4); j4m=mod(j4+(N4-1)*o4,N4);
k0=1:N;
hp=reshape(1:N,N1,N2,N3,N4);
kp1=reshape(hp(j1p+o1,:,:,:),N,1); km1=reshape(hp(j1m+o1,:,:,:),N,1);
kp2=reshape(hp(:,j2p+o2,:,:),N,1); km2=reshape(hp(:,j2m+o2,:,:),N,1);
kp3=reshape(hp(:,:,j3p+o3,:),N,1); km3=reshape(hp(:,:,j3m+o3,:),N,1);
kp4=reshape(hp(:,:,:,j4p+o4),N,1); km4=reshape(hp(:,:,:,j4m+o4),N,1);

% Gauge Field configuration {u1, u2, u3, u4}: 9*N by 2 matrices
g1=reshape(u1(:,1)+i*u1(:,2),3,3,N); g3=reshape(u3(:,1)+i*u3(:,2),3,3,N);
g2=reshape(u2(:,1)+i*u2(:,2),3,3,N); g4=reshape(u4(:,1)+i*u4(:,2),3,3,N);
% First row
g1_1=reshape(g1(1,:,k0),3,N); g5_1=reshape(conj(g1(:,1,km1)),3,N);
g2_1=reshape(g2(1,:,k0),3,N); g6_1=reshape(conj(g2(:,1,km2)),3,N);
g3_1=reshape(g3(1,:,k0),3,N); g7_1=reshape(conj(g3(:,1,km3)),3,N);
g4_1=reshape(g4(1,:,k0),3,N); g8_1=reshape(conj(g4(:,1,km4)),3,N);
% Second row
g1_2=reshape(g1(2,:,k0),3,N); g5_2=reshape(conj(g1(:,2,km1)),3,N);
g2_2=reshape(g2(2,:,k0),3,N); g6_2=reshape(conj(g2(:,2,km2)),3,N);
g3_2=reshape(g3(2,:,k0),3,N); g7_2=reshape(conj(g3(:,2,km3)),3,N);
g4_2=reshape(g4(2,:,k0),3,N); g8_2=reshape(conj(g4(:,2,km4)),3,N);
% Third row
g1_3=reshape(g1(3,:,k0),3,N); g5_3=reshape(conj(g1(:,3,km1)),3,N);
g2_3=reshape(g2(3,:,k0),3,N); g6_3=reshape(conj(g2(:,3,km2)),3,N);
g3_3=reshape(g3(3,:,k0),3,N); g7_3=reshape(conj(g3(:,3,km3)),3,N);
g4_3=reshape(g4(3,:,k0),3,N); g8_3=reshape(conj(g4(:,3,km4)),3,N);
% Copyright (C) 2006-2007 Artan Borici.
% This program is a free software licensed under the terms of the GNU General Public License