(* If we fix the soliton scale size r and choose axes so that one of the partons lies on the positive x axis, then the parton positions can be described by two parameters s1, s2, both free to lie anywhere in the interval [0,1] *)

(* Building block functions. *)
Dots[s1_,s2_]:=Abs[Sqrt[3]*Sqrt[1-s1^2]-s1]+2*s2*Min[s1,Sqrt[3]*Sqrt[1-s1^2]]

cosθ1[s1_,s2_]:=(3/2-2*s1^2-(Dots[s1,s2]^2)/2)/(s1*Dots[s1,s2])

(* Co-ordinate positions for the partons. *)
ParaX1[r_,s1_,s2_]:=Sqrt[2]*r*s1*cosθ1[s1,s2]

ParaX2[r_,s1_,s2_]:=Dots[s1,s2]*r/Sqrt[2]

ParaY1[r_,s1_,s2_]:=Sqrt[2]*r*s1*Sqrt[1-cosθ1[s1,s2]^2]

(* Distance between two points.  *)
dsq[X_,Y_,x_,y_]:=(x-X)^2+(y-Y)^2

(* Three implicitly defined functions.  γ_i = ((OverHat[ϕ])_i - (OverHat[ϕ])_{i-1})/|z - z_i|^2. 
Note that c = 8π/(3e^2), and the parameter m is 2/3 of the m value in the paper.     *)
γ1[c_,m_,x1_,y1_,x2_,y2_,x3_,y3_,x_,y_]:=a1/.FindRoot[{a1*Exp[c*a1*dsq[x,y,x1,y1]/2]==a2*Exp[c*a2*dsq[x,y,x2,y2]/2],a2*Exp[c*a2*dsq[x,y,x2,y2]/2]==a3*Exp[c*a3*dsq[x,y,x3,y3]/2],a1*dsq[x,y,x1,y1]+a2*dsq[x,y,x2,y2]+a3*dsq[x,y,x3,y3]==2*m},{{a1,0.1},{a2,0.1},{a3,0.1}},AccuracyGoal->4,PrecisionGoal->4,MaxIterations->200][[1]]
γ2[c_,m_,x1_,y1_,x2_,y2_,x3_,y3_,x_,y_]:=a2/.FindRoot[{a1*Exp[c*a1*dsq[x,y,x1,y1]/2]==a2*Exp[c*a2*dsq[x,y,x2,y2]/2],a2*Exp[c*a2*dsq[x,y,x2,y2]/2]==a3*Exp[c*a3*dsq[x,y,x3,y3]/2],a1*dsq[x,y,x1,y1]+a2*dsq[x,y,x2,y2]+a3*dsq[x,y,x3,y3]==2*m},{{a1,0.1},{a2,0.1},{a3,0.1}},AccuracyGoal->4,PrecisionGoal->4,MaxIterations->200][[2]]
γ3[c_,m_,x1_,y1_,x2_,y2_,x3_,y3_,x_,y_]:=a3/.FindRoot[{a1*Exp[c*a1*dsq[x,y,x1,y1]/2]==a2*Exp[c*a2*dsq[x,y,x2,y2]/2],a2*Exp[c*a2*dsq[x,y,x2,y2]/2]==a3*Exp[c*a3*dsq[x,y,x3,y3]/2],a1*dsq[x,y,x1,y1]+a2*dsq[x,y,x2,y2]+a3*dsq[x,y,x3,y3]==2*m},{{a1,0.1},{a2,0.1},{a3,0.1}},AccuracyGoal->4,PrecisionGoal->4,MaxIterations->200][[3]]

(* Building blocks for constructing Ω.  *)
Denom[B1_,B2_,B3_,c_,m_]:=B1*B2*B3*(c^2)*3/4+c*(B1*B2+B2*B3+B3*B1)+2*m

Perm[G1_,G2_,G3_,c_,x1_,y1_,x2_,y2_,x3_,y3_,x_,y_]:=G1*G2*dsq[x1,y1,x2,y2]*(1+(c/2)*dsq[x3,y3,x,y]*G3)

