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...

In fact, the matrix you are working on is one of the 4 versions of the matrix of the DST (Discrete Sine Transform). 
Please refer in the following document 
https://pdfs.semanticscholar.org/6088/357d5bf6099f83308023c59d1dc18b030e6a.pdf to definition (4.2a) with $N=n+1$ and a (harmless) normalization by $\sqrt{\dfrac{2}{N}}=\sqrt{\dfrac{2}{n+1}}$ allowing the author to give your property $P^{-1}=\dfrac{2}{n+1}P$ under the plain form $S^{-1}=S$ (property (4.4a) of this text).
