diff --git a/.appveyor/install_numpy.py b/.appveyor/install_numpy_scipy.py similarity index 60% rename from .appveyor/install_numpy.py rename to .appveyor/install_numpy_scipy.py index 12bdcb5..2bd1dc7 100644 --- a/.appveyor/install_numpy.py +++ b/.appveyor/install_numpy_scipy.py @@ -10,15 +10,14 @@ from six.moves.urllib.request import urlretrieve -def download_numpy_wheel(): - base_url = os.getenv('NUMPY_URL') +def download_wheel(name, version): + base_url = os.getenv('WHEEL_URL') if base_url is None: - raise ValueError('NUMPY_URL environment variable is missing.') + raise ValueError('WHEEL_URL environment variable is missing.') - version = '1.10.4+mkl' py = 'cp{0[0]}{0[1]}'.format(sys.version_info) if py not in {'cp27', 'cp34', 'cp35'}: - print('NumPy wheel not available for {}'.format(py)) + print('{} wheel not available for {}'.format(name, py)) return None bits = struct.calcsize('P') * 8 @@ -29,26 +28,38 @@ def download_numpy_wheel(): else: raise ValueError("Couldn't determine 32/64 bits.") - filename = 'numpy-{}-{}-none-{}.whl'.format(version, py, arch) + filename = '{}-{}-{}-none-{}.whl'.format(name, version, py, arch) - directory = 'astrodynamics-numpy-wheels' - os.mkdir(directory) + directory = 'astrodynamics-wheels' + try: + os.mkdir(directory) + except OSError: + pass filepath = os.path.join(directory, filename) url = base_url + filename - # Disable SSL. Shouldn't do this ever. This is just a script. - ssl._create_default_https_context = ssl._create_unverified_context urlretrieve(url, filepath) return filepath def install_numpy(): - filepath = download_numpy_wheel() + filepath = download_wheel('numpy', '1.10.4+mkl') if filepath: pip.main(['install', filepath]) else: pip.main(['install', 'numpy']) + +def install_scipy(): + filepath = download_wheel('scipy', '0.17.0') + if filepath: + pip.main(['install', filepath]) + else: + pip.main(['install', 'scipy']) + if __name__ == '__main__': + # Disable SSL. Shouldn't do this ever. This is just a script. + ssl._create_default_https_context = ssl._create_unverified_context install_numpy() + install_scipy() diff --git a/.travis.yml b/.travis.yml index 8689be2..c95e8e1 100644 --- a/.travis.yml +++ b/.travis.yml @@ -5,6 +5,16 @@ cache: - $HOME/.cache/pip # OS X cache not yet enabled: https://github.com/travis-ci/travis-ci/issues/4011 - $HOME/Library/Caches/pip +addons: + apt: + packages: + - libatlas-dev + - libatlas-base-dev + - liblapack-dev + - gfortran +env: + global: + - TOX_TESTENV_PASSENV=CI matrix: include: - python: 2.7 @@ -24,6 +34,10 @@ matrix: addons: apt: packages: + - libatlas-dev + - libatlas-base-dev + - liblapack-dev + - gfortran - libenchant-dev - graphviz - language: generic diff --git a/.travis/install.sh b/.travis/install.sh index 4f55b14..949f757 100755 --- a/.travis/install.sh +++ b/.travis/install.sh @@ -22,12 +22,12 @@ if [[ "$(uname -s)" == 'Darwin' ]]; then pyenv global 3.3.6 ;; py34) - pyenv install 3.4.3 - pyenv global 3.4.3 + pyenv install 3.4.4 + pyenv global 3.4.4 ;; py35) - pyenv install 3.5.0 - pyenv global 3.5.0 + pyenv install 3.5.1 + pyenv global 3.5.1 ;; esac pyenv rehash @@ -38,4 +38,9 @@ fi python -m virtualenv ~/.venv source ~/.venv/bin/activate -pip install tox codecov +pip install -U tox codecov setuptools wheel + +# Build wheels for cache, but do it here so pip output prevents Travis timing out. +# We install numpy since it has to be installed before we can build scipy wheel. +pip install numpy +pip install scipy diff --git a/appveyor.yml b/appveyor.yml index f1599f4..7256c1b 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -7,7 +7,7 @@ environment: # /E:ON and /V:ON options are not enabled in the batch script intepreter # See: http://stackoverflow.com/a/13751649/163740 BUILD: "cmd /E:ON /V:ON /C .\\.appveyor\\build.cmd" - NUMPY_URL: + WHEEL_URL: secure: drCgBIMN3vdVjPRXNyiBAHa4g9VeRmnylFuhd7aUUA7XYhA8FfyJ5bHdXtzkFhgT6GzUSXqp19RcXW6DOX9eow== matrix: @@ -41,7 +41,7 @@ install: - python -m virtualenv C:\projects\venv - C:\projects\venv\Scripts\activate.bat - pip install -U codecov coverage setuptools six wheel - - "%BUILD% python .appveyor\\install_numpy.py" + - "%BUILD% python .appveyor\\install_numpy_scipy.py" - "%BUILD% pip install .[test]" build: off diff --git a/dev-requirements.txt b/dev-requirements.txt new file mode 100644 index 0000000..8a28aec --- /dev/null +++ b/dev-requirements.txt @@ -0,0 +1,18 @@ +-e .[test] +doc8 +flake8 +flake8-coding +flake8-future-import +isort +pep8-naming +plumbum>=1.6.0 +pyenchant +pytest-cov +shovel +sphinx>=1.4.0 +sphinx_rtd_theme +sphinxcontrib-spelling +tabulate +tox +twine +watchdog diff --git a/docs/astrodynamics-docs.py b/docs/astrodynamics-docs.py index 1e1bb7e..6d73746 100644 --- a/docs/astrodynamics-docs.py +++ b/docs/astrodynamics-docs.py @@ -1,7 +1,7 @@ # coding: utf-8 # This file is dual licensed under the terms of the Apache License, Version -# 2.0, and the BSD License. See the LICENSE file in the root of this repository -# for complete details. +# 2.0, and the BSD License. See the licenses/CRYPTOGRAPHY.txt file in this +# repository for complete details. # # This file is adapted from cryptography-docs.py by the The Cryptography # Developers. diff --git a/docs/conf.py b/docs/conf.py index ee23798..05b1403 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -41,6 +41,8 @@ # ones. extensions = [ 'sphinx.ext.autodoc', + 'sphinx.ext.autosectionlabel', + 'sphinx.ext.doctest', 'sphinx.ext.inheritance_diagram', 'sphinx.ext.intersphinx', 'sphinx.ext.mathjax', diff --git a/docs/development/getting-started.rst b/docs/development/getting-started.rst index 4f9c8d2..f061f11 100644 --- a/docs/development/getting-started.rst +++ b/docs/development/getting-started.rst @@ -4,15 +4,15 @@ Getting started Working on ``astrodynamics`` requires the installation of some development dependencies. These can be installed in a `virtualenv`_ -using `pip`_ with the ``dev`` extras. You should install -``astrodynamics`` in ``editable`` mode. +using `pip`_ with the ``dev-requirements.txt`` file. This installs the +dependencies and installs ``astrodynamics`` in editable mode. For example: .. code-block:: console $ # Create a virtualenv and activate it - $ pip install --editable .[dev] + $ pip install -r dev-requirements.txt You are now ready to run the tests and build the documentation. diff --git a/docs/module.rst b/docs/module.rst index 17c71b5..0b25c61 100644 --- a/docs/module.rst +++ b/docs/module.rst @@ -7,4 +7,7 @@ Module Reference modules/bodies/index modules/constants + modules/frames/index + modules/lowlevel/index + modules/rotation modules/utils/index diff --git a/docs/modules/constants.rst b/docs/modules/constants.rst index 2a9694f..0c03e40 100644 --- a/docs/modules/constants.rst +++ b/docs/modules/constants.rst @@ -136,4 +136,4 @@ WGS84_MU WGS84 geocentric gravitational constant WGS84_ANGULAR_VELOCITY WGS84 nominal earth mean angular velocity World Geodetic System 1984 ================================= =============================================== =============================================================== -.. _`license`: https://raw.githubusercontent.com/python-astrodynamics/astrodynamics/master/licenses/ASTROPY_LICENSE.txt +.. _`license`: https://raw.githubusercontent.com/python-astrodynamics/astrodynamics/master/licenses/ASTROPY.txt diff --git a/docs/modules/frames/frame.rst b/docs/modules/frames/frame.rst new file mode 100644 index 0000000..821fd94 --- /dev/null +++ b/docs/modules/frames/frame.rst @@ -0,0 +1,7 @@ +******************** +Frame Implementation +******************** + +.. automodule:: astrodynamics.frames.frame + :members: + :show-inheritance: diff --git a/docs/modules/frames/index.rst b/docs/modules/frames/index.rst new file mode 100644 index 0000000..727fbaf --- /dev/null +++ b/docs/modules/frames/index.rst @@ -0,0 +1,11 @@ +****** +Frames +****** + +.. py:module:: astrodynamics.frames + +.. toctree:: + :maxdepth: 2 + + frame + mod diff --git a/docs/modules/frames/mod.rst b/docs/modules/frames/mod.rst new file mode 100644 index 0000000..6c0cd68 --- /dev/null +++ b/docs/modules/frames/mod.rst @@ -0,0 +1,7 @@ +************************** +Mean Equinox of Date Frame +************************** + +.. automodule:: astrodynamics.frames.mod + :members: + :show-inheritance: diff --git a/docs/modules/lowlevel/index.rst b/docs/modules/lowlevel/index.rst new file mode 100644 index 0000000..1af840f --- /dev/null +++ b/docs/modules/lowlevel/index.rst @@ -0,0 +1,8 @@ +********* +Low Level +********* + +.. py:module:: astrodynamics.frames + +.. toctree:: + :maxdepth: 2 diff --git a/docs/modules/rotation.rst b/docs/modules/rotation.rst new file mode 100644 index 0000000..6a1c8a3 --- /dev/null +++ b/docs/modules/rotation.rst @@ -0,0 +1,58 @@ +********* +Rotations +********* + +.. currentmodule:: astrodynamics.rotation + +.. autoclass:: astrodynamics.rotation.Rotation + :members: __init__, from_axis_angle, from_matrix, from_euler_angles, angle, + get_axis, get_angles, matrix, apply_to, compose, distance_to + :show-inheritance: + +.. autoclass:: astrodynamics.rotation.RotationOrder + :members: + :show-inheritance: + +Rotation Convention +=================== + +Vector Operator +--------------- + +Pass ``convention='vector'`` to use the semantics of a vector operator. +According to this convention, the rotation moves vectors with respect to a fixed +reference frame. + +This means that if we define rotation r is a 90 degrees rotation around the Z +axis, the image of vector I would be J, the image of vector J would be -I, the +image of vector K would be K, and the image of vector with coordinates (1, 2, 3) +would be vector (-2, 1, 3). This means that the vector rotates counterclockwise. + +The difference with the frame transform convention is only the semantics of the +sign of the angle. It is always possible to create or use a rotation using +either convention to really represent a rotation that would have been best +created or used with the other convention, by changing accordingly the sign of +the rotation angle. + +Frame Transform +--------------- + +Pass ``convention='frame'`` to use the semantics of a frame conversion. + +According to this convention, the rotation considered vectors to be fixed, but +their coordinates change as they are converted from an initial frame to a +destination frame rotated with respect to the initial frame. + +This means that if we define rotation r is a 90 degrees rotation around the Z +axis, the image of vector I would be -J, the image of vector J would be I, the +image of vector K would be K, and the image of vector with coordinates (1, 2, 3) +would be vector (2, -1, 3). This means that the coordinates of the vector +rotates clockwise, because they are expressed with respect to a destination +frame that is rotated counterclockwise. + +The difference with vector operator convention is only the semantics of the sign +of the angle. It is always possible to create or use a rotation using either +convention to really represent a rotation that would have been best created or +used with the other convention, by changing accordingly the sign of the rotation +angle. + diff --git a/licenses/ASTROPY_LICENSE.txt b/licenses/ASTROPY.txt similarity index 100% rename from licenses/ASTROPY_LICENSE.txt rename to licenses/ASTROPY.txt diff --git a/licenses/COMMONS-MATH.txt b/licenses/COMMONS-MATH.txt new file mode 100644 index 0000000..d97b49a --- /dev/null +++ b/licenses/COMMONS-MATH.txt @@ -0,0 +1,457 @@ + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + + 1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + + 2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + + 3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + + 4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + + 5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + + 6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + + 7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + + 8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + + 9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + + END OF TERMS AND CONDITIONS + + APPENDIX: How to apply the Apache License to your work. + + To apply the Apache License to your work, attach the following + boilerplate notice, with the fields enclosed by brackets "[]" + replaced with your own identifying information. (Don't include + the brackets!) The text should be enclosed in the appropriate + comment syntax for the file format. We also recommend that a + file or class name and description of purpose be included on the + same "printed page" as the copyright notice for easier + identification within third-party archives. + + Copyright [yyyy] [name of copyright owner] + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. + + +Apache Commons Math includes the following code provided to the ASF under the +Apache License 2.0: + + - The inverse error function implementation in the Erf class is based on CUDA + code developed by Mike Giles, Oxford-Man Institute of Quantitative Finance, + and published in GPU Computing Gems, volume 2, 2010 (grant received on + March 23th 2013) + - The LinearConstraint, LinearObjectiveFunction, LinearOptimizer, + RelationShip, SimplexSolver and SimplexTableau classes in package + org.apache.commons.math3.optimization.linear include software developed by + Benjamin McCann (http://www.benmccann.com) and distributed with + the following copyright: Copyright 2009 Google Inc. (grant received on + March 16th 2009) + - The class "org.apache.commons.math3.exception.util.LocalizedFormatsTest" which + is an adapted version of "OrekitMessagesTest" test class for the Orekit library + - The "org.apache.commons.math3.analysis.interpolation.HermiteInterpolator" + has been imported from the Orekit space flight dynamics library. + +=============================================================================== + + + +APACHE COMMONS MATH DERIVATIVE WORKS: + +The Apache commons-math library includes a number of subcomponents +whose implementation is derived from original sources written +in C or Fortran. License terms of the original sources +are reproduced below. + +=============================================================================== +For the lmder, lmpar and qrsolv Fortran routine from minpack and translated in +the LevenbergMarquardtOptimizer class in package +org.apache.commons.math3.optimization.general +Original source copyright and license statement: + +Minpack Copyright Notice (1999) University of Chicago. All rights reserved + +Redistribution and use in source and binary forms, with or +without modification, are permitted provided that the +following conditions are met: + +1. Redistributions of source code must retain the above +copyright notice, this list of conditions and the following +disclaimer. + +2. Redistributions in binary form must reproduce the above +copyright notice, this list of conditions and the following +disclaimer in the documentation and/or other materials +provided with the distribution. + +3. The end-user documentation included with the +redistribution, if any, must include the following +acknowledgment: + + "This product includes software developed by the + University of Chicago, as Operator of Argonne National + Laboratory. + +Alternately, this acknowledgment may appear in the software +itself, if and wherever such third-party acknowledgments +normally appear. + +4. WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS" +WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE +UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND +THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES +OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE +OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY +OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR +USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF +THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4) +DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION +UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL +BE CORRECTED. + +5. LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT +HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF +ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT, +INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF +ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF +PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER +SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT +(INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE, +EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE +POSSIBILITY OF SUCH LOSS OR DAMAGES. +=============================================================================== + +Copyright and license statement for the odex Fortran routine developed by +E. Hairer and G. Wanner and translated in GraggBulirschStoerIntegrator class +in package org.apache.commons.math3.ode.nonstiff: + + +Copyright (c) 2004, Ernst Hairer + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + +- Redistributions of source code must retain the above copyright +notice, this list of conditions and the following disclaimer. + +- Redistributions in binary form must reproduce the above copyright +notice, this list of conditions and the following disclaimer in the +documentation and/or other materials provided with the distribution. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A +PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +=============================================================================== + +Copyright and license statement for the original Mersenne twister C +routines translated in MersenneTwister class in package +org.apache.commons.math3.random: + + Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, + All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions + are met: + + 1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + + 3. The names of its contributors may not be used to endorse or promote + products derived from this software without specific prior written + permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +=============================================================================== + +The initial code for shuffling an array (originally in class +"org.apache.commons.math3.random.RandomDataGenerator", now replaced by +a method in class "org.apache.commons.math3.util.MathArrays") was +inspired from the algorithm description provided in +"Algorithms", by Ian Craw and John Pulham (University of Aberdeen 1999). +The textbook (containing a proof that the shuffle is uniformly random) is +available here: + http://citeseerx.ist.psu.edu/viewdoc/download;?doi=10.1.1.173.1898&rep=rep1&type=pdf + +=============================================================================== +License statement for the direction numbers in the resource files for Sobol sequences. + +----------------------------------------------------------------------------- +Licence pertaining to sobol.cc and the accompanying sets of direction numbers + +----------------------------------------------------------------------------- +Copyright (c) 2008, Frances Y. Kuo and Stephen Joe +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + + * Neither the names of the copyright holders nor the names of the + University of New South Wales and the University of Waikato + and its contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY +EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY +DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND +ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +=============================================================================== + +The initial commit of package "org.apache.commons.math3.ml.neuralnet" is +an adapted version of code developed in the context of the Data Processing +and Analysis Consortium (DPAC) of the "Gaia" project of the European Space +Agency (ESA). +=============================================================================== + +The initial commit of the class "org.apache.commons.math3.special.BesselJ" is +an adapted version of code translated from the netlib Fortran program, rjbesl +http://www.netlib.org/specfun/rjbesl by R.J. Cody at Argonne National +Laboratory (USA). There is no license or copyright statement included with the +original Fortran sources. +=============================================================================== + + +The BracketFinder (package org.apache.commons.math3.optimization.univariate) +and PowellOptimizer (package org.apache.commons.math3.optimization.general) +classes are based on the Python code in module "optimize.py" (version 0.5) +developed by Travis E. Oliphant for the SciPy library (http://www.scipy.org/) +Copyright © 2003-2009 SciPy Developers. + +SciPy license +Copyright © 2001, 2002 Enthought, Inc. +All rights reserved. + +Copyright © 2003-2013 SciPy Developers. +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + + * Neither the name of Enthought nor the names of the SciPy Developers may + be used to endorse or promote products derived from this software without + specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY +EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR ANY +DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND +ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +=============================================================================== + diff --git a/licenses/CRYPTOGRAPHY_LICENSE.txt b/licenses/CRYPTOGRAPHY.txt similarity index 100% rename from licenses/CRYPTOGRAPHY_LICENSE.txt rename to licenses/CRYPTOGRAPHY.txt diff --git a/licenses/OREKIT.txt b/licenses/OREKIT.txt new file mode 100644 index 0000000..ef51da2 --- /dev/null +++ b/licenses/OREKIT.txt @@ -0,0 +1,202 @@ + + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + +TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + +1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + +2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + +3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + +4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + +5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + +6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + +7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + +8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + +9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + +END OF TERMS AND CONDITIONS + +APPENDIX: How to apply the Apache License to your work. + + To apply the Apache License to your work, attach the following + boilerplate notice, with the fields enclosed by brackets "[]" + replaced with your own identifying information. (Don't include + the brackets!) The text should be enclosed in the appropriate + comment syntax for the file format. We also recommend that a + file or class name and description of purpose be included on the + same "printed page" as the copyright notice for easier + identification within third-party archives. + +Copyright [yyyy] [name of copyright owner] + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. diff --git a/setup.py b/setup.py index 13983bf..54cbed2 100644 --- a/setup.py +++ b/setup.py @@ -3,7 +3,6 @@ import re import sys -from collections import defaultdict from setuptools import setup, find_packages from setuptools.command.test import test as TestCommand # noqa @@ -23,6 +22,7 @@ requires = { 'appdirs', 'astropy>=1.0.5', + 'boltons', 'colorama', 'docopt', 'jplephem>=2.0', @@ -31,57 +31,21 @@ 'progress', 'represent>=1.4.0', 'requests', + 'scipy', 'six', } +extras_require = dict() -def add_to_extras(extras_require, dest, source): - """Add dependencies from `source` extra to `dest` extra, handling - conditional dependencies. - """ - for key, deps in list(extras_require.items()): - extra, _, condition = key.partition(':') - if extra == source: - if condition: - extras_require[dest + ':' + condition] |= deps - else: - extras_require[dest] |= deps - -extras_require = defaultdict(set) - -extras_require[':python_version<"3.4"'] = {'pathlib'} +extras_require[':python_version<"3.4"'] = {'enum34', 'pathlib'} extras_require['test'] = { - 'pytest>=2.7.3', + 'pytest>=2.9.0', 'responses', } extras_require['test:python_version<"3.3"'] = {'mock'} -extras_require['dev'] = { - 'doc8', - 'flake8', - 'flake8-coding', - 'flake8-future-import', - 'isort', - 'pep8-naming', - 'plumbum>=1.6.0', - 'pyenchant', - 'pytest-cov', - 'shovel', - 'sphinx', - 'sphinx_rtd_theme', - 'sphinxcontrib-spelling', - 'tabulate', - 'tox', - 'twine', - 'watchdog', -} - -add_to_extras(extras_require, 'dev', 'test') - -extras_require = dict(extras_require) - class PyTest(TestCommand): user_options = [('pytest-args=', 'a', "Arguments to pass to py.test")] diff --git a/shovel/code.py b/shovel/code.py index 76edf80..f8a60c8 100644 --- a/shovel/code.py +++ b/shovel/code.py @@ -18,7 +18,7 @@ @task def format_imports(): """Sort imports into a consistent style.""" - astrodynamics_dir = Path('astrodynamics') + astrodynamics_dir = Path('src', 'astrodynamics') constants_dir = astrodynamics_dir / 'constants' for initfile in astrodynamics_dir.glob('**/__init__.py'): if constants_dir in initfile.parents: diff --git a/shovel/constants.py b/shovel/constants.py index 93302bb..2c380aa 100644 --- a/shovel/constants.py +++ b/shovel/constants.py @@ -67,7 +67,7 @@ {details_table} -.. _`license`: https://raw.githubusercontent.com/python-astrodynamics/astrodynamics/master/licenses/ASTROPY_LICENSE.txt +.. _`license`: https://raw.githubusercontent.com/python-astrodynamics/astrodynamics/master/licenses/ASTROPY.txt """ # noqa diff --git a/src/astrodynamics/compat/abc.py b/src/astrodynamics/compat/abc.py new file mode 100644 index 0000000..a5d288d --- /dev/null +++ b/src/astrodynamics/compat/abc.py @@ -0,0 +1,33 @@ +# coding: utf-8 +from __future__ import absolute_import, division, print_function + +try: + from abc import abstractstaticmethod +except ImportError: + abstractstaticmethod = None + + +class _abstractstaticmethod(staticmethod): # noqa + """ + A decorator indicating abstract staticmethods. + + Similar to abstractmethod. + + Usage: + + class C(metaclass=ABCMeta): + @abstractstaticmethod + def my_abstract_staticmethod(...): + ... + + 'abstractstaticmethod' is deprecated. Use 'staticmethod' with + 'abstractmethod' instead. + """ + + __isabstractmethod__ = True + + def __init__(self, callable): + callable.__isabstractmethod__ = True + super(_abstractstaticmethod, self).__init__(callable) + +abstractstaticmethod = abstractstaticmethod or _abstractstaticmethod diff --git a/src/astrodynamics/frames/__init__.py b/src/astrodynamics/frames/__init__.py new file mode 100644 index 0000000..086d661 --- /dev/null +++ b/src/astrodynamics/frames/__init__.py @@ -0,0 +1,22 @@ +# coding: utf-8 +# The frames implementation in astrodynamics is a derivative work of the +# implementation in Orekit. It retains the original license (see +# licenses/OREKIT.txt) +from __future__ import absolute_import, division, print_function + +from .cirf import CIRF_CONVENTIONS_2010_SIMPLE_EOP +from .eme2000 import EME2000 +from .frame import Frame, FrameProxy +from .gcrf import GCRF +from .mod import MOD_CONVENTIONS_2010_SIMPLE_EOP +from .tod import TOD_CONVENTIONS_2010_SIMPLE_EOP + +__all__ = ( + 'CIRF_CONVENTIONS_2010_SIMPLE_EOP', + 'EME2000', + 'Frame', + 'FrameProxy', + 'GCRF', + 'MOD_CONVENTIONS_2010_SIMPLE_EOP', + 'TOD_CONVENTIONS_2010_SIMPLE_EOP', +) diff --git a/src/astrodynamics/frames/cirf.py b/src/astrodynamics/frames/cirf.py new file mode 100644 index 0000000..4987047 --- /dev/null +++ b/src/astrodynamics/frames/cirf.py @@ -0,0 +1,87 @@ +# coding: utf-8 +# The frames implementation in astrodynamics is a derivative work of the +# implementation in Orekit. It retains the original license (see +# licenses/OREKIT.txt) +from __future__ import absolute_import, division, print_function + +import errno +from math import atan2, cos, sin, sqrt + +import astropy.units as u +import numpy as np +from astropy._erfa import bpn2xy, pnm06a, s06 +from astropy.utils.iers import IERS_A, IERS_A_URL +from scipy.interpolate import pchip_interpolate + +from ..rotation import Rotation +from ..utils.helper import get_neighbors_index +from ..utils.iers import IERSConventions2010 +from .frame import Frame, FrameProxy +from .gcrf import GCRF +from .transform import AbstractTransformProvider, Transform + + +class CIRFSimpleEOPTransformProvider(AbstractTransformProvider): + def __init__(self, conventions): + self.conventions = conventions + """Implementation of :class:`AbstractIERSConventions` to use.""" + + def get_transform(self, time): + # TODO: File management somewhere else + try: + iers_a = IERS_A.open() + except OSError as e: + if e.errno != errno.ENOENT: + raise + iers_a = IERS_A.open(IERS_A_URL) + + idx = get_neighbors_index(iers_a['MJD'], time.mjd, 4) + + x, y = self.conventions.get_cip_position(time) + x, y = x.si.value, y.si.value + + dx = pchip_interpolate( + iers_a['MJD'][idx], + iers_a['dX_2000A_A'][idx].to(u.rad).value, time.mjd) + + dy = pchip_interpolate( + iers_a['MJD'][idx], + iers_a['dY_2000A_A'][idx].to(u.rad).value, time.mjd) + + # Position of the Celestial Intermediate Pole (CIP) + xc = x + dx + yc = y + dy + + # Position of the Celestial Intermediate Origin (CIO) + sc = self.conventions.get_cio_locator(time, xc, yc).si.value + + # Set up the bias, precession and nutation rotation + x2py2 = xc ** 2 + yc ** 2 + zp1 = 1 + sqrt(1 - x2py2) + r = sqrt(x2py2) + spe2 = 0.5 * (sc + atan2(yc, xc)) + spe2_cos = cos(spe2) + spe2_sin = sin(spe2) + xpr = xc + r + xpr_cos = xpr * spe2_cos + xpr_sin = xpr * spe2_sin + y_cos = yc * spe2_cos + y_sin = yc * spe2_sin + bpn = Rotation(zp1 * (xpr_cos + y_sin), -r * (y_cos + xpr_sin), + r * (xpr_cos - y_sin), zp1 * (y_cos - xpr_sin)) + + return Transform(time, rot=bpn) + + +CIRF_CONVENTIONS_2010_SIMPLE_EOP = FrameProxy() + + +@CIRF_CONVENTIONS_2010_SIMPLE_EOP.register_factory +def _(): + cirf_conventions_2010_simple_eop = Frame( + parent=GCRF, + transform_provider=CIRFSimpleEOPTransformProvider(IERSConventions2010()), + name='CIRF_CONVENTIONS_2010_SIMPLE_EOP', + pseudo_intertial=True) + + return cirf_conventions_2010_simple_eop diff --git a/src/astrodynamics/frames/eme2000.py b/src/astrodynamics/frames/eme2000.py new file mode 100644 index 0000000..5ed5037 --- /dev/null +++ b/src/astrodynamics/frames/eme2000.py @@ -0,0 +1,58 @@ +# coding: utf-8 +# The frames implementation in astrodynamics is a derivative work of the +# implementation in Orekit. It retains the original license (see +# licenses/OREKIT.txt) +from __future__ import absolute_import, division, print_function + +import astropy.units as u +import numpy as np +from astropy._erfa import bi00 +from astropy.time import Time + +from ..rotation import Rotation +from .frame import Frame, FrameProxy +from .gcrf import GCRF +from .transform import FixedTransformProvider, Transform + +EME2000 = FrameProxy() + + +@EME2000.register_factory +def _(): + # Obliquity of the ecliptic. + epsilon_0 = 84381.448 * u.arcsec + + # We get the angles from bi00 directly rather than use the rb matrix from + # bp00 which returns additional matrices we don't need. + d_psi_b, d_epsilon_b, alpha_0 = bi00() + + # Longitude correction + d_psi_b = d_psi_b * u.rad + + # Obliquity correction + d_epsilon_b = d_epsilon_b * u.rad + + # the ICRS right ascension of the J2000.0 mean equinox + alpha_0 = alpha_0 * u.rad + + j2000_epoch = Time('J2000', scale='tt') + + i = np.array([1, 0, 0]) + j = np.array([0, 1, 0]) + k = np.array([0, 0, 1]) + + # Obliquity correction + r1 = Rotation.from_axis_angle(axis=i, angle=d_epsilon_b) + r2 = Rotation.from_axis_angle(axis=j, angle=-d_psi_b * np.sin(epsilon_0)) + r3 = Rotation.from_axis_angle(axis=k, angle=-alpha_0) + + transform_provider = FixedTransformProvider( + Transform(date=j2000_epoch, rot=r1.compose(r2.compose(r3)))) + + eme2000 = Frame( + parent=GCRF, + transform_provider=transform_provider, + name='EME2000', + pseudo_intertial=True) + + return eme2000 diff --git a/src/astrodynamics/frames/frame.py b/src/astrodynamics/frames/frame.py new file mode 100644 index 0000000..8f6d3f6 --- /dev/null +++ b/src/astrodynamics/frames/frame.py @@ -0,0 +1,326 @@ +# coding: utf-8 +# The frames implementation in astrodynamics is a derivative work of the +# implementation in Orekit. It retains the original license (see +# licenses/OREKIT.txt) +from __future__ import absolute_import, division, print_function + +import weakref +from abc import ABCMeta, abstractmethod, abstractproperty +from collections import Sequence + +import six +from represent import ReprHelper + +from ..compat.abc import abstractstaticmethod +from ..utils import inherit_doc, read_only_property +from .transform import Transform + + +class AncestorsDescriptor(object): + """Uses :class:`AncestorsProxy` to add an attribute (typically called + ``ancestors``) to a class. + + Parameters: + doc (str): docstring for attribute. + + The owning class must adhere to the following protocol: + + * A ``parent`` attribute referencing its parent. + * A ``depth`` attribute such that ``self.depth == parent.depth + 1``. + * An object without a parent has ``parent=None`` and ``depth=0``. + + .. doctest:: + + >>> from astrodynamics.frames.frame import AncestorsDescriptor + >>> class Node: + ... ancestors = AncestorsDescriptor("Docstring") + ... def __init__(self, parent=None): + ... self.parent = parent + ... self.depth = parent.depth + 1 if parent is not None else 0 + >>> a = Node() + >>> b = Node(a) + >>> len(b.ancestors) + 2 + >>> b.ancestors[0] == b + True + >>> b.ancestors[1] == a + True + """ + def __init__(self, doc=None): + self.instances = weakref.WeakKeyDictionary() + if doc is not None: + self.__doc__ = doc + + def __get__(self, instance, owner): + if instance is None: + return self + + if instance not in self.instances: + self.instances[instance] = AncestorsProxy(instance) + + return self.instances[instance] + + def __set__(self, instance, owner): + # By providing __set__, we are a data descriptor and show up better + # when calling help() on the owning class. This is because the docstring + # of non-data descriptors isn't shown. + raise AttributeError("can't set attribute") + + +class AncestorsProxy(Sequence): + """Given an object with ``parent`` and ``depth`` attributes, provide lazy + access to its ancestors using item access. + + Example: + .. code-block:: python + + obj = MyObject(...) + ancestors = AncestorsProxy(obj) + assert ancestors[0] == obj + assert ancestors[1] == obj.parent + + .. note:: + + The recommended usage is adding :class:`AncestorsDescriptor` to a + class which has a ``parent`` attribute. + """ + def __init__(self, instance): + self._ancestor_refs = [weakref.ref(instance)] + + def __getitem__(self, item): + # Calculate index of last item we need to return + if isinstance(item, int): + if item >= 0: + last = item + else: + last = len(self) + item + elif isinstance(item, slice): + _, stop, _ = item.indices(len(self)) + last = stop - 1 + else: + raise TypeError('Ancestor indices must be integers or slices.') + + # Try to lazily load any parents we don't have weak references to. + if last > len(self._ancestor_refs) - 1: + ref = self._ancestor_refs[-1] + current = ref() + + for i in range(len(self._ancestor_refs), last + 1): + if current.parent is None: + raise IndexError('{} has no parent.'.format(current)) + self._ancestor_refs.append(weakref.ref(current.parent)) + current = current.parent + + if isinstance(item, int): + ref = self._ancestor_refs[item] + return ref() + elif isinstance(item, slice): + return [ref() for ref in self._ancestor_refs[item]] + + def __len__(self): + ref = self._ancestor_refs[0] + return ref().depth + 1 + + +@six.add_metaclass(ABCMeta) +class AbstractFrame(object): + """Abstract Frame class.""" + @abstractproperty + def parent(self): + raise NotImplementedError + + @abstractproperty + def depth(self): + raise NotImplementedError + + @abstractproperty + def transform_provider(self): + raise NotImplementedError + + @abstractproperty + def name(self): + raise NotImplementedError + + @abstractproperty + def pseudo_intertial(self): + raise NotImplementedError + + @abstractproperty + def ancestors(self): + raise NotImplementedError + + @abstractmethod + def get_transform_to(self, destination_frame, time): + # TODO: docstring + raise NotImplementedError + + @abstractmethod + def find_common_ancestor(self, other): + # TODO: docstring + raise NotImplementedError + + +@inherit_doc.resolve +class Frame(AbstractFrame): + """Reference frame class.""" + parent = read_only_property('_parent', 'Parent frame.') + depth = read_only_property('_depth', 'Depth of frame in tree.') + transform_provider = read_only_property( + '_transform_provider', 'Provides transform from parent to instance.') + + name = read_only_property('_name', 'Name of frame.') + pseudo_intertial = read_only_property( + '_pseudo_intertial', 'Whether frame is pseudo-intertial.') + + ancestors = AncestorsDescriptor( + "An immutable sequence providing access to this frames' ancestors") + + def __init__(self, parent, transform_provider, name, pseudo_intertial): + self._parent = parent + self._depth = parent.depth + 1 if parent is not None else 0 + self._transform_provider = transform_provider + self._name = name + self._pseudo_intertial = pseudo_intertial + + @inherit_doc.mark + def get_transform_to(self, destination_frame, time): + if destination_frame == self: + return Transform(time) + + common = self.find_common_ancestor(destination_frame) + + common_to_instance = Transform(time) + for frame in self.ancestors: + if frame == common: + break + common_to_instance = ( + frame.transform_provider.get_transform(time) + + common_to_instance) + + common_to_destination = Transform(time) + for frame in destination_frame.ancestors: + if frame == common: + break + common_to_destination = ( + frame.transform_provider.get_transform(time) + + common_to_destination) + + return ~common_to_instance + common_to_destination + + @inherit_doc.mark + def find_common_ancestor(self, other): + if self.depth > other.depth: + current_from = self.ancestors[self.depth - other.depth] + current_to = other + else: + current_from = self + current_to = other.ancestors[other.depth - self.depth] + + while current_from != current_to: + current_from = current_from.parent + current_to = current_to.parent + + return current_from + + def __repr__(self): + r = ReprHelper(self) + r.parantheses = ('<', '>') + r.keyword_from_attr('name') + return str(r) + + +@inherit_doc.resolve +class FrameProxy(AbstractFrame): + """Proxy a :class:`Frame` using a factory function for lazy initialisation. + + Parameters: + factory (Frame): Optional factory function that returns the frame to be + proxied. Alternatively, the :meth:`register_factory` + method may be used as a decorator. + """ + + def __init__(self, factory=None): + self._factory = factory + self._frame = None + + def register_factory(self, factory): + """Decorator to provide :class:`FrameProxy` instance with factory + function. This is an alternative to passing a factory + + .. code-block:: python + + MyFrame = FrameProxy() + + @MyFrame.register_factory + def _(): + return Frame(...) + """ + if self._factory is not None: + raise RuntimeError('This FrameProxy already has a registered factory.') + self._factory = factory + return factory + + def _lazy_init(self): + if self._frame is None: + if self._factory is None: + raise RuntimeError('Factory not registered with this FrameProxy') + self._frame = self._factory() + + @property + def frame(self): + self._lazy_init() + return self._frame + + @property + @inherit_doc.mark + def parent(self): + self._lazy_init() + return self._frame.parent + + @property + @inherit_doc.mark + def depth(self): + self._lazy_init() + return self._frame.depth + + @property + @inherit_doc.mark + def transform_provider(self): + self._lazy_init() + return self._frame.transform_provider + + @property + @inherit_doc.mark + def name(self): + self._lazy_init() + return self._frame.name + + @property + @inherit_doc.mark + def pseudo_intertial(self): + self._lazy_init() + return self._frame.pseudo_intertial + + @property + @inherit_doc.mark + def ancestors(self): + self._lazy_init() + return self._frame.ancestors + + @inherit_doc.mark + def get_transform_to(self, destination_frame, time): + self._lazy_init() + return self._frame.get_transform_to(destination_frame, time) + + @inherit_doc.mark + def find_common_ancestor(self, other): + self._lazy_init() + return self._frame.find_common_ancestor(self, other) + + def __eq__(self, other): + if isinstance(other, FrameProxy): + return self.frame == other.frame + elif isinstance(other, Frame): + return self.frame == other + else: + return NotImplemented diff --git a/src/astrodynamics/frames/gcrf.py b/src/astrodynamics/frames/gcrf.py new file mode 100644 index 0000000..8fb49a5 --- /dev/null +++ b/src/astrodynamics/frames/gcrf.py @@ -0,0 +1,28 @@ +# coding: utf-8 +# The frames implementation in astrodynamics is a derivative work of the +# implementation in Orekit. It retains the original license (see +# licenses/OREKIT.txt) +from __future__ import absolute_import, division, print_function + +from astropy.time import Time + +from .frame import Frame, FrameProxy +from .transform import FixedTransformProvider, Transform + +GCRF = FrameProxy() + + +@GCRF.register_factory +def _(): + j2000_epoch = Time('J2000', scale='tt') + + transform_provider = FixedTransformProvider( + transform=Transform(date=j2000_epoch)) + + gcrf = Frame( + parent=None, + transform_provider=transform_provider, + name='GCRF', + pseudo_intertial=True) + + return gcrf diff --git a/src/astrodynamics/frames/mod.py b/src/astrodynamics/frames/mod.py new file mode 100644 index 0000000..2cffc7d --- /dev/null +++ b/src/astrodynamics/frames/mod.py @@ -0,0 +1,61 @@ +# coding: utf-8 +# The frames implementation in astrodynamics is a derivative work of the +# implementation in Orekit. It retains the original license (see +# licenses/OREKIT.txt) +from __future__ import absolute_import, division, print_function + +import numpy as np +from astropy.time import Time + +from ..rotation import Rotation, RotationOrder +from ..utils import inherit_doc +from ..utils.iers import IERSConventions2010 +from .eme2000 import EME2000 +from .frame import Frame, FrameProxy +from .transform import AbstractTransformProvider, Transform + + +@inherit_doc.resolve +class MODTransformProvider(AbstractTransformProvider): + """Mean Equinox of Date transform provider.""" + def __init__(self, conventions): + self.conventions = conventions + """Implementation of :class:`AbstractIERSConventions` to use.""" + + i = np.array([1, 0, 0]) + + j2000_epoch = Time('J2000', scale='tt') + + eps0 = conventions.get_mean_obliquity(j2000_epoch) + r = Rotation.from_axis_angle(i, eps0, convention='frame') + + self.ecliptic_equator_pole_rotation = r + + @inherit_doc.mark + def get_transform(self, time): + pa = self.conventions.get_precession_angles(time) + + precession = self.ecliptic_equator_pole_rotation.compose( + Rotation.from_euler_angles( + RotationOrder.ZXZ, + -pa.psi_a, + -pa.omega_a, + pa.chi_a, + convention='frame'), + convention='frame') + + return Transform(time, rot=precession) + + +MOD_CONVENTIONS_2010_SIMPLE_EOP = FrameProxy() + + +@MOD_CONVENTIONS_2010_SIMPLE_EOP.register_factory +def _(): + mod_conventions_2010_simple_eop = Frame( + parent=EME2000, + transform_provider=MODTransformProvider(IERSConventions2010()), + name='MOD_CONVENTIONS_2010_SIMPLE_EOP', + pseudo_intertial=True) + + return mod_conventions_2010_simple_eop diff --git a/src/astrodynamics/frames/tod.py b/src/astrodynamics/frames/tod.py new file mode 100644 index 0000000..66896ff --- /dev/null +++ b/src/astrodynamics/frames/tod.py @@ -0,0 +1,84 @@ +# coding: utf-8 +# The frames implementation in astrodynamics is a derivative work of the +# implementation in Orekit. It retains the original license (see +# licenses/OREKIT.txt) +from __future__ import absolute_import, division, print_function + +import errno +from collections import namedtuple + +import astropy.units as u +import numpy as np +from astropy._erfa import nut06a, obl06, p06e, xy06 +from astropy.time import Time +from astropy.utils.iers import IERS_A, IERS_A_URL +from scipy.interpolate import pchip_interpolate + +from astrodynamics.rotation import Rotation + +from ..utils.helper import get_neighbors_index, inherit_doc +from ..utils.iers import IERSConventions2010 +from .frame import Frame, FrameProxy +from .mod import MOD_CONVENTIONS_2010_SIMPLE_EOP +from .transform import AbstractTransformProvider, Transform + +TOD_CONVENTIONS_2010_SIMPLE_EOP = FrameProxy() + +I = np.array([1, 0, 0]) +K = np.array([0, 0, 1]) + + +@inherit_doc.resolve +class TODTransformProvider(AbstractTransformProvider): + def __init__(self, conventions): + self.conventions = conventions + + @inherit_doc.mark + def get_transform(self, time): + moe = self.conventions.get_mean_obliquity(time) + dpsi, deps = self.conventions.get_nutation_angles(time) + + # TODO: File management somewhere else + try: + iers_a = IERS_A.open() + except OSError as e: + if e.errno != errno.ENOENT: + raise + iers_a = IERS_A.open(IERS_A_URL) + + idx = get_neighbors_index(iers_a['MJD'], time.mjd, 4) + + mjd = iers_a['MJD'][idx] + + equinox = self.conventions.to_equinox( + time, iers_a['dX_2000A_A'][idx], iers_a['dY_2000A_A'][idx]) + + ddpsi = pchip_interpolate( + mjd, equinox.dd_psi.si.value, time.mjd) * u.rad + + ddeps = pchip_interpolate( + mjd, equinox.dd_epsilon.si.value, time.mjd) * u.rad + + dpsi += ddpsi + deps += ddeps + + toe = moe + deps + + # set up the elementary rotations for nutation + r1 = Rotation.from_axis_angle(I, toe) + r2 = Rotation.from_axis_angle(K, dpsi) + r3 = Rotation.from_axis_angle(I, -moe) + + nutation = r1.compose(r2.compose(r3)) + return Transform(time, rot=nutation) + + +@TOD_CONVENTIONS_2010_SIMPLE_EOP.register_factory +def _(): + tod_conventions_2010_simple_eop = Frame( + parent=MOD_CONVENTIONS_2010_SIMPLE_EOP, + transform_provider=TODTransformProvider(IERSConventions2010()), + name='TOD_CONVENTIONS_2010_SIMPLE_EOP', + pseudo_intertial=True) + + return tod_conventions_2010_simple_eop diff --git a/src/astrodynamics/frames/transform.py b/src/astrodynamics/frames/transform.py new file mode 100644 index 0000000..d95bbe5 --- /dev/null +++ b/src/astrodynamics/frames/transform.py @@ -0,0 +1,262 @@ +# coding: utf-8 +from __future__ import absolute_import, division, print_function + +from abc import ABCMeta, abstractmethod + +import astropy.units as u +import numpy as np +import six + +from ..rotation import Rotation +from ..utils import read_only_property, verify_unit + + +@six.add_metaclass(ABCMeta) +class AbstractTransformProvider(object): + @abstractmethod + def get_transform(self, time): + # TODO: docstring + raise NotImplementedError + + +def check_shape(v, shape): + if v.shape != shape: + raise ValueError('vector should be an array with shape {}'.format(shape)) + + +class Transform(object): + date = read_only_property('_date') + translation = read_only_property('_translation') + velocity = read_only_property('_velocity') + acceleration = read_only_property('_acceleration') + rotation = read_only_property('_rotation') + angular_velocity = read_only_property('_angular_velocity') + angular_acceleration = read_only_property('_angular_acceleration') + + def __init__(self, date, trans=None, vel=None, acc=None, rot=None, + ang_vel=None, ang_acc=None): + self._date = date + + if trans is None: + trans = np.array([0, 0, 0]) * u.m + + if vel is None: + vel = np.array([0, 0, 0]) * u.m / u.s + + if acc is None: + acc = np.array([0, 0, 0]) * u.m / u.s ** 2 + + if rot is None: + rot = Rotation(1, 0, 0, 0, normalized=True) + + if ang_vel is None: + ang_vel = np.array([0, 0, 0]) * u.rad / u.s + + if ang_acc is None: + ang_acc = np.array([0, 0, 0]) * u.rad / u.s ** 2 + + check_shape(trans, (3,)) + check_shape(vel, (3,)) + check_shape(acc, (3,)) + check_shape(ang_vel, (3,)) + check_shape(ang_acc, (3,)) + + self._translation = verify_unit(trans, 'm') + self._velocity = verify_unit(vel, 'm / s') + self._acceleration = verify_unit(acc, 'm / s2') + self._rotation = rot + self._angular_velocity = verify_unit(ang_vel, 'rad / s') + self._angular_acceleration = verify_unit(ang_acc, 'rad / s2') + + def transform_position(self, position): + verify_unit(position, 'm') + return self.rotation.apply_to(self.translation.si.value + position.si.value) * u.m + + def transform(self, position, velocity, acceleration=None): + p1 = self.translation.si.value + v1 = self.velocity.si.value + a1 = self.acceleration.si.value + o1 = self.angular_velocity.si.value + o1_dot = self.angular_acceleration.si.value + + p2 = verify_unit(position, 'm').si.value + v2 = verify_unit(velocity, 'm / s').si.value + + if acceleration is not None: + a2 = verify_unit(acceleration, 'm / s2').si.value + else: + a2 = np.zeros(3) + + p2 = vector_linear_combination(1, p2, 1, p1) + v2 = vector_linear_combination(1, v2, 1, v1) + a2 = vector_linear_combination(1, a2, 1, a1) + + pt = self.rotation.apply_to(p2) + cross_p = np.cross(o1, pt) + vt = self.rotation.apply_to(v2) - cross_p + cross_v = np.cross(o1, cross_p) + cross_cross_p = np.cross(o1, cross_p) + cross_dot_p = np.cross(o1_dot, pt) + + u1 = self.rotation.apply_to(a2) + u2 = cross_v + u3 = cross_cross_p + u4 = cross_dot_p + + at = vector_linear_combination(1, u1, -2, u2, -1, u3, -1, u4) + return pt * u.m, vt * u.m / u.s, at * u.m / u.s ** 2 + + def __add__(self, other): + if not isinstance(other, Transform): + return NotImplemented + + assert self.date == other.date + + p1 = self.translation.to(u.m).value + v1 = self.velocity.to(u.m / u.s).value + a1 = self.acceleration.to(u.m / u.s ** 2).value + r1 = self.rotation + o1 = self.angular_velocity.to(u.rad / u.s).value + o1_dot = self.angular_acceleration.to(u.rad / u.s ** 2).value + + p2 = other.translation.to(u.m).value + v2 = other.velocity.to(u.m / u.s).value + a2 = other.acceleration.to(u.m / u.s ** 2).value + r2 = other.rotation + o2 = other.angular_velocity.to(u.rad / u.s).value + o2_dot = other.angular_acceleration.to(u.rad / u.s ** 2).value + + trans = p1 + (~r1).apply_to(p2) + + o1_x_p2 = np.cross(o1, p2) + + vel = v1 + (~r1).apply_to(v2 + o1_x_p2) + + o1_x_o1_x_p2 = np.cross(o1, o1_x_p2) + o1_x_v2 = np.cross(o1, v2) + o1_dot_x_p2 = np.cross(o1_dot, p2) + + u1 = a2 + u2 = o1_x_v2 + u3 = o1_x_o1_x_p2 + u4 = o1_dot_x_p2 + + vec1 = vector_linear_combination(1, u1, 2, u2, 1, u3, 1, u4) + + acc = a1 + (~r1).apply_to(vec1) + + rot = r1.compose(r2, convention='frame') + + ang_vel = o2 + r2.apply_to(o1) + + u1 = o2_dot + u2 = r2.apply_to(o1_dot) + u3 = np.cross(o2, r2.apply_to(o1)) + + ang_acc = vector_linear_combination(1, u1, 1, u2, -1, u3) + + return Transform( + date=self.date, + trans=trans * u.m, + vel=vel * u.m / u.s, + acc=acc * u.m / u.s ** 2, + rot=rot, + ang_vel=ang_vel * u.rad / u.s, + ang_acc=ang_acc * u.rad / u.s ** 2) + + def __invert__(self): + p = self.translation.to(u.m).value + v = self.velocity.to(u.m / u.s).value + a = self.acceleration.to(u.m / u.s ** 2).value + r = self.rotation + o = self.angular_velocity.to(u.rad / u.s).value + o_dot = self.angular_acceleration.to(u.rad / u.s ** 2).value + + rp = r.apply_to(p) + rv = r.apply_to(v) + ra = r.apply_to(a) + + o_x_rp = np.cross(o, rp) + + trans = -rp + vel = o_x_rp - rv + + o_x_rv = np.cross(o, rv) + o_dot_x_rp = np.cross(o_dot, rp) + o_x_o_x_rp = np.cross(o, o_x_rp) + + u1 = ra + u2 = o_x_rv + u3 = o_dot_x_rp + u4 = o_x_o_x_rp + + acc = vector_linear_combination(-1, u1, 2, u2, 1, u3, -1, u4) + + rot = ~self.rotation + ang_vel = -rot.apply_to(o) + ang_acc = -rot.apply_to(o_dot) + + return Transform( + date=self.date, + trans=trans * u.m, + vel=vel * u.m / u.s, + acc=acc * u.m / u.s ** 2, + rot=rot, + ang_vel=ang_vel * u.rad / u.s, + ang_acc=ang_acc * u.rad / u.s ** 2) + + +def vector_linear_combination(a1, u1, a2, u2, a3=None, u3=None, a4=None, u4=None): + if u3 is None: + u3 = (None, None, None) + + if u4 is None: + u4 = (None, None, None) + + return np.array([ + linear_combination(a1, u1[0], a2, u2[0], a3, u3[0], a4, u4[0]), + linear_combination(a1, u1[1], a2, u2[1], a3, u3[1], a4, u4[1]), + linear_combination(a1, u1[2], a2, u2[2], a3, u3[2], a4, u4[2]), + ]) + + +def linear_combination(a1, b1, a2, b2, a3=None, b3=None, a4=None, b4=None): + if a4 is not None or b4 is not None: + if a4 is None or b4 is None: + raise ValueError('a4 and b4 cannot be None if either is passed.') + + if a3 is None or b3 is None: + raise ValueError('a3 and b3 cannot be None if a4 and b4 are passed.') + else: + a4 = 0 + b4 = 1 + + if a3 is not None or b3 is not None: + if a3 is None or b3 is None: + raise ValueError('a3 and b3 cannot be None if either is passed.') + else: + a3 = 0 + b3 = 1 + + x = np.array([ + [b1, 0, 0, 0], + [0, b2, 0, 0], + [0, 0, b3, 0], + [0, 0, 0, b4], + ]) + y = np.array([a1, a2, a3, a4]) + try: + result = np.sum(np.linalg.solve(x, y)) + except np.linalg.LinAlgError: + result = a1 * b1 + a2 * b2 + a3 * b3 + a4 * b4 + return result + + +class FixedTransformProvider(AbstractTransformProvider): + def __init__(self, transform): + self.transform = transform + + def get_transform(self, date): + self.transform._date = date + # TODO: 'replace' method for new date + return self.transform diff --git a/src/astrodynamics/rotation.py b/src/astrodynamics/rotation.py new file mode 100644 index 0000000..46d9ed8 --- /dev/null +++ b/src/astrodynamics/rotation.py @@ -0,0 +1,759 @@ +# coding: utf-8 +# This file is a port of some functionality from the Rotation class in +# commons-math. It retains the original license (see licenses/COMMONS-MATH.txt) +from __future__ import absolute_import, division, print_function + +from enum import Enum, unique +from math import acos, asin, atan2, cos, pi, sin, sqrt + +import astropy.units as u +import numpy as np +from represent import ReprHelperMixin +from scipy.linalg import inv, sqrtm + +from .compat.math import isclose +from .utils import read_only_property, verify_unit + +I = np.array([1, 0, 0]) +J = np.array([0, 1, 0]) +K = np.array([0, 0, 1]) + +TWOPI = 2 * pi + + +def check_convention(convention): + if convention not in ('vector', 'frame'): + raise TypeError( + "convention must be 'vector' for vector operator semantics, or " + "'frame' for frame transform semantics.") + + +class Rotation(ReprHelperMixin): + """This class implements rotations in a three-dimensional space. + + Rotations can be represented by several different mathematical entities + (matrices, axis and angle, Cardan or Euler angles, quaternions). This class + presents an higher level abstraction, more user-oriented and hiding this + implementation details. Well, for the curious, we use quaternions for the + internal representation. The user can build a rotation from any of these + representations, and any of these representations can be retrieved from a + :class:`Rotation` instance (see the various constructors and getters). In + addition, a rotation can also be built implicitly from a set of vectors and + their image. + + This implies that this class can be used to convert from one representation + to another one. For example, converting a rotation matrix into a set of + Cardan angles can be done using the following single line of code: + + .. code-block:: python + + angles = Rotation.from_matrix(matrix).get_angles(RotationOrder.XYZ) + + Focus is oriented on what a rotation *does* rather than on its underlying + representation. Once it has been built, and regardless of its internal + representation, a rotation is an *operator* which basically transforms three + dimensional vectors into other three dimensional vectors. Depending on the + application, the meaning of these vectors may vary and the semantics of the + rotation also. + + Since a rotation is basically a vectorial operator, several rotations can be + composed together and the composite operation :math:`r = r_{1} \circ r_{2}` + (which means that for each vector :math:`u`, :math:`r(u) = r_{1}(r_{2}(u))`) + is also a rotation. Hence we can consider that in addition to vectors, a + rotation can be applied to other rotations as well (or to itself). With our + previous notations, we would say we can apply :math:`r_{1}` to :math:`r_{2}` + and the result we get is :math:`r = r_{1} \circ r_{2}`. For this purpose, + the class provides the :meth:`compose` method. + + Rotations are immutable. + """ + q0 = read_only_property('_q0') + q1 = read_only_property('_q1') + q2 = read_only_property('_q2') + q3 = read_only_property('_q3') + + def __init__(self, q0, q1, q2, q3, normalized=False): + """Build a rotation from the quaternion coordinates. A rotation can be + built from a *normalized* quaternion, i.e. a quaternion for which + :math:`q_{0}^{2} + q_{1}^{2} + q_{2}^{2} + q_{3}^{2} = 1`. By default, + the constructor normalizes the quaternion. This step can be skipped + using the ``normalized`` parameter. + + Note that some conventions put the scalar part of the quaternion as the + 4th component and the vector part as the first three components. This is + *not* our convention. We put the scalar part as the first component. + + Parameters: + q0: Scalar part of the quaternion. + q1: First coordinate of the vectorial part of the quaternion. + q2: Second coordinate of the vectorial part of the quaternion. + q3: Third coordinate of the vectorial part of the quaternion. + normalized: If the coordinates are already normalized, you may pass + ``True`` to skip the normalization step. + """ + if not normalized: + inv = 1 / sqrt(q0 ** 2 + q1 ** 2 + q2 ** 2 + q3 ** 2) + q0 *= inv + q1 *= inv + q2 *= inv + q3 *= inv + + self._q0 = q0 + self._q1 = q1 + self._q2 = q2 + self._q3 = q3 + + @classmethod + def from_axis_angle(cls, axis, angle, convention='vector'): + """Build a roation from an axis and an angle. + + Parameters: + axis: Axis around which to rotate. + angle: rotation angle [rad] + + Raises: + TypeError: When the axis isn't a 1-D array with length 3. + ValueError: When the axis norm is zero. + + :type angle: :class:`~astropy.units.quantity.Quantity` + """ + check_convention(convention) + angle = verify_unit(angle, 'rad') + angle = angle.si.value + + if axis.shape != (3,): + raise TypeError('axis should be a 1-D array with length 3') + + norm = np.linalg.norm(axis) + if norm == 0: + raise ValueError('Zero norm for rotation axis.') + + half_angle = -0.5 * angle if convention == 'vector' else 0.5 * angle + coeff = sin(half_angle) / norm + + q0 = cos(half_angle) + q1 = coeff * axis[0] + q2 = coeff * axis[1] + q3 = coeff * axis[2] + + return cls(q0, q1, q2, q3, normalized=True) + + @classmethod + def from_matrix(cls, matrix): + """Build a rotation from a 3X3 matrix. + + Rotation matrices are orthogonal matrices, i.e. unit matrices + (which are matrices for which :math:`m \cdot m^{T} = I`) with real + coefficients. The module of the determinant of unit matrices is + 1, among the orthogonal 3×3 matrices, only the ones having a + positive determinant (+1) are rotation matrices. + + When a rotation is defined by a matrix with truncated values (typically + when it is extracted from a technical sheet where only four to five + significant digits are available), the matrix is not orthogonal anymore. + This constructor handles this case transparently by using a copy of the + given matrix and applying a correction to the copy in order to perfect + its orthogonality. If the Frobenius norm of the correction needed is + above the given threshold, then the matrix is considered to be too far + from a true rotation matrix and an exception is thrown. + + Parameters: + matrix: Rotation matrix. + + Raises: + TypeError: When the matrix isn't a 3×3 matrix. + ValueError: When the determinant of the computed orthogonal matrix + is negative. + + """ + if matrix.shape != (3, 3): + raise TypeError('matrix should have shape (3, 3)') + + ort = matrix.dot(inv(sqrtm(matrix.T.dot(matrix)))) + + if np.linalg.det(ort) < 0: + raise ValueError('Not a rotation matrix.') + + return cls(*_quaternion_from_matrix(ort), normalized=True) + + @classmethod + def from_euler_angles(cls, order, alpha1, alpha2, alpha3, convention='vector'): + """Build a rotation from three Euler elementary rotations. + + Cardan rotations are three successive rotations around the + canonical axes X, Y and Z, each axis being used once. There are + 6 such sets of rotations (XYZ, XZY, YXZ, YZX, ZXY and ZYX). Euler + rotations are three successive rotations around the canonical + axes X, Y and Z, the first and last rotations being around the + same axis. There are 6 such sets of rotations (XYX, XZX, YXY, + YZY, ZXZ and ZYZ), the most popular one being ZXZ. + + Beware that many people routinely use the term Euler angles even + for what really are Cardan angles (this confusion is especially + widespread in the aerospace business where Roll, Pitch and Yaw angles + are often wrongly tagged as Euler angles). + + Parameters: + order: Order of rotations to compose, from left to right. + alpha1: Angle of the first elementary rotation [rad]. + alpah2: Angle of the second elementary rotation [rad]. + alpha3: Angle of the third elementary rotation [rad]. + convention: Convention to use for the semantics of the angle. + + :type order: :class:`RotationOrder` + :type alpha1: :class:`~astropy.units.quantity.Quantity` + :type alpha2: :class:`~astropy.units.quantity.Quantity` + :type alpha3: :class:`~astropy.units.quantity.Quantity` + + .. seealso:: :ref:`Rotation Convention` + """ + check_convention(convention) + + r1 = Rotation.from_axis_angle(order.axis1, alpha1, convention) + r2 = Rotation.from_axis_angle(order.axis2, alpha2, convention) + r3 = Rotation.from_axis_angle(order.axis3, alpha3, convention) + + composed = r1.compose(r2.compose(r3, convention), convention) + + q0 = composed.q0 + q1 = composed.q1 + q2 = composed.q2 + q3 = composed.q3 + + return Rotation(q0, q1, q2, q3, normalized=True) + + @property + def angle(self): + """Angle of the rotation (between 0 and :math:`\pi`)""" + if self.q0 < -0.1 or self.q0 > 0.1: + angle = 2 * asin(sqrt(self.q1 ** 2 + self.q2 ** 2 + self.q3 ** 2)) + elif self.q0 < 0: + angle = 2 * acos(-self.q0) + else: + angle = 2 * acos(self.q0) + + return angle * u.rad + + def get_axis(self, convention='vector'): + """Get the normalized axis of the rotation. + + Note that as :attr:`Rotation.angle` always returns an angle between 0 + and :math:`\pi`, changing the convention changes the direction of the + axis, not the sign of the angle. + + Parameters: + convention: Convention to use for the semantics of the angle. + + Returns: + Normalized axis of the rotation. + + .. seealso:: :ref:`Rotation Convention` + """ + check_convention(convention) + + squared_sine = self.q1 ** 2 + self.q2 ** 2 + self.q3 ** 2 + if squared_sine == 0: + return I if convention == 'vector' else -I + else: + sign = 1 if convention == 'vector' else -1 + + if self.q0 < 0: + inverse = sign / sqrt(squared_sine) + return np.array([ + self.q1 * inverse, self.q2 * inverse, self.q3 * inverse]) + + inverse = -sign / sqrt(squared_sine) + return np.array([ + self.q1 * inverse, self.q2 * inverse, self.q3 * inverse]) + + def get_angles(self, order, convention='vector'): + """Get the Euler angles corresponding to the rotation. + + The equations show that each rotation can be defined by two different + values of the Cardan or Euler angles set. For example if Cardan angles + are used, the rotation defined by the angles :math:`a_{1}`, + :math:`a_{2}` and :math:`a_{3}` is the same as the rotation defined by + the angles :math:`\pi + a_{1}`, :math:`\pi - a_{2}` and + :math:`\pi + a_{3}`. This method implements the following arbitrary + choices: + + * for Cardan angles, the chosen set is the one for which the second + angle is between :math:`-\\frac{\pi}{2}` and :math:`\\frac{pi}{2}` + (i.e its cosine is positive), + * for Proper Euler angles, the chosen set is the one for which the + second angle is between 0 and :math:`\pi` (i.e its sine is positive). + + Cardan and Proper Euler angles have a very disappointing drawback: all + of them have singularities. This means that if the instance is too close + to the singularities corresponding to the given rotation order, it will + be impossible to retrieve the angles. For Cardan angles, this is often + called gimbal lock. There is *nothing* to do to prevent this, it is an + intrinsic problem with Cardan and Proper Euler representation (but not a + problem with the rotation itself, which is perfectly well defined). For + Cardan angles, singularities occur when the second angle is close to + :math:`\pm\\frac{\pi}{2}`, for Proper Euler angle singularities occu + when the second angle is close to 0 or :math:`\pi`, this implies that + the identity rotation is always singular for Euler angles! + + Parameters: + order: Rotation order to use. + convention: Convention to use for the semantics of the angle. + + Returns: + An tuple of three angles, in the specified order. + + Raises: + ValueError: When the rotation is singular with respect to the + specified order. + + :type order: :class:`RotationOrder` + + .. seealso:: :ref:`Rotation Convention` + """ + check_convention(convention) + + cardan_angles = { + RotationOrder.XYZ, + RotationOrder.XZY, + RotationOrder.YXZ, + RotationOrder.YZX, + RotationOrder.ZXY, + RotationOrder.ZYX, + } + + a1y_sign = 1 + a1x_sign = 1 + a2_sign = 1 + a3y_sign = 1 + a3x_sign = 1 + + if order in cardan_angles: + if order in (RotationOrder.XYZ, RotationOrder.YZX, + RotationOrder.ZXY): + a1y_sign = -1 + a3y_sign = -1 + elif order in (RotationOrder.XZY, RotationOrder.YXZ, + RotationOrder.ZYX): + a2_sign = -1 + + if convention == 'vector': + v1 = self.apply_to(order.axis3) + v2 = (~self).apply_to(order.axis1) + a1_v = v1 + a2_index = order.index3 + a3_v = v2 + else: + v1 = self.apply_to(order.axis1) + v2 = (~self).apply_to(order.axis3) + a1_v = v2 + a2_index = order.index1 + a3_v = v1 + + if isclose(abs(v2[a2_index]), 1, abs_tol=1e-15): + raise ValueError('Cardan angles singularity.') + + return ( + atan2(a1y_sign * a1_v[order.index2], + a1x_sign * a1_v[order.index3]) * u.rad, + a2_sign * asin(v2[a2_index]) * u.rad, + atan2(a3y_sign * a3_v[order.index2], + a3x_sign * a3_v[order.index1]) * u.rad + ) + else: + if order in (RotationOrder.XZX, RotationOrder.YXY, + RotationOrder.ZYZ): + a3x_sign = -1 + elif order in (RotationOrder.XYX, RotationOrder.YZY, + RotationOrder.ZXZ): + a1x_sign = -1 + + other_index = {0, 1, 2} - {order.index1, order.index2, order.index3} + other_index, = other_index + + v1 = self.apply_to(order.axis3) + v2 = (~self).apply_to(order.axis1) + + if isclose(abs(v2[order.index3]), 1, abs_tol=1e-15): + raise ValueError('Proper Euler angles singularity.') + + if convention == 'vector': + a1_v = v1 + a3_v = v2 + else: + a1_v = v2 + a3_v = v1 + + return ( + atan2(a1y_sign * a1_v[order.index2], + a1x_sign * a1_v[other_index]) * u.rad, + a2_sign * acos(v2[order.index3]) * u.rad, + atan2(a3y_sign * a3_v[order.index2], + a3x_sign * a3_v[other_index]) * u.rad + ) + + @property + def matrix(self): + """Get the 3×3 matrix corresponding to the rotation.""" + q0q0 = self.q0 * self.q0 + q0q1 = self.q0 * self.q1 + q0q2 = self.q0 * self.q2 + q0q3 = self.q0 * self.q3 + q1q1 = self.q1 * self.q1 + q1q2 = self.q1 * self.q2 + q1q3 = self.q1 * self.q3 + q2q2 = self.q2 * self.q2 + q2q3 = self.q2 * self.q3 + q3q3 = self.q3 * self.q3 + + m = np.zeros((3, 3)) + m[0][0] = 2 * (q0q0 + q1q1) - 1 + m[1][0] = 2 * (q1q2 - q0q3) + m[2][0] = 2 * (q1q3 + q0q2) + + m[0][1] = 2 * (q1q2 + q0q3) + m[1][1] = 2 * (q0q0 + q2q2) - 1 + m[2][1] = 2 * (q2q3 - q0q1) + + m[0][2] = 2 * (q1q3 - q0q2) + m[1][2] = 2 * (q2q3 + q0q1) + m[2][2] = 2 * (q0q0 + q3q3) - 1 + + return m + + def apply_to(self, vector): + """Apply the rotation to a vector. + + Parameters: + vector: Vector to apply the rotation to. + + Returns: + A new vector which is the image of `vector` by the rotation. + + Raises: + TypeError: When the vector isn't a 1-D array with length 3. + """ + + if vector.shape != (3,): + raise TypeError('vector should be a 1-D array with length 3') + + x = vector[0] + y = vector[1] + z = vector[2] + + r = self.q1 * x + self.q2 * y + self.q3 * z + s = self + + return np.array([ + 2 * (s.q0 * (x * s.q0 - (s.q2 * z - s.q3 * y)) + r * s.q1) - x, + 2 * (s.q0 * (y * s.q0 - (s.q3 * x - s.q1 * z)) + r * s.q2) - y, + 2 * (s.q0 * (z * s.q0 - (s.q1 * y - s.q2 * x)) + r * s.q3) - z + ]) + + def compose(self, rotation, convention='vector'): + """Compose the instance with another rotation. + + If the semantics of the rotation's composition corresponds to the + :ref:`Vector Operator` convention, applying the rotation to another + rotation is computing the composition in an order compliant such that + the following is true: + + .. testsetup:: + + import astropy.units as un + import numpy as np + from astrodynamics.rotation import Rotation + + r1 = Rotation.from_axis_angle(np.array([1, 0, 0]), 5 * un.deg) + r2 = Rotation.from_axis_angle(np.array([0, 1, 0]), 15 * un.deg) + + .. doctest:: + + >>> import numpy as np + >>> u = np.array([4, 5, 6]) + >>> v = r1.apply_to(u) + >>> w = r2.apply_to(v) + >>> comp = r2.compose(r1, convention='vector') + >>> np.allclose(comp.apply_to(u), w) + True + + If the semantics of the rotations' composition corresponds to the + :ref:`Frame Transform` convention, the application order will be + reversed: + + .. doctest:: + + >>> comp = r1.compose(r2, convention='frame') + >>> np.allclose(comp.apply_to(u), w) + True + + Parameters: + rotation: Rotation to apply the rotation to. + convention: Convention to use for the semantics of the angle. + + :type rotation: :class:`Rotation` + + Returns: + A new rotation which is the composition of `rotation` by the + instance. + + .. seealso:: :ref:`Rotation Convention` + """ + check_convention(convention) + + a = self + b = rotation + + if convention == 'frame': + a, b = b, a + + q0 = b.q0 * a.q0 - (b.q1 * a.q1 + b.q2 * a.q2 + b.q3 * a.q3) + q1 = b.q1 * a.q0 + b.q0 * a.q1 + (b.q2 * a.q3 - b.q3 * a.q2) + q2 = b.q2 * a.q0 + b.q0 * a.q2 + (b.q3 * a.q1 - b.q1 * a.q3) + q3 = b.q3 * a.q0 + b.q0 * a.q3 + (b.q1 * a.q2 - b.q2 * a.q1) + + return Rotation(q0, q1, q2, q3, normalized=True) + + def distance_to(self, other): + """Compute the *distance* between two rotations. + + The *distance* is intended here as a way to check if two rotations are + almost similar (i.e. they transform vectors the same way) or very + different. It is mathematically defined as the angle of the rotation + :math:`r` that prepended to one of the rotations gives the other one: + + .. math:: + + r_{1}(r) = r_{2} + + This distance is an angle between 0 and :math:`\pi` radians. Its value + is the smallest possible upper bound of the angle between + :math:`r_{1}(v)` and :math:`r_{2}(v)` for all possible vectors :math:`v`. + This upper bound is reached for some :math:`v`. The distance is equal + to 0 if an donly if the two rotations are identical. Here, :math:`r_{1}` + is the instance and :math:`r_{2}` is `other`. + + Comparing two rotations should always be done using this value rather + than for example comparing the components of the quaternions. It is much + more stable, and has a geometric meaning. Also comparing quaternions + components is error prone since for example quaternions + (0.36, 0.48, -0.48, -0.64) and (-0.36, -0.48, 0.48, 0.64) represent + exactly the same rotation despite their components are different (they + are exact opposites). + """ + return (~self).compose(other).angle + + def __invert__(self): + return Rotation(-self.q0, self.q1, self.q2, self.q3, normalized=True) + + def __eq__(self, other): + if isinstance(other, Rotation): + return all([ + self.q0 == other.q0, + self.q1 == other.q1, + self.q2 == other.q2, + self.q3 == other.q3, + ]) + else: + return NotImplemented + + def _repr_helper_(self, r): + r.keyword_from_attr('q0') + r.keyword_from_attr('q1') + r.keyword_from_attr('q2') + r.keyword_from_attr('q3') + + +def _quaternion_from_matrix(ort): + quat = np.zeros(4) + + if ort.shape != (3, 3): + raise TypeError('matrix should have shape (3, 3)') + + # There are different ways to compute the quaternions elements + # from the matrix. They all involve computing one element from + # the diagonal of the matrix, and computing the three other ones + # using a formula involving a division by the first element, + # which unfortunately can be zero. Since the norm of the + # quaternion is 1, we know at least one element has an absolute + # value greater or equal to 0.5, so it is always possible to + # select the right formula and avoid division by zero and even + # numerical inaccuracy. Checking the elements in turn and using + # the first one greater than 0.45 is safe (this leads to a simple + # test since qi = 0.45 implies 4 qi^2 - 1 = -0.19) + s = ort[0][0] + ort[1][1] + ort[2][2] + if s > -0.19: + # compute q0 and deduce q1, q2 and q3 + quat[0] = 0.5 * sqrt(s + 1) + inv = 0.25 / quat[0] + quat[1] = inv * (ort[1][2] - ort[2][1]) + quat[2] = inv * (ort[2][0] - ort[0][2]) + quat[3] = inv * (ort[0][1] - ort[1][0]) + + return quat + + s = ort[0][0] - ort[1][1] - ort[2][2] + if s > -0.19: + # compute q1 and deduce q0, q2 and q3 + quat[1] = 0.5 * sqrt(s + 1) + inv = 0.25 / quat[1] + quat[0] = inv * (ort[1][2] - ort[2][1]) + quat[2] = inv * (ort[0][1] + ort[1][0]) + quat[3] = inv * (ort[0][2] + ort[2][0]) + + return quat + + s = ort[1][1] - ort[0][0] - ort[2][2] + if s > -0.19: + # compute q2 and deduce q0, q1 and q3 + quat[2] = 0.5 * sqrt(s + 1) + inv = 0.25 / quat[2] + quat[0] = inv * (ort[2][0] - ort[0][2]) + quat[1] = inv * (ort[0][1] + ort[1][0]) + quat[3] = inv * (ort[2][1] + ort[1][2]) + + return quat + + # compute q3 and deduce q0, q1 and q2 + s = ort[2][2] - ort[0][0] - ort[1][1] + quat[3] = 0.5 * sqrt(s + 1) + inv = 0.25 / quat[3] + quat[0] = inv * (ort[0][1] - ort[1][0]) + quat[1] = inv * (ort[0][2] + ort[2][0]) + quat[2] = inv * (ort[2][1] + ort[1][2]) + + return quat + + +@unique +class RotationOrder(Enum): + """This class is a utility representing a rotation order specification for + Euler angles. + + **Cardan Angles** + + .. attribute:: XYZ + + This ordered set of rotations is around X, then around Y, then around Z. + + .. attribute:: XZY + + This ordered set of rotations is around X, then around Z, then around Y. + + .. attribute:: YXZ + + This ordered set of rotations is around Y, then around X, then around Z. + + .. attribute:: YZX + + This ordered set of rotations is around Y, then around Z, then around X. + + .. attribute:: ZXY + + This ordered set of rotations is around Z, then around X, then around Y. + + .. attribute:: ZYX + + This ordered set of rotations is around Z, then around Y, then around X. + + **Proper Euler Angles** + + .. attribute:: XYX + + This ordered set of rotations is around X, then around Y, then around X. + + .. attribute:: XZX + + This ordered set of rotations is around X, then around Z, then around X. + + .. attribute:: YXY + + This ordered set of rotations is around Y, then around X, then around Y. + + .. attribute:: YZY + + This ordered set of rotations is around Y, then around Z, then around Y. + + .. attribute:: ZXZ + + This ordered set of rotations is around Z, then around X, then around Z. + + .. attribute:: ZYZ + + This ordered set of rotations is around Z, then around Y, then around Z. + + **Instance Attributes** + + .. attribute:: axis1 + + The axis of the first rotation. + + .. attribute:: axis2 + + The axis of the second rotation. + + .. attribute:: axis3 + + The axis of the third rotation. + + .. attribute:: index1 + + The index of the first rotation. + + .. attribute:: index2 + + The index of the second rotation. + + .. attribute:: index3 + + The index of the third rotation. + + .. doctest:: + + >>> from astrodynamics.rotation import RotationOrder + >>> xyz = RotationOrder.XYZ + >>> xyz.axis1 + array([1, 0, 0]) + >>> xyz.index1, xyz.index2, xyz.index3 + (0, 1, 2) + + """ + XYZ = 1 + XZY = 2 + YXZ = 3 + YZX = 4 + ZXY = 5 + ZYX = 6 + XYX = 7 + XZX = 8 + YXY = 9 + YZY = 10 + ZXZ = 11 + ZYZ = 12 + + def __init__(self, value): + # Use the name to derive the axes. + axis1, axis2, axis3 = self.name + + # Map axis names to array indices, e.g. v[0] is x coordinate of v + indices = {'X': 0, 'Y': 1, 'Z': 2} + + # The axes in array form. + arrays = [I, J, K] + + index1 = indices[axis1] + index2 = indices[axis2] + index3 = indices[axis3] + + self.axis1 = np.array(arrays[index1]) + self.axis2 = np.array(arrays[index2]) + self.axis3 = np.array(arrays[index3]) + + self.index1 = index1 + self.index2 = index2 + self.index3 = index3 + + +def normalize_angle(a, center=pi * u.rad): + a = verify_unit(a, 'rad') + center = verify_unit(center, 'rad') + return a - TWOPI * np.floor((a + (pi * u.rad) - center) / TWOPI) diff --git a/src/astrodynamics/utils/__init__.py b/src/astrodynamics/utils/__init__.py index 913846e..772a59d 100644 --- a/src/astrodynamics/utils/__init__.py +++ b/src/astrodynamics/utils/__init__.py @@ -4,6 +4,7 @@ from .compat import PY2, PY3, PY33, WINDOWS from .helper import ( format_size, + inherit_doc, prefix, qisclose, read_only_property, @@ -28,6 +29,7 @@ 'DownloadProgressBar', 'DownloadProgressSpinner', 'format_size', + 'inherit_doc', 'InvalidCategoryError', 'KernelNotFoundError', 'prefix', diff --git a/src/astrodynamics/utils/helper.py b/src/astrodynamics/utils/helper.py index e7c3782..63aec6e 100644 --- a/src/astrodynamics/utils/helper.py +++ b/src/astrodynamics/utils/helper.py @@ -4,8 +4,10 @@ import errno from contextlib import contextmanager +import numpy as np from astropy import units as u from astropy.units import Unit, UnitBase +from boltons.typeutils import make_sentinel from ..compat.contextlib import suppress from ..compat.math import isclose @@ -147,3 +149,87 @@ def prefix(prefix, iterable): """Prepend items from `iterable` with `prefix` string.""" for x in iterable: yield '{prefix}{x}'.format(prefix=prefix, x=x) + + +inherit_doc_sentinel = make_sentinel(var_name='inherit_doc_sentinel') + + +class _InheritDoc(object): + @staticmethod + def mark(fn): + if fn.__doc__ is not None: + raise TypeError( + 'docstring for {} is not None, inherit_doc would overwrite.' + .format(fn)) + + fn.__doc__ = inherit_doc_sentinel + return fn + + @staticmethod + def resolve(cls): + sentinel_found = False + for name, func in vars(cls).items(): + + # Only inherit docstrings for functions explicitly marked + if func.__doc__ is inherit_doc_sentinel: + # We will raise exceptions so that when nothing happens, it + # doesn't do so silently. + sentinel_found = True + parfunc_found = False + + for parent in cls.__bases__: + parfunc = getattr(parent, name, None) + parfunc_found = bool(parfunc) + + if parfunc and getattr(parfunc, '__doc__', None): + func.__doc__ = parfunc.__doc__ + break + else: + # If we don't break, then replace inherit_doc_sentinel with + # None + func.__doc__ = None + pass + + if not parfunc_found: + raise TypeError( + "Cannot inherit docstring, '{}' not found in any parent" + " classes.".format(name)) + + if not sentinel_found: + raise TypeError( + 'Could not resolve docstring inheritance for {}, you must use' + ' the inherit_doc decorator.'.format(cls)) + + return cls + +inherit_doc = _InheritDoc() + + +def find_nearest_index(array, value): + idx_sorted = np.argsort(array) + sorted_array = np.array(array[idx_sorted]) + idx = np.searchsorted(sorted_array, value, side='left') + if idx >= len(array): + idx_nearest = idx_sorted[len(array) - 1] + return idx_nearest + elif idx == 0: + idx_nearest = idx_sorted[0] + return idx_nearest + else: + if (abs(value - sorted_array[idx - 1]) < + abs(value - sorted_array[idx])): + idx_nearest = idx_sorted[idx - 1] + return idx_nearest + else: + idx_nearest = idx_sorted[idx] + return idx_nearest + + +def get_neighbors_index(array, central, num_neighbors): + idx = find_nearest_index(array, central) + + start = max(0, idx - (num_neighbors - 1) // 2) + end = min(len(array), start + num_neighbors) + start = end - num_neighbors + + return slice(start, end) diff --git a/src/astrodynamics/utils/iers.py b/src/astrodynamics/utils/iers.py new file mode 100644 index 0000000..d23ddb9 --- /dev/null +++ b/src/astrodynamics/utils/iers.py @@ -0,0 +1,129 @@ +# coding: utf-8 +from __future__ import absolute_import, division, print_function + +from abc import ABCMeta, abstractmethod +from collections import namedtuple + +import numpy as np +from astropy import units as u +from astropy._erfa import bpn2xy, nut06a, obl06, p06e, pnm06a, s06 +from astropy.time import Time +from six import add_metaclass + +from .helper import inherit_doc, verify_unit + +P06eResults = namedtuple( + 'P06eResults', [ + 'eps0', + 'psia', + 'oma', + 'bpa', + 'bqa', + 'pia', + 'bpia', + 'epsa', + 'chia', + 'za', + 'zetaa', + 'thetaa', + 'pa', + 'gam', + 'phi', + 'psi']) + +PrecessionAngles = namedtuple('PrecessionAngles', ['psi_a', 'omega_a', 'chi_a']) + +Position = namedtuple('Position', ['x', 'y']) + +NutationAngles = namedtuple('NutationAngles', ['d_psi', 'd_epsilon']) + +EquinoxNutationCorrection = namedtuple( + 'EquinoxNutationCorrection', ['dd_psi', 'dd_epsilon']) + + +@add_metaclass(ABCMeta) +class AbstractIERSConventions(object): + """Abstract class for implementing IERS Conventions.""" + @abstractmethod + def get_mean_obliquity(self, time): + # TODO: docstring + raise NotImplementedError + + @abstractmethod + def get_precession_angles(self, time): + # TODO: docstring + raise NotImplementedError + + @abstractmethod + def get_cip_position(self, time): + # TODO: docstring + raise NotImplementedError + + @abstractmethod + def get_cio_locator(self, time, x, y): + # TODO: docstring + raise NotImplementedError + + @abstractmethod + def get_nutation_angles(self, time): + # TODO: docstring + raise NotImplementedError + + @abstractmethod + def to_equinox(self, time, dx, dy): + # TODO: docstring + raise NotImplementedError + + +@inherit_doc.resolve +class IERSConventions2010(AbstractIERSConventions): + """IERS Conventions 2010""" + def __init__(self): + self.nutation_reference_epoch = Time('J2000', scale='tt') + self.epsilon0 = self.get_mean_obliquity(self.nutation_reference_epoch) + + @inherit_doc.mark + def get_mean_obliquity(self, time): + return obl06(time.jd1, time.jd2) * u.rad + + @inherit_doc.mark + def get_precession_angles(self, time): + p = P06eResults(*p06e(time.jd1, time.jd2)) + return PrecessionAngles( + psi_a=p.psia * u.rad, omega_a=p.oma * u.rad, chi_a=p.chia * u.rad) + + @inherit_doc.mark + def get_cip_position(self, time): + x, y = bpn2xy(pnm06a(time.jd1, time.jd2)) + return Position(x=x * u.one, y=y * u.one) + + @inherit_doc.mark + def get_cio_locator(self, time, x, y): + x = verify_unit(x, '') + y = verify_unit(y, '') + return s06(time.jd1, time.jd2, x.value, y.value) * u.rad + + @inherit_doc.mark + def get_nutation_angles(self, time): + dpsi, deps = nut06a(time.jd1, time.jd2) + return NutationAngles(d_psi=dpsi * u.rad, d_epsilon=deps * u.rad) + + @inherit_doc.mark + def to_equinox(self, time, dx, dy): + dx = verify_unit(dx, 'rad').si.value + dy = verify_unit(dy, 'rad').si.value + + angles = self.get_precession_angles(time) + + psi_a = angles.psi_a.si.value + chi_a = angles.chi_a.si.value + + sin_ea = np.sin(self.get_mean_obliquity(time)) + cos_e0 = np.cos(self.epsilon0) + c = psi_a * cos_e0 - chi_a + o_p_c2 = 1 + c ** 2 + + ddpsi = (dx - c * dy) / (sin_ea * o_p_c2) * u.rad + ddeps = (dy + c * dx) / o_p_c2 * u.rad + + return EquinoxNutationCorrection(dd_psi=ddpsi, dd_epsilon=ddeps) diff --git a/tests/check.py b/tests/check.py new file mode 100644 index 0000000..1892c78 --- /dev/null +++ b/tests/check.py @@ -0,0 +1,21 @@ +# coding: utf-8 +from __future__ import absolute_import, division, print_function + +import astropy.units as u +import numpy as np + +from astrodynamics.compat.math import isclose +from astrodynamics.rotation import normalize_angle + + +def check_vector(v1, v2, rtol=0, atol=1e-15): + assert np.allclose(v1, v2, rtol=rtol, atol=atol) + + +def check_rotation(r1, r2, rtol=0, atol=1e-15 * u.rad): + assert np.isclose(0, r1.distance_to(r2), rtol=rtol, atol=atol) + + +def check_angle(a1, a2, rtol=0, atol=1e-12 * u.rad): + assert np.isclose( + a1, normalize_angle(a2, center=a1), rtol=rtol, atol=atol) diff --git a/tests/frames/__init__.py b/tests/frames/__init__.py new file mode 100644 index 0000000..393def6 --- /dev/null +++ b/tests/frames/__init__.py @@ -0,0 +1,2 @@ +# coding: utf-8 +from __future__ import absolute_import, division, print_function diff --git a/tests/frames/test_cirf.py b/tests/frames/test_cirf.py new file mode 100644 index 0000000..97a04c3 --- /dev/null +++ b/tests/frames/test_cirf.py @@ -0,0 +1,28 @@ +# coding: utf-8 +from __future__ import absolute_import, division, print_function + +import astropy.units as u +import numpy as np +import pytest +from astropy.time import Time +from numpy.linalg import norm + +from astrodynamics.frames import CIRF_CONVENTIONS_2010_SIMPLE_EOP + + +@pytest.mark.skip( + reason='CIRF not fully implemented yet, requires IERS download each time.') +def test_rotation_rate(): + t_min = Time('2009-04-07 02:56:33.816', scale='utc') + t_max = Time('2043-12-16 10:47:20', scale='utc') + + t1 = CIRF_CONVENTIONS_2010_SIMPLE_EOP.transform_provider.get_transform(t_min) + t2 = CIRF_CONVENTIONS_2010_SIMPLE_EOP.transform_provider.get_transform(t_max) + + rad_s = u.rad / u.s + + min_rate = norm(t_min.angular_velocity) + max_rate = norm(t_max.angular_velocity) + + assert np.isclose(min_rate, 1.1e-15 * rad_s, rtol=0, atol=1.0e-16 * rad_s) + assert np.isclose(max_rate, 8.6e-12 * rad_s, rtol=0, atol=1.0e-13 * rad_s) diff --git a/tests/frames/test_eme2000.py b/tests/frames/test_eme2000.py new file mode 100644 index 0000000..053a319 --- /dev/null +++ b/tests/frames/test_eme2000.py @@ -0,0 +1,84 @@ +# coding: utf-8 +from __future__ import absolute_import, division, print_function + +import astropy.units as u +import numpy as np +from astropy.time import Time +from astropy._erfa import bp00 + +from astrodynamics.frames import GCRF, EME2000 +from astrodynamics.rotation import Rotation + +from tests.check import check_rotation, check_vector + + +def test_aas_reference_leo(): + """This reference test has been extracted from the following paper: + Implementation Issues Surrounding the New IAU Reference Systems for + Astrodynamics, by David A. Vallado, John H. Seago, P. Kenneth Seidelmann + + http://www.centerforspace.com/downloads/files/pubs/AAS-06-134.pdf + """ + t0 = Time('2004-04-06 07:51:28.386009', scale='utc') + t = GCRF.get_transform_to(EME2000, t0) + + p_gcrf_iau_2000_a = np.array([5102508.9579, 6123011.4038, 6378136.9252]) * u.m + v_gcrf_iau_2000_a = np.array([-4743.220156, 790.536497, 5533.755728]) * u.m / u.s + p_eme2000_eq_a = np.array([5102509.0383, 6123011.9758, 6378136.3118]) * u.m + v_eme2000_eq_a = np.array([-4743.219766, 790.536344, 5533.756084]) * u.m / u.s + + p, v, _ = t.transform(p_gcrf_iau_2000_a, v_gcrf_iau_2000_a) + + check_vector(p_eme2000_eq_a, p, rtol=0, atol=1.1e-4 * u.m) + check_vector(v_eme2000_eq_a, v, rtol=0, atol=2.6e-7 * u.m / u.s) + + p_gcrf_iau_2000_b = np.array([5102508.9579, 6123011.4012, 6378136.9277]) * u.m + v_gcrf_iau_2000_b = np.array([-4743.220156, 790.536495, 5533.755729]) * u.m / u.s + + p_eme2000_eq_b = np.array([5102509.0383, 6123011.9733, 6378136.3142]) * u.m + v_eme2000_eq_b = np.array([-4743.219766, 790.536342, 5533.756085]) * u.m / u.s + + p, v, _ = t.transform(p_gcrf_iau_2000_b, v_gcrf_iau_2000_b) + + check_vector(p_eme2000_eq_b, p, rtol=0, atol=7.4e-5 * u.m) + check_vector(v_eme2000_eq_b, v, rtol=0, atol=2.6e-7 * u.m / u.s) + + +def test_aas_reference_geo(): + """This reference test has been extracted from the following paper: + Implementation Issues Surrounding the New IAU Reference Systems for + Astrodynamics, by David A. Vallado, John H. Seago, P. Kenneth Seidelmann + + http://www.centerforspace.com/downloads/files/pubs/AAS-06-134.pdf + """ + t0 = Time('2004-06-01', scale='utc') + t = GCRF.get_transform_to(EME2000, t0) + + p_gcrf_iau_2000_a = np.array([-40588150.3617, -11462167.0397, 27143.1974]) * u.m + v_gcrf_iau_2000_a = np.array([834.787458, -2958.305691, -1.172993]) * u.m / u.s + p_eme2000_eq_a = np.array([-40588149.5482, -11462169.9118, 27146.8462]) * u.m + v_eme2000_eq_a = np.array([834.787667, -2958.305632, -1.172963]) * u.m / u.s + + p, v, _ = t.transform(p_gcrf_iau_2000_a, v_gcrf_iau_2000_a) + + check_vector(p_eme2000_eq_a, p, rtol=0, atol=5.8e-5 * u.m) + check_vector(v_eme2000_eq_a, v, rtol=0, atol=6.4e-7 * u.m / u.s) + + p_gcrf_iau_2000_b = np.array([-40588150.3617, -11462167.0397, 27143.2125]) * u.m + v_gcrf_iau_2000_b = np.array([834.787458, -2958.305691, -1.172999]) * u.m / u.s + + p_eme2000_eq_b = np.array([-40588149.5481, -11462169.9118, 27146.8613]) * u.m + v_eme2000_eq_b = np.array([834.787667, -2958.305632, -1.172968]) * u.m / u.s + + p, v, _ = t.transform(p_gcrf_iau_2000_b, v_gcrf_iau_2000_b) + + check_vector(p_eme2000_eq_b, p, rtol=0, atol=1.1e-4 * u.m) + check_vector(v_eme2000_eq_b, v, rtol=0, atol=5.5e-7 * u.m / u.s) + + +def test_erfa_bp00(): + t = Time('2004-02-14', scale='utc') + rb, rp, rbp = bp00(t.jd1, t.jd2) + r1 = Rotation.from_matrix(rb) + r2 = GCRF.get_transform_to(EME2000, t).rotation + check_rotation(r1, r2, atol=2.5e-16 * u.rad) diff --git a/tests/frames/test_mod.py b/tests/frames/test_mod.py new file mode 100644 index 0000000..cf2369b --- /dev/null +++ b/tests/frames/test_mod.py @@ -0,0 +1,20 @@ +# coding: utf-8 +from __future__ import absolute_import, division, print_function + +import astropy.units as u +import numpy as np +from astropy.time import Time + +from astrodynamics.frames import GCRF, MOD_CONVENTIONS_2010_SIMPLE_EOP +from astrodynamics.rotation import Rotation + +from tests.check import check_rotation, check_vector + + +def test_erfa_bp06(): + t = Time('2004-02-14', scale='utc') + from astropy._erfa import bp06 + rb, rp, rbp = bp06(t.jd1, t.jd2) + m1 = rbp + m2 = GCRF.get_transform_to(MOD_CONVENTIONS_2010_SIMPLE_EOP, t).rotation.matrix + assert np.allclose(m1, m2, rtol=0, atol=1.1e-12) diff --git a/tests/frames/test_tod.py b/tests/frames/test_tod.py new file mode 100644 index 0000000..2a22a9a --- /dev/null +++ b/tests/frames/test_tod.py @@ -0,0 +1,3 @@ +# coding: utf-8 +from __future__ import absolute_import, division, print_function + diff --git a/tests/test_constants.py b/tests/test_constants.py index f5ddf7a..d9ddab0 100644 --- a/tests/test_constants.py +++ b/tests/test_constants.py @@ -1,6 +1,6 @@ # coding: utf-8 # These tests are taken from astropy, as with the astrodynamics.constant.Constant -# class. It retains the original license (see licenses/ASTROPY_LICENSE.txt) +# class. It retains the original license (see licenses/ASTROPY.txt) from __future__ import absolute_import, division, print_function import copy diff --git a/tests/test_rotation.py b/tests/test_rotation.py new file mode 100644 index 0000000..e73f194 --- /dev/null +++ b/tests/test_rotation.py @@ -0,0 +1,249 @@ +# coding: utf-8 +from __future__ import absolute_import, division, print_function + +from math import pi, sqrt + +import astropy.units as u +import itertools +import numpy as np +import pytest +from astrodynamics.compat.math import isclose +from astrodynamics.rotation import Rotation, RotationOrder + +from .check import check_angle, check_rotation, check_vector + +I = np.array([1, 0, 0]) +J = np.array([0, 1, 0]) +K = np.array([0, 0, 1]) + + +def test_axis_angle_vector_convention(): + r = Rotation.from_axis_angle( + np.array([10, 10, 10]), 2 * pi / 3 * u.rad, convention='vector') + check_vector(r.apply_to(I), J) + check_vector(r.apply_to(J), K) + check_vector(r.apply_to(K), I) + + s = 1 / sqrt(3) + check_vector(r.get_axis(convention='vector'), np.array([s, s, s])) + check_vector(r.get_axis(convention='frame'), np.array([-s, -s, -s])) + + check_angle(r.angle, 2 * pi / 3 * u.rad) + + with pytest.raises(ValueError): + Rotation.from_axis_angle( + np.array([0, 0, 0]), 2 * pi / 3, convention='vector') + + r = Rotation.from_axis_angle(K, 1.5 * pi * u.rad, convention='vector') + check_vector(r.get_axis(convention='vector'), -K) + check_vector(r.get_axis(convention='frame'), K) + check_angle(r.angle, 0.5 * pi * u.rad) + + r = Rotation.from_axis_angle(J, pi * u.rad, convention='vector') + check_vector(r.get_axis(convention='vector'), J) + check_vector(r.get_axis(convention='frame'), -J) + check_angle(r.angle, pi * u.rad) + + +def test_axis_angle_frame_convention(): + r = Rotation.from_axis_angle( + np.array([10, 10, 10]), 2 * pi / 3 * u.rad, convention='frame') + check_vector(r.apply_to(I), K) + check_vector(r.apply_to(J), I) + check_vector(r.apply_to(K), J) + + s = 1 / sqrt(3) + check_vector(r.get_axis(convention='frame'), np.array([s, s, s])) + check_vector(r.get_axis(convention='vector'), np.array([-s, -s, -s])) + check_angle(r.angle, 2 * pi / 3 * u.rad) + + with pytest.raises(ValueError): + Rotation.from_axis_angle( + np.array([0, 0, 0]), 2 * pi / 3 * u.rad, convention='frame') + + r = Rotation.from_axis_angle(K, 1.5 * pi * u.rad, convention='frame') + check_vector(r.get_axis(convention='frame'), -K) + check_vector(r.get_axis(convention='vector'), K) + check_angle(r.angle, 0.5 * pi * u.rad) + + r = Rotation.from_axis_angle(J, pi * u.rad, convention='frame') + check_vector(r.get_axis(convention='frame'), J) + check_vector(r.get_axis(convention='vector'), -J) + check_angle(r.angle, pi * u.rad) + + +@pytest.mark.parametrize('convention', ['vector', 'frame']) +def test_invert(convention): + r = Rotation(0.001, 0.36, 0.48, 0.8) + inv = ~r + check_rotation(r.compose(inv, convention=convention), Rotation(1, 0, 0, 0)) + check_rotation(inv.compose(r, convention=convention), Rotation(1, 0, 0, 0)) + check_angle(r.angle, inv.angle) + x = np.dot(r.get_axis(convention=convention), inv.get_axis(convention=convention)) + assert isclose(-1, x, rel_tol=0, abs_tol=1e-15) + + +@pytest.mark.parametrize('matrix, exception', [ + (np.array([ + [0.0, 1.0, 0.0], + [1.0, 0.0, 0.0] + ]), TypeError), + (np.array([ + [0.445888, 0.797184, -0.407040], + [0.821760, -0.184320, 0.539200], + [-0.354816, 0.574912, 0.737280] + ]), ValueError) +]) +def test_matrix_error(matrix, exception): + with pytest.raises(exception): + Rotation.from_matrix(matrix) + + +@pytest.mark.parametrize('matrix, q0, q1, q2, q3', [ + ( + np.array([ + [0.445888, 0.797184, -0.407040], + [-0.354816, 0.574912, 0.737280], + [0.821760, -0.184320, 0.539200] + ]), + 0.8, 0.288, 0.384, 0.36 + ), + ( + np.array([ + [0.539200, 0.737280, 0.407040], + [0.184320, -0.574912, 0.797184], + [0.821760, -0.354816, -0.445888] + ]), + 0.36, 0.8, 0.288, 0.384 + ), + ( + np.array([ + [-0.445888, 0.797184, -0.407040], + [0.354816, 0.574912, 0.737280], + [0.821760, 0.184320, -0.539200] + ]), + 0.384, 0.36, 0.8, 0.288 + ), + ( + np.array([ + [-0.539200, 0.737280, 0.407040], + [-0.184320, -0.574912, 0.797184], + [0.821760, 0.354816, 0.445888] + ]), + 0.288, 0.384, 0.36, 0.8 + ) +]) +def test_matrix(matrix, q0, q1, q2, q3): + r = Rotation.from_matrix(matrix) + assert isclose(r.q0, q0, rel_tol=0, abs_tol=1e-15) + assert isclose(r.q1, q1, rel_tol=0, abs_tol=1e-15) + assert isclose(r.q2, q2, rel_tol=0, abs_tol=1e-15) + assert isclose(r.q3, q3, rel_tol=0, abs_tol=1e-15) + + +def test_matrix_other(): + m1 = np.array([ + [0, 1, 0], + [0, 0, 1], + [1, 0, 0], + ]) + r = Rotation.from_matrix(m1) + check_vector(r.apply_to(I), K) + check_vector(r.apply_to(J), I) + check_vector(r.apply_to(K), J) + + m2 = np.array([ + [0.83203, -0.55012, -0.07139], + [0.48293, 0.78164, -0.39474], + [0.27296, 0.29396, 0.91602] + ]) + r = Rotation.from_matrix(m2) + m3 = r.matrix + assert np.allclose(m2, m3, rtol=0, atol=6e-6) + + for i in range(3): + for j in range(3): + m3tm3 = m3[i][0] * m3[j][0] + m3[i][1] * m3[j][1] + m3[i][2] * m3[j][2] + if i == j: + assert isclose(m3tm3 - 1, 0, rel_tol=0, abs_tol=1e-14) + else: + assert isclose(m3tm3, 0, rel_tol=0, abs_tol=1e-15) + + +@pytest.mark.parametrize('convention', ['vector', 'frame']) +@pytest.mark.parametrize('order', [ + RotationOrder.XZX, + RotationOrder.YXY, + RotationOrder.ZYZ, + RotationOrder.XYX, + RotationOrder.YZY, + RotationOrder.ZXZ, +]) +def test_euler_angles_proper_euler(convention, order): + it = itertools.product( + np.arange(0.1, 6.2, 0.3) * u.rad, + np.arange(0.05, 3.1, 0.3) * u.rad, + np.arange(0.1, 6.2, 0.3) * u.rad, + ) + for alpha1, alpha2, alpha3 in it: + r = Rotation.from_euler_angles( + order, alpha1, alpha2, alpha3, + convention=convention) + + angles = r.get_angles(order, convention=convention) + check_angle(angles[0], alpha1) + check_angle(angles[1], alpha2) + check_angle(angles[2], alpha3) + + +@pytest.mark.parametrize('convention', ['vector', 'frame']) +@pytest.mark.parametrize('order', [ + RotationOrder.XYZ, + RotationOrder.XZY, + RotationOrder.YXZ, + RotationOrder.YZX, + RotationOrder.ZXY, + RotationOrder.ZYX, +]) +def test_euler_angles_cardan(convention, order): + it = itertools.product( + np.arange(0.1, 6.2, 0.3) * u.rad, + np.arange(-1.55, 1.55, 0.3) * u.rad, + np.arange(0.1, 6.2, 0.3) * u.rad, + ) + for alpha1, alpha2, alpha3 in it: + r = Rotation.from_euler_angles( + order, alpha1, alpha2, alpha3, + convention=convention) + + angles = r.get_angles(order, convention=convention) + check_angle(angles[0], alpha1) + check_angle(angles[1], alpha2) + check_angle(angles[2], alpha3) + + +@pytest.mark.parametrize('convention', ['vector', 'frame']) +def test_compose(convention): + r1 = Rotation.from_axis_angle( + np.array([2, -3, 5]), 1.7 * u.rad, convention=convention) + + r2 = Rotation.from_axis_angle( + np.array([-1, 3, 2]), 0.3 * u.rad, convention=convention) + + r3 = r2.compose(r1, convention=convention) + + float_range = np.arange(-0.9, 0.9, 0.2) + + it = itertools.product(float_range, float_range, float_range) + + for x, y, z in it: + w = np.array([x, y, z]) + + if convention == 'vector': + v1 = r2.apply_to(r1.apply_to(w)) + v2 = r3.apply_to(w) + elif convention == 'frame': + v1 = r1.apply_to(r2.apply_to(w)) + v2 = r3.apply_to(w) + + check_vector(v1, v2) diff --git a/tox.ini b/tox.ini index 3867ca8..b7c886b 100644 --- a/tox.ini +++ b/tox.ini @@ -44,6 +44,7 @@ deps= sphinxcontrib-spelling commands= sphinx-build -W -b html -d {envtmpdir}/doctrees docs docs/_build/html + sphinx-build -W -b doctest docs docs/_build/html sphinx-build -W -b spelling docs docs/_build/html doc8 docs/