Fix surface density, add cumulative surface density

[savchenkoyana]

No description provided.

[savchenkoyana]

Dear maintainers,
thank you for this project! It is very useful!

I think there is a bug with surface density because of the order of operations: Pi * i? smth : anything <-- here multiplication executes before the ternary operator, so Pi gets lost.
I fixed this one and also introduced cumulative surface density

[teuben]

thank you for the fix, we sadly don't have a full regression test, but I will check your PR this weekend and report here how it went.

[teuben]

For the record, this is how I compared the results:

mkplum p1M 1000000
gyrfalcON p1M . tstop=0.001 eps=0.05 kmax=7 manipname=projprof manipfile=junk1M_new.txt manippars=1,0,0,100

the dot for the output file simplifies re-running it, as data is not written to a file.

[teuben]

Indeed, a snippet of C code confirmed that the new expression is smaller by a factor of Pi.

[teuben]

If the new third column is a cumulative surface brightness, should it not converge to the total mass? Right now it seems it is printing the contribution in a shell, i.e. Pi.M.r^2

[savchenkoyana]

If the new third column is a cumulative surface brightness, should it not converge to the total mass? Right now it seems it is printing the contribution in a shell, i.e. Pi.M.r^2

PP.Mr(i) is M(<r), and PP.rad(i) is radius of bin:

https://github.com/teuben/nemo/blob/master/usr/dehnen/falcON/inc/public/profile.h#L119C1-L122C41

I made a small python script that compares Plummer density and parameters from projprof:

import numpy as np
import matplotlib.pyplot as plt

if __name__ == "__main__":
    data = np.loadtxt("junk1M_new.txt").T
    r_edges, sigma, sigma_cum = data[0], data[1], data[2]

    def plummer_sigma(r):
        return 1 / np.pi / (1 + r**2)**2

    def plummer_sigma_cum(r):
        return 1 / np.pi / (1 + r**2)

    # Compute bin centers from edges
    r_edges = np.hstack((np.zeros(1,), r_edges))
    r_centers = (r_edges[:-1] + r_edges[1:]) / 2.
    bin_widths = r_edges[1:] - r_edges[:-1]

    plt.bar(r_centers, sigma_cum, width=bin_widths, align="center", label="snapshot cum. surf. density", color="g")
    plt.plot(r_centers, plummer_sigma_cum(r_centers), label="theory cum. surf. density", color="y")
    plt.bar(r_centers, sigma, width=bin_widths, align="center", label="snapshot surf. density", color="r")
    plt.plot(r_centers, plummer_sigma(r_centers), label="theory surf. density", color="b")
    plt.yscale("log")
    plt.legend()

    plt.show()

Did I miss something? Or was the question about something else (cumulative mass, not cumulative density)?

[savchenkoyana]

I'm not very good with terms, what I wanted to add is M(<r) / (pi * r^2), where r is upper bound of a bin.

I don't know why I decided to put it in the same PR with bugfix. Maybe it is better to remove it for now? And then add it in another PR if it's useful

[teuben]

Cumulative mass makes more sense, since it's something that converges to the total mass. But because of binning and where one places the "r", it probably is not mathematically the same. Ideally it should be.

Your plummer_sigma_cum() function above is indeed the cumulative mass function, although I think there doesn't need to be a PI in there. It converges to 1 for R->large

Cumulative density is not something I've thought about being useful for anything, can you elaborate?

[savchenkoyana]

No, I think you're right, I removed it