General Neovius Surfaces (1|-1,-1)

d90[v_]:={v[[2]],-v[[1]],0}

nrml[v_]:=v/Sqrt[v.v]

For two points no vertically abov each other, find vertical plane through these points :

VerticalPlane[p1_,p2_]:=Module[{n},n=nrml[d90[p2-p1]];Plane[n,n.p2]]

Initialization

Choose triangle group

{r,s,t}={3,2,6};

{r,s,t}={6,2,3};

{r,s,t}={2,3,6};

{r,s,t}={2,4,4};

{r,s,t}={4,2,4};

{r,s,t}={3,3,3};

{a,b1,b2}={((-1+r)/r),-((-1+s)/s),-((-1+t)/t)};

Choose lattice constant

τ=.7I;

Period Problem

dh[z_]:=1

G[z_,p_,q1_,q2_]:=θ11[z-p,τ]^(-a)θ11[z+p,τ]^(a)θ11[z-(q1+τ/2),τ]^(-b1)θ11[z-(-q1+τ/2),τ]^(b1)θ11[z-(q2+τ/2),τ]^(-b2)θ11[z-(-q2+τ/2),τ]^(b2)

G[z_,p_,q1_]:=G[z,p,q1,((-2 q1 r (-1+s)+(1+2 p (-1+r)) s)/(2 (r+s)))]/G[0,p,q1,((-2 q1 r (-1+s)+(1+2 p (-1+r)) s)/(2 (r+s)))]

phi1[z_,a0_,b0_]:=G[z,a0,b0 ] dh[z];
phi2[z_,a0_,b0_]:=1/G[z,a0,b0]dh[z];

ω2[z_,a0_,b0_] :=I( phi2[z,a0,b0]+phi1[z,a0,b0] )/2;
ω1[z_,a0_,b0_] :=( phi2[z,a0,b0]-phi1[z,a0,b0] )/2;

test[a0_,b0_]:=Block[{z},{Re[NIntegrate[Evaluate[ω2[z,a0,b0]],{z,0,.5+τ/2}]],Re[NIntegrate[Evaluate[{ω1[z,a0,b0],ω2[z,a0,b0]}].cis[Pi/2-Pi/r],{z,τ/2,.5}]]}]

{p,q1}={a0,b0}/.FindRoot[ww=test[a0,b0],{a0,.25},{b0,.34},EvaluationMonitor:>Print[ww,"--  a0:",a0,"  b0:",b0]]

q2=((-2 q1 r (-1+s)+(1+2 p (-1+r)) s)/(2 (r+s)))

Surface Generation

G0[z_]:=θ11[z-p,τ]^(-a)θ11[z+p,τ]^(a)θ11[z-(q1+τ/2),τ]^(-b1)θ11[z-(-q1+τ/2),τ]^(b1)θ11[z-(q2+τ/2),τ]^(-b2)θ11[z-(-q2+τ/2),τ]^(b2)

g0=G0[0]

G[z_]:=G0[z]/g0

dh[z_]:=1

phi1[z_]:=G[z ] dh[z];
phi2[z_]:=1/G[z]dh[z];

ω1[z_] :=( phi2[z]-phi1[z] )/2;
ω2[z_] :=I( phi2[z]+phi1[z])/2;
ω3[z_]:=dh[z]

ϵ=10^-9.;

cp=.25+τ/4.;

f0[w_]:=Re[{
		NIntegrate[ω1[z],{z,cp,w}],NIntegrate[ω2[z],{z,cp,w}],
	w}]

f0[w_]:=Re[{
		NIntegrate[ω1[z],{z,cp,w},WorkingPrecision->20],NIntegrate[ω2[z],{z,cp,w},WorkingPrecision->20],
	w}]

δ0=f0[q2+τ/2-ϵ I]

f[w_]:=f0[w]-δ0

n[w_]:=StereographicProjection[G[w]]

Surface Parametrization

xspec=Union[{0,p,q1,q2,.5}];

y1=Im[τ]/2.;
nx=5;(*number of r-lines for catenoidal ends*)

ny=12;(*number of t-lines for all ends *)

fff[{x_,y_}]:=NRange[x,y,nx];
xr=Union[Flatten[fff/@Partition[xspec,2,1]],SameTest->AlmostEqual];
domain=RectangularDomain[xr,NRange[ϵ,y1-ϵ,ny]];

sp1=VerticalPlane[f[q1+τ/2-ϵ I],f[q2+τ/2-ϵ I]];

sp2=VerticalPlane[f[p+ϵ I],f[q1+τ/2-ϵ I]]

fr1=MeshPlot3D[{f[x+I y],n[x+I y]},{x,y},domain];

fr2=MeshJoin[fr1,MeshReflect[fr1,Plane[{0,0,1},.5]]];

fr3=MeshJoin[fr2,MeshReflect[fr2,sp1]];

fr4=MeshJoin@@Table[MeshRotateZ[fr3,2Pi i/t,0],{i,0,t-1}];

Show[gr=Mesh3DToGraphics3D[fr4],Axes->True,AxesLabel->{x,y,z}]