Quiet[Remove["Global`*"],{Remove::rmnsm}];Print["Mathematica $Version = “",$Version,"”"];Print["Execution time = ",DateString[DateList[],{"Hour",":","Minute"," on ","DayNameShort"," ","Day"," ","MonthNameShort"," ","Year"}]];

h=297(360/127)-30-24;
w=210(360/127)-24-24;
r=w/4;
Clear["yy"];yy=yy/.Solve[(yy-r)^2==(h-r-yy)^2+(w/2-r)^2,yy][[1]];
Print["h=",h,"   w=",w,"   r=",r,"   yy=",yy];

Print[N[circleCentres={
	{w/2,h-r},
	{r,yy},
	{w-r,yy},
	{r,r},
	{w-r,r}
}]//MatrixForm];
charges=Transpose[Join[Transpose[circleCentres],{RandomSample[Floor[Length[circleCentres]/2]+Range[Length[circleCentres]]]}]];
p1=Map[Graphics[{Red,Thick,Circle[{#[[1]],#[[2]]},r]}]&,circleCentres];
p2=Graphics[Map[Text[Style[#[[3]],FontSize->60,Bold,Green],{#[[1]],#[[2]]},{0,0}]&,charges]];
pot[x_,y_,qx_,qy_,q_]:=q/Sqrt[(x-qx)^2+(y-qy)^2];
potgridx[x_,y_]=D[Total[Table[pot[x,y,charges[[i,1]],charges[[i,2]],charges[[i,3]]],{i,1,Length[charges]}]],x];
potgridy[x_,y_]=D[Total[Table[pot[x,y,charges[[i,1]],charges[[i,2]],charges[[i,3]]],{i,1,Length[charges]}]],y];
p3=StreamPlot[{potgridx[x,y],potgridy[x,y]},{x,0,w},{y,0,h}];
Show[p1,p2,p3,AspectRatio->h/w]