  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.
