  The standard approach for photoacoustic imaging with variable speed of sound
is time reversal, which consists in solving a well-posed final-boundary value
problem for the wave equation backwards in time. This paper investigates the
iterative Landweber regularization algorithm, where convergence is guaranteed
by standard regularization theory, notably also in cases of trapping sound
speed or for short measurement times. We formulate and solve the direct and
inverse problem on the whole Euclidean space, what is common in standard
photoacoustic imaging, but not for time-reversal algorithms, where the problems
are considered on a domain enclosed by the measurement devices. We formulate
both the direct and adjoint photoacoustic operator as the solution of an
interior and an exterior differential equation which are coupled by
transmission conditions. The prior is solved numerically using a Galerkin
scheme in space and finite difference discretization in time, while the latter
consists in solving a boundary integral equation. We therefore use a BEM-FEM
approach for numerical solution of the forward operators. We analyze this
method, prove convergence, and provide numerical tests. Moreover, we compare
the approach to time reversal.
