W1 = ellperiods([1,I]);
W2 = ellperiods([1,exp(2/3*Pi*I)]);
hW1 = [1/2, I/2, 1/2+I/2];
hW2 = [1/2, exp(2/3*Pi*I)/2, 1/2+exp(2/3*Pi*I)/2];
epsilon = 0.0000000000000000000000000000000000000000001+0.000000000000000000000000000000000001*I;

print("lattice 1, at w=0.3+0.3*I, z = ", hW1[1], ", f(z) = ",ellwp(W1,hW1[1],1)[2]/(ellwp(W1,hW1[1],1)[1]-ellwp(W1,0.3+0.3*I)))
print("lattice 1, at w=0.3+0.3*I, z = ", hW1[2], ", f(z) = ",ellwp(W1,hW1[2],1)[2]/(ellwp(W1,hW1[2],1)[1]-ellwp(W1,0.3+0.3*I)))
print("lattice 1, at w=0.3+0.3*I, z = ", hW1[3], ", f(z) = ",ellwp(W1,hW1[3],1)[2]/(ellwp(W1,hW1[3],1)[1]-ellwp(W1,0.3+0.3*I)))

print("lattice 2, at w=0.3+0.3*I, z = ", hW2[1], ", f(z) = ",ellwp(W2,hW2[1],1)[2]/(ellwp(W2,hW2[1],1)[1]-ellwp(W2,0.3+0.3*I)))
print("lattice 2, at w=0.3+0.3*I, z = ", hW2[2], ", f(z) = ",ellwp(W2,hW2[2],1)[2]/(ellwp(W2,hW2[2],1)[1]-ellwp(W2,0.3+0.3*I)))
print("lattice 2, at w=0.3+0.3*I, z = ", hW2[3], ", f(z) = ",ellwp(W2,hW2[3],1)[2]/(ellwp(W2,hW2[3],1)[1]-ellwp(W2,0.3+0.3*I)))

print("lattice 1, at w=0.3+0.3*I, z near 0, 1/f(z) = ",(ellwp(W1,epsilon,1)[1]-ellwp(W1,0.3+0.3*I))/ellwp(W1,epsilon,1)[2])
print("lattice 2, at w=0.3+0.3*I, z near 0, 1/f(z) = ",(ellwp(W1,epsilon,1)[1]-ellwp(W1,0.3+0.3*I))/ellwp(W1,epsilon,1)[2])