Surface Density Profile
Firstly try and do as few computation as possible - e.g. the
np.logspace call is essentially constant, no reason to call it over
and over for every single element.

Also the dtype of object for surface_density makes using NumPy
somewhat useless for that array, try instead one dimension more and the
right datatype / dimensions of the array to actually save some space /
make things faster.

With NumPy you really want to avoid these loops and try to find some
operations that work on the whole array instead of per-element.  If
there's no way to do that easily, consider some of the apply functions
to at least make it look better and then look into processing parts of
the array in multiple threads / processes to speed things up as it looks
like to be a highly parallel computation.  Also Cython, but it's not clear to me if it would help that much here, even though trying it is not very complicated, so you could just do that and see.

Also avoid slower language features like the variable arguments in often called functions, e.g. here pass a fixed number of arguments instead of doing the *args thing.



Unrelated to NumPy and won't help here, but related to general Python, just in case for your next scripts:


There's also a typo in surface_denstiy_profile, it should be
density.
Don't use the range(..., len(...)) pattern, better would be to have
enumerate since you get both the index and value delivered
directly.