Genus One with four ends

This surface has four embedded catenoidal ends and genus one. It is, however, not embedded.
There is a 1-parameter family with varying growth rates of the catenoidal ends.
All surfaces are symmetric with respect to the coordinate planes.

<<Own`Mesh`

Solving the period problem

a0[b_]:=(2 (Sqrt[2]) b+2 (Sqrt[2]) (b^3)+2 I b (Sqrt[1+(b^2)]) EllipticE[I ArcCsch[b],(b^2)]-2 I (b^2) (Sqrt[1+(b^2)]) EllipticE[-I ArcSinh[b],(1/(b^2))]-I b (Sqrt[1+(b^2)]) EllipticF[I ArcCsch[b],(b^2)]+I (b^3) (Sqrt[1+(b^2)]) EllipticF[I ArcCsch[b],(b^2)]-I (Sqrt[1+(b^2)]) EllipticF[-I ArcSinh[b],(1/(b^2))]+I (b^2) (Sqrt[1+(b^2)]) EllipticF[-I ArcSinh[b],(1/(b^2))]+(1+b) √(-(1+(b^2)) (-2 b-2 (b^3)+b ((EllipticE[I ArcCsch[b],(b^2)])^2)+(b^3) ((EllipticE[-I ArcSinh[b],(1/(b^2))])^2)+2 I (Sqrt[2]) b (Sqrt[1+(b^2)]) EllipticF[I ArcCsch[b],(b^2)]-2 I (Sqrt[2]) (b^2) (Sqrt[1+(b^2)]) EllipticF[I ArcCsch[b],(b^2)]+b ((EllipticF[I ArcCsch[b],(b^2)])^2)-(b^2) ((EllipticF[I ArcCsch[b],(b^2)])^2)-(b^3) ((EllipticF[I ArcCsch[b],(b^2)])^2)+(b^4) ((EllipticF[I ArcCsch[b],(b^2)])^2)+2 I (Sqrt[2]) b (Sqrt[1+(b^2)]) EllipticF[-I ArcSinh[b],(1/(b^2))]-2 I (Sqrt[2]) (b^2) (Sqrt[1+(b^2)]) EllipticF[-I ArcSinh[b],(1/(b^2))]+((EllipticF[-I ArcSinh[b],(1/(b^2))])^2)-b ((EllipticF[-I ArcSinh[b],(1/(b^2))])^2)-(b^2) ((EllipticF[-I ArcSinh[b],(1/(b^2))])^2)+(b^3) ((EllipticF[-I ArcSinh[b],(1/(b^2))])^2)+2 (b^2) EllipticE[-I ArcSinh[b],(1/(b^2))] (I (Sqrt[2]) (Sqrt[1+(b^2)])-(-1+b) EllipticF[I ArcCsch[b],(b^2)]-(-1+b) EllipticF[-I ArcSinh[b],(1/(b^2))])-2 b EllipticE[I ArcCsch[b],(b^2)] (I (Sqrt[2]) (Sqrt[1+(b^2)])+b EllipticE[-I ArcSinh[b],(1/(b^2))]-(-1+b) EllipticF[I ArcCsch[b],(b^2)]+EllipticF[-I ArcSinh[b],(1/(b^2))]-b EllipticF[-I ArcSinh[b],(1/(b^2))]))))/((-1+b) ((Sqrt[2])+(Sqrt[2]) (b^2)+I (Sqrt[1+(b^2)]) EllipticE[I ArcCsch[b],(b^2)]-I b (Sqrt[1+(b^2)]) EllipticE[-I ArcSinh[b],(1/(b^2))]-I (Sqrt[1+(b^2)]) EllipticF[I ArcCsch[b],(b^2)]-I b (Sqrt[1+(b^2)]) EllipticF[I ArcCsch[b],(b^2)]+I (Sqrt[1+(b^2)]) EllipticF[-I ArcSinh[b],(1/(b^2))]+I b (Sqrt[1+(b^2)]) EllipticF[-I ArcSinh[b],(1/(b^2))]));

Weierstrass Data

G[x_]:=(((Sqrt[1+x]) (-a+x) (Sqrt[-b+x]))/((Sqrt[-1+x]) (a+x) (Sqrt[b+x])))

om3[x_]:=(((-a+x) (a+x))/((-1+x) (1+x) (-b+x) (b+x)))

Numerical Integration (not used)

phi1[x_]:=(((x-1))^(-3/2)) (((x-a))^2) (1/(Sqrt[x-b])) (1/(Sqrt[x+1])) (((x+b))^(-3/2));

