Combining an identity for harmonic numbers $H_n$, see the penultimate paragraph of Exercise 10 (page 8) of [1], and inspired in an integral, see identity $(22)$ of [2], I wrote using uniform convergence $$\int_0^1\frac{(1-x)\log^2(1-x)}{\sin (\pi x)}dx=2\sum_{n=1}^\infty\frac{H_{n}}{n+1}\int_0^1\frac{(1-x)x^{n+1}}{\sin (\pi x)}dx.\tag{1}$$
I know (using a CAS) some integral occuring in RHS of $(1)$, that are evaluated in terms of particular values of the zeta function $\zeta(s)$. But I think that evaluate for each integer $n\geq 1$ $$\int_0^1\frac{(1-x)x^{n+1}}{\sin (\pi x)}dx\tag{2}$$ is very difficult for me.
Question. Is it possible to evaluate $(2)$ with the purpose to write my identity $(1)$ as a statement for our integral in LHS of $(1)$ expressed as a series involving harmonic numbers and particular values of the Riemann's Zeta function? If these integrals in $(2)$ are known answer this question as a reference request, and I try to search such literature and read it. Thanks in advance.
References:
[1] Jack D’Aurizio, Superior Mathematics from an Elementary point of view, course notes, University of Pisa (2017-2018).
[2] Zurab Silagadze, Sums of Generalized Harmonic Series. For Kids from Five to Fifteen, RESONANCE (September 2015).
Given the ordinary generating function for the sequence $\{\zeta(2k)\}_{k\geq 1}$ and the identity $\frac{1}{\sin x}=\cot\frac{x}{2}-\cot x$ we have $$ \frac{\pi x}{\sin(\pi x)} = 1+2\sum_{m\geq 1}(-1)^{m+1}\zeta(2m)\left(\frac{2}{4^m}-1\right) x^{2m} \tag{1}$$ and since $\int_{0}^{1}\frac{(1-x)\log^2(1-x)}{\sin(\pi x)}\,dx = I = \int_{0}^{1}\frac{x\log^2(x)}{\sin(\pi x)}\,dx$ and $\int_{0}^{1}x^{2m}\log^2(x)\,dx = \frac{2}{(2m+1)^3}$ it follows that $$ I=\int_{0}^{1}\frac{x\log^2(x)}{\sin(\pi x)}\,dx = \frac{2}{\pi}+\frac{2}{\pi}\sum_{m\geq 1}(-1)^{m+1}\frac{\zeta(2m)}{(2m+1)^3}\left(\frac{4}{4^m}-2\right).\tag{2} $$ By exploiting $\zeta(2m)=\int_{0}^{+\infty}\frac{x^{2m-1}}{(2m-1)!}\cdot\frac{dx}{e^x-1}$ the RHS of $(2)$ can also be expressed as an integral over $\mathbb{R}^+$ of a function involving the product of $\frac{x}{e^x-1}$ and the difference of two $\phantom{}_2 F_3$ hypergeometric functions. Numerically $I\approx\frac{34}{49}$.