Elliptic Fourier series are a good example to kick the tires on array programming systems; J and Python, however, are both dynamically typed.

A modern statically approach makes things significantly easier on the programmer; here I will use my Apple language, which aspires to the interactivity and conciseness of the APL family.

One has a principled yet safe version of Haskell's zip or J's " (re-rank on a dyadic verb):

 > :ty \f.\x.\y. fx y
((a → (b → c)) → (Arr (i Cons Nil) a → (Arr (i Cons Nil) b → Arr (i Cons Nil) c)))

This avoids the pitfalls of zip-ing different-length lists (one cannot write eqList = zipWith (==) in Haskell) and avoids deferring failures to runtime as in J.

One can inspect the type signature of cons:

 > :ty (<|)
(a → (Arr (i Cons Nil) a → Arr ((i + 1) Cons Nil) a))

There is a safe last:

 > :ty }.
(Arr ((i + 1) Cons Nil) a → a)

These work together:

 > :ty \x.\xs. }.(x<|xs)
(b → (Arr (i Cons Nil) b → b))

I believe this combination of interactivity and type inference in array programming will be quite productive.

To address the elliptic Fourier series, we focus on the offsets $$A_0$$ and $$C_0$$. This is straightforward to express with mathematical notation:

$$\displaystyle A_0 = \frac{1}{T} \sum_{p=1}^{K} \left(\frac{\Delta x_p}{2\Delta t_p}(t_p^2-t_{p-1}^2) + \xi_p(t_p-t_{p-1}) \right) \\ \displaystyle C_0 = \frac{1}{T} \sum_{p=1}^{K} \left(\frac{\Delta y_p}{2\Delta t_p}(t_p^2-t_{p-1}^2) + \delta_p(t_p-t_{p-1}) \right)$$

where

$$\xi_p = \displaystyle\sum_{j=1}^{p-1}\Delta x_j - \frac{\Delta x_p}{\Delta t_p}\sum_{j=1}^{p-1}\Delta t_j \\ \delta_p = \displaystyle\sum_{j=1}^{p-1}\Delta y_j - \frac{\Delta y_p}{\Delta t_p}\sum_{j=1}^{p-1}\Delta t_j$$

In Apple:

λxs.λys.
{ sum ⇐ [(+)/1 0 x]
; dxs ⟜ (-)\~xs; dys ⟜ (-)\~ys
; dts ⟜ [√(x^2+y^2)]dxs dys
; pxs ← (+)Λ 0 dxs; pys ← (+)Λ 0 dys; pts ⟜ (+)Λ 0 dts
; dtss ⟜ 0<|(-)\~((^2)'1 pts)
; T ⟜}.pts
; 𝜉 ← (-)pxs ((*)((%)dxs dts) pts)
; 𝛿 ← (-)pys ((*)((%)dys dts) pts)
; A ← ((sum ((*)((%)dxs dts) dtss))%2 + (sum ((*)𝜉 dts)))%T
; C ← ((sum ((*)((%)dys dts) dtss))%2 + (sum ((*)𝛿 dts)))%T
; (A,C)
}

This all typechecks! So we can be sure that there will not be off-by-one errors.