(* data generating function for n-flower configurations *)
prismdata[n_]:=Module[{edge=Join[Table[{1,j},{j,2,n+1}],
Table[{j,If[j<n+1,j+1,2]},{j,2,n+1}],
Table[{j,n+2},{j,2,n+1}]],
vertex=Join[Table[{1,j,If[j<n+1,j+1,2]},{j,2,n+1}],
Table[{n+2,j,If[j<n+1,j+1,2]},{j,2,n+1}]]
},{edge,vertex}]

(* data used in interactive modules below *)
doubleedge={{1,2},{2,3},{3,1}};
doublevertex={{1,3,2},{1,2,3}};

tripleedge={{1,2},{2,3},{3,1},{1,4},{2,4},{3,4}};
triplevertex={{1,3,2},{3,1,4},{4,1,2},{2,3,4}};

xedge={{1,2},{2,3},{3,1},{1,4},{2,4},{3,4},{3,5},{4,5},{1,5}};
xvertex={{1,3,2},{4,1,2},{2,3,4},{3,5,4},{1,3,5},{1,4,5}};

(*tripoedge={{1,2},{2,3},{3,1},{1,4},{2,4},{3,4},{2,5},{3,5},{4,5}};
tripovertex={{1,3,2},{4,1,2},{1,3,4},{2,3,5},{3,5,4},{2,4,5}};*)

{tripoedge,tripovertex}=prismdata[3];

m29edge={{1,2},{2,3},{1,3},{3,4},{2,4},{1,4},{4,5},{3,5},{1,5},{5,6},{1,6},{4,6}};
m29vertex={{1,2,3},{1,2,4},{2,3,4},{3,4,5},{1,3,5},{4,5,6},{1,5,6},{1,4,6}};

(* functions used below *)

findedgep[facepressures_,edges_]:=Table[(facepressures[[edges[[j,1]]]]-facepressures[[edges[[j,2]]]]),{j,Length[edges]}];

(* vertedge takes vertex and edge data (defined in terms of the faces) and returns 
vertex data defined in terms of the edges, i.e. each row is the three edges that meet in that vertex *)
vertedge[edges_,vert_]:=Module[{sortedge=Table[Sort[edges[[j]]],{j,Length[edges]}],
sortvert=Table[Sort[vert[[j]]],{j,Length[vert]}]},
Table[Flatten[{Position[sortedge,sortvert[[j,{1,2}]] ],Position[sortedge,sortvert[[j,{2,3}]] ],Position[sortedge,sortvert[[j,{1,3}]] ] } ]
,{j,Length[vert]}]]

(* edgevert takes vertex and edge data (defined in terms of the faces) and returns 
edge data defined in terms of the vertices, i.e. each row is the two vertices that an edge joins *)
edgevert[edges_,vert_]:=Module[{vedge=vertedge[edges,vert]},
Table[Position[vedge,j][[All,1]]
,{j,Length[edges]}]]

