Inverse of a structured matrix of sines. Suppose I have a matrix $P$ defined by
$$
P =\begin{pmatrix}
\sin(\frac{\pi}{n+1}) & \sin(\frac{2\pi}{n+1}) & \cdots & \sin(\frac{n\pi}{n+1}) \\
\sin(\frac{2\pi}{n+1}) & \sin(\frac{4\pi}{n+1}) & \cdots & \\
\vdots & & \ddots & \\
\sin(\frac{n\pi}{n+1}) & \cdots & & \sin(\frac{n^2\pi}{n+1})
\end{pmatrix} 
$$
and suppose I wish to find its inverse. 
I claim that $P^2 = \frac{n+1}{2} I_n$ and hence $P^{-1} = \frac{2}{n+1} P$ and indeed computing this in Matlab for a variety of $n$ it appears to be true, but I've been struggling to show this rigourously. I suspect I'm just very rusty with my trig manipulations but any help would be appreciated.
So far, I have that
$$ 
(P^2)_{kl} = \sum_{j=1}^n \sin(\frac{kj\pi}{n+1})\sin(\frac{lj\pi}{n+1})
$$
which looks somewhat similar to the stuff you get in fourier anyalysis where $\sin$ and $\cos$ form orthogonal polynomials, but I've been struggling to adapt it to this sum situation.
From here I've been trying to use that
$$ 
(P^2)_{kl} = \frac{1}{2}\sum_{j=1}^n \cos\left(\frac{(k-l)j}{n+1} \pi\right) - \cos\left(\frac{(k+l)j}{n+1} \pi\right)
$$
but I really haven't made much progress. Even the diagonal case where $k = l$ seems to not work out as immediately as I hoped, since I'm a bit lost on what to do with the second term. I suspect I am missing something obvious...

I managed to solve it. It was a bit of a headache but here is my proof.
For $k, l \in \mathbb{N}$ and $n \in \mathbb{N}$ consider the sum
\begin{equation}
A_{kl} =  \sum_{j=1}^n \sin \left(\frac{k \pi}{n+1}j\right) \sin \left(\frac{l \pi}{n+1}j\right).
\end{equation}
Claim:
\begin{equation}
A_{kl} = \frac{n+1}{2}\delta_{kl}.
\end{equation}
We start by noting that
\begin{equation}
A_{kl} = \frac{1}{2}  \sum_{j=1}^n \cos \left(\frac{(k-l)\pi}{n+1} j\right) - \cos \left(\frac{(k+l)\pi}{n+1} j\right)
\end{equation}
and letting $\theta_\pm = \frac{(k\pm l)\pi}{n+1}$
we have
\begin{equation}
A_{kl} = \frac{1}{2} \sum_{j=1}^n \cos( j\theta_- ) - \cos(j\theta_+ ).
\end{equation}
We now recall the following trigonometric identity due to Lagrange:
\begin{equation}
\sum_{j=1}^n \cos(j \theta) = \frac{1}{2} \left[\frac{\sin((n+\frac{1}{2})\theta)}{\sin(\frac{1}{2}\theta)}-1 \right]
\end{equation}
which holds when $\theta \neq 0$.
Note as well that
\begin{align*}
\sin((n+\frac{1}{2})\theta) &=  \sin((n+1)\theta - \frac{1}{2}\theta) \\
&= \sin((n+1)\theta)\cos(\frac{1}{2}\theta) - \cos((n+1)\theta)\sin(\frac{1}{2}\theta).
\end{align*}
and hence the Lagrange identity reduces to
\begin{equation}
\sum_{j=1}^n \cos(j \theta) =\frac{1}{2} \left[ \sin((n+1)\theta)\cot(\frac{1}{2}\theta) - \cos((n+1)\theta) - 1\right].
\end{equation}
We start by considering the diagonal case, when $k = l$. In this situation $\theta_- = 0$ and $\theta_+ = 2\pi k$. Then
\begin{align*}
A_{kk} &= \frac{1}{2} \sum_{j=1}^n (1 - \cos(j\theta_+)) \\
&= \frac{n}{2} - \frac{1}{4} \left[\frac{\sin((n+\frac{1}{2})\theta_+)}{\sin(\frac{1}{2}\theta_+)}-1 \right] \\
&= \frac{n}{2}-\frac{1}{4}\left[\sin(2k\pi)\cot(\frac{k\pi}{n+1}) - \cos(2 k \pi) - 1\right] \\
&= \frac{n}{2} - \frac{1}{4}\left[ -2 \right] \\
&= \frac{n+1}{2},
\end{align*}
since $\sin(2k\pi) = 0$ and $\cos(2k\pi) = 1$ for all $k$.
If we now consider $k \neq l$, then the same method follows but we have to expand both terms.
\begin{align*}
A_{kl} &= \frac{1}{4} \left[\left[\frac{\sin((n+\frac{1}{2})\theta_-)}{\sin(\frac{1}{2}\theta_-)}-1 \right] - \left[\frac{\sin((n+\frac{1}{2})\theta_+)}{\sin(\frac{1}{2}\theta_+)}-1 \right] \right] \\
&= \frac{1}{4}\left[\sin((k-l)\pi)\cot(\frac{1}{2}\theta_-) - \cos((k-l)\pi) - \sin((k+l)\pi)\cot(\frac{1}{2}\theta_+) + \cos((k+l)\pi)\right] \\
&= \frac{1}{4}\left[\cos((k+l)\pi) - \cos((k-l)\pi) \right]\\
&= \frac{1}{4} \left[\cos{k\pi}\cos{l\pi} - \sin{k\pi}\sin{l\pi} - \cos{k\pi}\cos{l\pi} - \sin{k\pi}\sin{l\pi}\right] \\
&= -\frac{1}{2}\sin{k\pi}\sin{l\pi} \\
&=0
\end{align*}
since $\sin(m \pi) = 0$ for all $m \in \mathbb{Z}$.
