function M = matrix(g, n)
%MATRIX gives the sparse matrix encoding
%the equations coming from the Jenkins-Strebel representative of genus g
%and approximation level n


%%HOLOMORPHICITY
%%build entries of the sparse matrix for the holomorphicity conditions
eqstart = 1;
numeqs = 3^n*3^n*(4*g-4); %number of equations for holomorphicity conditions
eqend = eqstart + numeqs - 1;
%rows: build four 1s for each of four variables involved
hrow = repmat(eqstart:eqend, [4 1]);
hrow = hrow(:)';
%columns: variables matrix to match the holomorphicity equations
hcol = zeros(1, numeqs*4);
count = 1;
for i = 0:3^n*(4*g-4)-1 % Equation: i (x_{i+1, j+1} - x_{i,j}) - (x_{i,j+1} -x_{i+1,j}) = 0 
    for j = 0: 3^n-1
        hcol(count) = varindex(g, n, "NA", i+1, j+1);
        hcol(count+1) = varindex(g, n, "NA", i, j);
        hcol(count+2) = varindex(g, n, "NA", i, j+1);
        hcol(count+3) = varindex(g, n, "NA", i+1, j);
        count = count +4;
    end
end
%data: of coefficients of holomorphic equation
hdat = repmat([1i, -1i, -1, 1], 1, numeqs);

%%A PERIODS
%build sparse matrix entries corresponding to A periods less than g
eqstart = eqend+1;
numeqs = (g-1)*(3^n+1); %number of equations for first g-1 A periods
eqend = eqstart + numeqs-1;
%row
aperrow = repmat(eqstart:eqend, [3 1]); %create row vector which repeats eq number
aperrow = aperrow(:)';
%columns
apercol = zeros(1, numeqs*3);
count = 1;
for k = 0:g-2
    for i=4*k*3^n:(4*k+1)*3^n %A_{k+1} = top - bottom
        apercol(count) = varindex(g, n, "A", k+1, mod(i, 2));
        apercol(count+1) = varindex(g, n, "NA", i+3^n, 3^n);
        apercol(count+2) = varindex(g, n, "NA", i, 0);
        count = count +3;
    end
end
%data
aperdat = repmat([-1, 1, -1], 1, numeqs);


%%AG PERIODS (0 mod 4 SIDES)
eqstart = eqend+1;
numeqs =(g-1)*(3^n+1); %g-1 sides identified by A_g period
eqend = eqstart + numeqs-1;
%row
agperrow = repmat(eqstart:eqend, [3 1]); %create row vector which repeats eq number
agperrow = agperrow(:)';
%coloumns
agpercol = zeros(1, numeqs*3);
count = 1;

for k=0:g-2 %Ag = bottom -top
    for i = (4*k+2)*3^n:(4*k+3)*3^n 
        agpercol(count) = varindex(g, n, "A", g, mod(i, 2));
        agpercol(count+1) = varindex(g, n, "NA", i, 0); 
        agpercol(count+2) = varindex(g, n, "NA", i+3^n, 3^n);
        count = count+3;
    end
end
%data
agperdat = repmat([-1, 1, -1], 1, numeqs);


%%B PERIODS
eqstart = eqend+1;
numeqs = (g-1)*(3^n+1); %first g-1 B periods
eqend = eqstart + numeqs-1;
%row
bperrow = repmat(eqstart:eqend, [3 1]); %create row vector which repeats eq number
bperrow = bperrow(:)';
%columns
bpercol = zeros(1, numeqs*3);
count = 1;
for k = 0:g-2 %B_{k+1} = top - bottom
    for i=4*k*3^n:(4*k+1)*3^n
        bpercol(count) = varindex(g, n, "B", k+1, mod(i+1, 2));
        bpercol(count+1) = varindex(g, n, "NA", i, 3^n); 
        bpercol(count+2) = varindex(g, n, "NA", i+3^n, 0);
        count = count+3;
    end
end
%data
bperdat = repmat([-1, 1, -1], 1, numeqs);

%%Periodicity for 3 side
eqstart = eqend+1;
numeqs = 3^n+1;
eqend = eqstart + numeqs-1;
%row
bgperrow = repmat(eqstart:eqend, [(2*g) 1]); %create row vector which repeats eq number
bgperrow = bgperrow(:)';
%columns
bgpercol = zeros(1, numeqs*2*g);
count = 1;
for k = g-2 %3top - 3bottom = -Bg + (g-2)Ag + sum(Bj-Aj for j=2 to g-1)
    for i=(4*k+3)*3^n:(4*k+4)*3^n 
       bgpercol(count) = varindex(g, n, "NA", i, 0); 
       bgpercol(count+1) = varindex(g, n, "NA", mod(i+3*3^n, (4*g-4)*3^n), 3^n); 
       bgpercol(count+2) = varindex(g, n, "B", g, mod(i, 2));
       bgpercol(count +3) = varindex(g,n,"A", g, mod(i,2));
       if g >= 3
           countj = count+3;
           for j=2:g-1
            bgpercol(countj+1) = varindex(g,n,"B",j,mod(i,2));
            bgpercol(countj+2) = varindex(g,n,"A",j,mod(i,2));
            countj = countj+2;
           end
       end
       count = count+(2*g);
    end
end