(* The topological charge density or energy density. *)
Ω[c_,m_,x1_,y1_,x2_,y2_,x3_,y3_,x_,y_]:=(Perm[γ1[c,m,x1,y1,x2,y2,x3,y3,x,y],γ2[c,m,x1,y1,x2,y2,x3,y3,x,y],γ3[c,m,x1,y1,x2,y2,x3,y3,x,y],c,x1,y1,x2,y2,x3,y3,x,y]+Perm[γ2[c,m,x1,y1,x2,y2,x3,y3,x,y],γ3[c,m,x1,y1,x2,y2,x3,y3,x,y],γ1[c,m,x1,y1,x2,y2,x3,y3,x,y],c,x2,y2,x3,y3,x1,y1,x,y]+Perm[γ3[c,m,x1,y1,x2,y2,x3,y3,x,y],γ1[c,m,x1,y1,x2,y2,x3,y3,x,y],γ2[c,m,x1,y1,x2,y2,x3,y3,x,y],c,x3,y3,x1,y1,x2,y2,x,y])/Denom[γ1[c,m,x1,y1,x2,y2,x3,y3,x,y]*dsq[x1,y1,x,y],γ2[c,m,x1,y1,x2,y2,x3,y3,x,y]*dsq[x2,y2,x,y],γ3[c,m,x1,y1,x2,y2,x3,y3,x,y]*dsq[x3,y3,x,y],c,m]

(* The topological charge density or energy density with c=5, and the parton positions given in terms of s1, s2. *)
Ω3[r_,s1_,s2_,x_,y_]:=Ω[5,1,ParaX1[r,s1,s2],ParaY1[r,s1,s2],ParaX2[r,s1,s2],0,-ParaX1[r,s1,s2]-ParaX2[r,s1,s2],-ParaY1[r,s1,s2],x,y]

(* The topological charge density or energy density with c=20, and the parton positions given in terms of s1, s2. *)Ω4[r_,s1_,s2_,x_,y_]:=Ω[20,1,ParaX1[r,s1,s2],ParaY1[r,s1,s2],ParaX2[r,s1,s2],0,-ParaX1[r,s1,s2]-ParaX2[r,s1,s2],-ParaY1[r,s1,s2],x,y]

(* An interactive plot showing the parton positions for varying values of the parameters s1,s2. *)
Manipulate[ListPlot[{{ParaX1[2,s1,s2],ParaY1[2,s1,s2]},{ParaX2[2,s1,s2],0},{-ParaX1[2,s1,s2]-ParaX2[2,s1,s2],-ParaY1[2,s1,s2]}},AxesOrigin->{0,0},PlotMarkers->{Automatic,8},PlotRange->{{-5,5},{-5,5}},AspectRatio->1],{{s1,1/2},0,1},{{s2,1/2},0,1}]

GraphicsRow[{Plot3D[Ω[5,1,2,0,-1,Sqrt[3],-1,-Sqrt[3],x,y],{x,-4,4},{y,-4,4},BoxRatios->{1,1,1},AxesLabel->{"x","y","Ω"},BaseStyle->{FontFamily->"Helvetica",FontSize->12},Mesh->False,PlotPoints->25,PlotRange->All,ViewPoint->{1.3,-7,1.9}],Plot3D[Ω3[3,0.8,1,x,y],{x,-5,5},{y,-5,5},BoxRatios->{1,1,1},AxesLabel->{"x","y","Ω"},BaseStyle->{FontFamily->"Helvetica",FontSize->12},Mesh->False,PlotPoints->25,PlotRange->All,ViewPoint->{1.3,-7,1.9}],Plot3D[Ω3[3,0.7,0.3,x,y],{x,-5,5},{y,-5,5},BoxRatios->{1,1,1},AxesLabel->{"x","y","Ω"},BaseStyle->{FontFamily->"Helvetica",FontSize->12},Mesh->False,PlotPoints->25,PlotRange->All,ViewPoint->{1.3,-7,1.9}]},Spacings->Scaled[-0.1],ImageSize->1000]