Summation over Weierstrass $\wp$ functions

I've been trying to prove the following closed expression for a summation over Weierstrass $\wp$-functions:

\begin{equation} \sum_{k=1}^{N-1} \wp_N(k) = \frac{2}{\omega}\left(\zeta\left(\frac{\omega}{2}\right)-N\zeta_N\left(\frac{\omega}{2}\right)\right), \end{equation}

where $\wp_N$ is the usual WeierstrassP function with periods $(N,\omega)$, $\zeta_N$ is the Weierstrass zeta-function with periods $(N,\omega)$ and $\zeta$ is the Weierstrass zeta-function with periods $(1,\omega)$. Also, $N \in \mathbb{N}$, $N\geq 2$ and $\omega = i\pi/\kappa $ for some $\kappa\in \mathbb{R_{>0}}$.

My usual approach to such a problem is to find a quasi-periodic function $F$, i.e. satisfying $$F(z+1) = \alpha F(z),\qquad F(z+\omega)=F(z), $$ where $\alpha \in \mathbb{C}$ and $|\alpha| >0$. This $F$ should then have the left hand side of the expression above in its Laurent series. Secondly, I postulate a different function which has the exact same pole structure as $F$. Due to the Liouville theorem for elliptic functions, one can then conclude that they must be equal (quasi-periodicity dictates that their difference is not just a constant, but must be zero). Equating the Laurent-coefficients then should yield the equation given above.

In this case, one could use an appropriate modification of the function $F(z) = \sum_{k=0}^{N-1} \wp_N(z+k)$, which has as its zeroth order Laurent coefficient precisely $\sum_{k=1}^{N-1} \wp_N(k)$. This $F$ is doubly-periodic with periods $(1,\omega)$ and therefore does not satisfy quasi-periodicity yet. Its Laurent expansion equals $$ F(z) = \frac{1}{z^2} + \sum_{k=1}^{N-1} \wp_N(k)+ O(z) $$ and has therefore the exact same pole structure as the $\wp(z)$ (with periods $(1,\omega)$).

Does this approach work at all and if yes, which function $F$ should one use? If not, how could one prove the statement above?

