By considering any one-dimensional time-homogeneous solvable diffusion process,this paper develops a complete analytical framework for computing the distribution of the last hitting time, to any level, and its joint distribution with the process value on any finite time horizon. Our formalism allows for regular diffusions with any type of endpoint boundaries. We exploit the inherent link between last and first hitting times. The simpler known formula for the marginal distribution of the last hitting time on an infinite time horizon is easily recovered as a special limit. Furthermore, we derive general formulae for each component of the joint distribution, i.e., the jointly continuous, the partly continuous (defective) and the jointly defective portions. By employing spectral expansions of the transition densities and the first hitting time distributions, our derivations culminate in novel general spectral expansions for both marginal and joint distributions of the last hitting time and the process value on any finite time horizon. An additional main contribution of this paper lies in the application of our general formulae, giving rise to newly closed-form analytical formulae for several solvable diffusions. In particular, we systematically derive analytical expressions for each portion of the marginal and joint distributions of the last hitting time and the process value on any finite time horizon, without and with imposed killing at one or two interior points, for Brownian motion, Brownian motion with drift (geometric Brownian motion), the squared Bessel , squared radial Ornstein-Uhlenbeck (CIR) and Ornstein-Uhlenbeck processes. Most of our formulae are given in terms of spectral series that are rapidly convergent and efficiently implemented. We demonstrate this by presenting some numerical calculations of marginal and joint distributions using accurately truncated series.