I have a mass density, which comes from measurement and is represented by a cartesian grid and the average density over each grid cell.

I need to do some computations with this density. The final computation output of the computation represents again cell averages over the same grid. However for the internals of the computation, I need to pass to a spherical coordinate system.

Numerical there are many variants of transforming a density. A property I would really like to have is that
transforming a density to spherical and back is very close to the original density.

The kind of algorithms I tried so far is to just transform cell centers and then use Interpolations.jl (e.g. Gridded(Linear()), Gridded(Constant())). However with these the error between a density and a density transformed back and forth is too big.

Any recommendations? Julia packages or algorithms?

Given a diffeomorphism f : X -> Y between manifolds, we want numerical approximations F, G of the pullbacks along f and g = f⁻¹. One hacky way to guarantee G ∘ F = id is to make the grid on X fine enough and using nearest neighbor interpolation. Of course this has a lot of other bad properties. In particular F∘G is a projection and not the identity.

It might be necessary to understand why the grid points must be different in the spherical coordinate system and the Cartesian system. From outside one could imagine that there might be solutions in which we keep the same points, and define interpolation functions for each of the coordinate systems if intermediate points are required, but without propagating the interpolation errors.

Yes good question. Some of the spherical processing is done by libraries, that work best on just an array of values and a uniform spherical grid. So if there is a better pair of transforms, that would be the path of least work for me.

Imagine there is a source of radiation at one point and I am interested in how much radiation arrives at each point in space. There are many places, where a spherical grid is advantageous, but here is one example. We need to know the amount of matter between the origin of radiation and each point in space. In cartesian coordinates, I do not know an efficient way to do this. In spherical coordinates, this is just a cumsum along the r axis.