I want to test the hypothesis $H_0: X\sim N(\mu_0,\sigma_0^2)$ against $H_1:X\sim N(\mu_1,\sigma_1^2)$. By using Bayesian method or Neyman-Pearson, the rejection region has the form: $$\frac{n-1}{2}(\frac{1}{\sigma_1^2}-\frac{1}{\sigma_0^2})S^2+\frac{n}{2}(\frac{1}{\sigma_1^2}-\frac{1}{\sigma_0^2})\bar{X}^2-n(\frac{\mu_1}{\sigma_1^2}-\frac{\mu_0}{\sigma_0^2})\bar{X}<k,$$ where $\bar{X}=\frac{X_1+\dots+X_n}{n}$ and $S^2=\frac{1}{n-1}\sum_{i=1}^n(X_i-\bar{X})^2$. So to find the probability of type I error, I need to find the distribution of the LHS when $H_0$ is true. I know that if $X\sim N(\mu_0,\sigma_0^2)$ then $\frac{(n-1)S^2}{\sigma_0^2}\sim \chi^2(n-1)$, $\bar{X}\sim N(\mu_0,\frac{\sigma_0^2}{n})$, $S^2$ and $\bar{X}$ are independent. But I don't know how to compute the distribution of the LHS effectively. I tried to calculate the CDF of the LHS from the joint distribution of $(\bar{X},S^2)$ or use Jacobian matrix but it seems to be complicated.
Is there any way to know the exact distribution of the LHS? Or is there any way to approximate this distribution?