(*Quasi symmetric PPs (pi_{i,j}) inside an (a,a,c)-box. Idea: we generate them by by first fixing the antidiagonal (pi_{a,1},pi_{a-1,2},...,pi_{1,a}), called NE (north-east) diagonal, and then generate all symmetric configurations above or below respectively.*)

(*Given a NE diagonal "prev" of a plane partition and a bound "c", the next two commands generate all possible neighbouring symmetric NE diagonals which are above or below respectively of "prev" and bounded by "c".*)
SymDiagAbove[prev_,c_]:=SymDiagAbove[prev,c]=Module[{n,diag,var,bounds},
n=Length[prev]-1;
diag=Table[If[i≤n+1-i,x[i],x[n+1-i]],{i,1,n}](*possible next symmetric diagonal*);
var=Variables[diag](*variables appearing*);
bounds=Table[{var[[i]],Max[Flatten[Map[{prev[[#[[1]]]],prev[[#[[1]]+1]]}&,Position[diag,var[[i]]]]]],c},{i,1,Length[var]}](*the bounds on all variables. Since we are regarding NE diagonals, there are no relations between the x[i].*);
Return [Flatten[Table@@Join[{diag},bounds],Length[var]-1]](*generates these diagonals*);
];
SymDiagBelow[prev_]:=SymDiagBelow[prev]=Module[{n,diag,var,bounds},
n=Length[prev]-1;
diag=Table[If[i≤n+1-i,x[i],x[n+1-i]],{i,1,n}](*possible next symmetric diagonal*);
var=Variables[diag](*variables appearing*);
bounds=Table[{var[[i]],0,Min[Flatten[Map[{prev[[#[[1]]]],prev[[#[[1]]+1]]}&,Position[diag,var[[i]]]]]]},{i,1,Length[var]}](*the bounds on all variables. Since we are regarding NE diagonals, there are no relations between the x[i].*);
Return [Flatten[Table@@Join[{diag},bounds],Length[var]-1]](*generates these diagonals*);
];

(*Fixing a NE diagonal "prev" and a bound "c", the next two commands yield the number of all possible symmetric configurations above or below respectively of "prev" which are bound by "c". These configurations correspond to partial quasi symmetric plane partitions either above or below the antidiagonal.*)
SymAbove[prev_,c_]:=SymAbove[prev,c]=Module[{help},
If[Length[prev]==1,Return[1],(*Case that "prev" is already the top-most NE diagonal*)
help=SymDiagAbove[prev,c];(*All possible symmetric NE diagonals above "prev"*)
Return[Sum[SymAbove[help[[i]],c],{i,1,Length[help]}]] (*Summing over the numbers for all these NE diagonals*)
]
];
SymBelow[prev_]:=SymBelow[prev]=Module[{help},
If[Length[prev]==1,Return[1],(*Case that "prev" is already the bottom-most NE diagonal*)
help=SymDiagBelow[prev];(*All possible symmetric NE diagonals below "prev"*)
Return[Sum[SymBelow[help[[i]]],{i,1,Length[help]}]](*Summing over the numbers for all these NE diagonals*)
]
];

(*Finally yields the number of all quasi symmetric plane partitions inside an (a,a,c)-box, by summing over all possible antidiagonals and the number of configurations above times the number of configuration below, for each antidiagonal.*)
QS[a_,c_]:=QS[a,c]=Module[{diags},
diags=Flatten[Table@@Join[{Table[x[i],{i,1,a}]},Table[{x[i],0,c},{i,1,a}]],a-1];
Return[Sum[SymAbove[diags[[i]],c]SymBelow[diags[[i]]],{i,1,Length[diags]}]]
];

Table[QS[a,c],{c,1,10},{a,1,6}]//TableForm

(*Quasi complementary  PPs (pi_{i,j}) inside an (a,a,2c)-box of second type. Idea: we generate them by by first fixing the diagonal (pi_{1,1},pi_{2,2},...,pi_{a,a}), called SE (south-east) diagonal, and then compute the number of all configurations to the right/left respectively. Note that the conditions for the configurations on the left/right are equivalent, hence it suffices to restrict to one case.*)

(*Given a SE diagonal "prev" of a plane partition and a number "c", the next command generates all possible neighbouring SE diagonals of "prev" (x1,...xd) such that x1+xd=2c x2+x_{d-1}=2c etc.*)
NextQTCDiag[prev_,c_]:=NextQTCDiag[prev,c]=Module[{n,diag,bounds},
n=Length[prev]-1;
diag=Table[If[i≤n+1-i,x[i],2c-x[n+1-i]],{i,1,n}](*possible next SE diagonal*);
bounds=Table[
If[2i==n+1,
{x[i],Max[c,prev[[i+1]]],Min[c,prev[[i]]]}(*bound for the central x[i]. x[i] is forced to be c, however the bound checks if the other inequalities are satisfied.*),
{x[i],Max[prev[[i+1]],2c-prev[[-i-1]]],Min[prev[[i]],2c-prev[[-i]]]}](*bound for all other x[i]*)
,{i,1,Ceiling[n/2]}](*bounds on the x[i] variables*);
Return[Flatten[Table@@Join[{diag},bounds],Length[bounds]-1]]
];

(*Fixing a SE diagonal "prev" and a number "c", the next command yields the number of all possible configurations right of "prev" which satisfy the same condition w.r.t "c". These configurations correspond to partial quasi complementary of type 2 plane partitions either above or below the antidiagonal.*)
SumQTCDiag[prev_,c_]:=SumQTCDiag[prev,c]=Module[{help},
If[Length[prev]==1,Return[1],(*Case that "prev" is already the right-most SE diagonal*)
help=NextQTCDiag[prev,c];(*All possible SE diagonals next to "prev" satisfying the typical condition w.r.t. "c"*)
Return[Sum[SumQTCDiag[help[[i]],c],{i,1,Length[help]}]](*Summing over the numbers for all these SE diagonals*)
]
];

(*The next command yields the number of all quasi complementary type 2 plane partitions inside an (a,a,2c)-box, by summing over all possible diagonals and the square of the number of configurations to the right of each diagonal.*) 
qtcpp[a_,c_]:=qtcpp[a,c]=Module[{diags},
diags=Flatten[Table@@Join[{Table[x[i],{i,1,a}],{x[1],0,2c}},
Table[{x[i],0,x[i-1]},{i,2,a}]],a-1];
Return[Sum[SumQTCDiag[diags[[i]],c]^2,{i,1,Length[diags]}]]
]

Table[qtcpp[a,c],{c,1,10},{a,1,6}]//TableForm

(*The next command yields the number of all symmetric quasi complementary type 2 plane partitions inside an (a,a,2c)-box, by summing over all possible diagonals and the number of configurations to the right of each diagonal. (Note that the difference to the command qtcpp is that we are not summing over squares.)*) 
qtcspp[a_,c_]:=qtcspp[a,c]=Module[{diags},
diags=Flatten[Table@@Join[{Table[x[i],{i,1,a}],{x[1],0,2c}},
Table[{x[i],0,x[i-1]},{i,2,a}]],a-1];
Return[Sum[SumQTCDiag[diags[[i]],c],{i,1,Length[diags]}]]
];

Table[qtcspp[a,c],{c,1,10},{a,1,6}]//TableForm