In this article, we study several reconstruction methods for the inverse source problem of photoacoustic tomography (PAT) with spatially variable sound speed and damping. The backbone of these methods is the adjoint operators, which we thoroughly analyze in both the $L^2$- and $H^1$-settings. They are casted in the form of a nonstandard wave equation. We derive the well-pawedness of the aforementioned wave equation in a natural functional space, and also prove the finite speed of propagation. Under the uniqueness and visibility condition, our formulations of the standard iterative reconstruction methods, such as Landweber's and conjugate gradients (CG), achieve a linear rate of convergence in either $L^2$- or $H^1$-norm. When the visibility condition is not satisfied, the problem is severely ill-posed and one must apply a regularization technique to stabilize the solutions. To that end, we study two classes of regularization methods: (i) iterative, and (ii) variational regularization. In the case of full data, our simulations show that the CG method works best; it is very fast and robust. In the ill-posed case, the CG method behaves unstably. Total variation regularization method (TV), in this case, significantly improves the reconstruction quality.