Hi Pradeep,
Sorry for the delay in replying. To calculate the time series directly from the photoacoustic initial value p0 using 'Model II', as that paper calls it, it is necessary to deal somehow with the singularity due to the denominator in Eqs. (30) and (33). That's the hard part. In the paper we do so by only considering cylindrically symmetric p0 functions, and solve it in cylindrical coordinates (in which the singularity disappears).
That didn't answer your question on interpolation though. If you define a grid in x and t as the positions and time at which you want to know the photoacoustic signals, then you can calculate the values of kx and omega that they correspond to. From those you can calculate the values of (kx, zeta) at which you need to know the initial pressure, and you can interpolate from the values of p(kx,kz) that you know (because you know p(x,z)). Hope that helps.
(For others reading this post, the model we're referring to isn't implemented in k-Wave.)
Kind regards,
Ben