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. f`x 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.