diff --git a/MIGRATION.md b/MIGRATION.md new file mode 100644 index 0000000..1451570 --- /dev/null +++ b/MIGRATION.md @@ -0,0 +1,69 @@ +# Migration from pkg_resources to importlib + +## Summary + +This package has been updated to use modern Python packaging standards, replacing the deprecated `pkg_resources` API with `importlib.metadata` and `importlib.resources`. + +## Changes Made + +### 1. Version Detection +**Before:** +```python +from pkg_resources import get_distribution, DistributionNotFound +__version__ = get_distribution("pygedm").version +``` + +**After:** +```python +from importlib.metadata import version, PackageNotFoundError +__version__ = version("pygedm") +``` + +### 2. Resource File Access +**Before:** +```python +from pkg_resources import resource_filename +DATAPATH = os.path.dirname(resource_filename("pygedm", "spiral.txt")) +``` + +**After:** +```python +from importlib.resources import files +DATAPATH = str(files("pygedm")) +``` + +### 3. Updated Dependencies +- Added `importlib-metadata` for Python < 3.8 (backport) +- Added `importlib-resources` for Python < 3.9 (backport) +- Minimum Python version remains 3.8+ + +### 4. Modern Packaging +- Added `pyproject.toml` with PEP 621 metadata +- Updated `setup.py` with new dependencies +- Updated `requirements.txt` + +## Benefits + +1. **No more pkg_resources errors**: Works with setuptools 81.0+ and modern Python environments +2. **Standard library approach**: Uses built-in modules for Python 3.8+ +3. **Future-proof**: Follows current Python packaging best practices +4. **Faster imports**: `importlib.metadata` is more efficient than `pkg_resources` + +## Installation + +The package now requires Python 3.8 or later. Install as usual: + +```bash +pip install pygedm +``` + +For development: +```bash +pip install -e . +``` + +## Notes + +- The `pkg_resources` API was deprecated and is being removed from setuptools +- All functionality remains the same; only internal implementation changed +- No API changes for end users diff --git a/README.md b/README.md index cbe34c8..ff0ffee 100644 --- a/README.md +++ b/README.md @@ -3,35 +3,37 @@ [![Coverage](https://codecov.io/gh/FRBs/pygedm/branch/master/graph/badge.svg?token=TlBiPzD7DP)](https://codecov.io/gh/FRBs/pygedm)[![Documentation Status](https://readthedocs.org/projects/pygedm/badge/?version=latest)](https://pygedm.readthedocs.io/en/latest/?badge=latest) # PyGEDM -_Python bindings for the YMW16, NE2001 and YT2020 electron density models_ +_Python bindings for the YMW16, NE2001, NE2025 and YT2020 electron density models_ -This package is a Python interface to the YMW16 and NE2001 electron density models, and YT2020 halo model. +This package is a Python interface to the YMW16, NE2001, NE2025 electron density models, and YT2020 halo model. The Yao, Manchester and Wang (2017, [Astrophys. J., 835, 29](https://iopscience.iop.org/article/10.3847/1538-4357/835/1/29/meta); -[arXiv:1610.09448](https://arxiv.org/abs/1610.09448)) YMW16 electron density model, is written in C++, and the Cordes and Lazio +[arXiv:1610.09448](https://arxiv.org/abs/1610.09448)) YMW16 electron density model, is written in C++, and the Cordes and Lazio (2001, [arXiv:0207156](https://arxiv.org/abs/astro-ph/)) NE2001 model is written in FORTRAN. This package, PyGEDM, wraps these -two codes using [pybind11](https://pybind11.readthedocs.io/en/stable/intro.html) to make them usable from Python. Here, we have converted NE2001 to C++ using `f2c`. +two codes using [pybind11](https://pybind11.readthedocs.io/en/stable/intro.html) to make them usable from Python. Here, we have converted NE2001 to C++ using `f2c`. This package also provides an interface to the NE2025 package [arxiv:2602.11838](https://arxiv.org/abs/2602.11838) ### PyGEDM paper PyGEDM is detailed in [Price, D. C., Flynn, C., and Deller, A.](https://ui.adsabs.harvard.edu/abs/2021PASA...38...38P/exportcitation) (2021): Price, D. C., Flynn, C., and Deller, A., “A comparison of Galactic electron density models using PyGEDM”, _Publications of the Astronomical Society of Australia_, vol. 38, 2021. doi:[10.1017/pasa.2021.33](https://ui.adsabs.harvard.edu/abs/2021PASA...38...38P/exportcitation). -If you use PyGEDM, be sure to [reference](https://github.com/FRBs/pygedm#references) the NE2001 and YMW16 codes. +If you use PyGEDM, be sure to [reference](https://github.com/FRBs/pygedm#references) the underlying NE2001, NE2025, YMW16 and YT2020 codes. ### Web app -We provide a web app at https://apps.datacentral.org.au/pygedm/ +We provide a web app at https://apps.datacentral.org.au/pygedm/ The pygedm web app is kindly hosted by [Data Central](https://datacentral.org.au/). ### Usage -Some usage examples can be found in the [examples directory](https://github.com/telegraphic/pygedm/tree/master/examples). +Some usage examples can be found in the [examples directory](https://github.com/telegraphic/pygedm/tree/master/examples). ```python import pygedm # calculate DM at a given distance +DM, tau_sc = pygedm.dist_to_dm(204.0, -6.5, 200, method='ne2025') +DM, tau_sc = pygedm.dist_to_dm(204.0, -6.5, 200, method='ne2001p') DM, tau_sc = pygedm.dist_to_dm(204.0, -6.5, 200, method='ne2001') DM, tau_sc = pygedm.dist_to_dm(204.0, -6.5, 200, method='ymw16') @@ -54,6 +56,7 @@ The methods return astropy [Quantities](http://docs.astropy.org/en/stable/units/ import pygedm import astropy.units as u import astropy.coordinates as c + DM = u.Quantity(10.0, unit='pc cm^-3') ra, dec = c.Angle(23.0, unit='hourangle'), c.Angle('-43:00:02', unit='degree') sky_coords = c.SkyCoord(ra, dec, frame='icrs') @@ -68,7 +71,7 @@ print(tau_sc.to('ns')) ### Installation -Requires `pybind11`, `astropy`, `numpy`, `scipy`, a newish C compiler with C++11 support (Ubuntu 16.04+ default gcc will work), plus `f2c`. +Requires `pybind11`, `astropy`, `numpy`, `scipy`, a newish C compiler with C++11 support (Ubuntu 16.04+ default gcc will work), plus `f2c`. Basic installation is via `pip`: @@ -77,7 +80,7 @@ Basic installation is via `pip`: pip install pygedm ``` -or +or ``` pip install git+https://github.com/telegraphic/pygedm @@ -94,7 +97,7 @@ pip install . `f2c` can be installed via `apt-get f2c` in Ubuntu, or via `conda install -c conda-forge f2c` if you use conda. #### Mac OSX -For MacOS, you are best off using `conda` and getting `f2c` via `conda install -c conda-forge f2c`. +For MacOS, you are best off using `conda` and getting `f2c` via `conda install -c conda-forge f2c`. Some users have had success using `f2c` in macports, which installs into `/opt/local/include` and `/opt/local/lib`. If using macports, the install command is: @@ -108,31 +111,35 @@ Windows is not currently supported. #### Unit tests -To run unit tests, run `python setup.py test`. Note that these tests only check the Python bindings, +To run unit tests, run `python setup.py test`. Note that these tests only check the Python bindings, not the underlying C/FORTRAN source code (which is not supplied with unit tests). ### References If using PyGEDM in a journal article, please remember to cite the underlying electron density models: -[Cordes, J. M., & Lazio, T. J. W. (2002)](https://ui.adsabs.harvard.edu/abs/2002astro.ph..7156C/abstract), -_NE2001.I. A New Model for the Galactic Distribution of Free Electrons and its Fluctuations_, +[Cordes, J. M., & Lazio, T. J. W. (2002)](https://ui.adsabs.harvard.edu/abs/2002astro.ph..7156C/abstract), +_NE2001.I. A New Model for the Galactic Distribution of Free Electrons and its Fluctuations_, arXiv e-prints, astro-ph/0207156. -[Cordes, J. M., & Lazio, T. J. W. (2003)](https://ui.adsabs.harvard.edu/abs/2003astro.ph..1598C/abstract), -_NE2001. II. Using Radio Propagation Data to Construct a Model for the Galactic Distribution of Free Electrons_, +[Cordes, J. M., & Lazio, T. J. W. (2003)](https://ui.adsabs.harvard.edu/abs/2003astro.ph..1598C/abstract), +_NE2001. II. Using Radio Propagation Data to Construct a Model for the Galactic Distribution of Free Electrons_, arXiv e-prints, astro-ph/0301598. -[Yao, J. M., Manchester, R. N., & Wang, N. (2017)](https://ui.adsabs.harvard.edu/abs/2017ApJ...835...29Y/abstract), -_A New Electron-density Model for Estimation of Pulsar and FRB Distances_, +[Yao, J. M., Manchester, R. N., & Wang, N. (2017)](https://ui.adsabs.harvard.edu/abs/2017ApJ...835...29Y/abstract), +_A New Electron-density Model for Estimation of Pulsar and FRB Distances_, The Astrophysical Journal, Volume 888, Issue 2, id.105, Colume 835, id.29 -[Yamasaki, S., & Totani, T. (2020)](https://ui.adsabs.harvard.edu/abs/2019arXiv190900849Y/abstract), -_The Galactic Halo Contribution to the Dispersion Measure of Extragalactic Fast Radio Bursts_, +[Yamasaki, S., & Totani, T. (2020)](https://ui.adsabs.harvard.edu/abs/2019arXiv190900849Y/abstract), +_The Galactic Halo Contribution to the Dispersion Measure of Extragalactic Fast Radio Bursts_, The Astrophysical Journal, Volume 888, Issue 2, id.105 +[S.K. Ocker, J.M. Cordes (2026)](https://ui.adsabs.harvard.edu/abs/2026arXiv260211838O/abstract), +_NE2025: An Updated Electron Density Model for the Galactic Interstellar Medium_ +arXiv e-prints, astro-ph/2602.11838 + These are available in bibtex format in [references.bib](https://github.com/telegraphic/pygedm/references.bib), -and also as an [ADS library](https://ui.adsabs.harvard.edu/public-libraries/Ci6_0-TlSySPMLrHxTvhhw). +and also as an [ADS library](https://ui.adsabs.harvard.edu/public-libraries/Ci6_0-TlSySPMLrHxTvhhw). ## YMW16 C README @@ -199,5 +206,3 @@ Jumei Yao (yaojumei _@_ xao.ac.cn), Richard N Manchester 07 July 2002 To compile and execute the code, see [code.pdf](https://github.com/telegraphic/pygedm/blob/master/ne2001_src/code.pdf). - - diff --git a/app/app.py b/app/app.py index 3843a08..f3febc0 100644 --- a/app/app.py +++ b/app/app.py @@ -1,35 +1,47 @@ import dash -import dash_labs as dl +import dash_bootstrap_components as dbc import dash_core_components as dcc import dash_html_components as html -import dash_bootstrap_components as dbc - +import dash_labs as dl +import h5py +import numpy as np import plotly.express as px import plotly.graph_objects as pgo - -import numpy as np -from astropy.coordinates import Angle, SkyCoord -from astropy import units as u import xarray as xr -import h5py +from astropy import units as u +from astropy.coordinates import Angle, SkyCoord import pygedm # Load skymap data -dskymap = h5py.File('data/skymap.h5', mode='r') -skymap_dist = dskymap['dist'] +dskymap = h5py.File("data/skymap.h5", mode="r") +skymap_dist = dskymap["dist"] -skymap_data_ne = xr.DataArray(dskymap['ne2001/ddm'][:, ::2, ::2], dims=('distance_kpc', 'gb', 'gl'), - coords={'distance_kpc': skymap_dist, 'gl': dskymap['gl'][::2], 'gb': dskymap['gb'][::2]}, - attrs={'units': 'DM pc/cm3'}) -skymap_data_ymw = xr.DataArray(dskymap['ymw16/ddm'][:, ::2, ::2], dims=('distance_kpc', 'gb', 'gl'), - coords={'distance_kpc': skymap_dist, 'gl': dskymap['gl'][::2], 'gb': dskymap['gb'][::2]}, - attrs={'units': 'DM pc/cm3'}) +skymap_data_ne = xr.DataArray( + dskymap["ne2001/ddm"][:, ::2, ::2], + dims=("distance_kpc", "gb", "gl"), + coords={ + "distance_kpc": skymap_dist, + "gl": dskymap["gl"][::2], + "gb": dskymap["gb"][::2], + }, + attrs={"units": "DM pc/cm3"}, +) +skymap_data_ymw = xr.DataArray( + dskymap["ymw16/ddm"][:, ::2, ::2], + dims=("distance_kpc", "gb", "gl"), + coords={ + "distance_kpc": skymap_dist, + "gl": dskymap["gl"][::2], + "gb": dskymap["gb"][::2], + }, + attrs={"units": "DM pc/cm3"}, +) # APP SETUP app = dash.Dash(__name__, plugins=[dl.plugins.FlexibleCallbacks()]) app.title = "PyGEDM" -server = app.server +server = app.server theme_name = "minty" css_url = f"https://bootswatch.com/4/{theme_name}/bootstrap.css" @@ -38,16 +50,27 @@ app, tab_locations=["Plot", "Output", "Skymap", "Notes"], title="PyGEDM: Galactic Electron Density Models", - theme=css_url, figure_template=True, sidebar_columns=3, + theme=css_url, + figure_template=True, + sidebar_columns=3, ) + @app.callback( args=dict( model=tpl.new_dropdown(["NE2001", "YMW16"], label="Model"), - method=tpl.new_dropdown(["DM (pc/cm3) to Distance", "Distance (kpc) to DM"], label="Method", kind=dl.State), + method=tpl.new_dropdown( + ["DM (pc/cm3) to Distance", "Distance (kpc) to DM"], + label="Method", + kind=dl.State, + ), dmord=tpl.new_textbox(10, label=None, kind=dl.State), - nu=tpl.new_textbox(1.0, label='Frequency (GHz)', kind=dl.State), - coords=tpl.new_dropdown(["Galactic (gl, gb)", "Celestial (RA, DEC)"], label="Coordinates", kind=dl.State), + nu=tpl.new_textbox(1.0, label="Frequency (GHz)", kind=dl.State), + coords=tpl.new_dropdown( + ["Galactic (gl, gb)", "Celestial (RA, DEC)"], + label="Coordinates", + kind=dl.State, + ), x0=tpl.new_textbox("00:00:00.00", label=None, kind=dl.State), x1=tpl.new_textbox("00:00:00.00", label=None, kind=dl.State), go=tpl.new_button("Calculate", label=None), @@ -57,24 +80,24 @@ tpl.new_graph(location="Plot"), tpl.new_div(location="Output"), tpl.new_graph(location="Skymap"), - tpl.new_div(location="Notes") + tpl.new_div(location="Notes"), ], template=tpl, ) def callback(model, method, dmord, nu, coords, x0, x1, go, tab): print(f"{tab}") - + # Setup some error handling coord_error = False - freq_error = False + freq_error = False dmord_error = False - + try: nu = float(nu) except ValueError: nu = 1.0 freq_error = True - + if method == "DM (pc/cm3) to Distance": f = pygedm.dm_to_dist units = 1.0 * u.pc / (u.cm**3) @@ -83,10 +106,10 @@ def callback(model, method, dmord, nu, coords, x0, x1, go, tab): except ValueError: dmord = 10 * units dmord_error = True - xt = 'DM (pc / cm3)' - yt = 'Distance (kpc)' - xl = 'Dispersion measure' - yl = 'Distance' + xt = "DM (pc / cm3)" + yt = "Distance (kpc)" + xl = "Dispersion measure" + yl = "Distance" else: f = pygedm.dist_to_dm units = 1.0 * u.kpc @@ -95,124 +118,169 @@ def callback(model, method, dmord, nu, coords, x0, x1, go, tab): except ValueError: dmord = 10 * units dmord_error = True - #print(f, dmord) - yt = 'DM (pc / cm3)' - xt = 'Distance (kpc)' - xl = 'Distance' - yl = 'Dispersion measure' + # print(f, dmord) + yt = "DM (pc / cm3)" + xt = "Distance (kpc)" + xl = "Distance" + yl = "Dispersion measure" if coords == "Galactic (gl, gb)": try: - gl = Angle(x0, unit='degree') - gb = Angle(x1, unit='degree') - sc = SkyCoord(gl, gb, frame='galactic') + gl = Angle(x0, unit="degree") + gb = Angle(x1, unit="degree") + sc = SkyCoord(gl, gb, frame="galactic") except ValueError: - sc = SkyCoord(0*u.deg, 0*u.deg, frame='galactic') + sc = SkyCoord(0 * u.deg, 0 * u.deg, frame="galactic") coord_error = True else: try: - ra = Angle(x0, unit='hourangle') - dec = Angle(x1, unit='degree') + ra = Angle(x0, unit="hourangle") + dec = Angle(x1, unit="degree") sc = SkyCoord(ra, dec) except: - sc = SkyCoord(0*u.deg, 0*u.deg, frame='galactic') + sc = SkyCoord(0 * u.deg, 0 * u.deg, frame="galactic") coord_error = True - + print(sc.galactic.l, sc.galactic.b, dmord, f) - dout = f(sc.galactic.l, sc.galactic.b, dmord, method=model, nu=nu) - dout_ne = f(sc.galactic.l, sc.galactic.b, dmord, method='ne2001', nu=nu) - dout_ymw = f(sc.galactic.l, sc.galactic.b, dmord, method='ymw16', nu=nu) + dout = f(sc.galactic.l, sc.galactic.b, dmord, method=model, nu=nu) + dout_ne = f(sc.galactic.l, sc.galactic.b, dmord, method="ne2001", nu=nu) + dout_ymw = f(sc.galactic.l, sc.galactic.b, dmord, method="ymw16", nu=nu) # Make plots D = np.linspace(0.1, dmord.value) y_ne21 = np.zeros_like(D) - y_ymw = np.zeros_like(D) + y_ymw = np.zeros_like(D) for ii, d in enumerate(D): - d_ne21 = f(sc.galactic.l, sc.galactic.b, d * units, method='ne2001', nu=nu) - d_ymw = f(sc.galactic.l, sc.galactic.b, d * units, method='ymw16', nu=nu) + d_ne21 = f(sc.galactic.l, sc.galactic.b, d * units, method="ne2001", nu=nu) + d_ymw = f(sc.galactic.l, sc.galactic.b, d * units, method="ymw16", nu=nu) if method == "DM (pc/cm3) to Distance": - y_ne21[ii] = d_ne21[0].to('kpc').value - y_ymw[ii] = d_ymw[0].to('kpc').value + y_ne21[ii] = d_ne21[0].to("kpc").value + y_ymw[ii] = d_ymw[0].to("kpc").value else: y_ne21[ii] = d_ne21[0].value - y_ymw[ii] = d_ymw[0].value + y_ymw[ii] = d_ymw[0].value - #print(d, y) + # print(d, y) fig = pgo.Figure() - fig.add_trace(pgo.Scatter(x=D, y=y_ne21, mode='lines', name='NE2001')) - fig.add_trace(pgo.Scatter(x=D, y=y_ymw , mode='lines', name='YMW16')) - fig.update_layout(xaxis_title=xt, yaxis_title=yt, - title=f'ICRS: {sc.icrs.ra:2.2f} {sc.icrs.dec:2.2f} \t Galactic: {sc.galactic.l:2.2f} {sc.galactic.b:2.2f}') + fig.add_trace(pgo.Scatter(x=D, y=y_ne21, mode="lines", name="NE2001")) + fig.add_trace(pgo.Scatter(x=D, y=y_ymw, mode="lines", name="YMW16")) + fig.update_layout( + xaxis_title=xt, + yaxis_title=yt, + title=f"ICRS: {sc.icrs.ra:2.2f} {sc.icrs.dec:2.2f} \t Galactic: {sc.galactic.l:2.2f} {sc.galactic.b:2.2f}", + ) # SKYMAP - if model == 'NE2001': + if model == "NE2001": skymap_data = skymap_data_ne else: skymap_data = skymap_data_ymw - + print(skymap_data.shape) - skymap = px.imshow(animation_frame=0, img=skymap_data, origin='upper' ) + skymap = px.imshow(animation_frame=0, img=skymap_data, origin="upper") skymap["layout"].pop("updatemenus") - skymap.update_layout(xaxis_title='Galactic longitude (deg)', - yaxis_title='Galactic latitude (Deg)', - title=f'{model}: All-sky DM map' - ) - - l_wrapped = sc.galactic.l.wrap_at(Angle(180, unit='deg')).value - + skymap.update_layout( + xaxis_title="Galactic longitude (deg)", + yaxis_title="Galactic latitude (Deg)", + title=f"{model}: All-sky DM map", + ) + + l_wrapped = sc.galactic.l.wrap_at(Angle(180, unit="deg")).value + skymap.add_shape( - type='rect', - x0=l_wrapped-2, x1=l_wrapped+2, y0=sc.galactic.b.value-2, y1=sc.galactic.b.value+2, - xref='x', yref='y', - line_color='cyan' + type="rect", + x0=l_wrapped - 2, + x1=l_wrapped + 2, + y0=sc.galactic.b.value - 2, + y1=sc.galactic.b.value + 2, + xref="x", + yref="y", + line_color="cyan", ) ## TEXT OUTPUT hdr = [html.H2(method), html.H4(f"Sky coordinates: {sc}")] if coord_error: - hdr.append(dbc.Alert('Input coordinates invalid, please check', color='danger')) + hdr.append(dbc.Alert("Input coordinates invalid, please check", color="danger")) if freq_error: - hdr.append(dbc.Alert('Could not parse frequency input, please check.', color='danger')) + hdr.append( + dbc.Alert("Could not parse frequency input, please check.", color="danger") + ) if dmord_error: - hdr.append(dbc.Alert('Could not parse DM/distance input, please check.', color='danger')) + hdr.append( + dbc.Alert( + "Could not parse DM/distance input, please check.", color="danger" + ) + ) hdr = html.Div(hdr) table_header = [ html.Thead(html.Tr([html.Th(""), html.Th("YMW16"), html.Th("NE2001")])) ] row1 = html.Tr([html.Th(f"{xl}"), html.Td(f"{dmord}"), html.Td(f"{dmord}")]) - row2 = html.Tr([html.Th(f"{yl}"), html.Td(f"{dout_ymw[0]:2.4f}"), html.Td(f"{dout_ne[0]:2.4f}")]) - row3 = html.Tr([html.Th(f"Scattering timescale @ {nu} GHz"), html.Td(f"{dout_ymw[1]:2.4e}"), html.Td(f"{dout_ne[1]:2.4e}")]) + row2 = html.Tr([ + html.Th(f"{yl}"), + html.Td(f"{dout_ymw[0]:2.4f}"), + html.Td(f"{dout_ne[0]:2.4f}"), + ]) + row3 = html.Tr([ + html.Th(f"Scattering timescale @ {nu} GHz"), + html.Td(f"{dout_ymw[1]:2.4e}"), + html.Td(f"{dout_ne[1]:2.4e}"), + ]) table_body = [html.Tbody([row1, row2, row3])] gedm_out = html.Div([hdr, dbc.Table(table_header + table_body, bordered=True)]) - - notes = html.Div([html.H2("PyGEDM"), html.P(f"Version: {pygedm.__version__}"), - html.P(["Authors: Danny C. Price, Chris Flynn, Adam Deller"]), - html.P(["PyGEDM Documentation: ", html.A("https://pygedm.readthedocs.io", href='https://pygedm.readthedocs.io')]), - html.P(["Github: ", html.A("https://github.com/FRBs/pygedm", href='https://github.com/FRBs/pygedm')]), - html.P(["Journal article: A comparison of Galactic electron density models using PyGEDM ", - html.A("ADS", href='https://ui.adsabs.harvard.edu/abs/2021arXiv210615816P/abstract'), " | ", - html.A("arXiv", href='https://arxiv.org/abs/2106.15816') - ]), - html.H4("NE2001"), - html.P("Cordes, J. M., & Lazio, T. J. W. (2002),"), - html.P("NE2001.I. A New Model for the Galactic Distribution of Free Electrons and its Fluctuations, arXiv e-prints, astro-ph/0207156."), - - html.H4("YMW16"), - html.P("Yao, J. M., Manchester, R. N., & Wang, N. (2017),"), - html.P("A New Electron-density Model for Estimation of Pulsar and FRB Distances, ApJ, Volume 888, Issue 2, id.105, Colume 835, id.29"), - html.P(["Web app hosted by ", html.A("Data Central", href='https://datacentral.org.au/')]), - ]) - - return fig, gedm_out, skymap, notes + + notes = html.Div([ + html.H2("PyGEDM"), + html.P(f"Version: {pygedm.__version__}"), + html.P(["Authors: Danny C. Price, Chris Flynn, Adam Deller"]), + html.P([ + "PyGEDM Documentation: ", + html.A( + "https://pygedm.readthedocs.io", href="https://pygedm.readthedocs.io" + ), + ]), + html.P([ + "Github: ", + html.A( + "https://github.com/FRBs/pygedm", href="https://github.com/FRBs/pygedm" + ), + ]), + html.P([ + "Journal article: A comparison of Galactic electron density models using PyGEDM ", + html.A( + "ADS", + href="https://ui.adsabs.harvard.edu/abs/2021arXiv210615816P/abstract", + ), + " | ", + html.A("arXiv", href="https://arxiv.org/abs/2106.15816"), + ]), + html.H4("NE2001"), + html.P("Cordes, J. M., & Lazio, T. J. W. (2002),"), + html.P( + "NE2001.I. A New Model for the Galactic Distribution of Free Electrons and its Fluctuations, arXiv e-prints, astro-ph/0207156." + ), + html.H4("YMW16"), + html.P("Yao, J. M., Manchester, R. N., & Wang, N. (2017),"), + html.P( + "A New Electron-density Model for Estimation of Pulsar and FRB Distances, ApJ, Volume 888, Issue 2, id.105, Colume 835, id.29" + ), + html.P([ + "Web app hosted by ", + html.A("Data Central", href="https://datacentral.org.au/"), + ]), + ]) + + return fig, gedm_out, skymap, notes app.layout = dbc.Container(fluid=True, children=tpl.children) if __name__ == "__main__": - app.run_server(host='0.0.0.0', debug=False) \ No newline at end of file + app.run_server(host="0.0.0.0", debug=False) diff --git a/docs/conf.py b/docs/conf.py index b39102e..7140214 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -12,18 +12,27 @@ # import os import sys + +import sphinx from m2r import MdInclude from recommonmark.transform import AutoStructify -import sphinx + def monkeypatch(cls): - """ decorator to monkey-patch methods """ + """decorator to monkey-patch methods""" + def decorator(f): method = f.__name__ old_method = getattr(cls, method) - setattr(cls, method, lambda self, *args, **kwargs: f(old_method, self, *args, **kwargs)) + setattr( + cls, + method, + lambda self, *args, **kwargs: f(old_method, self, *args, **kwargs), + ) + return decorator + # workaround until https://github.com/miyakogi/m2r/pull/55 is merged @monkeypatch(sphinx.registry.SphinxComponentRegistry) def add_source_parser(_old_add_source_parser, self, *args, **kwargs): @@ -33,15 +42,16 @@ def add_source_parser(_old_add_source_parser, self, *args, **kwargs): args = args[1:] return _old_add_source_parser(self, *args, **kwargs) -sys.path.insert(0, os.path.abspath('..')) -autoclass_content = 'both' +sys.path.insert(0, os.path.abspath("..")) + +autoclass_content = "both" # -- Project information ----------------------------------------------------- -project = 'PyGEDM' -copyright = '2020' -author = 'D. C. Price' +project = "PyGEDM" +copyright = "2020" +author = "D. C. Price" # -- General configuration --------------------------------------------------- @@ -49,15 +59,20 @@ def add_source_parser(_old_add_source_parser, self, *args, **kwargs): # Add any Sphinx extension module names here, as strings. They can be # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom # ones. -extensions = ['recommonmark', 'sphinx.ext.autodoc', 'sphinx.ext.coverage', 'sphinx.ext.napoleon'] +extensions = [ + "recommonmark", + "sphinx.ext.autodoc", + "sphinx.ext.coverage", + "sphinx.ext.napoleon", +] # Add any paths that contain templates here, relative to this directory. -templates_path = ['_templates'] +templates_path = ["_templates"] # List of patterns, relative to source directory, that match files and # directories to ignore when looking for source files. # This pattern also affects html_static_path and html_extra_path. -exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store'] +exclude_patterns = ["_build", "Thumbs.db", ".DS_Store"] # -- Options for HTML output ------------------------------------------------- @@ -66,36 +81,36 @@ def add_source_parser(_old_add_source_parser, self, *args, **kwargs): # html_theme = "sphinx_rtd_theme" html_theme_options = { - 'logo_only': False, - 'display_version': False, - 'prev_next_buttons_location': 'bottom', - 'style_external_links': False, + "logo_only": False, + "display_version": False, + "prev_next_buttons_location": "bottom", + "style_external_links": False, # Toc options - 'collapse_navigation': True, - 'sticky_navigation': True, - 'navigation_depth': 4, - 'includehidden': True, - 'titles_only': False + "collapse_navigation": True, + "sticky_navigation": True, + "navigation_depth": 4, + "includehidden": True, + "titles_only": False, } # Add any paths that contain custom static files (such as style sheets) here, # relative to this directory. They are copied after the builtin static files, # so a file named "default.css" will overwrite the builtin "default.css". -html_static_path = ['_static'] +html_static_path = ["_static"] def setup(app): config = { # 'url_resolver': lambda url: github_doc_root + url, - 'auto_toc_tree_section': 'Contents', - 'enable_eval_rst': True, + "auto_toc_tree_section": "Contents", + "enable_eval_rst": True, } - app.add_config_value('recommonmark_config', config, True) + app.add_config_value("recommonmark_config", config, True) app.add_transform(AutoStructify) # from m2r to make `mdinclude` work - app.add_config_value('no_underscore_emphasis', False, 'env') - app.add_config_value('m2r_parse_relative_links', False, 'env') - app.add_config_value('m2r_anonymous_references', False, 'env') - app.add_config_value('m2r_disable_inline_math', False, 'env') - app.add_directive('mdinclude', MdInclude) + app.add_config_value("no_underscore_emphasis", False, "env") + app.add_config_value("m2r_parse_relative_links", False, "env") + app.add_config_value("m2r_anonymous_references", False, "env") + app.add_config_value("m2r_disable_inline_math", False, "env") + app.add_directive("mdinclude", MdInclude) diff --git a/examples/compare_to_ne2001.py b/examples/compare_to_ne2001.py index dc5af71..b4c618a 100644 --- a/examples/compare_to_ne2001.py +++ b/examples/compare_to_ne2001.py @@ -3,16 +3,19 @@ Compares output to that of pyne2001, plots galaxic contribution """ -import pyymw16 -import pyne2001 + +import os import numpy as np import pylab as plt +import pyne2001 +import pyymw16 + import hickle as hkl -import os -def plot_gplane(d, title='', cbar_label="DM [pc cm$^{-3}$]"): - plt.imshow(d.T[::-1, ::-1], extent=(-180, 180, -90, 90), cmap='magma') + +def plot_gplane(d, title="", cbar_label="DM [pc cm$^{-3}$]"): + plt.imshow(d.T[::-1, ::-1], extent=(-180, 180, -90, 90), cmap="magma") plt.ylim(-60, 60) cb = plt.colorbar() cb.set_label(cbar_label) @@ -22,13 +25,14 @@ def plot_gplane(d, title='', cbar_label="DM [pc cm$^{-3}$]"): plt.ylabel("gb [deg]") plt.minorticks_on() -if not os.path.exists('gedm.hkl'): + +if not os.path.exists("gedm.hkl"): N = 256 d_2001 = np.zeros([N, N]) - d_2016= np.zeros([N, N]) + d_2016 = np.zeros([N, N]) for ii in range(N): - print("%i / %i" % (ii+1, N)) + print("%i / %i" % (ii + 1, N)) for jj in range(N): l = float(ii) / N * 360 - 180 b = float(jj) / N * 90 - 45 @@ -36,20 +40,20 @@ def plot_gplane(d, title='', cbar_label="DM [pc cm$^{-3}$]"): dm, tau = pyymw16.dist_to_dm(l, b, 30000) d_2016[ii, jj] = dm.value - hkl.dump({'NE2001':d_2001, 'YMW16':d_2016}, 'gedm.hkl') + hkl.dump({"NE2001": d_2001, "YMW16": d_2016}, "gedm.hkl") else: - d = hkl.load('gedm.hkl') + d = hkl.load("gedm.hkl") plt.figure(figsize=(9, 9)) - plt.subplot(3,1,1) - plot_gplane(d['NE2001'], 'NE2001') + plt.subplot(3, 1, 1) + plot_gplane(d["NE2001"], "NE2001") - plt.subplot(3,1,2) - plot_gplane(d['YMW16'], 'YMW16') + plt.subplot(3, 1, 2) + plot_gplane(d["YMW16"], "YMW16") - plt.subplot(3,1,3) - d_delta = (d['YMW16'] - d['NE2001']) - plot_gplane(d_delta, 'Difference') + plt.subplot(3, 1, 3) + d_delta = d["YMW16"] - d["NE2001"] + plot_gplane(d_delta, "Difference") plt.xlabel("gl [deg]") plt.tight_layout() - plt.savefig('compare_to_ne2001.png') + plt.savefig("compare_to_ne2001.png") plt.show() diff --git a/examples/create_healpix_map.py b/examples/create_healpix_map.py index b5da537..ab5a3f4 100644 --- a/examples/create_healpix_map.py +++ b/examples/create_healpix_map.py @@ -1,15 +1,17 @@ """ # create_healpix_map.py -Create a healpix map showing the Galactic DM contribution. +Create a healpix map showing the Galactic DM contribution. """ + +import os + import healpy as hp import numpy as np -import pyymw16 import pylab as plt -import os +import pyymw16 -filename = 'ymw16_allsky.fits' +filename = "ymw16_allsky.fits" if os.path.exists(filename): print("File %s exists, skipping write" % filename) @@ -19,13 +21,13 @@ NSIDE = 128 pix_id = np.arange(hp.nside2npix(NSIDE)) gl, gb = hp.pix2ang(NSIDE, pix_id, lonlat=True) - d_2016 = np.zeros_like(pix_id, dtype='float32') + d_2016 = np.zeros_like(pix_id, dtype="float32") for ii in pix_id: dm, tau = pyymw16.dist_to_dm(gl[ii], gb[ii], 30000) d_2016[ii] = dm.value print("Writing to %s" % filename) - hp.write_map(filename, d_2016, coord='G') + hp.write_map(filename, d_2016, coord="G") -hp.mollzoom(np.log(d_2016), title='YMW16', cmap='magma') +hp.mollzoom(np.log(d_2016), title="YMW16", cmap="magma") plt.show() diff --git a/examples/plot_galaxy_ne.py b/examples/plot_galaxy_ne.py index 9fcf56d..e11d7ed 100644 --- a/examples/plot_galaxy_ne.py +++ b/examples/plot_galaxy_ne.py @@ -5,44 +5,57 @@ Plot electron density of galaxy using ymw16.ne_crd """ -import pygedm -import pylab as plt import numpy as np +import pylab as plt + +import pygedm # Create empty arrays for NE2001 and YMW16 model output d_2016 = np.zeros((100, 100)) d_2001 = np.zeros((100, 100)) -# Loop through pixels +# Loop through pixels for ii in range(100): for jj in range(100): # generate X, Y, Z in PC from indexes x, y, z = (ii - 50) * 400, (jj - 50) * 400, 0 - + # Compute EDM - d_2016[ii, jj] = pygedm.calculate_electron_density_xyz(x, y, z, method='ymw16').value - d_2001[ii, jj] = pygedm.calculate_electron_density_xyz(x, y, z, method='ne2001').value + d_2016[ii, jj] = pygedm.calculate_electron_density_xyz( + x, y, z, method="ymw16" + ).value + d_2001[ii, jj] = pygedm.calculate_electron_density_xyz( + x, y, z, method="ne2001" + ).value + - # Plot output images -plt.rcParams['font.size'] = 14 +plt.rcParams["font.size"] = 14 -plt.figure(figsize=(4*3*2, 3*3)) -plt.subplot(1,2,1) +plt.figure(figsize=(4 * 3 * 2, 3 * 3)) +plt.subplot(1, 2, 1) plt.title("NE2001 electron density model, Z=0 plane") -plt.imshow(np.log10(d_2001), extent=(-50*400/1e3, 50*400/1e3, -50*400/1e3, 50*400/1e3), clim=(-5, 0)) +plt.imshow( + np.log10(d_2001), + extent=(-50 * 400 / 1e3, 50 * 400 / 1e3, -50 * 400 / 1e3, 50 * 400 / 1e3), + clim=(-5, 0), +) plt.xlabel("X [kpc]") plt.ylabel("Y [kpc]") cbar = plt.colorbar() cbar.set_label("log10($N_e$)") -plt.subplot(1,2,2) +plt.subplot(1, 2, 2) plt.title("YMW16 electron density model, Z=0 plane") -plt.imshow(np.log10(d_2016), extent=(-50*400/1e3, 50*400/1e3, -50*400/1e3, 50*400/1e3), clim=(-5, 0)) +plt.imshow( + np.log10(d_2016), + extent=(-50 * 400 / 1e3, 50 * 400 / 1e3, -50 * 400 / 1e3, 50 * 400 / 1e3), + clim=(-5, 0), +) plt.xlabel("X [kpc]") plt.ylabel("Y [kpc]") cbar = plt.colorbar() cbar.set_label("log10($N_e$)") plt.savefig("plot_galaxy_ne.png") -plt.show() \ No newline at end of file +plt.show() diff --git a/examples/plot_los_dm.py b/examples/plot_los_dm.py index 6ad8a5d..b19ec93 100644 --- a/examples/plot_los_dm.py +++ b/examples/plot_los_dm.py @@ -6,17 +6,17 @@ of 10,000 parsecs. """ -import pyymw16 -import pylab as plt import numpy as np +import pylab as plt +import pyymw16 d = np.zeros((100, 100)) for ii in range(100): for jj in range(100): - gl, gb, dist = (ii-50) * 360./100, (jj-50) * 180./100, 10000 - dm, tau = pyymw16.dist_to_dm(gl, gb, dist) - d[ii, jj] = dm.value + gl, gb, dist = (ii - 50) * 360.0 / 100, (jj - 50) * 180.0 / 100, 10000 + dm, tau = pyymw16.dist_to_dm(gl, gb, dist) + d[ii, jj] = dm.value plt.figure(figsize=(8, 5)) plt.title("YMW16 Galactic DM (to 10,000 pc)") @@ -24,7 +24,7 @@ plt.ylim(-45, 45) plt.xlabel("gl [deg]") plt.ylabel("gb [deg]") -cbar = plt.colorbar(orientation='horizontal') +cbar = plt.colorbar(orientation="horizontal") cbar.set_label("DM [pc cm$^{-3}$]") plt.savefig("plot_los_dm.png") plt.show() diff --git a/ne21c/ne2001_src/compile_ne2001.py b/ne21c/ne2001_src/compile_ne2001.py index d2d5780..95f73e7 100755 --- a/ne21c/ne2001_src/compile_ne2001.py +++ b/ne21c/ne2001_src/compile_ne2001.py @@ -1,30 +1,46 @@ #!/usr/bin/env python import os + NE_SRC_PATH = os.path.dirname(os.path.abspath(__file__)) -FORTRAN_SRCS = ['dmdsm.NE2001.f', 'density.NE2001.f', 'neLISM.NE2001.f', - 'neclumpN.f', 'nevoidN.f', 'scattering98.f'] +FORTRAN_SRCS = [ + "dmdsm.NE2001.f", + "density.NE2001.f", + "neLISM.NE2001.f", + "neclumpN.f", + "nevoidN.f", + "scattering98.f", +] + +DATA_FILES = [ + "gal01.inp", + "ne_arms_log_mod.inp", + "ne_gc.inp", + "nelism.inp", + "neclumpN.NE2001.dat", + "nevoidN.NE2001.dat", +] -DATA_FILES = ['gal01.inp', 'ne_arms_log_mod.inp', 'ne_gc.inp', - 'nelism.inp', 'neclumpN.NE2001.dat', 'nevoidN.NE2001.dat'] def runcmd(cmd): - """ Run a command with os.system, printing command to screen""" + """Run a command with os.system, printing command to screen""" print("\n> " + cmd) os.system(cmd) + def list_data_files(): return [os.path.join(NE_SRC_PATH, df) for df in DATA_FILES] + def cleanup(): """ """ print("---- Running cleanup ----") - runcmd('rm *.o') - if os.path.exists('sgnFile.pyf'): - runcmd('rm sgnFile.pyf') + runcmd("rm *.o") + if os.path.exists("sgnFile.pyf"): + runcmd("rm sgnFile.pyf") -def compile_ne2001(fcomp='gfortran', ar_flags='rc'): +def compile_ne2001(fcomp="gfortran", ar_flags="rc"): orig_cwd = os.getcwd() try: os.chdir(NE_SRC_PATH) @@ -35,31 +51,36 @@ def compile_ne2001(fcomp='gfortran', ar_flags='rc'): print("Working directory: " + NE_SRC_PATH) print("Fortran compiler: " + fcomp) - AR_CMD = 'ar {ar_flags} libNE2001.a '.format(ar_flags=ar_flags) - RANLIB_CMD = 'ranlib libNE2001.a' + AR_CMD = "ar {ar_flags} libNE2001.a ".format(ar_flags=ar_flags) + RANLIB_CMD = "ranlib libNE2001.a" for fsrc in FORTRAN_SRCS: - fobj = fsrc.replace('.f', '.o') - runcmd("{fcomp} -O -std=gnu -c -o {fobj} {fsrc}".format(fcomp=fcomp, fobj=fobj, fsrc=fsrc)) - AR_CMD += '{fobj} '.format(fobj=fobj) + fobj = fsrc.replace(".f", ".o") + runcmd( + "{fcomp} -O -std=gnu -c -o {fobj} {fsrc}".format( + fcomp=fcomp, fobj=fobj, fsrc=fsrc + ) + ) + AR_CMD += "{fobj} ".format(fobj=fobj) runcmd(AR_CMD) runcmd(RANLIB_CMD) print("\n---- Generating F2PY shared object for DMDSM ----") - runcmd('f2py -m dmdsm -h sgnFile.pyf dmdsm.NE2001.f --overwrite-signature') - runcmd('f2py -c sgnFile.pyf dmdsm.NE2001.f -L./ -lNE2001 -m dmdsm') - runcmd('rm sgnFile.pyf') + runcmd("f2py -m dmdsm -h sgnFile.pyf dmdsm.NE2001.f --overwrite-signature") + runcmd("f2py -c sgnFile.pyf dmdsm.NE2001.f -L./ -lNE2001 -m dmdsm") + runcmd("rm sgnFile.pyf") print("\n---- Generating F2PY shared object for density ----") - runcmd('f2py -m density -h sgnFile.pyf density.NE2001.f --overwrite-signature') - runcmd('f2py -c sgnFile.pyf density.NE2001.f -L./ -lNE2001 -m density') - runcmd('rm sgnFile.pyf') + runcmd("f2py -m density -h sgnFile.pyf density.NE2001.f --overwrite-signature") + runcmd("f2py -c sgnFile.pyf density.NE2001.f -L./ -lNE2001 -m density") + runcmd("rm sgnFile.pyf") except: raise finally: os.chdir(orig_cwd) + if __name__ == "__main__": try: compile_ne2001() @@ -68,6 +89,3 @@ def compile_ne2001(fcomp='gfortran', ar_flags='rc'): raise finally: cleanup() - - - diff --git a/ne21c/test_ne21c.py b/ne21c/test_ne21c.py index 5815c06..159d739 100644 --- a/ne21c/test_ne21c.py +++ b/ne21c/test_ne21c.py @@ -1,16 +1,17 @@ -import ne21c import pprint + +import ne21c import numpy as np print(ne21c.__doc__) print(ne21c.dm_to_dist.__doc__) a = ne21c.dm_to_dist(0, 0, 100) -assert np.isclose(a['dist'], 2.068159818649292) +assert np.isclose(a["dist"], 2.068159818649292) pprint.pprint(a) print(ne21c.dist_to_dm.__doc__) -a = ne21c.dist_to_dm(0, 0, a['dist']) +a = ne21c.dist_to_dm(0, 0, a["dist"]) pprint.pprint(a) print(ne21c.density_xyz.__doc__) @@ -18,5 +19,4 @@ pprint.pprint(a) # Compare against pre-computed value -assert np.isclose(a['ne'], 10.04799611523049) - +assert np.isclose(a["ne"], 10.04799611523049) diff --git a/pygedm/__init__.py b/pygedm/__init__.py index 3a79f56..54761d4 100644 --- a/pygedm/__init__.py +++ b/pygedm/__init__.py @@ -1,9 +1,20 @@ -from pkg_resources import resource_filename, get_distribution, DistributionNotFound +try: + from importlib.metadata import version, PackageNotFoundError +except ImportError: + from importlib_metadata import version, PackageNotFoundError try: - __version__ = get_distribution('pygedm').version -except DistributionNotFound: # pragma: no cover - __version__ = '0.0.1 - manual install' + __version__ = version("pygedm") +except PackageNotFoundError: # pragma: no cover + # Package is not installed, use a fallback version + __version__ = "4.0.0" -from .pygedm import dm_to_dist, dist_to_dm, calculate_electron_density_lbr, calculate_electron_density_xyz, \ - calculate_halo_dm, generate_healpix_dm_map, convert_lbr_to_xyz \ No newline at end of file +from .pygedm import ( + calculate_electron_density_lbr, + calculate_electron_density_xyz, + calculate_halo_dm, + convert_lbr_to_xyz, + dist_to_dm, + dm_to_dist, + generate_healpix_dm_map, +) diff --git a/pygedm/healpix_utils.py b/pygedm/healpix_utils.py index 418c7f4..0834304 100644 --- a/pygedm/healpix_utils.py +++ b/pygedm/healpix_utils.py @@ -1,9 +1,9 @@ - def check_for_healpix_install(): - """ Returns true if python can import healpix """ + """Returns true if python can import healpix""" try: import healpy as hp + HAS_HEALPIX = True except ImportError: # pragma: no cover HAS_HEALPIX = False - return HAS_HEALPIX \ No newline at end of file + return HAS_HEALPIX diff --git a/pygedm/ne2001_wrapper.py b/pygedm/ne2001_wrapper.py index 974a253..b9daa80 100644 --- a/pygedm/ne2001_wrapper.py +++ b/pygedm/ne2001_wrapper.py @@ -14,22 +14,24 @@ arXiv e-prints, astro-ph/0301598. """ -import numpy as np import os from functools import wraps -from astropy import units as u + import ne21c +import numpy as np +from astropy import units as u DATA_PATH = os.path.dirname(os.path.abspath(__file__)) def run_from_pkgdir(f): - """ Decorator function to chdir() into package directory when running + """Decorator function to chdir() into package directory when running NE2001 code doesn't know the relative path to its data files. This wraps the function call, changing into the right directory first, calling it, then changing back to original directory. """ + @wraps(f) def wrapped(*args, **kwargs): cwdpath = os.getcwd() @@ -37,15 +39,16 @@ def wrapped(*args, **kwargs): os.chdir(DATA_PATH) r = f(*args, **kwargs) return r - except: # pragma: no cover + except: # pragma: no cover raise finally: os.chdir(cwdpath) + return wrapped def TAUISS(d, sm, nu): - """ Convert a scattering measure (SM) to scattering timescale at given frequency. + """Convert a scattering measure (SM) to scattering timescale at given frequency. Direct port from FORTRAN code scattering98.f @@ -74,13 +77,34 @@ def TAUISS(d, sm, nu): tauiss = 1000. * (sm / 292.)**1.2 * d * nu**(-4.4) end """ - tauiss = 1. * (sm / 292.)**1.2 * d * nu**(-4.4) + tauiss = 1.0 * (sm / 292.0) ** 1.2 * d * nu ** (-4.4) return tauiss +def TRANSITION_FREQUENCY(d_eff: float, sm: float, xi: float = None) -> float: + """The transition frequency between weak and strong scattering + + Python implementation of Equation 17 in C&L 2002 [1]. + + Computes nu = (318 GHz) * xi^(10/17) * sm^(6/17) * d_eff^(5/17) + + Args: + d_eff (float): effective distance to the scattering medium (kpc) + sm (float): Scattering measure (kpc m^{-20/3}) + xi (float): Fresnel scale scaling factor (~1). If not supplied, + sqrt(2 pi) is used as preferred definition. + + Returns: + nu (float): Transition frequency in GHz + """ + xi = np.sqrt(np.pi * 2) if xi is None else xi + nu = 318.0 * xi ^ (10 / 17) * sm ** (6 / 17) * d_eff ** (5 / 17) + return nu + + @run_from_pkgdir def dm_to_dist(l, b, dm, nu=1.0, full_output=False): - """ Convert DM to distance and compute scattering timescale + """Convert DM to distance and compute scattering timescale Args: l (float): galactic longitude in degrees @@ -99,18 +123,18 @@ def dm_to_dist(l, b, dm, nu=1.0, full_output=False): return 0.0 * u.pc, 0.0 * u.s else: d = ne21c.dm_to_dist(l_rad, b_rad, dm) - tau_sc = TAUISS(float(d['dist']), d['smtau'], nu=nu) + tau_sc = TAUISS(float(d["dist"]), d["smtau"], nu=nu) if not full_output: - return (float(d['dist']) * u.kpc).to('pc'), tau_sc * u.s + return (float(d["dist"]) * u.kpc).to("pc"), tau_sc * u.s else: - d['tau_sc'] = tau_sc + d["tau_sc"] = tau_sc return d @run_from_pkgdir def dist_to_dm(l, b, dist, nu=1.0, full_output=False): - """ Convert distance to DM and compute scattering timescale + """Convert distance to DM and compute scattering timescale Args: l (float): galactic longitude in degrees @@ -125,21 +149,22 @@ def dist_to_dm(l, b, dist, nu=1.0, full_output=False): l_rad = np.deg2rad(l) b_rad = np.deg2rad(b) - if np.isclose(dist, 0): # Catch infinite timeout bug + if np.isclose(dist, 0): # Catch infinite timeout bug return 0.0 * u.pc / u.cm**3, 0.0 * u.s else: d = ne21c.dist_to_dm(l_rad, b_rad, dist) - tau_sc = TAUISS(float(dist), d['smtau'], nu=nu) + tau_sc = TAUISS(float(dist), d["smtau"], nu=nu) if not full_output: - return float(d['dm']) * u.pc / u.cm**3, tau_sc * u.s + return float(d["dm"]) * u.pc / u.cm**3, tau_sc * u.s else: - d['tau_sc'] = tau_sc + d["tau_sc"] = tau_sc return d + @run_from_pkgdir def calculate_electron_density_xyz(x, y, z): - """ Compute electron density at Galactocentric X, Y, Z coordinates + """Compute electron density at Galactocentric X, Y, Z coordinates x,y,z are Galactocentric Cartesian coordinates, measured in kpc (NOT pc!) with the axes parallel to (l, b) = (90, 0), (180, 0), and (0, 90) degrees @@ -153,4 +178,4 @@ def calculate_electron_density_xyz(x, y, z): ne_out (astropy.Quantity): Electron density in cm-3 """ ne_out = ne21c.density_xyz(x, y, z) - return ne_out['ne'] / u.cm**3 + return ne_out["ne"] / u.cm**3 diff --git a/pygedm/ne2025_wrapper.py b/pygedm/ne2025_wrapper.py new file mode 100644 index 0000000..4402552 --- /dev/null +++ b/pygedm/ne2025_wrapper.py @@ -0,0 +1,118 @@ +""" +Wrapper from mwprop / ne2025 + +References: + [1] `S.K. Ocker, J.M. Cordes (2026) `_, + *NE2025: An Updated Electron Density Model for the Galactic Interstellar Medium*, + https://arxiv.org/abs/2602.11838 +""" + +import os +import numpy as np + +from astropy import units as u +from astropy.coordinates import Angle +from astropy.units import Quantity, Unit + +from mwprop.nemod.NE2025 import ne2025 +from mwprop.nemod.NE2001 import ne2001 +from mwprop.nemod.density import density_2001 + + +def dm_to_dist(gl, gb, dm, nu=1.0, model='ne2025', full_output=False): + """Convert a DM to a distance + + Args: + gl (float in deg): galactic longitude + gb (float in deg): galactic latitude + dm (float in pc/cm3): dispersion measure (pc cm^-3) + nu (float in GHz): Observing frequency + model (str): One of 'ne2001' or 'ne2025' + + Notes: + In this wrapper, 'ne2001' == 'ne2001p'. + + Returns: + dist (astropy.Quantity), tau_sc (astropy.Quantity): distance (pc) and scattering time scale (s) + """ + if np.isclose(dm, 0): # WAR Catch spline failure + return 0.0 * u.pc, 0.0 * u.s + else: + if model.lower() == 'ne2025': + _Dk,Dv,_Du,_Dd = ne2025(gl, gb, dm, ndir=1, dmd_only=False, classic=False) + else: + _Dk,Dv,_Du,_Dd = ne2001(gl, gb, dm, ndir=1, dmd_only=False, classic=False) + dist, sm_tau = Dv['DIST'], Dv['SMtau'] + + # Convert a scattering measure (SM) to scattering timescale at given frequency. + tau_sc = 1.0 * (sm_tau / 292.0) ** 1.2 * dist * nu ** (-4.4) + + if not full_output: + return (float(dist) * u.kpc).to("pc"), tau_sc * u.s + else: + Dv["tau_sc"] = tau_sc + return Dv + +def dist_to_dm(gl, gb, dist, nu=1.0, model='ne2025', full_output=False): + """Convert a DM to a distance + + Args: + gl (float in deg): galactic longitude + gb (float in deg): galactic latitude + dist (float in pc/cm3): dispersion measure (pc cm^-3) + nu (float in GHz): Observing frequency + model (str): One of 'ne2001' or 'ne2025' + + Notes: + In this wrapper, 'ne2001' == 'ne2001p'. + + Returns: + dist (kpc), tau_sc (astropy.Quantity): distance (pc) and scattering time scale (s) + """ + if np.isclose(dist, 0): # WAR Catch spline failure + return 0.0 * u.pc, 0.0 * u.s + else: + if model.lower() == 'ne2025': + _Dk,Dv,_Du,_Dd = ne2025(gl, gb, dist, ndir=-1, dmd_only=False, classic=False) + else: + _Dk,Dv,_Du,_Dd = ne2001(gl, gb, dist, ndir=-1, dmd_only=False, classic=False) + + + dm, sm_tau = Dv['DM'], Dv['SMtau'] + + # Convert a scattering measure (SM) to scattering timescale at given frequency. + tau_sc = 1.0 * (sm_tau / 292.0) ** 1.2 * dist * nu ** (-4.4) + + if not full_output: + return float(dm) * u.pc / u.cm**3, tau_sc * u.s + else: + Dv["tau_sc"] = tau_sc + return Dv + + +def calculate_electron_density_xyz(x, y, z): + """Compute electron density at Galactocentric X, Y, Z coordinates + + x,y,z are Galactocentric Cartesian coordinates, measured in kpc (NOT pc!) + with the axes parallel to (l, b) = (90, 0), (180, 0), and (0, 90) degrees + + Args: + x (float): Galactocentric coordinates in kpc + y (float): Galactocentric coordinates in kpc + z (float): Galactocentric coordinates in kpc + + Returns: + ne_out (astropy.Quantity): Electron density in cm-3 + """ + density_list = density_2001(x, y, z) + + keys = ('ne1','ne2','nea','negc','nelism','necN','nevN', + 'F1', 'F2', 'Fa', 'Fgc', 'Flism', 'FcN', 'FvN', + 'whicharm', 'wlism', 'wldr', 'wlhb', 'wlsb', 'wloopI', + 'hitclump', 'hitvoid', 'wvoid') + + dd = dict(zip(keys, density_list)) + ne = dd['ne1'] + dd['ne2'] + dd['nea'] + dd['negc'] + dd['nelism'] + dd['necN'] + dd['nevN'] + ne = float(ne) + + return ne / u.cm**3 diff --git a/pygedm/pygedm.py b/pygedm/pygedm.py index dccff37..1446ea6 100644 --- a/pygedm/pygedm.py +++ b/pygedm/pygedm.py @@ -20,23 +20,25 @@ The Astrophysical Journal, Volume 888, Issue 2, id.105 """ -from . import ymw16_wrapper -from . import ne2001_wrapper -from . import yt2020 -from . import healpix_utils +import warnings + +import numpy as np from astropy import units as u -from astropy.coordinates import Angle, SkyCoord, Galactocentric +from astropy.coordinates import Angle, Galactocentric, SkyCoord from astropy.units import Quantity, Unit -import numpy as np -import warnings + +from tqdm import tqdm + +from . import healpix_utils, ne2001_wrapper, ymw16_wrapper, yt2020, ne2025_wrapper HAS_HEALPIX = healpix_utils.check_for_healpix_install() if HAS_HEALPIX: import healpy as hp -def _gl_gb_convert(gl, gb, unit='rad'): - """ Convert gl, gb astropy.Angle to floats/arrays with correct units + +def _gl_gb_convert(gl, gb, unit="rad"): + """Convert gl, gb astropy.Angle to floats/arrays with correct units If not an astropy quantity, returns value unchanged. """ @@ -48,7 +50,7 @@ def _gl_gb_convert(gl, gb, unit='rad'): def _unit_convert(q, unit): - """ Convert astropy.unit into float/array with given unit + """Convert astropy.unit into float/array with given unit If not an astropy quantity, returns value unchanged. """ @@ -57,8 +59,8 @@ def _unit_convert(q, unit): return q -def dm_to_dist(gl, gb, dm, dm_host=0, mode='gal', method='ymw16', nu=1.0): - """ Convert a DM to a distance +def dm_to_dist(gl, gb, dm, dm_host=0, mode="gal", method="ymw16", nu=1.0): + """Convert a DM to a distance Args: gl (float in deg or astropy.Angle): galactic longitude @@ -71,22 +73,27 @@ def dm_to_dist(gl, gb, dm, dm_host=0, mode='gal', method='ymw16', nu=1.0): Returns: dist (astropy.Quantity), tau_sc (astropy.Quantity): Distance (pc), scattering timescale at 1 GHz (s) """ - gl, gb = _gl_gb_convert(gl, gb, 'deg') - dm = _unit_convert(dm, 'pc cm^(-3)') - nu = _unit_convert(nu, 'GHz') - - if method.lower() == 'ymw16': + gl, gb = _gl_gb_convert(gl, gb, "deg") + dm = _unit_convert(dm, "pc cm^(-3)") + nu = _unit_convert(nu, "GHz") + + if method.lower() in ("ne2001", "ne2001p", "ne2025"): + if mode != "gal": + raise RuntimeError("NE2001 only supports Galactic (gal) mode.") + if method.lower() == 'ne2025': + return ne2025_wrapper.dm_to_dist(gl, gb, dm - dm_host, nu=nu, model='ne2025') + elif method.lower() == 'ne2001p': + return ne2025_wrapper.dm_to_dist(gl, gb, dm - dm_host, nu=nu, model='ne2001p') + else: + return ne2001_wrapper.dm_to_dist(gl, gb, dm - dm_host, nu=nu) + elif method.lower() == "ymw16": return ymw16_wrapper.dm_to_dist(gl, gb, dm, dm_host=0, mode=mode, nu=nu) - elif method.lower() == 'ne2001': - if mode != 'gal': - raise RuntimeError('NE2001 only supports Galactic (gal) mode.') - return ne2001_wrapper.dm_to_dist(gl, gb, dm - dm_host, nu=nu) else: - raise RuntimeError("Only ymw16 and ne2001 models supported.") + raise RuntimeError("Only ymw16, ne2001, ne2001p and ne2025 models supported.") -def dist_to_dm(gl, gb, dist, mode='gal', method='ymw16', nu=1.0): - """ Convert a distance to a DM +def dist_to_dm(gl, gb, dist, mode="gal", method="ymw16", nu=1.0): + """Convert a distance to a DM Args: gl (float in deg or astropy.Angle): galactic longitude @@ -99,33 +106,38 @@ def dist_to_dm(gl, gb, dist, mode='gal', method='ymw16', nu=1.0): Returns: dm (astropy.Quantity), tau_sc (astropy.Quantity): Dispersion measure (pc / cm3), scattering timescale at 1 GHz (s) """ - gl, gb = _gl_gb_convert(gl, gb, 'deg') + gl, gb = _gl_gb_convert(gl, gb, "deg") - if mode == 'igm': - dist = _unit_convert(dist, 'Mpc') + if mode == "igm": + dist = _unit_convert(dist, "Mpc") else: - dist = _unit_convert(dist, 'pc') - + dist = _unit_convert(dist, "pc") + # Catch NE2001 crash if dist too large - if _unit_convert(dist, 'pc') >= 100000: + if _unit_convert(dist, "pc") >= 100000: dist = 50000 warnings.warn("Distance too large, setting to 50 kpc.", UserWarning) - nu = _unit_convert(nu, 'GHz') + nu = _unit_convert(nu, "GHz") - if method.lower() == 'ymw16': + if method.lower() == "ymw16": return ymw16_wrapper.dist_to_dm(gl, gb, dist, mode=mode, nu=nu) - elif method.lower() == 'ne2001': - if mode != 'gal': - raise RuntimeError('NE2001 only supports Galactic (gal) mode.') + elif method.lower() in ("ne2001", "ne2001p", "ne2025"): + if mode != "gal": + raise RuntimeError("NE2001 only supports Galactic (gal) mode.") dist_kpc = dist / 1000.0 - return ne2001_wrapper.dist_to_dm(gl, gb, dist_kpc, nu=nu) + if method.lower() == 'ne2001': + return ne2001_wrapper.dist_to_dm(gl, gb, dist_kpc, nu=nu) + elif method.lower() == "ne2001p": + return ne2025_wrapper.dist_to_dm(gl, gb, dist_kpc, nu=nu, model='ne2001p') + else: + return ne2025_wrapper.dist_to_dm(gl, gb, dist_kpc, nu=nu, model='ne2025') else: - raise RuntimeError("Only ymw16 and ne2001 models supported.") + raise RuntimeError("Only ymw16, ne2001, ne2001p and ne2025 models supported.") -def calculate_electron_density_xyz(x, y, z, method='ymw16'): - """ Calculate electron density at a point with galactocentric coords (x, y, z) +def calculate_electron_density_xyz(x, y, z, method="ymw16"): + """Calculate electron density at a point with galactocentric coords (x, y, z) Args: x (float or Quantity): galactocentric X coordinate in pc @@ -135,20 +147,22 @@ def calculate_electron_density_xyz(x, y, z, method='ymw16'): Returns: N_e (astropy.quantity): electron density in cm^-3 """ - x = _unit_convert(x, 'pc') - y = _unit_convert(y, 'pc') - z = _unit_convert(z, 'pc') + x = _unit_convert(x, "pc") + y = _unit_convert(y, "pc") + z = _unit_convert(z, "pc") - if method.lower() == 'ymw16': + if method.lower() == "ymw16": return ymw16_wrapper.calculate_electron_density_xyz(x, y, z) - elif method.lower() == 'ne2001': - return ne2001_wrapper.calculate_electron_density_xyz(x/1e3, y/1e3, z/1e3) + elif method.lower() == "ne2001": + return ne2001_wrapper.calculate_electron_density_xyz(x / 1e3, y / 1e3, z / 1e3) + elif method.lower() == "ne2025": + return ne2025_wrapper.calculate_electron_density_xyz(x / 1e3, y / 1e3, z / 1e3) else: - raise RuntimeError("Only ymw16 and ne2001 models supported.") + raise RuntimeError("Only ymw16, ne2001, and ne2025 models supported.") -def calculate_electron_density_lbr(gl, gb, dist, method='ymw16'): - """ Calculate electron density at a point with Galactic coords (ga, gl) at given distance +def calculate_electron_density_lbr(gl, gb, dist, method="ymw16"): + """Calculate electron density at a point with Galactic coords (ga, gl) at given distance Args: gl (float, Angle, or Quantity): Galatic longitude in degrees (or astropy Angle) @@ -158,20 +172,22 @@ def calculate_electron_density_lbr(gl, gb, dist, method='ymw16'): Returns: N_e (astropy.Quantity): electron density in cm^-3 """ - gl, gb = _gl_gb_convert(gl, gb, 'deg') - dist = _unit_convert(dist, 'pc') + gl, gb = _gl_gb_convert(gl, gb, "deg") + dist = _unit_convert(dist, "pc") - if method.lower() == 'ymw16': + if method.lower() == "ymw16": return ymw16_wrapper.calculate_electron_density_lbr(gl, gb, dist) - elif method.lower() == 'ne2001': - x, y, z = convert_lbr_to_xyz(gl, gb, dist, method='ne2001') - return ne2001_wrapper.calculate_electron_density_xyz(x.to('kpc').value, y.to('kpc').value, z.to('kpc').value) + elif method.lower() == "ne2001": + x, y, z = convert_lbr_to_xyz(gl, gb, dist, method="ne2001") + return ne2001_wrapper.calculate_electron_density_xyz( + x.to("kpc").value, y.to("kpc").value, z.to("kpc").value + ) else: raise RuntimeError("Only ymw16 and ne2001 models supported.") -def convert_lbr_to_xyz(gl, gb, dist, method='ymw16'): - """ Convert Galactic (l,b,r) coords to Galactocentric (x,y,z) coords +def convert_lbr_to_xyz(gl, gb, dist, method="ymw16"): + """Convert Galactic (l,b,r) coords to Galactocentric (x,y,z) coords Args: gl (float, Angle, or Quantity): Galatic longitude in degrees (or astropy Angle) @@ -205,24 +221,26 @@ def convert_lbr_to_xyz(gl, gb, dist, method='ymw16'): (, , ) """ - gl, gb = _gl_gb_convert(gl, gb, 'deg') + gl, gb = _gl_gb_convert(gl, gb, "deg") gl = np.deg2rad(gl) gb = np.deg2rad(gb) - dist = _unit_convert(dist, 'pc') # Unit conversion and strip quantity - dist = dist * u.pc # Make sure distance is in units of pc - if method == 'astropy': - lbr = SkyCoord(gl, gb, distance=dist, frame='galactic', unit=('rad', 'rad', 'pc')) + dist = _unit_convert(dist, "pc") # Unit conversion and strip quantity + dist = dist * u.pc # Make sure distance is in units of pc + if method == "astropy": + lbr = SkyCoord( + gl, gb, distance=dist, frame="galactic", unit=("rad", "rad", "pc") + ) xyz = lbr.transform_to(Galactocentric()) return xyz.x, xyz.y, xyz.z - elif method.lower() == 'ymw16': + elif method.lower() == "ymw16": Rsun = 8300 * u.pc Zsun = 6.0 * u.pc x = dist * np.sin(gl) * np.cos(gb) y = Rsun - dist * np.cos(gl) * np.cos(gb) z = Zsun + dist * np.sin(gb) return x, y, z - elif method.lower() == 'ne2001': + elif method.lower() == "ne2001": Rsun = 8500 * u.pc x = dist * np.sin(gl) * np.cos(gb) y = Rsun - dist * np.cos(gl) * np.cos(gb) @@ -232,13 +250,13 @@ def convert_lbr_to_xyz(gl, gb, dist, method='ymw16'): raise RuntimeError("method must be astropy, ne2001, or ymw16") -def generate_healpix_dm_map(dist=1, nside=64, method='ymw16'): - """ Generate an all-sky healpix map for a given distance. +def generate_healpix_dm_map(dist=1, nside=64, method="ymw16"): + """Generate an all-sky healpix map for a given distance. Args: dist (float or Quantity): Distance to integrate EDM out to. 30 kpc will go to edge nside (int): The NSIDE parameter for the healpix map (power of 2, larger => higher resolution) - method (str): one of 'ymw16', 'ne2001', 'yt2020' or 'yt2020_analytic' + method (str): one of 'ymw16', 'ne2001', 'ne2025', 'yt2020' or 'yt2020_analytic' Notes: This method takes a fair amount of time to run -- tens of seconds for NSIDE=32. @@ -250,18 +268,18 @@ def generate_healpix_dm_map(dist=1, nside=64, method='ymw16'): if HAS_HEALPIX: pix_id = np.arange(hp.nside2npix(nside)) gl, gb = hp.pix2ang(nside, pix_id, lonlat=True) - d = np.zeros_like(pix_id, dtype='float32') - if 'yt2020' not in method: - for ii in pix_id: + d = np.zeros_like(pix_id, dtype="float32") + if "yt2020" not in method: + for ii in tqdm(pix_id): dm, tau_sc = dist_to_dm(gl[ii], gb[ii], dist, method=method) d[ii] = dm.value - elif 'yt2020' in method: + elif "yt2020" in method: # gl is in (0, 360), wrap to (-180, 180) as reqd for yt2020 function # Then call yt2020.compute_halo_dm function directly gl_u = np.copy(gl) gl_u[gl >= 180] -= 360 - if method == 'yt2020': + if method == "yt2020": for ii in pix_id: dm_halo = yt2020.calculate_halo_dm(gl_u[ii], gb[ii]) d[ii] = dm_halo.value @@ -274,8 +292,8 @@ def generate_healpix_dm_map(dist=1, nside=64, method='ymw16'): raise RuntimeError("Healpy installation not found.") -def calculate_halo_dm(gl, gb, method='yt2020', component='both'): - """ Compute halo DM +def calculate_halo_dm(gl, gb, method="yt2020", component="both"): + """Compute halo DM Args: gl (float, Angle, or Quantity): Galatic latitude in degrees (or astropy Angle) @@ -287,13 +305,13 @@ def calculate_halo_dm(gl, gb, method='yt2020', component='both'): DM (float): Dispersion measure in (pc/cm3) """ - gl, gb = _gl_gb_convert(gl, gb, 'deg') - gl = Angle(gl, unit='deg').wrap_at('180d').value ## Ensure in correct range + gl, gb = _gl_gb_convert(gl, gb, "deg") + gl = Angle(gl, unit="deg").wrap_at("180d").value ## Ensure in correct range - if method == 'yt2020': + if method == "yt2020": return yt2020.calculate_halo_dm(gl, gb, component) - elif method == 'yt2020_analytic': - if component != 'both': + elif method == "yt2020_analytic": + if component != "both": raise RuntimeError("Analytic formula requires component='both'") return yt2020.calculate_halo_dm_analytic(gl, gb) else: diff --git a/pygedm/ymw16_wrapper.py b/pygedm/ymw16_wrapper.py index 34eabfa..ebba697 100644 --- a/pygedm/ymw16_wrapper.py +++ b/pygedm/ymw16_wrapper.py @@ -8,23 +8,26 @@ ApJ, 835, 29. """ +import os + import ymw16 from astropy import units as u from astropy.coordinates import Angle from astropy.units import Quantity, Unit -import os -from pkg_resources import resource_filename, get_distribution, DistributionNotFound -DATAPATH = os.path.dirname(resource_filename("pygedm", "spiral.txt")) +try: + from importlib.resources import files +except ImportError: + from importlib_resources import files +DATAPATH = str(files("pygedm")) -MODE_IDS = {'gal': 1, - 'mc': 0, - 'igm': -1} +MODE_IDS = {"gal": 1, "mc": 0, "igm": -1} -def dm_to_dist(gl, gb, dm, dm_host=0, mode='gal', nu=1.0): - """ Convert a DM to a distance + +def dm_to_dist(gl, gb, dm, dm_host=0, mode="gal", nu=1.0): + """Convert a DM to a distance Args: gl (float in deg or astropy.Angle): galactic longitude @@ -38,21 +41,21 @@ def dm_to_dist(gl, gb, dm, dm_host=0, mode='gal', nu=1.0): """ mode_id = MODE_IDS.get(mode.lower().strip()) - ndir, vbs, txt = 1, 0, '' + ndir, vbs, txt = 1, 0, "" r = ymw16.dmdtau(gl, gb, dm, dm_host, ndir, mode_id, vbs, DATAPATH, txt) - if mode == 'igm': - r['dist'] *= u.Mpc - r['tau_sc'] = r['tau_FRB'] + if mode == "igm": + r["dist"] *= u.Mpc + r["tau_sc"] = r["tau_FRB"] else: - r['dist'] *= u.pc - dist = r['dist'] - tau_sc = r['tau_sc'] * nu ** -4 * u.s + r["dist"] *= u.pc + dist = r["dist"] + tau_sc = r["tau_sc"] * nu**-4 * u.s return (dist, tau_sc) -def dist_to_dm(gl, gb, dist, mode='gal', nu=1.0): - """ Convert a distance to a DM +def dist_to_dm(gl, gb, dist, mode="gal", nu=1.0): + """Convert a distance to a DM Args: gl (float in deg or astropy.Angle): galactic longitude @@ -65,19 +68,19 @@ def dist_to_dm(gl, gb, dist, mode='gal', nu=1.0): dm (astropy.Quantity), tau_sc (astropy.Quantity): dispersion measure (pc/cm3) and scattering time scale (s) """ mode_id = MODE_IDS.get(mode.lower().strip()) - ndir, dm_host, vbs, txt = 2, 0, 0, '' + ndir, dm_host, vbs, txt = 2, 0, 0, "" r = ymw16.dmdtau(gl, gb, dist, dm_host, ndir, mode_id, vbs, DATAPATH, txt) - if mode == 'igm': - r['DM'] = r['DM_IGM'] - r['tau_sc'] = r['tau_FRB'] - dm = r['DM'] * u.pc / u.cm**3 - tau_sc = r['tau_sc'] * nu ** -4 * u.s + if mode == "igm": + r["DM"] = r["DM_IGM"] + r["tau_sc"] = r["tau_FRB"] + dm = r["DM"] * u.pc / u.cm**3 + tau_sc = r["tau_sc"] * nu**-4 * u.s return dm, tau_sc def calculate_electron_density_xyz(x, y, z): - """ Calculate electron density at a point with galactocentric coords (x, y, z) + """Calculate electron density at a point with galactocentric coords (x, y, z) Args: x (float or Quantity): galactocentric X coordinate in pc @@ -87,13 +90,13 @@ def calculate_electron_density_xyz(x, y, z): Returns: N_e (astropy.quantity): electron density in cm^-3 """ - ndir, vbs, txt = 1, 0, '' + ndir, vbs, txt = 1, 0, "" ne = ymw16.ne_crd(x, y, z, 0, 0, 0, ndir, vbs, DATAPATH, txt) return ne / u.cm**3 def calculate_electron_density_lbr(gl, gb, dist): - """ Calculate electron density at a point with Galactic coords (ga, gl) + """Calculate electron density at a point with Galactic coords (ga, gl) at a given distance in pc Args: @@ -104,6 +107,6 @@ def calculate_electron_density_lbr(gl, gb, dist): Returns: N_e (astropy.Quantity): electron density in cm^-3 """ - ndir, vbs, txt = 2, 0, '' + ndir, vbs, txt = 2, 0, "" ne = ymw16.ne_crd(0, 0, 0, gl, gb, dist, ndir, vbs, DATAPATH, txt) - return ne / u.cm**3 + return ne / u.cm**3 diff --git a/pygedm/yt2020.py b/pygedm/yt2020.py index 45fd828..81b9275 100644 --- a/pygedm/yt2020.py +++ b/pygedm/yt2020.py @@ -9,44 +9,62 @@ Notes: Adapted from S. Yamasaki's ``DM_halo_yt2020_numerical.py`` command-line python code """ + import numpy as np -from scipy import integrate -from numpy import cos, sin, sqrt, fabs from astropy import units as u +from numpy import cos, fabs, sin, sqrt +from scipy import integrate ### numerical constants ### -U_G = 6.67 # G = 6.67259e-8 cm^3/g/s^2 -U_Msun = 1.99 # 1 M_sun = 1.99e+33 g -U_pc = 3.09 # 1 pc = 3.086e+18 cm -U_mp = 1.67 # m_p = 1.6726231e-24 g -U_eV = 1.60 # 1 eV = 1.602177e-12 erg -c_scaling = U_Msun/U_pc**3*pow(10,-18) # scaling factor [10^{12} M_sun/kpc^3] -> [g/cm^3] -m_p = 1.6726231*pow(10,-24) # proton mass [g] +U_G = 6.67 # G = 6.67259e-8 cm^3/g/s^2 +U_Msun = 1.99 # 1 M_sun = 1.99e+33 g +U_pc = 3.09 # 1 pc = 3.086e+18 cm +U_mp = 1.67 # m_p = 1.6726231e-24 g +U_eV = 1.60 # 1 eV = 1.602177e-12 erg +c_scaling = ( + U_Msun / U_pc**3 * pow(10, -18) +) # scaling factor [10^{12} M_sun/kpc^3] -> [g/cm^3] +m_p = 1.6726231 * pow(10, -24) # proton mass [g] ### key parameters ### -T = 0.3 # halo temperature [keV] -M_b = 0.036 # total halo gas mass [10^{12} Z^{-1} M_sun] (such that M_b = 0.12 for Z = Z_halo) -Z_halo = 0.3 # halo metallicity [Z_sun] -xi_H = 28.0/34.0 # abundance ratio of hydrogen to electron -mu_e = 1.18 # mean molecular mass per electron -mu = 0.62 # mean particle mass -D_sun = 8.5 # Galactocentoric distance of the Sun [kpc] -c_NFW = 12 # NFW concentration -r_vir = 260 # virial radius of MW [kpc] -M_vir = 1.0 # virial mass of MW [10^{12} M_sun] -r_s = r_vir/c_NFW # NFW scale radius [kpc] -f_NFW = np.log(1.0+c_NFW)-c_NFW/(1.0+c_NFW) # f(x)=np.log(1.0+x)-x/(1.0+x) and x = c_NFW -rho_s = M_vir/4.0/np.pi/r_s**(3.0)/f_NFW # NFW scale density [10^{12} M_sun/kpc^3] +T = 0.3 # halo temperature [keV] +M_b = 0.036 # total halo gas mass [10^{12} Z^{-1} M_sun] (such that M_b = 0.12 for Z = Z_halo) +Z_halo = 0.3 # halo metallicity [Z_sun] +xi_H = 28.0 / 34.0 # abundance ratio of hydrogen to electron +mu_e = 1.18 # mean molecular mass per electron +mu = 0.62 # mean particle mass +D_sun = 8.5 # Galactocentoric distance of the Sun [kpc] +c_NFW = 12 # NFW concentration +r_vir = 260 # virial radius of MW [kpc] +M_vir = 1.0 # virial mass of MW [10^{12} M_sun] +r_s = r_vir / c_NFW # NFW scale radius [kpc] +f_NFW = np.log(1.0 + c_NFW) - c_NFW / ( + 1.0 + c_NFW +) # f(x)=np.log(1.0+x)-x/(1.0+x) and x = c_NFW +rho_s = ( + M_vir / 4.0 / np.pi / r_s ** (3.0) / f_NFW +) # NFW scale density [10^{12} M_sun/kpc^3] ### calculate n_0^{sphe} (eq.[3]) by normalization (eq.[4]) ### -r_array = np.linspace(r_s/100,r_vir,100000) # r = 0-r_vir [kpc] -Upsilon = 4*np.pi*10*U_G*U_Msun/U_pc*U_mp/U_eV*r_s**2*rho_s*mu/T # Upsilon appearing in right-hand side of eq.(3) +r_array = np.linspace(r_s / 100, r_vir, 100000) # r = 0-r_vir [kpc] +Upsilon = ( + 4 * np.pi * 10 * U_G * U_Msun / U_pc * U_mp / U_eV * r_s**2 * rho_s * mu / T +) # Upsilon appearing in right-hand side of eq.(3) # integrated function excluding n_0^{sphe} (left-hand side of eq.[4]) -y_sphe = 4.0*np.pi*r_array**2*mu_e*m_p/c_scaling*np.exp(-Upsilon*(1-np.log(1+r_array/r_s)*r_s/r_array)) +y_sphe = ( + 4.0 + * np.pi + * r_array**2 + * mu_e + * m_p + / c_scaling + * np.exp(-Upsilon * (1 - np.log(1 + r_array / r_s) * r_s / r_array)) +) + I_sphe = integrate.simpson(y_sphe,r_array) # integration n_0_sphe = M_b/Z_halo/I_sphe # central electron number density for spherical comp. [cm^{-3}] @@ -56,8 +74,10 @@ ### key functions ### -def s_max(l, b): # maximum integration limit corresponsing to r = r_vir as function of sky coordinates [kpc] - """ Compute integration limit s_max for given sky coordinates +def s_max( + l, b +): # maximum integration limit corresponsing to r = r_vir as function of sky coordinates [kpc] + """Compute integration limit s_max for given sky coordinates Args: l (float): Galactic longitude, in radians (-pi to +pi) @@ -66,11 +86,13 @@ def s_max(l, b): # maximum integration limit corresponsing to r = r_vir as funct Returns: s_max (float), maximum integration limit corresponsing to r = r_vir """ - return D_sun*cos(b)*cos(l)+sqrt(D_sun**2*(cos(b)**2*cos(l)**2-1)+(1.0*r_vir)**2) + return D_sun * cos(b) * cos(l) + sqrt( + D_sun**2 * (cos(b) ** 2 * cos(l) ** 2 - 1) + (1.0 * r_vir) ** 2 + ) def ne_sphe(l, b, s): - """ Compute electron density for spherical component for (l, b) at distance s + """Compute electron density for spherical component for (l, b) at distance s Args: l (float): Galactic longitude, in radians (-pi to +pi) @@ -80,15 +102,19 @@ def ne_sphe(l, b, s): Returns: ne (float): electron density in cm^{-3} """ - R = sqrt(D_sun**2+(s*cos(b))**2-2*D_sun*s*cos(b)*cos(l)) - z = s*sin(b) - r = sqrt(R**2+z**2) - Upsilon = 4*np.pi*10*U_G*U_Msun/U_pc*U_mp/U_eV*r_s**2*rho_s*mu/T - return n_0_sphe*np.exp(-Upsilon*(1-np.log(1+r/r_s)*r_s/r)) # [cm^{-3}] + R = sqrt(D_sun**2 + (s * cos(b)) ** 2 - 2 * D_sun * s * cos(b) * cos(l)) + z = s * sin(b) + r = sqrt(R**2 + z**2) + Upsilon = ( + 4 * np.pi * 10 * U_G * U_Msun / U_pc * U_mp / U_eV * r_s**2 * rho_s * mu / T + ) + return n_0_sphe * np.exp( + -Upsilon * (1 - np.log(1 + r / r_s) * r_s / r) + ) # [cm^{-3}] def ne_disk(l, b, s): - """ Compute electron density for spherical component for (l, b) at distance s + """Compute electron density for spherical component for (l, b) at distance s Args: l (float): Galactic longitude, in radians (-pi to +pi) @@ -98,20 +124,21 @@ def ne_disk(l, b, s): Returns: ne (float): electron density in cm^{-3} """ - R = sqrt(D_sun**2+(s*cos(b))**2-2*D_sun*s*cos(b)*cos(l)) - z = fabs(s*sin(b)) - n_0_disk = 7.4*pow(10,-3) # [Z^{-1} cm^{-3}] - z_0 = 2.4 # [kpc] - R_0 = 4.9 # [kpc] - return n_0_disk/np.exp(R/R_0+z/z_0)/Z_halo # [cm^{-3}] + R = sqrt(D_sun**2 + (s * cos(b)) ** 2 - 2 * D_sun * s * cos(b) * cos(l)) + z = fabs(s * sin(b)) + n_0_disk = 7.4 * pow(10, -3) # [Z^{-1} cm^{-3}] + z_0 = 2.4 # [kpc] + R_0 = 4.9 # [kpc] + return n_0_disk / np.exp(R / R_0 + z / z_0) / Z_halo # [cm^{-3}] ### vectorize density profile functions ### vfunc_sphe = np.vectorize(ne_sphe) vfunc_disk = np.vectorize(ne_disk) -def calculate_halo_dm(l, b, component='both'): - """ Compute halo DM + +def calculate_halo_dm(l, b, component="both"): + """Compute halo DM Args: l (float): Galactic longitude, in degrees (-180 to +180) @@ -133,23 +160,23 @@ def calculate_halo_dm(l, b, component='both'): y_DM_sphe = vfunc_sphe(l_FRB, b_FRB, xx) y_DM_disk = vfunc_disk(l_FRB, b_FRB, xx) - if component == 'spherical' or component == 'both': - DM_sphe = 1000 * integrate.simps(y_DM_sphe, xx) * u.pc / u.cm**3 # [pc/cm^3] - if component == 'disk' or component == 'both': - DM_disk = 1000 * integrate.simps(y_DM_disk, xx) * u.pc / u.cm**3 # [pc/cm^3] + if component == "spherical" or component == "both": + DM_sphe = 1000 * integrate.simpson(y_DM_sphe, xx) * u.pc / u.cm**3 # [pc/cm^3] + if component == "disk" or component == "both": + DM_disk = 1000 * integrate.simpson(y_DM_disk, xx) * u.pc / u.cm**3 # [pc/cm^3] - if component == 'both': + if component == "both": return DM_sphe + DM_disk - elif component == 'spherical': + elif component == "spherical": return DM_sphe - elif component == 'disk': + elif component == "disk": return DM_disk else: raise RuntimeError("Component must be 'both', 'spherical', or 'disk' ") def calculate_halo_dm_analytic(l, b): - """ Calculate the DM contribution of the Galactic halo. + """Calculate the DM contribution of the Galactic halo. Uses an analytical formula for speed. Useful for all-sky mapping. @@ -157,18 +184,20 @@ def calculate_halo_dm_analytic(l, b): l (float): Galactic longitude, in degrees (-180 to +180) b (float): Galactic latitude, in degrees (-90 to 90) """ - YT19_coeffs = np.array([[250.12, -871.06, 1877.5, -2553.0, 2181.3, -1127.5, 321.72, -38.905], - [-154.82, 783.43, -1593.9, 1727.6, -1046.5, 332.09, -42.815, 0], - [-116.72, -76.815, 428.49, -419.00, 174.60, -27.610, 0, 0], - [216.67, -193.30, 12.234, 32.145, -8.3602, 0, 0, 0], - [-129.95, 103.80, -22.800, 0.44171, 0, 0, 0, 0], - [39.652, -21.398, 2.7694, 0, 0, 0, 0, 0], - [-6.1926, 1.6162, 0, 0, 0, 0, 0, 0], - [0.39346, 0, 0, 0, 0, 0, 0, 0]]) + YT19_coeffs = np.array([ + [250.12, -871.06, 1877.5, -2553.0, 2181.3, -1127.5, 321.72, -38.905], + [-154.82, 783.43, -1593.9, 1727.6, -1046.5, 332.09, -42.815, 0], + [-116.72, -76.815, 428.49, -419.00, 174.60, -27.610, 0, 0], + [216.67, -193.30, 12.234, 32.145, -8.3602, 0, 0, 0], + [-129.95, 103.80, -22.800, 0.44171, 0, 0, 0, 0], + [39.652, -21.398, 2.7694, 0, 0, 0, 0, 0], + [-6.1926, 1.6162, 0, 0, 0, 0, 0, 0], + [0.39346, 0, 0, 0, 0, 0, 0, 0], + ]) l_abs = np.abs(np.deg2rad(l)) b_abs = np.abs(np.deg2rad(b)) DM_halo = 0 for ii in range(8): for jj in range(8): - DM_halo += YT19_coeffs[ii, jj] * l_abs ** ii * b_abs ** jj - return DM_halo * u.pc / u.cm**3 \ No newline at end of file + DM_halo += YT19_coeffs[ii, jj] * l_abs**ii * b_abs**jj + return DM_halo * u.pc / u.cm**3 diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..ae7ef27 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,93 @@ +[build-system] +requires = ["setuptools>=64.0", "wheel", "pybind11>=2.2"] +build-backend = "setuptools.build_meta" + +[project] +name = "pygedm" +version = "4.0.0" +description = "Python/C++ version of NE2001, YMW16, and YT2020 electron density models" +readme = "README.md" +requires-python = ">=3.8" +license = {text = "MIT"} +authors = [ + {name = "D. C. Price", email = "dan@thetelegraphic.com"} +] +keywords = ["astronomy", "radio", "pulsars", "dispersion measure"] +classifiers = [ + "Development Status :: 4 - Beta", + "Intended Audience :: Science/Research", + "License :: OSI Approved :: MIT License", + "Programming Language :: Python :: 3", + "Programming Language :: Python :: 3.8", + "Programming Language :: Python :: 3.9", + "Programming Language :: Python :: 3.10", + "Programming Language :: Python :: 3.11", + "Programming Language :: Python :: 3.12", + "Programming Language :: Python :: 3.13", + "Topic :: Scientific/Engineering :: Astronomy", +] +dependencies = [ + "pybind11>=2.2", + "astropy", + "numpy", + "scipy", + "tqdm", + "importlib-metadata>=1.0; python_version < '3.8'", + "importlib-resources>=1.3; python_version < '3.9'" +] + +[project.optional-dependencies] +test = [ + "pytest>=6.0", + "pytest-cov", + "healpy", +] +dev = [ + "pytest>=6.0", + "pytest-cov", + "healpy", + "build", + "twine", +] + +[project.urls] +Homepage = "https://github.com/frbs/pygedm" +Repository = "https://github.com/frbs/pygedm" +"Bug Tracker" = "https://github.com/frbs/pygedm/issues" +Documentation = "https://pygedm.readthedocs.io" + +[tool.setuptools] +packages = ["pygedm"] +zip-safe = false +include-package-data = true + +[tool.setuptools.package-data] +pygedm = [ + "spiral.txt", + "ymw16par.txt", + "gal01.inp", + "ne_arms_log_mod.inp", + "ne_gc.inp", + "nelism.inp", + "neclumpN.NE2001.dat", + "nevoidN.NE2001.dat", +] + +[tool.pytest.ini_options] +addopts = ["--verbose", "--cov=pygedm"] +testpaths = ["tests"] +minversion = "6.0" + +[tool.coverage.run] +source = ["pygedm"] +omit = ["*/tests/*", "*/test_*.py"] + +[tool.coverage.report] +exclude_lines = [ + "pragma: no cover", + "def __repr__", + "raise AssertionError", + "raise NotImplementedError", + "if __name__ == .__main__.:", + "if TYPE_CHECKING:", +] diff --git a/references.bib b/references.bib index 499813d..66f1ad3 100644 --- a/references.bib +++ b/references.bib @@ -70,3 +70,19 @@ @ARTICLE{Yamasaki:2020 adsnote = {Provided by the SAO/NASA Astrophysics Data System} } +@ARTICLE{Ocker:2026, + author = {{Ocker}, S.~K. and {Cordes}, J.~M.}, + title = "{NE2025: An Updated Electron Density Model for the Galactic Interstellar Medium}", + journal = {arXiv e-prints}, + keywords = {Astrophysics of Galaxies, High Energy Astrophysical Phenomena}, + year = 2026, + month = feb, + eid = {arXiv:2602.11838}, + pages = {arXiv:2602.11838}, + doi = {10.48550/arXiv.2602.11838}, +archivePrefix = {arXiv}, + eprint = {2602.11838}, + primaryClass = {astro-ph.GA}, + adsurl = {https://ui.adsabs.harvard.edu/abs/2026arXiv260211838O}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} +} diff --git a/requirements.txt b/requirements.txt index 6b68ace..fa7cce4 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,5 +1,7 @@ -pybind11 +pybind11>=2.2 astropy numpy scipy -healpy \ No newline at end of file +healpy +importlib-metadata>=1.0; python_version < "3.8" +importlib-resources>=1.3; python_version < "3.9" diff --git a/setup.cfg b/setup.cfg index 9003394..3b478a1 100644 --- a/setup.cfg +++ b/setup.cfg @@ -2,7 +2,7 @@ universal=1 [metadata] -description-file=README.md +description_file=README.md [aliases] test=pytest diff --git a/setup.py b/setup.py index bd617cf..5590d3d 100644 --- a/setup.py +++ b/setup.py @@ -1,100 +1,99 @@ -# To increment version -# Check you have ~/.pypirc filled in -# git tag x.y.z -# git push --tags -# python setup.py sdist bdist_wheel -# twine upload dist/* - -from setuptools import setup, Extension -from setuptools.command.build_ext import build_ext -import sys +# Modern build configuration for pygedm with C++ extensions +# All package metadata is now in pyproject.toml following PEP 621 +# +# To build and publish: +# python -m build +# twine upload dist/* + import os -import glob -import setuptools +import sys +import tempfile +from setuptools import Extension, setup +from setuptools.command.build_ext import build_ext -__version__ = '3.3.0' __here__ = os.path.abspath(os.path.dirname(__file__)) + class get_pybind_include(object): """Helper class to determine the pybind11 include path The purpose of this class is to postpone importing pybind11 until it is actually installed, so that the ``get_include()`` - method can be invoked. """ + method can be invoked.""" def __init__(self, user=False): self.user = user def __str__(self): import pybind11 + return pybind11.get_include(self.user) +prefix = os.environ.get("CONDA_PREFIX") + ext_modules = [ Extension( - 'ymw16', + "ymw16", sources=[ - 'ymw16_src/main.cpp', - 'ymw16_src/dora.cpp', - 'ymw16_src/fermibubble.cpp', - 'ymw16_src/frb_d.cpp', - 'ymw16_src/galcen.cpp', - 'ymw16_src/gum.cpp', - 'ymw16_src/lmc.cpp', - 'ymw16_src/localbubble.cpp', - 'ymw16_src/ne_crd.cpp', - 'ymw16_src/nps.cpp', - 'ymw16_src/smc.cpp', - 'ymw16_src/spiral.cpp', - 'ymw16_src/thick.cpp', - 'ymw16_src/thin.cpp', - 'ymw16_src/ymw16par.cpp', - 'ymw16_src/dmdtau2.cpp', + "ymw16_src/main.cpp", + "ymw16_src/dora.cpp", + "ymw16_src/fermibubble.cpp", + "ymw16_src/frb_d.cpp", + "ymw16_src/galcen.cpp", + "ymw16_src/gum.cpp", + "ymw16_src/lmc.cpp", + "ymw16_src/localbubble.cpp", + "ymw16_src/ne_crd.cpp", + "ymw16_src/nps.cpp", + "ymw16_src/smc.cpp", + "ymw16_src/spiral.cpp", + "ymw16_src/thick.cpp", + "ymw16_src/thin.cpp", + "ymw16_src/ymw16par.cpp", + "ymw16_src/dmdtau2.cpp", ], include_dirs=[ # Path to pybind11 headers get_pybind_include(), get_pybind_include(user=True), - os.path.join(__here__, 'ymw16_src'), + os.path.join(__here__, "ymw16_src"), ], - - extra_link_args=['-lm'], - language='c++' + extra_link_args=["-lm"], + language="c++", ), Extension( - 'ne21c', + "ne21c", sources=[ - 'ne21c/main.cpp', + "ne21c/main.cpp", ], include_dirs=[ # Path to pybind11 headers get_pybind_include(), get_pybind_include(user=True), - os.path.join(__here__, 'ne21c'), + os.path.join(__here__, "ne21c"), + f"{prefix}/include" ], - extra_compile_args=['-Wno-write-strings'], - extra_link_args=['-lm', '-lf2c'], - language='c++', + library_dirs=[f"{prefix}/lib"], + extra_compile_args=["-Wno-write-strings"], + extra_link_args=["-lm", "-lf2c"], + language="c++", ), ] -ymw16_data_files = ['spiral.txt', 'ymw16par.txt'] -ne2001_data_files = ['gal01.inp', 'ne_arms_log_mod.inp', 'ne_gc.inp', - 'nelism.inp', 'neclumpN.NE2001.dat', 'nevoidN.NE2001.dat'] -data_files = ymw16_data_files + ne2001_data_files - -# As of Python 3.6, CCompiler has a `has_flag` method. -# cf http://bugs.python.org/issue26689 +ymw16_data_files = ["spiral.txt", "ymw16par.txt"] +# C++ extensions def has_flag(compiler, flagname): """Return a boolean indicating whether a flag name is supported on the specified compiler. """ import tempfile - with tempfile.NamedTemporaryFile('w', suffix='.cpp') as f: - f.write('int main (int argc, char **argv) { return 0; }') + + with tempfile.NamedTemporaryFile("w", suffix=".cpp") as f: + f.write("int main (int argc, char **argv) { return 0; }") try: compiler.compile([f.name], extra_postargs=[flagname]) - except setuptools.distutils.errors.CompileError: + C++ extensionsr: return False return True @@ -102,64 +101,46 @@ def has_flag(compiler, flagname): def cpp_flag(compiler): """Return the -std=c++[11/14] compiler flag. - The c++14 is prefered over c++11 (when it is available). + The c++14 is prefered over c++11 (when it is availa, delete=False) as f: + f.write("int main (int argc, char **argv) { return 0; }") + fname = f.name + try: + compiler.compile([fname], extra_postargs=[flagname]) + return True + except Exception: + return False + finally: + try: + os.unlink(fname)4/17] compiler flag. + + C++17 is preferred, then C++14, then C++11. """ - if has_flag(compiler, '-std=c++14'): - return '-std=c++14' - elif has_flag(compiler, '-std=c++11'): - return '-std=c++11' - else: - raise RuntimeError('Unsupported compiler -- at least C++11 support ' - 'is needed!') - - -class BuildExt(build_ext): - """A custom build extension for adding compiler-specific options.""" - c_opts = { - 'msvc': ['/EHsc'], - 'unix': [], + for flag in ["-std=c++17", "-std=c++14", "-std=c++11"]: + if has_flag(compiler, flag): + return flag + raise RuntimeError("Unsupported compiler -- at least C++11 support is needed!"ts = { + "msvc": ["/EHsc"], + "unix": [], } - if sys.platform == 'darwin': - c_opts['unix'] += ['-stdlib=libc++', '-mmacosx-version-min=10.7'] + if sys.platform == "darwin": + c_opts["unix"] += ["-stdlib=libc++", "-mmacosx-version-min=10.7"] def build_extensions(self): ct = self.compiler.compiler_type opts = self.c_opts.get(ct, []) - if ct == 'unix': - opts.append('-DVERSION_INFO="%s"' % self.distribution.get_version()) - opts.append(cpp_flag(self.compiler)) - # This seems to break libf2c :/ - #if has_flag(self.compiler, '-fvisibility=hidden'): - # opts.append('-fvisibility=hidden') - elif ct == 'msvc': + if ct == "unix": + elif ct == "msvc": opts.append('/DVERSION_INFO=\\"%s\\"' % self.distribution.get_version()) for ext in self.extensions: ext.extra_compile_args += opts build_ext.build_extensions(self) -with open("README.md", "r") as fh: - long_description = fh.read() - +# Run setup with configuration from pyproject.toml setup( - name='pygedm', - version=__version__, - author='D. C. Price', - author_email='dancpr@berkeley.edu', - description='Python/C++ version of NE2001, YMW16, and YT2020 electron density models', - long_description=long_description, - long_description_content_type='text/markdown', - url='https://github.com/frbs/pygedm', - download_url='https://github.com/telegraphic/pygedm/archive/%s.tar.gz' % __version__, - python_requires='>=3.6', - install_requires=['pybind11>=2.2', 'astropy'], - tests_require= ['pytest', 'astropy', 'numpy', 'healpy'], - setup_requires= ['pytest-runner', 'pytest-cov', 'pybind11>=2.2'], - ext_modules=ext_modules, - packages=['pygedm'], - package_data={'pygedm': data_files}, + ext_modules=ext_modules"pygedm": data_files}, include_package_data=True, zip_safe=False, - cmdclass={'build_ext': BuildExt}, + cmdclass={"build_ext": BuildExt}, ) diff --git a/tests/test_basic_ne2001.py b/tests/test_basic_ne2001.py index 7fdd91a..00eccf2 100644 --- a/tests/test_basic_ne2001.py +++ b/tests/test_basic_ne2001.py @@ -1,125 +1,144 @@ -import pygedm -import numpy as np -from astropy.coordinates import Angle -from astropy.units import Unit, Quantity import astropy.units as u +import numpy as np import pytest +from astropy.coordinates import Angle +from astropy.units import Quantity, Unit + +import pygedm def test_dm_to_dist(): - """ Test that astropy units / angles work with dm_to_dist """ - a = pygedm.dm_to_dist(204, -6.5, 200, method='ne2001') - b = pygedm.dm_to_dist(Angle(204, unit='degree'), Angle(-6.5, unit='degree'), 200, method='ne2001') - c = pygedm.dm_to_dist(204, -6.5, 200 * Unit('pc cm^-3'), method='ne2001') + """Test that astropy units / angles work with dm_to_dist""" + a = pygedm.dm_to_dist(204, -6.5, 200, method="ne2001") + b = pygedm.dm_to_dist( + Angle(204, unit="degree"), Angle(-6.5, unit="degree"), 200, method="ne2001" + ) + c = pygedm.dm_to_dist(204, -6.5, 200 * Unit("pc cm^-3"), method="ne2001") assert a[0] == b[0] == c[0] assert a[1] == b[1] == c[1] def test_dist_to_dm(): - """ Test that astropy units / angles work with dist_to_dm """ - a = pygedm.dist_to_dm(204, -6.5, 200, method='ne2001') - b = pygedm.dist_to_dm(Angle(204, unit='degree'), Angle(-6.5, unit='degree'), 200, method='ne2001') - c = pygedm.dist_to_dm(204, -6.5, 200 * Unit('pc'), method='ne2001') + """Test that astropy units / angles work with dist_to_dm""" + a = pygedm.dist_to_dm(204, -6.5, 200, method="ne2001") + b = pygedm.dist_to_dm( + Angle(204, unit="degree"), Angle(-6.5, unit="degree"), 200, method="ne2001" + ) + c = pygedm.dist_to_dm(204, -6.5, 200 * Unit("pc"), method="ne2001") assert a[0] == b[0] == c[0] assert a[1] == b[1] == c[1] def test_calculate_electron_density_xyz(): - pc = Unit('pc') - a = pygedm.calculate_electron_density_xyz(1, 2, 3, method='ne2001') - b = pygedm.calculate_electron_density_xyz(1 * pc, 2, 3, method='ne2001') - c = pygedm.calculate_electron_density_xyz(1, 2 * pc, 3, method='ne2001') - d = pygedm.calculate_electron_density_xyz(1, 2, 3 * pc, method='ne2001') + pc = Unit("pc") + a = pygedm.calculate_electron_density_xyz(1, 2, 3, method="ne2001") + b = pygedm.calculate_electron_density_xyz(1 * pc, 2, 3, method="ne2001") + c = pygedm.calculate_electron_density_xyz(1, 2 * pc, 3, method="ne2001") + d = pygedm.calculate_electron_density_xyz(1, 2, 3 * pc, method="ne2001") assert a == b == c == d def test_calculate_electron_density_lbr(): - ed_gc = pygedm.calculate_electron_density_xyz(0, 0, 0, method='ne2001') - ed_gc_lbr = pygedm.calculate_electron_density_lbr(0, 0, 8500, method='ne2001') + ed_gc = pygedm.calculate_electron_density_xyz(0, 0, 0, method="ne2001") + ed_gc_lbr = pygedm.calculate_electron_density_lbr(0, 0, 8500, method="ne2001") assert ed_gc == ed_gc_lbr def test_basic(): - """ Basic tests of YMW16 model + """Basic tests of YMW16 model Note: tested against online NE2001 interface https://www.nrl.navy.mil/rsd/RORF/ne2001/ """ # No access to actual model via web interface - #a = pygedm.calculate_electron_density_xyz(1, 2, 3) - #assert np.isclose(a.value, 5.220655, atol=0.0001) + # a = pygedm.calculate_electron_density_xyz(1, 2, 3) + # assert np.isclose(a.value, 5.220655, atol=0.0001) # FRB180301 value - dm, tau = pygedm.dist_to_dm(204, -6.5, 25*u.kpc, method='ne2001') + dm, tau = pygedm.dist_to_dm(204, -6.5, 25 * u.kpc, method="ne2001") assert np.isclose(dm.value, 150.80, atol=0.01) # Loop through distances and check round trip - for dist in (10.*u.pc, 100.*u.pc, 1000.*u.pc): - dm, tau = pygedm.dist_to_dm(0, 0, dist, method='ne2001') - dist_out, tau = pygedm.dm_to_dist(0, 0, dm, method='ne2001') + for dist in (10.0 * u.pc, 100.0 * u.pc, 1000.0 * u.pc): + dm, tau = pygedm.dist_to_dm(0, 0, dist, method="ne2001") + dist_out, tau = pygedm.dm_to_dist(0, 0, dm, method="ne2001") print(dist, dm, dist_out) - assert np.isclose(dist_out.to('pc').value, dist.to('pc').value, rtol=0.1) + assert np.isclose(dist_out.to("pc").value, dist.to("pc").value, rtol=0.1) def test_igm(): - """ Test that IGM mode FAILS as expected """ + """Test that IGM mode FAILS as expected""" with pytest.raises(RuntimeError): - pygedm.dm_to_dist(100, 10, 100, mode='igm', method='ne2001') + pygedm.dm_to_dist(100, 10, 100, mode="igm", method="ne2001") def test_dm_wrapper(): - """ Run test against known values + """Run test against known values ## Test data from https://www.nrl.navy.mil/rsd/RORF/ne2001_src/ """ test_data = { - 'l': [0, 2, 97.5, ], - 'b': [0, 7.5, 85.2, ], - 'dm': [10, 20, 11.1], - 'dist': [0.461, 0.781, 0.907] + "l": [ + 0, + 2, + 97.5, + ], + "b": [ + 0, + 7.5, + 85.2, + ], + "dm": [10, 20, 11.1], + "dist": [0.461, 0.781, 0.907], } - for ii in range(len(test_data['l'])): - dist, tau_sc = pygedm.ne2001_wrapper.dm_to_dist(test_data['l'][ii], test_data['b'][ii], test_data['dm'][ii]) - assert np.allclose(dist.to('kpc').value, test_data['dist'][ii], atol=2) + for ii in range(len(test_data["l"])): + dist, tau_sc = pygedm.ne2001_wrapper.dm_to_dist( + test_data["l"][ii], test_data["b"][ii], test_data["dm"][ii] + ) + assert np.allclose(dist.to("kpc").value, test_data["dist"][ii], atol=2) - dm, tau_sc = pygedm.ne2001_wrapper.dist_to_dm(test_data['l'][ii], test_data['b'][ii], test_data['dist'][ii]) - assert np.allclose(dm.value, test_data['dm'][ii], atol=2) + dm, tau_sc = pygedm.ne2001_wrapper.dist_to_dm( + test_data["l"][ii], test_data["b"][ii], test_data["dist"][ii] + ) + assert np.allclose(dm.value, test_data["dm"][ii], atol=2) def test_zero_dm(): - """ Check that zero DM doesn't cause timeout bug + """Check that zero DM doesn't cause timeout bug Fortran code hangs if DM or D == 0 """ dist, tau_sc = pygedm.ne2001_wrapper.dm_to_dist(0, 0, 0) - dm, tau_sc = pygedm.ne2001_wrapper.dist_to_dm(0, 0, 0) + dm, tau_sc = pygedm.ne2001_wrapper.dist_to_dm(0, 0, 0) assert dist.value == dm.value == 0 def test_dm_wrapper_b0353(): - """ Test against NE2001 online values for PSR B0353+52 + """Test against NE2001 online values for PSR B0353+52 l, b, dm = 149.0993, -0.5223, 102.5 D = 2.746 kpc log_sm = -2.25 kpc m-20/3 pulse broad = 6.57 us """ - l = 149.0993 - b = -0.5223 + l = 149.0993 + b = -0.5223 dm = 102.50 dist, tau_sc = pygedm.ne2001_wrapper.dm_to_dist(l, b, dm) - assert np.isclose(dist.to('kpc').value, 2.746, atol=0.01, rtol=0) - assert np.isclose(tau_sc.to('us').value, 6.57, atol=0.01, rtol=0) + assert np.isclose(dist.to("kpc").value, 2.746, atol=0.01, rtol=0) + assert np.isclose(tau_sc.to("us").value, 6.57, atol=0.01, rtol=0) + def test_full_output(): - """ Make sure full_output arg works """ + """Make sure full_output arg works""" a = pygedm.ne2001_wrapper.dist_to_dm(0, 0, 0.1, full_output=True) b = pygedm.ne2001_wrapper.dm_to_dist(0, 0, 10, full_output=True) assert isinstance(a, dict) assert isinstance(b, dict) + if __name__ == "__main__": test_basic() test_dm_to_dist() diff --git a/tests/test_basic_ne2025.py b/tests/test_basic_ne2025.py new file mode 100644 index 0000000..89ce7e8 --- /dev/null +++ b/tests/test_basic_ne2025.py @@ -0,0 +1,163 @@ +import astropy.units as u +import numpy as np +import pytest +from astropy.coordinates import Angle +from astropy.units import Quantity, Unit + +import pygedm + + +def test_dm_to_dist(): + """Test that astropy units / angles work with dm_to_dist""" + a = pygedm.dm_to_dist(204, -6.5, 200, method="ne2025") + b = pygedm.dm_to_dist( + Angle(204, unit="degree"), Angle(-6.5, unit="degree"), 200, method="ne2025" + ) + c = pygedm.dm_to_dist(204, -6.5, 200 * Unit("pc cm^-3"), method="ne2025") + assert a[0] == b[0] == c[0] + assert a[1] == b[1] == c[1] + + +def test_dm_to_dist_ne2001p(): + """Test that dm_to_dist works with ne2001p model""" + a = pygedm.dm_to_dist(204, -6.5, 200, method="ne2001p") + b = pygedm.dm_to_dist( + Angle(204, unit="degree"), Angle(-6.5, unit="degree"), 200, method="ne2001p" + ) + c = pygedm.dm_to_dist(204, -6.5, 200 * Unit("pc cm^-3"), method="ne2001p") + assert a[0] == b[0] == c[0] + assert a[1] == b[1] == c[1] + + +def test_dist_to_dm(): + """Test that astropy units / angles work with dist_to_dm""" + a = pygedm.dist_to_dm(204, -6.5, 200, method="ne2025") + b = pygedm.dist_to_dm( + Angle(204, unit="degree"), Angle(-6.5, unit="degree"), 200, method="ne2025" + ) + c = pygedm.dist_to_dm(204, -6.5, 200 * Unit("pc"), method="ne2025") + assert a[0] == b[0] == c[0] + assert a[1] == b[1] == c[1] + + +def test_dist_to_dm_ne2001p(): + """Test that dist_to_dm works with ne2001p model""" + a = pygedm.dist_to_dm(204, -6.5, 200, method="ne2001p") + b = pygedm.dist_to_dm( + Angle(204, unit="degree"), Angle(-6.5, unit="degree"), 200, method="ne2001p" + ) + c = pygedm.dist_to_dm(204, -6.5, 200 * Unit("pc"), method="ne2001p") + assert a[0] == b[0] == c[0] + assert a[1] == b[1] == c[1] + + +def test_calculate_electron_density_xyz(): + pc = Unit("pc") + a = pygedm.calculate_electron_density_xyz(1, 2, 3, method="ne2025") + b = pygedm.calculate_electron_density_xyz(1 * pc, 2, 3, method="ne2025") + c = pygedm.calculate_electron_density_xyz(1, 2 * pc, 3, method="ne2025") + d = pygedm.calculate_electron_density_xyz(1, 2, 3 * pc, method="ne2025") + assert a == b == c == d + + +def test_basic(): + """Basic tests of NE2025 model """ + + # FRB180301 value + dm, tau = pygedm.dist_to_dm(204, -6.5, 25 * u.kpc, method="ne2025") + assert np.isclose(dm.value, 154.7895926, atol=0.01) + + # Loop through distances and check round trip + for dist in (50.0 * u.pc, 500.0 * u.pc, 5000.0 * u.pc): + dm, tau = pygedm.dist_to_dm(0, 0, dist, method="ne2025") + dist_out, tau = pygedm.dm_to_dist(0, 0, dm, method="ne2025") + print(dist, dm, dist_out) + assert np.isclose(dist_out.to("pc").value, dist.to("pc").value, rtol=0.1) + + +def test_igm(): + """Test that IGM mode FAILS as expected""" + with pytest.raises(RuntimeError): + pygedm.dm_to_dist(100, 10, 100, mode="igm", method="ne2025") + + +def test_dm_wrapper(): + """Run test against known values + + """ + test_data = { + "l": [ + 0, + 2, + 97.5, + ], + "b": [ + 0, + 7.5, + 85.2, + ], + "dm": [10, 20, 11.1], + "dist": [0.5637, 0.9096, 1.1436], + } + + for ii in range(len(test_data["l"])): + dist, tau_sc = pygedm.ne2025_wrapper.dm_to_dist( + test_data["l"][ii], test_data["b"][ii], test_data["dm"][ii] + ) + assert np.allclose(dist.to("kpc").value, test_data["dist"][ii], atol=2) + + dm, tau_sc = pygedm.ne2025_wrapper.dist_to_dm( + test_data["l"][ii], test_data["b"][ii], test_data["dist"][ii] + ) + assert np.allclose(dm.value, test_data["dm"][ii], atol=2) + + +def test_zero_dm(): + """Check that zero DM doesn't cause timeout bug + + Fortran code hangs if DM or D == 0 + """ + dist, tau_sc = pygedm.ne2025_wrapper.dm_to_dist(0, 0, 0) + dm, tau_sc = pygedm.ne2025_wrapper.dist_to_dm(0, 0, 0) + assert dist.value == dm.value == 0 + + +def test_dm_wrapper_b0353(): + """Test against precomputed values for PSR B0353+52 + + Values from NE2001 old online tool: + l, b, dm = 149.0993, -0.5223, 102.5 + D = 2.746 kpc + log_sm = -2.25 kpc m-20/3 + pulse broad = 6.57 us + + Expected values from NE2025: + Note pulse broadening estimate much larger + (, ) + """ + l = 149.0993 + b = -0.5223 + dm = 102.50 + dist, tau_sc = pygedm.ne2025_wrapper.dm_to_dist(l, b, dm) + + assert np.isclose(dist.to("kpc").value, 3.452, atol=0.01, rtol=0) + assert np.isclose(tau_sc.to("us").value, 53.63, atol=0.01, rtol=0) + + +def test_full_output(): + """Make sure full_output arg works""" + a = pygedm.ne2025_wrapper.dist_to_dm(0, 0, 0.1, full_output=True) + b = pygedm.ne2025_wrapper.dm_to_dist(0, 0, 10, full_output=True) + assert isinstance(a, dict) + assert isinstance(b, dict) + + +if __name__ == "__main__": + test_basic() + test_dm_to_dist() + test_dist_to_dm() + test_calculate_electron_density_xyz() + test_igm() + test_dm_wrapper() + test_dm_wrapper_b0353() + test_full_output() diff --git a/tests/test_basic_raises.py b/tests/test_basic_raises.py index 132273a..64d5934 100644 --- a/tests/test_basic_raises.py +++ b/tests/test_basic_raises.py @@ -1,26 +1,29 @@ -import pygedm -import numpy as np -from astropy.coordinates import Angle -from astropy.units import Unit, Quantity import astropy.units as u +import numpy as np import pytest +from astropy.coordinates import Angle +from astropy.units import Quantity, Unit + +import pygedm + def test_raises(): - """ Test that IGM mode FAILS as expected """ + """Test that IGM mode FAILS as expected""" with pytest.raises(RuntimeError): - pygedm.dm_to_dist(100, 10, 100, method='ymw1066') + pygedm.dm_to_dist(100, 10, 100, method="ymw1066") with pytest.raises(RuntimeError): - pygedm.dist_to_dm(100, 10, 100, method='ne2020') + pygedm.dist_to_dm(100, 10, 100, method="ne2020") with pytest.raises(RuntimeError): - pygedm.calculate_electron_density_xyz(100, 10, 100, method='tc93') + pygedm.calculate_electron_density_xyz(100, 10, 100, method="tc93") with pytest.raises(RuntimeError): - pygedm.calculate_electron_density_lbr(100, 10, 100, method='ymwPaleolithic') + pygedm.calculate_electron_density_lbr(100, 10, 100, method="ymwPaleolithic") with pytest.raises(RuntimeError): - pygedm.dist_to_dm(100, 10, 100, mode='igm', method='ne2001') + pygedm.dist_to_dm(100, 10, 100, mode="igm", method="ne2001") with pytest.raises(RuntimeError): - pygedm.convert_lbr_to_xyz(0, 0, 0, method='chicken') + pygedm.convert_lbr_to_xyz(0, 0, 0, method="chicken") with pytest.warns(UserWarning): - pygedm.dist_to_dm(0, 0, 555555, method='ne2001') + pygedm.dist_to_dm(0, 0, 555555, method="ne2001") + if __name__ == "__main__": test_raises() diff --git a/tests/test_basic_ymw.py b/tests/test_basic_ymw.py index 969b5ad..70f93d9 100644 --- a/tests/test_basic_ymw.py +++ b/tests/test_basic_ymw.py @@ -1,43 +1,49 @@ -import pygedm +import astropy.units as u import numpy as np from astropy.coordinates import Angle -from astropy.units import Unit, Quantity -import astropy.units as u +from astropy.units import Quantity, Unit + +import pygedm + def test_dm_to_dist(): - """ Test that astropy units / angles work with dm_to_dist """ + """Test that astropy units / angles work with dm_to_dist""" a = pygedm.dm_to_dist(204, -6.5, 200) - b = pygedm.dm_to_dist(Angle(204, unit='degree'), Angle(-6.5, unit='degree'), 200) - c = pygedm.dm_to_dist(204, -6.5, 200 * Unit('pc cm^-3')) + b = pygedm.dm_to_dist(Angle(204, unit="degree"), Angle(-6.5, unit="degree"), 200) + c = pygedm.dm_to_dist(204, -6.5, 200 * Unit("pc cm^-3")) assert a[0] == b[0] == c[0] assert a[1] == b[1] == c[1] + def test_dist_to_dm(): - """ Test that astropy units / angles work with dist_to_dm """ + """Test that astropy units / angles work with dist_to_dm""" a = pygedm.dist_to_dm(204, -6.5, 200) - b = pygedm.dist_to_dm(Angle(204, unit='degree'), Angle(-6.5, unit='degree'), 200) - c = pygedm.dist_to_dm(204, -6.5, 200 * Unit('pc')) + b = pygedm.dist_to_dm(Angle(204, unit="degree"), Angle(-6.5, unit="degree"), 200) + c = pygedm.dist_to_dm(204, -6.5, 200 * Unit("pc")) assert a[0] == b[0] == c[0] assert a[1] == b[1] == c[1] + def test_calculate_electron_density_xyz(): - pc = Unit('pc') + pc = Unit("pc") a = pygedm.calculate_electron_density_xyz(1, 2, 3) b = pygedm.calculate_electron_density_xyz(1 * pc, 2, 3) c = pygedm.calculate_electron_density_xyz(1, 2 * pc, 3) d = pygedm.calculate_electron_density_xyz(1, 2, 3 * pc) assert a == b == c == d + def test_calculate_electron_density_lbr(): - pc = Unit('pc') + pc = Unit("pc") a = pygedm.calculate_electron_density_lbr(1, 2, 3) - b = pygedm.calculate_electron_density_lbr(Angle(1, unit='deg'), 2, 3) - c = pygedm.calculate_electron_density_lbr(1, Angle(2, unit='deg'), 3) + b = pygedm.calculate_electron_density_lbr(Angle(1, unit="deg"), 2, 3) + c = pygedm.calculate_electron_density_lbr(1, Angle(2, unit="deg"), 3) d = pygedm.calculate_electron_density_lbr(1, 2, 3 * pc) assert a == b == c == d + def test_basic(): - """ Basic tests of YMW16 model + """Basic tests of YMW16 model Note: tested against online YMW16 interface http://www.atnf.csiro.au/research/pulsar/ymw16/index.php @@ -47,21 +53,21 @@ def test_basic(): assert np.isclose(a.value, 5.220655, atol=0.0001) a = pygedm.calculate_electron_density_lbr(0, 0, 4000) - assert np.isclose(a.value, 0.388407, atol=0.0001) + assert np.isclose(a.value, 0.388407, atol=0.0001) # FRB180301 value dm, tau = pygedm.dist_to_dm(204, -6.5, 25000) assert np.isclose(dm.value, 252.0501, atol=0.01) # Loop through distances and check round trip - for dist in (10., 100., 1000.): + for dist in (10.0, 100.0, 1000.0): dm, tau = pygedm.dist_to_dm(0, 0, dist) dist_out, tau = pygedm.dm_to_dist(0, 0, dm.value) assert np.isclose(dist_out.value, dist, rtol=0.1) def test_igm(): - """ Test that IGM mode works as expected + """Test that IGM mode works as expected Note: tested against YMW16 code with: # CMD: ./ ymw16 -d data -v IGM 204 -6.5 2000 100 1 @@ -69,21 +75,21 @@ def test_igm(): # z: 2.311 Dist: 5336.4 log(tau_sc): -2.218 """ - dist, tau = pygedm.dm_to_dist(204, -6.5, 2000, dm_host=100, mode='igm') + dist, tau = pygedm.dm_to_dist(204, -6.5, 2000, dm_host=100, mode="igm") assert np.isclose(dist.value, 5336.4, rtol=0.1) assert np.isclose(np.log10(tau.value), -2.218, rtol=0.1) - dm, tau = pygedm.dist_to_dm(204, -6.5, 5336.4, mode='igm') - dm_total = dm.value + 252.05 + 100 # Add galactic and host contribution + dm, tau = pygedm.dist_to_dm(204, -6.5, 5336.4, mode="igm") + dm_total = dm.value + 252.05 + 100 # Add galactic and host contribution assert np.isclose(dm_total, 2000, rtol=0.1) - dm, tau = pygedm.dist_to_dm(204, -6.5, 5336.4 * u.Mpc, mode='igm') + dm, tau = pygedm.dist_to_dm(204, -6.5, 5336.4 * u.Mpc, mode="igm") dm_total = dm.value + 252.05 + 100 # Add galactic and host contribution assert np.isclose(dm_total, 2000, rtol=0.1) def test_magellanic_cloud(): - """ Test that MC mode agrees with YMW16 + """Test that MC mode agrees with YMW16 Note: tested against YMW16 code with: CMD: ./ymw16 -d data -v MC 280.46 -32.88 50000 2 @@ -91,10 +97,11 @@ def test_magellanic_cloud(): dtest=50000.000000, nstep=10000.000000, dstep=5.000000 DM_Gal: 58.03 DM_MC: 78.87 DM: 136.90 log(tau_sc): -5.399 """ - dm, tau = pygedm.dist_to_dm(280.46, -32.88, 50000, mode='mc') + dm, tau = pygedm.dist_to_dm(280.46, -32.88, 50000, mode="mc") assert np.isclose(dm.value, 136.90) assert np.isclose(np.log10(tau.value), -5.399, rtol=0.1) + if __name__ == "__main__": test_basic() test_dm_to_dist() diff --git a/tests/test_basic_yt2020.py b/tests/test_basic_yt2020.py index 7cde310..2cbc226 100644 --- a/tests/test_basic_yt2020.py +++ b/tests/test_basic_yt2020.py @@ -1,10 +1,12 @@ -from pygedm import yt2020 -import pygedm -import numpy as np -from astropy.coordinates import Angle -from astropy.units import Unit, Quantity import astropy.units as u +import numpy as np import pytest +from astropy.coordinates import Angle +from astropy.units import Quantity, Unit + +import pygedm +from pygedm import yt2020 + def test_basic(): """ @@ -39,16 +41,18 @@ def test_basic(): DM_halo_sphe(-1.8,0.6) = 21.169654 [pc cm^{-3}] DM_halo(-1.8,0.6) = 35.989149 [pc cm^{-3} """ - test_list = [(0, 0, 220.362398, 24.857097, 245.219494), - (10, 45, 24.843628, 23.325799, 48.169427), - (-100, 32, 14.876829, 21.177893, 36.054722), - (-100.3, 32.1, 14.819495, 21.169654, 35.989149)] + test_list = [ + (0, 0, 220.362398, 24.857097, 245.219494), + (10, 45, 24.843628, 23.325799, 48.169427), + (-100, 32, 14.876829, 21.177893, 36.054722), + (-100.3, 32.1, 14.819495, 21.169654, 35.989149), + ] for line in test_list: l, b, dm_disk, dm_sphe, dm_tot = line - dm_disk_out = yt2020.calculate_halo_dm(l, b, component='disk') - dm_sphe_out = yt2020.calculate_halo_dm(l, b, component='spherical') - dm_tot_out = yt2020.calculate_halo_dm(l, b, component='both') + dm_disk_out = yt2020.calculate_halo_dm(l, b, component="disk") + dm_sphe_out = yt2020.calculate_halo_dm(l, b, component="spherical") + dm_tot_out = yt2020.calculate_halo_dm(l, b, component="both") assert np.isclose(dm_disk, dm_disk_out.value) assert np.isclose(dm_sphe, dm_sphe_out.value) @@ -56,16 +60,18 @@ def test_basic(): # Tests for higher-level wrapper in __init__.py with pytest.raises(RuntimeError): - pygedm.calculate_halo_dm(l, b, component='nonexistent_component') + pygedm.calculate_halo_dm(l, b, component="nonexistent_component") with pytest.raises(RuntimeError): - pygedm.calculate_halo_dm(l, b, method='nonexistent_component') + pygedm.calculate_halo_dm(l, b, method="nonexistent_component") for line in test_list: l, b, dm_disk, dm_sphe, dm_tot = line - dm_disk_out = pygedm.calculate_halo_dm(l, b, method='yt2020', component='disk') - dm_sphe_out = pygedm.calculate_halo_dm(l, b, method='yt2020', component='spherical') - dm_tot_out = pygedm.calculate_halo_dm(l, b, method='yt2020', component='both') + dm_disk_out = pygedm.calculate_halo_dm(l, b, method="yt2020", component="disk") + dm_sphe_out = pygedm.calculate_halo_dm( + l, b, method="yt2020", component="spherical" + ) + dm_tot_out = pygedm.calculate_halo_dm(l, b, method="yt2020", component="both") assert np.isclose(dm_disk, dm_disk_out.value) assert np.isclose(dm_sphe, dm_sphe_out.value) @@ -74,20 +80,22 @@ def test_basic(): lbd_list = ((0, 0, 250.12), (0, 30, 68.07), (15, 0, 204.9)) for line in lbd_list: l, b, dm = line - dm_out = pygedm.calculate_halo_dm(l, b, method='yt2020_analytic') + dm_out = pygedm.calculate_halo_dm(l, b, method="yt2020_analytic") assert np.isclose(dm, dm_out.value, atol=0.1) with pytest.raises(RuntimeError): - pygedm.calculate_halo_dm(l, b, method='yt2020_analytic', component='disk') + pygedm.calculate_halo_dm(l, b, method="yt2020_analytic", component="disk") + def compute_dm_halo_analytic(): - """ Test the analytical version (runs faster) """ + """Test the analytical version (runs faster)""" lbd_list = ((0, 0, 250.12), (0, 30, 68.07), (15, 0, 204.9)) for line in lbd_list: l, b, dm = line dm_out = yt2020.calculate_halo_dm_analytic(l, b) assert np.isclose(dm, dm_out.value, atol=0.1) + if __name__ == "__main__": test_basic() compute_dm_halo_analytic() diff --git a/tests/test_compare.py b/tests/test_compare.py index aff7308..6d8467d 100644 --- a/tests/test_compare.py +++ b/tests/test_compare.py @@ -1,13 +1,13 @@ import pygedm print("\n--- DM to dist ---") -print("YMW16:", pygedm.dm_to_dist(100, 0, 250, dm_host=0, method='ymw16')) -print("NE2001:", pygedm.dm_to_dist(100, 0, 250, dm_host=0, method='ne2001')) +print("YMW16:", pygedm.dm_to_dist(100, 0, 250, dm_host=0, method="ymw16")) +print("NE2001:", pygedm.dm_to_dist(100, 0, 250, dm_host=0, method="ne2001")) print("\n--- dist to DM ---") -print("YMW16:", pygedm.dist_to_dm(100, 0, 250, method='ymw16')) -print("NE2001:", pygedm.dist_to_dm(100, 0, 250, method='ne2001')) +print("YMW16:", pygedm.dist_to_dm(100, 0, 250, method="ymw16")) +print("NE2001:", pygedm.dist_to_dm(100, 0, 250, method="ne2001")) print("\n--- Electron density ---") -print("YMW16:", pygedm.calculate_electron_density_xyz(0, 0, 0, method='ymw16')) -print("NE2001:", pygedm.calculate_electron_density_xyz(0, 0, 0, method='ne2001')) \ No newline at end of file +print("YMW16:", pygedm.calculate_electron_density_xyz(0, 0, 0, method="ymw16")) +print("NE2001:", pygedm.calculate_electron_density_xyz(0, 0, 0, method="ne2001")) diff --git a/tests/test_healpix_map.py b/tests/test_healpix_map.py index 524f975..bcf3a0a 100644 --- a/tests/test_healpix_map.py +++ b/tests/test_healpix_map.py @@ -1,7 +1,9 @@ -import pygedm import numpy as np import pytest +import pygedm + + def test_healpix_map(): hpmap = pygedm.generate_healpix_dm_map(dist=30000, nside=32) @@ -10,10 +12,10 @@ def test_healpix_map(): assert np.isclose(hpmap[0], 19.607964) assert np.isclose(hpmap[1729], 29.772997) - hpmap = pygedm.generate_healpix_dm_map(nside=2, method='yt2020') + hpmap = pygedm.generate_healpix_dm_map(nside=2, method="yt2020") assert np.isclose(hpmap[0], 35.190342) - hpmap = pygedm.generate_healpix_dm_map(nside=2, method='yt2020_analytic') + hpmap = pygedm.generate_healpix_dm_map(nside=2, method="yt2020_analytic") assert np.isclose(hpmap[0], 35.01475) # Check that error is raised if healpy isn't installed @@ -22,6 +24,5 @@ def test_healpix_map(): pygedm.generate_healpix_dm_map(dist=30000, nside=32) - if __name__ == "__main__": - test_healpix_map() \ No newline at end of file + test_healpix_map() diff --git a/tests/test_lbr_to_xyz.py b/tests/test_lbr_to_xyz.py index cec026c..5b57db1 100644 --- a/tests/test_lbr_to_xyz.py +++ b/tests/test_lbr_to_xyz.py @@ -1,36 +1,42 @@ -import pygedm -import numpy as np -from astropy.coordinates import Angle, Galactocentric -from astropy.units import Unit, Quantity import astropy.units as u +import numpy as np import pytest +from astropy.coordinates import Angle, Galactocentric +from astropy.units import Quantity, Unit +import pygedm def test_xyz_to_lbr(): - x,y,z = pygedm.convert_lbr_to_xyz(0, 0, 0, method='ymw16') + x, y, z = pygedm.convert_lbr_to_xyz(0, 0, 0, method="ymw16") assert x == 0 assert y == 8300 * u.pc assert z == 6 * u.pc - x,y,z = pygedm.convert_lbr_to_xyz(0, 0, 0, method='ne2001') + x, y, z = pygedm.convert_lbr_to_xyz(0, 0, 0, method="ne2001") assert x == 0 assert y == 8500 * u.pc assert z == 0 - # Who even knows where the centre of the Galaxy is? + # Who even knows where the centre of the Galaxy is? # https://docs.astropy.org/en/stable/api/astropy.coordinates.galactocentric_frame_defaults.html from astropy.coordinates import galactocentric_frame_defaults - _ = galactocentric_frame_defaults.set('v4.0') - x,y,z = pygedm.convert_lbr_to_xyz(0, 0, 0, method='astropy') - assert np.isclose(x.to('pc').value, -1 * Galactocentric().galcen_distance.to('pc').value) + + _ = galactocentric_frame_defaults.set("v4.0") + x, y, z = pygedm.convert_lbr_to_xyz(0, 0, 0, method="astropy") + assert np.isclose( + x.to("pc").value, -1 * Galactocentric().galcen_distance.to("pc").value + ) assert y == 0 - assert np.isclose(z.to('pc').value, Galactocentric().z_sun.to('pc').value) + assert np.isclose(z.to("pc").value, Galactocentric().z_sun.to("pc").value) - x,y,z = pygedm.convert_lbr_to_xyz(Angle(0, unit='degree'), Angle(0, unit='rad'), 0 * u.pc, method='ymw16') + x, y, z = pygedm.convert_lbr_to_xyz( + Angle(0, unit="degree"), Angle(0, unit="rad"), 0 * u.pc, method="ymw16" + ) assert x == 0 assert y == 8300 * u.pc assert z == 6 * u.pc + if __name__ == "__main__": - test_xyz_to_lbr() \ No newline at end of file + test_xyz_to_lbr() diff --git a/tests/test_shared_obj.py b/tests/test_shared_obj.py index ac8ab81..29176f9 100644 --- a/tests/test_shared_obj.py +++ b/tests/test_shared_obj.py @@ -1,8 +1,13 @@ -from pygedm.ymw16_wrapper import ymw16 import os -from pkg_resources import resource_filename -DATAPATH = os.path.dirname(resource_filename("pygedm", "spiral.txt")) -gl, gb, dist, dm_host, ndir, mode_id, vbs, txt = 204.0, -6.5, 2000, 0, 2, -1, 0, '' +try: + from importlib.resources import files +except ImportError: + from importlib_resources import files + +from pygedm.ymw16_wrapper import ymw16 + +DATAPATH = str(files("pygedm")) +gl, gb, dist, dm_host, ndir, mode_id, vbs, txt = 204.0, -6.5, 2000, 0, 2, -1, 0, "" a = ymw16.dmdtau(gl, gb, dist, dm_host, ndir, mode_id, vbs, DATAPATH, txt) print(a) diff --git a/tests/test_tau_sc.py b/tests/test_tau_sc.py index ed4d036..7604e07 100644 --- a/tests/test_tau_sc.py +++ b/tests/test_tau_sc.py @@ -1,42 +1,44 @@ -import pygedm +import astropy.units as u import numpy as np from astropy.coordinates import Angle -from astropy.units import Unit, Quantity -import astropy.units as u +from astropy.units import Quantity, Unit + +import pygedm + def test_tau_sc_nu(): - """ Test that tau_sc changes with nu """ - dm, tau_sc = pygedm.dist_to_dm(0, 0, 100, method='ymw16', nu=1) - dm_, tau_sc_ = pygedm.dist_to_dm(0, 0, 100, method='ymw16', nu=1000*u.MHz) + """Test that tau_sc changes with nu""" + dm, tau_sc = pygedm.dist_to_dm(0, 0, 100, method="ymw16", nu=1) + dm_, tau_sc_ = pygedm.dist_to_dm(0, 0, 100, method="ymw16", nu=1000 * u.MHz) assert dm == dm_ assert tau_sc == tau_sc_ - dm, tau_sc = pygedm.dist_to_dm(0, 0, 100, method='ne2001', nu=1) - dm_, tau_sc_ = pygedm.dist_to_dm(0, 0, 100, method='ne2001', nu=1000*u.MHz) + dm, tau_sc = pygedm.dist_to_dm(0, 0, 100, method="ne2001", nu=1) + dm_, tau_sc_ = pygedm.dist_to_dm(0, 0, 100, method="ne2001", nu=1000 * u.MHz) assert dm == dm_ - assert tau_sc == tau_sc_ - - dist, tau_sc = pygedm.dist_to_dm(0, 0, 100, method='ymw16', nu=1) - dist_, tau_sc_ = pygedm.dist_to_dm(0, 0, 100, method='ymw16', nu=1000*u.MHz) + assert tau_sc == tau_sc_ + + dist, tau_sc = pygedm.dist_to_dm(0, 0, 100, method="ymw16", nu=1) + dist_, tau_sc_ = pygedm.dist_to_dm(0, 0, 100, method="ymw16", nu=1000 * u.MHz) assert dist == dist_ - assert tau_sc == tau_sc_ + assert tau_sc == tau_sc_ - dist, tau_sc = pygedm.dist_to_dm(0, 0, 100, method='ne2001', nu=1) - dist_, tau_sc_ = pygedm.dist_to_dm(0, 0, 100, method='ne2001', nu=1000*u.MHz) + dist, tau_sc = pygedm.dist_to_dm(0, 0, 100, method="ne2001", nu=1) + dist_, tau_sc_ = pygedm.dist_to_dm(0, 0, 100, method="ne2001", nu=1000 * u.MHz) assert dist == dist_ - assert tau_sc == tau_sc_ - - dm, tau_sc_1GHz = pygedm.dm_to_dist(0, 0, 1000, method='ymw16', nu=1.0) - dm, tau_sc_100MHz = pygedm.dm_to_dist(0, 0, 1000, method='ymw16', nu=0.1) + assert tau_sc == tau_sc_ + + dm, tau_sc_1GHz = pygedm.dm_to_dist(0, 0, 1000, method="ymw16", nu=1.0) + dm, tau_sc_100MHz = pygedm.dm_to_dist(0, 0, 1000, method="ymw16", nu=0.1) assert np.isclose(tau_sc_1GHz.value, 0.31681767) assert np.isclose(tau_sc_100MHz.value, 3168.17671061) - - assert np.isclose((0.1/1.0)**(-4) * tau_sc_1GHz.value, tau_sc_100MHz.value) - - dm, tau_sc_1GHz = pygedm.dm_to_dist(0, 0, 1000, method='ne2001', nu=1.0) - dm, tau_sc_100MHz = pygedm.dm_to_dist(0, 0, 1000, method='ne2001', nu=0.1) - assert np.isclose(tau_sc_1GHz.value, 198.57881596) - assert np.isclose(tau_sc_100MHz.value, 4988074.33385041) + + assert np.isclose((0.1 / 1.0) ** (-4) * tau_sc_1GHz.value, tau_sc_100MHz.value) + + dm, tau_sc_1GHz = pygedm.dm_to_dist(0, 0, 1000, method="ne2001", nu=1.0) + dm, tau_sc_100MHz = pygedm.dm_to_dist(0, 0, 1000, method="ne2001", nu=0.1) + assert np.isclose(tau_sc_1GHz.value, 198.55289306) + assert np.isclose(tau_sc_100MHz.value, 4987423.18015693) if __name__ == "__main__": test_tau_sc_nu()