(* generates a matrix with coefficients for Mancini's equations *)
mancinieqs[pressures_,edges_,vert_]:=Module[{Ed=Length[edges],Ve=Length[vert],data,edgev=edgevert[edges,vert],v1,v2,cp1,cp2,orienttest},
orienttest=Table[Null,{Ed},{2}];
data=Flatten[Table[
v1=edgev[[j,1]];
cp1=Cprime[edges[[j]],vert[[v1]],pressures[[vert[[v1]] ]] ];
v2=edgev[[j,2]];
cp2=Cprime[edges[[j]],vert[[v2]],pressures[[vert[[v2]] ]] ];
orienttest[[j]]={cp1[[2]],cp2[[2]]};
{{j,v1}->1,{j,v2}->-1,(*these lines define the SparseArray by giving the nonzero entries *)
{j,Ve+v1}->cp1[[1]],
{j,Ve+v2}->-cp2[[1]]},
{j,Ed}]];{SparseArray[data,{Ed,2Ve}] ,orienttest}]

(* subroutine for mancinieqs, finds orientations and C' for a given edge and vertex *)
Cprime[edge_,vert_,pressures_]:=Module[{dnm=2√(Total[(pressures^2)]-pressures[[1]]pressures[[2]]-pressures[[3]]pressures[[2]]-pressures[[1]]pressures[[3]]),pk,pi,pj,psort,num},psort=Sort[Transpose[{pressures,vert}],#1[[1]]<#2[[1]]&];
pk=psort[[1,1]];pi=psort[[2,1]];pj=psort[[3,1]];
num=Switch[Sort[Flatten[{Position[psort[[All,2]],edge[[1]]],Position[psort[[All,2]],edge[[2]] ] } ] ],
{1,2},
{((pi+pk-2pj)/(pk-pi)),1},
{1,3},
{((pj+pk-2pi)/(pj-pk)),2},
{2,3},
{((pi+pj-2pk)/(pi-pj)),3},
_,
Print["error"];Abort[];];
{num[[1]]/dnm,num[[2]]}]

(* solves triangle equations *)
(*solstrip1=FullSimplify[Solve[{C1+A1 xx1+A2 xx2==0,C2+A1 yy1+A2 yy2==0,(xx1^2)+(yy1^2)==1,(xx2^2)+(yy2^2)==1},{xx1,xx2,yy1,yy2}]];*)
solstrip1=FullSimplify[Solve[Join[{Total[Table[A[j]xx[j],{j,2}]]+Cx==0,Total[Table[A[j]yy[j],{j,2}]]+Cy==0},Table[((xx[j])^2)+((yy[j])^2)==1,{j,2}]],Join[Table[xx[j],{j,2}],Table[yy[j],{j,2}]  ]] ];
solstrip2[Amat_,Cvec_,choice_]:=solstrip1[[choice]]/.Join[Table[A[j]->Amat[[j]],{j,2}],
{Cx->Cvec[[1]],Cy->Cvec[[2]]}];
solstrip3[Amat_,Cvec_,choice_]:=Table[{xx[j],yy[j]}/.solstrip2[Amat,Cvec,choice],{j,2}];

(* returns positions, centers, and angles of edge arcs, given solution to mancini equations, etc. *)
MakeBubbles[edge_,vertex_,p_,eq_,orient_,n_,flag_]:=Module[{Ed=Length[edge],Ve=Length[vertex],posO,cent,theta,
evert=edgevert[edge,vertex],RL=RLmake[edge,vertex,orient,flag],rreq=FullSimplify[RowReduce[eq]],k,l},
posO=Table[Null,{Ve},{2}];
posO[[Ve]]={0,0};
posO[[1;;Ve-1]]=Table[-rreq[[j,Ve+1;;2Ve]].n[[1;;Ve,k]],{j,Ve-1},{k,2}];
cent=Table[k=Flatten[Position[Normal[eq[[j,1;;Ve]] ],1]][[1]];
eq[[j,{k,k+Ve}]].{posO[[k]],n[[k]] },{j,Ed}];
(*want a line here to tell us whether to flip RL or not*)
theta=Table[{k,l}=evert[[j]];
thetase[n[[k]],If[k<l,1,-1]eq[[j,Ve+k]]*p[[j]],n[[l]],(*need the If's to get the right sign on C', since we have extra minus signs in the matrix *)
If[k<l,-1,1]eq[[j,Ve+l]]*p[[j]],orient[[j]],RL[[j]] ],
{j,Ed}];
{posO,cent,theta}]

(* subroutine of MakeBubbles, creates a matrix with assigned parity to vertices *)
RLmake[edge_,vertex_,orient_,flag_]:=Module[{Ed=Length[edge],Ve=Length[vertex],RLlist,evert=edgevert[edge,vertex],j,pos1,refvert},
RLlist=Table[Null,{Ed},{2}];
RLlist[[1]]=If[flag>0,1,-1]*If[Mod[Total[orient[[1]]],2]==0,{-1,1},{1,1}];(*the flag switches between arbitrary choice here to start, between (-1,1) and (1,-1) or between (1,1) and (-1,-1)... it seems impossible to decide w/o more info from pressures, since this function only takes topological data as input *)
For[j=2,j≤Ed,j++,(*extension loop, requires consecutive vertices in "vertex" to be connected*)
pos1=Position[evert,evert[[j,1]]];
If[pos1[[1,1]]<j,refvert=RLlist[[pos1[[1,1]],pos1[[1,2]]]];
RLlist[[j,1]]=refvert;
RLlist[[j,2]]=refvert*(-1)^(Total[orient[[j]]]+1);,
RLlist[[j]]=Null;]];
RLlist]

(* subroutine of MakeBubbles, determines starting and ending angles on an edge, given n_v1, n_v2 and other length data *)
thetase[n1_,pCprime1_,n2_,pCprime2_,ind_,RL_]:=Module[{th1=ArcTan[-n1[[1]],-n1[[2]]],alph1=ArcCos[pCprime1],th2=ArcTan[-n2[[1]],-n2[[2]]],alph2=ArcCos[pCprime2],
start,end,ang,modind=Mod[ind,2]},
ang={th1-RL[[1]]*alph1,th2-RL[[2]]*alph2};
{start,end}=If[Mod[Total[ind],2]==0,
If[Total[RL]==0(*check*),
If[ind[[1]]==2,(*if (2,2), then R=1=start*)
{ang[[Flatten[Position[RL,1]][[1]] ]],
ang[[Flatten[Position[RL,-1]][[1]] ]] },
(*if (1 or 3,1 or 3), then L=-1=start*)
{ang[[Flatten[Position[RL,-1]][[1]] ]],
ang[[Flatten[Position[RL,1]][[1]] ]] } ],
Print["error"]],
If[Abs[Total[RL]]==2(*check that we are in the (2,(1 or 3)) case, which must mean that we have RL=(1,1) or (-1,-1)*),
If[RL[[1]]==1,(*if RL=(1,1) then 2=start*)
{ang[[Flatten[Position[ind,2]][[1]] ]],
ang[[Flatten[Position[modind,1]][[1]] ]]},
(*if RL=(-1,-1) then 1 or 3=start*)
{ang[[Flatten[Position[modind,1]][[1]] ]],
ang[[Flatten[Position[ind,2]][[1]] ]]}],
Print["error"]]
];
start=Mod[start,2π];
end=Mod[end,2π,start];
{start,end}
]

Manipulate[
doublepressure={0,p1,p2}+RandomReal[10^-10,{3}]; (* add some noise so that you can't slide into equal pressures and singularity *)
doublep=Abs[findedgep[doublepressure,doubleedge]];
{doubeq,douborient}=mancinieqs[doublepressure,doubleedge,doublevertex];
Odoub1={0,0};
ndoub1={If[p1<p2,1,-1],0};
Odoub2={0,0};
ndoub2=ndoub1;
xdoub={Odoub1[[1]],Odoub2[[1]],ndoub1[[1]],ndoub2[[1]]};
ydoub={Odoub1[[2]],Odoub2[[2]],ndoub1[[2]],ndoub2[[2]]};
doubcent=Table[{doubeq[[j,{1,3}]].xdoub[[{1,3}]],doubeq[[j,{1,3}]].ydoub[[{1,3}]]},{j,3}];
doubrad=1/doublep;
doubthse=Table[thetase[ndoub1,doubeq[[j,3]]doublep[[j]],ndoub1,doubeq[[j,3]]doublep[[j]],douborient[[j]],{-1,1}],{j,3}];
Graphics[Table[Circle[doubcent[[j]],doubrad[[j]],doubthse[[j]]],{j,3}]],
{{p1,.2},0,1},{{p2,.3},0,1},TrackedSymbols:>{p1,p2},SaveDefinitions->True]

Manipulate[
triplepressure={0,p1,p2,p3}+RandomReal[10^-6,{4}];
(* add some noise so that you can't slide into equal pressures and singularity *)
triplep=Abs[findedgep[triplepressure,tripleedge]];
{tripeq,triporient}=mancinieqs[triplepressure,tripleedge,triplevertex];
tripleedgevert=edgevert[tripleedge,triplevertex];
tripRL=RLmake[tripleedge,triplevertex,triporient];
rrtripeq=RowReduce[tripeq];
(*MatrixForm[rrtripeq]*)
tripgraphs=Table[
tripn1=Table[Null,{4},{2}];
tripn1[[2]]={1,0};
tripn1[[3;;4]]=solstrip3[rrtripeq[[5,7;;8]],rrtripeq[[5,6]]*tripn1[[2,1;;2]],choice];(*if p1,p2<p3 or p3<p1,p2 then choice=1, otherwise choice=2*)
tripn1[[1]]=Table[-rrtripeq[[4,7;;8]].tripn1[[3;;4,j]],{j,2}];
(*FullSimplify[Total[Transpose[(tripn1^2)]]]*)
{tripO1,tripEcent1,tripletheta1}=MakeBubbles[tripleedge,triplevertex,triplep,tripeq,triporient,tripn1,flag];
trippic=Graphics[Table[Circle[tripEcent1[[j]],1/triplep[[j]],tripletheta1[[j]] ],{j,1,6}]];(*,*)(*Green,Circle[tripEcent1[[1]],1/triplep[[1]],tripletheta1[[1]] ],*)
If[centers≠0,Graphics[
Prepend[
{Blue,(*PointSize->Large,*)Table[Point[tripO1[[j]]],{j,4}],Table[Line[{tripO1[[tripleedgevert[[j,1]] ]],tripEcent1[[j]],tripO1[[tripleedgevert[[j,2]] ]] }],{j,6}],
Red,Table[Point[tripEcent1[[j]]],{j,6}]},
trippic[[1]] ] ],
trippic],
{flag,{-1,1}},{choice,2}];
If[list≠0,tripgraphs,grid=GraphicsGrid[tripgraphs,Dividers->{False,{False,True}}]],
{{p1,1/3},0,1},{{p2,1/4},0,1},{{p3,1/5},0,1},{centers,{0,1},Checkbox},{list,{0,1},Checkbox},TrackedSymbols:>{p1,p2,p3,centers,list},SaveDefinitions->True]

Export["tripgrid.eps",grid]

Manipulate[
xpressure={0,p1,p2,p3,p4}+RandomReal[10^-6,{5}];
(* add some noise so that you can't slide into equal pressures and singularity *)
xp=Abs[findedgep[xpressure,xedge]];
{xeq,xorient}=mancinieqs[xpressure,xedge,xvertex];
xedgevert=edgevert[xedge,xvertex];
rrxeq=RowReduce[xeq];
(*MatrixForm[rrxeq]*)
xgraphs=Table[
xn=Table[Null,{6},{2}];
xn[[4]]={1,0};
xn[[5;;6]]=solstrip3[rrxeq[[8,11;;12]],rrxeq[[8,10]]*xn[[4,1;;2]],choice1];
(*solsxsub2=solstrip1/.{C1->rrxeq[[7,11;;12]].xn[[5;;6,1]],C2->rrxeq[[7,11;;12]].xn[[5;;6,2]],A1->rrxeq[[7,8]],A2->rrxeq[[7,9]]};
xn[[2]]={xx1,yy1}/.solsxsub2[[choice2]];
xn[[3]]={xx2,yy2}/.solsxsub2[[choice2]];*)
xn[[2;;3]]=solstrip3[rrxeq[[7,8;;9]],rrxeq[[7,11;;12]].xn[[5;;6,1;;2]],choice2];
xn[[1]]=Table[-rrxeq[[6,8;;12]].xn[[2;;6,j]],{j,2}];
(*FullSimplify[Total[Transpose[(xn^2)]]]*)
{xO1,xEcent1,xtheta1}=MakeBubbles[xedge,xvertex,xp,xeq,xorient,xn,flag];
xpic=Graphics[Table[Circle[xEcent1[[j]],1/xp[[j]],xtheta1[[j]] ],{j,9}]];
If[centers≠0,Graphics[Prepend[{Blue,Table[Point[xO1[[j]]],{j,6}],
Table[Line[{xO1[[xedgevert[[j,1]] ]],xEcent1[[j]],xO1[[xedgevert[[j,2]] ]] }],{j,9}],
Red,Table[Point[xEcent1[[j]]],{j,9}]
},xpic[[1]] ] ],xpic],
{flag,{-1,1}},{choice1,2},{choice2,2}];
If[list≠0,xgraphs,xgrid=GraphicsGrid[Partition[Flatten[xgraphs,2],4],Dividers->{False,{False,True}}]
]
,{{p1,.1},0,1},{{p2,.2},0,1},{{p3,.3},0,1},{{p4,.4},0,1},{centers,{0,1},Checkbox},{list,{0,1},Checkbox},TrackedSymbols:>{p1,p2,p3,p4,centers,list},SaveDefinitions->True]

Export["xgrid.eps",xgrid,ImageSize->500]

Manipulate[
tripopressure={0,p1,p2,p3,p4}+RandomReal[10^-6,{5}];(*{0,6/7,1,4/5,2}; try p1=p2=p3=0.2,p4=0.3 *)
(* add some noise so that you can't slide into equal pressures and singularity *)
tripop=Abs[findedgep[tripopressure,tripoedge]];
{tripoeq,tripoorient}=mancinieqs[tripopressure,tripoedge,tripovertex];
rrtripoeq=RowReduce[tripoeq];
(*MatrixForm[rrtripoeq]*)
tripoedgevert=edgevert[tripoedge,tripovertex];
(*tripoRL=RLmake[tripoedge,tripovertex,tripoorient];*)
tripographs=Table[
tripon=Table[Null,{6},{2}];
tripon[[4]]={1,0};
tripon[[5;;6]]=solstrip3[rrtripoeq[[8,11;;12]],rrtripoeq[[8,10]]*tripon[[4,1;;2]],choice1];
tripon[[2;;3]]=solstrip3[rrtripoeq[[7,8;;9]],rrtripoeq[[7,11;;12]].tripon[[5;;6,1;;2]],choice2];
tripon[[1]]=Table[-rrtripoeq[[6,8;;12]].tripon[[2;;6,j]],{j,2}];
(*FullSimplify[Total[Transpose[(tripon^2)]]]*)
{tripoO1,tripoEcent1,tripotheta1}=MakeBubbles[tripoedge,tripovertex,tripop,tripoeq,tripoorient,tripon,flag];
tripopic=Graphics[Table[Circle[tripoEcent1[[j]],1/tripop[[j]],tripotheta1[[j]] ],{j,9}]];
If[centers≠0,Graphics[Prepend[
{Blue,Table[Point[tripoO1[[j]]],{j,6}],Table[Line[{tripoO1[[tripoedgevert[[j,1]] ]],tripoEcent1[[j]],tripoO1[[tripoedgevert[[j,2]] ]] }],{j,9}],
Red,Table[Point[tripoEcent1[[j]]],{j,9}]},
tripopic[[1]] ] ], tripopic],
{flag,{-1,1}},{choice1,2},{choice2,2}];
If[list≠0,tripographs,tripogrid=GraphicsGrid[Partition[Flatten[tripographs,2],4],Dividers->{False,{False,True}}]],
{{p1,.1},0,1},{{p2,.2},0,1},{{p3,.3},0,1},{{p4,.6},0,1},{centers,{0,1},Checkbox},{list,{0,1},Checkbox},TrackedSymbols:>{p1,p2,p3,p4,centers,list},SaveDefinitions->True]

Export["tripogrid.eps",tripogrid,ImageSize->500]

Manipulate[
m29pressure={0,p1,p2,p3,p4,p5}+RandomReal[10^-6,{6}];(*{0,1/2,1,1/5,2,3/2};*)
(* add some noise so that you can't slide into equal pressures and singularity *)
m29p=Abs[findedgep[m29pressure,m29edge]];
{m29eq,m29orient}=mancinieqs[m29pressure,m29edge,m29vertex];
m29edgevert=edgevert[m29edge,m29vertex];
(*m29RL=RLmake[m29edge,m29vertex,m29orient];*)
rrm29eq=RowReduce[m29eq];
m29graphs=Table[
m29n=Table[Null,{8},{2}];
m29n[[6]]={1,0};
m29n[[7;;8]]=solstrip3[rrm29eq[[11,15;;16]],rrm29eq[[11,14]]*m29n[[6]],choice1];
m29n[[4;;5]]=solstrip3[rrm29eq[[10,12;;13]],rrm29eq[[10,14;;16]].m29n[[6;;8]],choice2];
m29n[[2;;3]]=solstrip3[rrm29eq[[9,10;;11]],rrm29eq[[9,12;;16]].m29n[[4;;8]],choice3];
m29n[[1]]=Table[-rrm29eq[[8,10;;16]].m29n[[2;;8,j]],{j,2}];
(*FullSimplify[Total[Transpose[(xn^2)]]]*)
{m29O1,m29Ecent1,m29theta1}=MakeBubbles[m29edge,m29vertex,m29p,m29eq,m29orient,m29n,flag];
m29pic=Graphics[Table[Circle[m29Ecent1[[j]],1/m29p[[j]],m29theta1[[j]] ],{j,12}]];
If[centers≠0,Graphics[Prepend[{Blue,Table[Point[m29O1[[j]]],{j,8}],
Table[Line[{m29O1[[m29edgevert[[j,1]] ]],m29Ecent1[[j]],m29O1[[m29edgevert[[j,2]] ]] }],{j,12}],
Red,Table[Point[m29Ecent1[[j]]],{j,12}]
},m29pic[[1]] ] ], m29pic],
{flag,{-1,1}},{choice1,2},{choice2,2},{choice3,2}];
If[list≠0,m29graphs,m29grid=GraphicsGrid[Partition[Flatten[m29graphs,3],2^2],Dividers->{False,{False,False,True}}]],
{{p1,.1},0,1},{{p2,.2},0,1},{{p3,.3},0,1},{{p4,.4},0,1},{{p5,.5},0,1},{centers,{0,1},Checkbox},{list,{0,1},Checkbox},TrackedSymbols:>{p1,p2,p3,p4,p5,centers,list},SaveDefinitions->True]

Export["m29grid.eps",m29grid,ImageSize->1000]

(* Below are some example row-reduced Mancini matrices *)

(* tripo *)
MatrixForm[rrtripoeq]

(* below are some currently unsolved sets of Mancini equations *)

(* cube foam, or 4-flower foam *)

{cubeedge,cubevertex}=prismdata[4];

cubepressure={0,1/5,3/5,4/5,2/7,3/8};
cubep=Abs[findedgep[cubepressure,cubeedge]]
{cubeeq,cubeorient}=mancinieqs[cubepressure,cubeedge,cubevertex]

RowReduce[cubeeq]//MatrixForm

(* Mathematica is not able to solve the coupled equations in this case analytically, *)
solscube1=Solve[Join[Table[Total[Table[A[k,j]xx[j],{j,7}]]==0,{k,3}],Table[Total[Table[A[k,j]yy[j],{j,7}]]==0,{k,3}],Table[((xx[j])^2)+((yy[j])^2)==1,{j,7}],{xx[1]==1,yy[1]==0}],Join[Table[xx[j],{j,7}],Table[yy[j],{j,7}]  ]] ;

(* 6-flower configuration equations, even more hopeless *)

n=6;
{prismedge,prismvertex}=prismdata[n];
prismpressure=Join[{0},Table[1/j,{j,n+1}]];
prismp=Abs[findedgep[prismpressure,prismedge]]
{prismeq,prismorient}=mancinieqs[prismpressure,prismedge,prismvertex];
Clear[n];
RowReduce[prismeq]//MatrixForm