In the mean time, I have found a function $F$ which might do the trick. Using the argument above, I proved that $$ F(z)=\sum_{k=0}^{N-1} e^{ipk}\wp_N(z+k) $$ equals $$ G(z) = -\frac{\sigma(z+r)}{\sigma(z-r)}e^{\frac{p}{\pi}\zeta(\omega/2)z}\left( \wp(z) -\wp(r) +\Delta(r)\left(\frac{\wp'(z)-\wp'(r)}{\wp(z)-\wp(r)} -\frac{\wp''(r)}{\wp'(r)} \right) \right) $$ where $r=-ip/(4\kappa)$, $\Delta(r) = \zeta(r) +\frac{p}{2\pi}\zeta(\omega/2)$ and $\wp$ and $\zeta$ are defined on the lattice $(1,\omega)$. Here $\sigma$ is the usual Weierstrass $\sigma$ function and is also defined on the lattice $(1,\omega)$. $p$ is a non-zero complex number. The expansion of $G$ around the point $z=0$ equals, writing $\delta=\frac{p}{\pi}\zeta\left(\frac{\omega}{2}\right)$ \begin{eqnarray}\label{laurent0} G(z) &=& \left( 1+2 \zeta(r)z+2\zeta^2(r) z^2 +O(z^3)\right)\left(1+\delta z +\frac{1}{2}(\delta z)^2\right) \nonumber \\ & &\times \left\{\frac{1}{z^2}-\wp(r) +\Delta\left( -\frac{2}{z} -\wp(r)z + \wp'(r) z^2 -\frac{\wp''(r)}{\wp'(r)}\right) \right\} \nonumber \\ &=& \frac{1}{z^2} + \frac{2\zeta(r) +\delta -2\Delta}{z} + \left(-\wp(r) -\Delta\frac{\wp''(r)}{\wp'(r)} +\left(\zeta(r) +\delta/2\right)\left(-4\Delta+2(\zeta(r) +\delta/2) \right) \right) +O(z). \nonumber \\ \end{eqnarray} One sees that the term going as $\frac{1}{z}$ vanishes due to our definition of $\Delta(r)$. By equating their respective Laurent series, I found $$ \sum_{k=1}^{N-1} e^{ipk}\wp_N(k) = -\wp(r) -\Delta(r) \frac{\wp''(r)}{\wp'(r)}-2\Delta(r)^2. $$ However, when I take the limit $p\rightarrow 0$ -- by expanding the right hand side in $p$ and taking the zeroth order term -- I find $$ \sum_{k=1}^{N-1}\wp_N(k) = \frac{2}{\omega}\zeta\left(\frac{\omega}{2}\right), $$ which is not exactly the answer I know to be true found on the top of this question. Also, numerical analysis of this expression shows conclusively that the equation above cannot be true, whereas the equation on the top of this page yields correct results every time.

For $\alpha, \beta \in \mathbb{C}$, let $L(\alpha,\beta)$ be the lattice $\{\; m\alpha+n\beta : m, n \in \mathbb{Z}\;\}$ and $[\alpha, \beta]$ be the line segment $\{\; t\alpha + (1-t)\beta : t\in\mathbb{R}, 0 \le t \le 1\;\}$ joining $\alpha$ and $\beta$. Consider the function $\sum\limits_{k=0}^{N-1}\wp_N(z+k)$, it has double poles on the lattice $L( 1, \omega)$ and for any $z$ near $\lambda \in L(1,\omega)$, we have

$$\sum_{k=0}^{N-1} \wp_N(z+k) \sim \frac{1}{(z-\lambda)^2} + O(1)$$

i.e. it has the same set of poles as $\wp(z)$ and the singular behavior at each pole is also the same. This means their difference is a doubly periodic entire function and hence is a constant. Let $K$ be that constant:

$$\sum_{k=0}^{N-1} \wp_N(z+k) - \wp(z) = K\tag{*1}$$

Notice for $z$ near $0$, we have:

$$\wp_{N}(z) = \frac{1}{z^2} + O(z^2)\quad\text{ and }\quad\wp(z) = \frac{1}{z^2} + O(z^2)$$ By taking the limit of $z \to 0$ in $(*1)$, we obtain:

$$\sum_{k=1}^{N-1} \wp_N(k) = K\tag{*2}$$

To fix $K$, pick a $\tau$ such that $[\tau-\frac{\omega}{2},\tau+\frac{\omega}{2}] \bigcap L(1,\omega) = \emptyset$ and integrate $(*1)$ along the line segment $[\tau-\frac{\omega}{2},\tau+\frac{\omega}{2}]$, we get

$$K\omega = \int_{\tau-\frac{\omega}{2}}^{\tau+\frac{\omega}{2}}\left( \sum_{k=0}^{N-1} \wp_N(z+k) - \wp(z)\right) dz$$

Using $\zeta'_N(z) = -\wp_N(z)$ and $\zeta'(z) = -\wp(z)$, this reduces to:

$$K\omega = -\sum_{k=0}^{N-1} \left(\zeta_N\left(\tau+k+\frac{\omega}{2}\right) - \zeta_N\left(\tau+k - \frac{\omega}{2}\right)\right) + \left(\zeta\left(\tau+\frac{\omega}{2}\right) - \zeta\left(\tau - \frac{\omega}{2}\right)\right)$$

As a function of $\tau$, $\zeta\left(\tau+\frac{\omega}{2}\right) - \zeta\left(\tau-\frac{\omega}{2}\right)$ is a constant because

$$\zeta'\left(\tau+\frac{\omega}{2}\right) - \zeta'\left(\tau-\frac{\omega}{2}\right) = -\wp\left(\tau+\frac{\omega}{2}\right) + \wp\left(\tau-\frac{\omega}{2}\right) = 0.$$

Taking the limit $\tau \to 0$ and using the fact $\zeta(\tau)$ is an odd function, we get

$$\zeta\left(\tau+\frac{\omega}{2}\right) - \zeta\left(\tau-\frac{\omega}{2}\right) = 2\zeta\left(\frac{\omega}{2}\right)\tag{*3a}$$

The same thing happens to $\zeta_N(\cdot)$ and we have

$$\zeta_N\left(\tau+\frac{\omega}{2}\right) - \zeta_N\left(\tau-\frac{\omega}{2}\right) = 2\zeta_N\left(\frac{\omega}{2}\right)\tag{*3b}$$ Substitute $(*3a)$ and $(*3b)$ into above expression for $K\omega$ and combine with $(*2)$, we get

$$K\omega = 2\left(\zeta\left(\frac{\omega}{2}\right) - N\zeta_N\left(\frac{\omega}{2}\right)\right) \quad\implies\quad \sum_{k=1}^{N-1} \wp_N(k) = \frac{2}{\omega}\left(\zeta\left(\frac{\omega}{2}\right) - N\zeta_N\left(\frac{\omega}{2}\right)\right) $$