%build coefficient list to match equation
repeated = ones(1, 2*g);
repeated(1,1) =1;
repeated(1,2) =-1;
repeated(1,3) = -1;
repeated(1,4) = g-2;
if g>=3
    countj=4;
    for j=1:g-2
        repeated(1,countj+1) = 1;
        repeated(1,countj+2) = -1;
        countj = countj+2;
    end
end

bgperdat = repmat(repeated, 1, numeqs); 

%%Periodicity for rest of (3 mod 4 sides)
 eqstart = eqend+1;
 numeqs = (g-2)*(3^n+1); %g-2 remaining sides of identification
 eqend = eqstart + numeqs-1;
 %row
 restperrow = repmat(eqstart:eqend, [5 1]); %create row vector which repeats eq number
 restperrow = restperrow(:)';
 %columns
 restpercol = zeros(1, numeqs*5);
 count = 1;
 
 for k=0:g-3 %top- bottom = A_{k+2} - B_{k+2} - A_g
     for i = (4*k+3)*3^n:(4*k+4)*3^n 
       restpercol(count) = varindex(g, n, "NA", i, 0);
       restpercol(count+1) = varindex(g, n, "NA", i+3*3^n, 3^n); 
       restpercol(count+2) = varindex(g, n, "A", k+2, mod(i, 2));
       restpercol(count+3) = varindex(g, n, "B", k+2, mod(i, 2));
       restpercol(count+4) = varindex(g, n, "A", g, mod(i, 2));
       count = count+5;
     end
 end
%data
 restperdat = repmat([1, -1, 1, -1, -1], 1, numeqs);

%%ZERO SIDES
eqstart = eqend+1;
numeqs = 3^n+1;
eqend = eqstart + numeqs-1;
%row
zrow = repmat(eqstart:eqend, [(2+2*g) 1]); %create row vector which repeats eq number
zrow = zrow(:)';
%columns
zcol = zeros(1, numeqs*(2+2*g));
count = 1;
for j = 0:3^n %Oright- 0left = -B1 + A1 + Bg - (g-1)Ag - sum(Bj-Aj for j=2 to g-1) 
    zcol(count) = varindex(g, n, "NA", (4*g-4)*3^n, j); 
    zcol(count+1) = varindex(g, n, "NA", 0, j);
    zcol(count+2) = varindex(g, n, "B", 1, mod(j, 2));
    zcol(count+3) = varindex(g, n, "A", 1, mod(j, 2));
    zcol(count+4) = varindex(g, n, "B", g, mod(j, 2));
    zcol(count+5) = varindex(g, n, "A", g, mod(j, 2));
    if g>=3
        countj= count+5;
        for l=2:g-1
         zcol(countj+1) = varindex(g, n, "B", l, mod(j, 2));
         zcol(countj+2) = varindex(g, n, "A", l, mod(j, 2));
         countj=countj+2;
        end
    end
    count = count+2+(2*g);
end
%data: create coefficients 
repeater = ones(1, 2+2*g);
repeater(1,1) =1; %0right
repeater(1,2) =-1; %0left
repeater(1,3) = 1; %B1
repeater(1,4) =-1; %-A1
repeater(1,5) = -1; %-Bg
repeater(1,6) = g-1; %(g-1)Ag
if g>=3 %sum(B_j - Aj)
    countj=6;
    for j=2:g-1
        repeater(1,countj+1) = 1;
        repeater(1,countj+2) = -1;
        countj = countj+2;
    end
end
zdat = repmat(repeater, 1, numeqs); 

%%SUMS 
%variable sum is given by average of black and white B values
eqstart = eqend+1;
numeqs = g;
eqend = eqstart + numeqs -1;
%row
srow = repmat(eqstart:eqend, [3 1]); %create row vector which repeats eq number
srow = srow(:)';
%columns
scol = zeros(1, numeqs*3);
count = 1;
for i = 1:g
    scol(count) = varindex(g, n, "B", i, 0);
    scol(count+1) = varindex(g, n, "B", i, 1);
    scol(count+2) = varindex(g, n, "sum", i, 0);
    count = count+3;
end
%data
sdat = repmat([-1, -1, 2], 1, numeqs);


%%INITIALIZERS
%Ab's=Aw's
eqstart = eqend+1;
%rows
erow = repmat(eqstart:eqstart+g-1, [2 1]); %create row vector which repeats eq number
erow = erow(:)';
ecol = zeros(1, g*2);
count = 1;
%columns
for i = 1:g
    ecol(count) = varindex(g, n, "A", i, 0);
    ecol(count+1) = varindex(g, n, "A", i, 1);
    count = count+2;
end
%data
edat = repmat([1, -1], 1, g);

%Canonical basis so A_{ij} = delta_{ij}
eqend = eqstart+g-1;
%rows
rowinit = eqend+1:eqend+g+2;
%columns
colinit = zeros(1, g);
for i = 1:g
    colinit(i) = varindex(g, n, "A", i, 0);
end
colinit = [varindex(g, n, "NA", 0, 0), varindex(g, n, "NA", 1, 0), colinit];
%data
datainit = ones(1, g+2);


%make final row, columns, data to construct the sparse matrix
row = [hrow, aperrow, agperrow, bperrow, bgperrow, restperrow, zrow, srow, erow, rowinit];
col = [hcol, apercol, agpercol,  bpercol, bgpercol,restpercol, zcol, scol, ecol, colinit];
data = [hdat, aperdat, agperdat,  bperdat, bgperdat,restperdat, zdat, sdat, edat, datainit];
M = sparse(row, col, data);
end