3D disk galaxy models

[rmandelb]

For some Euclid simulations, the group at Oxford (Lance Miller, Niraj Welikala) has been building 3D bulge+disk galaxies. The bulge can be represented by GalSim already, but the 3D disk cannot be. They model it as a disk of finite thickness when viewed edge-on, with a sech^2(z) surface brightness distribution as observed for real edge-on galaxies (the radial part of the profile is indeed exponential as you move along the disk at fixed z, so it's the sech^2(z) aspect that would be new to GalSim). Initially they did this via photon-shooting but they recently began representing it using DFTs. Lance has offered some C code that could be used to implement the 3D disk models in GalSim. He does not plan to implement it in GalSim himself.

Unlike many of our profiles, the input parameter would be the disk inclination, with the ellipticity being a derived parameter. I'm mentioning this not as a problem, but rather as an interesting feature.

Any thoughts, questions, concerns, etc.?

[brgillis]

Unless someone already familiar with working on GalSim is keen to pick this up, I can try implementing this.

The use of inclination as an input parameter is actually more convenient for simulations like what I'm working on than ellipticity is. It might also be a good idea to implement a similar class for bulges, so inclination can be passed directly to GalSim rather than requiring a calculation of the resulting ellipticity.

[rmandelb]

If you want to give this a try, then my suggestion is that you go through the C++ layer and look at the examples of base class implementations, come up with a proposal for how/where to add this in the C++ and python layers, and post it here. I'm happy to iterate with you and make suggestions about code design before you start trying to actually incorporate it.

[brgillis]

Alright, I've given this some thought about how to implement it now. What we have is a profile with an analytic form in Fourier space, but no analytic form in real space (aside from of course the Fourier transform integral). I can think of a couple ways to implement this. In both ways, implement in the C++ layer (like SBExponential.h/.cpp) and use the analytic form for kValue and fillKValue.

1. For xValue, perform the Fourier integral each time the method is called to calculate for that point. For fillXValue, perform an FFT on an array of k values.
2. Store a lookup table of x values (if fillXValue is called first, store a copy of the matrix made here, otherwise create one with an FFT of k values), and interpolate when xValue is called.

Option 2 will be much faster (but less accurate) if there are frequent calls to xValue. The question is, how often is xValue called? I suspect it isn't that often, and so for most cases, we'd want to err on the side of better accuracy at a slower speed. What are your thoughts on this, Rachel?

[rmandelb]

I would like some feedback from @rmjarvis and @jmeyers314 on this. But for now, here are my thoughts:

• yes, analytic form for kValue and fillKValue. This is useful as it means we get something exact when we're using DFT rendering.
• the way I normally think of kValue and fillKValue vs. xValue and fillXValue is that the former is required for DFT rendering of basic profiles, while the latter is needed for photon-shooting of basic profiles. (By "basic profiles" I mean one term of a convolution, like a profile for a galaxy or for a PSF. i.e., convolutions of a galaxy with a PSF don't tend to be analytic in real-space, but that's not a problem for photon-shooting, because you determine the photon offsets for each element of the convolution and then just add them.) The reason we have a lot of infrastructure built around tabulating the Hankel transforms needed to write a Sersic profile in k-space is that most of our image rendering is done via DFT, which means that we needed a fast yet accurate way to do this. That infrastructure includes the inclusion in the galsim.GSParams of tunable parameters that set the speed vs. accuracy tradeoff. It also includes the possibility of caching the table so later Sersic profiles that get initialized with the same n can access the table and don't have to rebuild it. This option is similar to your option (2) if you swap x<->k, but more powerful given the tunability of the speed vs. accuracy tradeoff.
• If I haven't over-simplified the situation regarding the need for kValue and xValue (hopefully @rmjarvis or @jmeyers314 will comment), then it seems to me that for the default drawing method, we can simply have xValue be undefined and set isAnalyticX to False. If the profile is initialized and then the user tells it to render the image via photon-shooting, it can at that point build the lookup table needed for xValue and fillXValue (and cache it).

[brgillis]

For reference of others reading it, the analytic form in Fourier space is simply the FT of the exponential profile multiplied by: pihk_ysin(i) / (2 * sinh( pihk_ysin(i)/2 )), where h is the scale height, i is the inclination angle (0 = edge-on), and k_y is the k component in the direction of the minor axis.

[brgillis]

I have some free time this afternoon, so I'm going to go ahead with starting up a first pass of implementing it (or at least some of the code for it), using the Sersic profile as a base. I'll let you know once I have something substantial, so we can make sure it's working properly.

[rmjarvis]

If I haven't over-simplified the situation regarding the need for kValue and xValue (hopefully @rmjarvis or @jmeyers314 will comment), then it seems to me that for the default drawing method, we can simply have xValue be undefined and set isAnalyticX to False.

I agree with this. The easiest path forward would be to have this profile be defined only in k-space. (cf. SBDeconvolve and SBFourierSqrt for two others like this.) If sometime later we decide we want to be able to photon shoot this, we can investigate options to enable that efficiently then.

If the profile is initialized and then the user tells it to render the image via photon-shooting, it can at that point build the lookup table needed for xValue and fillXValue (and cache it).

This would likely be the way to do so, but I'd leave that off at this point. It's fine to have this profile not be photon-shootable I think. At least for now. Enabling that can be a separate project (and separate issue) once the fft rendering is working.

[brgillis]

I've written up a description of my proposal for how to calculate the k-space expression of this profile, which you can see here: https://www.dropbox.com/s/z6h14bgd199czsi/Inclined_Exponential.pdf?dl=0 This includes Lance's initial work, plus my proposal for a coordinate transform which will allow the profile to be cached based on only one value. It will be a bit trickier than a Sersic profile to cache though, since it isn't axisymmetric, so it might indeed be a good idea to wait on implementing that for the future.

[rmandelb]

OK. Sounds like we're converging on just defining it in Fourier space for now. I think this is fine - and definitely a much simpler way to start. This will not be the first time we've introduced a new base class and only implemented photon shooting later on.

Is Lance willing to provide some example images that could be used as the basis for unit tests? In addition to internal unit tests (e.g., ones that make sure keyword arguments are being parsed properly), we often have unit tests that involve comparison of rendered images with some precomputed ones, either from an earlier version of GalSim or from some external software that we believe is a useful reference. You can see examples of this in tests/test_base.py (the tests that read in from tests/SBProfile_comparison_images). I think that if he can provide some images, it will not only be useful for unit testing, but will also provide basic sanity checks for you as you are working on the implementation of this profile.

Finally, do you have any questions about how to work on this in the GalSim repo? This wiki page will be useful:
https://github.com/GalSim-developers/GalSim/wiki/Developer-Information
and it points you to other documents within the repo that explain our branch naming scheme (all development on branches named for the issue number), etc. Another useful wiki page is here:
https://github.com/GalSim-developers/GalSim/wiki/Checklist-for-making-a-pull-request
It lists what you'd need to do to get the code ready to be reviewed.

[brgillis]

I don't have any questions right now. I've already read through the info you've linked. The only thing to note is that I don't have the ability to push my working branch to the server (unless I were to fork it). I'll either need to be added to the group I think, or just work locally.

I'll check with Lance to see if he has any samples that can be used for unit tests.

[rmandelb]

I just added you. You should have gotten the invitation from github.