  We present a new algorithm for reconstructing an unknown source in
Thermoacoustic and Photoacoustic Tomography based on the recent advances in
understanding the theoretical nature of the problem. We work with variable
sound speeds that might be also discontinuous across some surface. The latter
problem arises in brain imaging. The new algorithm is based on an explicit
formula in the form of a Neumann series. We present numerical examples with
non-trapping, trapping and piecewise smooth speeds, as well as examples with
data on a part of the boundary. These numerical examples demonstrate the robust
performance of the new algorithm.
