Fourier series of elliptic F integral 
Question. How to expand $$\mathbf{F}(\theta\mid k)=\int_0^\theta\frac{dt}{\sqrt{1-k^2\sin^2t}} \text{ (where $|k|<1$ and $|\theta|<\pi$)}$$as a Fourier series w.r.t. $\theta$?

I noticed that the question boils down to finding the Fourier expansion of $$f(\theta)=\frac1{\sqrt{1-k^2\sin^2\theta}}$$ as we can change the order of integration and differentiation. We can see clearly that $f$ is even, so $$\int_{-\pi}^\pi f(\theta)\sin n\theta d\theta=0.$$
Also, $$\int_{-\pi}^{\pi}f(\theta)\cos(2n+1)\theta d\theta=2\int_0^\pi f(\theta)\cos(2n+1)\theta d\theta\\=2\int_0^\pi f(\theta)\cos((2n+1)(\pi-\theta))d\theta=0$$
So the only hard part left is $$\int_0^{\pi/2}\frac{\cos 2n\theta}{\sqrt{1-k^2\sin^2\theta}}d\theta.$$
For $n=0$ I can clearly see it's $\mathbf{K}(k)$, the elliptic K integral.
Clearly there is a polynomial $P_n(x)$ s.t. $\cos 2nx=P_n(\sin^2x)$ with $\deg P_n=n$, but note that this does not make the question boil down to evaluating $\displaystyle\int_0^{\pi/2}\frac{\sin^{2n}(t)}{\sqrt{1-k^2\sin^2 t}}dt$ because after evaluating this elliptic-like integral, which I ensure that it evaluates to a form of $a\mathbf{K}(k)+b\mathbf{E}(k)$ by partial fraction expansion, where $a,b\in\mathbb Q$, there is a hard finite summation left involving the coefficients of $P$.
When $n=1$, $a_n$ equals $2\mathbf{E}/k^2+(1-2/k^2)\mathbf{K}$.
I don't have any further thoughts.
If we cannot find the general term with $k$ varying, can we at least find the coefficients when $k^2=-1$?
 A: $$\int_0^\theta\dfrac{dt}{\sqrt{1-k^2\sin^2t}}$$
$$=\int_0^\theta\sum\limits_{n=0}^\infty\dfrac{(2n)!k^{2n}\sin^{2n}t}{4^n(n!)^2}~dt$$
$$=\int_0^\theta\sum\limits_{n=0}^\infty\dfrac{(2n)!k^{2n}C_n^{2n}}{16^n(n!)^2}~dt+\int_0^\theta\sum\limits_{n=1}^\infty\sum\limits_{m=1}^n\dfrac{(-1)^m(2n)!k^{2n}C_{n+m}^{2n}\cos2mt}{2^{4n-1}(n!)^2}~dt$$
(according to https://en.wikipedia.org/wiki/List_of_trigonometric_identities#Power-reduction_formulae)
$$=\left[\sum\limits_{n=0}^\infty\dfrac{((2n)!)^2k^{2n}t}{16^n(n!)^4}\right]_0^\theta+\left[\sum\limits_{n=1}^\infty\sum\limits_{m=1}^n\dfrac{(-1)^m((2n)!)^2k^{2n-1}\sin2mt}{16^n(n!)^2(n+m)!(n-m)!}\right]_0^\theta$$
$$=\sum\limits_{n=0}^\infty\dfrac{\left(\tfrac{1}{2}\right)_n\left(\tfrac{1}{2}\right)_nk^{2n}\theta}{(n!)^2}+\sum\limits_{m=1}^\infty\sum\limits_{n=m}^\infty\dfrac{(-1)^m\left(\tfrac{1}{2}\right)_n\left(\tfrac{1}{2}\right)_nk^{2n-1}\sin2m\theta}{(n+m)!(n-m)!}$$
$$=\sum\limits_{n=0}^\infty\dfrac{\left(\tfrac{1}{2}\right)_n\left(\tfrac{1}{2}\right)_nk^{2n}\theta}{(n!)^2}+\sum\limits_{m=1}^\infty\sum\limits_{n=0}^\infty\dfrac{(-1)^m\left(\tfrac{1}{2}\right)_{m+n}\left(\tfrac{1}{2}\right)_{m+n}k^{2m+2n-1}\sin2m\theta}{(2m+n)!n!}$$
