summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorlim <lim@stsci.edu>2012-03-29 11:31:13 -0400
committerlim <lim@stsci.edu>2012-03-29 11:31:13 -0400
commitaf99439c057871ea8974f5e4cdd674b8c6110dfe (patch)
treef82170a77c8cb63691f1f235e46ee14664895eb7
parent69e439b92cece8d8a2def7d38c2cb38b639eab0d (diff)
downloadstsci.sphere-af99439c057871ea8974f5e4cdd674b8c6110dfe.tar.gz
lim added files to sphere branch
git-svn-id: http://svn.stsci.edu/svn/ssb/stsci_python/stsci_python/branches/sphere@15878 fe389314-cf27-0410-b35b-8c050e845b92 Former-commit-id: b51412b5b545598ae101152eeaed470b08616176
-rw-r--r--doc/Makefile153
-rw-r--r--doc/make.bat190
-rw-r--r--doc/source/api.rst30
-rw-r--r--doc/source/conf.py246
-rw-r--r--doc/source/cutlines.pngbin0 -> 54023 bytes
-rw-r--r--doc/source/index.rst23
-rw-r--r--doc/source/inside.pngbin0 -> 86567 bytes
-rw-r--r--doc/source/user.rst160
-rw-r--r--lib/__init__.py0
-rw-r--r--lib/graph.py792
-rw-r--r--lib/great_circle_arc.py318
-rw-r--r--lib/polygon.py679
-rw-r--r--lib/skyline.py48
-rw-r--r--lib/test/__init__.py1
-rw-r--r--lib/test/benchmarks.py39
-rw-r--r--lib/test/data/1904-66_TAN.fits375
-rw-r--r--lib/test/data/2chipA.fits.REMOVED.git-id1
-rw-r--r--lib/test/data/2chipB.fits.REMOVED.git-id1
-rw-r--r--lib/test/data/2chipC.fits.REMOVED.git-id1
-rw-r--r--lib/test/test.py260
-rw-r--r--lib/test/test_intersection.py167
-rw-r--r--lib/test/test_union.py169
-rw-r--r--lib/test/test_util.py14
-rw-r--r--lib/vector.py213
-rwxr-xr-xsetup.py45
-rw-r--r--src/math_util.c374
26 files changed, 4299 insertions, 0 deletions
diff --git a/doc/Makefile b/doc/Makefile
new file mode 100644
index 0000000..6dd7fc4
--- /dev/null
+++ b/doc/Makefile
@@ -0,0 +1,153 @@
+# Makefile for Sphinx documentation
+#
+
+# You can set these variables from the command line.
+SPHINXOPTS =
+SPHINXBUILD = sphinx-build
+PAPER =
+BUILDDIR = build
+
+# Internal variables.
+PAPEROPT_a4 = -D latex_paper_size=a4
+PAPEROPT_letter = -D latex_paper_size=letter
+ALLSPHINXOPTS = -d $(BUILDDIR)/doctrees $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) source
+# the i18n builder cannot share the environment and doctrees with the others
+I18NSPHINXOPTS = $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) source
+
+.PHONY: help clean html dirhtml singlehtml pickle json htmlhelp qthelp devhelp epub latex latexpdf text man changes linkcheck doctest gettext
+
+help:
+ @echo "Please use \`make <target>' where <target> is one of"
+ @echo " html to make standalone HTML files"
+ @echo " dirhtml to make HTML files named index.html in directories"
+ @echo " singlehtml to make a single large HTML file"
+ @echo " pickle to make pickle files"
+ @echo " json to make JSON files"
+ @echo " htmlhelp to make HTML files and a HTML help project"
+ @echo " qthelp to make HTML files and a qthelp project"
+ @echo " devhelp to make HTML files and a Devhelp project"
+ @echo " epub to make an epub"
+ @echo " latex to make LaTeX files, you can set PAPER=a4 or PAPER=letter"
+ @echo " latexpdf to make LaTeX files and run them through pdflatex"
+ @echo " text to make text files"
+ @echo " man to make manual pages"
+ @echo " texinfo to make Texinfo files"
+ @echo " info to make Texinfo files and run them through makeinfo"
+ @echo " gettext to make PO message catalogs"
+ @echo " changes to make an overview of all changed/added/deprecated items"
+ @echo " linkcheck to check all external links for integrity"
+ @echo " doctest to run all doctests embedded in the documentation (if enabled)"
+
+clean:
+ -rm -rf $(BUILDDIR)/*
+
+html:
+ $(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html
+ @echo
+ @echo "Build finished. The HTML pages are in $(BUILDDIR)/html."
+
+dirhtml:
+ $(SPHINXBUILD) -b dirhtml $(ALLSPHINXOPTS) $(BUILDDIR)/dirhtml
+ @echo
+ @echo "Build finished. The HTML pages are in $(BUILDDIR)/dirhtml."
+
+singlehtml:
+ $(SPHINXBUILD) -b singlehtml $(ALLSPHINXOPTS) $(BUILDDIR)/singlehtml
+ @echo
+ @echo "Build finished. The HTML page is in $(BUILDDIR)/singlehtml."
+
+pickle:
+ $(SPHINXBUILD) -b pickle $(ALLSPHINXOPTS) $(BUILDDIR)/pickle
+ @echo
+ @echo "Build finished; now you can process the pickle files."
+
+json:
+ $(SPHINXBUILD) -b json $(ALLSPHINXOPTS) $(BUILDDIR)/json
+ @echo
+ @echo "Build finished; now you can process the JSON files."
+
+htmlhelp:
+ $(SPHINXBUILD) -b htmlhelp $(ALLSPHINXOPTS) $(BUILDDIR)/htmlhelp
+ @echo
+ @echo "Build finished; now you can run HTML Help Workshop with the" \
+ ".hhp project file in $(BUILDDIR)/htmlhelp."
+
+qthelp:
+ $(SPHINXBUILD) -b qthelp $(ALLSPHINXOPTS) $(BUILDDIR)/qthelp
+ @echo
+ @echo "Build finished; now you can run "qcollectiongenerator" with the" \
+ ".qhcp project file in $(BUILDDIR)/qthelp, like this:"
+ @echo "# qcollectiongenerator $(BUILDDIR)/qthelp/SphericalGeometryToolkit.qhcp"
+ @echo "To view the help file:"
+ @echo "# assistant -collectionFile $(BUILDDIR)/qthelp/SphericalGeometryToolkit.qhc"
+
+devhelp:
+ $(SPHINXBUILD) -b devhelp $(ALLSPHINXOPTS) $(BUILDDIR)/devhelp
+ @echo
+ @echo "Build finished."
+ @echo "To view the help file:"
+ @echo "# mkdir -p $$HOME/.local/share/devhelp/SphericalGeometryToolkit"
+ @echo "# ln -s $(BUILDDIR)/devhelp $$HOME/.local/share/devhelp/SphericalGeometryToolkit"
+ @echo "# devhelp"
+
+epub:
+ $(SPHINXBUILD) -b epub $(ALLSPHINXOPTS) $(BUILDDIR)/epub
+ @echo
+ @echo "Build finished. The epub file is in $(BUILDDIR)/epub."
+
+latex:
+ $(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex
+ @echo
+ @echo "Build finished; the LaTeX files are in $(BUILDDIR)/latex."
+ @echo "Run \`make' in that directory to run these through (pdf)latex" \
+ "(use \`make latexpdf' here to do that automatically)."
+
+latexpdf:
+ $(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex
+ @echo "Running LaTeX files through pdflatex..."
+ make -C $(BUILDDIR)/latex all-pdf
+ @echo "pdflatex finished; the PDF files are in $(BUILDDIR)/latex."
+
+text:
+ $(SPHINXBUILD) -b text $(ALLSPHINXOPTS) $(BUILDDIR)/text
+ @echo
+ @echo "Build finished. The text files are in $(BUILDDIR)/text."
+
+man:
+ $(SPHINXBUILD) -b man $(ALLSPHINXOPTS) $(BUILDDIR)/man
+ @echo
+ @echo "Build finished. The manual pages are in $(BUILDDIR)/man."
+
+texinfo:
+ $(SPHINXBUILD) -b texinfo $(ALLSPHINXOPTS) $(BUILDDIR)/texinfo
+ @echo
+ @echo "Build finished. The Texinfo files are in $(BUILDDIR)/texinfo."
+ @echo "Run \`make' in that directory to run these through makeinfo" \
+ "(use \`make info' here to do that automatically)."
+
+info:
+ $(SPHINXBUILD) -b texinfo $(ALLSPHINXOPTS) $(BUILDDIR)/texinfo
+ @echo "Running Texinfo files through makeinfo..."
+ make -C $(BUILDDIR)/texinfo info
+ @echo "makeinfo finished; the Info files are in $(BUILDDIR)/texinfo."
+
+gettext:
+ $(SPHINXBUILD) -b gettext $(I18NSPHINXOPTS) $(BUILDDIR)/locale
+ @echo
+ @echo "Build finished. The message catalogs are in $(BUILDDIR)/locale."
+
+changes:
+ $(SPHINXBUILD) -b changes $(ALLSPHINXOPTS) $(BUILDDIR)/changes
+ @echo
+ @echo "The overview file is in $(BUILDDIR)/changes."
+
+linkcheck:
+ $(SPHINXBUILD) -b linkcheck $(ALLSPHINXOPTS) $(BUILDDIR)/linkcheck
+ @echo
+ @echo "Link check complete; look for any errors in the above output " \
+ "or in $(BUILDDIR)/linkcheck/output.txt."
+
+doctest:
+ $(SPHINXBUILD) -b doctest $(ALLSPHINXOPTS) $(BUILDDIR)/doctest
+ @echo "Testing of doctests in the sources finished, look at the " \
+ "results in $(BUILDDIR)/doctest/output.txt."
diff --git a/doc/make.bat b/doc/make.bat
new file mode 100644
index 0000000..d9ed58d
--- /dev/null
+++ b/doc/make.bat
@@ -0,0 +1,190 @@
+@ECHO OFF
+
+REM Command file for Sphinx documentation
+
+if "%SPHINXBUILD%" == "" (
+ set SPHINXBUILD=sphinx-build
+)
+set BUILDDIR=build
+set ALLSPHINXOPTS=-d %BUILDDIR%/doctrees %SPHINXOPTS% source
+set I18NSPHINXOPTS=%SPHINXOPTS% source
+if NOT "%PAPER%" == "" (
+ set ALLSPHINXOPTS=-D latex_paper_size=%PAPER% %ALLSPHINXOPTS%
+ set I18NSPHINXOPTS=-D latex_paper_size=%PAPER% %I18NSPHINXOPTS%
+)
+
+if "%1" == "" goto help
+
+if "%1" == "help" (
+ :help
+ echo.Please use `make ^<target^>` where ^<target^> is one of
+ echo. html to make standalone HTML files
+ echo. dirhtml to make HTML files named index.html in directories
+ echo. singlehtml to make a single large HTML file
+ echo. pickle to make pickle files
+ echo. json to make JSON files
+ echo. htmlhelp to make HTML files and a HTML help project
+ echo. qthelp to make HTML files and a qthelp project
+ echo. devhelp to make HTML files and a Devhelp project
+ echo. epub to make an epub
+ echo. latex to make LaTeX files, you can set PAPER=a4 or PAPER=letter
+ echo. text to make text files
+ echo. man to make manual pages
+ echo. texinfo to make Texinfo files
+ echo. gettext to make PO message catalogs
+ echo. changes to make an overview over all changed/added/deprecated items
+ echo. linkcheck to check all external links for integrity
+ echo. doctest to run all doctests embedded in the documentation if enabled
+ goto end
+)
+
+if "%1" == "clean" (
+ for /d %%i in (%BUILDDIR%\*) do rmdir /q /s %%i
+ del /q /s %BUILDDIR%\*
+ goto end
+)
+
+if "%1" == "html" (
+ %SPHINXBUILD% -b html %ALLSPHINXOPTS% %BUILDDIR%/html
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Build finished. The HTML pages are in %BUILDDIR%/html.
+ goto end
+)
+
+if "%1" == "dirhtml" (
+ %SPHINXBUILD% -b dirhtml %ALLSPHINXOPTS% %BUILDDIR%/dirhtml
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Build finished. The HTML pages are in %BUILDDIR%/dirhtml.
+ goto end
+)
+
+if "%1" == "singlehtml" (
+ %SPHINXBUILD% -b singlehtml %ALLSPHINXOPTS% %BUILDDIR%/singlehtml
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Build finished. The HTML pages are in %BUILDDIR%/singlehtml.
+ goto end
+)
+
+if "%1" == "pickle" (
+ %SPHINXBUILD% -b pickle %ALLSPHINXOPTS% %BUILDDIR%/pickle
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Build finished; now you can process the pickle files.
+ goto end
+)
+
+if "%1" == "json" (
+ %SPHINXBUILD% -b json %ALLSPHINXOPTS% %BUILDDIR%/json
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Build finished; now you can process the JSON files.
+ goto end
+)
+
+if "%1" == "htmlhelp" (
+ %SPHINXBUILD% -b htmlhelp %ALLSPHINXOPTS% %BUILDDIR%/htmlhelp
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Build finished; now you can run HTML Help Workshop with the ^
+.hhp project file in %BUILDDIR%/htmlhelp.
+ goto end
+)
+
+if "%1" == "qthelp" (
+ %SPHINXBUILD% -b qthelp %ALLSPHINXOPTS% %BUILDDIR%/qthelp
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Build finished; now you can run "qcollectiongenerator" with the ^
+.qhcp project file in %BUILDDIR%/qthelp, like this:
+ echo.^> qcollectiongenerator %BUILDDIR%\qthelp\SphericalGeometryToolkit.qhcp
+ echo.To view the help file:
+ echo.^> assistant -collectionFile %BUILDDIR%\qthelp\SphericalGeometryToolkit.ghc
+ goto end
+)
+
+if "%1" == "devhelp" (
+ %SPHINXBUILD% -b devhelp %ALLSPHINXOPTS% %BUILDDIR%/devhelp
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Build finished.
+ goto end
+)
+
+if "%1" == "epub" (
+ %SPHINXBUILD% -b epub %ALLSPHINXOPTS% %BUILDDIR%/epub
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Build finished. The epub file is in %BUILDDIR%/epub.
+ goto end
+)
+
+if "%1" == "latex" (
+ %SPHINXBUILD% -b latex %ALLSPHINXOPTS% %BUILDDIR%/latex
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Build finished; the LaTeX files are in %BUILDDIR%/latex.
+ goto end
+)
+
+if "%1" == "text" (
+ %SPHINXBUILD% -b text %ALLSPHINXOPTS% %BUILDDIR%/text
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Build finished. The text files are in %BUILDDIR%/text.
+ goto end
+)
+
+if "%1" == "man" (
+ %SPHINXBUILD% -b man %ALLSPHINXOPTS% %BUILDDIR%/man
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Build finished. The manual pages are in %BUILDDIR%/man.
+ goto end
+)
+
+if "%1" == "texinfo" (
+ %SPHINXBUILD% -b texinfo %ALLSPHINXOPTS% %BUILDDIR%/texinfo
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Build finished. The Texinfo files are in %BUILDDIR%/texinfo.
+ goto end
+)
+
+if "%1" == "gettext" (
+ %SPHINXBUILD% -b gettext %I18NSPHINXOPTS% %BUILDDIR%/locale
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Build finished. The message catalogs are in %BUILDDIR%/locale.
+ goto end
+)
+
+if "%1" == "changes" (
+ %SPHINXBUILD% -b changes %ALLSPHINXOPTS% %BUILDDIR%/changes
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.The overview file is in %BUILDDIR%/changes.
+ goto end
+)
+
+if "%1" == "linkcheck" (
+ %SPHINXBUILD% -b linkcheck %ALLSPHINXOPTS% %BUILDDIR%/linkcheck
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Link check complete; look for any errors in the above output ^
+or in %BUILDDIR%/linkcheck/output.txt.
+ goto end
+)
+
+if "%1" == "doctest" (
+ %SPHINXBUILD% -b doctest %ALLSPHINXOPTS% %BUILDDIR%/doctest
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Testing of doctests in the sources finished, look at the ^
+results in %BUILDDIR%/doctest/output.txt.
+ goto end
+)
+
+:end
diff --git a/doc/source/api.rst b/doc/source/api.rst
new file mode 100644
index 0000000..5cce592
--- /dev/null
+++ b/doc/source/api.rst
@@ -0,0 +1,30 @@
+API Documentation
+=================
+
+Vectors
+-------
+
+.. automodule:: vector
+ :members:
+
+Great circle arcs
+-----------------
+
+.. automodule:: great_circle_arc
+ :members:
+
+Spherical polygons
+------------------
+
+.. currentmodule:: polygon
+
+.. automodule:: polygon
+ :members:
+
+Graph operations on polygons
+----------------------------
+
+.. currentmodule:: graph
+
+.. automodule:: graph
+ :members:
diff --git a/doc/source/conf.py b/doc/source/conf.py
new file mode 100644
index 0000000..5ee1814
--- /dev/null
+++ b/doc/source/conf.py
@@ -0,0 +1,246 @@
+# -*- coding: utf-8 -*-
+#
+# Spherical Geometry Toolkit documentation build configuration file, created by
+# sphinx-quickstart on Tue Oct 11 09:04:48 2011.
+#
+# This file is execfile()d with the current directory set to its containing dir.
+#
+# Note that not all possible configuration values are present in this
+# autogenerated file.
+#
+# All configuration values have a default; values that are commented out
+# serve to show the default.
+
+from stsci.sphinxext.conf import *
+
+import sys, os
+
+# If extensions (or modules to document with autodoc) are in another directory,
+# add these directories to sys.path here. If the directory is relative to the
+# documentation root, use os.path.abspath to make it absolute, like shown here.
+#sys.path.insert(0, os.path.abspath('.'))
+
+# -- General configuration -----------------------------------------------------
+
+# If your documentation needs a minimal Sphinx version, state it here.
+#needs_sphinx = '1.0'
+
+# Add any Sphinx extension module names here, as strings. They can be extensions
+# coming with Sphinx (named 'sphinx.ext.*') or your custom ones.
+#extensions += ['sphinx.ext.autodoc', 'sphinx.ext.intersphinx', 'sphinx.ext.mathjax', 'sphinx.ext.viewcode']
+extensions += ['sphinx.ext.autodoc', 'sphinx.ext.intersphinx', 'sphinx.ext.pngmath', 'sphinx.ext.viewcode']
+
+# Add any paths that contain templates here, relative to this directory.
+templates_path = ['_templates']
+
+# The suffix of source filenames.
+source_suffix = '.rst'
+
+# The encoding of source files.
+#source_encoding = 'utf-8-sig'
+
+# The master toctree document.
+master_doc = 'index'
+
+# General information about the project.
+project = u'Spherical Geometry Toolkit'
+copyright = u'2011, Michael Droettboom, STScI'
+
+# The version info for the project you're documenting, acts as replacement for
+# |version| and |release|, also used in various other places throughout the
+# built documents.
+#
+# The short X.Y version.
+version = '0.1'
+# The full version, including alpha/beta/rc tags.
+release = '0.1'
+
+# The language for content autogenerated by Sphinx. Refer to documentation
+# for a list of supported languages.
+#language = None
+
+# There are two options for replacing |today|: either, you set today to some
+# non-false value, then it is used:
+#today = ''
+# Else, today_fmt is used as the format for a strftime call.
+#today_fmt = '%B %d, %Y'
+
+# List of patterns, relative to source directory, that match files and
+# directories to ignore when looking for source files.
+exclude_patterns = []
+
+# The reST default role (used for this markup: `text`) to use for all documents.
+#default_role = None
+
+# If true, '()' will be appended to :func: etc. cross-reference text.
+#add_function_parentheses = True
+
+# If true, the current module name will be prepended to all description
+# unit titles (such as .. function::).
+#add_module_names = True
+
+# If true, sectionauthor and moduleauthor directives will be shown in the
+# output. They are ignored by default.
+#show_authors = False
+
+# The name of the Pygments (syntax highlighting) style to use.
+pygments_style = 'sphinx'
+
+# A list of ignored prefixes for module index sorting.
+#modindex_common_prefix = []
+
+
+# -- Options for HTML output ---------------------------------------------------
+
+# The theme to use for HTML and HTML Help pages. See the documentation for
+# a list of builtin themes.
+#html_theme = 'default'
+
+# Theme options are theme-specific and customize the look and feel of a theme
+# further. For a list of options available for each theme, see the
+# documentation.
+#html_theme_options = {}
+
+# Add any paths that contain custom themes here, relative to this directory.
+#html_theme_path = []
+
+# The name for this set of Sphinx documents. If None, it defaults to
+# "<project> v<release> documentation".
+#html_title = None
+
+# A shorter title for the navigation bar. Default is the same as html_title.
+#html_short_title = None
+
+# The name of an image file (relative to this directory) to place at the top
+# of the sidebar.
+#html_logo = None
+
+# The name of an image file (within the static path) to use as favicon of the
+# docs. This file should be a Windows icon file (.ico) being 16x16 or 32x32
+# pixels large.
+#html_favicon = None
+
+# Add any paths that contain custom static files (such as style sheets) here,
+# relative to this directory. They are copied after the builtin static files,
+# so a file named "default.css" will overwrite the builtin "default.css".
+html_static_path = ['_static']
+
+# If not '', a 'Last updated on:' timestamp is inserted at every page bottom,
+# using the given strftime format.
+#html_last_updated_fmt = '%b %d, %Y'
+
+# If true, SmartyPants will be used to convert quotes and dashes to
+# typographically correct entities.
+#html_use_smartypants = True
+
+# Custom sidebar templates, maps document names to template names.
+#html_sidebars = {}
+
+# Additional templates that should be rendered to pages, maps page names to
+# template names.
+#html_additional_pages = {}
+
+# If false, no module index is generated.
+#html_domain_indices = True
+
+# If false, no index is generated.
+#html_use_index = True
+
+# If true, the index is split into individual pages for each letter.
+#html_split_index = False
+
+# If true, links to the reST sources are added to the pages.
+#html_show_sourcelink = True
+
+# If true, "Created using Sphinx" is shown in the HTML footer. Default is True.
+#html_show_sphinx = True
+
+# If true, "(C) Copyright ..." is shown in the HTML footer. Default is True.
+#html_show_copyright = True
+
+# If true, an OpenSearch description file will be output, and all pages will
+# contain a <link> tag referring to it. The value of this option must be the
+# base URL from which the finished HTML is served.
+#html_use_opensearch = ''
+
+# This is the file name suffix for HTML files (e.g. ".xhtml").
+#html_file_suffix = None
+
+# Output file base name for HTML help builder.
+htmlhelp_basename = 'SphericalGeometryToolkitdoc'
+
+
+# -- Options for LaTeX output --------------------------------------------------
+
+# The paper size ('letter' or 'a4').
+#latex_paper_size = 'letter'
+
+# The font size ('10pt', '11pt' or '12pt').
+#latex_font_size = '10pt'
+
+# Grouping the document tree into LaTeX files. List of tuples
+# (source start file, target name, title, author, documentclass [howto/manual]).
+latex_documents = [
+ ('index', 'SphericalGeometryToolkit.tex', u'Spherical Geometry Toolkit Documentation',
+ u'Michael Droettboom, STScI', 'manual'),
+]
+
+# The name of an image file (relative to this directory) to place at the top of
+# the title page.
+#latex_logo = None
+
+# For "manual" documents, if this is true, then toplevel headings are parts,
+# not chapters.
+#latex_use_parts = False
+
+# If true, show page references after internal links.
+#latex_show_pagerefs = False
+
+# If true, show URL addresses after external links.
+#latex_show_urls = False
+
+# Additional stuff for the LaTeX preamble.
+#latex_preamble = ''
+
+# Documents to append as an appendix to all manuals.
+#latex_appendices = []
+
+# If false, no module index is generated.
+#latex_domain_indices = True
+
+
+# -- Options for manual page output --------------------------------------------
+
+# One entry per manual page. List of tuples
+# (source start file, name, description, authors, manual section).
+man_pages = [
+ ('index', 'sphericalgeometrytoolkit', u'Spherical Geometry Toolkit Documentation',
+ [u'Michael Droettboom, STScI'], 1)
+]
+
+# If true, show URL addresses after external links.
+#man_show_urls = False
+
+
+# -- Options for Texinfo output ------------------------------------------------
+
+# Grouping the document tree into Texinfo files. List of tuples
+# (source start file, target name, title, author,
+# dir menu entry, description, category)
+texinfo_documents = [
+ ('index', 'SphericalGeometryToolkit', u'Spherical Geometry Toolkit Documentation', u'Michael Droettboom, STScI',
+ 'SphericalGeometryToolkit', 'One line description of project.', 'Miscellaneous'),
+]
+
+# Documents to append as an appendix to all manuals.
+#texinfo_appendices = []
+
+# If false, no module index is generated.
+#texinfo_domain_indices = True
+
+# How to display URL addresses: 'footnote', 'no', or 'inline'.
+#texinfo_show_urls = 'footnote'
+
+
+# Example configuration for intersphinx: refer to the Python standard library.
+intersphinx_mapping = {'http://docs.python.org/': None}
diff --git a/doc/source/cutlines.png b/doc/source/cutlines.png
new file mode 100644
index 0000000..6766848
--- /dev/null
+++ b/doc/source/cutlines.png
Binary files differ
diff --git a/doc/source/index.rst b/doc/source/index.rst
new file mode 100644
index 0000000..680a41a
--- /dev/null
+++ b/doc/source/index.rst
@@ -0,0 +1,23 @@
+.. Spherical Geometry Toolkit documentation master file, created by
+ sphinx-quickstart on Tue Oct 11 09:04:48 2011.
+ You can adapt this file completely to your liking, but it should at least
+ contain the root `toctree` directive.
+
+Welcome to Spherical Geometry Toolkit's documentation!
+======================================================
+
+Contents:
+
+.. toctree::
+ :maxdepth: 2
+
+ user
+ api
+
+Indices and tables
+==================
+
+* :ref:`genindex`
+* :ref:`modindex`
+* :ref:`search`
+
diff --git a/doc/source/inside.png b/doc/source/inside.png
new file mode 100644
index 0000000..a56246f
--- /dev/null
+++ b/doc/source/inside.png
Binary files differ
diff --git a/doc/source/user.rst b/doc/source/user.rst
new file mode 100644
index 0000000..e17db6f
--- /dev/null
+++ b/doc/source/user.rst
@@ -0,0 +1,160 @@
+User documentation
+==================
+
+.. currentmodule:: sphere
+
+The `sphere` library is a pure Python package for handling spherical
+polygons that represent arbitrary regions of the sky.
+
+Requirements
+------------
+
+- Python 2.7
+
+- Numpy 1.4 or later
+
+Coordinate representation
+-------------------------
+
+Coordinates in world space are traditionally represented by right
+ascension and declination (*ra* and *dec*), or longitude and latitude.
+While these representations are convenient, they have discontinuities
+at the poles, making operations on them trickier at arbitrary
+locations on the sky sphere. Therefore, all internal operations of
+this library are done in 3D vector space, where coordinates are
+represented as (*x*, *y*, *z*) vectors. The `sphere.vector` module
+contains functions to convert between (*ra*, *dec*) and (*x*, *y*,
+*z*) representations.
+
+While any (*x*, *y*, *z*) triple represents a vector and therefore a
+location on the sky sphere, a distinction must be made between
+normalized coordinates fall exactly on the unit sphere, and
+unnormalized coordinates which do not. A normalized coordinate is
+defined as a vector whose length is 1, i.e.:
+
+.. math::
+
+ \sqrt{x^2 + y^2 + z^2} = 1
+
+To prevent unnecessary recomputation, many methods in this library
+assume that the vectors passed in are already normalized. If this is
+not the case, `sphere.vector.normalize_vector` can be used to
+normalize an array of vectors.
+
+The library allows the user to work in either degrees or radians. All
+methods that require or return an angular value have a `degrees`
+keyword argument. When `degrees` is `True`, these measurements are in
+degrees, otherwise they are in radians.
+
+Spherical polygons
+------------------
+
+Spherical polygons are arbitrary areas on the sky sphere enclosed by
+great circle arcs. They are represented by the
+`~sphere.polygon.SphericalPolygon` class.
+
+Representation
+``````````````
+
+The points defining the polygon are available from the
+`~polygon.SphericalPolygon.points` property. It is a Nx3 array where
+each row is an (*x*, *y*, *z*) vector, normalized. The polygon points
+are explicitly closed, i.e., the first and last points are the same.
+
+Where is the inside?
+^^^^^^^^^^^^^^^^^^^^
+
+The edges of a polygon serve to separate the “inside” from the
+“outside” area. On a traditional 2D planar surface, the “inside” is
+defined as the finite area and the “outside” is the infinite area.
+However, since the surface of a sphere is cyclical, i.e., it wraps
+around on itself, the a spherical polygon actually defines two finite
+areas. To specify which should be considered the “inside” vs. the
+“outside”, the definition of the polygon also has an “inside point”
+which is just any point that should be considered inside of the
+polygon.
+
+In the following image, the inside point (marked with the red dot)
+declares that the area of the polygon is the green region, and not the
+white region.
+
+.. image:: inside.png
+
+The inside point of the the polygon can be obtained from the
+`~polygon.SphericalPolygon.inside` property.
+
+Cut lines
+^^^^^^^^^
+
+If the polygon represents two disjoint areas or the polygon has holes,
+those areas will be connected by cut lines. The following image shows
+a polygon made from the union of a number of cone areas which has both
+a hole and a disjoint region connected by cut lines.
+
+.. image:: cutlines.png
+
+Creating spherical polygons
+```````````````````````````
+
+.. currentmodule:: sphere.polygon
+
+`SphericalPolygon` objects have 4 different constructors:
+
+ - `SphericalPolygon`: Takes an array of (*x*, *y*, *z*)
+ points and an inside point.
+
+ - `SphericalPolygon.from_radec`: Takes an array of (*ra*, *dec*)
+ points and an inside point.
+
+ - `SphericalPolygon.from_cone`: Creates a polygon from a cone on the
+ sky shere. Takes (*ra*, *dec*, *radius*).
+
+ - `SphericalPolygon.from_wcs`: Creates a polygon from the footprint
+ of a FITS image using its WCS header keywords. Takes a FITS
+ filename or a `pyfits.Header` object.
+
+Operations on Spherical Polygons
+````````````````````````````````
+
+Once one has a `SphericalPolygon` object, there are a number of
+operations available:
+
+ - `~SphericalPolygon.contains_point`: Determines if the given point is inside the polygon.
+
+ - `~SphericalPolygon.intersects_poly`: Determines if one polygon intersects with another.
+
+ - `~SphericalPolygon.area`: Determine the area of a polygon.
+
+ - `~SphericalPolygon.union` and `~SphericalPolygon.multi_union`:
+ Return a new polygon that is the union of two or more polygons.
+
+ - `~SphericalPolygon.intersection` and
+ `~SphericalPolygon.multi_intersection`: Return a new polygon that
+ is the intersection of two or more polygons.
+
+ - `~SphericalPolygon.overlap`: Determine how much a given polygon
+ overlaps another.
+
+ - `~SphericalPolygon.draw`: Plots the polygon using matplotlib’s
+ Basemap toolkit. This feature is rather bare and intended
+ primarily for debugging purposes.
+
+Great circle arcs
+-----------------
+
+.. currentmodule:: sphere.great_circle_arc
+
+As seen above, great circle arcs are used to define the edges of the
+polygon. The `sphere.great_circle_arc` module contains a number of
+functions that are useful for dealing with them.
+
+- `length`: Returns the angular distance between two points on the sphere.
+
+- `intersection`: Returns the intersection point between two great
+ circle arcs.
+
+- `intersects`: Determines if two great circle arcs intersect.
+
+- `angle`: Calculate the angle between two great circle arcs.
+
+- `midpoint`: Calculate the midpoint along a great circle arc.
diff --git a/lib/__init__.py b/lib/__init__.py
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/lib/__init__.py
diff --git a/lib/graph.py b/lib/graph.py
new file mode 100644
index 0000000..e9f03ba
--- /dev/null
+++ b/lib/graph.py
@@ -0,0 +1,792 @@
+# -*- coding: utf-8 -*-
+
+# Copyright (C) 2011 Association of Universities for Research in
+# Astronomy (AURA)
+#
+# 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 name of AURA and its representatives may not be used to
+# endorse or promote products derived from this software without
+# specific prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY AURA ``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 AURA 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.
+
+"""
+This contains the code that does the actual unioning of regions.
+"""
+# TODO: Weak references for memory management problems?
+
+# STDLIB
+import itertools
+import weakref
+
+# THIRD-PARTY
+import numpy as np
+
+# LOCAL
+#from . import great_circle_arc
+#from . import vector
+import great_circle_arc
+import vector
+
+# Set to True to enable some sanity checks
+DEBUG = False
+
+
+class Graph:
+ """
+ A graph of nodes connected by edges. The graph is used to build
+ unions between polygons.
+
+ .. note::
+ This class is not meant to be used directly. Instead, use
+ `sphere.polygon.SphericalPolygon.union` and
+ `sphere.polygon.SphericalPolygon.intersection`.
+ """
+
+ class Node:
+ """
+ A `Node` represents a single point, connected by an arbitrary
+ number of `Edge` objects to other `Node` objects.
+ """
+ def __init__(self, point, source_polygons=[]):
+ """
+ Parameters
+ ----------
+ point : 3-sequence (*x*, *y*, *z*) coordinate
+
+ source_polygon : `~sphere.polygon.SphericalPolygon` instance, optional
+ The polygon(s) this node came from. Used for bookkeeping.
+ """
+ self._point = np.asanyarray(point)
+ self._source_polygons = set(source_polygons)
+ self._edges = weakref.WeakSet()
+
+ def __repr__(self):
+ return "Node(%s %d)" % (str(self._point), len(self._edges))
+
+ def follow(self, edge):
+ """
+ Follows from one edge to another across this node.
+
+ Parameters
+ ----------
+ edge : `Edge` instance
+ The edge to follow away from.
+
+ Returns
+ -------
+ other : `Edge` instance
+ The other edge.
+ """
+ edges = list(self._edges)
+ assert len(edges) == 2
+ try:
+ return edges[not edges.index(edge)]
+ except IndexError:
+ raise ValueError("Following from disconnected edge")
+
+ def equals(self, other):
+ """
+ Returns `True` if the location of this and the *other*
+ `Node` are the same.
+ """
+ return np.array_equal(self._point, other._point)
+
+
+ class Edge:
+ """
+ An `Edge` represents a connection between exactly two `Node`
+ objects. This `Edge` class has no direction.
+ """
+ def __init__(self, A, B, source_polygons=[]):
+ """
+ Parameters
+ ----------
+ A, B : `Node` instances
+
+ source_polygon : `~sphere.polygon.SphericalPolygon` instance, optional
+ The polygon this edge came from. Used for bookkeeping.
+ """
+ self._nodes = [A, B]
+ for node in self._nodes:
+ node._edges.add(self)
+ self._source_polygons = set(source_polygons)
+
+ def __repr__(self):
+ nodes = self._nodes
+ return "Edge(%s -> %s)" % (nodes[0]._point, nodes[1]._point)
+
+ def follow(self, node):
+ """
+ Follow along the edge from the given *node* to the other
+ node.
+
+ Parameters
+ ----------
+ node : `Node` instance
+
+ Returns
+ -------
+ other : `Node` instance
+ """
+ nodes = self._nodes
+ try:
+ return nodes[not nodes.index(node)]
+ except IndexError:
+ raise ValueError("Following from disconnected node")
+
+
+ def __init__(self, polygons):
+ """
+ Parameters
+ ----------
+ polygons : sequence of `~sphere.polygon.SphericalPolygon` instances
+ Build a graph from this initial set of polygons.
+ """
+ self._nodes = set()
+ self._edges = set()
+ self._source_polygons = set()
+
+ self.add_polygons(polygons)
+
+ def add_polygons(self, polygons):
+ """
+ Add more polygons to the graph.
+
+ .. note::
+ Must be called before `union` or `intersection`.
+
+ Parameters
+ ----------
+ polygons : sequence of `~sphere.polygon.SphericalPolygon` instances
+ Set of polygons to add to the graph
+ """
+ for polygon in polygons:
+ self.add_polygon(polygon)
+
+ def add_polygon(self, polygon):
+ """
+ Add a single polygon to the graph.
+
+ .. note::
+ Must be called before `union` or `intersection`.
+
+ Parameters
+ ----------
+ polygon : `~sphere.polygon.SphericalPolygon` instance
+ Polygon to add to the graph
+ """
+ points = polygon._points
+
+ if len(points) < 3:
+ raise ValueError("Too few points in polygon")
+
+ self._source_polygons.add(polygon)
+
+ start_node = nodeA = self.add_node(points[0], [polygon])
+ for i in range(1, len(points) - 1):
+ # Don't create self-pointing edges
+ if not np.array_equal(points[i], nodeA._point):
+ nodeB = self.add_node(points[i], [polygon])
+ self.add_edge(nodeA, nodeB, [polygon])
+ nodeA = nodeB
+ # Close the polygon
+ self.add_edge(nodeA, start_node, [polygon])
+
+ def add_node(self, point, source_polygons=[]):
+ """
+ Add a node to the graph. It will be disconnected until used
+ in a call to `add_edge`.
+
+ Parameters
+ ----------
+ point : 3-sequence (*x*, *y*, *z*) coordinate
+
+ source_polygon : `~sphere.polygon.SphericalPolygon` instance, optional
+ The polygon this node came from. Used for bookkeeping.
+
+ Returns
+ -------
+ node : `Node` instance
+ The new node
+ """
+ node = self.Node(point, source_polygons)
+ self._nodes.add(node)
+ return node
+
+ def remove_node(self, node):
+ """
+ Removes a node and all of the edges that touch it.
+
+ .. note::
+ It is assumed that *Node* is already a part of the graph.
+
+ Parameters
+ ----------
+ node : `Node` instance
+ """
+ assert node in self._nodes
+
+ for edge in list(node._edges):
+ nodeB = edge.follow(node)
+ nodeB._edges.remove(edge)
+ self._edges.remove(edge)
+ self._nodes.remove(node)
+
+ def add_edge(self, A, B, source_polygons=[]):
+ """
+ Add an edge between two nodes.
+
+ .. note::
+ It is assumed both nodes already belong to the graph.
+
+ Parameters
+ ----------
+ A, B : `Node` instances
+
+ source_polygons : `~sphere.polygon.SphericalPolygon` instance, optional
+ The polygon(s) this edge came from. Used for bookkeeping.
+
+ Returns
+ -------
+ edge : `Edge` instance
+ The new edge
+ """
+ assert A in self._nodes
+ assert B in self._nodes
+
+ edge = self.Edge(A, B, source_polygons)
+ self._edges.add(edge)
+ return edge
+
+ def remove_edge(self, edge):
+ """
+ Remove an edge from the graph. The nodes it points to remain intact.
+
+ .. note::
+ It is assumed that *edge* is already a part of the graph.
+
+ Parameters
+ ----------
+ edge : `Edge` instance
+ """
+ assert edge in self._edges
+
+ A, B = edge._nodes
+ A._edges.remove(edge)
+ if len(A._edges) == 0:
+ self.remove_node(A)
+ if A is not B:
+ B._edges.remove(edge)
+ if len(B._edges) == 0:
+ self.remove_node(B)
+ self._edges.remove(edge)
+
+ def split_edge(self, edge, node):
+ """
+ Splits an `Edge` *edge* at `Node` *node*, removing *edge* and
+ replacing it with two new `Edge` instances. It is intended
+ that *E* is along the original edge, but that is not enforced.
+
+ Parameters
+ ----------
+ edge : `Edge` instance
+ The edge to split
+
+ node : `Node` instance
+ The node to insert
+
+ Returns
+ -------
+ edgeA, edgeB : `Edge` instances
+ The two new edges on either side of *node*.
+ """
+ assert edge in self._edges
+ assert node in self._nodes
+
+ A, B = edge._nodes
+ edgeA = self.add_edge(A, node, edge._source_polygons)
+ edgeB = self.add_edge(node, B, edge._source_polygons)
+ self.remove_edge(edge)
+ return [edgeA, edgeB]
+
+ def _sanity_check(self, node_is_2=False):
+ """
+ For debugging purposes: assert that edges and nodes are
+ connected to each other correctly and there are no orphaned
+ edges or nodes.
+ """
+ if not DEBUG:
+ return
+
+ try:
+ unique_edges = set()
+ for edge in self._edges:
+ for node in edge._nodes:
+ assert edge in node._edges
+ assert node in self._nodes
+ edge_repr = [tuple(x._point) for x in edge._nodes]
+ edge_repr.sort()
+ edge_repr = tuple(edge_repr)
+ # assert edge_repr not in unique_edges
+ unique_edges.add(edge_repr)
+
+ for node in self._nodes:
+ if node_is_2:
+ assert len(node._edges) == 2
+ else:
+ assert len(node._edges) >= 2
+ for edge in node._edges:
+ assert node in edge._nodes
+ assert edge in self._edges
+ except AssertionError as e:
+ import traceback
+ traceback.print_exc()
+ self._dump_graph()
+ raise
+
+ def _dump_graph(self, title=None, lon_0=0, lat_0=90,
+ projection='ortho', func=lambda x: len(x._edges)):
+ from mpl_toolkits.basemap import Basemap
+ from matplotlib import pyplot as plt
+ fig = plt.figure()
+ m = Basemap(projection="ortho",
+ lon_0=lon_0,
+ lat_0=lat_0)
+ m.drawparallels(np.arange(-90., 90., 20.))
+ m.drawmeridians(np.arange(0., 420., 20.))
+ m.drawmapboundary(fill_color='white')
+
+ counts = {}
+ for node in self._nodes:
+ count = func(node)
+ counts.setdefault(count, [])
+ counts[count].append(list(node._point))
+
+ for k, v in counts.items():
+ v = np.array(v)
+ ra, dec = vector.vector_to_radec(v[:, 0], v[:, 1], v[:, 2])
+ x, y = m(ra, dec)
+ m.plot(x, y, 'o', label=str(k))
+
+ for edge in list(self._edges):
+ A, B = [x._point for x in edge._nodes]
+ r0, d0 = vector.vector_to_radec(A[0], A[1], A[2])
+ r1, d1 = vector.vector_to_radec(B[0], B[1], B[2])
+ m.drawgreatcircle(r0, d0, r1, d1, color='blue')
+
+ if title:
+ plt.title(title)
+ plt.legend()
+ plt.show()
+
+ def union(self):
+ """
+ Once all of the polygons have been added to the graph,
+ join the polygons together.
+
+ Returns
+ -------
+ points : Nx3 array of (*x*, *y*, *z*) points
+ This is a list of points outlining the union of the
+ polygons that were given to the constructor. If the
+ original polygons are disjunct or contain holes, cut lines
+ will be included in the output.
+ """
+ self._remove_cut_lines()
+ self._sanity_check()
+ self._find_all_intersections()
+ self._sanity_check()
+ self._remove_interior_nodes()
+ self._sanity_check()
+ self._remove_3ary_edges()
+ self._sanity_check(True)
+ return self._trace()
+
+ def intersection(self):
+ """
+ Once all of the polygons have been added to the graph,
+ calculate the intersection.
+
+ Returns
+ -------
+ points : Nx3 array of (*x*, *y*, *z*) points
+ This is a list of points outlining the intersection of the
+ polygons that were given to the constructor. If the
+ resulting polygons are disjunct or contain holes, cut lines
+ will be included in the output.
+ """
+ self._remove_cut_lines()
+ self._sanity_check()
+ self._find_all_intersections()
+ self._sanity_check()
+ self._remove_exterior_edges()
+ self._sanity_check()
+ self._remove_3ary_edges(large_first=True)
+ self._sanity_check(True)
+ return self._trace()
+
+ def _remove_cut_lines(self):
+ """
+ Removes any cutlines that may already have existed in the
+ input polygons. This is so any cutlines in the final result
+ will be optimized to be as short as possible and won't
+ intersect each other.
+
+ This works by finding coincident edges that are reverse to
+ each other, and then splicing around them.
+ """
+ # As this proceeds, edges are removed from the graph. It
+ # iterates over a static list of all edges that exist at the
+ # start, so each time one is selected, we need to ensure it
+ # still exists as part of the graph.
+
+ # This transforms the following (where = is the cut line)
+ #
+ # \ /
+ # A' + + B'
+ # | |
+ # A +====================+ B
+ #
+ # D +====================+ C
+ # | |
+ # D' + + C'
+ # / \
+ #
+ # to this:
+ #
+ # \ /
+ # A' + + B'
+ # | |
+ # A + + C
+ # | |
+ # D' + + C'
+ # / \
+ #
+
+ edges = list(self._edges)
+ for i in xrange(len(edges)):
+ AB = edges[i]
+ if AB not in self._edges:
+ continue
+ A, B = AB._nodes
+ for j in xrange(i + 1, len(edges)):
+ CD = edges[j]
+ if CD not in self._edges:
+ continue
+ C, D = CD._nodes
+ # To be a cutline, the candidate edges need to run in
+ # the opposite direction, hence A == D and B == C, not
+ # A == C and B == D.
+ if (A.equals(D) and B.equals(C)):
+ # Create new edges A -> D' and C -> B'
+ self.add_edge(
+ A, D.follow(CD).follow(D),
+ AB._source_polygons | CD._source_polygons)
+ self.add_edge(
+ C, B.follow(AB).follow(B),
+ AB._source_polygons | CD._source_polygons)
+
+ # Remove B and D which are identical to C and A
+ # respectively. We do not need to remove AB and
+ # CD because this will remove it for us.
+ self.remove_node(D)
+ self.remove_node(B)
+
+ break
+
+ def _find_all_intersections(self):
+ """
+ Find all the intersecting edges in the graph. For each
+ intersecting pair, four new edges are created around the
+ intersection point.
+ """
+ def get_edge_points(edges):
+ return (np.array([x._nodes[0]._point for x in edges]),
+ np.array([x._nodes[1]._point for x in edges]))
+
+ # For speed, we want to vectorize all of the intersection
+ # calculations. Therefore, there is a list of edges, and two
+ # arrays containing the end points of those edges. They all
+ # need to have things added and removed from them at the same
+ # time to keep them in sync, but of course the interface for
+ # doing so is different between Python lists and numpy arrays.
+
+ edges = list(self._edges)
+ starts, ends = get_edge_points(edges)
+
+ while len(edges) > 1:
+ AB = edges.pop(0)
+ A = starts[0]; starts = starts[1:] # numpy equiv of "pop(0)"
+ B = ends[0]; ends = ends[1:] # numpy equiv of "pop(0)"
+
+ # Calculate the intersection points between AB and all
+ # other remaining edges
+ with np.errstate(invalid='ignore'):
+ intersections = great_circle_arc.intersection(
+ A, B, starts, ends)
+ # intersects is `True` everywhere intersections has an
+ # actual intersection
+ intersects = np.isfinite(intersections[..., 0])
+
+ intersection_indices = np.nonzero(intersects)[0]
+
+ # Iterate through the candidate intersections, if any --
+ # we want to eliminate intersections that only intersect
+ # at the end points
+ for j in intersection_indices:
+ C = starts[j]
+ D = ends[j]
+ if (np.array_equal(C, A) or np.array_equal(C, B) or
+ np.array_equal(D, A) or np.array_equal(D, B)):
+ continue
+
+ CD = edges[j]
+ E = intersections[j]
+
+ # This is a bona-fide intersection, and E is the
+ # point at which the two lines intersect. Make a
+ # new node for it -- this must belong to the all
+ # of the source polygons of both of the edges that
+ # crossed.
+
+ # A
+ # |
+ # C--E--D
+ # |
+ # B
+
+ E = self.add_node(
+ E, AB._source_polygons | CD._source_polygons)
+ newA, newB = self.split_edge(AB, E)
+ newC, newD = self.split_edge(CD, E)
+
+ new_edges = [newA, newB, newC, newD]
+
+ # Delete CD, and push the new edges to the
+ # front so they will be tested for intersection
+ # against all remaining edges.
+ del edges[j] # CD
+ edges = new_edges + edges
+ new_starts, new_ends = get_edge_points(new_edges)
+ starts = np.vstack(
+ (new_starts, starts[:j], starts[j+1:]))
+ ends = np.vstack(
+ (new_ends, ends[:j], ends[j+1:]))
+ break
+
+ def _remove_interior_nodes(self):
+ """
+ Removes any nodes that are contained inside other polygons.
+ What's left is the (possibly disjunct) outline.
+ """
+ polygons = self._source_polygons
+
+ # node._count is the number of polygons that the node is
+ # inside, other than the polygon that it came from
+ for node in self._nodes:
+ node._count = 0
+ for polygon in polygons:
+ if polygon in node._source_polygons:
+ continue
+ if polygon.contains_point(node._point):
+ node._count += 1
+
+ for node in list(self._nodes):
+ # Nodes with a count >= 1 are completely contained in the
+ # outer polygon, so they should be removed. What's left
+ # are only outer nodes.
+ if node._count >= 1:
+ self.remove_node(node)
+
+ def _remove_exterior_edges(self):
+ """
+ Removes any edges that are not contained in all of the source
+ polygons. What's left is the (possibly disjunct) outline.
+ """
+ polygons = self._source_polygons
+
+ for edge in self._edges:
+ edge._count = 0
+ for polygon in polygons:
+ if (polygon in edge._source_polygons or
+ polygon.intersects_arc(
+ edge._nodes[0]._point, edge._nodes[1]._point)):
+ edge._count += 1
+
+ for edge in list(self._edges):
+ if edge._count < len(polygons):
+ self.remove_edge(edge)
+
+ def _remove_3ary_edges(self, large_first=False):
+ """
+ Remove edges between pairs of nodes that have more than 3
+ edges. This removes triangles that can't be traced.
+ """
+ if large_first:
+ max_ary = 0
+ for node in self._nodes:
+ max_ary = max(len(node._edges), max_ary)
+ order = range(max_ary + 1, 2, -1)
+ else:
+ order = [3]
+
+ for i in order:
+ removals = []
+ for edge in list(self._edges):
+ if (len(edge._nodes[0]._edges) >= i and
+ len(edge._nodes[1]._edges) >= i):
+ removals.append(edge)
+
+ for edge in removals:
+ if edge in self._edges:
+ self.remove_edge(edge)
+
+ def _trace(self):
+ """
+ Given a graph that has had cutlines removed and all
+ intersections found, traces it to find a resulting single
+ polygon.
+ """
+ polygons = []
+ edges = set(self._edges) # copy
+ seen_nodes = set()
+ while len(edges):
+ polygon = []
+ edge = edges.pop()
+ start_node = node = edge._nodes[0]
+ while True:
+ # TODO: Do we need this if clause any more?
+ if len(polygon):
+ if not np.array_equal(polygon[-1], node._point):
+ polygon.append(node._point)
+ else:
+ polygon.append(node._point)
+ edge = node.follow(edge)
+ edges.discard(edge)
+ node = edge.follow(node)
+ if node is start_node:
+ polygon.append(node._point)
+ break
+
+ polygons.append(np.asarray(polygon))
+
+ if len(polygons) == 1:
+ return polygons[0]
+ elif len(polygons) == 0:
+ return []
+ else:
+ return self._join(polygons)
+
+ def _join(self, polygons):
+ """
+ If the graph is disjunct, joins the parts with cutlines.
+
+ The closest nodes between each pair that don't intersect
+ any other edges are used as cutlines.
+
+ TODO: This is not optimal, because the closest distance
+ between two polygons may not in fact be between two vertices,
+ but may be somewhere along an edge.
+ """
+ def do_join(polygons):
+ all_polygons = polygons[:]
+
+ skipped = 0
+
+ polyA = polygons.pop(0)
+ while len(polygons):
+ polyB = polygons.pop(0)
+
+ # If fewer than 3 edges, it's not a polygon,
+ # just throw it out
+ if len(polyB) < 4:
+ continue
+
+ # Find the closest set of vertices between polyA and
+ # polyB that don't cross any of the edges in *any* of
+ # the polygons
+ closest = np.inf
+ closest_pair_idx = (None, None)
+ for a in xrange(len(polyA) - 1):
+ A = polyA[a]
+ distances = great_circle_arc.length(A, polyB[:-1])
+ b = np.argmin(distances)
+ distance = distances[b]
+ if distance < closest:
+ B = polyB[b]
+ # Does this candidate line cross other edges?
+ crosses = False
+ for poly in all_polygons:
+ if np.any(
+ great_circle_arc.intersects(
+ A, B, poly[:-1], poly[1:])):
+ crosses = True
+ break
+ if not crosses:
+ closest = distance
+ closest_pair_idx = (a, b)
+
+ if not np.isfinite(closest):
+ # We didn't find a pair of points that don't cross
+ # something else, so we want to try to join another
+ # polygon. Defer the current polygon until later.
+ if len(polygons) in (0, skipped):
+ return None
+ polygons.append(polyB)
+ skipped += 1
+ else:
+ # Splice the two polygons together using a cut
+ # line
+ a, b = closest_pair_idx
+ new_poly = np.vstack((
+ # polyA up to and including the cut point
+ polyA[:a+1],
+
+ # polyB starting with the cut point and
+ # wrapping around back to the cut point.
+ # Ignore the last point in polyB, because it
+ # is the same as the first
+ np.roll(polyB[:-1], -b, 0),
+
+ # The cut point on polyB
+ polyB[b:b+1],
+
+ # the rest of polyA, starting with the cut
+ # point
+ polyA[a:]
+ ))
+
+ skipped = 0
+ polyA = new_poly
+
+ return polyA
+
+ for permutation in itertools.permutations(polygons):
+ poly = do_join(list(permutation))
+ if poly is not None:
+ return poly
+
+ raise RuntimeError("Could not find cut points")
diff --git a/lib/great_circle_arc.py b/lib/great_circle_arc.py
new file mode 100644
index 0000000..0fff174
--- /dev/null
+++ b/lib/great_circle_arc.py
@@ -0,0 +1,318 @@
+# -*- coding: utf-8 -*-
+
+# Copyright (C) 2011 Association of Universities for Research in
+# Astronomy (AURA)
+#
+# 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 name of AURA and its representatives may not be used to
+# endorse or promote products derived from this software without
+# specific prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY AURA ``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 AURA 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 `sphere.great_circle_arc` module contains functions for computing
+the length, intersection, angle and midpoint of great circle arcs.
+
+Great circles are circles on the unit sphere whose center is
+coincident with the center of the sphere. Great circle arcs are the
+section of those circles between two points on the unit sphere.
+"""
+
+from __future__ import with_statement, division, absolute_import
+
+# THIRD-PARTY
+import numpy as np
+
+
+try:
+ #from . import math_util
+ import math_util
+ HAS_C_UFUNCS = True
+except ImportError:
+ HAS_C_UFUNCS = False
+
+
+__all__ = ['angle', 'intersection', 'intersects', 'length', 'midpoint']
+
+
+if not HAS_C_UFUNCS:
+ from numpy.core.umath_tests import inner1d
+
+ def _fast_cross(a, b):
+ """
+ This is a reimplementation of `numpy.cross` that only does 3D x
+ 3D, and is therefore faster since it doesn't need any
+ conditionals.
+ """
+ cp = np.empty(np.broadcast(a, b).shape)
+ aT = a.T
+ bT = b.T
+ cpT = cp.T
+
+ cpT[0] = aT[1]*bT[2] - aT[2]*bT[1]
+ cpT[1] = aT[2]*bT[0] - aT[0]*bT[2]
+ cpT[2] = aT[0]*bT[1] - aT[1]*bT[0]
+
+ return cp
+
+ def _cross_and_normalize(A, B):
+ T = _fast_cross(A, B)
+ # Normalization
+ l = np.sqrt(np.sum(T ** 2, axis=-1))
+ l = np.expand_dims(l, 2)
+ # Might get some divide-by-zeros, but we don't care
+ with np.errstate(invalid='ignore'):
+ TN = T / l
+ return TN
+
+ def intersection(A, B, C, D):
+ ur"""
+ Returns the point of intersection between two great circle arcs.
+ The arcs are defined between the points *AB* and *CD*. Either *A*
+ and *B* or *C* and *D* may be arrays of points, but not both.
+
+ Parameters
+ ----------
+ A, B : (*x*, *y*, *z*) triples or Nx3 arrays of triples
+ Endpoints of the first great circle arc.
+
+ C, D : (*x*, *y*, *z*) triples or Nx3 arrays of triples
+ Endpoints of the second great circle arc.
+
+ Returns
+ -------
+ T : (*x*, *y*, *z*) triples or Nx3 arrays of triples
+ If the given arcs intersect, the intersection is returned. If
+ the arcs do not intersect, the triple is set to all NaNs.
+
+ Notes
+ -----
+ The basic intersection is computed using linear algebra as follows
+ [1]_:
+
+ .. math::
+
+ T = \lVert(A × B) × (C × D)\rVert
+
+ To determine the correct sign (i.e. hemisphere) of the
+ intersection, the following four values are computed:
+
+ .. math::
+
+ s_1 = ((A × B) × A) · T
+
+ s_2 = (B × (A × B)) · T
+
+ s_3 = ((C × D) × C) · T
+
+ s_4 = (D × (C × D)) · T
+
+ For :math:`s_n`, if all positive :math:`T` is returned as-is. If
+ all negative, :math:`T` is multiplied by :math:`-1`. Otherwise
+ the intersection does not exist and is undefined.
+
+ References
+ ----------
+
+ .. [1] Method explained in an `e-mail
+ <http://www.mathworks.com/matlabcentral/newsreader/view_thread/276271>`_
+ by Roger Stafford.
+
+ http://www.mathworks.com/matlabcentral/newsreader/view_thread/276271
+ """
+ A = np.asanyarray(A)
+ B = np.asanyarray(B)
+ C = np.asanyarray(C)
+ D = np.asanyarray(D)
+
+ A, B = np.broadcast_arrays(A, B)
+ C, D = np.broadcast_arrays(C, D)
+
+ ABX = _fast_cross(A, B)
+ CDX = _fast_cross(C, D)
+ T = _cross_and_normalize(ABX, CDX)
+ T_ndim = len(T.shape)
+
+ if T_ndim > 1:
+ s = np.zeros(T.shape[0])
+ else:
+ s = np.zeros(1)
+ s += np.sign(inner1d(_fast_cross(ABX, A), T))
+ s += np.sign(inner1d(_fast_cross(B, ABX), T))
+ s += np.sign(inner1d(_fast_cross(CDX, C), T))
+ s += np.sign(inner1d(_fast_cross(D, CDX), T))
+ if T_ndim > 1:
+ s = np.expand_dims(s, 2)
+
+ cross = np.where(s == -4, -T, np.where(s == 4, T, np.nan))
+
+ # If they share a common point, it's not an intersection. This
+ # gets around some rounding-error/numerical problems with the
+ # above.
+ equals = (np.all(A == C, axis = -1) |
+ np.all(A == D, axis = -1) |
+ np.all(B == C, axis = -1) |
+ np.all(B == D, axis = -1))
+
+ equals = np.expand_dims(equals, 2)
+
+ return np.where(equals, np.nan, cross)
+else:
+ inner1d = math_util.inner1d
+ _fast_cross = math_util.cross
+ _cross_and_normalize = math_util.cross_and_norm
+ intersection = math_util.intersection
+
+
+def length(A, B, degrees=True):
+ ur"""
+ Returns the angular distance between two points (in vector space)
+ on the unit sphere.
+
+ Parameters
+ ----------
+ A, B : (*x*, *y*, *z*) triples or Nx3 arrays of triples
+ The endpoints of the great circle arc, in vector space.
+
+ degrees : bool, optional
+ If `True` (default) the result is returned in decimal degrees,
+ otherwise radians.
+
+ Returns
+ -------
+ length : scalar or array of scalars
+ The angular length of the great circle arc.
+
+ Notes
+ -----
+ The length is computed using the following:
+
+ .. math::
+
+ \Delta = \arccos(A · B)
+ """
+ A = np.asanyarray(A)
+ B = np.asanyarray(B)
+
+ dot = inner1d(A, B)
+ with np.errstate(invalid='ignore'):
+ result = np.arccos(dot)
+
+ if degrees:
+ return np.rad2deg(result)
+ else:
+ return result
+
+
+def intersects(A, B, C, D):
+ """
+ Returns `True` if the great circle arcs between *AB* and *CD*
+ intersect. Either *A* and *B* or *C* and *D* may be arrays of
+ points, but not both.
+
+ Parameters
+ ----------
+ A, B : (*x*, *y*, *z*) triples or Nx3 arrays of triples
+ Endpoints of the first great circle arc.
+
+ C, D : (*x*, *y*, *z*) triples or Nx3 arrays of triples
+ Endpoints of the second great circle arc.
+
+ Returns
+ -------
+ intersects : bool or array of bool
+ If the given arcs intersect, the intersection is returned as
+ `True`.
+ """
+ with np.errstate(invalid='ignore'):
+ intersections = intersection(A, B, C, D)
+ return np.isfinite(intersections[..., 0])
+
+
+def angle(A, B, C, degrees=True):
+ """
+ Returns the angle at *B* between *AB* and *BC*.
+
+ This always returns the shortest angle < π.
+
+ Parameters
+ ----------
+ A, B, C : (*x*, *y*, *z*) triples or Nx3 arrays of triples.
+
+ degrees : bool, optional
+ If `True` (default) the result is returned in decimal degrees,
+ otherwise radians.
+
+ Returns
+ -------
+ angle : float or array of floats
+ The angle at *B* between *AB* and *BC*.
+
+ References
+ ----------
+
+ .. [1] Miller, Robert D. Computing the area of a spherical
+ polygon. Graphics Gems IV. 1994. Academic Press.
+ """
+ A = np.asanyarray(A)
+ B = np.asanyarray(B)
+ C = np.asanyarray(C)
+
+ A, B, C = np.broadcast_arrays(A, B, C)
+
+ ABX = _fast_cross(A, B)
+ ABX = _cross_and_normalize(B, ABX)
+ BCX = _fast_cross(C, B)
+ BCX = _cross_and_normalize(B, BCX)
+ with np.errstate(invalid='ignore'):
+ angle = np.arccos(inner1d(ABX, BCX))
+
+ if degrees:
+ angle = np.rad2deg(angle)
+ return angle
+
+
+def midpoint(A, B):
+ """
+ Returns the midpoint on the great circle arc between *A* and *B*.
+
+ Parameters
+ ----------
+ A, B : (*x*, *y*, *z*) triples or Nx3 arrays of triples
+ The endpoints of the great circle arc. It is assumed that
+ these points are already normalized.
+
+ Returns
+ -------
+ midpoint : (*x*, *y*, *z*) triple or Nx3 arrays of triples
+ The midpoint between *A* and *B*, normalized on the unit
+ sphere.
+ """
+ P = (A + B) / 2.0
+ # Now normalize...
+ l = np.sqrt(np.sum(P * P, axis=-1))
+ l = np.expand_dims(l, 2)
+ return P / l
diff --git a/lib/polygon.py b/lib/polygon.py
new file mode 100644
index 0000000..b9b5d8c
--- /dev/null
+++ b/lib/polygon.py
@@ -0,0 +1,679 @@
+# -*- coding: utf-8 -*-
+
+# Copyright (C) 2011 Association of Universities for Research in
+# Astronomy (AURA)
+#
+# 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 name of AURA and its representatives may not be used to
+# endorse or promote products derived from this software without
+# specific prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY AURA ``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 AURA 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 `polygon` module defines the `SphericalPolygon` class for managing
+polygons on the unit sphere.
+"""
+from __future__ import division, print_function
+
+# STDLIB
+from copy import copy
+
+# THIRD-PARTY
+import numpy as np
+
+# LOCAL
+#from . import great_circle_arc
+#from . import vector
+import great_circle_arc
+import vector
+
+__all__ = ['SphericalPolygon']
+
+
+class SphericalPolygon(object):
+ ur"""
+ Polygons are represented by both a set of points (in
+ Cartesian (*x*, *y*, *z*) normalized on the unit sphere),
+ and an inside point. The inside point is necessary, because
+ both the inside and outside of the polygon are finite areas
+ on the great sphere, and therefore we need a way of
+ specifying which is which.
+ """
+
+ def __init__(self, points, inside):
+ ur"""
+ Parameters
+ ----------
+ points : An Nx3 array of (*x*, *y*, *z*) triples in vector space
+ These points define the boundary of the polygon. It must
+ be "closed", i.e., the last point is the same as the first.
+
+ It may contain zero points, in which it defines the null
+ polygon. It may not contain one, two or three points.
+ Four points are needed to define a triangle, since the
+ polygon must be closed.
+
+ inside : An (*x*, *y*, *z*) triple
+ This point must be inside the polygon.
+ """
+ if len(points) == 0:
+ pass
+ elif len(points) < 3:
+ raise ValueError("Polygon made of too few points")
+ else:
+ assert np.array_equal(points[0], points[-1]), 'Polygon is not closed'
+ self._points = np.asanyarray(points)
+ self._inside = np.asanyarray(inside)
+
+ # TODO: Detect self-intersection and fix
+
+ def __copy__(self):
+ return self.__class__(np.copy(self._points), np.copy(self._inside))
+
+ @property
+ def points(self):
+ """
+ The points defining the polygon. It is an Nx3 array of
+ (*x*, *y*, *z*) vectors. The polygon will be explicitly
+ closed, i.e., the first and last points are the same.
+ """
+ return self._points
+
+ @property
+ def inside(self):
+ """
+ Get the inside point of the polygon.
+ """
+ return self._inside
+
+ @classmethod
+ def from_radec(cls, ra, dec, center=None, degrees=True):
+ ur"""
+ Create a new `SphericalPolygon` from a list of (*ra*, *dec*)
+ points.
+
+ Parameters
+ ----------
+ ra, dec : 1-D arrays of the same length
+ The vertices of the polygon in right-ascension and
+ declination. It must be \"closed\", i.e., that is, the
+ last point is the same as the first.
+
+ center : (*ra*, *dec*) pair, optional
+ A point inside of the polygon to define its inside. If no
+ *center* point is provided, the mean of the polygon's
+ points in vector space will be used. That approach may
+ not work for concave polygons.
+
+ degrees : bool, optional
+ If `True`, (default) *ra* and *dec* are in decimal degrees,
+ otherwise in radians.
+
+ Returns
+ -------
+ polygon : `SphericalPolygon` object
+ """
+ # Convert to Cartesian
+ x, y, z = vector.radec_to_vector(ra, dec, degrees=degrees)
+
+ if center is None:
+ xc = x.mean()
+ yc = y.mean()
+ zc = z.mean()
+ center = vector.normalize_vector(xc, yc, zc)
+ else:
+ center = vector.radec_to_vector(*center, degrees=degrees)
+
+ return cls(np.dstack((x, y, z))[0], center)
+
+ @classmethod
+ def from_cone(cls, ra, dec, radius, degrees=True, steps=16.0):
+ ur"""
+ Create a new `SphericalPolygon` from a cone (otherwise known
+ as a "small circle") defined using (*ra*, *dec*, *radius*).
+
+ The cone is not represented as an ideal circle on the sphere,
+ but as a series of great circle arcs. The resolution of this
+ conversion can be controlled using the *steps* parameter.
+
+ Parameters
+ ----------
+ ra, dec : float scalars
+ This defines the center of the cone
+
+ radius : float scalar
+ The radius of the cone
+
+ degrees : bool, optional
+ If `True`, (default) *ra*, *dec* and *radius* are in
+ decimal degrees, otherwise in radians.
+
+ steps : int, optional
+ The number of steps to use when converting the small
+ circle to a polygon.
+
+ Returns
+ -------
+ polygon : `SphericalPolygon` object
+ """
+ u, v, w = vector.radec_to_vector(ra, dec, degrees=degrees)
+ if degrees:
+ radius = np.deg2rad(radius)
+
+ # Get an arbitrary perpendicular vector. This be be obtained
+ # by crossing (u, v, w) with any unit vector that is not itself.
+ which_min = np.argmin([u, v, w])
+ if which_min == 0:
+ perp = np.cross([u, v, w], [1., 0., 0.])
+ elif which_min == 1:
+ perp = np.cross([u, v, w], [0., 1., 0.])
+ else:
+ perp = np.cross([u, v, w], [0., 0., 1.])
+
+ # Rotate by radius around the perpendicular vector to get the
+ # "pen"
+ x, y, z = vector.rotate_around(
+ u, v, w, perp[0], perp[1], perp[2], radius, degrees=False)
+
+ # Then rotate the pen around the center point all 360 degrees
+ C = np.linspace(0, np.pi * 2.0, steps)
+ # Ensure that the first and last elements are exactly the
+ # same. 2π should equal 0, but with rounding error that isn't
+ # always the case.
+ C[-1] = 0
+ X, Y, Z = vector.rotate_around(x, y, z, u, v, w, C, degrees=False)
+
+ return cls(np.dstack((X, Y, Z))[0], (u, v, w))
+
+ @classmethod
+ def from_wcs(cls, fitspath, steps=1, crval=None):
+ ur"""
+ Create a new `SphericalPolygon` from the footprint of a FITS
+ WCS specification.
+
+ This method requires having `pywcs` and `pyfits` installed.
+
+ Parameters
+ ----------
+ fitspath : path to a FITS file, `pyfits.Header`, or `pywcs.WCS`
+ Refers to a FITS header containing a WCS specification.
+
+ steps : int, optional
+ The number of steps along each edge to convert into
+ polygon edges.
+
+ Returns
+ -------
+ polygon : `SphericalPolygon` object
+ """
+ import pywcs
+ import pyfits
+
+ if isinstance(fitspath, pyfits.Header):
+ header = fitspath
+ wcs = pywcs.WCS(header)
+ elif isinstance(fitspath, pywcs.WCS):
+ wcs = fitspath
+ else:
+ wcs = pywcs.WCS(fitspath)
+ if crval is not None:
+ wcs.wcs.crval = crval
+ xa, ya = [wcs.naxis1, wcs.naxis2]
+
+ length = steps * 4 + 1
+ X = np.empty(length)
+ Y = np.empty(length)
+
+ # Now define each of the 4 edges of the quadrilateral
+ X[0 :steps ] = np.linspace(1, xa, steps, False)
+ Y[0 :steps ] = 1
+ X[steps :steps*2] = xa
+ Y[steps :steps*2] = np.linspace(1, ya, steps, False)
+ X[steps*2:steps*3] = np.linspace(xa, 1, steps, False)
+ Y[steps*2:steps*3] = ya
+ X[steps*3:steps*4] = 1
+ Y[steps*3:steps*4] = np.linspace(ya, 1, steps, False)
+ X[-1] = 1
+ Y[-1] = 1
+
+ # Use wcslib to convert to (ra, dec)
+ ra, dec = wcs.all_pix2sky(X, Y, 1)
+
+ # Convert to Cartesian
+ x, y, z = vector.radec_to_vector(ra, dec)
+
+ # Calculate an inside point
+ ra, dec = wcs.all_pix2sky(xa / 2.0, ya / 2.0, 1)
+ xc, yc, zc = vector.radec_to_vector(ra, dec)
+
+ return cls(np.dstack((x, y, z))[0], (xc[0], yc[0], zc[0]))
+
+ def contains_point(self, point):
+ ur"""
+ Determines if this `SphericalPolygon` contains a given point.
+
+ Parameters
+ ----------
+ point : an (*x*, *y*, *z*) triple
+ The point to test.
+
+ Returns
+ -------
+ contains : bool
+ Returns `True` if the polygon contains the given *point*.
+ """
+ P = self._points
+ r = self._inside
+ point = np.asanyarray(point)
+
+ intersects = great_circle_arc.intersects(P[:-1], P[1:], r, point)
+ crossings = np.sum(intersects)
+
+ return (crossings % 2) == 0
+
+ def intersects_poly(self, other):
+ ur"""
+ Determines if this `SphericalPolygon` intersects another
+ `SphericalPolygon`.
+
+ This method is much faster than actually computing the
+ intersection region between two polygons.
+
+ Parameters
+ ----------
+ other : `SphericalPolygon`
+
+ Returns
+ -------
+ intersects : bool
+ Returns `True` if this polygon intersects the *other*
+ polygon.
+
+ Notes
+ -----
+
+ The algorithm proceeds as follows:
+
+ 1. Determine if any single point of one polygon is contained
+ within the other.
+
+ 2. Deal with the case where only the edges overlap as in::
+
+ : o---------o
+ : o----+---------+----o
+ : | | | |
+ : o----+---------+----o
+ : o---------o
+
+ In this case, an edge from one polygon must cross an
+ edge from the other polygon.
+ """
+ assert isinstance(other, SphericalPolygon)
+
+ # The easy case is in which a point of one polygon is
+ # contained in the other polygon.
+ for point in other._points:
+ if self.contains_point(point):
+ return True
+ for point in self._points:
+ if other.contains_point(point):
+ return True
+
+ # The hard case is when only the edges overlap, as in:
+ #
+ # o---------o
+ # o----+---------+----o
+ # | | | |
+ # o----+---------+----o
+ # o---------o
+ #
+ for i in range(len(self._points) - 1):
+ A = self._points[i]
+ B = self._points[i+1]
+ if np.any(great_circle_arc.intersects(
+ A, B, other._points[:-1], other._points[1:])):
+ return True
+ return False
+
+ def intersects_arc(self, a, b):
+ """
+ Determines if this `SphericalPolygon` intersects or contains
+ the given arc.
+ """
+ return self.contains_point(great_circle_arc.midpoint(a, b))
+
+ if self.contains_point(a) and self.contains_point(b):
+ return True
+
+ P = self._points
+ intersects = great_circle_arc.intersects(P[:-1], P[1:], a, b)
+ if np.any(intersects):
+ return True
+ return False
+
+ def area(self):
+ ur"""
+ Returns the area of the polygon on the unit sphere.
+
+ The algorithm is not able to compute the area of polygons
+ that are larger than half of the sphere. Therefore, the
+ area will always be less than 2π.
+
+ The area is computed based on the sum of the radian angles:
+
+ .. math::
+
+ S = \sum^n_{i=0} \theta_i - (n - 2) \pi
+ """
+ if len(self._points) < 3:
+ return np.array(0.0)
+
+ A = self._points[:-1]
+ B = np.roll(A, 1, 0)
+ C = np.roll(B, 1, 0)
+
+ angles = great_circle_arc.angle(A, B, C, degrees=False)
+
+ midpoints = great_circle_arc.midpoint(A, C)
+ interior = [self.contains_point(x) for x in midpoints]
+ angles = np.where(interior, angles, 2.*np.pi - angles)
+
+ sum = np.sum(angles)
+ area = sum - float(len(angles) - 2) * np.pi
+
+ return area
+
+ def union(self, other):
+ """
+ Return a new `SphericalPolygon` that is the union of *self*
+ and *other*.
+
+ If the polygons are disjoint, they result will be connected
+ using cut lines. For example::
+
+ : o---------o
+ : | |
+ : o---------o=====o----------o
+ : | |
+ : o----------o
+
+ Parameters
+ ----------
+ other : `SphericalPolygon`
+
+ Returns
+ -------
+ polygon : `SphericalPolygon` object
+
+ See also
+ --------
+ multi_union
+
+ Notes
+ -----
+ For implementation details, see the :mod:`~sphere.graph`
+ module.
+ """
+ import graph
+ g = graph.Graph([self, other])
+
+ polygon = g.union()
+
+ return self.__class__(polygon, self._inside)
+
+ @classmethod
+ def multi_union(cls, polygons, method='parallel'):
+ """
+ Return a new `SphericalPolygon` that is the union of all of the
+ polygons in *polygons*.
+
+ Parameters
+ ----------
+ polygons : sequence of `SphericalPolygon`
+
+ method : 'parallel' or 'serial', optional
+ Specifies the method that is used to perform the unions:
+
+ - 'parallel' (default): A graph is built using all of
+ the polygons, and the union operation is computed on
+ the entire thing globally.
+
+ - 'serial': The polygon is built in steps by adding one
+ polygon at a time and unioning at each step.
+
+ This option is provided because one may be faster than the
+ other depending on context, but it primarily exposed for
+ testing reasons. Both modes should theoretically provide
+ equivalent results.
+
+ Returns
+ -------
+ polygon : `SphericalPolygon` object
+
+ See also
+ --------
+ union
+ """
+ assert len(polygons)
+ for polygon in polygons:
+ assert isinstance(polygon, SphericalPolygon)
+
+ import graph
+
+ if method.lower() == 'parallel':
+ g = graph.Graph(polygons)
+ polygon = g.union()
+ return cls(polygon, polygons[0]._inside)
+ elif method.lower() == 'serial':
+ result = copy(polygons[0])
+ for polygon in polygons[1:]:
+ result = result.union(polygon)
+ return result
+ else:
+ raise ValueError("method must be 'parallel' or 'serial'")
+
+ @staticmethod
+ def _find_new_inside(points, polygons):
+ """
+ Finds an acceptable inside point inside of *points* that is
+ also inside of *polygons*. Used by the intersection
+ algorithm, and is really only useful in that context because
+ it requires existing polygons with known inside points.
+ """
+ if len(points) < 4:
+ return np.array([0, 0, 0])
+
+ # Special case for a triangle
+ if len(points) == 4:
+ return np.sum(points[:3]) / 3.0
+
+ for i in range(len(points) - 1):
+ A = points[i]
+ # Skip the adjacent point, since it is by definition on
+ # the edge of the polygon, not potentially running through
+ # the middle.
+ for j in range(i + 2, len(points) - 1):
+ B = points[j]
+ C = great_circle_arc.midpoint(A, B)
+ in_all = True
+ for polygon in polygons:
+ if not polygon.contains_point(C):
+ in_all = False
+ break
+ if in_all:
+ return C
+
+ raise RuntimeError("Suitable inside point could not be found")
+
+ def intersection(self, other):
+ """
+ Return a new `SphericalPolygon` that is the intersection of
+ *self* and *other*.
+
+ If the intersection is empty, a `SphericalPolygon` with zero
+ points will be returned.
+
+ If the result is disjoint, the pieces will be connected using
+ cut lines. For example::
+
+ : o---------o
+ : | |
+ : o---------o=====o----------o
+ : | |
+ : o----------o
+
+ Parameters
+ ----------
+ other : `SphericalPolygon`
+
+ Returns
+ -------
+ polygon : `SphericalPolygon` object
+
+ Notes
+ -----
+ For implementation details, see the :mod:`~sphere.graph`
+ module.
+ """
+ # if not self.intersects_poly(other):
+ # return self.__class__([], [0, 0, 0])
+
+ import graph
+ g = graph.Graph([self, other])
+
+ polygon = g.intersection()
+
+ inside = self._find_new_inside(polygon, [self, other])
+
+ return self.__class__(polygon, inside)
+
+ @classmethod
+ def multi_intersection(cls, polygons, method='parallel'):
+ """
+ Return a new `SphericalPolygon` that is the intersection of
+ all of the polygons in *polygons*.
+
+ Parameters
+ ----------
+ polygons : sequence of `SphericalPolygon`
+
+ method : 'parallel' or 'serial', optional
+ Specifies the method that is used to perform the
+ intersections:
+
+ - 'parallel' (default): A graph is built using all of
+ the polygons, and the intersection operation is computed on
+ the entire thing globally.
+
+ - 'serial': The polygon is built in steps by adding one
+ polygon at a time and computing the intersection at
+ each step.
+
+ This option is provided because one may be faster than the
+ other depending on context, but it primarily exposed for
+ testing reasons. Both modes should theoretically provide
+ equivalent results.
+
+ Returns
+ -------
+ polygon : `SphericalPolygon` object
+ """
+ assert len(polygons)
+ for polygon in polygons:
+ assert isinstance(polygon, SphericalPolygon)
+
+ # for i in range(len(polygons)):
+ # polyA = polygons[i]
+ # for j in range(i + 1, len(polygons)):
+ # polyB = polygons[j]
+ # if not polyA.intersects_poly(polyB):
+ # return cls([], [0, 0, 0])
+
+ import graph
+
+ if method.lower() == 'parallel':
+ g = graph.Graph(polygons)
+ polygon = g.intersection()
+ inside = cls._find_new_inside(polygon, polygons)
+ return cls(polygon, inside)
+ elif method.lower() == 'serial':
+ result = copy(polygons[0])
+ for polygon in polygons[1:]:
+ result = result.intersection(polygon)
+ return result
+ else:
+ raise ValueError("method must be 'parallel' or 'serial'")
+
+ def overlap(self, other):
+ ur"""
+ Returns the fraction of *self* that is overlapped by *other*.
+
+ Let *self* be *a* and *other* be *b*, then the overlap is
+ defined as:
+
+ .. math::
+
+ \frac{S_a}{S_{a \cap b}}
+
+ Parameters
+ ----------
+ other : `SphericalPolygon`
+
+ Returns
+ -------
+ frac : float
+ The fraction of *self* that is overlapped by *other*.
+ """
+ s1 = self.area()
+ intersection = self.intersection(other)
+ s2 = intersection.area()
+ return s2 / s1
+
+ def draw(self, m, **plot_args):
+ """
+ Draws the polygon in a matplotlib.Basemap axes.
+
+ Parameters
+ ----------
+ m : Basemap axes object
+
+ **plot_args : Any plot arguments to pass to basemap
+ """
+ if not len(self._points):
+ return
+ if not len(plot_args):
+ plot_args = {color: 'blue'}
+ points = self._points
+ ra, dec = vector.vector_to_radec(
+ points[:, 0], points[:, 1], points[:, 2],
+ degrees=True)
+ for r0, d0, r1, d1 in zip(ra[0:-1], dec[0:-1], ra[1:], dec[1:]):
+ m.drawgreatcircle(r0, d0, r1, d1, **plot_args)
+ ra, dec = vector.vector_to_radec(
+ *self._inside, degrees=True)
+ x, y = m(ra, dec)
+ m.scatter(x, y, 1, **plot_args)
+
diff --git a/lib/skyline.py b/lib/skyline.py
new file mode 100644
index 0000000..05fd279
--- /dev/null
+++ b/lib/skyline.py
@@ -0,0 +1,48 @@
+"""
+skyline -- Manage outlines on the sky
+
+This module provides support for working with footprints on the sky.
+Primary use case would use the following generalized steps::
+
+1. Initialize SkyLine objects for each input image. This object would be the
+ union of all the input image's individual chips WCS footprints.
+2. Determine overlap between all images. The determination would employ a
+ recursive operation to return the extended list of all overlap values
+ computed as [img1 vs [img2,img3,...,imgN],img2 vs [img3,...,imgN],...]
+3. Select the pair with the largest overlap, or the pair which produces the
+ largest overlap with the first input image. This defines the initial
+ reference SkyLine object.
+4. Perform some operation on the 2 images: for example, match sky in intersecting
+ regions, or aligning second image with the first (reference) image.
+5. Update the second image, either apply the sky value or correct the WCS, then
+ generate a new SkyLine object for that image.
+6. Create a new reference SkyLine object as the union of the initial reference
+ object and the newly updated SkyLine object.
+7. Repeat Steps 2-6 for all remaining input images.
+
+This process will work reasonably fast as most operations are performed using
+the SkyLine objects and WCS information solely, not image data itself.
+"""
+import pyfits
+
+import sphere
+
+class SkyLine(object):
+ def __init__(self,fname):
+ pass
+
+ def add_image(self,skyline):
+ pass
+
+ def compute_overlap(self, skyline):
+ pass
+
+ def find_intersection(self, skyline):
+ pass
+
+ def within(self, pos):
+ pass
+
+ def create_wcs(self):
+ pass
+
diff --git a/lib/test/__init__.py b/lib/test/__init__.py
new file mode 100644
index 0000000..8b13789
--- /dev/null
+++ b/lib/test/__init__.py
@@ -0,0 +1 @@
+
diff --git a/lib/test/benchmarks.py b/lib/test/benchmarks.py
new file mode 100644
index 0000000..9711d7b
--- /dev/null
+++ b/lib/test/benchmarks.py
@@ -0,0 +1,39 @@
+import os
+import sys
+import time
+
+import numpy as np
+from sphere import *
+from test_util import *
+
+def point_in_poly_lots():
+ poly1 = SphericalPolygon.from_wcs(
+ os.path.join(ROOT_DIR, '1904-66_TAN.fits'), 64, crval=[0, 87])
+ poly2 = SphericalPolygon.from_wcs(
+ os.path.join(ROOT_DIR, '1904-66_TAN.fits'), 64, crval=[20, 89])
+ poly3 = SphericalPolygon.from_wcs(
+ os.path.join(ROOT_DIR, '1904-66_TAN.fits'), 64, crval=[180, 89])
+
+ points = get_point_set(density=25)
+
+ count = 0
+ for point in points:
+ if poly1.contains_point(point) or poly2.contains_point(point) or \
+ poly3.contains_point(point):
+ count += 1
+
+ assert count == 5
+ assert poly1.intersects_poly(poly2)
+ assert not poly1.intersects_poly(poly3)
+ assert not poly2.intersects_poly(poly3)
+
+if __name__ == '__main__':
+ for benchmark in [point_in_poly_lots]:
+ t = time.time()
+ sys.stdout.write(benchmark.__name__)
+ sys.stdout.write('...')
+ sys.stdout.flush()
+ benchmark()
+ sys.stdout.write(' %.03fs\n' % (time.time() - t))
+ sys.stdout.flush()
+
diff --git a/lib/test/data/1904-66_TAN.fits b/lib/test/data/1904-66_TAN.fits
new file mode 100644
index 0000000..09c3929
--- /dev/null
+++ b/lib/test/data/1904-66_TAN.fits
@@ -0,0 +1,375 @@
+SIMPLE = T BITPIX = -32 / IEEE (big-endian) 32-bit floating point data NAXIS = 2 NAXIS1 = 192 NAXIS2 = 192 BUNIT = 'JY/BEAM ' CTYPE1 = 'RA---TAN' CRPIX1 = -2.680658087122E+02 CDELT1 = -6.666666666667E-02 CRVAL1 = 0.000000000000E+00 CTYPE2 = 'DEC--TAN' CRPIX2 = -5.630437201085E-01 CDELT2 = 6.666666666667E-02 CRVAL2 = -9.000000000000E+01 LONPOLE = 1.800000000000E+02 / Native longitude of celestial pole LATPOLE = -9.000000000000E+01 / Native latitude of celestial pole EQUINOX = 2.000000000000E+03 / Equinox of equatorial coordinates BMAJ = 2.399999936422E-01 / Beam major axis in degrees BMIN = 2.399999936422E-01 / Beam minor axis in degrees BPA = 0.000000000000E+00 / Beam position angle in degrees RESTFRQ = 1.420405750000E+09 / Line rest frequency, Hz HISTORY Parkes Multibeam continuum map HISTORY Formed on Mon 2004/02/09 01:23:37 GMT by "pksgridzilla" which was HISTORY compiled on Feb 9 2004 12:08:02 (local time) within HISTORY AIPS++ version 19.405.00 dated . HISTORY Polarization mode: A and B aggregated HISTORY Gridding parameters: HISTORY Method: WGTMED HISTORY Clip fraction: 0.000 HISTORY Tsys weighting: applied HISTORY Beam weight order: 1 HISTORY Beam FWHM: 14.4 arcmin HISTORY Beam normalization: applied HISTORY Smoothing kernel type: TOP-HAT HISTORY Kernel FWHM: 12.0 arcmin HISTORY Cutoff radius: 6.0 arcmin HISTORY Beam RSS cutoff: 0.0 HISTORY Input data sets: HISTORY 97-10-09_0356_193558-66_206a.sdfits HISTORY 97-10-12_0142_182123-66_193a.sdfits HISTORY 97-10-12_0151_182707-66_194a.sdfits HISTORY 97-10-12_0200_183252-66_195a.sdfits HISTORY 97-11-07_0510_183836-66_196a.sdfits HISTORY 97-11-07_0519_184420-66_197a.sdfits HISTORY 97-11-07_0528_185004-66_198a.sdfits HISTORY 97-11-07_0537_185548-66_199a.sdfits HISTORY 97-11-07_0546_190132-66_200a.sdfits HISTORY 97-11-07_0556_190717-66_201a.sdfits HISTORY 97-11-07_0645_191301-66_202a.sdfits HISTORY 97-11-07_0654_191845-66_203a.sdfits HISTORY 97-11-07_0703_192429-66_204a.sdfits HISTORY 97-11-07_0712_193013-66_205a.sdfits HISTORY 97-11-07_0724_194142-66_207a.sdfits HISTORY 97-11-18_0256_193815-66_206c.sdfits HISTORY 97-11-18_0306_194359-66_207c.sdfits HISTORY 97-11-19_0447_182341-66_193c.sdfits HISTORY 97-11-19_0456_182925-66_194c.sdfits HISTORY 97-11-19_0507_190350-66_200c.sdfits HISTORY 97-11-19_0516_190934-66_201c.sdfits HISTORY 97-11-19_0525_191519-66_202c.sdfits HISTORY 97-11-19_0534_192103-66_203c.sdfits HISTORY 97-11-19_0544_192647-66_204c.sdfits HISTORY 97-11-19_0553_193231-66_205c.sdfits HISTORY 97-11-19_0602_183509-66_195c.sdfits HISTORY 97-11-19_0612_184053-66_196c.sdfits HISTORY 97-11-19_0622_184638-66_197c.sdfits HISTORY 97-11-19_0631_185222-66_198c.sdfits HISTORY 97-11-19_0640_185806-66_199c.sdfits HISTORY 98-03-24_2107_193706-66_206b.sdfits HISTORY 98-03-24_2116_194251-66_207b.sdfits HISTORY 98-03-25_2020_190826-66_201b.sdfits HISTORY 98-03-25_2029_191410-66_202b.sdfits HISTORY 98-03-25_2038_191954-66_203b.sdfits HISTORY 98-03-25_2047_192538-66_204b.sdfits HISTORY 98-03-25_2056_193122-66_205b.sdfits HISTORY 98-03-26_2048_190459-66_200d.sdfits HISTORY 98-03-27_2034_191627-66_202d.sdfits HISTORY 98-03-27_2043_192212-66_203d.sdfits HISTORY 98-03-27_2052_192756-66_204d.sdfits HISTORY 98-03-27_2102_193340-66_205d.sdfits HISTORY 98-03-27_2111_193924-66_206d.sdfits HISTORY 98-03-27_2120_194508-66_207d.sdfits HISTORY 98-03-27_2130_191043-66_201d.sdfits HISTORY 98-05-10_2123_182232-66_193b.sdfits HISTORY 98-05-10_2133_182816-66_194b.sdfits HISTORY 98-05-10_2142_183400-66_195b.sdfits HISTORY 98-05-10_2151_183945-66_196b.sdfits HISTORY 98-05-10_2200_184529-66_197b.sdfits HISTORY 98-05-10_2209_185113-66_198b.sdfits HISTORY 98-05-10_2219_185657-66_199b.sdfits HISTORY 98-05-10_2228_190241-66_200b.sdfits HISTORY 98-05-13_2132_182450-66_193d.sdfits HISTORY 98-05-13_2151_183034-66_194d.sdfits HISTORY 98-05-13_2200_183618-66_195d.sdfits HISTORY 98-05-13_2210_184202-66_196d.sdfits HISTORY 98-05-13_2219_184746-66_197d.sdfits HISTORY 98-05-13_2228_185331-66_198d.sdfits HISTORY 98-05-13_2237_185915-66_199d.sdfits HISTORY 98-05-25_1711_182559-66_193e.sdfits HISTORY 98-05-25_1720_183143-66_194e.sdfits HISTORY 98-05-25_1729_183727-66_195e.sdfits HISTORY 98-05-25_1738_184311-66_196e.sdfits HISTORY 98-05-25_1747_184855-66_197e.sdfits HISTORY 98-05-25_1756_185439-66_198e.sdfits HISTORY 98-05-25_1806_190024-66_199e.sdfits HISTORY 98-05-25_1815_190608-66_200e.sdfits HISTORY 98-05-25_1824_191152-66_201e.sdfits HISTORY 98-05-25_1833_191736-66_202e.sdfits HISTORY 98-05-25_1842_192320-66_203e.sdfits HISTORY 98-05-25_1851_192905-66_204e.sdfits HISTORY 98-05-25_1901_193449-66_205e.sdfits HISTORY 98-05-25_1910_194033-66_206e.sdfits HISTORY 98-05-25_1919_194617-66_207e.sdfits HISTORY Original FITS filename "1904-66_TAN.continuum.fits". HISTORY Noise level of continuum map: 59 mJy (RMS) END sPwľ'O$ `zģz@
+żF*p,E־F/[ k>V詤Q"UXҨHp^K+F5;*nҽD`:HtCʍmaO<߽X<P=tY=8;B=yH=ѯf=>k>68½@ѽj1wXLϽ -`n˔Yۋм߽g<wH=й=-<c ڽ5 a<<H={Ƚ`K _WojA׽C;r+p?;4sn<k<L=0L*=>μj<С< <#=
+̂;5===u>/N>H>?y?.z)?,,>_>>W[$T>?1
+CfTv¼=[=.={<(;< =Uf'=~=cMq3 O+l,c{7yt<2o=AJ== >8=1N=y< <$=]==T>
+5>5>GuC>W#?S?v?K>>,=Ӽw2rT| <\uFGk=!5=zO=Jͺ=Vq=HK{|k+4~7:=`d=Q=.P=<=><=t;el>&t=@<ޠ=M=%>y> >o0>==g3==&=Gi<>j>l=Ħ=17}I<cIͽ*. |V z\sЍK=0=X`>= =g; *J=f~=m>}6>>9>p>??S?=> m=\t='󿽑n䱽`A;^=S: =&i=/C<ܽo7m#c oL@9&ɽ؃OV<ټdh=kb=zoi0{y6"=OJ= 9==&&-<=<;r==<Һm <=؈=>Cv>X%>c s>j>> x==>!s>ɝ>r>Zko=-<O,<yќ>&on75*$%y=N=n>=l<Z[ϲASQ-=e>=>90`>>>û-f=qS&P!ƽ kT弨ۚ< &=t='=S<eZ$£-)mviأ܈=VS=;1u=KR8(h򑽃: z Wqgս7G<Q&K_c ==@"=T ;=&R==P<T=Xh=p=4=QF=Q=;/==l>-΅>H>>b->:<=>>;>0)T>2>,j=V5<Ԭ`<'G1e1Y罖7ԽRAfv l=%XE>(>P=m=K$Ž{Qo佩ýo8+ =W_ <:AϽq{ F%ɽ?=11]=!=>Λ>,=<H,G= AO!nԾOC<KJs-7{<Žc˽CuM}p⽚L< ¼]Añq,;Iyp_ k Cjkaw= =;>x= == 9=p6= J
+6l<s==ae=)=a== |=מ<=>>>> =&H=v
+=> w>E=<l&;L; =PbJl}V˽늽9=7><䗣XHą`=7=='A=>k==B酐𫽟޾"xoỼ3: ʰTxhҫ/ֽ?Jս+:y=3=sV>p>T
+D>,v> A=R9w;Jl<=]='i<?;7ˌuzD>%zv
+ȽV9pɽQAU~𽃭4rUY>R~e}:9gzO϶D;"bwμ.= W=ކ==z=+=.<D̽Fi|ݽ =
+jHO<y\ ̽GG弧e(R8J2m=^q>Z>Je> m>n==Lt==X='k=
+K
+ٻ/==xϽUF7ae]>:====D?( F=^=^<>
+M^>I=<ϼ&R\Fw/C ܾZƾ(IHF$WF`[6crG@kyAἢGvA7==h=ʙ= =y=lɻɻJ;r'۾JQ =}<=oɼASitzn0[pؽ>߼HS<n<=[Ѽ*׼Z@=>g>UN%>u>zA>a>Z?>b=> '2>&L>2^==M$>3+k%Z]<=7===놡=Y>!
+t4ۼh=="C;};⋯9wG߫,?z-cd=c<=!뼚\=>
+[=){=>l=o:<_<-<<Za< ;.el:LhӖ^Y@$m<<Ӡ=D=9>Ic>vx>vc>3=V2=c_<,<&h{F<[<wP;悽P]p4P"ƽ\Tl<B-!<&7%#K!w>|Qn<u(Y5^>U&p:
+/A=Y=3r=4 <<Nl1<=:*=r==sc=YG5;Z\1W .QͲjGffHS],CNGfn'þDA!iE9
+~.m˜<=.= H=y=8;0<=v=9=@=5Y= =ϛz>x֯>>>>b->T>Q>N>A6===
+N=Gi=Y=2=hD=2<C=HJ=\u=<w<=dR3<R=U=#=y=\=CU> :==&3Rp>i?#z>h>U>;g>"ql=П=0=
+=Ze9=>M?>Yn>y.>s2>Z>"]v=؊(<}弇:_ۮH$u _;Fcͼe\`F R{p<ѭ'=y=]> |>e>I>=ƌ; ׀Dz /AyƋ
+ڽn]rS9Ge[ѽ\~Erټ61弒*T~:H#Cv ƾT"ݼ7Ǣ-
+$=L7=<66<W<A ==)1===!<6;kL:%YaF-
+тx㽽e
+Gݽ#b'ΑP꽃u7뢽N{ƽ9ý#!,rǽ(,.ν'@f'"Q|
+(ҽjKl)XKz8vD⽋%3u^pg;WXusykKU;=O=B> 9=M==Fy<g<`=1=)v=&,0:Zo 콡[9y}BHuo+6<FdOO澽#` ޽>*9z7k@˽Umýq%ݽ/:(# T; rz;N<rS<=(===~\=79qE3n҆c-:Ǘͽ䧽 쎽G$pYtڽjT:(;~<= ^=w=$n<iS8O
+"׽)\#õ 9Q=!>&ڸ>,>t>>?+? 0>> S>k>>>9`>)> %E>0
+w ė½ޓݽ@~[G0wd:FgmCU8ǽtW0{:a<x=5= =4*p=ݷ=o%==o=V=7>Jsżȼħ}5`ܽ|$.^4ǽJxk}d>(*)㟽QC7f7*)"e1 <+;:`<z<=$Z=h{<#<ѕxQ$pW>H_EϽٽ;[1y4ܽ:-kvlN*z <6<V]$<s=/ލ=`X<;/V;Z!?@JoM'֯jڻ3<<=3`>>@,>0>">To@>c\>'=1=^=
+>%I>>cC>c>;=͐
+=n;=ཟ=4=#== K;EK=#'{<H2{^27H"ֽg [ GG˪}N<CRɾ}|a=uX>3>Ll\=<ϻ>,t>6V>2[>1v>1=@f=aA<(=;><>v}>@0>A
+놽Ŗ6˅ƎhɃн
+M=$>;>_G>r>? 91?Vl>_>[>y==
+żԷ&<VOWDx9)3ҽ%5M; <] <üaV"4>)]> 4>>aO>.=/Pq=>N>-t><={>==C<j=n?==}$=qI=)= ټ(@(׽iŽ/Hʽ#.ְK9T౥3eklu9oMiWf@ɽж޻2XP7.}:-~Tǻ<GQ<0=j<tGux5u.OT;;]*4`S,O/ڽ+ ]{J(D]tq#2L'35,g mcTw3fL<M2= =&%=*< <-h/ge|5:Ƚ/㽜NkλC:<3:ǻu&VLYkܯxbh8ChjٽD̞ <;`ٺ<!#91N<j:$Oe<'-=җ=n>3nj><>)?Y@5@x@@U)@/?a>ٰ>\=!;>D0>G=y<m<d9=/<SoG<;&.;<Nw4<yV=N<w= W;*k=YȻ_Ύ<I<= f=,<u_l=}m=G|=>=+_=b=jE=1=>b>2`}>2b2> O=/>=d=Q=m=\ =^==;F5<=P& <=E~=u==);OգmYƽ܇~!?aýS4l
+?p6ս֑ } jdEYyӱ:xͼ'S8}=0*Lj-<I<oaU;!{xSd2 ʜL;k; `<0O;H(:ߵ<!<Ѧw=*~=d=pD<>ƻXB! /#Xs hVϽ gX̼{ȼ
+ͻ =4:= %=E=vO<}(?Hd0˽x0^5<@>z;B<=#*==e;dϽ扽!cڼ n<Tƽ_{h<9mO;қ<<TM<L <g<1=,=Ż>=$>k>?~@?r@^,@Ͱ@_@' @
+@6W?/>>& >Hi>J>19=Q=|AW=&=z/<;yۼUmi2gټdx;Yo_<ym
+H=(O=\r=ki>9T=c=;m>"|~<=<*<R;޻;e<#һUؼWO廓#<$<<x;:B_R\yBQ?,/kuH?ǘʽg~> X!VxJ.y'鵽p@򻸥޻mּ⽊V7JNnOһM5/xd<485`%;RdB<1=n
+>^R>?T@Y(@!A
+t5ht"V:DR:)8 &7S`> <>O=R=$R=P]=9;z==S=K=JE<F;: P!
+GLlAV,dPB/Ni¼#5(:$KAF
+pe"<#=2Y=(Q=5^=<==@=>{>m?x?T?1?\*>x=jM4C˺ٽVYϽ}ltos,꽞2Dۦ<Ƽ*q\ dm\ofyK:;<g=o=[+=<&=T=->n9>>Ņ>f=bw<ٽLGEK轀NF,Iz>g8;rX=4=Ķ>"">j>(?yVz@'@@͞AfA-q@~W@ز?`>f>b>,=ya> s>*> t>/ܳ>*=]:nĽr૽ھ u߾2BX=>p/'eJP"پc,^>L CJؽ[`Rtm ^<n=#=Z2=?!=*=jD=SEsi+ 9.%Bݽ5#L5#Fw.:G=^e>I<NL<(N<@<.<f}*ul7$QG
+;(meWD;ܼt_޽mhZRS??oF 3+_eW0ڲ༼<>=9=́=8==l=ņ>>?p???p?q?"K`>W9]˼A+ƽˍ̝XzIT䢽Aʽ+@Ҡp}>6qNAPK1l#rażϠ<<)<O<'=o>UO>֙>I>>>I?>;fB E+q)"3;,+̽L'b?7t;I<[%=K>m>?.޺?d@<@ckM@@Q@X@<@E?|>=P="ʼ;g-=[Y>5>C>.>N>*< 彅UѧF'.<[KF+9p6 uƆr_Aa<|+y6v^׽dn3~{zRg=P:=i=)=3=ȹ= p(;Y޽ eRyyry竽@QX <a=EN=OV=2===A|꼮X[l!l1<=]=̬=WZM<[20ĽK*潔65ڷ|ޠ4>;Z369[1DCgw96lf~6v<B=z==ev=3=[==>>Uh?n$?3@N@?X?a;>."wLSG G W
+<ɒ6!PqٽĽL
+>  >+?J?X?Q?)}>?Q>?p=1<.+1Ent9󮽸 ý\7!h[q29a=^&>ׄ?d ??R@'@/A@O}@@@Q?>T=Zl=18!9ý(\1<<=sT>4N>:&> =`O &.ăwdlǠF9I3 7NqоR1v9vm +W"h?-*+Mǽ(WJU4>o>0>o<dഽL?5~;9zFͼ<=L=LR=6=O1==T=]I='d}R;&u=>)>ӌ=Ja< 7̽7>Fo;P<{G<(ҽJ;GVF
+SϼA?wyir;i\Y%Ii<=F7=;A< <O=a=>s?[d???̟??dY>qû]佾 #}⹾)k DʽҨ3>:. Ui׽9vmI ὧzu-=h>X%>?&?ch?Zy?U>>#s=Ѿ=<*11ؽ⽡}彅k@G
+Dc;}>~fO?4?ǒ(?E@@,@ Z?Z?ːr???Z>}'K
+i`.w\ۧ A4 8y
+Rl]={6>=?)?V?+?B>}>6.><rTlWνH@Yײ'VDC`ڽD뽔 |e^R:OψI̽nھоt"Ӿ S W`Zf|<T =l>R
+Fj S0͟ =ݫ}oֽ6%-,s}e sQ<4J#=>]>!? t?5?
+۹>>>>N=h=GSսhǽg5,SjL<kB>#?.V?)~@܁@@ ??U?X> >'=8VTԼ|,ܗ½.o'Mft׼ ;d<ח</&+4
+d<i<Ļߒ0rļx0rн~׷4&ܾ}F]ʼ8C'6-Z^p瀼k[>8T{==(V=k=3=->u=s>&=^= .<=9Ι= ɺN_/N¼s%B#<G4<t*< <3=' =#%<6μ{ш)ޛ:^<[===S=
+1ѽϼ5<O'< =
+x;6żFпs?֟%߽?|1F\{0I7^>|˾Ӿ
+,)aUE;=l%>>>^F>>Z>6L><ۼHм{nqGG
+JC&}3]4ʽ>½S?cf1 )k}f޽jNԼf{!}<{=?>l'>@>&>Z>a->|>>>Lp>9| ߺr2~ؽHf"$Ƃ`< >mz? l)?GG?e??4%>=`=.x+q)սmG]{󇽨gH+\i@/ :z;=,<7<No8;?\=[Tgcr dyx_,*& 2iAYsۼ6Y>C8b>>">-O==>=> $>tE>XR>SH>#>s==^,E=qK=n=;C:@;ť<Zhb="Lz=X=Pk=N={=b;ҽ \e-%w<]V<¼n<>x;7 S3x۠B!pNɮs1\2guTսlRJRk؈}t8F@&ƅ52=@+>>e>>?>[9>%*=f;gm^ԽF|;b;LD( " t*:
+B~d gZp|h+:;;.r;I;;+ͼjF?8<=N$>&{>Y>}>2L>>>U>r>Yz>5H=x<m=7
+Xٽ
+41o'}sͽkS]QvT=>,>$? >3>~> =W^)?3k*<ˏ~ھ&;߽Ű۽4*Ľ;/%XؼC<jT;B<H8Hfh5hKp/y
+><GP#a' <FK<<fO ⼪nFF 8vu_3.l⼘k"<1==`=oT=@s8=<<m=& =>'>,>+>#<>>3R>,ݸ=R4=ʭ<Q Fzf 좽Cy[~_nY潿(F no'=s_>#>`= B&zlͽHȯJ;
+`==S4=j<u ;t/gls:FUخƼoOZ}JKEu;MT-kGb󽰔cWz>nk%n[%Խp
+i <R=9N==@:rļdp̽R^1I
+.g½Eb=J$:g:?},?[?<?EA>4>E>>}F>7a>F>"~==Q=<}5<<w=B<A=;;ف=6E=8=kGH==*+U<I<k4J=<:\`
+<;>;]7>f~hhW/ZI R 伿gQ$Ӽ0!˟
+ȼB(dt턼26.lV/Q?YֽE<|ֽ!RB]H+<D =1>=Ԙ=:;KTȼ~ݽ C)^Vn*{-g=7t
+=d>(>0>=`=j=-"=mo=9=Jwe==Ax=>1>8>"=q>t=hq=#z<\ =MZ<u;1AwU'+˼{5dk/&;_T;%Z<y;;SνD^~b{cR񈽄нc_GVk6{O5ȁ <=4 <m<LZug$2żQl~=g= > O=;*~EǶ𼍹뽥!LϾ1 @u|P࿻)l!<=p?y?@W> 9>"f>ط>p%>%=C=+==o=m2<<U<4<<k=4>=x
+=XS_=`=A3=?=
+6=¼>#♼=7\;:t];<%E<{cw" _\ɹϭ_'c;[EOP.8>a.̷MCHEˉ~X
+f
+λ2<]<4<|;s_a<7Z½zνm[DApR{~⼢)]νg@r+7z<
+$Z=@
+<.;; ;= <J <B6 `FIG3!4dAQϼ' y1
+O㲽WA""216n٭k5=M<->}> >};)>y>U.>B"><?=7={:=4s<E4'⚼~ey'&;<=8< 86=Bu>uJZSUܼ/W!|<R%=@P=:<D;{7Pp=gZ/9.JhHdzؽ\ǴƻI~4ػ,;?byzn>6)ǽG)gM/ ּ|DihӼq;=7<<;jvf;]4佑11sS>м:=W>>uX>f>*J>#>>V
+>2֥>!w>8>~>/> >>?>`=Ѹ=;;][<z:̼TV}2U( K@9bcY+ս.vTK4ᢻx;u;I/y"ƽ Ƚ>($$ʽ8)ؽ.`9/; zbzrCyO:@{T.5]1D;2^ܼQ:sm;K;%:|8;<=NNo=ANY<tťO⽴'<{޼ĽZӷOXo=t<A=>@w>:XF=,>MrP>CF->x>2d=tn=T=C@="HNA 85/M{QoDj*onV*.<K`<١=)<t7RF 1Xi4QO_-8mO<=R(D=J==V;ּB
+ܼGph_fX!pPa埐]땼oɻyuT¼ޭHDԩx)t}EMۼqFIQ`_ռL`;<uemӽBitW+GloV=5=q[>4\>w>j>>1>{X>D>1>YR>U>
+%t;tEk<ɲ<#=U<$<=EO q+q)uQ<[
+<kl(<]<"=`ޫ==[>Ŝ>ف==fSݸ
+==]̼1
+[TA*3y7*-N<H9<=5>ej> y=<z¼ĵpin=ƛ2{սo6Bń<=T=($=&8<9Ƽ>x:'XhYCy ü"-2n /Q̪(KԛRpS=(oiݿiDr&; 9ۼǃUZʽU9 3°q;ͼyM;Ƹ<= =Q>!n>X>DG>>u>MS>"h>R8>)
+>Q>??>~n>]=t0<}Fz49Dշ#w;@b㽼FyQ彤Ľn_"Yսms 㨽y^;!D<#<];WJ4m7輧μz11;s;]-t$6R.Pڽ7=v ( ս<N5.IѽOp)G/3;J֚<-
+>tG>H>L>>U>;=z{;J/ D2Ahü3½IUݽ}Ž@?/
+vm<La=7=nF*=Evt<=*Q<ۡu_#xq"s)[8T @֗])q}Ijuj\n;7ԻiK<P<\<Y&ʻ%y <AF<֠<oFJR1ܩCʽ:m*]}c揽ݽVm%ֽ7bY= >Ć???
+>Y=C~νh&D߼6
+`
+==Q<΍;a81n
+MAXG0'-;d==0b>^*>bdy>c=I=q`κ ڼ<#JG=17=ͫ>8|>v>Q>2o=y====.==yf<qW8۽ !\Y0.]sp?ݼt_SNc:WY.hKg2ud3нý@D !;D[<l<l<P<=ڑ;Ԗp;el<D3<^d<ghM;<8=e=8x===?6;|ϼpW-I\}kE7҂1 6R ZpyQ;ze=F-мGgr=ds=;σ;K:2z=
+n^̐FCYC;=0y=s(=4f=JQN<>;K:ÓI⼞'Ͻ.ݽ0,naoW=>z>>h?B>B>A+=m.G;Rt-
+R;m=>>B=<w5><3;я:=.=y³==z<w<Op<oe1uLS
+R8r!nooPI⽑%;ǽB^ BDH಼ 45<>h<<D<;ld#<F <;7@<<r=$J<
+(<:#=B'=t`=Q>=C=3<Y H;lнfS=Buk<g*Rem]mFק}6Hu+y(5
+߻M뷿;<|!=#o?=M=&<⹙+<0<B=)ө<io<;sɻ~`q E<=y=< <:wglAdDJt뽏 ½6>(~?;?v??-P?>>T=u`T·ʽlr#2j:t 2`;=y=K=RV=a=p(<GȽH&_ܳ
+{j::;T<?=D=+=9S=eD=Sb=c=?=s=|<^;ٽ?c;ZQ;Ov<<Fɼu᝻U:B~mJٰ{ F$V[ؽ$Ԍ8<">N$?Q%?G?ֶ?ћ?v>?= 6 ukнٽٽϽp:W+}`=<>=> >%>>=[R'GY|)< =5n=qC=jԸT;;go9۹ʽ߼蔧tr(EY4r{5]4~f<%<i~u<k;_;H;n<pOv=|=;T=l)=1AF=<Cf=O=PH=3;3@k@t' 1,\`gqqj8ܭ ޽틽S߈OO87µ&$=<^<='-=&=.=Z=2f=P5'Uoj=<8޼~,0-)$ 8ܾ | Hž/9&},"6ts(w7սOm
+^ M=Ͻ|ӽPySmB="¼*; G; ;<g<g=8O='=U=5====Eh=<G</< <n<J;}zU #ƼK9A8)V<9eS=eb<+ڼQ8ȼglQ
+#8]T3J)KY6BS2ozsb8,)C0sjڼEź[)j<{=6e=b=?֣ }<TJsٽc$CqIWgG
+;=|I=o=B==Ҫ=GI=!=
+<Uo<-ltR<ROB,b##'"k;CF<<UW<)=K=Ah=
+='=K=,V:Ps˻̻`h1'JR.
+Rν"$lX>I\^q嶼5_$t<<̜d<L/=2; =U==,=^,<Qe*Iٽ]9iAPC<<B=<j黟hv>+<֡i KsC
+c>"$>ua>>
+J>U`>\P=?=H<|8<=="r<(x{l"ׯ#;y;w!<Im= <c+Vؼ!i bB;uuJ6G/{ڃӽ6-4YDBneo}Ԗ>u=Ӻ=/>->a/>_9>"=$v{<h=Ė>(>$>
+BidD5;(<0*<Q;ͼ_ƼѽBd$#PK˼۔f~Hb<FѺC۹@m;d-t`A
+(͟C:QܽSͽ%I{<=;\hg)'=Dp |~[7`<:3<<s<ypƼ箼SՏ׽#>#T>r=ʯ=<8:M|яSu ӽF[1ν牽FȔr"ҽKHh @ɬP-=z%=>>C>*> έ>Q=u=W=8a=%=hA=y=sv==^=5(==w=H%=9===)> =>c>M5>]1>|>4>kV>1g=j~ <BJ<.G<J<<\ǺLϼ+伄üi:.si;A==d<eK#L̼{و?PӽWve܆=-EMBCѻbGoi+сcxǽI$<u>}>>>>y>0=H;4< MH< =B_>
+\._dY4155<= I=U8=b==۩<'g
+(<7,#Lhz,0__=5ڽ~hNxM5%*{;!$-Bg|k2:_=d=>s> ?>rD>va=7=4==n=|y>,>vJb>j>G=+<݄߼d2
+h>@%=(=s=A}C>a</ةLs<`$QJ,yhFs'C
+ݽWp1/Y_ӽM*R<'h<=n
+9R98H^J2$Ӽ>Z.5Oͽsܽq\˽,|rp!@:PFh:P;<}[b
+R>\>Vׂ>-
+=E<~(gLZlnO2]&sv5s{E8Q齆;H ͼ7@P芼k^SIGT\d9UJmD!Q8wݼƽ0cw̲(";:; 5:ț;0<=8\= 1<Nh<=4y=W=>,g>h&8>zZ>Sl> 4!=1V=<7<6y˼d<rNǜv۽8b=7L=jbѪq<<*= ==v(<YVD/溿Zh<=4=G=Q5=*<=
+Z@=O=r>P>d>{S?> >>j= ===Β== =E=
+%<ɎU=^=2h=e.=X=E=<;tQQA:;<E<Ӓ_<;ǻ]}N!9j }V(SY2xki~bFr k9-KF0}@|9<VH ;\T<N<<<%1^)I驱mй =m>; > >>>Ey==@R<<r(=!F= n>e>S==<raݶ@tv!;"MfKRSu̽aCHi꽚Z
+_lv/9½9E%[S7kSM ;׻kkSCʗ;U䭝*mBP_;i<<X<a<hKQ< <=b=Sb>>b>u>b)e>>!Kp=>=SՇNցƼ3ýuvt̽4mqCӽ@ڼ<<*<[<<a<+<I̻7jzڼ&M<
+=)=W}=3*=Q6:ļ< =q==3=a:T<g<E3}VθYFv|OY>0EݽAc<
+R;#: -6.lg3:zJ[XN;L;3Y|= eC=+W=>>ET3>?8>>9==+g얽6 ?0ɓ5Rὔ罖D%!½HO{A9h=
+=b=HI<ͪ/<nD9<.+$ dE @N^pt<b
+'<=v=Wn=g\>\??N?Wg?:
+ D¼>7=Ub=>
+L>
+>#=|<<\*P;7VF橼:39l8!27VD~L] @~PNfY",ڽ|'#kVؽϽ?[P;M@@1Vݽ^~O?ڽI/f'Cn{;ڽ
+9
+|q(Co
+T=\#TS^6׽HxM^;=>2%>=^>=5>9R> |t=:*==~=~E<Y?#$?Lb95M0żHe_w:qd;{$V<}b<H#.{FvIUw}Ժ=i
+۽ے;_< 9==[=y4=p>,->h{>uk>U> Z=fv:;e)U:\KQѦM(|wUk8l,e<<<9T; ;VZHqL֪;UctȔJp
+.UYҽ==>==<=\w;琽5w n c>&sD^i<oTPo2];~P<YU<<ե<橫=ס=`2=!.<g=X(=u=D<*˼*:ޥ޹]h̼#>ּܔ*racEFzTJȞڽZ彷ƥ
+<!C<Ȗ'</z;6Xq35ټ:N !!`^
+, =Mؤi<P
+mLW`qUԷ> G= =:=>Ў=L <j= = ===m:=R=|Q=/<<}@E-O8B1<A,;Jq<*E<67<o<|=l==6=!G5<6)< ;ݼ`( r~һ蓻c_$u8o<~L=>4Z>gگ>n!>#= r~K#a!!M6)G,K[kT>5 S<:<z<㯦<翤<z<></{
+=^="=> a=.==ݩ=S<jF<߻M/SǿQ齂zBQ b wFF%0zNt0׽k~u=7n=w=5=
+V==
+O%s;B$=5{=
+S,yܔcH;R<=gK=&=q@F鼊¼r=u}<4̡<!<el<:!lӒjP<(<U==F=T~=X==v=K>
+>#}">,G>==<!٪;fwJ Z9o:E憼 I.gLѽ(@e{޼w)& ![Wʽ"ˆ>g۽Ru=A~<f=WΈ>̏>SD>E<>94>=s>> >v>p===&[<^<~<N=,;^{Cx[ԽP34+`;<Mh<f{u<<M<3E WƼL'.4ݠ&i"3*a2]==ԩ=x=j= ==9==黠==P==]==v=<E:;WVmX <L=
+=
+C==T==ż==B=>y<Z59̻;<T1=i=q:=r=C=%;$<(1<E= G`==='=~;WGͽj2NAj`zw
+;T=9s=R=NH=S=2;P4p2FO %Q"巼*yB[*U;%h;9Vvgk꼳Pj ~4Gs <%=L=>g5<CAy;Hn<+ۥ<1=s=) =0[=a==m=J=fn=c=+=a= mc<"|;p}&Rmω|`׼)<U/xX}i$6ȼ6H
+D3ŽLEO3s;p5\۽׽s;a=\P> &>O>[.>Oi>>`>v=ܳ=(^=a=Z:8^Zؼ \6'-`vY7ⲽzVӼ$q;Y=$"=UM="==bʒ<y<<O#(;ʍ׼<F ;;h㢒p=$==T=1>t=ذ==J=<J>ϧ>V>= D=<&<Z:<(L===Q==g(7=W=G=t{=U=I=_<`+<I'='=_>ߓ>E|>1~>=K6<3&>=< =C!j=]=8? m4ڵqB
+JYL(ֽ&⽇AԼ%<=
+Q<mлi}y }
+6h 1bÕ3;Vl2=P=9=i8=6>7=§=_=3L==G=Z=J]= E=kx fٶU6n#<==_=
+<)Լ 9Ƽ-x ZZҼ&?缣=NqKu%fýgܽ=E$Fͼ.80Tz;
+;4 <Fz y;%[<fU_<,<z5<AM;O<Z<c͠἖pPh2v߹աjp#XD :pK}#ͽqf|Jf+bkj>y$U=Aku½ iT ;̻0z<=h=R=c=ƶ=J==h=(><`^j24B ?M\5+J\ʼBkݽ0%a3%伊ٺ
+s=r>>\>b=.&=O|====Z>!jy>*`齛&<S<G= <=dy=====>qo><>
+"=K%:\Xc50Z<u2Ĩ(
+׽k銽k"Dq e3FlȽ$OMg޼iL&=<====8=;=@Ϝ=; ;;)<T<5<Ssl;"㻢
+0hμ@xѩfϼ'v0)9޽Cp Q;;i6`ټEJƼKΙ<Z^<{<;E;4]"~] ؝+j?Ua7MuOhh`.ɽwRoo4ܽUA?P"=QC̽C鞼g8;RҼü5#Iu䕼绉<5 <B=3$8=
+\;gvW&L?Ɔ͙< XA.r%arӼt+FW5-XicX<<\=J1=#0>@> *=O&=l^=h=Jy=J}=7= [>P= <e1f<<<=8 =k*/=@=b=Z=ni=
+%=5s=Q\=e
+=P==e;寧V);g<k<D<^<<7~<_4<<V<u<+_C;N;IX;X < >(>Ma2>~=>>P> u=;YIC26ƽ(=
+O\<ʟ̄ƻsmQt*n<
+9콳b\ǔd$1;ɽPI4½罶t큽rXϋڽný<j=A>R>,>>i>>n֥> <=Ń=xٹ==Ys=w=e="=XK<;Tȼzr69&_軁!E `MO½% ռj.:i8;*;ϼ;!;< =Cn=^W=OK<gA<Կ;ڃFgz{Y5{)뙽eؽ`9Ǽ?>_]
+<b<G;Ǽ@Вpdm48ʼ + %T|"E; jS&I4q\u_ulSI I8=R>>
+<X<b<8*N<_=`==<2=
+97<Ky%w(B/N ]o<<<ງ @
+s96
+7>c>/ =^;==/t`[8krh`-'мKD*o9mx)`<<d<Q<h<j<j:< <=JM==zI>M*=F=(<< ڛ<!b<:A Gdo5;QIetR켆*<:ǹ׺*2:<<H=3=>>>;=26<Kj༡*< <<<¼|tHN}KH:ùo: o;N;A;W<
+$.5\GuC=dU弰i߼#'üy`aAټWII%dnex
+_HAk=qLCh=qYV/QۘLQl{i2sG1̔uʛY!罠ѻY=>@>#>>>0>=ﲕ<)>1̑¼b1F)>ha⽈;sc&\ϼ_:C;1<<@<<s;Y.9:jcvī:Źv<d<e=#=A=w7==,=;;̧7;;Ý;Q%8k,&:_^9?g;զ;<%j:ܛ<:y;;լiz<=y>[>9w=e<(PfT|޼Q) &^A42R4 ]Վr?Lļ1&[KiJ'`(rFY3/FJ׼_zܽ\jT<lHS~_3݈)4-,V@2v׽Q;C[ݼ #ŻEĭ<<<\c^A㖼ү;߽XF1'_¼μ,ۢ;$3:<XG<</9.]Zm<ҽqi?l8*%ƽ Q,@"Q;|`\:Ė aď֩ q]#sͽqнI۽P@Hi mn5) A>][5ҽhսh>
+>;>AEf>5=Bm<>ERcC٢--:~޾V2(팽N
+g;N%.q*Q8~J'3(ڼW:_I5s :Q=;^fO1 yyýr-rDqkL|K<.Ǡ<Z<<;,r_hm`)",~Zx׽:
+L57'Eؽš||fU~>⽡9Ͻ񽢥Q2l'н3Z F
+=zzr=<Aϼ^;-;*jŽ@! %C<@Ľd=cS#>fz9>*><ԕ>]=6=Hr=(=<=L<,<o =';UμP<SH<w=HN={<փm }~Җ^J׹{zgd}>I "v G#7࠻O;;H,^H8 :lٻ<q<=ʽ<yw:#\j:쐾M/v!J; ==>&>h>O>}@= =b=h==x.XԼt+q)~޽XkO;0heIa+DԽXps1
+[幷*_?c\Dt==>1HF>>>g>R}>W<KDļʻ_<=$>:p>TQ>l>j>e>=U=0eꩡCK*»oǻT6Aqj0xeh sC
+>Zu>>Z=G==_={ =B<<<lzk ټ⛽1u1u<N z=%"齴X^%o i 0 iG齂UaQ
+=3<=k= ==<=,v=N<<A3o! н odǼF@]D <pe?="_=(>'٪>Iō>Z>;ۿ>g=R=SS;]*4.A&"wmhO (@r UPq=MRj~+^&9:92;1<*=w> g>=>g >{UH>1G=n<;Z$Z7ҽMBf;i-:@#ߴ
+<M;=ò>>2?F?Z?.>b>X=C1<<
+LikW,T;(;<{j'<$߻Q༟9μm6ؼżxȼgҼdkZ(+}v{"<<Ȳ=n
+=ʞX>J>`>s~>>v6><01==t=V0=1<Ƶ:*[j:!l^~Ԣ;Z;Rb<>'d>o?<}?h[?"??:>=W;_|3 <N=t>2>OZ>a>,>ik>3j> =:==5<><=7$2==T^U<N-#ؽZHsڽ`YXd>d=ʢ)=~==*=#[<<!;8 t@髛g"'fN%_i:-,x@`GNZ@zؽ\ ?;+=\>N>Wa>r4>s>: >>N>* >r=LU=k@=Bˏ<s<Cn9b`kJp驼,Ѽ%E<5S=S=y>44>==$'= ;#O5]L:o<<C<|8LIh[Q 7<x,;':q[<</=>`?$!?P?`?}?2>e=Zn=<P<;<
+,׺+<N]<=!<J<żzwz{D{?pVxX0dhB$I"L[;'<V=4=D=;>Rn>=և>x>Ю>B,>=GZ=.1=og3=9R#<<<@:e%<,W<l<<?PF< vL<k>}>y?!?u??ږ?:,>(>D ;6lT;&^ļo<E.=>/2>8>=V>====M@\;˭*@<0,=)=w={=E}<X^'5yl+WN=ݿ=$=*!<<3-v Y>A>]ҽQ䨻 ϻ<:@=/==Ͳn=°;𓽉2fP2h #^B=-%>9C>{>>~>>>K>2>A
+>N4]>YӜ>7==AJ+<<<fx:79i/uuX:PCۼM<y <ğ==T;F== =g< ]r]OyR;<x<‹<;b,Ի55 ; <<=U;nu/<LpW<]K=i>>U?o3\?b?X?]?c>[>7&=Q<Ƕ<>޺ *ѻ;>S;P<
+c==>\MO>.?)?;K?1
+>?>j=:=阽H,xὧ>Cyif<Is==Z+=Ee=_;=O}=i<<WR<2;KϼZ<r]=Y=z^>mE==6N<w'Br2 >p =!|<a;2;7DcW4Dpnq>?di@PCr՟3<8=f>%$>,>=|CJ`󂽶Y}x6# =4>XP>0H>>T>X>f>Z>V&1>]Ŝ>P>![==7<|^
+> >!``<\<-)<<`jd;Ph<D<U
+;~J;(1(jZݐ}r k+Źbӽڗy<%]<$<<;< <[=R=A7= ==|\=ڼ=R}=o==~Z=6*<<@<%
+<w=2x=%F<]E< {=.&>|>[>8>G> f=_ ;[rN۽~$v}&4yI< <Lv<^|J;vɼ1o:҂;x;==R!=>3yk>BO>'$=|x<&;(J)~ۅǽ]=xkP=_<Ϲ<Mb<9C轑Ľg59}˽Z͕3;=" K=>S">B>2U=
+! glR6>;9=Ts> D>{>\>z>>9>zG>G >7AE>5=0=a=!z;s_ ȦT7eWCD|5+4Bc 4QJ-HQ)災}q-n\<Y`<`5<a<U";;^:ٹ<G==:L=Ic=g=M<!<g<g|<;Z<W/<C?w<<e:=D>H>_??!IY? >>=y<E<;4<5\0❼' :(~߼ļ0W} 轖Iet v#59&I[<G=21=5<<,';<Nd֭<)<@$=%=W-=J= W=ι=p==x== ,<= =oCV=w=L\=@]ؽ)KH;= =%1==d=G<Cɺ_ *ogY5.vs7bsIm  >[>#}6<<S<= =@C=9>ax>j>op>I
+ʾv8M|=Wh=>R>;>) )<!;q Neв2̻0sZ]n4Ǎ3,+i<߫=x> \B>=פ=Z). ؼ=ʼ?XN=S=>3=>oFO>5>:>z7>N><ј=R=5<;%<ũ<'9;:`hM^8mt9[% OY$k9`o$<=M=]=^L=N=v 6=~M=f=A?=5
++=X[=\Cf<%:ʼn1
+^XĽ1
+kAteI;<f<=UN=Nl=N==p=0ȴb. 7j»gFYx3̐E=Vug&cTT8~UN~gu=Ul4)=/;=~=\K=o;DOⰼ,.W_ZW/:<n<3</@=MJ=8|=%=u=;=i=3=
+3<^
+;< <`;J! <<#6;5{䘄:;52B;XcT>!;AdҽgDF
+F))~ʃ輟T8:iS幝Ch6<x=O_'=,ϼR~G$lݮ 7[Yu}j;';ұ<=)==Q=;=_K='<<n[R<<A~I伪
+= k =KQ#=$=V=]-<hqWߐU7oCq~_I&mr;J< =! ={p=ɀr=R=݆=v=15<}y?* Dq<= ;=>=!{<N]BG?]t]7hٽt̽LTo&rJ>ӽ8ξ\<8=.|2=w=[==L.=l<;;YKѼ}wC+=(+ԅ;λ# 9;;<<$;qr:|;n<9ݷ<zXM<3f
+]ż_j;gT9޺@~ < 6%;T:$ڍãN2hNZ:?fM-/4R'gϽ9é)ΛŘ޼ #C158 4+<!=#=K=xУ=~C=L<<V8_<xrےwka}۽PyZAd
+׽|NV&G:s]џ:)govp5*\μ睬`nvq:w$J4<+)<A<<X< <~84笽4x v4==â=~YF=Hٕ= ?]<==^=8뻖r y̳W 4Q(R*?𽠫l
+~޼wi;+I<k=A=?oZ={=J=<9"dȻT+鼑:˽P˛DH8v<i%<+;K
+q`.| Tcg;h<׻ﻉ?R:ArMF;eV;g<<<+<N<I;T9B@h%ݶ bPP'ǽ.н Qbbܼ <fO<;t5;U:[ViD6q
+pkսa1_p̽B~y=AKw ԼBۼ樽Լm4B74<]o<Xo=F=t#=u=Z<<1<<A<v i$'X iz,,/T;v-g%?r)p kp"7kۑFBH/"ۘ½AʡН(:;<<%< <<]<<Uy퉽ʽĻ`=P=@=y=ĥ;d=I|===9unǐ(<< *8]J Yoz: 2=%]=4=/<BT<u!<&< ;l;O6G&<GɽR}P ) B%F6 JsEJճl|;G=
+*ʽrT)Jh-s{Ӈ0弼vWPlr߽ c0K/]
+g
+8 5饼}MHi]զ)ϼ𴐽GpD;R<=E =V3=kC<<o.;mt ( _Խ=7v;<=\s" Kg3p=KqҢC#*9<O<ϼ
+̄qz+fZvb;< @Zzɜ׸$G9; _:SXjoF>?6񱭼`r
+7qP3!Q
+%=Me=:=8^>X=i=e<CD
+ڼϸ؋ l? ༳8׼jSƼTfM˼ 53F`3M%뼎UyzV9\<E=5==o<:1ꩽǽC & =a20=!=1<L0^wc(/rR i 8{aI:~K;-޽6@~iKKӟxw;*%o)%
+MkJDJ k"Qc]`뻭3;P_(O颽<: w su]n;/<K<+S<<=j=@>f(>0> rc=g=
+繽༸ռb;;6b<ш%<<=i5=j=e=a=<}y$crKÆ\$h8亽RDj$
+`ՅL/U!;<}#ۑj:|23<p<z_='=-s< `= <%]<4<_?9ԩ:W]o=!^ 99N;z5[=*Υ=U=>J=ˬr= 佼S9p&"#+C3nH<lo<h =L=kR==dZ=]=&='=A%<:}Ͻ EcpHA2;7ꈼ E󼮐.;&⼱V}Q3/t
+k==j>8>?!3v?D>ŚI>sPb=1=< 8ͅZ߂vtJu#;#W8񃻒ܼ8*׻鏮tT:н$ ݽPEoh,Qh̽wT>XRb7y<!e<2<}<e{m7 ཉZݽ\< 9< :<=&C{>۽Oƽw{#W<
+=:f=N=T@r
+c<l<;Y<+;J<b;[I͐d9`Ү9%:#S0=#>>&2>(G!> @=AAeżn_&Ѽ#8;1<e<CqN:hg< =Y=,==h= =H=cm=5h<z;njP|dBk4$f
+,p2<G^,݂Ѽ!BqHu=N=Uȼv7?'D˼nĽ3fS<=]W==c@=',<:͇˼>GS,ORݼ5l1軆嬼RgN5
+Uُ+|ؼ;HwHPF~pMהļV&d ޼ ּ5=og;֊<~c<$<ߛ<<a9A;0<:p;c;::` ֆ,ƪ&y<9=->x>D>C+>g=l
+=$Vk9BF+7k5P|Lk
+3>==Cu:
+HڗK(!DT<L<c9<]::<RU<=vw=q=1=N=r<;N`Z ʗ:Ç3:-s\|<<lo<Ɵ?<u)UÖ0r{-CH>*mϿ*׽IOüN2<=^4>'>>>>]? L>&8>]>=G'\<+ZlShxʼmݼ91h0; f<Ir<vs<G;
+L OhR);v<Ր=CS=,
+/iMc$<=s<z=j<=J<ν=XD<<;;~C='QjOx~;Dt<.;Ɋ<=V<=5==>F>j>@>>_D:>
+7;RD<ڝ<6:Y;sF<Q<l+<&4<~(=
+<<<2c)`;Q;;;'(3_P(7'PU2ͽN/@;XݽkUd\g}<
+>iE>hu>o>,>>MA>3==h=v=԰(s+eRXb;/K$O"ev:S;ovjr$Xӽ&J!m𼟔߼Xl(9K;;L; qBetba+(L^%c:OL3I_:9x;a|L[!=6Ď
+Gj{8m::O;є;!l;)v$3%/#z+*B8LPaQz) RYI:k½2伅?:]/<C< ]ǻfP6vʎD5e(Gw1`9Y/M ȼugB<\<P<><*<j<s<c<k<d<~"<m=̺<S`;i_U%ޮ&8}$ֻ-s
+ử;qWwWNu^̽cE署O޽d,ȧV<!\<NH< a<<ɸ=%mi=Ic=El<<zC;`;=ڠYrջ,d0Ŭ d%ѽ2GrѼ¯A:q
+>z~t߽_ڽ `ƒ„ ;t=-> C>y&s>a'a>o>K=,Y>{>t>0p>*e>#M== [=4iTռX zϼnV
+ݽ!
+4<M<{H< <o<<w<ĻM%mz 8
+N(\r.I 9μ+J՞.:-;A<F2=g=#<<>M<Y<Rc<=B<0ܼ̂a༴sDNX~*`YԼWPad.<'ߔܯꀫ볽[`ogk~|ܽ\GCGl/-#y2A6(QE.5IB%rJ_/ǽILm; LJ4n @Mb$֢Z*z7m;
+Z#.k;k[;Q<,=%=(<x<i<H7<><B#;0K<<}k0=T<Jt<=}ܻTae'?rre63L0HDͽ%I\ *$>&?* 쇽w3➼"r;8:i YRM=@)@)]M<x<T<`F=kv=u=# =CG=^<-=- =R=_<P)->&ͼvJƽ0(:f 6:Q(̱;
+e|CT|[Pc%#Pv<}=.==kz=4bz<=>==뛐>'>"2=)M=9:ݽ~H IŽi;IyBǽ&ډŦFbfSP/rս&Ԝ˼7{/
+:=G< ;;ek;;JAh"fwJ)I=)+L琽^P|2Kt.:z-6~R]2쵽pD0b(6>rǼr>} zr 7Z*-EBI)jE4necO+15(%C*WVo)O8E7&t|=)K s[܌v=ļ̽7q蘼*"?nHX;K%<^0= ?~=&n(= #2=6=1=<<Mg;j4<<n%<4><<Mv?Y[~"EXw"}BtF:DiW!?龽
+4H
+Q^'aY:E
+l;0<aD<yS<ǺoVPԻ<ro*ֽ
+~DEZ7WPVOȽG ?m h:J9'G?M#dۻ2'XG &ҽ
+;q<Is,;elN:z伒2)T$Q
+Mⲽl;"-Lˈ<:=ds<Õ<`<oY< >3F>~sX=>==r=X2==o&=#=@`=.S Gquh*f80z7Hœ<-uX:t:ek ZqzBo}콆ͽD%Sλ5<5<
+<h<sNyƽfܼuv
+<U
+1>?>j>;ai>e=ۋD=R:wz]q::襻^PLIcmQ0]˽WTż}򒼂Sn *f~bcּ)<r<y;B<ފ<Wz==M}>/K7>[6>]q>>=ŦE==l> =ڙ={l=<VJw~zi؊s#iŽ*TU4g<5`=_)=i==<QksTuؽ2"l˽߻;<s=`*H=A8=%K=;Z < ;w#l/ABRҽOFwüDֻgX#wZ^a1i磽k^1q1,vwx޽R׺1%ҽ*? Ja-y;d"<=WFQ>8 >>1>>yX>?@=w=:<<;&!:`F+/$◼?꼑uax <l:=QD=_=<"<k;a;;><<մ<;qBFOc&T<<P<<;p;;Nq켔8;S6Udܦ
+waf3k𧽄$rEӽR\e(7|K2t`q\kw/़LuF"x:;6fUjMH)ϼT*<N=_
+=,>==>D5>o=;==C|du'ǯ;Jw;yL6<NZ8c ż[ż aQfż:
+
+&K8m`{*ŽH:<U=*ڄ>:n>:|>}q="=y;WKXʼ;V_"躅c<Ƣ7=4= ==K=9'<?f;;3dRt((ࣽ>#[L0t~)̽`ji9ýa_ㇽet) + NȽYq;˽u2|'aO Sfa~/p꼝Ru<DW=K>T?eC?'?>Y>W>>=g<;=Mɼ4O<#o=8on}K<8l=@=v?=Զ=="</<V<L7<#< ȼ_9u%:<T='<(w<&;R;z;:o<b;ok޼/ ~c[qb2h't޼zX ;ݾ0;-׻H)owqY<_T
+^{3;O;K;"bԼ]TL)S?;4<y=7=C==c=1==Gb;ƭ q+F.<er=,=b^=y<9}BA**ʼ.1F@k=gO=#C<I;%;;9<?<=T7=Q>}>*=ۈ>AO>>^)>0I>" =H=<F==J=" <uZwaToŤ׼Ġ< ?k=:>(>T>>4>W8==Se<y;́λ+ |:=>=_=o>=3=;ӻe*<>St<<*2AeUJnؽ jŽ
+|O(ag%7A*7O'4ͽm.wu6<j=˝=f$=*,="X=95-=
+6igHDBTP@^Bfy>3C98~&;͘<%=O= &==m=;==h==X=h=y{==l<<sEԼ3^VEj%=='<;c E;O?;<K;a<
+=B=
+>81&>D==+=>5Β==l=V=<<|;<n<mՑаQ"(>f3-=#=>I>>T>>.> $=bv<fa<;u<6z<G~=+V=\l=:h=i=r=7u<7;n;/(<4<< s1ٽ.ԽgܽPս;ܽ#ƽN80yeVx^fsD9>WY?˧=<R<=$>)>>6>>KA>>< =e=?S= 9\伹漗M nK2&Uh6Ƚ|Xa޼k;BB^=*3==䝌=!9=y=v< = < <D<=|Gӫ{9W<T2<I==Jd=J<<)<fC<<C<*z=G=!=&<)eT<=U=5===h<:%G^(
+"y
+p̼v6((:Qլмх\i};ͷ+&;$ ˤ;+1<]=]=+>p>|.=ޥ`=&=AN=pW={=j= =
+< ;3<N̼ޘ꼨μ<B<ݥ0=+:ս&};"Q=^<|L<^<,s=>
+<8= Yd<K<<,=&;"[ؼY켯k1<!}=3>>1h>=#<\wM~;Յ<<:2BJbc3ROt;n;TZ;|36A;;P<2<=C=s>!R>XA>N9>)d== N=%%=
+u<1=&=4$<V<>`uF.jMѼ>cͺnN^{;Z߼H<4;ʈ=+̽<<<1=v>,Vۼƻo;NMbVo].Z6GمTf^dzڻ= =H=>"D>m>>>
+<I<GUּkpֽ괽6YRrʼ2⻵޻w[Ya½wSUY[%8C([5]x}1Q<_;=Xǻ>mIط7ٻ"Q<Y7*N׼h 9;3:K;e<F=^
+bcס@K۽0ν὎@i69B3 Ƚ2"_غ R<=i=_B=l -==>x>F4>=Im<B= D0<Oݹ||
+\:=ϼ20u7Sw c{|h1WX悔 {헫=z07+ʽZƂ;ἙU;!+<H <9t:G=b=>=J<hL<$<<<x=n=>P=?:YX0'R;c
+ ȼYkK<x"=6h=$Q=<x<Y<u=Ky=2 <Ʀ<<;ƻ_ +꼜sL <@0<:<'<'<ž<<ސ <=!=k<[;ՙ<7
+;U;*{;;S<'ē=5=-fs=*g=%<λ4i;[x")bAy;|;jqb;[5;8;ɻ= K:<a͈=w
+> s>D==n<z<5 <8< <[ϐR2/=2͖T<9<Z;@ʻۻ3;7<f<:Z<=)7=dn>?p>c>>>;>]@L>q5=
+ҏ<0;$<
+<:͵<q;cs<':j7{мGR=[û
+=<B<Q<%=,-=$= =3i=l<R<pUL]~,(x.۽^f.%x^߆o㸮<.=^=l((=`1@=;:e<<A[=;<O<k<0ٸ_W/4)ql2~{nȂʽ ?ZxUҽrTֹ_C3Q4Lxf=l=!Z=~<zw WY-j<|=Ѕ=J=D<=Jn= Q5<=<M===^O?<ɸ;W<;v4M9|+<c<ņ<%<B;l:\<o<B<o<\<};ΓR$hhN<9<h<=\/=<C\<<:;=I<7<9(<<|;Y <q)<Rm;R;ҩ;s;<0<tf;E7;U.eWxG6X aMw'! :#;Șc:u߮ džӼ)<q=}=="=L<S<e<! ;k;}:zٻ* uMbJˉ&o
+>>j>N>M>$=tM;WbJ92ځt O$y3=;&_
+yU;u<<=?<ٿ=Q=-?<l<=/=!z*cpD]l]#=R8ݽ߽π";r<oa=K 3=2:<Ŵ/8t; <O<53QM\C +
+E\VB߼.  !,PϽ쌽Nl^:.0t[x=^=zx=<Cɽ);m΀%<QL=yq===Q=:[=n=O=>==1=c<2<<,cL<]<?-<n#;pF;;!:*.c9;8';\< I7y1Wzֻ!Z+<= 9==3=a< ;PKiqm^G< i6 Eeۃ?}n;R
+9 Ļ7A))ؼb;a)<=f=[T=F<G<֮<qHd;%4M+F0<nP< 0Ɨ|s>ǽ8
+jRApRv C2=>J>3>}E>V[>>{9
+nO:<7:u~/= =ȼ
+Yzsiɼ6)7[:¼a?PTIŻU#;$<o<6ϼIY0A,y' H0C9}K&U 5H0=޴;ѻռY!7="=6Xw<ٍ<T8j=PxjwZ==V==j<<=,K=3=x=i=4<ݫX<<&<i<~R<;u޼4ἎBλ;KC8;co)ݼ^.B;x;׿;@=+M==>|== :`lP kҊ$ Ƽ<։.hq; #<I'<]4\<9;o>'Tg+mۻr0u֩;ά:<><=h==?=#S<c;mD;;Ӥ8Hթ;
+:z;<-x=1=&2O<n{<Aia<*Qo/ =n<P,/ۣW_:C8ֽ, 5m̼Gn !޼2nA=)>>@>|7>w>H$V>M`==ǐ6 -暽 ټy o&;I7L E@[
+1)~(FΫ;}
+l=2=L=7>; 4=1>k> *=6=>=1S<=D=#==v=aB$OR-HAiq-=^Խؽaʽ6+4F~</= H<[<ƃl?ڵ"ͽ[;+3;C~<7M%<2N9#ʼ񋍽H=J&_%v<:!h;9yoL ;?SS"rȼ9m8#Cdg$ȼDjcܻ6FX=& =;,p<[=Wj=W=;=zl<P <I<,V;;L<;;Dċdy+/-m :İ 1o rT;;=ٳ;5<33K<=s=n>(%==A=2.:x9jzAV+1+
+s<A;h:|;Sq;^K<1\<'Hu<.562
++òѼVf?=+6:ԽIGG^*`TE ݼTF!<0J=K>n4>u>J%>*#;> ={=<3\< <5029\wyNOepż
+Z}3K(x/,<-<<<::9{#r";!g
+‚4;<<J_=U=<;]"46N-AB" eX?sѼ>JU`
+:ӼD@ۅIA ߼ ػMӻμF'V_<aW<v<~=ٹ=S==R=|%=R=I<`=D='Ej=:H;wcܼ@-FqoD};<!<w<?= = a<%w<EB<2==C>6>0>'S==@~;XTI(Fͽ:h<l{0<I;!uQ69<s_;A<E<<i<Z<d"<YF<<ȶ<z<s=D,==np=^]="<bI):ݼ+zf;3G;_ARZ
+yD<;{;O`Իм&WEDyԔX;h4Nk^׬mf%0nin=w=b=(==f=u=7==<;=I;<;<9oknjJ:Z:O.<j&;[f;O;BL<;)<<=$>=]= /=x8=</<v
+O(=R==@<fϼfv x- o( Suj'|{;2V};;;ڔ<tZؼD׽` =(8F<X=*=)?=h&$Pٽ#Z⽉ ܳ{\oC8z~X"
+k=z==e= ==P<e5 <
+={=b;׼0˼/
+Ye"\;J<~<<E;S8h;<>d=Zc=z=w=uA=/7= <w `<h;<<Sc<Lj=W='w= i<%NI㼸8`a_dw!TI2⼾1f{߼?\|6kټ*X09c۽*;$Aٽs>`m~s} ONڻ(=Wff=*=
+r=q==A=r<<<fz;@<X<cm:o-;_`y=7=*V%;;=;9B<
+=#=<QgN=I<x/yU]+itؽֽ퍱o޽|ғʽn9
+=3=&=Y=rIT<8 So!=ʠO'u<[H<ܙ|=<=
+=, ,> 4"$ DXw60]39=(zA˖w欏奻0q; k<N;[ˉ]YlvD;<Y=J=8<<H<K<= =\>$*>Y
+E>K7>=sA=s~<Xz;g;|;Bؼ.'3^<{<V<5<AJ):J<O=TO=>1>>5W>1>y={<S#F#c(ٽQ!C<<<Bu~: =?===z,==5y<X;xq ;0u3;w;7;Ԗ<<{5,V`lZٽ_ټż5Y{˽^׾{(|I|ѽv+ؽ
+w\(p@uTm=bp=<0<ܮ===^b=j<n+ؼIH;6$5<!<g,7=̻qi⛽*
+~ȼ#q NzG*6+&/==U)<&) ޤx[J׺̎Ҽ\M`oΙʀH˶ᅟ罿׽ƽQ;߶>0Y>Kg>)>ӳ>PH=Y=;k*<&3c<F
+<ˣ\=.>2>t1>zf>==  =$&==04<Քcduۼ':d4o"_<=VH=> ~>=§=<;(;Žd\˽ ⼐'&g;zg;tA*4;
+
+% spɻq&ݼuXսPr;xl<ٜ<'<<0.<=~=Ҽ=p=CI^<CJ<~c<,<3#PP;SX+Dҽsーp.5X`wVPȽ\*}=!=
+ *(l7p̼:⻥Ӽ-#¼0Oڨ| 4R
+0< <Y<
+<ػ[\0- Ty7:2{"c,(F0dNB@_}]s4u^h(~`)+;<l-8:d=ƻ 弢ct` &<vaٺ)Ey0= 5=Ll=<'<~=@@=zQ==/=RН=?<R<<Yp<_;<%>(8t6$מo3z
+オb=]j<:J
+ͻ<;o<<tr8=
+FL>!
+Pټ/ռ;W<Q< $<{]f
+;齜:̯<D=K=a# =>>Z>}>2>K(==n;I;)b<@?=<J<t6bͺYؼծu{R*s-''-;1<
+;J)V7k$=>@p>1=l=/;M`ں:e<<2ڼPuO ý
+X+=<>=fp=UL=<;p';+==4>$f(>P_>pj>~>[>}><A=m=
+mPh1>>Z =c<;׹<<;<u<(<,]y&p,?x]i0b
+r缇ɼ*Z=G18=r-=٫==q>'=s=Y=t==}>/:>=[!U=%=ݷ=B|==SZ=h=<c</$;q;Ɯ;k
+l7<ou=<
+=2=u=],u<;+/޻d==N
+=);i2<p<;k<I<0ON D[VKbNJjSx kupV<F<L<Q<s<a==b= =!<VջYP̐x@żwo׻7<Y+1=LA=lW=O<XȻ[c%{]<ý+̬}s\]/:׽YpD|Ei U><
+2 A㼽7;@<
+<_-;A8o2'ցQL #KXKȽ/;ѽҳ?7)ơ.ɽcU;1Hv&Լw˼"E;<n=y=ˌ>6= ,=<9v,;=J=>S>i>g>M],>;==כU=<[L<+<F<}=6<t=y<-Ҿȃ)6q;s|ƽ<%$F* tν|ts=uսI/r
+鼅~:g8u_+<Z<0Ի/I~ʽ
+U="=Bi<(Xս0BY~&;k;&0[͜Xl]^gXsR.j~nݽ/L}T@"<s>>_z%>&= =bn,wQȽ޼
+=*=<˸ɽ3D'z)648j؟|(mZ׼Ycp;z<w<`3=5=a|{=>=`$<d<1);~ ;><Z<CѼJj1!#E('8Ƚ9
+Sy\,m;6ԺI9nH9)_֚iɼӷ(DQ'L.~s7WIf;+DI<Un<=ld=V=Fe=L*8=V1;͘CC:Ӽ)Gm+<"6=F=q ==DMm<%
+]3Mo:v,%
+ځm k;7)Qr:ڼDӫ<A1<=z =|B=A#=U1<`"P;&if׻fn<;3{> ?(1hlt 8#,şqt-9૽ cר/dT潛p O H㼄a{R#Ӽ
+MD";h=z> >Bw8>1>e>8 =]&;;Cw˼!='===7=Wl=,=o,N~:e `̽ ڽYk27e=@yqzΰd,JXM:LP녽7.(N4V<j;`&>q>d>=i =wB=`O#y<OIWX,i*ý9+u<tcq#JJCRQVʛ;i=</)=gr==U=Y=<A;R8xwx{.x_Rb٪4B 5
+ЧD7 b 
+Il#}Nl\Ҽ9&[$p1ͼ<!<yg(=N=
+=ÌT=x=8H<y<<Z< <<0<$<F<x,&F
+p۴ V&BR3ⷽVUFy9S@AA^<t=2$= 7s=<<5\BȀy\<=RI>t>Fo>z>z>1$=\;%cHI=J=n=v=ʖ=}Q=J<l;;C)q\~+ нia0٫{{N`|jcO<m~y)i0G=AYӽR쳽{"GE>7>?>q==W>#P=J;[vI84:~( +U6g7`4~,r տ*߼tC<AW$=y*==_<i<:P<QEnW.k'-]C3>. e)\Sj^
+;
+<T0<@
+hz4Ѽ^SC=@=-D@<C=3=!<q>ZpR1J= !>&>]>,j>]=>hy>!0oT兽IqK:U/Yz"w^a𖽋󧽘uI'"쓽-;mq< <9%i .R* čT?=AN튽<Pqr=&{MjѽcS%)e|
+*B'<=˽6~,7T? j1]ɯ<ܻ<E="HP<:<#<Ui<Il<)<;MJ<<_<R!<l<6`=Z<פ<C:[ܽ"t.%l!%3yx*8W2(нgwwq/ ;=x=N<>
+4=C==*"<ݼt:1s-_9=u4=\ҼD\6h<bC<=JA=n =`<@^:\ӼQ+e]vSn舽p9ǽ8zGa,S8޼ RüA93=S={–=ͣ=.a=g%=C;澽z%wӽ >b0==Ec>,4>>X>^M><e=1p<򽉫8ݸk>?zIǽY_e#ņ}9;9Pd(j F9䝼s.N꾽a{l:zyѽj a&`uzUz{R^FݽIǽ<:}9样)<=B >p8><==U=O|H
+ǼV;V< ;3Z;G69C~,L"kVEX!x+f2$Հ;
+St +Qνѽכ/`-2VDV:<<g6=E=J؄=W{<
+Ⱥ;)#{;<;8KARX&
+g~ebd!ż)s-@<.9!dV,
+3l}F 4 ý#)uĽ9hgBrdof$'9<^H=EN=I=9=h8=Kb;_;!:'.:t; <rp='}o=\=T===!<Jc̡,Hټ 锻 rű<1<<<`n|JhW"(k;<Z=<=+>3>?h>_e{>g==<s8zl)ܽXZl֜NXVI"ټ)WNJ*O;W;]kj nPao}QD xwe:<g!<Je֔;;SA?m<a<66<ܵ<pZ<^:=<F=T=9==U=4= <Y>5K=ú>DvqEoƾG}̼Žy1";H;/t:[/$A; < <[FmF!=
+"ڙ=nf ;]`F{?g$t"7Ts.D~w=D&t>V>N>w>>>=S+溝7MBcM#@("?|k;8<<<<tT\4]ݠ7熼S~T:!;}JmKNSW&!rNջw7ni:M:̓A=Y=a===T=U=o; ;B:Wv<Ou<=-?===<S=mr=I";\N:_P;"<{ϻམ=43=5W=.)=$9n<@"ݻ~8u5P<Ć=\=P>|>>/?>|>wۜ>=o;ѼR'N!ǽ҂ ڇ?ɲj2_LSJ0_. e T;U;h:;UN~9Mm KEӽmսCݼVp듽iQaɼ?qb< <p</7ѻ<$<_<<iZ<ȑ==N=~=.t=r=}=l=e<QIཤ⡼pbnRYXOja㽒.v`C>m<i;lO:ʏ;+'ż+c2UUٽQn{=μ@~Q2uc]
+F{/7<*F=?3<Q=V<Fm<C;cȺW3-vC彀 #W|dc6n9>$=m='A=y=65A=*2=<a<<=W=w<%<&=%=l=v2=qe]=#IGK
+p_}ԼkϽ
+\1ٽ2(;S</mY<<O7(<#<[<J\Ye.RQz;<-;<q:,-<3k!"Bx:=/I>>>->>L>$> =>%<dcjAl Ub׽xJ<=aC=ƅ>Z=Kn==v<<G=Rp<nj:#x2 u\kp%M6-Zͽ4Zi4=G= /=o#==SZ=H`=='s=M=t=x=};F*<h*=
+4<<c<f<fO]/V]$p"ϼPؼhfARv\8ϓ䶴˒֨@'ꗻ
+=+E=)=R=]x<] <6> == D<۪<eM<vT"(<Dc==Q6ݽܸn7r <n%Yo4ֽ75z<<=|EӼ—Nֿe=ف"഼k|y;}:I<=
+†q(v
+)Ѽ.,q:@<`K=AdU=fc==S=ϫ=J=$N;sl0+?l7 'J׽EȽkwٶZB,'K;rM<.:=eO= r#<=t=4
+'=n=oX<Ῥ<<:l#:;<<t<c;76:xMgwÈͽN'B屽<ٽ@q-uJڋ9UͽYн?zH;мҜ2g ٽQcʬ<c=>f`>}>xJ>]C>m>g;=C=8<kPgDov0[=<@=E@i= Z==[<ֲ<<:al=n>+>"b=]==K:k;q=.=P=v>==Ǭ=Fo"\?:(<hn0Ld<y;m-=Jս'f <$=<B<y==w$===y='a<=M< Ӽ0Ҽ~<: +2<I=K=I=R8=q=?H=7q=o=
+e5Nj0Gh.ҽNd^jp\琽>a9';==n$={j=Ȫ===D= v==[-=t=°o=f=L3=xOK5c<K=|5={<;f&iTpغ]y==V=Q=B==LM<:=⬐<&g=<Ưo=5s<F<nd=Y<A<`F+9򽎲Ὄabɨ, 1)g^&Qy\NiI+exEo%D~x60K& .<'=ȟ>*s>FR>>s=<= M=q<-S磼iW缳ۭ:87"J<
+=c=^=<;Ѽ*"<֞=Z=F<Ek<>>
+>={=?b=ѧ>>MO>N
+$A'p3rJk^hT6;,U;ӎ=!:Pwr̙gqL< $*==;>._">Taq>@> v>V==
+=P=ӔS=F=h;#L7-&,*񦼣 ̼@0*wGa<<=s=Fs==x8==|=j=Z<?V<====:>:j>>&>$>6a=ӻkT_[\\o¼ݽ5罽#Q̽O\<_=jz>>$>E>K>J>=< a[L~&&c_@<۶=&=_"=3=4;=
+;;9q<o={=+==R=L"=q=-<|m<Eg8p/:yg;<s=f:> >6>mi>|>b =CА?FL^|nrk~M7x.C=9Y>%97>?4?|?oɸ?,>>=e=}n<L Q1<=>A>I-B=== F}ĬI%^Ѽ-59lLnFAKPM/>4+pN$:ھ$4=D^==2=5< \ɔ;n<\f<rq>=pK====Ê=S=<:8< d;<g(<
+M<\ <M<Ą<ɘ==< <(/:(P<J</<ă==<<N%<';QˇF`o]@be;le=8=Wo=_i=a@>)>/w> l=2Y#rнѼTٻq3= t=R>=-=tJ= =S=V; =s+;Qq`Lv%_
+Հ=gy==a=Ա=<<;'H7߻ua;L 7<=K>>$\> =ڏ<<o>ҽZ\D0Q0FlѭD+;t=>9X?OC??"??b̫?<>O==GLɼ;
+=0>;">>">v>6u'=6ɐ* zu5K꼽V MCz=w}p}S4!U,M?ɾ*XVyP9t<ż#:T =t=G=zD;!;6S=_=S=
+=}=-o<a<xo?<<l<J'<gL<<e8n;:<]j<`=pUU=\<=/<ht,< qrcl[Όr5N?;(=+==>!>">L>n=*Z ۽)|}/ : 2=UZ=/>>ؾ=s=<R:"wyq^jbEm"IQV' :%5K>KL){GU+:}E<@ =-=;==\3<;;lȼع7
+dPӼ"5&ֽ+aɽB+߼$[w?gSϼSa;d"<mq==9*=,;x':Y] aI|V;W<O<<+s4_I ؽ>|2u{B6WnQq B꽃ڽH}ν)뿜=H>{??;<?ŝ?&~?@?r>s= <PлkJ#3;.< =e>/)>#>V>u A>S=Xp<iffIV1 Lǽy0/[_ ^M
+=8Cs95޽ =J7=|=U`=f7=-<&=Y$=J=P=C=\$LVI7SL;[<= GI<5Ћ<S<=ߊ=/=O<p:M볼p40[3,9rl:<' =W=[=@ ='< )<<Y.J<&Pd||I]a=!A=ߙ>>#>=w==<9T`[ѽo8b+[B-ǽ<
+7t̽XX1o
+<`Q<C<|J&<<1b=/"<ھ3;i Q*7!F+eԼ5ytX-zZ I=QcVL@p?X0ͽI&wk7Mwٯ%W{/+XW]*eż8νOf*5",j 'P[OýcͽsнbJʖ#ۻlV!X=9+>>>az>'K#=s<<a-פM@ed>y/;l<y<<nE;<(Xpo<^=L==<O\=t 02ͽA$,7˼9V^<Jѽľ'Ym=$===,|==c=p[Y4*/G<iub`d K%[=d=
+<uD<lt1U=q)=W<]@<=/%=<\<9;GucK/tNҕ,1Qn <;<X\<fl<s;
+tYlWʽf@Aܽ7 e2 ;j҆Ӻ+Ż{Y>e#-5iٽYjJ /3F͊B08w$/W!Y
+Ľh奮>=!j==N=<;ν0 kֽ%;<xM<[W<v< ;ݶ=,o6nIF3<VG<<<i<<=<6E;^;j4+H)Y',Iea国V>DS;:< ) OÚ꽑B
+=,=T<޺
+{3Y1==X=;"<]b IнQuo=58=={> ==
+뼒6n+;<=h<<:;C<} ǁ^G
+ټtĘZyN W)~:=񼰱׼wٔU;fF7y|韼B.'y&b彡j_^dB%rR%=>3Z>+W>%>I>>Bk>>k>>N==s
+s;c:<a>=%=˕=C=Ԙ=;1ͽ'YxxMXsҽ|ͽфƮ=D(_uF;^6<F<x<<zAH:#vr#ZνrDqEo̽JR⯽w^⽍g;3ӽ
+kvKDR<f<_$?kfL$E2VD2}/=NQ`>n# <*0ݽ@z*FD߽)߽SrIyq1A@׼a':5;_]<+<9lPL], kOؼ ߼=ӼW~kQɺ5._;3+<$B;x; ;y((~){ +6rSWІ!87B`+$<)>.CH>;b>Wf?I?yc0?^J?A$? i>޳>[~>>m=;u<V=,=y=e=΄=L'a;̽>tmMjsyV C1G @֦|ڢ=c^P<pl<=<`=<w &+4v<*ݽ㺽 !58I^aL2gPld~X~2,wҽ<Z[Q9Lϑ}'qRJJ5;#@+¾:⾅9
+~ O`$:N
+4bDͶ]K8 ݱt'/Z* ;-<ںT
+=hs>?@5@@@@"]@??>c>^3=<Q$Nؽ@M<
+#,aּa Fp=ɓ<?[̧$uj<]_<;=݅M__5A3cTT}a ʼҝ [ Ƚw}GЅT=}Gbd<";µ=W==y<%;ڵ<
+:,Md~[$=x<=q^e=Qq?e@@ԳA|CA-A"@@.3?/>;> ~=(;./ ]ҽ56ĽÑvld Iy_:[6Ͻ,w8촚@E;=?<a#<1e}E1Wڼs*ksސ ;<R^;F<<Þ<r&hս?Լ̼`= =KV=->a=(=$=4սVս2" |y^ZD
+i=Xm=(6f@¢n2$<:&uf==\= t&,%ݡ<^|=
+R=6=#]<` 9h'漗u9dHۼz* I;;4p<ID<ۯ<@<'F$q(
+43*9$W=ڼ=XY<#ED0e;
+X==<iX<6<O;}㒼;%]+5j65]3W;vk;/|'Fg;@<
+!L]R^$ !J*CV43z~@n*1`:¹=rjo=> T>Q>?ܜ?ɾ@7?9?>}>[== =O<
+˽S;~ٽV\[wG½ĵȽd걼ë"{= =<)<@U<NY<%=M=O=23-=<yU^/3).Ob:PJ:;v[;5;A<u<_<J}:kD <wZ<f=VA=ս>>t>xb>':%MOF,$־!@<ߦ?a]ZֻEQ=sj=&e=_;b0 o==eQ=
+Ľ'
+ռ|.ɽ νH̽m=ɼûxxlN< <2=s> >]>q
+>x>|>v>><>= =T<pWҽyZ *(y"UB[uŽz%<r<b9V<Yx<O</<~<<+,=Q<`r;獑5Tշ1` ><>j<i<<<Q<<b-;Z;);<$:AUp<)3:茌=U==V*>|P>N=>-V.A6ݖKV򽋄=7>#=Z}!;٬<:D<=U,=V7<=,c=
+O= =m`?=:=-y=2=}=eM=w3%=<A<?!!O&pveݿӼ˧¼4܁ἝWThć wq786i; O=
+=Us>
+3 "p= V ΠH8s):ww:<ڹ= .=u<;[WfGR+´@&E߽{+hr|q"}9<cѼ;G<vC=
+'<<'ջu-Ӽ΀SoS;M3%;<=#{q==r=>\=A=Ϡ=>!%>l>f>=X0غ<=Ct{==
+:j;;=_=2h=1 "/^SuW~Pyw=6=x<A;ėWl<FQ* TL;;0peνld*05_&%KC;
+Sρϼ,ȼrX6¥qWmd`ri$bUϐE#ժm9jG-v b/K <^<=;l;J=+<X<-===
+5D/`p5 :84I]׼qCrż~gsWG`ۼּmkqVƅ`Riӽwн e BQ'齟IUfmXi=S+*
+ҽ!~`
+ֽv!x}n<>
+ߺ>B=ؽ$!eN(<);<A=i?=W+l=*=~=a=i=w=7@<5
+=/=d=9r|=oV==-=kv>
+=0dRqn<U1`dMf+N|<J=!:<c<5=3= Q>>>C
diff --git a/lib/test/data/2chipA.fits.REMOVED.git-id b/lib/test/data/2chipA.fits.REMOVED.git-id
new file mode 100644
index 0000000..50242ca
--- /dev/null
+++ b/lib/test/data/2chipA.fits.REMOVED.git-id
@@ -0,0 +1 @@
+fe2343477504847a56ac369c96fbe2d1b9d835db \ No newline at end of file
diff --git a/lib/test/data/2chipB.fits.REMOVED.git-id b/lib/test/data/2chipB.fits.REMOVED.git-id
new file mode 100644
index 0000000..1d70532
--- /dev/null
+++ b/lib/test/data/2chipB.fits.REMOVED.git-id
@@ -0,0 +1 @@
+8a32e201693484546f631c52e3774dfbf278151b \ No newline at end of file
diff --git a/lib/test/data/2chipC.fits.REMOVED.git-id b/lib/test/data/2chipC.fits.REMOVED.git-id
new file mode 100644
index 0000000..7afd241
--- /dev/null
+++ b/lib/test/data/2chipC.fits.REMOVED.git-id
@@ -0,0 +1 @@
+167e2502832be1609cdde29c52e9dc8c571cc3ae \ No newline at end of file
diff --git a/lib/test/test.py b/lib/test/test.py
new file mode 100644
index 0000000..c6c1fcd
--- /dev/null
+++ b/lib/test/test.py
@@ -0,0 +1,260 @@
+import os
+import random
+
+import numpy as np
+from numpy.testing import assert_almost_equal, assert_array_less
+
+from sphere import graph
+from sphere import great_circle_arc
+from sphere import polygon
+from sphere import vector
+
+from .test_util import *
+
+
+graph.DEBUG = True
+
+
+def test_normalize_vector():
+ x, y, z = np.ogrid[-100:100:11,-100:100:11,-100:100:11]
+ xn, yn, zn = vector.normalize_vector(x.flatten(), y.flatten(), z.flatten())
+ l = np.sqrt(xn ** 2 + yn ** 2 + zn ** 2)
+ assert_almost_equal(l, 1.0)
+
+
+def test_radec_to_vector():
+ npx, npy, npz = vector.radec_to_vector(np.arange(-360, 360, 1), 90)
+ assert_almost_equal(npx, 0.0)
+ assert_almost_equal(npy, 0.0)
+ assert_almost_equal(npz, 1.0)
+
+ spx, spy, spz = vector.radec_to_vector(np.arange(-360, 360, 1), -90)
+ assert_almost_equal(spx, 0.0)
+ assert_almost_equal(spy, 0.0)
+ assert_almost_equal(spz, -1.0)
+
+ eqx, eqy, eqz = vector.radec_to_vector(np.arange(-360, 360, 1), 0)
+ assert_almost_equal(eqz, 0.0)
+
+
+def test_vector_to_radec():
+ ra, dec = vector.vector_to_radec(0, 0, 1)
+ assert_almost_equal(dec, 90)
+
+ ra, dec = vector.vector_to_radec(0, 0, -1)
+ assert_almost_equal(dec, -90)
+
+ ra, dec = vector.vector_to_radec(1, 1, 0)
+ assert_almost_equal(ra, 45.0)
+ assert_almost_equal(dec, 0.0)
+
+
+def test_intersects_poly_simple():
+ ra1 = [-10, 10, 10, -10, -10]
+ dec1 = [30, 30, 0, 0, 30]
+ poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1)
+
+ ra2 = [-5, 15, 15, -5, -5]
+ dec2 = [20, 20, -10, -10, 20]
+ poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2)
+
+ assert poly1.intersects_poly(poly2)
+
+ # Make sure it isn't order-dependent
+ ra1 = ra1[::-1]
+ dec1 = dec1[::-1]
+ poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1)
+
+ ra2 = ra2[::-1]
+ dec2 = dec2[::-1]
+ poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2)
+
+ assert poly1.intersects_poly(poly2)
+
+
+def test_intersects_poly_fully_contained():
+ ra1 = [-10, 10, 10, -10, -10]
+ dec1 = [30, 30, 0, 0, 30]
+ poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1)
+
+ ra2 = [-5, 5, 5, -5, -5]
+ dec2 = [20, 20, 10, 10, 20]
+ poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2)
+
+ assert poly1.intersects_poly(poly2)
+
+ # Make sure it isn't order-dependent
+ ra1 = ra1[::-1]
+ dec1 = dec1[::-1]
+ poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1)
+
+ ra2 = ra2[::-1]
+ dec2 = dec2[::-1]
+ poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2)
+
+ assert poly1.intersects_poly(poly2)
+
+
+def test_hard_intersects_poly():
+ ra1 = [-10, 10, 10, -10, -10]
+ dec1 = [30, 30, 0, 0, 30]
+ poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1)
+
+ ra2 = [-20, 20, 20, -20, -20]
+ dec2 = [20, 20, 10, 10, 20]
+ poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2)
+
+ assert poly1.intersects_poly(poly2)
+
+ # Make sure it isn't order-dependent
+ ra1 = ra1[::-1]
+ dec1 = dec1[::-1]
+ poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1)
+
+ ra2 = ra2[::-1]
+ dec2 = dec2[::-1]
+ poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2)
+
+ assert poly1.intersects_poly(poly2)
+
+
+def test_not_intersects_poly():
+ ra1 = [-10, 10, 10, -10, -10]
+ dec1 = [30, 30, 5, 5, 30]
+ poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1)
+
+ ra2 = [-20, 20, 20, -20, -20]
+ dec2 = [-20, -20, -10, -10, -20]
+ poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2)
+
+ assert not poly1.intersects_poly(poly2)
+
+ # Make sure it isn't order-dependent
+ ra1 = ra1[::-1]
+ dec1 = dec1[::-1]
+ poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1)
+
+ ra2 = ra2[::-1]
+ dec2 = dec2[::-1]
+ poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2)
+
+ assert not poly1.intersects_poly(poly2)
+
+
+def test_point_in_poly():
+ point = np.asarray([-0.27475449, 0.47588873, -0.83548781])
+ points = np.asarray([[ 0.04821217, -0.29877206, 0.95310589],
+ [ 0.04451801, -0.47274119, 0.88007608],
+ [-0.14916503, -0.46369786, 0.87334649],
+ [-0.16101648, -0.29210164, 0.94273555],
+ [ 0.04821217, -0.29877206, 0.95310589]])
+ inside = np.asarray([-0.03416009, -0.36858623, 0.9289657])
+ poly = polygon.SphericalPolygon(points, inside)
+ assert not poly.contains_point(point)
+
+
+def test_point_in_poly_lots():
+ import pyfits
+ fits = pyfits.open(os.path.join(ROOT_DIR, '1904-66_TAN.fits'))
+ header = fits[0].header
+
+ poly1 = polygon.SphericalPolygon.from_wcs(
+ header, 1, crval=[0, 87])
+ poly2 = polygon.SphericalPolygon.from_wcs(
+ header, 1, crval=[20, 89])
+ poly3 = polygon.SphericalPolygon.from_wcs(
+ header, 1, crval=[180, 89])
+
+ points = get_point_set()
+ count = 0
+ for point in points:
+ if poly1.contains_point(point) or poly2.contains_point(point) or \
+ poly3.contains_point(point):
+ count += 1
+
+ assert count == 5
+ assert poly1.intersects_poly(poly2)
+ assert not poly1.intersects_poly(poly3)
+ assert not poly2.intersects_poly(poly3)
+
+
+def test_great_circle_arc_intersection():
+ A = [-10, -10]
+ B = [10, 10]
+
+ C = [-25, 10]
+ D = [15, -10]
+
+ E = [-20, 40]
+ F = [20, 40]
+
+ correct = [0.99912414, -0.02936109, -0.02981403]
+
+ A = vector.radec_to_vector(*A)
+ B = vector.radec_to_vector(*B)
+ C = vector.radec_to_vector(*C)
+ D = vector.radec_to_vector(*D)
+ E = vector.radec_to_vector(*E)
+ F = vector.radec_to_vector(*F)
+
+ assert great_circle_arc.intersects(A, B, C, D)
+ r = great_circle_arc.intersection(A, B, C, D)
+ assert r.shape == (3,)
+ assert_almost_equal(r, correct)
+
+ assert np.all(great_circle_arc.intersects([A], [B], [C], [D]))
+ r = great_circle_arc.intersection([A], [B], [C], [D])
+ assert r.shape == (1, 3)
+ assert_almost_equal(r, [correct])
+
+ assert np.all(great_circle_arc.intersects([A], [B], C, D))
+ r = great_circle_arc.intersection([A], [B], C, D)
+ assert r.shape == (1, 3)
+ assert_almost_equal(r, [correct])
+
+ assert not np.all(great_circle_arc.intersects([A, E], [B, F], [C], [D]))
+ r = great_circle_arc.intersection([A, E], [B, F], [C], [D])
+ assert r.shape == (2, 3)
+ assert_almost_equal(r[0], correct)
+ assert np.all(np.isnan(r[1]))
+
+ # Test parallel arcs
+ r = great_circle_arc.intersection(A, B, A, B)
+ assert np.all(np.isnan(r))
+
+
+def test_great_circle_arc_length():
+ A = [90, 0]
+ B = [-90, 0]
+ A = vector.radec_to_vector(*A)
+ B = vector.radec_to_vector(*B)
+ assert great_circle_arc.length(A, B) == 180.0
+
+ A = [135, 0]
+ B = [-90, 0]
+ A = vector.radec_to_vector(*A)
+ B = vector.radec_to_vector(*B)
+ assert great_circle_arc.length(A, B) == 135.0
+
+ A = [0, 0]
+ B = [0, 90]
+ A = vector.radec_to_vector(*A)
+ B = vector.radec_to_vector(*B)
+ assert great_circle_arc.length(A, B) == 90.0
+
+
+def test_great_circle_arc_angle():
+ A = [1, 0, 0]
+ B = [0, 1, 0]
+ C = [0, 0, 1]
+ assert great_circle_arc.angle(A, B, C) == 90.0
+
+ # TODO: More angle tests
+
+
+def test_cone():
+ random.seed(0)
+ for i in range(50):
+ ra = random.randrange(-180, 180)
+ dec = random.randrange(20, 90)
+ cone = polygon.SphericalPolygon.from_cone(ra, dec, 8, steps=64)
diff --git a/lib/test/test_intersection.py b/lib/test/test_intersection.py
new file mode 100644
index 0000000..81512d1
--- /dev/null
+++ b/lib/test/test_intersection.py
@@ -0,0 +1,167 @@
+from __future__ import print_function
+
+# STDLIB
+import functools
+import itertools
+import math
+import os
+import random
+import sys
+
+# THIRD-PARTY
+import numpy as np
+from numpy.testing import assert_array_almost_equal
+
+# LOCAL
+from sphere import polygon
+
+GRAPH_MODE = False
+ROOT_DIR = os.path.join(os.path.dirname(__file__), 'data')
+
+
+class intersection_test:
+ def __init__(self, lon_0, lat_0, proj='ortho'):
+ self._lon_0 = lon_0
+ self._lat_0 = lat_0
+ self._proj = proj
+
+ def __call__(self, func):
+ @functools.wraps(func)
+ def run(*args, **kwargs):
+ polys = func(*args, **kwargs)
+
+ intersections = []
+ num_permutations = math.factorial(len(polys))
+ step_size = int(max(float(num_permutations) / 20.0, 1.0))
+ if GRAPH_MODE:
+ print("%d permutations" % num_permutations)
+ for method in ('parallel', 'serial'):
+ for i, permutation in enumerate(
+ itertools.islice(
+ itertools.permutations(polys),
+ None, None, step_size)):
+ filename = '%s_%s_intersection_%04d.svg' % (
+ func.__name__, method, i)
+ print(filename)
+
+ intersection = polygon.SphericalPolygon.multi_intersection(
+ permutation, method=method)
+ intersections.append(intersection)
+ areas = [x.area() for x in permutation]
+ intersection_area = intersection.area()
+ assert np.all(intersection_area < areas)
+
+ if GRAPH_MODE:
+ fig = plt.figure()
+ m = Basemap(projection=self._proj,
+ lon_0=self._lon_0,
+ lat_0=self._lat_0)
+ m.drawparallels(np.arange(-90., 90., 20.))
+ m.drawmeridians(np.arange(0., 420., 20.))
+ m.drawmapboundary(fill_color='white')
+
+ intersection.draw(m, color='red', linewidth=3)
+ for poly in permutation:
+ poly.draw(m, color='blue', alpha=0.5)
+ plt.savefig(filename)
+ fig.clear()
+
+ lengths = np.array([len(x._points) for x in intersections])
+ assert np.all(lengths == [lengths[0]])
+ areas = np.array([x.area() for x in intersections])
+ assert_array_almost_equal(areas, areas[0])
+
+ return run
+
+
+@intersection_test(0, 90)
+def test1():
+ import pyfits
+ fits = pyfits.open(os.path.join(ROOT_DIR, '1904-66_TAN.fits'))
+ header = fits[0].header
+
+ poly1 = polygon.SphericalPolygon.from_wcs(
+ header, 1, crval=[0, 87])
+ poly2 = polygon.SphericalPolygon.from_wcs(
+ header, 1, crval=[20, 89])
+
+ return [poly1, poly2]
+
+
+@intersection_test(0, 90)
+def test2():
+ poly1 = polygon.SphericalPolygon.from_cone(0, 60, 8, steps=16)
+ poly2 = polygon.SphericalPolygon.from_cone(0, 68, 8, steps=16)
+ poly3 = polygon.SphericalPolygon.from_cone(12, 66, 8, steps=16)
+ return [poly1, poly2, poly3]
+
+
+@intersection_test(0, 90)
+def test3():
+ import pyfits
+ fits = pyfits.open(os.path.join(ROOT_DIR, '1904-66_TAN.fits'))
+ header = fits[0].header
+
+ poly1 = polygon.SphericalPolygon.from_wcs(
+ header, 1, crval=[0, 87])
+ poly3 = polygon.SphericalPolygon.from_wcs(
+ header, 1, crval=[175, 89])
+
+ return [poly1, poly3]
+
+
+def test4():
+ import pyfits
+ import pywcs
+
+ A = pyfits.open(os.path.join(ROOT_DIR, '2chipA.fits.gz'))
+ B = pyfits.open(os.path.join(ROOT_DIR, '2chipB.fits.gz'))
+
+ wcs = pywcs.WCS(A[1].header, fobj=A)
+ chipA1 = polygon.SphericalPolygon.from_wcs(wcs)
+ wcs = pywcs.WCS(A[4].header, fobj=A)
+ chipA2 = polygon.SphericalPolygon.from_wcs(wcs)
+ wcs = pywcs.WCS(B[1].header, fobj=B)
+ chipB1 = polygon.SphericalPolygon.from_wcs(wcs)
+ wcs = pywcs.WCS(B[4].header, fobj=B)
+ chipB2 = polygon.SphericalPolygon.from_wcs(wcs)
+
+ Apoly = chipA1.union(chipA2)
+ Bpoly = chipB1.union(chipB2)
+
+ X = Apoly.intersection(Bpoly)
+
+
+def test5():
+ import pyfits
+ import pywcs
+
+ A = pyfits.open(os.path.join(ROOT_DIR, '2chipA.fits.gz'))
+ B = pyfits.open(os.path.join(ROOT_DIR, '2chipB.fits.gz'))
+
+ wcs = pywcs.WCS(A[1].header, fobj=A)
+ chipA1 = polygon.SphericalPolygon.from_wcs(wcs)
+ wcs = pywcs.WCS(A[4].header, fobj=A)
+ chipA2 = polygon.SphericalPolygon.from_wcs(wcs)
+ wcs = pywcs.WCS(B[1].header, fobj=B)
+ chipB1 = polygon.SphericalPolygon.from_wcs(wcs)
+ wcs = pywcs.WCS(B[4].header, fobj=B)
+ chipB2 = polygon.SphericalPolygon.from_wcs(wcs)
+
+ Apoly = chipA1.union(chipA2)
+ Bpoly = chipB1.union(chipB2)
+
+ Apoly.overlap(chipB1)
+
+
+
+if __name__ == '__main__':
+ if '--profile' not in sys.argv:
+ GRAPH_MODE = True
+ from mpl_toolkits.basemap import Basemap
+ from matplotlib import pyplot as plt
+
+ functions = [(k, v) for k, v in globals().items() if k.startswith('test')]
+ functions.sort()
+ for k, v in functions:
+ v()
diff --git a/lib/test/test_union.py b/lib/test/test_union.py
new file mode 100644
index 0000000..66abda8
--- /dev/null
+++ b/lib/test/test_union.py
@@ -0,0 +1,169 @@
+from __future__ import print_function
+
+# STDLIB
+import functools
+import itertools
+import math
+import os
+import random
+import sys
+
+# THIRD-PARTY
+import numpy as np
+from numpy.testing import assert_array_almost_equal
+
+# LOCAL
+from sphere import polygon
+
+GRAPH_MODE = False
+ROOT_DIR = os.path.join(os.path.dirname(__file__), 'data')
+
+
+class union_test:
+ def __init__(self, lon_0, lat_0, proj='ortho'):
+ self._lon_0 = lon_0
+ self._lat_0 = lat_0
+ self._proj = proj
+
+ def __call__(self, func):
+ @functools.wraps(func)
+ def run(*args, **kwargs):
+ polys = func(*args, **kwargs)
+
+ unions = []
+ num_permutations = math.factorial(len(polys))
+ step_size = int(max(float(num_permutations) / 20.0, 1.0))
+ if GRAPH_MODE:
+ print("%d permutations" % num_permutations)
+ for method in ('parallel', 'serial'):
+ for i, permutation in enumerate(
+ itertools.islice(
+ itertools.permutations(polys),
+ None, None, step_size)):
+ filename = '%s_%s_union_%04d.svg' % (
+ func.__name__, method, i)
+ print(filename)
+
+ union = polygon.SphericalPolygon.multi_union(
+ permutation, method=method)
+ unions.append(union)
+ areas = [x.area() for x in permutation]
+ union_area = union.area()
+ assert np.all(union_area >= areas)
+
+ if GRAPH_MODE:
+ fig = plt.figure()
+ m = Basemap(projection=self._proj,
+ lon_0=self._lon_0,
+ lat_0=self._lat_0)
+ m.drawmapboundary(fill_color='white')
+ m.drawparallels(np.arange(-90., 90., 20.))
+ m.drawmeridians(np.arange(0., 420., 20.))
+
+ union.draw(m, color='red', linewidth=3)
+ for poly in permutation:
+ poly.draw(m, color='blue', alpha=0.5)
+ plt.savefig(filename)
+ fig.clear()
+
+ lengths = np.array([len(x._points) for x in unions])
+ assert np.all(lengths == [lengths[0]])
+ areas = np.array([x.area() for x in unions])
+ # assert_array_almost_equal(areas, areas[0], 1)
+
+ return run
+
+
+@union_test(0, 90)
+def test1():
+ import pyfits
+ fits = pyfits.open(os.path.join(ROOT_DIR, '1904-66_TAN.fits'))
+ header = fits[0].header
+
+ poly1 = polygon.SphericalPolygon.from_wcs(
+ header, 1, crval=[0, 87])
+ poly2 = polygon.SphericalPolygon.from_wcs(
+ header, 1, crval=[20, 89])
+ poly3 = polygon.SphericalPolygon.from_wcs(
+ header, 1, crval=[175, 89])
+ poly4 = polygon.SphericalPolygon.from_cone(
+ 90, 70, 10, steps=50)
+
+ return [poly1, poly2, poly3, poly4]
+
+
+@union_test(0, 90)
+def test2():
+ poly1 = polygon.SphericalPolygon.from_cone(0, 60, 7, steps=16)
+ poly2 = polygon.SphericalPolygon.from_cone(0, 72, 7, steps=16)
+ poly3 = polygon.SphericalPolygon.from_cone(20, 60, 7, steps=16)
+ poly4 = polygon.SphericalPolygon.from_cone(20, 72, 7, steps=16)
+ poly5 = polygon.SphericalPolygon.from_cone(35, 55, 7, steps=16)
+ poly6 = polygon.SphericalPolygon.from_cone(60, 60, 3, steps=16)
+ return [poly1, poly2, poly3, poly4, poly5, poly6]
+
+
+@union_test(0, 90)
+def test3():
+ random.seed(0)
+ polys = []
+ for i in range(10):
+ polys.append(polygon.SphericalPolygon.from_cone(
+ random.randrange(-180, 180),
+ random.randrange(20, 90),
+ random.randrange(5, 16),
+ steps=16))
+ return polys
+
+
+@union_test(0, 15)
+def test4():
+ random.seed(64)
+ polys = []
+ for i in range(10):
+ polys.append(polygon.SphericalPolygon.from_cone(
+ random.randrange(-30, 30),
+ random.randrange(-15, 60),
+ random.randrange(5, 16),
+ steps=16))
+ return polys
+
+
+def test5():
+ import pyfits
+ import pywcs
+
+ A = pyfits.open(os.path.join(ROOT_DIR, '2chipA.fits.gz'))
+
+ wcs = pywcs.WCS(A[1].header, fobj=A)
+ chipA1 = polygon.SphericalPolygon.from_wcs(wcs)
+ wcs = pywcs.WCS(A[4].header, fobj=A)
+ chipA2 = polygon.SphericalPolygon.from_wcs(wcs)
+
+ null_union = chipA1.union(chipA2)
+
+
+def test6():
+ import pyfits
+ import pywcs
+
+ A = pyfits.open(os.path.join(ROOT_DIR, '2chipC.fits.gz'))
+
+ wcs = pywcs.WCS(A[1].header, fobj=A)
+ chipA1 = polygon.SphericalPolygon.from_wcs(wcs)
+ wcs = pywcs.WCS(A[4].header, fobj=A)
+ chipA2 = polygon.SphericalPolygon.from_wcs(wcs)
+
+ null_union = chipA1.union(chipA2)
+
+
+if __name__ == '__main__':
+ if '--profile' not in sys.argv:
+ GRAPH_MODE = True
+ from mpl_toolkits.basemap import Basemap
+ from matplotlib import pyplot as plt
+
+ functions = [(k, v) for k, v in globals().items() if k.startswith('test')]
+ functions.sort()
+ for k, v in functions:
+ v()
diff --git a/lib/test/test_util.py b/lib/test/test_util.py
new file mode 100644
index 0000000..47a3444
--- /dev/null
+++ b/lib/test/test_util.py
@@ -0,0 +1,14 @@
+import os
+
+import numpy as np
+from sphere import vector
+
+ROOT_DIR = os.path.join(os.path.dirname(__file__), 'data')
+
+def get_point_set(density=25):
+ points = []
+ for i in np.linspace(-85, 85, density, True):
+ for j in np.linspace(-180, 180, np.cos(np.deg2rad(i)) * density):
+ points.append([j, i])
+ points = np.asarray(points)
+ return np.dstack(vector.radec_to_vector(points[:,0], points[:,1]))[0]
diff --git a/lib/vector.py b/lib/vector.py
new file mode 100644
index 0000000..7323c92
--- /dev/null
+++ b/lib/vector.py
@@ -0,0 +1,213 @@
+# -*- coding: utf-8 -*-
+
+# Copyright (C) 2011 Association of Universities for Research in
+# Astronomy (AURA)
+#
+# 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 name of AURA and its representatives may not be used to
+# endorse or promote products derived from this software without
+# specific prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY AURA ``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 AURA 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 `sphere.vector` module contains the basic operations for handling
+vectors and converting them to and from other representations.
+"""
+
+# THIRD-PARTY
+import numpy as np
+
+
+__all__ = ['radec_to_vector', 'vector_to_radec', 'normalize_vector',
+ 'rotate_around']
+
+
+def radec_to_vector(ra, dec, degrees=True):
+ ur"""
+ Converts a location on the unit sphere from right-ascension and
+ declination to an *x*, *y*, *z* vector.
+
+ Parameters
+ ----------
+ ra, dec : scalars or 1-D arrays
+
+ degrees : bool, optional
+
+ If `True`, (default) *ra* and *dec* are in decimal degrees,
+ otherwise in radians.
+
+ Returns
+ -------
+ x, y, z : tuple of scalars or 1-D arrays of the same length
+
+ Notes
+ -----
+ Where right-ascension is *α* and declination is *δ*:
+
+ .. math::
+ x = \cos\alpha \cos\delta
+
+ y = \sin\alpha \cos\delta
+
+ z = \sin\delta
+ """
+ ra = np.asanyarray(ra)
+ dec = np.asanyarray(dec)
+
+ if degrees:
+ ra_rad = np.deg2rad(ra)
+ dec_rad = np.deg2rad(dec)
+ else:
+ ra_rad = ra
+ dec_rad = dec
+
+ cos_dec = np.cos(dec_rad)
+
+ return (
+ np.cos(ra_rad) * cos_dec,
+ np.sin(ra_rad) * cos_dec,
+ np.sin(dec_rad))
+
+
+def vector_to_radec(x, y, z, degrees=True):
+ ur"""
+ Converts a vector to right-ascension and declination.
+
+ Parameters
+ ----------
+ x, y, z : scalars or 1-D arrays
+ The input vectors
+
+ degrees : bool, optional
+ If `True` (default) the result is returned in decimal degrees,
+ otherwise radians.
+
+ Returns
+ -------
+ ra, dec : tuple of scalars or arrays of the same length
+
+ Notes
+ -----
+ Where right-ascension is *α* and declination is
+ *δ*:
+
+ .. math::
+ \alpha = \arctan2(y, x)
+
+ \delta = \arctan2(z, \sqrt{x^2 + y^2})
+ """
+ x = np.asanyarray(x)
+ y = np.asanyarray(y)
+ z = np.asanyarray(z)
+
+ result = (
+ np.arctan2(y, x),
+ np.arctan2(z, np.sqrt(x ** 2 + y ** 2)))
+
+ if degrees:
+ return np.rad2deg(result[0]), np.rad2deg(result[1])
+ else:
+ return result
+
+
+def normalize_vector(x, y, z, inplace=False):
+ ur"""
+ Normalizes a vector so it falls on the unit sphere.
+
+ *x*, *y*, *z* may be scalars or 1-D arrays
+
+ Parameters
+ ----------
+ x, y, z : scalars or 1-D arrays of the same length
+ The input vectors
+
+ inplace : bool, optional
+ When `True`, the original arrays will be normalized in place,
+ otherwise a normalized copy is returned.
+
+ Returns
+ -------
+ X, Y, Z : scalars of 1-D arrays of the same length
+ The normalized output vectors
+ """
+ x = np.asanyarray(x)
+ y = np.asanyarray(y)
+ z = np.asanyarray(z)
+
+ l = np.sqrt(x ** 2 + y ** 2 + z ** 2)
+
+ if inplace:
+ x /= l
+ y /= l
+ z /= l
+ return (x, y, z)
+ else:
+ return (x / l, y / l, z / l)
+
+
+def rotate_around(x, y, z, u, v, w, theta, degrees=True):
+ r"""
+ Rotates the vector (*x*, *y*, *z*) around the arbitrary axis defined by
+ vector (*u*, *v*, *w*) by *theta*.
+
+ It is assumed that both (*x*, *y*, *z*) and (*u*, *v*, *w*) are
+ already normalized.
+
+ Parameters
+ ----------
+ x, y, z : doubles
+ The normalized vector to rotate
+
+ u, v, w : doubles
+ The normalized vector to rotate around
+
+ theta : double, or array of doubles
+ The amount to rotate
+
+ degrees : bool, optional
+ When `True`, *theta* is given in degrees, otherwise radians.
+
+ Returns
+ -------
+ X, Y, Z : doubles
+ The rotated vector
+ """
+ if degrees:
+ theta = np.deg2rad(theta)
+
+ u2 = u**2
+ v2 = v**2
+ w2 = w**2
+ costheta = np.cos(theta)
+ sintheta = np.sin(theta)
+ icostheta = 1.0 - costheta
+
+ det = (-u*x - v*y - w*z)
+ X = (-u*det)*icostheta + x*costheta + (-w*y + v*z)*sintheta
+ Y = (-v*det)*icostheta + y*costheta + ( w*x - u*z)*sintheta
+ Z = (-w*det)*icostheta + z*costheta + (-v*x + u*y)*sintheta
+
+ return X, Y, Z
diff --git a/setup.py b/setup.py
new file mode 100755
index 0000000..b810728
--- /dev/null
+++ b/setup.py
@@ -0,0 +1,45 @@
+#!/usr/bin/env python
+
+CONTACT = "Michael Droettboom"
+EMAIL = "help@stsci.edu"
+VERSION = "0.1"
+
+from distutils.core import setup, Extension
+import sys
+
+try:
+ import numpy
+except ImportError:
+ print("numpy must be installed to build sphere.")
+ print("ABORTING.")
+ raise
+
+major, minor, rest = numpy.__version__.split(".", 2)
+if (int(major), int(minor)) < (1, 4):
+ print("numpy version 1.4 or later must be installed to build pywcs.")
+ print("ABORTING.")
+ raise ImportError
+
+try:
+ numpy_include = numpy.get_include()
+except AttributeError:
+ numpy_include = numpy.get_numpy_include()
+
+extensions = [
+ Extension('sphere.math_util',
+ ['src/math_util.c'],
+ include_dirs = [numpy_include],
+ libraries = ['m'])
+ ]
+
+setup(
+ name = 'sphere',
+ version = VERSION,
+ description = "Python based tools for spherical geometry",
+ author = CONTACT,
+ author_email = EMAIL,
+ packages = ['sphere', 'sphere.test'],
+ package_dir = {'sphere': 'lib', 'sphere.test': 'lib/test'},
+ package_data = {'sphere.test': ['data/*.fits', 'data/*.fits.gz']},
+ ext_modules = extensions
+ )
diff --git a/src/math_util.c b/src/math_util.c
new file mode 100644
index 0000000..f67bea5
--- /dev/null
+++ b/src/math_util.c
@@ -0,0 +1,374 @@
+#include "Python.h"
+#include "numpy/arrayobject.h"
+#include "numpy/ufuncobject.h"
+
+#include "numpy/npy_3kcompat.h"
+
+/*
+ *****************************************************************************
+ ** BASICS **
+ *****************************************************************************
+ */
+
+typedef npy_intp intp;
+
+#define INIT_OUTER_LOOP_1 \
+ intp dN = *dimensions++; \
+ intp N_; \
+ intp s0 = *steps++;
+
+#define INIT_OUTER_LOOP_2 \
+ INIT_OUTER_LOOP_1 \
+ intp s1 = *steps++;
+
+#define INIT_OUTER_LOOP_3 \
+ INIT_OUTER_LOOP_2 \
+ intp s2 = *steps++;
+
+#define INIT_OUTER_LOOP_4 \
+ INIT_OUTER_LOOP_3 \
+ intp s3 = *steps++;
+
+#define INIT_OUTER_LOOP_5 \
+ INIT_OUTER_LOOP_4 \
+ intp s4 = *steps++;
+
+#define BEGIN_OUTER_LOOP_3 \
+ for (N_ = 0; N_ < dN; N_++, args[0] += s0, args[1] += s1, args[2] += s2) {
+
+#define BEGIN_OUTER_LOOP_4 \
+ for (N_ = 0; N_ < dN; N_++, args[0] += s0, args[1] += s1, args[2] += s2, args[3] += s3) {
+
+#define BEGIN_OUTER_LOOP_5 \
+ for (N_ = 0; N_ < dN; N_++, args[0] += s0, args[1] += s1, args[2] += s2, args[3] += s3, args[4] += s4) {
+
+#define END_OUTER_LOOP }
+
+static inline void
+load_point(const char *in, const intp s, double* out) {
+ out[0] = (*(double *)in);
+ in += s;
+ out[1] = (*(double *)in);
+ in += s;
+ out[2] = (*(double *)in);
+}
+
+static inline void
+save_point(const double* in, char *out, const intp s) {
+ *(double *)out = in[0];
+ out += s;
+ *(double *)out = in[1];
+ out += s;
+ *(double *)out = in[2];
+}
+
+static inline void
+cross(const double *A, const double *B, double *C) {
+ C[0] = A[1]*B[2] - A[2]*B[1];
+ C[1] = A[2]*B[0] - A[0]*B[2];
+ C[2] = A[0]*B[1] - A[1]*B[0];
+}
+
+static inline void
+normalize(double *A) {
+ double l = sqrt(A[0]*A[0] + A[1]*A[1] + A[2]*A[2]);
+ A[0] /= l;
+ A[1] /= l;
+ A[2] /= l;
+}
+
+static inline double
+dot(const double *A, const double *B) {
+ return A[0]*B[0] + A[1]*B[1] + A[2]*B[2];
+}
+
+static inline double
+sign(const double A) {
+ return (A == 0.0) ? 0.0 : ((A < 0.0) ? -1.0 : 1.0);
+}
+
+static inline int
+equals(const double *A, const double *B) {
+ return A[0] == B[0] && A[1] == B[1] && A[2] == B[2];
+}
+
+static inline void
+mult(double *T, const double f) {
+ T[0] *= f;
+ T[1] *= f;
+ T[2] *= f;
+}
+
+/*
+ *****************************************************************************
+ ** UFUNC LOOPS **
+ *****************************************************************************
+ */
+
+/*///////////////////////////////////////////////////////////////////////////
+ inner1d
+*/
+
+char *inner1d_signature = "(i),(i)->()";
+
+static void
+DOUBLE_inner1d(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+ INIT_OUTER_LOOP_3
+ intp di = dimensions[0];
+ intp i;
+ intp is1=steps[0], is2=steps[1];
+ BEGIN_OUTER_LOOP_3
+ char *ip1=args[0], *ip2=args[1], *op=args[2];
+ double sum = 0;
+ for (i = 0; i < di; i++) {
+ sum += (*(double *)ip1) * (*(double *)ip2);
+ ip1 += is1;
+ ip2 += is2;
+ }
+ *(double *)op = sum;
+ END_OUTER_LOOP
+}
+
+static PyUFuncGenericFunction inner1d_functions[] = { DOUBLE_inner1d };
+static void * inner1d_data[] = { (void *)NULL };
+static char inner1d_signatures[] = { PyArray_DOUBLE, PyArray_DOUBLE, PyArray_DOUBLE };
+
+/*///////////////////////////////////////////////////////////////////////////
+ cross
+*/
+
+char *cross_signature = "(i),(i)->(i)";
+
+static void
+DOUBLE_cross(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+ double A[3];
+ double B[3];
+ double C[3];
+
+ INIT_OUTER_LOOP_3
+ intp is1=steps[0], is2=steps[1], is3=steps[2];
+ BEGIN_OUTER_LOOP_3
+ char *ip1=args[0], *ip2=args[1], *op=args[2];
+
+ load_point(ip1, is1, A);
+ load_point(ip2, is2, B);
+
+ cross(A, B, C);
+
+ save_point(C, op, is3);
+ END_OUTER_LOOP
+}
+
+static PyUFuncGenericFunction cross_functions[] = { DOUBLE_cross };
+static void * cross_data[] = { (void *)NULL };
+static char cross_signatures[] = { PyArray_DOUBLE, PyArray_DOUBLE, PyArray_DOUBLE };
+
+/*///////////////////////////////////////////////////////////////////////////
+ cross_and_norm
+*/
+
+char *cross_and_norm_signature = "(i),(i)->(i)";
+
+/*
+ * This implements the function
+ * out[n] = sum_i { in1[n, i] * in2[n, i] }.
+ */
+static void
+DOUBLE_cross_and_norm(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+ double A[3];
+ double B[3];
+ double C[3];
+
+ INIT_OUTER_LOOP_3
+ intp is1=steps[0], is2=steps[1], is3=steps[2];
+ BEGIN_OUTER_LOOP_3
+ char *ip1=args[0], *ip2=args[1], *op=args[2];
+
+ load_point(ip1, is1, A);
+ load_point(ip2, is2, B);
+
+ cross(A, B, C);
+ normalize(C);
+
+ save_point(C, op, is3);
+ END_OUTER_LOOP
+}
+
+static PyUFuncGenericFunction cross_and_norm_functions[] = { DOUBLE_cross_and_norm };
+static void * cross_and_norm_data[] = { (void *)NULL };
+static char cross_and_norm_signatures[] = { PyArray_DOUBLE, PyArray_DOUBLE, PyArray_DOUBLE };
+
+/*///////////////////////////////////////////////////////////////////////////
+ intersection
+*/
+
+char *intersection_signature = "(i),(i),(i),(i)->(i)";
+
+/*
+ * Finds the intersection of 2 great circle arcs AB and CD.
+ */
+static void
+DOUBLE_intersection(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+ double A[3];
+ double B[3];
+ double C[3];
+ double D[3];
+
+ double ABX[3];
+ double CDX[3];
+ double T[3];
+ double tmp[3];
+
+ double nans[3];
+
+ double s;
+ int match;
+
+ nans[0] = nans[1] = nans[2] = NPY_NAN;
+
+ INIT_OUTER_LOOP_5
+ intp is1=steps[0], is2=steps[1], is3=steps[2], is4=steps[3], is5=steps[4];
+ BEGIN_OUTER_LOOP_5
+ char *ip1=args[0], *ip2=args[1], *ip3=args[2], *ip4=args[3], *op=args[4];
+
+ load_point(ip1, is1, A);
+ load_point(ip2, is2, B);
+ load_point(ip3, is3, C);
+ load_point(ip4, is4, D);
+
+ match = !(equals(A, C) | equals(A, D) | equals(B, C) | equals(B, D));
+
+ if (match) {
+ cross(A, B, ABX);
+ cross(C, D, CDX);
+ cross(ABX, CDX, T);
+ normalize(T);
+
+ match = 0;
+ cross(ABX, A, tmp);
+ s = sign(dot(tmp, T));
+ cross(B, ABX, tmp);
+ if (s == sign(dot(tmp, T))) {
+ cross(CDX, C, tmp);
+ if (s == sign(dot(tmp, T))) {
+ cross(D, CDX, tmp);
+ if (s == sign(dot(tmp, T))) {
+ match = 1;
+ }
+ }
+ }
+ }
+
+ if (match) {
+ mult(T, s);
+ save_point(T, op, is5);
+ } else {
+ save_point(nans, op, is5);
+ }
+ END_OUTER_LOOP
+}
+
+static PyUFuncGenericFunction intersection_functions[] = { DOUBLE_intersection };
+static void * intersection_data[] = { (void *)NULL };
+static char intersection_signatures[] = { PyArray_DOUBLE, PyArray_DOUBLE, PyArray_DOUBLE, PyArray_DOUBLE, PyArray_DOUBLE };
+
+/*
+ *****************************************************************************
+ ** MODULE **
+ *****************************************************************************
+ */
+
+static void
+addUfuncs(PyObject *dictionary) {
+ PyObject *f;
+
+ f = PyUFunc_FromFuncAndDataAndSignature(
+ inner1d_functions, inner1d_data, inner1d_signatures, 2, 2, 1,
+ PyUFunc_None, "inner1d",
+ "inner on the last dimension and broadcast on the rest \n" \
+ " \"(i),(i)->()\" \n",
+ 0, inner1d_signature);
+ PyDict_SetItemString(dictionary, "inner1d", f);
+ Py_DECREF(f);
+
+ f = PyUFunc_FromFuncAndDataAndSignature(
+ cross_functions, cross_data, cross_signatures, 2, 2, 1,
+ PyUFunc_None, "cross",
+ "cross product of 3-vectors only \n" \
+ " \"(i),(i)->(i)\" \n",
+ 0, cross_signature);
+ PyDict_SetItemString(dictionary, "cross", f);
+ Py_DECREF(f);
+
+ f = PyUFunc_FromFuncAndDataAndSignature(
+ cross_and_norm_functions, cross_and_norm_data, cross_and_norm_signatures, 2, 2, 1,
+ PyUFunc_None, "cross_and_norm",
+ "cross_and_norm product of 3-vectors only \n" \
+ " \"(i),(i)->(i)\" \n",
+ 0, cross_and_norm_signature);
+ PyDict_SetItemString(dictionary, "cross_and_norm", f);
+ Py_DECREF(f);
+
+ f = PyUFunc_FromFuncAndDataAndSignature(
+ intersection_functions, intersection_data, intersection_signatures, 2, 4, 1,
+ PyUFunc_None, "intersection",
+ "intersection product of 3-vectors only \n" \
+ " \"(i),(i),(i),(i)->(i)\" \n",
+ 0, intersection_signature);
+ PyDict_SetItemString(dictionary, "intersection", f);
+ Py_DECREF(f);
+}
+
+#if defined(NPY_PY3K)
+static struct PyModuleDef moduledef = {
+ PyModuleDef_HEAD_INIT,
+ "math_util",
+ NULL,
+ -1,
+ NULL,
+ NULL,
+ NULL,
+ NULL,
+ NULL
+};
+#endif
+
+#if defined(NPY_PY3K)
+#define RETVAL m
+PyObject *PyInit_math_util(void)
+#else
+#define RETVAL
+PyMODINIT_FUNC
+initmath_util(void)
+#endif
+{
+ PyObject *m;
+ PyObject *d;
+
+#if defined(NPY_PY3K)
+ m = PyModule_Create(&moduledef);
+#else
+ m = Py_InitModule("math_util", NULL);
+#endif
+ if (m == NULL)
+ return RETVAL;
+
+ import_array();
+ import_ufunc();
+
+ d = PyModule_GetDict(m);
+
+ /* Load the ufunc operators into the module's namespace */
+ addUfuncs(d);
+
+ if (PyErr_Occurred()) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "cannot load umath_tests module.");
+ }
+
+ return RETVAL;
+}