Help with deriving a log likelihood contribution for a panel model with fixed effects

[TheS]

I would appreciate some help with the following problem:

We have the following panel model: $$y_{it} = h(x_{it}\beta + c_i) + \epsilon_{it}, \quad t=1,\ldots,T, \quad i=1,\ldots,N$$$$c_i | x_i \text{ is }\mathcal{N}(0, \sigma_c^2)$$
Also,
$$\{\epsilon_{it}\}_{t=1}^{T} \text{ are independently distributed } \mathcal{N}(0, \sigma_{\epsilon}^2) \text{ over time conditional on} (x_i, c_i)$$
$$y_i := (y_{i1}, \ldots, y_{iT}) \text{ is a vector of outcomes }$$
$$x_i := (x_{i1}, \ldots, x_{iT}) \text{ gathers the possibly time-varying regressors }$$
$$x_{it} = (1, x_{it2}, x_{it3}) \text{ with } x_{it2} \text{ binary and } x_{it3} \text{ continuous }$$

$$\text{The model parameters involve } (\beta, \sigma_{\epsilon}, \sigma_c, h), \text{ where } \sigma_c \in \mathbb{R}_{++} \text{ and } h: \mathbb{R} \rightarrow \mathbb{R} \text{ is some unknown function.}$$
$$\text{Show that the log-likelihood contribution }\text{ for the model viewed as a function of } (\beta, \sigma_{\epsilon}, \sigma_c, h) \text{ equals}:$$
$$l_i(\beta, \sigma_{\epsilon}, \sigma_c, h) = -T\log(\sigma_{\epsilon}) -$$
$$\log \left[ \int_{-\infty}^{\infty} \exp\left( -\frac{1}{\sigma_{\epsilon}^2} \sum_{t=1}^{T}[(y_{it} - h(x_{it}\beta + \sigma_cc)]^2 \right) \phi(c) \, dc \right] - \frac{T}{2}\log(2\pi)$$
Any help would be much appreciated. My guess is that the first steps involve recognizing that the cond densities: $$f(y_{it} = h(x_{it}\beta + c_i) + \epsilon_{it}|h,\sigma_{\epsilon})= f(\epsilon_{it} =y_{it} - h(x_{it}\beta + c_i) + |h,\sigma_{\epsilon})$$ This we can write as a normal pdf:

$$\frac{1}{\sqrt{2\pi}\sigma_{\epsilon}} \exp\left(-\frac{\epsilon_{it}^2}{2\sigma_{\epsilon}^2}\right).$$
Then taking logarithms etc.. But I am still nowhere near the answer. Any expert help available?