richteague / eddy

Python tools to recover precise rotation profiles of disks.

eddy - Extracting Disk Dynamics

What Is It?

eddy is a suite of Python tools to recover precise velocity profiles of protoplanetary disks from Doppler shifted line emission. eddy makes fitting of first moment maps and the inference of a rotation velocity from an annulus of spectra a breeze.

Installation

The most simple method is with pip,

pip install astro-eddy

The only real dependencies for this are numpy, scipy, matplotlib, and emcee, If you want to run the Gaussian Process method you will also need celerite which can be easily installed if you follow their installation guide.

If things have installed correctly you should be able to run the Jupyter Notebooks with no errors. If something goes wrong, please open an issue.

Useage

For guides on how to use eddy you will find extensive examples in the documents. We shamelessly recommend bettermoments to make the moment maps required for the fitting.

Citations

If you use eddy as part of your research, please cite the JOSS article:

@article{eddy,
    doi = {10.21105/joss.01220},
    url = {https://doi.org/10.21105/joss.01220},
    year = {2019},
    month = {feb},
    publisher = {The Open Journal},
    volume = {4},
    number = {34},
    pages = {1220},
    author = {Richard Teague},
    title = {eddy},
    journal = {The Journal of Open Source Software}
}

A full list of citations including dependencies can be found on the citations page.

Works Using eddy

A list of the works using (or at least citing) eddy can be found on ADS.

eddy's Issues

Weighting for v0 Map Fits

Include a 1/r^2 weighting for the fitting to account for larger number of pixels in the outer disk.

Needs tests.

Make some tests.

Plotting maxima in `plot_river` fails

Often (mostly) the plot_river function fails to properly find the correct maxima. Should check the gridding.

`fit_map` fails when only 1 parameter left free

As title.

Add optimize call before emcee

Need to run a quick optimize before the emcee run for the GP method.

Tutorial Notebooks

Tidy up and properly run the tutorial notebooks to include in the documentation.

Build models from percentiles

Currentlyevaluate_models only takes a posterior samples. This should be changed so it can also accept 16th, 50th and 84th percentiles to draw models from (assuming a simple Gaussian posterior).

Documentation

Make documentation available via ReadTheDocs.

'rotationmap' import error

Hello Dr. Teague! My name is Shreya Karthikeyan, and I am using the package eddy for my research at the University of Florida on protoplanetary disks. I've been able to utilize other packages of yours with no problem, but I've been running into an issue with the 'rotationmap' class in eddy. I installed eddy but rotationmap is simply not available when I list the available functions. Is this a common issue? Does it depend on the version of eddy that I installed?

Thank you in advance.

Priors are ignored when using `multiprocessing`.

As title.

Recommended change to default threading behavior

On my fork (commit 9f361f3), I've changed the default behavior of one of the emcee calls to use a specific multithreading package rather than be single-threaded.

Pros: runs faster if you have 'billiard' installed
Cons: eats up all your CPUs, though you can always 'nice' your python call. There may be a way to throttle it, but I haven't taken the effort

Feel free to pull the change if you like it.

Rotationmap fit_map method, error using "shadowed" parameter.

Hello Dr. Richard. I'm Francisco Urbina, an astronomy student from Chile, currently working as an intern for the NRAO. I was trying to use eddy for my research project in debris disks, but got an error I couldn't fix and I don't know if it is an usual one or maybe it's something different. I'll copy the code I've been using and the error on a txt file.

error.txt

In case you are wondering, I'm using python 3.9.7 and anaconda for package managenement.

Thanks in advance for your reply.

Multi-Ellipse Fitting

Implement multiple-ellipse fitting as is common in extra-galactic studies (e.g. https://ui.adsabs.harvard.edu/#abs/2006MNRAS.366..787K/abstract). This may allow us to distinguish between radial flows and warps.

"percentiles" return option currently returns 60th percentile instead of median

Line 267 of fit_cube.py currently reads
to_return += [np.percentile(samples, [16, 60, 84], axis=0)]

The 60 should presumably be 50.

Remove `tilt` Parameter

The parameter tilt should not be necessary as it is the difference between a positive and negative inclination. This should also help with the convergence of the fit (maybe).

Add self-gravity term

Use Eqn. 2 from https://arxiv.org/pdf/2104.09530.pdf to allow for self-gravity to be included in the fitting.

Allow for a GP to model covariance matrix

Posteriors are wildly underestimated, likely due to the poor (non-existent) handling of uncertainties. Try and implement a GP to model the covariance matrix and then marginalize over its hyper-parameters to see if that improves the returned posteriors.

gaussian_ensemble crashes when specifying an optical depth

When I try to create an ensemble of spectra using the gaussian_ensemble function, it works when when the optical depth tau is set to its default (None), but crashes whenever tau is set to a value. This is true regardless of the value. The error message that I get is

spectra,theta,velax = gaussian_ensemble(vrot=1500.,tau=5.0)

ValueError Traceback (most recent call last)
in
----> 1 spectra,theta,velax = gaussian_ensemble(vrot=1500.,tau=5.0)

/Volumes/Untitled2/eddy/eddy/modelling.py in gaussian_ensemble(vrot, Tb, dV, tau, rms, dV_chan, N, oversample, PAmin, PAmax, linear_sample, plot, return_ensemble)
68 spectra = np.array(spectra)
69 else:
---> 70 spectra = _thick_line(velax[None, :], v_los[:, None], dV, Tex, tau)
71 else:
72 if oversample:

/Volumes/Untitled2/eddy/eddy/modelling.py in _thick_line(x, x0, dx, Tex, tau, N)
232 """Optically thick line profile."""
233 if not N:
--> 234 return Tex * (1. - np.exp(-_gaussian(x, x0, dx, tau, N=1)))
235 xx = np.linspace(x[0], x[-1], x.size * int(N))
236 yy = Tex * (1. - np.exp(-_gaussian(xx, x0, dx, tau, N=False)))

/Volumes/Untitled2/eddy/eddy/modelling.py in _gaussian(x, x0, dx, A, N)
224 if not N:
225 return A * np.exp(-np.power((x - x0) / dx, 2))
--> 226 xx = np.linspace(x[0], x[-1], x.size * int(N))
227 yy = A * np.exp(-np.power((xx - x0) / dx, 2))
228 return np.array([np.average(yy[i*N:(i+1)*N]) for i in range(x.size)])

~/anaconda2/envs/python3/lib/python3.7/site-packages/numpy/core/function_base.py in linspace(start, stop, num, endpoint, retstep, dtype)
122 if num > 1:
123 step = delta / div
--> 124 if step == 0:
125 # Special handling for denormal numbers, gh-5437
126 y /= div

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

fit_method='SNR' and fit_method='doublegauss' can give NaN uncertainties

Something in scipy.optimize.minimize?

missing default_parameters.yml on pip install

Hi Richard, I'm having my class use eddy in a problem set for a planet formation seminar but they were having problems with it. I just uninstalled and installed it again to see if I could reproduce their problem and I now find it broken too! It is missing the default_parameters.yml file. (If I get this file from github directly and copy it to the eddy directory in site-packages, all is fine but this is probably not what you want.)

Anyway, thanks for a very cool package!

Documentation for inferring rotation velocity from real data

One thing that would be helpful in the Jupyter notebooks is an example, or more explicit description, of how to estimate the rotation velocity from actual data. The gaussian_ensemble class is very nice for testing the code, but once the testing is done the next step is to apply this to real data. The information on how to do this is contained throughout the various notebooks, but a section that explicitly lays out how to estimate rotation velocities from actual data would be a huge benefit for anyone trying to apply this code to their data.

If it is easy, I would suggest using the TW Hya CO data that is already used in the 'Fitting Rotation Maps' notebook as an example in the 'Inferring Rotation Maps...' notebooks to explicitly demonstrate how someone could create the necessary spectra, theta, and velax arrays, how to use these arrays to create an annulus object, and then how to use the annulus object to infer the rotation velocity. Barring this, a simple step-by-step list (e.g. (1) create an array of spectra at a range of azimuthal angles from your data with shape N x M, (2) make an array of the azimuthal angles with shape N, etc.) would help someone that is skimming the notebooks looking for instructions on how to apply this code to their data.