phi2[x_]:=(1/(Sqrt[x-1])) (((x-b))^(-3/2)) (((x+1))^(-3/2)) (((x+a))^2) (1/(Sqrt[x+b]));

om2[z_]:=(1/2) I (phi1[z]+phi2[z]);
om1[z_]:=-(1/2) (phi1[z]-phi2[z]);

f0[z_]:=Re[NIntegrate[{om1[x],om2[x],om3[x]},{x,I,z}]]

Symbolic Integration

phi1[x_]:=(-(b-x) (1+x) (b (-b (-2+x)+x)-2 a b (-1+b+2 x)+(a^2) (1+(b^2)-x+b x))-(b^2) (-(a^2) (-1+b)+4 a b+(-1+b) b)(Sqrt[1-x]) (Sqrt[1+x]) (Sqrt[1-(x/b)]) (Sqrt[1+(x/b)]) EllipticE[ArcSin[x],(1/(b^2))]+b (2 a-(a^2)+b) (-1+(b^2)) (Sqrt[1-x]) (Sqrt[1+x]) (Sqrt[1-(x/b)]) (Sqrt[1+(x/b)]) EllipticF[ArcSin[x],(1/(b^2))])/((-1+b) b (((1+b))^2) (Sqrt[-1+x]) (Sqrt[1+x]) (Sqrt[-b+x]) (Sqrt[b+x]))

phi2[x_]:=(-(-1+x) (b+x) (-2 a b (-1+b-2 x)+(a^2) (1+(b^2)+x-b x)+b (-x+b (2+x)))-(b^2) (-(a^2) (-1+b)+4 a b+(-1+b) b) (Sqrt[1-x]) (Sqrt[1+x]) (Sqrt[1-(x/b)]) (Sqrt[1+(x/b)])EllipticE[ArcSin[x],(1/(b^2))]+b (2 a-(a^2)+b) (-1+(b^2)) (Sqrt[1-x]) (Sqrt[1+x]) (Sqrt[1-(x/b)]) (Sqrt[1+(x/b)]) EllipticF[ArcSin[x],(1/(b^2))])/((-1+b) b (((1+b))^2) (Sqrt[-1+x]) (Sqrt[1+x]) (Sqrt[-b+x]) (Sqrt[b+x]));

f1[x_]:=(phi2[x]-phi1[x])/2;
f2[x_]:=I(phi2[x]+phi1[x])/2;

f3[x_]:=(1/(2 b (-1+(b^2))))((-1+(a^2)) b Log[-1+x]+(b-(a^2) b) Log[1+x]-((a^2)-(b^2)) (Log[-b+x]-Log[b+x]))

f0[x_]:=Re[{f1[x],f2[x],f3[x]}]

Choose parameter b>1

b=5.;
a=Re[a0[b]]

b=3.;
a=1.342358905985905`;

rr=-Sqrt[b];

f[x_]:=f0[rr(x+1)/(x-1)]

n[x_]:=StereographicProjection[G[rr(x+1)/(x-1)]]

Parametrization

ϵ=.00000000001;

Map polar coordinates in upper half plane to upper half disk Preimages of catenoidal ends:

{x1,x2}=((-rr-#)/(rr-#))& /@{1,b}

x1: handle catenoids
x2: big catenoids

{y1,y2}=(-(1/2) (#1+(1/(#1)))&)/@{x1,x2}

{i1,i2}=(((-y1+#1)/(-y2+#1))&)/@{1.,-1.}

ff[w_]:=-w+(Sqrt[-1+w]) (Sqrt[1+w]);
mo1[z_]:=((y2 z-1 y1)/(z-1));

mfac=2;

nx1=3*mfac;
nx2=2*mfac;
nx3=2*mfac;
nx4=2*mfac;
ny=4*mfac;

m1=PolarDomain[Exp[Union[
NRange[-8,Log[-i1],nx1],NRange[Log[-i1],0,nx2],NRange[0,Log[-i2],nx3],NRange[Log[-i2],3.5,nx4],SameTest->AlmostEqual]],NRange[ϵ,Pi-ϵ,ny]];

m2=MeshApply[ff[mo1[#]]&,m1];

mp1=MeshPlot3D[{f[x+I y],n[x+I y]},{x,y},m2];

mp2=MeshJoin[mp1,MeshReflect[mp1,Plane[{0,0,1},0]]];
mp3=MeshJoin[mp2,MeshReflect[mp2,Plane[{0,1,0},0]]];
mp4=MeshJoin[mp3,MeshReflect[mp3,Plane[{1,0,0},0]]];

Show[gl=Mesh3DToGraphics3D[mp4]]