diff --git a/.github/workflows/on_pr.yml b/.github/workflows/on_pr.yml new file mode 100644 index 0000000..b23d642 --- /dev/null +++ b/.github/workflows/on_pr.yml @@ -0,0 +1,42 @@ +name: Latest per-python environments + +on: + workflow_dispatch: + pull_request: + +jobs: + build: + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + python-version: ["3.8", "3.9", "3.10", "3.11", "3.12"] + os: [ubuntu-22.04, macos-13, macos-14] # https://docs.github.com/en/actions/using-github-hosted-runners/about-github-hosted-runners + toolchain: + - {compiler: gcc, version: 13} + + name: ${{ matrix.os }} python${{ matrix.python-version}} + steps: + - name: Checkout photodynam + uses: actions/checkout@v4 + + - name: Setup python ${{ matrix.python-version }} + uses: actions/setup-python@v5 + with: + python-version: ${{ matrix.python-version }} + + - name: Update pip + run: | + python -m pip install --upgrade pip + + - name: Install photodynam + run: | + pip install . + + - name: Test photodynam install + run: | + python -c "import photodynam" + + - name: Run tests + run: | + pytest --verbose --capture=no tests/ diff --git a/.gitignore b/.gitignore index 6dd8c2e..32cf55f 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,3 @@ .DS_Store build -photodynam +*.egg-info \ No newline at end of file diff --git a/Makefile b/Makefile index ce35e5b..7ca5aa3 100644 --- a/Makefile +++ b/Makefile @@ -4,5 +4,5 @@ CXXFLAGS=-O3 -Wall -std=c++11 -I./include all: photodynam -%: source/%.cpp +%: photodynam/%.cpp $(CXX) $< -o $@ $(CXXFLAGS) diff --git a/README.md b/README.md new file mode 100644 index 0000000..10e5203 --- /dev/null +++ b/README.md @@ -0,0 +1,263 @@ +# photodynam + +This code was written by Josh Carter (see citation requirements below) and +the public release is maintained—with the permission of the original author—by +Dan Foreman-Mackey. Any questions should be raised as issues on the GitHub +repository (https://github.com/dfm/photodynam). + +Contents +-------- + +-License +-Credit +-Overview +-Example use of functions +-photodynam + -Installing photodynam + -Documenation + - + - + -Example input file + -Example report file + -Running the example + +License +------- + + +The MIT License (MIT) + +Copyright (c) 2013 Joshua Carter + +Permission is hereby granted, free of charge, to any person obtaining a copy of +this software and associated documentation files (the "Software"), to deal in +the Software without restriction, including without limitation the rights to +use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of +the Software, and to permit persons to whom the Software is furnished to do so, +subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS +FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR +COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER +IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN +CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + +Credit +------ + +Please cite both + +Science 4 February 2011: Vol. 331 no. 6017 pp. 562-565 DOI:10.1126/science.1201274 +MNRAS (2012) 420 (2): 1630-1635. doi: 10.1111/j.1365-2966.2011.20151.x + +when using this code towards a publication. + +Overview +-------- + +The code included in this package facilitates so-called "photometric-dynamical" modeling. This model is quite simple and this is reflected in the code base. A N-body code provides coordinates and the photometric code produces light curves based on coordinates. + +The source code (in the directory source/) should be all one needs to produce forward model light curves. Briefly: + +A. Pal's code to compute overlap integrals: + elliptic.c + elliptic.h + icirc.c + icirc.h + mttr.c + scpolyint.c + scpolyint.h + +This code forms the base of the photometric code. When using that code alone or when using the full "photodynam" code please cite him appropriately: + + MNRAS (2012) 420 (2): 1630-1635. doi: 10.1111/j.1365-2966.2011.20151.x + +n_body.cpp, n_body.h: + This code performs the N-body integration with a Burlisch-Stoer integration scheme. Refer to the header for the description of the simple function "evolve." Alternatively, use the NBodyState object to access the integrator. + +n_body_state.cpp, n_body_state.h: + Defines a class (containing public member functions and constructors) that retains a N-body "state" which is most simply the ICs and masses at some time. You may operate on this object in a number of ways including N-body integration which is accomplished through the overloaded operator (). + +n_body_lc.cpp, n_body.h: + Contains code (relying on A. Pal's code, see above) to produce light curves (integrated light) for any number of spherical, quadratically limb-darkened bodies of arbitrary radius and flux. Will handle multiple overlaps ("mutual events"). + +kepcart.c: + Matt Holman's code to covert between cartesian coordinates of Keplerian elements. (Matt, credit?) + +photodynam.cpp: + User-end code to produce light curves and other information for a standardized input format. Easy to use, probably a good starting point to get some light curves produced. + +Example use of functions +------------------------ + +#include "n_body_state.h" +#include "n_body_lc.h" + +int main(int argc, char* argv[]) { + + // Kepler-16 + + int N = 3; + double t0 = 212.12316; + + // Object properties (two stars and one planet) AU, day, radians + + double mass[N] = {0.00020335520, 5.977884E-05, 9.320397E-08} + double radii[N] = {0.00301596700, 0.00104964500, 0.00035941463} + double flux[N] = {0.98474961000, 0.01525038700, 0.00000000000} + double u1[N] = {0.65139908000, 0.2, 0.0} + double u2[N] = {0.00587581200, 0.3, 0.0} + + // Now, the N-1 Jacobian Keplerian elements + + double a[N-1] = {2.240546E-01,7.040813E-01} + double e[N-1] = {1.595442E-01,7.893413E-03} + double inc[N-1] = {1.576745E+00,1.571379E+00} + double om[N-1] = {4.598385E+00,-5.374484E-01} + double ln[N-1] = {0,-8.486496E-06} + double ma[N-1] = {3.296652E+00,2.393066E+00} + + // Instantiate state. Time t0 is epoch of above coordinates + + NBodyState state(mass,a,e,inc,om,ln,ma,N,t0); + + int status; + double flux; + + // Evaluate the flux at time t0 using the getBaryLT() member method + // of NBodyState which returns NX3 array of barycentric, light-time + // corrected coordinates + + flux = occultn(state.getBaryLT(),radii,u1,u2,flux,N); + + // Now integrate forward in time to time t0+100 with stepsize 0.01 days orbit + // error tolerance of 1e-20 and minimum step size of 1e-10 days + + status = state(t0+100,0.01,1e-20,1e-10); + + // Now get the flux at the new time + + flux = occultn(state.getBaryLT(),radii,u1,u2,flux,N); + + // Print out the barycentric, light-time corrected x coordinate of body 0 + + cout << state.X_LT(0) << endl; + + return 0; +} + +// Refer to n_body_state.h for other access functions... + +photodynam +---------- + + Installing photodynam + --------------------- + + Unpack directory, change to that directory, type make. Program "photodynam" will be built + in top directory + + Documentation + ------------- + Call code as photodynam [> ]. + + Output is written to standard out unless redirected (shown in the optional listing above). + + file + ----------------- + + is file of initial coordinates and properties in + following format: + + + + + ... + ... + ... + ... + ... + + + ... + + + where the Keplerian coordinates (a = semimajor axis, e = eccentricity, i = inclination, + o = argument periapse, l = nodal longitude, m = mean anomaly) are the + N-1 Jacobian coordinates associated with the masses as ordered above. + Angles are assumed to be in radians. The observer is along the positive z axis. + Rotations are performed according to Murray and Dermott. + + file + ------------------- + + This file is a list of times to report the outputs. + The first line is a space-separated list of single character-defined + output fields according to: + + t = time + F = flux + a = semi-major axes + e = eccentricities + i = sky-plane inclinations + o = arguments of periapse + l = nodal longitudes + m = mean anomalies + K = full keplerian osculating elements + x = barycentric, light-time corrected coordinates + v = barycentric, light-time corrected velocities + M = masses + E = fractional energy change from t0 + L = fraction Lz change from t0 + + For example, the first line could be + + t F E + + and the output would have three columns of time flux and + fractional energy loss. + + + Example input file + ------------------ examples/kepler16_input.txt: + + 3 212.12316 + 0.01 1e-16 + + 0.00020335520 5.977884E-05 9.320397E-08 + 0.00301596700 0.00104964500 0.00035941463 + 0.98474961000 0.01525038700 0.00000000000 + 0.65139908000 0.2 0.0 + 0.00587581200 0.3 0.0 + + 2.240546E-01 1.595442E-01 1.576745E+00 4.598385E+00 0.000000E+00 3.296652E+00 + 7.040813E-01 7.893413E-03 1.571379E+00 -5.374484E-01 -8.486496E-06 2.393066E+00 + + Example report file + ------------------- examples/kepler16_report.txt: + + t F E e + -46.461114 -46.440679 -46.420245 -46.399811 -46.379377 + ... + + Running the example + ------------------- + + Run: + + ./photodynam examples/kepler16_input.txt examples/kepler16_report.txt + + to write output to standard out + + or + + ./photodynam examples/kepler16_input.txt examples/kepler16_report.txt > output.txt + + to dump the result into the file named output.txt + + Compare this file (whatever you call it) to examples/output.txt diff --git a/source/icirc.cpp b/photodynam/icirc.cpp similarity index 100% rename from source/icirc.cpp rename to photodynam/icirc.cpp diff --git a/source/photodynam.cpp b/photodynam/photodynam.cpp similarity index 100% rename from source/photodynam.cpp rename to photodynam/photodynam.cpp diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..590931a --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,64 @@ +# photodynam build specification +# +# Refer to the following document for specification: +# https://packaging.python.org/en/latest/specifications/ +# +# Key specification is given here: +# https://packaging.python.org/en/latest/specifications/declaring-project-metadata/#declaring-project-metadata +# +# Classifier strings are given here: +# https://pypi.org/classifiers/ + +[project] +name = "photodynam" +version = "1.0" +description = "Adapted and pythonized version of Josh Carter's photodynam code" +readme = "README.md" +requires-python = ">=3.7" +license = { text = "GPL-3.0-or-later" } +authors = [ + { name = "Josh Carter", email = "josh.a.carter@gmail.com" }, + { name = "Kyle Conroy", email = "kyleconroy@gmail.com" }, + { name = "Martin Horvat", email = "horvat77@gmail.com" }, +] +maintainers = [ + { name = "Andrej Prša", email = "aprsa@villanova.edu" }, +] +keywords = [ + "mathematics", + "science", +] +classifiers = [ + "Development Status :: 5 - Production/Stable", + "Environment :: Console", + "Intended Audience :: Developers", + "Intended Audience :: Science/Research", + "License :: OSI Approved :: GNU General Public License v3 (GPLv3)", + "Natural Language :: English", + "Operating System :: OS Independent", + "Programming Language :: C++", + "Programming Language :: Python :: 3", + "Topic :: Scientific/Engineering", + "Topic :: Software Development :: Libraries", + "Topic :: Utilities", +] +dependencies = [ + "numpy", "pytest" +] + +[project.urls] +repository = "https://github.com/phoebe-project/photodynam" + +[build-system] +requires = ["setuptools", "numpy", "wheel"] +build-backend = "setuptools.build_meta" + +[tool.setuptools] +packages = [ + "photodynam", +] + +[tool.pytest.ini_options] +addopts = [ + "--import-mode=importlib", +] diff --git a/setup.py b/setup.py index fd1bda8..5dbd3fc 100644 --- a/setup.py +++ b/setup.py @@ -1,4 +1,5 @@ -from numpy.distutils.core import setup, Extension +from setuptools import Extension, setup +import numpy # Set to true if you want to link against electric fence: CDEBUG = False @@ -8,14 +9,14 @@ libraries += ['efence'] ext_modules = [ - Extension('photodynam', - sources = ['./pywrapper/libphotodynam.cpp'], - extra_compile_args = ["-std=c++11", "-I./include"]), + Extension( + 'photodynam', + sources=[ + './pywrapper/libphotodynam.cpp' + ], + extra_compile_args=["-std=c++11", "-I./include"], + include_dirs=[numpy.get_include()] + ), ] -setup (name = 'photodynam', - version = 'devel', - description = "python wrapper of a modified version of Josh Carter's photodynam", - packages = [], - install_requires=['numpy'], - ext_modules = ext_modules) \ No newline at end of file +setup(ext_modules=ext_modules) diff --git a/tests/test_photodynam.py b/tests/test_photodynam.py new file mode 100644 index 0000000..8b03a08 --- /dev/null +++ b/tests/test_photodynam.py @@ -0,0 +1,67 @@ +import pytest + + +def test_import(): + try: + import photodynam + except ImportError: + raise 'cannot import photodynam' + + funclist = dir(photodynam) + assert 'do_dynamics' in funclist + + +def test_photodynam(): + times = (-46.461114, -46.358943, -46.256772, -46.154601, -46.052430) + masses = (0.00020335520, 5.977884E-05, 9.320397E-08) + + smas = (2.240546E-01, 7.040813E-01) + eccs = (1.595442E-01, 7.893413E-03) + incls = (1.576745E+00, 1.571379E+00) + pers = (4.598385E+00, -5.374484E-01) + longans = (0.000000E+00, -8.486496E-06) + meanans = (3.296652E+00, 2.393066E+00) + + t0 = 212.12316 + maxh = 0.01 + orbit_error = 1e-20 + ltte = 0 + return_keplerian = 0 + + expected_x = ( + (-0.0511502363666544, 0.1733783950574056, 0.40027488952282686), + (-0.05120885754146305, 0.17358032430034634, 0.3986636365083068), + (-0.05125528678276906, 0.17374078374549284, 0.3970494150725483), + (-0.051289602075730856, 0.17386003866299207, 0.3954322375092533), + (-0.05131188325061271, 0.17393836059978, 0.3938121160337223) + ) + + expected_y = ( + (2.365394719739024e-05, -7.950635374491489e-05, -0.0006153768301292659), + (2.8342166765658853e-05, -9.545552554840314e-05, -0.0006148493695065278), + (3.3023672311021896e-05, -0.00011138186510231049, -0.0006143171228620109), + (3.769740344450476e-05, -0.00012728176515662345, -0.0006137801074630697), + (4.236231708968865e-05, -0.0001431516773537286, -0.0006132383414235248) + ) + + expected_z = ( + (-0.004168119311124201, 0.013271320201116753, 0.5822134963370964), + (-0.004956660339538397, 0.01595198491740589, 0.5833603506253912), + (-0.005744066262074637, 0.018628795459728106, 0.5845025738455908), + (-0.006530158864043062, 0.021301145601418804, 0.5856401514370174), + (-0.007314762841696089, 0.0239684390182777, 0.5867730688113036) + ) + + import photodynam + res = photodynam.do_dynamics(times, masses, smas, eccs, incls, pers, longans, meanans, t0, maxh, orbit_error, ltte, return_keplerian) + + for x1, x2 in zip(res['x'], expected_x): + assert x1 == pytest.approx(x2) + for y1, y2 in zip(res['y'], expected_y): + assert y1 == pytest.approx(y2) + for z1, z2 in zip(res['z'], expected_z): + assert z1 == pytest.approx(z2) + + +if __name__ == '__main__': + test_photodynam()