diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..12d5d18 --- /dev/null +++ b/.gitignore @@ -0,0 +1,8 @@ +*.pyc + +.project + +*.pydevproject + + +_build \ No newline at end of file diff --git a/.travis.yml b/.travis.yml new file mode 100644 index 0000000..ad92833 --- /dev/null +++ b/.travis.yml @@ -0,0 +1,15 @@ +sudo: false + +language: python + +python: + - "2.7" + - "3.3" + - "3.4" + - "3.5" + +cache: pip + +script: + - make install + - make test \ No newline at end of file diff --git a/geomag/MANIFEST.in b/MANIFEST.in similarity index 95% rename from geomag/MANIFEST.in rename to MANIFEST.in index dadfdd2..1cffb43 100644 --- a/geomag/MANIFEST.in +++ b/MANIFEST.in @@ -1,2 +1,2 @@ -include geomag/*.COF -exclude geomag/.svn +include geomag/*.COF +exclude geomag/.svn diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..4606ff9 --- /dev/null +++ b/Makefile @@ -0,0 +1,12 @@ +help: + @echo "Usage:" + @echo "make help -- display this help" + @echo "make test -- run the tests" + @echo "make install -- install for development" + +install: + pip install -r requirements.txt + pip install -e . + +test: + pytest diff --git a/geomag/README.txt b/README.txt similarity index 97% rename from geomag/README.txt rename to README.txt index 7276621..b591b81 100644 --- a/geomag/README.txt +++ b/README.txt @@ -1,5 +1,5 @@ -Magnetic variation/declination ------------------------------- - -Calculates magnetic variation/declination for any latitude/longitude/altitude, -for any date. Uses the NOAA National Geophysical Data Center, epoch 2015 data. +Magnetic variation/declination +------------------------------ + +Calculates magnetic variation/declination for any latitude/longitude/altitude, +for any date. Uses the NOAA National Geophysical Data Center, epoch 2015 data. diff --git a/doc/Makefile b/doc/Makefile new file mode 100644 index 0000000..c5ff2a8 --- /dev/null +++ b/doc/Makefile @@ -0,0 +1,192 @@ +# Makefile for Sphinx documentation +# + +# You can set these variables from the command line. +SPHINXOPTS = +SPHINXBUILD = sphinx-build +PAPER = +BUILDDIR = _build + +# User-friendly check for sphinx-build +ifeq ($(shell which $(SPHINXBUILD) >/dev/null 2>&1; echo $$?), 1) +$(error The '$(SPHINXBUILD)' command was not found. Make sure you have Sphinx installed, then set the SPHINXBUILD environment variable to point to the full path of the '$(SPHINXBUILD)' executable. Alternatively you can add the directory with the executable to your PATH. If you don't have Sphinx installed, grab it from http://sphinx-doc.org/) +endif + +# Internal variables. +PAPEROPT_a4 = -D latex_paper_size=a4 +PAPEROPT_letter = -D latex_paper_size=letter +ALLSPHINXOPTS = -d $(BUILDDIR)/doctrees $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) . +# the i18n builder cannot share the environment and doctrees with the others +I18NSPHINXOPTS = $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) . + +.PHONY: help clean html dirhtml singlehtml pickle json htmlhelp qthelp devhelp epub latex latexpdf text man changes linkcheck doctest coverage gettext + +help: + @echo "Please use \`make ' where 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 " applehelp to make an Apple Help Book" + @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 " latexpdfja to make LaTeX files and run them through platex/dvipdfmx" + @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 " xml to make Docutils-native XML files" + @echo " pseudoxml to make pseudoxml-XML files for display purposes" + @echo " linkcheck to check all external links for integrity" + @echo " doctest to run all doctests embedded in the documentation (if enabled)" + @echo " coverage to run coverage check of 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/GeoMAg.qhcp" + @echo "To view the help file:" + @echo "# assistant -collectionFile $(BUILDDIR)/qthelp/GeoMAg.qhc" + +applehelp: + $(SPHINXBUILD) -b applehelp $(ALLSPHINXOPTS) $(BUILDDIR)/applehelp + @echo + @echo "Build finished. The help book is in $(BUILDDIR)/applehelp." + @echo "N.B. You won't be able to view it unless you put it in" \ + "~/Library/Documentation/Help or install it in your application" \ + "bundle." + +devhelp: + $(SPHINXBUILD) -b devhelp $(ALLSPHINXOPTS) $(BUILDDIR)/devhelp + @echo + @echo "Build finished." + @echo "To view the help file:" + @echo "# mkdir -p $$HOME/.local/share/devhelp/GeoMAg" + @echo "# ln -s $(BUILDDIR)/devhelp $$HOME/.local/share/devhelp/GeoMAg" + @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." + +latexpdfja: + $(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex + @echo "Running LaTeX files through platex and dvipdfmx..." + $(MAKE) -C $(BUILDDIR)/latex all-pdf-ja + @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." + +coverage: + $(SPHINXBUILD) -b coverage $(ALLSPHINXOPTS) $(BUILDDIR)/coverage + @echo "Testing of coverage in the sources finished, look at the " \ + "results in $(BUILDDIR)/coverage/python.txt." + +xml: + $(SPHINXBUILD) -b xml $(ALLSPHINXOPTS) $(BUILDDIR)/xml + @echo + @echo "Build finished. The XML files are in $(BUILDDIR)/xml." + +pseudoxml: + $(SPHINXBUILD) -b pseudoxml $(ALLSPHINXOPTS) $(BUILDDIR)/pseudoxml + @echo + @echo "Build finished. The pseudo-XML files are in $(BUILDDIR)/pseudoxml." diff --git a/doc/conf.py b/doc/conf.py new file mode 100644 index 0000000..cba11b5 --- /dev/null +++ b/doc/conf.py @@ -0,0 +1,290 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +# +# GeoMag documentation build configuration file, created by +# sphinx-quickstart on Fri Oct 16 16:45:07 2015. +# +# 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. + +import sys +import os +import shlex + +# 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('.')) +sys.path.append(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', +] + +# Add any paths that contain templates here, relative to this directory. +templates_path = ['_templates'] + +# The suffix(es) of source filenames. +# You can specify multiple suffix as a list of string: +# source_suffix = ['.rst', '.md'] +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 = 'GeoMag' +copyright = '2015, Christopher Weiss' +author = 'Christopher Weiss' + +# 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.0.1' +# The full version, including alpha/beta/rc tags. +release = '0.0.1' + +# The language for content autogenerated by Sphinx. Refer to documentation +# for a list of supported languages. +# +# This is also used if you do content translation via gettext catalogs. +# Usually you set "language" from the command line for these cases. +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 = ['_build'] + +# 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 = [] + +# If true, keep warnings as "system message" paragraphs in the built documents. +#keep_warnings = False + +# If true, `todo` and `todoList` produce output, else they produce nothing. +todo_include_todos = False + + +# -- 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 = 'alabaster' + +# 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 +# " v 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'] + +# Add any extra paths that contain custom files (such as robots.txt or +# .htaccess) here, relative to this directory. These files are copied +# directly to the root of the documentation. +#html_extra_path = [] + +# 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 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 + +# Language to be used for generating the HTML full-text search index. +# Sphinx supports the following languages: +# 'da', 'de', 'en', 'es', 'fi', 'fr', 'h', 'it', 'ja' +# 'nl', 'no', 'pt', 'ro', 'r', 'sv', 'tr' +#html_search_language = 'en' + +# A dictionary with options for the search language support, empty by default. +# Now only 'ja' uses this config value +#html_search_options = {'type': 'default'} + +# The name of a javascript file (relative to the configuration directory) that +# implements a search results scorer. If empty, the default will be used. +#html_search_scorer = 'scorer.js' + +# Output file base name for HTML help builder. +htmlhelp_basename = 'GeoMagdoc' + +# -- Options for LaTeX output --------------------------------------------- + +latex_elements = { +# The paper size ('letterpaper' or 'a4paper'). +#'papersize': 'letterpaper', + +# The font size ('10pt', '11pt' or '12pt'). +#'pointsize': '10pt', + +# Additional stuff for the LaTeX preamble. +#'preamble': '', + +# Latex figure (float) alignment +#'figure_align': 'htbp', +} + +# Grouping the document tree into LaTeX files. List of tuples +# (source start file, target name, title, +# author, documentclass [howto, manual, or own class]). +latex_documents = [ + (master_doc, 'GeoMag.tex', 'GeoMag Documentation', + 'Todd Dembrey', '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 + +# 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 = [ + (master_doc, 'geomag', 'GeoMag Documentation', + [author], 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 = [ + (master_doc, 'GeoMag', 'GeoMag Documentation', + author, 'GeoMag', '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' + +# If true, do not generate a @detailmenu in the "Top" node's menu. +#texinfo_no_detailmenu = False + +autoclass_content = 'both' diff --git a/doc/index.rst b/doc/index.rst new file mode 100644 index 0000000..55bf070 --- /dev/null +++ b/doc/index.rst @@ -0,0 +1,30 @@ +.. GeoMag documentation master file, created by + sphinx-quickstart on Fri Oct 16 16:45:07 2015. + You can adapt this file completely to your liking, but it should at least + contain the root `toctree` directive. + +Welcome to GeoMag's documentation! +================================== + +.. toctree:: + :maxdepth: 2 + + + +Indices and tables +================== + +* :ref:`genindex` +* :ref:`modindex` +* :ref:`search` + +Module Information +================== + +.. automodule:: geomag + +World Magnetic Model +================== + +.. autoclass:: geomag.WorldMagneticModel + :members: diff --git a/doc/make.bat b/doc/make.bat new file mode 100644 index 0000000..035556b --- /dev/null +++ b/doc/make.bat @@ -0,0 +1,263 @@ +@ECHO OFF + +REM Command file for Sphinx documentation + +if "%SPHINXBUILD%" == "" ( + set SPHINXBUILD=sphinx-build +) +set BUILDDIR=_build +set ALLSPHINXOPTS=-d %BUILDDIR%/doctrees %SPHINXOPTS% . +set I18NSPHINXOPTS=%SPHINXOPTS% . +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 ^` where ^ 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. xml to make Docutils-native XML files + echo. pseudoxml to make pseudoxml-XML files for display purposes + echo. linkcheck to check all external links for integrity + echo. doctest to run all doctests embedded in the documentation if enabled + echo. coverage to run coverage check of 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 +) + + +REM Check if sphinx-build is available and fallback to Python version if any +%SPHINXBUILD% 2> nul +if errorlevel 9009 goto sphinx_python +goto sphinx_ok + +:sphinx_python + +set SPHINXBUILD=python -m sphinx.__init__ +%SPHINXBUILD% 2> nul +if errorlevel 9009 ( + echo. + echo.The 'sphinx-build' command was not found. Make sure you have Sphinx + echo.installed, then set the SPHINXBUILD environment variable to point + echo.to the full path of the 'sphinx-build' executable. Alternatively you + echo.may add the Sphinx directory to PATH. + echo. + echo.If you don't have Sphinx installed, grab it from + echo.http://sphinx-doc.org/ + exit /b 1 +) + +:sphinx_ok + + +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\GeoMAg.qhcp + echo.To view the help file: + echo.^> assistant -collectionFile %BUILDDIR%\qthelp\GeoMAg.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" == "latexpdf" ( + %SPHINXBUILD% -b latex %ALLSPHINXOPTS% %BUILDDIR%/latex + cd %BUILDDIR%/latex + make all-pdf + cd %~dp0 + echo. + echo.Build finished; the PDF files are in %BUILDDIR%/latex. + goto end +) + +if "%1" == "latexpdfja" ( + %SPHINXBUILD% -b latex %ALLSPHINXOPTS% %BUILDDIR%/latex + cd %BUILDDIR%/latex + make all-pdf-ja + cd %~dp0 + echo. + echo.Build finished; the PDF 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 +) + +if "%1" == "coverage" ( + %SPHINXBUILD% -b coverage %ALLSPHINXOPTS% %BUILDDIR%/coverage + if errorlevel 1 exit /b 1 + echo. + echo.Testing of coverage in the sources finished, look at the ^ +results in %BUILDDIR%/coverage/python.txt. + goto end +) + +if "%1" == "xml" ( + %SPHINXBUILD% -b xml %ALLSPHINXOPTS% %BUILDDIR%/xml + if errorlevel 1 exit /b 1 + echo. + echo.Build finished. The XML files are in %BUILDDIR%/xml. + goto end +) + +if "%1" == "pseudoxml" ( + %SPHINXBUILD% -b pseudoxml %ALLSPHINXOPTS% %BUILDDIR%/pseudoxml + if errorlevel 1 exit /b 1 + echo. + echo.Build finished. The pseudo-XML files are in %BUILDDIR%/pseudoxml. + goto end +) + +:end diff --git a/examples/world_plot.py b/examples/world_plot.py new file mode 100644 index 0000000..a700c9a --- /dev/null +++ b/examples/world_plot.py @@ -0,0 +1,42 @@ +''' +Plots declination for +-70deg latitude and 260deg longditude +''' +import matplotlib.pyplot as plt +import matplotlib.cm as cm +from mpl_toolkits.basemap import Basemap +from datetime import date +import numpy as np + +from geomag import WorldMagneticModel + + +def _gen_2d_array(size_x, size_y, default=None): + return [[default] * size_x for _ in range(size_y)] + +wmm = WorldMagneticModel() +lat = range(-70, 70) +lon = range(360) +len_lat = len(lat) +len_lon = len(lon) +world = _gen_2d_array(len_lon, len_lat, date(2015, 1, 1)) +x = _gen_2d_array(len_lon, len_lat, date(2015, 1, 1)) +y = _gen_2d_array(len_lon, len_lat, date(2015, 1, 1)) +for i, _lon in enumerate(lon): + for j, _lat in enumerate(lat): + world[j][i] = wmm.calc_mag_field(_lat, _lon).declination + x[j][i] = _lon + y[j][i] = _lat + +plot_levels = np.arange(-90, 90, 5) + +map = Basemap(llcrnrlon=lon[0], llcrnrlat=lat[0], urcrnrlon=lon[-1], urcrnrlat=lat[-1], projection='mill') +map_x, map_y = map(np.asarray(x), np.asarray(y)) +contour_set = map.contour(map_x, map_y, world, plot_levels, linewidths=3, cmap=cm.jet) +map.drawcoastlines(linewidth=1.25) +map.fillcontinents(color='0.8') +map.drawparallels(np.arange(lat[0],lat[-1],10),labels=[1,1,0,0]) +map.drawmeridians(np.arange(lon[0],lon[-1],20),labels=[0,0,0,1]) + +plt.clabel(contour_set, fontsize='xx-small', inline_spacing=1.5, fmt='%1.0f') +plt.suptitle('Magnetic Declination (degrees)') +plt.show() diff --git a/geomag/__init__.py b/geomag/__init__.py new file mode 100644 index 0000000..542a92a --- /dev/null +++ b/geomag/__init__.py @@ -0,0 +1,24 @@ +""" +by Christopher Weiss cmweiss@gmail.com +and Todd Dembrey todd.dembrey@gmail.com + +Adapted from the geomagc software and World Magnetic Model of the NOAA +Satellite and Information Service, National Geophysical Data Center +http://www.ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml + +Suggestions for improvements are appreciated. + +""" +from .world_magnetic_model import WorldMagneticModel + +__all__ = ['WorldMagneticModel'] + + +# The following functions are for backwards compatability + +def declination(*args, **kargs): + return WorldMagneticModel().calc_mag_field(*args, **kargs).declination + + +def heading(hdg, *args, **kargs): + return WorldMagneticModel().calc_mag_field(*args, **kargs).mag_heading(hdg) diff --git a/geomag/geomag/__init__.py b/geomag/geomag/__init__.py deleted file mode 100644 index 0f65a81..0000000 --- a/geomag/geomag/__init__.py +++ /dev/null @@ -1,36 +0,0 @@ -"""geomag package -by Christopher Weiss cmweiss@gmail.com - -Adapted from the geomagc software and World Magnetic Model of the NOAA -Satellite and Information Service, National Geophysical Data Center -http://www.ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml - -Suggestions for improvements are appreciated. - -USAGE: ->>> import geomag ->>> geomag.declination(80,0) --6.1335150785195536 -""" - -from . import geomag - -__singleton__ = geomag.GeoMag() - -def declination(*args, **kargs): - """Calculate magnetic declination in degrees - dlat = latitude in degrees - dlon = longitude in degrees - h = altitude in feet, default=0 - calc_date = date for computing declination, default=today - """ - mag = __singleton__.GeoMag(*args, **kargs) - return mag.dec - -def mag_heading(hdg, *args, **kargs): - """Calculates the magnetic heading from a true heading. - hdg = true heading in degrees - All other parameters are the same as declination. - """ - dec = declination(*args, **kargs) - return (hdg - dec + 360.0) % 360 diff --git a/geomag/geomag/geomag.py b/geomag/geomag/geomag.py deleted file mode 100644 index d546d55..0000000 --- a/geomag/geomag/geomag.py +++ /dev/null @@ -1,309 +0,0 @@ -# geomag.py -# by Christopher Weiss cmweiss@gmail.com - -# Adapted from the geomagc software and World Magnetic Model of the NOAA -# Satellite and Information Service, National Geophysical Data Center -# http://www.ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml -# -# Suggestions for improvements are appreciated. - -# USAGE: -# -# >>> gm = geomag.GeoMag("WMM.COF") -# >>> mag = gm.GeoMag(80,0) -# >>> mag.dec -# -6.1335150785195536 -# >>> - -import math, os, unittest -from datetime import date - -class GeoMag: - - def GeoMag(self, dlat, dlon, h=0, time=date.today()): # latitude (decimal degrees), longitude (decimal degrees), altitude (feet), date - #time = date('Y') + date('z')/365 - time = time.year+((time - date(time.year,1,1)).days/365.0) - alt = h/3280.8399 - - otime = oalt = olat = olon = -1000.0 - - dt = time - self.epoch - glat = dlat - glon = dlon - rlat = math.radians(glat) - rlon = math.radians(glon) - srlon = math.sin(rlon) - srlat = math.sin(rlat) - crlon = math.cos(rlon) - crlat = math.cos(rlat) - srlat2 = srlat*srlat - crlat2 = crlat*crlat - self.sp[1] = srlon - self.cp[1] = crlon - - #/* CONVERT FROM GEODETIC COORDS. TO SPHERICAL COORDS. */ - if (alt != oalt or glat != olat): - q = math.sqrt(self.a2-self.c2*srlat2) - q1 = alt*q - q2 = ((q1+self.a2)/(q1+self.b2))*((q1+self.a2)/(q1+self.b2)) - ct = srlat/math.sqrt(q2*crlat2+srlat2) - st = math.sqrt(1.0-(ct*ct)) - r2 = (alt*alt)+2.0*q1+(self.a4-self.c4*srlat2)/(q*q) - r = math.sqrt(r2) - d = math.sqrt(self.a2*crlat2+self.b2*srlat2) - ca = (alt+d)/r - sa = self.c2*crlat*srlat/(r*d) - - if (glon != olon): - for m in range(2,self.maxord+1): - self.sp[m] = self.sp[1]*self.cp[m-1]+self.cp[1]*self.sp[m-1] - self.cp[m] = self.cp[1]*self.cp[m-1]-self.sp[1]*self.sp[m-1] - - aor = self.re/r - ar = aor*aor - br = bt = bp = bpp = 0.0 - for n in range(1,self.maxord+1): - ar = ar*aor - - #for (m=0,D3=1,D4=(n+m+D3)/D3;D4>0;D4--,m+=D3): - m=0 - D3=1 - #D4=(n+m+D3)/D3 - D4=(n+m+1) - while D4>0: - - # /* - # COMPUTE UNNORMALIZED ASSOCIATED LEGENDRE POLYNOMIALS - # AND DERIVATIVES VIA RECURSION RELATIONS - # */ - if (alt != oalt or glat != olat): - if (n == m): - self.p[m][n] = st * self.p[m-1][n-1] - self.dp[m][n] = st*self.dp[m-1][n-1]+ct*self.p[m-1][n-1] - - elif (n == 1 and m == 0): - self.p[m][n] = ct*self.p[m][n-1] - self.dp[m][n] = ct*self.dp[m][n-1]-st*self.p[m][n-1] - - elif (n > 1 and n != m): - if (m > n-2): - self.p[m][n-2] = 0 - if (m > n-2): - self.dp[m][n-2] = 0.0 - self.p[m][n] = ct*self.p[m][n-1]-self.k[m][n]*self.p[m][n-2] - self.dp[m][n] = ct*self.dp[m][n-1] - st*self.p[m][n-1]-self.k[m][n]*self.dp[m][n-2] - - # /* - # TIME ADJUST THE GAUSS COEFFICIENTS - # */ - if (time != otime): - self.tc[m][n] = self.c[m][n]+dt*self.cd[m][n] - if (m != 0): - self.tc[n][m-1] = self.c[n][m-1]+dt*self.cd[n][m-1] - - # /* - # ACCUMULATE TERMS OF THE SPHERICAL HARMONIC EXPANSIONS - # */ - par = ar*self.p[m][n] - - if (m == 0): - temp1 = self.tc[m][n]*self.cp[m] - temp2 = self.tc[m][n]*self.sp[m] - else: - temp1 = self.tc[m][n]*self.cp[m]+self.tc[n][m-1]*self.sp[m] - temp2 = self.tc[m][n]*self.sp[m]-self.tc[n][m-1]*self.cp[m] - - bt = bt-ar*temp1*self.dp[m][n] - bp = bp + (self.fm[m] * temp2 * par) - br = br + (self.fn[n] * temp1 * par) - # /* - # SPECIAL CASE: NORTH/SOUTH GEOGRAPHIC POLES - # */ - if (st == 0.0 and m == 1): - if (n == 1): - self.pp[n] = self.pp[n-1] - else: - self.pp[n] = ct*self.pp[n-1]-self.k[m][n]*self.pp[n-2] - parp = ar*self.pp[n] - bpp = bpp + (self.fm[m]*temp2*parp) - - D4=D4-1 - m=m+1 - - if (st == 0.0): - bp = bpp - else: - bp = bp/st - # /* - # ROTATE MAGNETIC VECTOR COMPONENTS FROM SPHERICAL TO - # GEODETIC COORDINATES - # */ - bx = -bt*ca-br*sa - by = bp - bz = bt*sa-br*ca - # /* - # COMPUTE DECLINATION (DEC), INCLINATION (DIP) AND - # TOTAL INTENSITY (TI) - # */ - bh = math.sqrt((bx*bx)+(by*by)) - ti = math.sqrt((bh*bh)+(bz*bz)) - dec = math.degrees(math.atan2(by,bx)) - dip = math.degrees(math.atan2(bz,bh)) - # /* - # COMPUTE MAGNETIC GRID VARIATION IF THE CURRENT - # GEODETIC POSITION IS IN THE ARCTIC OR ANTARCTIC - # (I.E. GLAT > +55 DEGREES OR GLAT < -55 DEGREES) - - # OTHERWISE, SET MAGNETIC GRID VARIATION TO -999.0 - # */ - gv = -999.0 - if (math.fabs(glat) >= 55.): - if (glat > 0.0 and glon >= 0.0): - gv = dec-glon - if (glat > 0.0 and glon < 0.0): - gv = dec+math.fabs(glon); - if (glat < 0.0 and glon >= 0.0): - gv = dec+glon - if (glat < 0.0 and glon < 0.0): - gv = dec-math.fabs(glon) - if (gv > +180.0): - gv = gv - 360.0 - if (gv < -180.0): - gv = gv + 360.0 - - otime = time - oalt = alt - olat = glat - olon = glon - - class RetObj: - pass - retobj = RetObj() - retobj.dec = dec - retobj.dip = dip - retobj.ti = ti - retobj.bh = bh - retobj.bx = bx - retobj.by = by - retobj.bz = bz - retobj.lat = dlat - retobj.lon = dlon - retobj.alt = h - retobj.time = time - - return retobj - - def __init__(self, wmm_filename=None): - if not wmm_filename: - wmm_filename = os.path.join(os.path.dirname(__file__), 'WMM.COF') - wmm=[] - with open(wmm_filename) as wmm_file: - for line in wmm_file: - linevals = line.strip().split() - if len(linevals) == 3: - self.epoch = float(linevals[0]) - self.model = linevals[1] - self.modeldate = linevals[2] - elif len(linevals) == 6: - linedict = {'n': int(float(linevals[0])), - 'm': int(float(linevals[1])), - 'gnm': float(linevals[2]), - 'hnm': float(linevals[3]), - 'dgnm': float(linevals[4]), - 'dhnm': float(linevals[5])} - wmm.append(linedict) - - z = [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] - self.maxord = self.maxdeg = 12 - self.tc = [z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13]] - self.sp = z[0:14] - self.cp = z[0:14] - self.cp[0] = 1.0 - self.pp = z[0:13] - self.pp[0] = 1.0 - self.p = [z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14]] - self.p[0][0] = 1.0 - self.dp = [z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13]] - self.a = 6378.137 - self.b = 6356.7523142 - self.re = 6371.2 - self.a2 = self.a*self.a - self.b2 = self.b*self.b - self.c2 = self.a2-self.b2 - self.a4 = self.a2*self.a2 - self.b4 = self.b2*self.b2 - self.c4 = self.a4 - self.b4 - - self.c = [z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14]] - self.cd = [z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14]] - - for wmmnm in wmm: - m = wmmnm['m'] - n = wmmnm['n'] - gnm = wmmnm['gnm'] - hnm = wmmnm['hnm'] - dgnm = wmmnm['dgnm'] - dhnm = wmmnm['dhnm'] - if (m <= n): - self.c[m][n] = gnm - self.cd[m][n] = dgnm - if (m != 0): - self.c[n][m-1] = hnm - self.cd[n][m-1] = dhnm - - #/* CONVERT SCHMIDT NORMALIZED GAUSS COEFFICIENTS TO UNNORMALIZED */ - self.snorm = [z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13]] - self.snorm[0][0] = 1.0 - self.k = [z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13]] - self.k[1][1] = 0.0 - self.fn = [0.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0,11.0,12.0,13.0] - self.fm = [0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0,11.0,12.0] - for n in range(1,self.maxord+1): - self.snorm[0][n] = self.snorm[0][n-1]*(2.0*n-1)/n - j=2.0 - #for (m=0,D1=1,D2=(n-m+D1)/D1;D2>0;D2--,m+=D1): - m=0 - D1=1 - D2=(n-m+D1)/D1 - while (D2 > 0): - self.k[m][n] = (((n-1)*(n-1))-(m*m))/((2.0*n-1)*(2.0*n-3.0)) - if (m > 0): - flnmj = ((n-m+1.0)*j)/(n+m) - self.snorm[m][n] = self.snorm[m-1][n]*math.sqrt(flnmj) - j = 1.0 - self.c[n][m-1] = self.snorm[m][n]*self.c[n][m-1] - self.cd[n][m-1] = self.snorm[m][n]*self.cd[n][m-1] - self.c[m][n] = self.snorm[m][n]*self.c[m][n] - self.cd[m][n] = self.snorm[m][n]*self.cd[m][n] - D2=D2-1 - m=m+D1 - -class GeoMagTest(unittest.TestCase): - - d1=date(2015,1,1) - d2=date(2017,7,2) - - test_values = ( - # date, alt, lat, lon, var - (d1, 0, 80, 0, -3.85), - (d1, 0, 0, 120, 0.57), - (d1, 0, -80, 240, 69.81), - (d1, 328083.99, 80, 0, -4.27), - (d1, 328083.99, 0, 120, 0.56), - (d1, 328083.99, -80, 240, 69.22), - (d2, 0, 80, 0, -2.75), - (d2, 0, 0, 120, 0.32), - (d2, 0, -80, 240, 69.58), - (d2, 328083.99, 80, 0, -3.17), - (d2, 328083.99, 0, 120, 0.32), - (d2, 328083.99, -80, 240, 69.00), - ) - - def test_declination(self): - gm = GeoMag() - for values in self.test_values: - calcval=gm.GeoMag(values[2], values[3], values[1], values[0]) - self.assertAlmostEqual(values[4], calcval.dec, 2, 'Expected %s, result %s' % (values[4], calcval.dec)) - -if __name__ == '__main__': - unittest.main() diff --git a/geomag/latlon.py b/geomag/latlon.py new file mode 100644 index 0000000..a3469b0 --- /dev/null +++ b/geomag/latlon.py @@ -0,0 +1,106 @@ +from __future__ import division + +import math + + +def normalise_plus_minus_range(value, norm_range): + """Normalise a value to within a positive and negative range + + **Parameters** + + value + input value to be limited + norm_range : {'lat','lon',*number*} + A number will limit the range between +/- that value. 'lat' + and 'lon' will limit within +/-90 and 180 respectively. + + """ + norm_range_dict = {'lat': 90, 'lon': 180} + try: + valid_range = norm_range_dict[norm_range] + value %= [-1, 1][value > 0] * valid_range * 2 + except KeyError: + valid_range = norm_range + if abs(value) > valid_range: + return value % ([1, -1][value > 0] * valid_range) + else: + return value + + +class LatLon(object): + """Implements a latitude Longitude class with conversion to spherical co-ords""" + + def __init__(self, latitude, longitude, unit='degrees'): + """Performs conversion and stored both degrees and radians. + The range is limited to within +-90 and +-180 degrees + """ + + if unit in ['deg', 'degrees']: + _latitude_deg = normalise_plus_minus_range(latitude, 'lat') + _longitude_deg = normalise_plus_minus_range(longitude, 'lon') + _latitude_rad = math.radians(_latitude_deg) + _longitude_rad = math.radians(_longitude_deg) + elif unit in ['rad', 'radians']: + _latitude_deg = normalise_plus_minus_range( + math.degrees(latitude), 'lat') + _longitude_deg = normalise_plus_minus_range( + math.degrees(longitude), 'lon') + _latitude_rad = math.radians(_latitude_deg) + _longitude_rad = math.radians(_longitude_deg) + else: + raise Exception('Incorrect unit specified, please use "deg", "degrees", "rad" or "radians"') + self._lat = _latitude_deg + self._lon = _longitude_deg + self._lat_rad = _latitude_rad + self._lon_rad = _longitude_rad + + self._default_unit = unit + + @property + def lat(self): + """Return Latitude""" + return self._lat + + @property + def lon(self): + """Return Longitude""" + return self._lon + + @property + def lon_rad(self): + """Return Longitude""" + return self._lon_rad + + @property + def lat_rad(self): + """Return Longitude""" + return self._lat_rad + + def __call__(self): + return self._lat_lon() + + def __str__(self): + return "({} N,{} W)".format(*self._lat_lon()) + + def _lat_lon(self): + return self.lat, self.lon + + def convert_spherical(self, altitude, unit='radians'): + """Converts to spherical co-ordinates + + Returns in radians as standard + """ + cur_lat_in_radians = self._lat_rad + equator_radius = 6378137 + flattening = 1 / 298.257223563 + eccentricity_squared = (flattening * (2 - flattening)) + prime_vertical = equator_radius / math.sqrt(1 - eccentricity_squared * math.sin(cur_lat_in_radians)**2) + p = (prime_vertical + altitude) * math.cos(cur_lat_in_radians) + z = (prime_vertical * (1 - eccentricity_squared) + altitude) * math.sin(cur_lat_in_radians) + r = math.sqrt(p**2 + z**2) + spherical_latitude = math.asin(z / r) + spherical_longitude = self.lon + if unit == 'degrees': + spherical_latitude = math.degrees(spherical_latitude) + spherical_longitude = math.degrees(spherical_longitude) + return spherical_latitude, spherical_longitude, r diff --git a/geomag/geomag/WMM.COF b/geomag/model_data/WMM.COF similarity index 100% rename from geomag/geomag/WMM.COF rename to geomag/model_data/WMM.COF diff --git a/geomag/geomag/WMM2010.COF b/geomag/model_data/WMM2010.COF similarity index 98% rename from geomag/geomag/WMM2010.COF rename to geomag/model_data/WMM2010.COF index eb78503..6bfacd6 100644 --- a/geomag/geomag/WMM2010.COF +++ b/geomag/model_data/WMM2010.COF @@ -1,93 +1,93 @@ - 2010.0 WMM-2010 11/20/2009 - 1 0 -29496.6 0.0 11.6 0.0 - 1 1 -1586.3 4944.4 16.5 -25.9 - 2 0 -2396.6 0.0 -12.1 0.0 - 2 1 3026.1 -2707.7 -4.4 -22.5 - 2 2 1668.6 -576.1 1.9 -11.8 - 3 0 1340.1 0.0 0.4 0.0 - 3 1 -2326.2 -160.2 -4.1 7.3 - 3 2 1231.9 251.9 -2.9 -3.9 - 3 3 634.0 -536.6 -7.7 -2.6 - 4 0 912.6 0.0 -1.8 0.0 - 4 1 808.9 286.4 2.3 1.1 - 4 2 166.7 -211.2 -8.7 2.7 - 4 3 -357.1 164.3 4.6 3.9 - 4 4 89.4 -309.1 -2.1 -0.8 - 5 0 -230.9 0.0 -1.0 0.0 - 5 1 357.2 44.6 0.6 0.4 - 5 2 200.3 188.9 -1.8 1.8 - 5 3 -141.1 -118.2 -1.0 1.2 - 5 4 -163.0 0.0 0.9 4.0 - 5 5 -7.8 100.9 1.0 -0.6 - 6 0 72.8 0.0 -0.2 0.0 - 6 1 68.6 -20.8 -0.2 -0.2 - 6 2 76.0 44.1 -0.1 -2.1 - 6 3 -141.4 61.5 2.0 -0.4 - 6 4 -22.8 -66.3 -1.7 -0.6 - 6 5 13.2 3.1 -0.3 0.5 - 6 6 -77.9 55.0 1.7 0.9 - 7 0 80.5 0.0 0.1 0.0 - 7 1 -75.1 -57.9 -0.1 0.7 - 7 2 -4.7 -21.1 -0.6 0.3 - 7 3 45.3 6.5 1.3 -0.1 - 7 4 13.9 24.9 0.4 -0.1 - 7 5 10.4 7.0 0.3 -0.8 - 7 6 1.7 -27.7 -0.7 -0.3 - 7 7 4.9 -3.3 0.6 0.3 - 8 0 24.4 0.0 -0.1 0.0 - 8 1 8.1 11.0 0.1 -0.1 - 8 2 -14.5 -20.0 -0.6 0.2 - 8 3 -5.6 11.9 0.2 0.4 - 8 4 -19.3 -17.4 -0.2 0.4 - 8 5 11.5 16.7 0.3 0.1 - 8 6 10.9 7.0 0.3 -0.1 - 8 7 -14.1 -10.8 -0.6 0.4 - 8 8 -3.7 1.7 0.2 0.3 - 9 0 5.4 0.0 -0.0 0.0 - 9 1 9.4 -20.5 -0.1 -0.0 - 9 2 3.4 11.5 0.0 -0.2 - 9 3 -5.2 12.8 0.3 0.0 - 9 4 3.1 -7.2 -0.4 -0.1 - 9 5 -12.4 -7.4 -0.3 0.1 - 9 6 -0.7 8.0 0.1 -0.0 - 9 7 8.4 2.1 -0.1 -0.2 - 9 8 -8.5 -6.1 -0.4 0.3 - 9 9 -10.1 7.0 -0.2 0.2 - 10 0 -2.0 0.0 0.0 0.0 - 10 1 -6.3 2.8 -0.0 0.1 - 10 2 0.9 -0.1 -0.1 -0.1 - 10 3 -1.1 4.7 0.2 0.0 - 10 4 -0.2 4.4 -0.0 -0.1 - 10 5 2.5 -7.2 -0.1 -0.1 - 10 6 -0.3 -1.0 -0.2 -0.0 - 10 7 2.2 -3.9 0.0 -0.1 - 10 8 3.1 -2.0 -0.1 -0.2 - 10 9 -1.0 -2.0 -0.2 0.0 - 10 10 -2.8 -8.3 -0.2 -0.1 - 11 0 3.0 0.0 0.0 0.0 - 11 1 -1.5 0.2 0.0 -0.0 - 11 2 -2.1 1.7 -0.0 0.1 - 11 3 1.7 -0.6 0.1 0.0 - 11 4 -0.5 -1.8 -0.0 0.1 - 11 5 0.5 0.9 0.0 0.0 - 11 6 -0.8 -0.4 -0.0 0.1 - 11 7 0.4 -2.5 -0.0 0.0 - 11 8 1.8 -1.3 -0.0 -0.1 - 11 9 0.1 -2.1 0.0 -0.1 - 11 10 0.7 -1.9 -0.1 -0.0 - 11 11 3.8 -1.8 -0.0 -0.1 - 12 0 -2.2 0.0 -0.0 0.0 - 12 1 -0.2 -0.9 0.0 -0.0 - 12 2 0.3 0.3 0.1 0.0 - 12 3 1.0 2.1 0.1 -0.0 - 12 4 -0.6 -2.5 -0.1 0.0 - 12 5 0.9 0.5 -0.0 -0.0 - 12 6 -0.1 0.6 0.0 0.1 - 12 7 0.5 -0.0 0.0 0.0 - 12 8 -0.4 0.1 -0.0 0.0 - 12 9 -0.4 0.3 0.0 -0.0 - 12 10 0.2 -0.9 0.0 -0.0 - 12 11 -0.8 -0.2 -0.1 0.0 - 12 12 0.0 0.9 0.1 0.0 -999999999999999999999999999999999999999999999999 -999999999999999999999999999999999999999999999999 + 2010.0 WMM-2010 11/20/2009 + 1 0 -29496.6 0.0 11.6 0.0 + 1 1 -1586.3 4944.4 16.5 -25.9 + 2 0 -2396.6 0.0 -12.1 0.0 + 2 1 3026.1 -2707.7 -4.4 -22.5 + 2 2 1668.6 -576.1 1.9 -11.8 + 3 0 1340.1 0.0 0.4 0.0 + 3 1 -2326.2 -160.2 -4.1 7.3 + 3 2 1231.9 251.9 -2.9 -3.9 + 3 3 634.0 -536.6 -7.7 -2.6 + 4 0 912.6 0.0 -1.8 0.0 + 4 1 808.9 286.4 2.3 1.1 + 4 2 166.7 -211.2 -8.7 2.7 + 4 3 -357.1 164.3 4.6 3.9 + 4 4 89.4 -309.1 -2.1 -0.8 + 5 0 -230.9 0.0 -1.0 0.0 + 5 1 357.2 44.6 0.6 0.4 + 5 2 200.3 188.9 -1.8 1.8 + 5 3 -141.1 -118.2 -1.0 1.2 + 5 4 -163.0 0.0 0.9 4.0 + 5 5 -7.8 100.9 1.0 -0.6 + 6 0 72.8 0.0 -0.2 0.0 + 6 1 68.6 -20.8 -0.2 -0.2 + 6 2 76.0 44.1 -0.1 -2.1 + 6 3 -141.4 61.5 2.0 -0.4 + 6 4 -22.8 -66.3 -1.7 -0.6 + 6 5 13.2 3.1 -0.3 0.5 + 6 6 -77.9 55.0 1.7 0.9 + 7 0 80.5 0.0 0.1 0.0 + 7 1 -75.1 -57.9 -0.1 0.7 + 7 2 -4.7 -21.1 -0.6 0.3 + 7 3 45.3 6.5 1.3 -0.1 + 7 4 13.9 24.9 0.4 -0.1 + 7 5 10.4 7.0 0.3 -0.8 + 7 6 1.7 -27.7 -0.7 -0.3 + 7 7 4.9 -3.3 0.6 0.3 + 8 0 24.4 0.0 -0.1 0.0 + 8 1 8.1 11.0 0.1 -0.1 + 8 2 -14.5 -20.0 -0.6 0.2 + 8 3 -5.6 11.9 0.2 0.4 + 8 4 -19.3 -17.4 -0.2 0.4 + 8 5 11.5 16.7 0.3 0.1 + 8 6 10.9 7.0 0.3 -0.1 + 8 7 -14.1 -10.8 -0.6 0.4 + 8 8 -3.7 1.7 0.2 0.3 + 9 0 5.4 0.0 -0.0 0.0 + 9 1 9.4 -20.5 -0.1 -0.0 + 9 2 3.4 11.5 0.0 -0.2 + 9 3 -5.2 12.8 0.3 0.0 + 9 4 3.1 -7.2 -0.4 -0.1 + 9 5 -12.4 -7.4 -0.3 0.1 + 9 6 -0.7 8.0 0.1 -0.0 + 9 7 8.4 2.1 -0.1 -0.2 + 9 8 -8.5 -6.1 -0.4 0.3 + 9 9 -10.1 7.0 -0.2 0.2 + 10 0 -2.0 0.0 0.0 0.0 + 10 1 -6.3 2.8 -0.0 0.1 + 10 2 0.9 -0.1 -0.1 -0.1 + 10 3 -1.1 4.7 0.2 0.0 + 10 4 -0.2 4.4 -0.0 -0.1 + 10 5 2.5 -7.2 -0.1 -0.1 + 10 6 -0.3 -1.0 -0.2 -0.0 + 10 7 2.2 -3.9 0.0 -0.1 + 10 8 3.1 -2.0 -0.1 -0.2 + 10 9 -1.0 -2.0 -0.2 0.0 + 10 10 -2.8 -8.3 -0.2 -0.1 + 11 0 3.0 0.0 0.0 0.0 + 11 1 -1.5 0.2 0.0 -0.0 + 11 2 -2.1 1.7 -0.0 0.1 + 11 3 1.7 -0.6 0.1 0.0 + 11 4 -0.5 -1.8 -0.0 0.1 + 11 5 0.5 0.9 0.0 0.0 + 11 6 -0.8 -0.4 -0.0 0.1 + 11 7 0.4 -2.5 -0.0 0.0 + 11 8 1.8 -1.3 -0.0 -0.1 + 11 9 0.1 -2.1 0.0 -0.1 + 11 10 0.7 -1.9 -0.1 -0.0 + 11 11 3.8 -1.8 -0.0 -0.1 + 12 0 -2.2 0.0 -0.0 0.0 + 12 1 -0.2 -0.9 0.0 -0.0 + 12 2 0.3 0.3 0.1 0.0 + 12 3 1.0 2.1 0.1 -0.0 + 12 4 -0.6 -2.5 -0.1 0.0 + 12 5 0.9 0.5 -0.0 -0.0 + 12 6 -0.1 0.6 0.0 0.1 + 12 7 0.5 -0.0 0.0 0.0 + 12 8 -0.4 0.1 -0.0 0.0 + 12 9 -0.4 0.3 0.0 -0.0 + 12 10 0.2 -0.9 0.0 -0.0 + 12 11 -0.8 -0.2 -0.1 0.0 + 12 12 0.0 0.9 0.1 0.0 +999999999999999999999999999999999999999999999999 +999999999999999999999999999999999999999999999999 diff --git a/geomag/scalar_potential.py b/geomag/scalar_potential.py new file mode 100644 index 0000000..065d4b3 --- /dev/null +++ b/geomag/scalar_potential.py @@ -0,0 +1,180 @@ +from __future__ import division + +import math +try: + import functools32 as functools +except ImportError: + import functools +import operator + + +def _gen_2d_array(size_x, size_y, default=None): + return [[default] * size_x for _ in range(size_y)] + + +def kronecker_delta(j, i): + """ Kronecker delta is defined as Iji = 1 if i = j and Iji = 0 otherwise""" + return 1 if j == i else 0 + + +def double_factorial(n): + """Returns the double factorial (n!!) + + Input must be an integer or greater than 0 otherwise ValueError is raised + + """ + if int(n) != n: + raise ValueError('n must be an integer') + elif n < 0: + raise ValueError('n must be greater than 0') + start_number = 2 - (n % 2) + return functools.reduce(operator.mul, range(start_number, n + 1, 2)) + + +def schmidt_quasi_normalisation(max_n): + """Returns an array of the Schmidt Quasi-normalised values + + Array is symmetrical about the diagonal + """ + schmidt = _gen_2d_array(max_n, max_n, 0.0) + for n in range(max_n): + for m in range(n + 1): + if n == 0: + # This is a bit of a hack to get round 2n-1 evaluating to + # -1 and erroring when it should be returning 1 + double_fact = 1.0 + else: + double_fact = double_factorial(2 * n - 1) + k_delta = kronecker_delta(0, m) + n_minus_m_fact = math.factorial(n - m) + n_plus_m_fact = math.factorial(n + m) + schmidt[m][n] = math.sqrt( + ((2 - k_delta) * n_minus_m_fact) / n_plus_m_fact) * double_fact / n_minus_m_fact + return schmidt + + +def tuplize_array(array): + return tuple(tuple(line) for line in array) + + +def recursion_constants(max_n): + """ Calculates the values useful for performing the scalar potential algorithm """ + k = _gen_2d_array(max_n, max_n, 0) + for n in range(1, max_n): + for m in range(n + 1): + k[m][n] = (((n - 1) * (n - 1)) - (m**2)) / ((2.0 * n - 1) * (2.0 * n - 3.0)) + return tuplize_array(k) + + +@functools.lru_cache(maxsize=None) +def associated_polynomials(cos_theta, sin_theta, max_n, max_m, k=None): + """Specific legendre legrende_polynomials for application in geomagnetic calculation + + Recursion Constants (k) can be pre-calculated to improve performance, otherwise the + function will perform this for you + + """ + if k is None: + k = recursion_constants(max_n) + associated_poly = _gen_2d_array(max_n, max_m, 0.0) + derivative_associated_poly = _gen_2d_array(max_n, max_m, 0.0) + associated_poly[0][0] = 1 + for n in range(1, max_n): + for m in range(n + 1): + if m == n: + previous_ass_poly = associated_poly[m - 1][n - 1] + associated_poly[m][n] = cos_theta * previous_ass_poly + derivative_associated_poly[m][n] = (cos_theta * derivative_associated_poly[m - 1][n - 1] + + sin_theta * previous_ass_poly) + else: + previous_ass_poly = associated_poly[m][n - 1] + current_k = k[m][n] + associated_poly[m][n] = sin_theta * previous_ass_poly - current_k * associated_poly[m][n - 2] + derivative_associated_poly[m][n] = (sin_theta * derivative_associated_poly[m][n - 1] + - cos_theta * previous_ass_poly + - current_k * derivative_associated_poly[m][n - 2]) + + return associated_poly, derivative_associated_poly + + +def legrende_polynomials(max_n, v): + """ Returns the legendre polynomials + Based on the algorithm here: + http://uk.mathworks.com/help/symbolic/mupad_ref/orthpoly-legendre.html + """ + poly = [1, v] + [0] * (max_n - 1) + for n in range(2, max_n + 1): + poly[n] = (2 * n - 1) / n * v * poly[n - 1] - (n - 1) / n * poly[n - 2] + return poly + + +def multiple_angle(phi, max_n): + """ Returns the values of sin(nx) and cos(nx) series + + """ + sin_phi = math.sin(phi) + cos_phi = math.cos(phi) + sin_n = [0] * max_n + cos_n = [1] * max_n + sin_n[1] = sin_phi + cos_n[1] = cos_phi + for i in range(2, max_n): + last_sin = sin_n[i - 1] + last_cos = cos_n[i - 1] + sin_n[i] = sin_phi * last_cos + cos_phi * last_sin + cos_n[i] = cos_phi * last_cos - sin_phi * last_sin + return tuple(sin_n), tuple(cos_n) + + +def scalar_potential(coeff, phi, theta, max_n, max_m, radial_alt, ref_radius=6371200, k=None): + """ + + ref_radius defaults to the radius of the earth in m. + """ + cos_theta = math.cos(theta) + if cos_theta is 0: + cos_theta += 0.00000001 + sin_theta = math.sin(theta) + sin_n, cos_n = multiple_angle(phi, max_n) + legendre_poly, legendre_poly_derivative = associated_polynomials(cos_theta, sin_theta, max_n, max_m, k) + a_over_r_pow = _altitude_ratios(ref_radius, radial_alt, max_n) + magnetic_field = _calc_scalar_potential( + coeff, + cos_theta, + sin_n, + cos_n, + legendre_poly, + legendre_poly_derivative, + a_over_r_pow, + max_n) + return magnetic_field + + +def _calc_scalar_potential(coeff, cos_theta, sin_n, cos_n, legendre_poly, legendre_poly_der, a_over_r_pow, max_n): + """ Calculates the partial scalar potential values""" + B_r = 0 + B_theta = 0 + B_phi = 0 + for n in range(1, max_n): + current_a_over_r_power = a_over_r_pow[n] + for m in range(n + 1): + current_cos_n = cos_n[m] + current_sin_n = sin_n[m] + g_cos = coeff[m][n] * current_cos_n + g_sin = coeff[m][n] * current_sin_n + h_sin = coeff[n][m - 1] * current_sin_n + h_cos = coeff[n][m - 1] * current_cos_n + B_r += current_a_over_r_power * (n + 1) * (g_cos + h_sin) * legendre_poly[m][n] + B_theta -= current_a_over_r_power * (g_cos + h_sin) * legendre_poly_der[m][n] + B_phi -= current_a_over_r_power * m * (-g_sin + h_cos) * legendre_poly[m][n] + try: + B_phi *= 1 / cos_theta + except ZeroDivisionError: + B_phi = B_phi + return B_r, B_theta, B_phi + + +@functools.lru_cache(maxsize=None) +def _altitude_ratios(ref_radius, altitude_radius, n): + a_over_r = ref_radius / altitude_radius + return tuple(a_over_r**(i + 2) for i in range(n)) diff --git a/geomag/world_magnetic_model.py b/geomag/world_magnetic_model.py new file mode 100644 index 0000000..fb4bff5 --- /dev/null +++ b/geomag/world_magnetic_model.py @@ -0,0 +1,287 @@ +from __future__ import division + +import math +import os +from datetime import date + +from .scalar_potential import scalar_potential, schmidt_quasi_normalisation, recursion_constants +from .latlon import LatLon + + +def _gen_square_array(size_x, default=None): + """Creates a square array with x by x elements""" + return [[default] * size_x for _ in range(size_x)] + + +def _calculate_decimal_year(date_of_year): + """ + .total_seconds() call makes the function 2/3 compliant. timedelta division was added in 3.2. + This does limit the module to 2.7 as the minimum supported + """ + year = date_of_year.year + start_of_this_year = date(year, 1, 1) + days_this_year = (date(year + 1, 1, 1) - start_of_this_year).total_seconds() + days_into_year = (date_of_year - start_of_this_year).total_seconds() + return year + days_into_year / days_this_year + + +def _convert_to_km(value, unit): + conversion_factor = {'ft': 3280.8399, 'm': 1000, 'km': 1} + try: + value_in_km = value / conversion_factor[unit] + except KeyError: + raise KeyError('Unknown unit: {unit}') + return value_in_km + +DEFAULT_PATH = os.path.join(os.path.dirname(__file__), 'model_data/WMM.COF') + + +class MagneticModelData(object): + """ Class to hold the data and calculations from the file.""" + _last_calculated_datetime = None + max_order = degree_of_expansion = 12 + array_size = max_order + 1 + coeff = {} + coefficient = _gen_square_array(array_size, 0.0) + coefficient_dot = _gen_square_array(array_size, 0.0) + time_adjusted_coefficients_cache = _gen_square_array(array_size, 0.0) + + def __init__(self, file): + with open(file) as world_magnetic_model_file: + for line in world_magnetic_model_file: + linevals = line.strip().split() + if len(linevals) == 3: + self.epoch = float(linevals[0]) + self.model = linevals[1] + self.modeldate = linevals[2] + elif len(linevals) == 6: + degree_n = int(float(linevals[0])) + order_m = int(float(linevals[1])) + gauss_g = float(linevals[2]) + gauss_h = float(linevals[3]) + gauss_g_dot = float(linevals[4]) + gauss_h_dot = float(linevals[5]) + self.coefficient[order_m][degree_n] = gauss_g + self.coefficient_dot[order_m][degree_n] = gauss_g_dot + if order_m != 0: + self.coefficient[degree_n][order_m - 1] = gauss_h + self.coefficient_dot[degree_n][order_m - 1] = gauss_h_dot + self._unnormalise_gauss_coefficients() + + def _unnormalise_gauss_coefficients(self): + """ Convert Schmidt normalized Gauss coefficients to unnormalized """ + schmidt_norm = schmidt_quasi_normalisation(self.array_size) + for n in range(self.array_size): + for m in range(self.array_size): + if m <= n: + self.coefficient[m][n] = schmidt_norm[m][n] * self.coefficient[m][n] + self.coefficient_dot[m][n] = schmidt_norm[m][n] * self.coefficient_dot[m][n] + else: + self.coefficient[m][n] = schmidt_norm[n + 1][m] * self.coefficient[m][n] + self.coefficient_dot[m][n] = schmidt_norm[n + 1][m] * self.coefficient_dot[m][n] + + def time_adjust_gauss(self, time): + """Time adjust the Gauss Coefficients + + There is a very basic cache happening here where if the time delta + hasn't changed then the previous calculation is returned + """ + current_delta_time = _calculate_decimal_year(time) - self.epoch + if self._last_calculated_datetime != current_delta_time: + self._update_time_coefficients(current_delta_time) + self._last_calculated_datetime = current_delta_time + return self.time_adjusted_coefficients_cache + + def _update_time_coefficients(self, delta_time): + for n in range(1, self.max_order + 1): + for m in range(0, n + 1): + self.time_adjusted_coefficients_cache[m][n] = (self.coefficient[m][n] + + delta_time * self.coefficient_dot[m][n]) + if m != 0: + self.time_adjusted_coefficients_cache[n][m - 1] = (self.coefficient[n][m - 1] + + delta_time * self.coefficient_dot[n][m - 1]) + + +class WorldMagneticModel(object): + """Class for calculating geomagnetic variation according to the world magnetic model + + Example Usage: + + >>> from geomag import WorldMagneticModel + >>> wmm = WorldMagneticModel() + >>> wmm.calc_mag_field(80,0).declination + + """ + + radius_earth = 6371200 + + _northerly_intensity = None + _easterly_intensity = None + _vertical_intensity = None + _horizontal_intensity = None + _total_intensity = None + _declination = None + _inclination = None + _grid_variation = None + + def __init__(self, world_magnetic_model_filename=DEFAULT_PATH): + """__init__(self,world_magnetic_model_filename='WMM.COF') + Loads a file containing the constants for the magnetic model. + + The coefficients for the model are included for the year 2015 if no + variable is provided. + File should be in the format of: + + ========= ======== ======== ======== ======== ======== + degree(n) order(m) g |mnt0| h |mnt0| g |mnt0| h |mnt0| + ========= ======== ======== ======== ======== ======== + + .. |mnt0| replace:: \ :sup:`m`:sub:`n`\ (t\ :sub:`0`\ ) + + """ + self.data = MagneticModelData(world_magnetic_model_filename) + self.k = recursion_constants(self.data.array_size) + + @property + def grid_variation(self): + return self._grid_variation + + @property + def dip(self): + return self._inclination + + @property + def declination(self): + return self._declination + + @property + def total_intensity(self): + return self._total_intensity + + @property + def Bh(self): + return self._horizontal_intensity + + @property + def Bx(self): + return self._northerly_intensity + + @property + def By(self): + return self._easterly_intensity + + @property + def Bz(self): + return self._vertical_intensity + + def calc_mag_field(self, dlat, dlon, altitude=0, date=date.today(), unit='ft'): + """calc_mag_field(self, dlat, dlon, altitude=0, date=date.today(), unit='ft') + + Calculates the magnetic field for a given latitude and longitude in decimal degrees. + + **Parameters** + + dlat + Latitude in degrees + dlon + Longitude in degrees + altitude : optional + Altitude at which to evaluate magnetic field + date : datetime.date - optional + Time will default to today + unit : {'ft','m','km'} - optional + Unit for altitude + + """ + altitude_in_km = _convert_to_km(altitude, unit) + + if not -1 <= altitude_in_km <= 850: + raise ValueError('World Magnetic Model is not valid outside the rage -1 to 850km') + + lat_lon = LatLon(dlat, dlon) + spherical_latitude, _, radial = lat_lon.convert_spherical(altitude_in_km * 1000) + + b_radius, b_theta, b_phi = scalar_potential(self.data.time_adjust_gauss(date), + lat_lon.lon_rad, + spherical_latitude, + self.data.array_size, + self.data.array_size, + radial, + k=self.k) + # Matching the method in the document + b_radius = -b_radius + b_theta = -b_theta + # /* + # ROTATE MAGNETIC VECTOR COMPONENTS FROM SPHERICAL TO + # GEODETIC COORDINATES + # */ + delta_latitude_radians = spherical_latitude - lat_lon.lat_rad + sin_delta_latitude = math.sin(delta_latitude_radians) + cos_delta_latitude = math.cos(delta_latitude_radians) + + self._northerly_intensity = b_theta * cos_delta_latitude - b_radius * sin_delta_latitude + self._easterly_intensity = b_phi + self._vertical_intensity = b_theta * sin_delta_latitude + b_radius * cos_delta_latitude + self._horizontal_intensity = math.sqrt(self._northerly_intensity ** 2 + self._easterly_intensity**2) + self._total_intensity = math.sqrt(self._horizontal_intensity**2 + self._vertical_intensity**2) + self._declination = math.degrees(math.atan2(self._easterly_intensity, self._northerly_intensity)) + self._inclination = math.degrees(math.atan2(self._vertical_intensity, self._horizontal_intensity)) + self._grid_variation = self._calculate_grid_variation(lat_lon) + return self + + def _calculate_grid_variation(self, lat_lon): + """Calculate the magnetic grid variation + + Compute magnetic grid variation if the current + geodetic position is in the arctic or antarctic + (i.e. glat > +55 degrees or glat < -55 degrees) + Otherwise, set magnetic grid variation to 0 + """ + cur_lat = lat_lon.lat + grid_variation = self._declination + if cur_lat > 55: + grid_variation -= lat_lon.lon + elif cur_lat < -55: + grid_variation += lat_lon.lon + grid_variation %= 360 + return grid_variation + + def mag_heading(self, hdg): + """Calculates the magnetic heading from a true heading. + hdg = true heading in degrees + Example Usage: + + >>> wmm.calc_mag_field(0,80).mag_heading(20) + + """ + return (hdg - self.declination + 360.0) % 360 + + def field_vectors(self): + """ Returns the main magnetic components: X, Y and Z""" + return self.return_array('X', 'Y', 'Z') + + def return_array(self, *variable_list): + """Populates a return array with the requested variables, where: + + | X = Northerly intensity + | Y = Easterly intensity + | Z = Vertical intensity + | H = Horizontal Intensity + | F = Total Intensity + | I = Inclination + | D = Dip + | GV = Grid Variation + + Example Usage: + + >>> wmm.calc_mag_field(0,80).return_array('GV','I') + + """ + variable_dict = { + 'X': self.Bx, 'Y': self.By, 'Z': self.Bz, 'H': self.Bh, + 'F': self.total_intensity, 'I': self.dip, 'D': self.dip, + 'GV': self.grid_variation} + return_list = [] + for variable in list(variable_list): + return_list.append(variable_dict[variable]) + return return_list diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..a5313f2 --- /dev/null +++ b/requirements.txt @@ -0,0 +1,2 @@ +py==1.4.31 +pytest==3.0.3 diff --git a/geomag/setup.py b/setup.py similarity index 88% rename from geomag/setup.py rename to setup.py index dec8558..9a314e0 100644 --- a/geomag/setup.py +++ b/setup.py @@ -1,4 +1,12 @@ from distutils.core import setup + +import sys + +if sys.version_info >= (3, 2): + install_requires = [] +else: + install_requires = ["functools32"] + setup( name = "geomag", packages = ["geomag"], @@ -11,6 +19,7 @@ url = "http://geomag.googlecode.com/", download_url = "//pypi.python.org/packages/source/g/geomag/geomag-0.9.2015.zip", keywords = ["magnetic", "variation", "declination"], + install_requires = install_requires, classifiers = [ "Programming Language :: Python", "Development Status :: 4 - Beta", diff --git a/tests/test_geomag.py b/tests/test_geomag.py new file mode 100644 index 0000000..74b632a --- /dev/null +++ b/tests/test_geomag.py @@ -0,0 +1,69 @@ +import pytest +from collections import namedtuple +from datetime import date + +from geomag import WorldMagneticModel + +DATE_IN_2015 = date(2015, 1, 1) +DATE_IN_2017 = date(2017, 7, 2) + +WMM = WorldMagneticModel() + +RESULTS = namedtuple('RESULTS', ['X', 'Y', 'Z', 'H', 'F', 'I', 'D', 'GV']) +CASES = [([80, 0, 0, DATE_IN_2015], + RESULTS(6627.1, -445.9, 54432.3, 6642.1, 54836, 83.04, -3.85, -3.85)), + ([0, 120, 0, DATE_IN_2015], + RESULTS(39518.2, 392.9, -11252.4, 39520.2, 41090.9, -15.89, 0.57, 0.57)), + ([-80, 240, 0, DATE_IN_2015], + RESULTS(5797.3, 15761.1, -52919.1, 16793.5, 55519.8, -72.39, 69.81, 309.81)), + ([80, 0, 328083.99, DATE_IN_2015], + RESULTS(6314.3, -471.6, 52269.8, 6331.9, 52652, 83.09, -4.27, -4.27)), + ([0, 120, 328083.99, DATE_IN_2015], + RESULTS(37535.6, 364.4, -10773.4, 37537.3, 39052.7, -16.01, 0.56, 0.56)), + ([-80, 240, 328083.99, DATE_IN_2015], + RESULTS(5613.1, 14791.5, -50378.6, 15820.7, 52804.4, -72.57, 69.22, 309.22)), + ([80, 0, 0, DATE_IN_2017], + RESULTS(6599.4, -317.1, 54459.2, 6607, 54858.5, 83.08, -2.75, -2.75)), + ([0, 120, 0, DATE_IN_2017], + RESULTS(39571.4, 222.5, -11030.1, 39572, 41080.5, -15.57, 0.32, 0.32)), + ([-80, 240, 0, DATE_IN_2017], + RESULTS(5873.8, 15781.4, -52687.9, 16839.1, 55313.4, -72.28, 69.58, 309.58)), + ([80, 0, 328083.99, DATE_IN_2017], + RESULTS(6290.5, -348.5, 52292.7, 6300.1, 52670.9, 83.13, -3.17, -3.17)), + ([0, 120, 328083.99, DATE_IN_2017], + RESULTS(37585.5, 209.5, -10564.2, 37586.1, 39042.5, -15.7, 0.32, 0.32)), + ([-80, 240, 328083.99, DATE_IN_2017], + RESULTS(5683.5, 14808.8, -50163, 15862, 52611.1, -72.45, 69, 309)), + ] + + +@pytest.fixture(scope="class", params=CASES) +def setup_class(request): + print(request) + request.cls.results = WMM.calc_mag_field(*request.param[0]) + request.cls.expected = request.param[1] + + +@pytest.mark.usefixtures("setup_class") +class TestGeoMag: + + def test_northerly_intensity(self): + assert(self.expected.X - self.results.Bx < 1) + + def test_easterly_intensity(self): + assert(self.expected.Y - self.results.By < 1) + + def test_down_intensity(self): + assert(self.expected.Z - self.results.Bz < 1) + + def test_total_intensity(self): + assert(self.expected.F - self.results.total_intensity < 1) + + def test_inclination(self): + assert(self.expected.I - self.results.dip < 0.01) + + def test_declination(self): + assert(self.expected.D - self.results.declination < 0.01) + + def test_grid_variation(self): + assert(self.expected.GV - self.results.grid_variation < 0.01) diff --git a/tests/test_latlon.py b/tests/test_latlon.py new file mode 100644 index 0000000..2d7a689 --- /dev/null +++ b/tests/test_latlon.py @@ -0,0 +1,42 @@ +import pytest + +from geomag.latlon import normalise_plus_minus_range +from geomag.latlon import LatLon + + +@pytest.mark.parametrize(("lat_lon", "expected_lat_lon"), [ + ((0, 0), (0, 0)), + ((-91, 181), (89, -179)), + ((0, 540), (0, 180)), +]) +def test_initilisation(lat_lon, expected_lat_lon): + location = LatLon(*lat_lon) + assert(location.lat == expected_lat_lon[0]) + assert(location.lon == expected_lat_lon[1]) + + +def test_error_on_wrong_unit(): + with pytest.raises(Exception): + LatLon((0, 0), (0, 0), 'wrong_unit') + + +@pytest.mark.parametrize(("value", "norm_range", "expected"), [ + (0, 'lat', 0), + (0, 'lon', 0), + (0, 10, 0), + (20, 10, 0), + (20, 10, 0), + (15, 10, -5), + (91, 'lat', -89), + (180, 'lon', 180), + (-180, 'lon', -180), + (181, 'lon', -179), + (-181, 'lon', 179), + (360, 'lon', 0), + (-360, 'lon', 0), + (-185, 'lon', 175), + (540, 'lon', 180), + (0, 0, 0), +]) +def test_norm_range(value, norm_range, expected): + assert expected == normalise_plus_minus_range(value, norm_range) diff --git a/tests/test_scalar_potential.py b/tests/test_scalar_potential.py new file mode 100644 index 0000000..a3d75ac --- /dev/null +++ b/tests/test_scalar_potential.py @@ -0,0 +1,33 @@ +import pytest + +from geomag.scalar_potential import legrende_polynomials, associated_polynomials, double_factorial + + +@pytest.mark.parametrize(("n", "expected"), [ + (1, 1), + (2, 2), + (3, 3), + (4, 8), + (5, 15), +]) +def test_double_factorial(n, expected): + assert double_factorial(n) == expected + + +@pytest.mark.parametrize(("in_array", "expected_poly"), [ + ((1, 1), [1, 1]), + ((2, 1), [1, 1, 1]), + ((3, 1), [1, 1, 1, 1]), + ((1, 2), [1, 2]), + ((2, 2), [1, 2, 5.5]), + ((3, 2), [1, 2, 5.5, 17]), +]) +def test_polynomials(in_array, expected_poly): + assert legrende_polynomials(*in_array) == pytest.approx(expected_poly) + + +@pytest.mark.parametrize(("in_array", "expected_poly"), [ + (({'cos_theta': 1, 'sin_theta': 1, 'max_n': 1, 'max_m': 1}), ([[1]], [[0]])), +]) +def test_associated_polynomials(in_array, expected_poly): + assert associated_polynomials(**in_array) == expected_poly