From ac58ab00f52deea11abd6773675fbbbb94dbaba4 Mon Sep 17 00:00:00 2001 From: Dimitrios Theodorakis Date: Thu, 26 Aug 2021 15:15:37 +0100 Subject: [PATCH] planetMagFields integration Integrate planetMagFields and require scipy --- MANIFEST.in | 5 +- README.md | 37 +- setup.cfg | 3 +- src/astroedu/__main__.py | 60 +- src/astroedu/planetmagfields/LICENSE | 674 ++++++++++++++++++ src/astroedu/planetmagfields/README.md | 312 ++++++++ src/astroedu/planetmagfields/__init__.py | 5 + src/astroedu/planetmagfields/data/IGRF13.dat | 195 +++++ .../planetmagfields/data/ganymede.dat | 10 + src/astroedu/planetmagfields/data/jupiter.dat | 120 ++++ src/astroedu/planetmagfields/data/mercury.dat | 4 + src/astroedu/planetmagfields/data/neptune.dat | 18 + src/astroedu/planetmagfields/data/saturn.dat | 14 + src/astroedu/planetmagfields/data/uranus.dat | 18 + src/astroedu/planetmagfields/libbfield.py | 93 +++ src/astroedu/planetmagfields/libgauss.py | 349 +++++++++ src/astroedu/planetmagfields/planet.py | 160 +++++ src/astroedu/planetmagfields/plotlib.py | 146 ++++ src/astroedu/planetmagfields/potextra.py | 120 ++++ 19 files changed, 2328 insertions(+), 15 deletions(-) create mode 100644 src/astroedu/planetmagfields/LICENSE create mode 100644 src/astroedu/planetmagfields/README.md create mode 100644 src/astroedu/planetmagfields/__init__.py create mode 100644 src/astroedu/planetmagfields/data/IGRF13.dat create mode 100644 src/astroedu/planetmagfields/data/ganymede.dat create mode 100644 src/astroedu/planetmagfields/data/jupiter.dat create mode 100644 src/astroedu/planetmagfields/data/mercury.dat create mode 100644 src/astroedu/planetmagfields/data/neptune.dat create mode 100644 src/astroedu/planetmagfields/data/saturn.dat create mode 100644 src/astroedu/planetmagfields/data/uranus.dat create mode 100644 src/astroedu/planetmagfields/libbfield.py create mode 100644 src/astroedu/planetmagfields/libgauss.py create mode 100644 src/astroedu/planetmagfields/planet.py create mode 100644 src/astroedu/planetmagfields/plotlib.py create mode 100644 src/astroedu/planetmagfields/potextra.py diff --git a/MANIFEST.in b/MANIFEST.in index 363c66d..481c25f 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -1,2 +1,5 @@ include src/astroedu/interactives/* -include src/astroedu/datasets/* \ No newline at end of file +include src/astroedu/datasets/* + +include src/astroedu/planetmagfields/* +include src/astroedu/planetmagfields/data/* \ No newline at end of file diff --git a/README.md b/README.md index 63851b9..3b6f387 100644 --- a/README.md +++ b/README.md @@ -3,9 +3,14 @@ [![PyPI version](https://badge.fury.io/py/astroedu.svg)](https://badge.fury.io/py/astroedu) -[![Contributor Covenant](https://img.shields.io/badge/Contributor%20Covenant-2.0-4baaaa.svg)](code_of_conduct.md) +[![Contributor Covenant](https://img.shields.io/badge/Contributor%20Covenant-2.0-4baaaa.svg)](code_of_conduct.md) **Alpha** -This package is in alpha. +A package containing small interactives, datasets, functions etc for teaching astronomy. + +Contributors: + +- Dimitrios Theodorakis +- Ankit Barik, planetMagFields module, [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.4706157.svg)](https://doi.org/10.5281/zenodo.4706157), [github.com/AnkitBarik/planetMagFields](github.com/AnkitBarik/planetMagFields) ## Installation @@ -13,14 +18,6 @@ To install **astroedu** run ``` >>> pip install astroedu ``` -then run -``` ->>> astroedu build -``` -which creates the configuration file **config.ini**. - -**config.ini** stores the location of your documents folder. -This is where an astroedu folder for interactive activites will me made. (not yet implemented) ## Current Functionality @@ -120,4 +117,22 @@ The class has built in methods. For instance to calculate the tides on the Earth >>> forces = earth.tides(moon, step=0.25, scale=5.972*10**24) ``` Documentation coming soon. -More methods will be added at a later date including calculating gravitational potentials and plotting tides & potentials. \ No newline at end of file +More methods will be added at a later date including calculating gravitational potentials and plotting tides & potentials. + +## Submodules + +### planetMagFields + +planetMagFields by Ankit Barik. +[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.4706157.svg)](https://doi.org/10.5281/zenodo.4706157), [github.com/AnkitBarik/planetMagFields](github.com/AnkitBarik/planetMagFields) + +See Ankit's GitHub for usage. Note: cartopy is required for some plots which requires these non-python packages to be installed, [GEOS](https://trac.osgeo.org/geos/) and [PROJ](https://proj.org/). Some functions also require other libraries, see Ankit's GitHub for more info. + +Since the dataset location is defined relative to the astroedu install there is no need to specify a datDir for instance: + +``` +>>> from astroedu.planetmagfields import * +>>> p = planet(name='jupiter') +>>> # not p = planet(name='jupiter',datDir='planetmagfields/data/') +>>> p.plot(r=0.85,proj='Mollweide') +``` \ No newline at end of file diff --git a/setup.cfg b/setup.cfg index 664a094..46edd0b 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,7 +1,7 @@ [metadata] # replace with your username: name = astroedu -version = 0.3.0a1 +version = 0.4.0a0 author = Dimitrios Theodorakis author_email = astrodimitrios@gmail.com description = A python package for astronomy educators @@ -25,6 +25,7 @@ python_requires = >=3.7 include_package_data = True install_requires = numpy >= 1.20.0 + scipy >= 1.6.0 matplotlib >= 3.3.0 pandas >= 1.1 jupyterlab >= 3.0.0 diff --git a/src/astroedu/__main__.py b/src/astroedu/__main__.py index 398ef03..2b9a8f4 100644 --- a/src/astroedu/__main__.py +++ b/src/astroedu/__main__.py @@ -7,7 +7,7 @@ from astroedu.__build__ import get_astroedu_path from astroedu.__build__ import build_path_config -commands = ['build', 'interactive'] +commands = ['build', 'interactive', 'magfield'] astroedu_path = get_astroedu_path() config_file_name = '/config.ini' config_file_path = astroedu_path + config_file_name @@ -46,6 +46,60 @@ def _load_interactive(*args): print(f'Loading interactive {interactive_name}!') subprocess.run(run_interactive) +def _planetmagfield_quick(*args): + ''' Utility function to make wuick plots using planetMagFields + See https://zenodo.org/record/5140421#.YSdp7t_TUuU + ''' + levels=20 + cmap='RdBu_r' + proj = 'Mollweide' + r = 1 + + import numpy as np + + if len(args) > 4: + raise ValueError("Too many arguments, exiting ...\n") + elif len(args) == 4: + planet = str(args[1]).lower() + r = np.float32(args[2]) + proj = str(args[3]) + elif len(args) == 3: + planet = str(args[1]).lower() + try: + r = np.float32(args[2]) + except: + proj = str(args[2]) + elif len(args) == 2: + if args[2] == '--help': + print('planetMagFields\nSubmodule by Ankit Barik - 10.5281/zenodo.5140421\nUsage see - github.com/AnkitBarik/planetMagFields') + else: + print("Radius not specified, using surface\n") + planet = str(args[1]).lower() + else: + print("Planet or radius not specified, plotting for Earth's surface\n") + planet="earth" + r=1. + + import matplotlib.pyplot as plt + from astroedu.planetmagfields.libbfield import getBr, plotAllFields, plotMagField + + if planet == 'all': + with files('astroedu.planetmagfields').joinpath('data/') as p: + plotAllFields(datDir=p,r=r,levels=levels,cmap=cmap,proj=proj) + plt.tight_layout() + plt.subplots_adjust(top=0.895, + bottom=0.035, + left=0.023, + right=0.976, + hspace=0.38, + wspace=0.109) + else: + with files('astroedu.planetmagfields').joinpath('data/') as p: + plotMagField(planet=planet,r=r,datDir=p,levels=levels,cmap=cmap,proj=proj) + plt.tight_layout() + + plt.show() + def startup(*args): ''' Startup function called by Fire cli function main() @@ -67,7 +121,7 @@ def startup(*args): ''' if not args: raise ValueError('No arguments passed to astroedu call.') - command = args[0] + command = args[0].lower() if command not in commands: raise ValueError('Command not recognised.') elif command == 'build': @@ -76,6 +130,8 @@ def startup(*args): raise ValueError(f'No arguments passed to {command}') elif command == 'interactive': _load_interactive(*args) + elif command == 'magfield': + _planetmagfield_quick(*args) if __name__ == '__main__': fire.Fire(startup) \ No newline at end of file diff --git a/src/astroedu/planetmagfields/LICENSE b/src/astroedu/planetmagfields/LICENSE new file mode 100644 index 0000000..f288702 --- /dev/null +++ b/src/astroedu/planetmagfields/LICENSE @@ -0,0 +1,674 @@ + GNU GENERAL PUBLIC LICENSE + Version 3, 29 June 2007 + + Copyright (C) 2007 Free Software Foundation, Inc. + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The GNU General Public License is a free, copyleft license for +software and other kinds of works. + + The licenses for most software and other practical works are designed +to take away your freedom to share and change the works. By contrast, +the GNU General Public License is intended to guarantee your freedom to +share and change all versions of a program--to make sure it remains free +software for all its users. We, the Free Software Foundation, use the +GNU General Public License for most of our software; it applies also to +any other work released this way by its authors. You can apply it to +your programs, too. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +them if you wish), that you receive source code or can get it if you +want it, that you can change the software or use pieces of it in new +free programs, and that you know you can do these things. + + To protect your rights, we need to prevent others from denying you +these rights or asking you to surrender the rights. Therefore, you have +certain responsibilities if you distribute copies of the software, or if +you modify it: responsibilities to respect the freedom of others. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must pass on to the recipients the same +freedoms that you received. You must make sure that they, too, receive +or can get the source code. And you must show them these terms so they +know their rights. + + Developers that use the GNU GPL protect your rights with two steps: +(1) assert copyright on the software, and (2) offer you this License +giving you legal permission to copy, distribute and/or modify it. + + For the developers' and authors' protection, the GPL clearly explains +that there is no warranty for this free software. For both users' and +authors' sake, the GPL requires that modified versions be marked as +changed, so that their problems will not be attributed erroneously to +authors of previous versions. + + Some devices are designed to deny users access to install or run +modified versions of the software inside them, although the manufacturer +can do so. This is fundamentally incompatible with the aim of +protecting users' freedom to change the software. The systematic +pattern of such abuse occurs in the area of products for individuals to +use, which is precisely where it is most unacceptable. Therefore, we +have designed this version of the GPL to prohibit the practice for those +products. If such problems arise substantially in other domains, we +stand ready to extend this provision to those domains in future versions +of the GPL, as needed to protect the freedom of users. + + Finally, every program is threatened constantly by software patents. +States should not allow patents to restrict development and use of +software on general-purpose computers, but in those that do, we wish to +avoid the special danger that patents applied to a free program could +make it effectively proprietary. To prevent this, the GPL assures that +patents cannot be used to render the program non-free. + + The precise terms and conditions for copying, distribution and +modification follow. + + TERMS AND CONDITIONS + + 0. Definitions. + + "This License" refers to version 3 of the GNU General Public License. + + "Copyright" also means copyright-like laws that apply to other kinds of +works, such as semiconductor masks. + + "The Program" refers to any copyrightable work licensed under this +License. Each licensee is addressed as "you". "Licensees" and +"recipients" may be individuals or organizations. + + To "modify" a work means to copy from or adapt all or part of the work +in a fashion requiring copyright permission, other than the making of an +exact copy. The resulting work is called a "modified version" of the +earlier work or a work "based on" the earlier work. + + A "covered work" means either the unmodified Program or a work based +on the Program. + + To "propagate" a work means to do anything with it that, without +permission, would make you directly or secondarily liable for +infringement under applicable copyright law, except executing it on a +computer or modifying a private copy. Propagation includes copying, +distribution (with or without modification), making available to the +public, and in some countries other activities as well. + + To "convey" a work means any kind of propagation that enables other +parties to make or receive copies. Mere interaction with a user through +a computer network, with no transfer of a copy, is not conveying. + + An interactive user interface displays "Appropriate Legal Notices" +to the extent that it includes a convenient and prominently visible +feature that (1) displays an appropriate copyright notice, and (2) +tells the user that there is no warranty for the work (except to the +extent that warranties are provided), that licensees may convey the +work under this License, and how to view a copy of this License. If +the interface presents a list of user commands or options, such as a +menu, a prominent item in the list meets this criterion. + + 1. Source Code. + + The "source code" for a work means the preferred form of the work +for making modifications to it. "Object code" means any non-source +form of a work. + + A "Standard Interface" means an interface that either is an official +standard defined by a recognized standards body, or, in the case of +interfaces specified for a particular programming language, one that +is widely used among developers working in that language. + + The "System Libraries" of an executable work include anything, other +than the work as a whole, that (a) is included in the normal form of +packaging a Major Component, but which is not part of that Major +Component, and (b) serves only to enable use of the work with that +Major Component, or to implement a Standard Interface for which an +implementation is available to the public in source code form. A +"Major Component", in this context, means a major essential component +(kernel, window system, and so on) of the specific operating system +(if any) on which the executable work runs, or a compiler used to +produce the work, or an object code interpreter used to run it. + + The "Corresponding Source" for a work in object code form means all +the source code needed to generate, install, and (for an executable +work) run the object code and to modify the work, including scripts to +control those activities. However, it does not include the work's +System Libraries, or general-purpose tools or generally available free +programs which are used unmodified in performing those activities but +which are not part of the work. For example, Corresponding Source +includes interface definition files associated with source files for +the work, and the source code for shared libraries and dynamically +linked subprograms that the work is specifically designed to require, +such as by intimate data communication or control flow between those +subprograms and other parts of the work. + + The Corresponding Source need not include anything that users +can regenerate automatically from other parts of the Corresponding +Source. + + The Corresponding Source for a work in source code form is that +same work. + + 2. Basic Permissions. + + All rights granted under this License are granted for the term of +copyright on the Program, and are irrevocable provided the stated +conditions are met. This License explicitly affirms your unlimited +permission to run the unmodified Program. The output from running a +covered work is covered by this License only if the output, given its +content, constitutes a covered work. This License acknowledges your +rights of fair use or other equivalent, as provided by copyright law. + + You may make, run and propagate covered works that you do not +convey, without conditions so long as your license otherwise remains +in force. You may convey covered works to others for the sole purpose +of having them make modifications exclusively for you, or provide you +with facilities for running those works, provided that you comply with +the terms of this License in conveying all material for which you do +not control copyright. Those thus making or running the covered works +for you must do so exclusively on your behalf, under your direction +and control, on terms that prohibit them from making any copies of +your copyrighted material outside their relationship with you. + + Conveying under any other circumstances is permitted solely under +the conditions stated below. Sublicensing is not allowed; section 10 +makes it unnecessary. + + 3. Protecting Users' Legal Rights From Anti-Circumvention Law. + + No covered work shall be deemed part of an effective technological +measure under any applicable law fulfilling obligations under article +11 of the WIPO copyright treaty adopted on 20 December 1996, or +similar laws prohibiting or restricting circumvention of such +measures. + + When you convey a covered work, you waive any legal power to forbid +circumvention of technological measures to the extent such circumvention +is effected by exercising rights under this License with respect to +the covered work, and you disclaim any intention to limit operation or +modification of the work as a means of enforcing, against the work's +users, your or third parties' legal rights to forbid circumvention of +technological measures. + + 4. Conveying Verbatim Copies. + + You may convey verbatim copies of the Program's source code as you +receive it, in any medium, provided that you conspicuously and +appropriately publish on each copy an appropriate copyright notice; +keep intact all notices stating that this License and any +non-permissive terms added in accord with section 7 apply to the code; +keep intact all notices of the absence of any warranty; and give all +recipients a copy of this License along with the Program. + + You may charge any price or no price for each copy that you convey, +and you may offer support or warranty protection for a fee. + + 5. Conveying Modified Source Versions. + + You may convey a work based on the Program, or the modifications to +produce it from the Program, in the form of source code under the +terms of section 4, provided that you also meet all of these conditions: + + a) The work must carry prominent notices stating that you modified + it, and giving a relevant date. + + b) The work must carry prominent notices stating that it is + released under this License and any conditions added under section + 7. This requirement modifies the requirement in section 4 to + "keep intact all notices". + + c) You must license the entire work, as a whole, under this + License to anyone who comes into possession of a copy. This + License will therefore apply, along with any applicable section 7 + additional terms, to the whole of the work, and all its parts, + regardless of how they are packaged. This License gives no + permission to license the work in any other way, but it does not + invalidate such permission if you have separately received it. + + d) If the work has interactive user interfaces, each must display + Appropriate Legal Notices; however, if the Program has interactive + interfaces that do not display Appropriate Legal Notices, your + work need not make them do so. + + A compilation of a covered work with other separate and independent +works, which are not by their nature extensions of the covered work, +and which are not combined with it such as to form a larger program, +in or on a volume of a storage or distribution medium, is called an +"aggregate" if the compilation and its resulting copyright are not +used to limit the access or legal rights of the compilation's users +beyond what the individual works permit. Inclusion of a covered work +in an aggregate does not cause this License to apply to the other +parts of the aggregate. + + 6. Conveying Non-Source Forms. + + You may convey a covered work in object code form under the terms +of sections 4 and 5, provided that you also convey the +machine-readable Corresponding Source under the terms of this License, +in one of these ways: + + a) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by the + Corresponding Source fixed on a durable physical medium + customarily used for software interchange. + + b) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by a + written offer, valid for at least three years and valid for as + long as you offer spare parts or customer support for that product + model, to give anyone who possesses the object code either (1) a + copy of the Corresponding Source for all the software in the + product that is covered by this License, on a durable physical + medium customarily used for software interchange, for a price no + more than your reasonable cost of physically performing this + conveying of source, or (2) access to copy the + Corresponding Source from a network server at no charge. + + c) Convey individual copies of the object code with a copy of the + written offer to provide the Corresponding Source. This + alternative is allowed only occasionally and noncommercially, and + only if you received the object code with such an offer, in accord + with subsection 6b. + + d) Convey the object code by offering access from a designated + place (gratis or for a charge), and offer equivalent access to the + Corresponding Source in the same way through the same place at no + further charge. You need not require recipients to copy the + Corresponding Source along with the object code. If the place to + copy the object code is a network server, the Corresponding Source + may be on a different server (operated by you or a third party) + that supports equivalent copying facilities, provided you maintain + clear directions next to the object code saying where to find the + Corresponding Source. Regardless of what server hosts the + Corresponding Source, you remain obligated to ensure that it is + available for as long as needed to satisfy these requirements. + + e) Convey the object code using peer-to-peer transmission, provided + you inform other peers where the object code and Corresponding + Source of the work are being offered to the general public at no + charge under subsection 6d. + + A separable portion of the object code, whose source code is excluded +from the Corresponding Source as a System Library, need not be +included in conveying the object code work. + + A "User Product" is either (1) a "consumer product", which means any +tangible personal property which is normally used for personal, family, +or household purposes, or (2) anything designed or sold for incorporation +into a dwelling. In determining whether a product is a consumer product, +doubtful cases shall be resolved in favor of coverage. For a particular +product received by a particular user, "normally used" refers to a +typical or common use of that class of product, regardless of the status +of the particular user or of the way in which the particular user +actually uses, or expects or is expected to use, the product. A product +is a consumer product regardless of whether the product has substantial +commercial, industrial or non-consumer uses, unless such uses represent +the only significant mode of use of the product. + + "Installation Information" for a User Product means any methods, +procedures, authorization keys, or other information required to install +and execute modified versions of a covered work in that User Product from +a modified version of its Corresponding Source. The information must +suffice to ensure that the continued functioning of the modified object +code is in no case prevented or interfered with solely because +modification has been made. + + If you convey an object code work under this section in, or with, or +specifically for use in, a User Product, and the conveying occurs as +part of a transaction in which the right of possession and use of the +User Product is transferred to the recipient in perpetuity or for a +fixed term (regardless of how the transaction is characterized), the +Corresponding Source conveyed under this section must be accompanied +by the Installation Information. But this requirement does not apply +if neither you nor any third party retains the ability to install +modified object code on the User Product (for example, the work has +been installed in ROM). + + The requirement to provide Installation Information does not include a +requirement to continue to provide support service, warranty, or updates +for a work that has been modified or installed by the recipient, or for +the User Product in which it has been modified or installed. Access to a +network may be denied when the modification itself materially and +adversely affects the operation of the network or violates the rules and +protocols for communication across the network. + + Corresponding Source conveyed, and Installation Information provided, +in accord with this section must be in a format that is publicly +documented (and with an implementation available to the public in +source code form), and must require no special password or key for +unpacking, reading or copying. + + 7. Additional Terms. + + "Additional permissions" are terms that supplement the terms of this +License by making exceptions from one or more of its conditions. +Additional permissions that are applicable to the entire Program shall +be treated as though they were included in this License, to the extent +that they are valid under applicable law. If additional permissions +apply only to part of the Program, that part may be used separately +under those permissions, but the entire Program remains governed by +this License without regard to the additional permissions. + + When you convey a copy of a covered work, you may at your option +remove any additional permissions from that copy, or from any part of +it. (Additional permissions may be written to require their own +removal in certain cases when you modify the work.) You may place +additional permissions on material, added by you to a covered work, +for which you have or can give appropriate copyright permission. + + Notwithstanding any other provision of this License, for material you +add to a covered work, you may (if authorized by the copyright holders of +that material) supplement the terms of this License with terms: + + a) Disclaiming warranty or limiting liability differently from the + terms of sections 15 and 16 of this License; or + + b) Requiring preservation of specified reasonable legal notices or + author attributions in that material or in the Appropriate Legal + Notices displayed by works containing it; or + + c) Prohibiting misrepresentation of the origin of that material, or + requiring that modified versions of such material be marked in + reasonable ways as different from the original version; or + + d) Limiting the use for publicity purposes of names of licensors or + authors of the material; or + + e) Declining to grant rights under trademark law for use of some + trade names, trademarks, or service marks; or + + f) Requiring indemnification of licensors and authors of that + material by anyone who conveys the material (or modified versions of + it) with contractual assumptions of liability to the recipient, for + any liability that these contractual assumptions directly impose on + those licensors and authors. + + All other non-permissive additional terms are considered "further +restrictions" within the meaning of section 10. If the Program as you +received it, or any part of it, contains a notice stating that it is +governed by this License along with a term that is a further +restriction, you may remove that term. If a license document contains +a further restriction but permits relicensing or conveying under this +License, you may add to a covered work material governed by the terms +of that license document, provided that the further restriction does +not survive such relicensing or conveying. + + If you add terms to a covered work in accord with this section, you +must place, in the relevant source files, a statement of the +additional terms that apply to those files, or a notice indicating +where to find the applicable terms. + + Additional terms, permissive or non-permissive, may be stated in the +form of a separately written license, or stated as exceptions; +the above requirements apply either way. + + 8. Termination. + + You may not propagate or modify a covered work except as expressly +provided under this License. Any attempt otherwise to propagate or +modify it is void, and will automatically terminate your rights under +this License (including any patent licenses granted under the third +paragraph of section 11). + + However, if you cease all violation of this License, then your +license from a particular copyright holder is reinstated (a) +provisionally, unless and until the copyright holder explicitly and +finally terminates your license, and (b) permanently, if the copyright +holder fails to notify you of the violation by some reasonable means +prior to 60 days after the cessation. + + Moreover, your license from a particular copyright holder is +reinstated permanently if the copyright holder notifies you of the +violation by some reasonable means, this is the first time you have +received notice of violation of this License (for any work) from that +copyright holder, and you cure the violation prior to 30 days after +your receipt of the notice. + + Termination of your rights under this section does not terminate the +licenses of parties who have received copies or rights from you under +this License. If your rights have been terminated and not permanently +reinstated, you do not qualify to receive new licenses for the same +material under section 10. + + 9. Acceptance Not Required for Having Copies. + + You are not required to accept this License in order to receive or +run a copy of the Program. Ancillary propagation of a covered work +occurring solely as a consequence of using peer-to-peer transmission +to receive a copy likewise does not require acceptance. However, +nothing other than this License grants you permission to propagate or +modify any covered work. These actions infringe copyright if you do +not accept this License. Therefore, by modifying or propagating a +covered work, you indicate your acceptance of this License to do so. + + 10. Automatic Licensing of Downstream Recipients. + + Each time you convey a covered work, the recipient automatically +receives a license from the original licensors, to run, modify and +propagate that work, subject to this License. You are not responsible +for enforcing compliance by third parties with this License. + + An "entity transaction" is a transaction transferring control of an +organization, or substantially all assets of one, or subdividing an +organization, or merging organizations. If propagation of a covered +work results from an entity transaction, each party to that +transaction who receives a copy of the work also receives whatever +licenses to the work the party's predecessor in interest had or could +give under the previous paragraph, plus a right to possession of the +Corresponding Source of the work from the predecessor in interest, if +the predecessor has it or can get it with reasonable efforts. + + You may not impose any further restrictions on the exercise of the +rights granted or affirmed under this License. For example, you may +not impose a license fee, royalty, or other charge for exercise of +rights granted under this License, and you may not initiate litigation +(including a cross-claim or counterclaim in a lawsuit) alleging that +any patent claim is infringed by making, using, selling, offering for +sale, or importing the Program or any portion of it. + + 11. Patents. + + A "contributor" is a copyright holder who authorizes use under this +License of the Program or a work on which the Program is based. The +work thus licensed is called the contributor's "contributor version". + + A contributor's "essential patent claims" are all patent claims +owned or controlled by the contributor, whether already acquired or +hereafter acquired, that would be infringed by some manner, permitted +by this License, of making, using, or selling its contributor version, +but do not include claims that would be infringed only as a +consequence of further modification of the contributor version. For +purposes of this definition, "control" includes the right to grant +patent sublicenses in a manner consistent with the requirements of +this License. + + Each contributor grants you a non-exclusive, worldwide, royalty-free +patent license under the contributor's essential patent claims, to +make, use, sell, offer for sale, import and otherwise run, modify and +propagate the contents of its contributor version. + + In the following three paragraphs, a "patent license" is any express +agreement or commitment, however denominated, not to enforce a patent +(such as an express permission to practice a patent or covenant not to +sue for patent infringement). To "grant" such a patent license to a +party means to make such an agreement or commitment not to enforce a +patent against the party. + + If you convey a covered work, knowingly relying on a patent license, +and the Corresponding Source of the work is not available for anyone +to copy, free of charge and under the terms of this License, through a +publicly available network server or other readily accessible means, +then you must either (1) cause the Corresponding Source to be so +available, or (2) arrange to deprive yourself of the benefit of the +patent license for this particular work, or (3) arrange, in a manner +consistent with the requirements of this License, to extend the patent +license to downstream recipients. "Knowingly relying" means you have +actual knowledge that, but for the patent license, your conveying the +covered work in a country, or your recipient's use of the covered work +in a country, would infringe one or more identifiable patents in that +country that you have reason to believe are valid. + + If, pursuant to or in connection with a single transaction or +arrangement, you convey, or propagate by procuring conveyance of, a +covered work, and grant a patent license to some of the parties +receiving the covered work authorizing them to use, propagate, modify +or convey a specific copy of the covered work, then the patent license +you grant is automatically extended to all recipients of the covered +work and works based on it. + + A patent license is "discriminatory" if it does not include within +the scope of its coverage, prohibits the exercise of, or is +conditioned on the non-exercise of one or more of the rights that are +specifically granted under this License. You may not convey a covered +work if you are a party to an arrangement with a third party that is +in the business of distributing software, under which you make payment +to the third party based on the extent of your activity of conveying +the work, and under which the third party grants, to any of the +parties who would receive the covered work from you, a discriminatory +patent license (a) in connection with copies of the covered work +conveyed by you (or copies made from those copies), or (b) primarily +for and in connection with specific products or compilations that +contain the covered work, unless you entered into that arrangement, +or that patent license was granted, prior to 28 March 2007. + + Nothing in this License shall be construed as excluding or limiting +any implied license or other defenses to infringement that may +otherwise be available to you under applicable patent law. + + 12. No Surrender of Others' Freedom. + + If conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot convey a +covered work so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you may +not convey it at all. For example, if you agree to terms that obligate you +to collect a royalty for further conveying from those to whom you convey +the Program, the only way you could satisfy both those terms and this +License would be to refrain entirely from conveying the Program. + + 13. Use with the GNU Affero General Public License. + + Notwithstanding any other provision of this License, you have +permission to link or combine any covered work with a work licensed +under version 3 of the GNU Affero General Public License into a single +combined work, and to convey the resulting work. The terms of this +License will continue to apply to the part which is the covered work, +but the special requirements of the GNU Affero General Public License, +section 13, concerning interaction through a network will apply to the +combination as such. + + 14. Revised Versions of this License. + + The Free Software Foundation may publish revised and/or new versions of +the GNU General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + + Each version is given a distinguishing version number. If the +Program specifies that a certain numbered version of the GNU General +Public License "or any later version" applies to it, you have the +option of following the terms and conditions either of that numbered +version or of any later version published by the Free Software +Foundation. If the Program does not specify a version number of the +GNU General Public License, you may choose any version ever published +by the Free Software Foundation. + + If the Program specifies that a proxy can decide which future +versions of the GNU General Public License can be used, that proxy's +public statement of acceptance of a version permanently authorizes you +to choose that version for the Program. + + Later license versions may give you additional or different +permissions. However, no additional obligations are imposed on any +author or copyright holder as a result of your choosing to follow a +later version. + + 15. Disclaimer of Warranty. + + THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY +APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT +HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY +OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, +THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM +IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF +ALL NECESSARY SERVICING, REPAIR OR CORRECTION. + + 16. Limitation of Liability. + + IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS +THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY +GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE +USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF +DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD +PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), +EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF +SUCH DAMAGES. + + 17. Interpretation of Sections 15 and 16. + + If the disclaimer of warranty and limitation of liability provided +above cannot be given local legal effect according to their terms, +reviewing courts shall apply local law that most closely approximates +an absolute waiver of all civil liability in connection with the +Program, unless a warranty or assumption of liability accompanies a +copy of the Program in return for a fee. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +state the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + + Copyright (C) + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + +Also add information on how to contact you by electronic and paper mail. + + If the program does terminal interaction, make it output a short +notice like this when it starts in an interactive mode: + + Copyright (C) + This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, your program's commands +might be different; for a GUI interface, you would use an "about box". + + You should also get your employer (if you work as a programmer) or school, +if any, to sign a "copyright disclaimer" for the program, if necessary. +For more information on this, and how to apply and follow the GNU GPL, see +. + + The GNU General Public License does not permit incorporating your program +into proprietary programs. If your program is a subroutine library, you +may consider it more useful to permit linking proprietary applications with +the library. If this is what you want to do, use the GNU Lesser General +Public License instead of this License. But first, please read +. diff --git a/src/astroedu/planetmagfields/README.md b/src/astroedu/planetmagfields/README.md new file mode 100644 index 0000000..0fba5ce --- /dev/null +++ b/src/astroedu/planetmagfields/README.md @@ -0,0 +1,312 @@ +# planetMagFields +[![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0) [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.4706157.svg)](https://doi.org/10.5281/zenodo.4706157) + +Routines to easily access information about magnetic fields of planets in our solar system and visualize them in both 2D and 3D. These require [NumPy](https://numpy.org/), [Matplotlib](https://matplotlib.org/) and [SciPy](https://www.scipy.org/) (pronounced "Sigh Pie"). Other than that, the following external libraries are used for a few different functions: + + - 2D plotting for map projections other than Hammer : [Cartopy](https://scitools.org.uk/cartopy/docs/latest/) library ( see more under [Projections](#projections) ) + - Potential extrapolation: [SHTns](https://bitbucket.org/nschaeff/shtns) library + - Writing vts files for 3D visualisation: [PyEVTK](https://github.com/paulo-herrera/PyEVTK) library + +# The `planet` class + +This gives access to all the relevant properties of a planet and has methods to plot +the field and write a `vts` file for 3D visualization. Usage: + +```python +from planetmagfields import * +p = planet(name='earth',datDir='planetmagfields/data/') +``` + +This displays the some information about the planet + +``` +Planet: Earth +l_max = 13 +Dipole tilt (degrees) = -9.410531 +``` + +and gives access to +variables associated with the planet such as: + + - `p.lmax` : maximum spherical harmonic degree till which data is available + - `p.glm`, `p.hlm`: the Gauss coefficients + - `p.Br` : computed radial magnetic field at surface + - `p.dipTheta` : dipole tilt with respect to the rotation axis + - `p.dipPhi` : dipole longitude ( in case zero longitude is known, applicable to Earth ) + - `p.idx` : indices to get values of Gauss coefficients + +Example: + +```python +In [1]: from planetmagfields import * + +In [2]: p = planet(name='jupiter') +Planet: Jupiter +l_max = 10 +Dipole tilt (degrees) = 10.307870 + +In [3]: p.glm[p.idx[2,0]] # g20 +Out[3]: 11670.4 + +In [4]: p.hlm[p.idx[4,2]] # h42 +Out[4]: 27811.2 +``` + +as well as the functions: + +## `planet.plot()` + +This function plots a 2D surface plot of the radial magnetic field at radius `r` given in terms of the surface radius. +For example, + +```python +from planetmagfields import * +p = planet(name='earth',datDir='planetmagfields/data/') +p.plot(r=1,proj='Mollweide') +``` + +produces the info mentioned above first and then the following plot of Earth's magnetic field using a Mollweide projection + + + +while + +```python +from planetmagfields import * +p = planet(name='jupiter',datDir='planetmagfields/data/') +p.plot(r=0.85,proj='Mollweide') +``` +produces the following info about Jupiter and then plot that follows + +``` +Planet: Jupiter +l_max = 10 +Dipole tilt (degrees) = 10.307870 +``` + + + +This can be compared with Fig. 1 g from [Moore et al. 2018](https://doi.org/10.1038/s41586-018-0468-5) + +## `planet.spec()` + +This function computes the Lowes spectrum of a planet at a given radius. It adds an array `p.emag_spec` which contains the energy at different spherical harmonic degrees and two variables `p.dipolarity` and `p.dip_tot` which provide the fraction of energies in the axial dipole and the total dipole with respect to the total energy at all degrees. Usage example: + +```python +from planetmagfields import * +p = planet(name='jupiter') +p.spec() +``` + +will provide variables + +```python +In [8]: p.dipolarity +Out[8]: 0.7472047129875864 + +In [9]: p.dip_tot +Out[9]: 0.7719205112704368 + +In [10]: p.emag_spec +Out[10]: +array([0.00000000e+00, 3.47735422e+11, 2.36340423e+10, 2.12851278e+10, + 1.75661779e+10, 1.92219842e+10, 9.91200756e+09, 3.34535475e+09, + 3.95317971e+09, 2.59333412e+09, 1.23423769e+09]) +``` + +and will produce Jupiter's surface spectrum: + + + +The plotting can be suppressed setting the logical `p.spec(iplot=False)`. + +## `planet.writeVtsFile()` + +This function writes a vts file that can be used to produce 3D visualizations of field lines with Paraview/VisIt. Usage: + +```python +p.writeVtsFile(potExtra=True, ratio_out=2, nrout=32) +``` +where, + + - `potExtra` : bool, whether to use potential extrapolation + - `ratio_out`: float, radius till which the field would be extrapolated in terms of the surface radius + - `nrout`: radial resolution for extrapolation + +Example of a 3D image produced using Paraview for Neptune's field, extrapolated till 5 times the surface radius. + + + +# Field filtering using `planet.plot_filt` + +The `planet` class also provides a function for producing a filtered view of the radial magnetic field using the function `plot_filt`. +This function can take in either an arbitrary array of spherical harmonic degrees and orders or cut-off values. This is illustrated +below with examples, assuming the user is in the repository directory. + +## Saturn's small-scale magnetic field at a depth of 0.75 planetary radius ( degree > 3 ) + +```python +from planetmagfields import * +p = planet(name='saturn') +p.plot_filt(r=0.75,lCutMin=4,proj='Mollweide') +``` + + + +Compare this with Fig. 20 B from [Cao et al. 2020](https://doi.org/10.1016/j.icarus.2019.113541) + +## Jupiter's surface field restricted to degrees 1,2,3 and order 3 + +```python +from planetmagfields import * +p = planet(name='jupiter') +p.plot_filt(r=1,larr=[1,2,3],marr=[3],proj='Mollweide') +``` + + + +## Earth's smaller scale surface field for degree > 4 and order > 3 + +```python +from planetmagfields import * +p = planet(name='earth') +p.plot_filt(r=1,lCutMin=5,mmin=4,proj='Mollweide') +``` + + + +# Potential extrapolation + +The `potextra` module provides a method for potential extrapolation of a planet's magnetic field. +This uses the [SHTns](https://bitbucket.org/nschaeff/shtns) library for spherical harmonic transforms. +Usage example: + +```python +import numpy as np +from planetmagfields import * +p = planet('saturn') +ratio_out = 5 # Ratio (> 1) in terms of surface radius to which to extrapolate +nrout = 32 # Number of grid points in radius between 1 and ratio_out +rout = np.linspace(1,ratio_out,nrout) +brout, btout, bpout = potextra.extrapot(p.lmax,1.,p.Br,rout) +``` + +# Quickplot using the `magField.py` script: + +``` +$ ./magField.py +``` + +This will plot the radial magnetic field of a planet (any of the names from the list +below, case insensitive) at a radius given in terms of the surface radius. The default +is the surface field. For example, + +``` +$ ./magField.py earth Mollweide +``` + +displays the same information as above about Earth's field and produces the surface field of Earth while + +``` +$ ./magField.py jupiter 0.85 Mollweide +``` + +produces the same plot of Jupiter's field as shown before. + +``` +$ ./magField.py all +``` + +would produce a table of information about dipole tilt for each planet and magnetic field maps of all different planets at the given radius in a single figure. + +For example: + +``` +$ ./magField all 0.9 Mollweide +``` + +would give + +``` +|=========|======|=======| +|Planet | Theta| Phi | +|=========|======|=======| +|Mercury | 0.0 | 0.0 | +|Earth | -9.4 | -72.7 | +|Jupiter | 10.3 | -16.6 | +|Saturn | 0.0 | 0.0 | +|Uranus | 58.6 | -53.6 | +|Neptune | 46.9 | -72.0 | +|Ganymede | -4.2 | 25.5 | +|---------|------|-------| +``` + +followed by the following plot + +![All fields r=0.9](planetmagfields/images/magField_all_09.png) + +# Spherical harmonic normalization and Condon-Shortley phase + +All the Gauss coefficients in the collected data are Schmidt semi-normalized. +Only the data for Earth uses a Condon-Shortley phase, the others do not. + +# Projections + +By default, the plot function tries to use the Mollweide projection. However, using the power of the [cartopy](https://scitools.org.uk/cartopy/docs/latest/) library, any projection from [this list](https://scitools.org.uk/cartopy/docs/latest/crs/projections.html) is supported. In the absence of the cartopy library, the 2D plots fall back to the internally written Hammer projection. Examples of Jupiter's radial magnetic field at r=0.85 with different projections are shown below: + +```python +In [1]: from planetmagfields import * + +In [2]: p = planet(name='jupiter') +Planet: Jupiter +l_max = 10 +Dipole tilt (degrees) = 10.307870 + +In [3]: p.plot(r=0.85) + +In [4]: projlist=['Mercator','Robinson','Stereographic','AzimuthalEquidistant'] + +In [5]: for k,proj in enumerate(projlist): + ...: p.plot(r=0.85,proj=proj) + ...: savefig('jup_r0_85'+proj+'.png',dpi=150,bbox_inches='tight') + ...: close() +``` + +`In[3]` produces the plot of Jupiter's field already shown above. `In[5]` produces the following figures with the Mercator, Robinson, Stereographic and azimuthal equidistant projections, respectively. + + + + + +This also works with the `magField.py` script for quick plotting. Examples: + +```bash +./magField.py earth 0.9 Robinson +``` + +or even with plots of all planets together + +```bash +./magField.py all 0.9 Robinson +``` + +Note that I have chosen to keep the projection information out of the plot titles to prevent too much text. + +❗ | The Orthographic projection often does not create correct plots, be cautious while using it +:---: | :--- + +# Data sources + +Mercury : [Anderson et al. 2012](https://doi.org/10.1029/2012JE004159) + +Earth : [IGRF 13](https://www.ngdc.noaa.gov/IAGA/vmod/igrf.html), [Alken et al. 2021](https://doi.org/10.1186/s40623-020-01288-x) + +Jupiter : [JRM09, Connerney et al. 2018](https://doi.org/10.1002/2018GL077312) + +Saturn : [Cassini 11+, Cao et al. 2020](https://doi.org/10.1016/j.icarus.2019.113541) + +Uranus : [Connerny et al. 1987](https://doi.org/10.1029/JA092iA13p15329) + +Neptune : [Connerny et al. 1991](https://doi.org/10.1029/91JA01165) + +Ganymede: [Kivelson et al. 2002](https://doi.org/10.1006/icar.2002.6834) diff --git a/src/astroedu/planetmagfields/__init__.py b/src/astroedu/planetmagfields/__init__.py new file mode 100644 index 0000000..8b8c597 --- /dev/null +++ b/src/astroedu/planetmagfields/__init__.py @@ -0,0 +1,5 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +from .planet import planet +from .potextra import extrapot \ No newline at end of file diff --git a/src/astroedu/planetmagfields/data/IGRF13.dat b/src/astroedu/planetmagfields/data/IGRF13.dat new file mode 100644 index 0000000..db67ed5 --- /dev/null +++ b/src/astroedu/planetmagfields/data/IGRF13.dat @@ -0,0 +1,195 @@ +g 1 0 -31543 -31464 -31354 -31212 -31060 -30926 -30805 -30715 -30654 -30594 -30554 -30500 -30421 -30334 -30220 -30100 -29992 -29873 -29775 -29692 -29619.4 -29554.63 -29496.57 -29441.46 -29404.8 5.7 +g 1 1 -2298 -2298 -2297 -2306 -2317 -2318 -2316 -2306 -2292 -2285 -2250 -2215 -2169 -2119 -2068 -2013 -1956 -1905 -1848 -1784 -1728.2 -1669.05 -1586.42 -1501.77 -1450.9 7.4 +h 1 1 5922 5909 5898 5875 5845 5817 5808 5812 5821 5810 5815 5820 5791 5776 5737 5675 5604 5500 5406 5306 5186.1 5077.99 4944.26 4795.99 4652.5 -25.9 +g 2 0 -677 -728 -769 -802 -839 -893 -951 -1018 -1106 -1244 -1341 -1440 -1555 -1662 -1781 -1902 -1997 -2072 -2131 -2200 -2267.7 -2337.24 -2396.06 -2445.88 -2499.6 -11.0 +g 2 1 2905 2928 2948 2956 2959 2969 2980 2984 2981 2990 2998 3003 3002 2997 3000 3010 3027 3044 3059 3070 3068.4 3047.69 3026.34 3012.20 2982.0 -7.0 +h 2 1 -1061 -1086 -1128 -1191 -1259 -1334 -1424 -1520 -1614 -1702 -1810 -1898 -1967 -2016 -2047 -2067 -2129 -2197 -2279 -2366 -2481.6 -2594.50 -2708.54 -2845.41 -2991.6 -30.2 +g 2 2 924 1041 1176 1309 1407 1471 1517 1550 1566 1578 1576 1581 1590 1594 1611 1632 1663 1687 1686 1681 1670.9 1657.76 1668.17 1676.35 1677.0 -2.1 +h 2 2 1121 1065 1000 917 823 728 644 586 528 477 381 291 206 114 25 -68 -200 -306 -373 -413 -458.0 -515.43 -575.73 -642.17 -734.6 -22.4 +g 3 0 1022 1037 1058 1084 1111 1140 1172 1206 1240 1282 1297 1302 1302 1297 1287 1276 1281 1296 1314 1335 1339.6 1336.30 1339.85 1350.33 1363.2 2.2 +g 3 1 -1469 -1494 -1524 -1559 -1600 -1645 -1692 -1740 -1790 -1834 -1889 -1944 -1992 -2038 -2091 -2144 -2180 -2208 -2239 -2267 -2288.0 -2305.83 -2326.54 -2352.26 -2381.2 -5.9 +h 3 1 -330 -357 -389 -421 -445 -462 -480 -494 -499 -499 -476 -462 -414 -404 -366 -333 -336 -310 -284 -262 -227.6 -198.86 -160.40 -115.29 -82.1 6.0 +g 3 2 1256 1239 1223 1212 1205 1202 1205 1215 1232 1255 1274 1288 1289 1292 1278 1260 1251 1247 1248 1249 1252.1 1246.39 1232.10 1225.85 1236.2 3.1 +h 3 2 3 34 62 84 103 119 133 146 163 186 206 216 224 240 251 262 271 284 293 302 293.4 269.72 251.75 245.04 241.9 -1.1 +g 3 3 572 635 705 778 839 881 907 918 916 913 896 882 878 856 838 830 833 829 802 759 714.5 672.51 633.73 581.69 525.7 -12.0 +h 3 3 523 480 425 360 293 229 166 101 43 -11 -46 -83 -130 -165 -196 -223 -252 -297 -352 -427 -491.1 -524.72 -537.03 -538.70 -543.4 0.5 +g 4 0 876 880 884 887 889 891 896 903 914 944 954 958 957 957 952 946 938 936 939 940 932.3 920.55 912.66 907.42 903.0 -1.2 +g 4 1 628 643 660 678 695 711 727 744 762 776 792 796 800 804 800 791 782 780 780 780 786.8 797.96 808.97 813.68 809.5 -1.6 +h 4 1 195 203 211 218 220 216 205 188 169 144 136 133 135 148 167 191 212 232 247 262 272.6 282.07 286.48 283.54 281.9 -0.1 +g 4 2 660 653 644 631 616 601 584 565 550 544 528 510 504 479 461 438 398 361 325 290 250.0 210.65 166.58 120.49 86.3 -5.9 +h 4 2 -69 -77 -90 -109 -134 -163 -195 -226 -252 -276 -278 -274 -278 -269 -266 -265 -257 -249 -240 -236 -231.9 -225.23 -211.03 -188.43 -158.4 6.5 +g 4 3 -361 -380 -400 -416 -424 -426 -422 -415 -405 -421 -408 -397 -394 -390 -395 -405 -419 -424 -423 -418 -403.0 -379.86 -356.83 -334.85 -309.4 5.2 +h 4 3 -210 -201 -189 -173 -153 -130 -109 -90 -72 -55 -37 -23 3 13 26 39 53 69 84 97 119.8 145.15 164.46 180.95 199.7 3.6 +g 4 4 134 146 160 178 199 217 234 249 265 304 303 290 269 252 234 216 199 170 141 122 111.3 100.00 89.40 70.38 48.0 -5.1 +h 4 4 -75 -65 -55 -51 -57 -70 -90 -114 -141 -178 -210 -230 -255 -269 -279 -288 -297 -297 -299 -306 -303.8 -305.36 -309.72 -329.23 -349.7 -5.0 +g 5 0 -184 -192 -201 -211 -221 -230 -237 -241 -241 -253 -240 -229 -222 -219 -216 -218 -218 -214 -214 -214 -218.8 -227.00 -230.87 -232.91 -234.3 -0.3 +g 5 1 328 328 327 327 326 326 327 329 334 346 349 360 362 358 359 356 357 355 353 352 351.4 354.41 357.29 360.14 363.2 0.5 +h 5 1 -210 -193 -172 -148 -122 -96 -72 -51 -33 -12 3 15 16 19 26 31 46 47 46 46 43.8 42.72 44.58 46.98 47.7 0.0 +g 5 2 264 259 253 245 236 226 218 211 208 194 211 230 242 254 262 264 261 253 245 235 222.3 208.95 200.26 192.35 187.8 -0.6 +h 5 2 53 56 57 58 58 58 60 64 71 95 103 110 125 128 139 148 150 150 154 165 171.9 180.25 189.01 196.98 208.3 2.5 +g 5 3 5 -1 -9 -16 -23 -28 -32 -33 -33 -20 -20 -23 -26 -31 -42 -59 -74 -93 -109 -118 -130.4 -136.54 -141.05 -140.94 -140.7 0.2 +h 5 3 -33 -32 -33 -34 -38 -44 -53 -64 -75 -67 -87 -98 -117 -126 -139 -152 -151 -154 -153 -143 -133.1 -123.45 -118.06 -119.14 -121.2 -0.6 +g 5 4 -86 -93 -102 -111 -119 -125 -131 -136 -141 -142 -147 -152 -156 -157 -160 -159 -162 -164 -165 -166 -168.6 -168.05 -163.17 -157.40 -151.2 1.3 +h 5 4 -124 -125 -126 -126 -125 -122 -118 -115 -113 -119 -122 -121 -114 -97 -91 -83 -78 -75 -69 -55 -39.3 -19.57 -0.01 15.98 32.3 3.0 +g 5 5 -16 -26 -38 -51 -62 -69 -74 -76 -76 -82 -76 -69 -63 -62 -56 -49 -48 -46 -36 -17 -12.9 -13.55 -8.03 4.30 13.5 0.9 +h 5 5 3 11 21 32 43 51 58 64 69 82 80 78 81 81 83 88 92 95 97 107 106.3 103.85 101.04 100.12 98.9 0.3 +g 6 0 63 62 62 61 61 61 60 59 57 59 54 47 46 45 43 45 48 53 61 68 72.3 73.60 72.78 69.55 66.0 -0.5 +g 6 1 61 60 58 57 55 54 53 53 54 57 57 57 58 61 64 66 66 65 65 67 68.2 69.56 68.69 67.57 65.5 -0.3 +h 6 1 -9 -7 -5 -2 0 3 4 4 4 6 -1 -9 -10 -11 -12 -13 -15 -16 -16 -17 -17.4 -20.33 -20.90 -20.61 -19.1 0.0 +g 6 2 -11 -11 -11 -10 -10 -9 -9 -8 -7 6 4 3 1 8 15 28 42 51 59 68 74.2 76.74 75.92 72.79 72.9 0.4 +h 6 2 83 86 89 93 96 99 102 104 105 100 99 96 99 100 100 99 93 88 82 72 63.7 54.75 44.18 33.30 25.1 -1.6 +g 6 3 -217 -221 -224 -228 -233 -238 -242 -246 -249 -246 -247 -247 -237 -228 -212 -198 -192 -185 -178 -170 -160.9 -151.34 -141.40 -129.85 -121.5 1.3 +h 6 3 2 4 5 8 11 14 19 25 33 16 33 48 60 68 72 75 71 69 69 67 65.1 63.63 61.54 58.74 52.8 -1.3 +g 6 4 -58 -57 -54 -51 -46 -40 -32 -25 -18 -25 -16 -8 -1 4 2 1 4 4 3 -1 -5.9 -14.58 -22.83 -28.93 -36.2 -1.4 +h 6 4 -35 -32 -29 -26 -22 -18 -16 -15 -15 -9 -12 -16 -20 -32 -37 -41 -43 -48 -52 -58 -61.2 -63.53 -66.26 -66.64 -64.5 0.8 +g 6 5 59 57 54 49 44 39 32 25 18 21 12 7 -2 1 3 6 14 16 18 19 16.9 14.58 13.10 13.14 13.5 0.0 +h 6 5 36 32 28 23 18 13 8 4 0 -16 -12 -12 -11 -8 -6 -4 -2 -1 1 1 0.7 0.24 3.02 7.35 8.9 0.0 +g 6 6 -90 -92 -95 -98 -101 -103 -104 -106 -107 -104 -105 -107 -113 -111 -112 -111 -108 -102 -96 -93 -90.4 -86.36 -78.09 -70.85 -64.7 0.9 +h 6 6 -69 -67 -65 -62 -57 -52 -46 -40 -33 -39 -30 -24 -17 -7 1 11 17 21 24 36 43.8 50.94 55.40 62.41 68.1 1.0 +g 7 0 70 70 71 72 73 73 74 74 74 70 65 65 67 75 72 71 72 74 77 77 79.0 79.88 80.44 81.29 80.6 -0.1 +g 7 1 -55 -54 -54 -54 -54 -54 -54 -53 -53 -40 -55 -56 -56 -57 -57 -56 -59 -62 -64 -72 -74.0 -74.46 -75.00 -75.99 -76.7 -0.2 +h 7 1 -45 -46 -47 -48 -49 -50 -51 -52 -52 -45 -35 -50 -55 -61 -70 -77 -82 -83 -80 -69 -64.6 -61.14 -57.80 -54.27 -51.5 0.6 +g 7 2 0 0 1 2 2 3 4 4 4 0 2 2 5 4 1 1 2 3 2 1 0.0 -1.65 -4.55 -6.79 -8.2 0.0 +h 7 2 -13 -14 -14 -14 -14 -14 -15 -17 -18 -18 -17 -24 -28 -27 -27 -26 -27 -27 -26 -25 -24.2 -22.57 -21.20 -19.53 -16.9 0.6 +g 7 3 34 33 32 31 29 27 25 23 20 0 1 10 15 13 14 16 21 24 26 28 33.3 38.73 45.24 51.82 56.5 0.7 +h 7 3 -10 -11 -12 -12 -13 -14 -14 -14 -14 2 0 -4 -6 -2 -4 -5 -5 -2 0 4 6.2 6.82 6.54 5.59 2.2 -0.8 +g 7 4 -41 -41 -40 -38 -37 -35 -34 -33 -31 -29 -40 -32 -32 -26 -22 -14 -12 -6 -1 5 9.1 12.30 14.00 15.07 15.8 0.1 +h 7 4 -1 0 1 2 4 5 6 7 7 6 10 8 7 6 8 10 16 20 21 24 24.0 25.35 24.96 24.45 23.5 -0.2 +g 7 5 -21 -20 -19 -18 -16 -14 -12 -11 -9 -10 -7 -11 -7 -6 -2 0 1 4 5 4 6.9 9.37 10.46 9.32 6.4 -0.5 +h 7 5 28 28 28 28 28 29 29 29 29 28 36 28 23 26 23 22 18 17 17 17 14.8 10.93 7.03 3.27 -2.2 -1.1 +g 7 6 18 18 18 19 19 19 18 18 17 15 5 9 17 13 13 12 11 10 9 8 7.3 5.42 1.64 -2.88 -7.2 -0.8 +h 7 6 -12 -12 -13 -15 -16 -17 -18 -19 -20 -17 -18 -20 -18 -23 -23 -23 -23 -23 -23 -24 -25.4 -26.32 -27.61 -27.50 -27.2 0.1 +g 7 7 6 6 6 6 6 6 6 6 5 29 19 18 8 1 -2 -5 -2 0 0 -2 -1.2 1.94 4.92 6.61 9.8 0.8 +h 7 7 -22 -22 -22 -22 -22 -21 -20 -19 -19 -22 -16 -18 -17 -12 -11 -12 -10 -7 -4 -6 -5.8 -4.64 -3.28 -2.32 -1.8 0.3 +g 8 0 11 11 11 11 11 11 11 11 11 13 22 11 15 13 14 14 18 21 23 25 24.4 24.80 24.41 23.98 23.7 0.0 +g 8 1 8 8 8 8 7 7 7 7 7 7 15 9 6 5 6 6 6 6 5 6 6.6 7.62 8.21 8.89 9.7 0.1 +h 8 1 8 8 8 8 8 8 8 8 8 12 5 10 11 7 7 6 7 8 10 11 11.9 11.20 10.84 10.04 8.4 -0.2 +g 8 2 -4 -4 -4 -4 -3 -3 -3 -3 -3 -8 -4 -6 -4 -4 -2 -1 0 0 -1 -6 -9.2 -11.73 -14.50 -16.78 -17.6 -0.1 +h 8 2 -14 -15 -15 -15 -15 -15 -15 -15 -14 -21 -22 -15 -14 -12 -15 -16 -18 -19 -19 -21 -21.5 -20.88 -20.03 -18.26 -15.3 0.6 +g 8 3 -9 -9 -9 -9 -9 -9 -9 -9 -10 -5 -1 -14 -11 -14 -13 -12 -11 -11 -10 -9 -7.9 -6.88 -5.59 -3.16 -0.5 0.4 +h 8 3 7 7 6 6 6 6 5 5 5 -12 0 5 7 9 6 4 4 5 6 8 8.5 9.83 11.83 13.18 12.8 -0.2 +g 8 4 1 1 1 2 2 2 2 1 1 9 11 6 2 0 -3 -8 -7 -9 -12 -14 -16.6 -18.11 -19.34 -20.56 -21.1 -0.1 +h 8 4 -13 -13 -13 -13 -14 -14 -14 -15 -15 -7 -21 -23 -18 -16 -17 -19 -22 -23 -22 -23 -21.5 -19.71 -17.41 -14.60 -11.7 0.5 +g 8 5 2 2 2 3 4 4 5 6 6 7 15 10 10 8 5 4 4 4 3 9 9.1 10.17 11.61 13.33 15.3 0.4 +h 8 5 5 5 5 5 5 5 5 5 5 2 -8 3 4 4 6 6 9 11 12 15 15.5 16.22 16.71 16.16 14.9 -0.3 +g 8 6 -9 -8 -8 -8 -7 -7 -6 -6 -5 -10 -13 -7 -5 -1 0 0 3 4 4 6 7.0 9.36 10.85 11.76 13.7 0.3 +h 8 6 16 16 16 16 17 17 18 18 19 18 17 23 23 24 21 18 16 14 12 11 8.9 7.61 6.96 5.69 3.6 -0.4 +g 8 7 5 5 5 6 6 7 8 8 9 7 5 6 10 11 11 10 6 4 2 -5 -7.9 -11.25 -14.05 -15.98 -16.5 -0.1 +h 8 7 -5 -5 -5 -5 -5 -5 -5 -5 -5 3 -4 -4 1 -3 -6 -10 -13 -15 -16 -16 -14.9 -12.76 -10.74 -9.10 -6.9 0.5 +g 8 8 8 8 8 8 8 8 8 7 7 2 -1 9 8 4 3 1 -1 -4 -6 -7 -7.0 -4.87 -3.54 -2.02 -0.3 0.4 +h 8 8 -18 -18 -18 -18 -19 -19 -19 -19 -19 -11 -17 -13 -20 -17 -16 -17 -15 -11 -10 -4 -2.1 -0.06 1.64 2.26 2.8 0.0 +g 9 0 8 8 8 8 8 8 8 8 8 5 3 4 4 8 8 7 5 5 4 4 5.0 5.58 5.50 5.33 5.0 0.0 +g 9 1 10 10 10 10 10 10 10 10 10 -21 -7 9 6 10 10 10 10 10 9 9 9.4 9.76 9.45 8.83 8.4 0.0 +h 9 1 -20 -20 -20 -20 -20 -20 -20 -20 -21 -27 -24 -11 -18 -22 -21 -21 -21 -21 -20 -20 -19.7 -20.11 -20.54 -21.77 -23.4 0.0 +g 9 2 1 1 1 1 1 1 1 1 1 1 -1 -4 0 2 2 2 1 1 1 3 3.0 3.58 3.45 3.02 2.9 0.0 +h 9 2 14 14 14 14 14 14 14 15 15 17 19 12 12 15 16 16 16 15 15 15 13.4 12.69 11.51 10.76 11.0 0.0 +g 9 3 -11 -11 -11 -11 -11 -11 -12 -12 -12 -11 -25 -5 -9 -13 -12 -12 -12 -12 -12 -10 -8.4 -6.94 -5.27 -3.22 -1.5 0.0 +h 9 3 5 5 5 5 5 5 5 5 5 29 12 7 2 7 6 7 9 9 11 12 12.5 12.67 12.75 11.74 9.8 0.0 +g 9 4 12 12 12 12 12 12 12 11 11 3 10 2 1 10 10 10 9 9 9 8 6.3 5.01 3.13 0.67 -1.1 0.0 +h 9 4 -3 -3 -3 -3 -3 -3 -3 -3 -3 -9 2 6 0 -4 -4 -4 -5 -6 -7 -6 -6.2 -6.72 -7.14 -6.74 -5.1 0.0 +g 9 5 1 1 1 1 1 1 1 1 1 16 5 4 4 -1 -1 -1 -3 -3 -4 -8 -8.9 -10.76 -12.38 -13.20 -13.2 0.0 +h 9 5 -2 -2 -2 -2 -2 -2 -2 -3 -3 4 2 -2 -3 -5 -5 -5 -6 -6 -7 -8 -8.4 -8.16 -7.42 -6.88 -6.3 0.0 +g 9 6 -2 -2 -2 -2 -2 -2 -2 -2 -2 -3 -5 1 -1 -1 0 -1 -1 -1 -2 -1 -1.5 -1.25 -0.76 -0.10 1.1 0.0 +h 9 6 8 8 8 8 9 9 9 9 9 9 8 10 9 10 10 10 9 9 9 8 8.4 8.10 7.97 7.79 7.8 0.0 +g 9 7 2 2 2 2 2 2 3 3 3 -4 -2 2 -2 5 3 4 7 7 7 10 9.3 8.76 8.43 8.68 8.8 0.0 +h 9 7 10 10 10 10 10 10 10 11 11 6 8 7 8 10 11 11 10 9 8 5 3.8 2.92 2.14 1.04 0.4 0.0 +g 9 8 -1 0 0 0 0 0 0 0 1 -3 3 2 3 1 1 1 2 1 1 -2 -4.3 -6.66 -8.42 -9.06 -9.3 0.0 +h 9 8 -2 -2 -2 -2 -2 -2 -2 -2 -2 1 -11 -6 0 -4 -2 -3 -6 -7 -7 -8 -8.2 -7.73 -6.08 -3.89 -1.4 0.0 +g 9 9 -1 -1 -1 -1 -1 -1 -2 -2 -2 -4 8 5 -1 -2 -1 -2 -5 -5 -6 -8 -8.2 -9.22 -10.08 -10.54 -11.9 0.0 +h 9 9 2 2 2 2 2 2 2 2 2 8 -7 5 5 1 1 1 2 2 2 3 4.8 6.01 7.01 8.44 9.6 0.0 +g 10 0 -3 -3 -3 -3 -3 -3 -3 -3 -3 -3 -8 -3 1 -2 -3 -3 -4 -4 -3 -3 -2.6 -2.17 -1.94 -2.01 -1.9 0.0 +g 10 1 -4 -4 -4 -4 -4 -4 -4 -4 -4 11 4 -5 -3 -3 -3 -3 -4 -4 -4 -6 -6.0 -6.12 -6.24 -6.26 -6.2 0.0 +h 10 1 2 2 2 2 2 2 2 2 2 5 13 -4 4 2 1 1 1 1 2 1 1.7 2.19 2.73 3.28 3.4 0.0 +g 10 2 2 2 2 2 2 2 2 2 2 1 -1 -1 4 2 2 2 2 3 2 2 1.7 1.42 0.89 0.17 -0.1 0.0 +h 10 2 1 1 1 1 1 1 1 1 1 1 -2 0 1 1 1 1 0 0 1 0 0.0 0.10 -0.10 -0.40 -0.2 0.0 +g 10 3 -5 -5 -5 -5 -5 -5 -5 -5 -5 2 13 2 0 -5 -5 -5 -5 -5 -5 -4 -3.1 -2.35 -1.07 0.55 1.7 0.0 +h 10 3 2 2 2 2 2 2 2 2 2 -20 -10 -8 0 2 3 3 3 3 3 4 4.0 4.46 4.71 4.55 3.6 0.0 +g 10 4 -2 -2 -2 -2 -2 -2 -2 -2 -2 -5 -4 -3 -1 -2 -1 -2 -2 -2 -2 -1 -0.5 -0.15 -0.16 -0.55 -0.9 0.0 +h 10 4 6 6 6 6 6 6 6 6 6 -1 2 -2 2 6 4 4 6 6 6 5 4.9 4.76 4.44 4.40 4.8 0.0 +g 10 5 6 6 6 6 6 6 6 6 6 -1 4 7 4 4 6 5 5 5 4 4 3.7 3.06 2.45 1.70 0.7 0.0 +h 10 5 -4 -4 -4 -4 -4 -4 -4 -4 -4 -6 -3 -4 -5 -4 -4 -4 -4 -4 -4 -5 -5.9 -6.58 -7.22 -7.92 -8.6 0.0 +g 10 6 4 4 4 4 4 4 4 4 4 8 12 4 6 4 4 4 3 3 3 2 1.0 0.29 -0.33 -0.67 -0.9 0.0 +h 10 6 0 0 0 0 0 0 0 0 0 6 6 1 1 0 0 -1 0 0 0 -1 -1.2 -1.01 -0.96 -0.61 -0.1 0.0 +g 10 7 0 0 0 0 0 0 0 0 0 -1 3 -2 1 0 1 1 1 1 1 2 2.0 2.06 2.13 2.13 1.9 0.0 +h 10 7 -2 -2 -2 -2 -2 -2 -2 -1 -1 -4 -3 -3 -1 -2 -1 -1 -1 -1 -2 -2 -2.9 -3.47 -3.95 -4.16 -4.3 0.0 +g 10 8 2 2 2 1 1 1 1 2 2 -3 2 6 -1 2 0 0 2 2 3 5 4.2 3.77 3.09 2.33 1.4 0.0 +h 10 8 4 4 4 4 4 4 4 4 4 -2 6 7 6 3 3 3 4 4 3 1 0.2 -0.86 -1.99 -2.85 -3.4 0.0 +g 10 9 2 2 2 2 3 3 3 3 3 5 10 -2 2 2 3 3 3 3 3 1 0.3 -0.21 -1.03 -1.80 -2.4 0.0 +h 10 9 0 0 0 0 0 0 0 0 0 0 11 -1 0 0 1 1 0 0 -1 -2 -2.2 -2.31 -1.97 -1.12 -0.1 0.0 +g 10 10 0 0 0 0 0 0 0 0 0 -2 3 0 0 0 -1 -1 0 0 0 0 -1.1 -2.09 -2.80 -3.59 -3.8 0.0 +h 10 10 -6 -6 -6 -6 -6 -6 -6 -6 -6 -2 8 -3 -7 -6 -4 -5 -6 -6 -6 -7 -7.4 -7.93 -8.31 -8.72 -8.8 0.0 +g 11 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.7 2.95 3.05 3.00 3.0 0.0 +g 11 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1.7 -1.60 -1.48 -1.40 -1.4 0.0 +h 11 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.1 0.26 0.13 0.00 0.0 0.0 +g 11 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1.9 -1.88 -2.03 -2.30 -2.5 0.0 +h 11 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.3 1.44 1.67 2.11 2.5 0.0 +g 11 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.5 1.44 1.65 2.08 2.3 0.0 +h 11 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.9 -0.77 -0.66 -0.60 -0.6 0.0 +g 11 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.1 -0.31 -0.51 -0.79 -0.9 0.0 +h 11 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -2.6 -2.27 -1.76 -1.05 -0.4 0.0 +g 11 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.1 0.29 0.54 0.58 0.3 0.0 +h 11 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.9 0.90 0.85 0.76 0.6 0.0 +g 11 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.7 -0.79 -0.79 -0.70 -0.7 0.0 +h 11 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.7 -0.58 -0.39 -0.20 -0.2 0.0 +g 11 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.7 0.53 0.37 0.14 -0.1 0.0 +h 11 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -2.8 -2.69 -2.51 -2.12 -1.7 0.0 +g 11 8 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.7 1.80 1.79 1.70 1.4 0.0 +h 11 8 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.9 -1.08 -1.27 -1.44 -1.6 0.0 +g 11 9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.1 0.16 0.12 -0.22 -0.6 0.0 +h 11 9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1.2 -1.58 -2.11 -2.57 -3.0 0.0 +g 11 10 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.2 0.96 0.75 0.44 0.2 0.0 +h 11 10 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1.9 -1.90 -1.94 -2.01 -2.0 0.0 +g 11 11 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4.0 3.99 3.75 3.49 3.1 0.0 +h 11 11 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.9 -1.39 -1.86 -2.34 -2.6 0.0 +g 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -2.2 -2.15 -2.12 -2.09 -2.0 0.0 +g 12 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.3 -0.29 -0.21 -0.16 -0.1 0.0 +h 12 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.4 -0.55 -0.87 -1.08 -1.2 0.0 +g 12 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.2 0.21 0.30 0.46 0.5 0.0 +h 12 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.3 0.23 0.27 0.37 0.5 0.0 +g 12 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.9 0.89 1.04 1.23 1.3 0.0 +h 12 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.5 2.38 2.13 1.75 1.4 0.0 +g 12 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.2 -0.38 -0.63 -0.89 -1.2 0.0 +h 12 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -2.6 -2.63 -2.49 -2.19 -1.8 0.0 +g 12 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.9 0.96 0.95 0.85 0.7 0.0 +h 12 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.7 0.61 0.49 0.27 0.1 0.0 +g 12 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.5 -0.30 -0.11 0.10 0.3 0.0 +h 12 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.3 0.40 0.59 0.72 0.8 0.0 +g 12 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.3 0.46 0.52 0.54 0.5 0.0 +h 12 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0 0.01 0.00 -0.09 -0.2 0.0 +g 12 8 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.3 -0.35 -0.39 -0.37 -0.3 0.0 +h 12 8 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0 0.02 0.13 0.29 0.6 0.0 +g 12 9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.4 -0.36 -0.37 -0.43 -0.5 0.0 +h 12 9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.3 0.28 0.27 0.23 0.2 0.0 +g 12 10 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.1 0.08 0.21 0.22 0.1 0.0 +h 12 10 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.9 -0.87 -0.86 -0.89 -0.9 0.0 +g 12 11 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.2 -0.49 -0.77 -0.94 -1.1 0.0 +h 12 11 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.4 -0.34 -0.23 -0.16 0.0 0.0 +g 12 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.4 -0.08 0.04 -0.03 -0.3 0.0 +h 12 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.8 0.88 0.87 0.72 0.5 0.0 +g 13 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.2 -0.16 -0.09 -0.02 0.1 0.0 +g 13 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.9 -0.88 -0.89 -0.92 -0.9 0.0 +h 13 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.9 -0.76 -0.87 -0.88 -0.9 0.0 +g 13 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.3 0.30 0.31 0.42 0.5 0.0 +h 13 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.2 0.33 0.30 0.49 0.6 0.0 +g 13 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.1 0.28 0.42 0.63 0.7 0.0 +h 13 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.8 1.72 1.66 1.56 1.4 0.0 +g 13 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.4 -0.43 -0.45 -0.42 -0.3 0.0 +h 13 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.4 -0.54 -0.59 -0.50 -0.4 0.0 +g 13 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.3 1.18 1.08 0.96 0.8 0.0 +h 13 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1.0 -1.07 -1.14 -1.24 -1.3 0.0 +g 13 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.4 -0.37 -0.31 -0.19 0.0 0.0 +h 13 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.1 -0.04 -0.07 -0.10 -0.1 0.0 +g 13 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.7 0.75 0.78 0.81 0.8 0.0 +h 13 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.7 0.63 0.54 0.42 0.3 0.0 +g 13 8 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.4 -0.26 -0.18 -0.13 0.0 0.0 +h 13 8 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.3 0.21 0.10 -0.04 -0.1 0.0 +g 13 9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.3 0.35 0.38 0.38 0.4 0.0 +h 13 9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.6 0.53 0.49 0.48 0.5 0.0 +g 13 10 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.1 -0.05 0.02 0.08 0.1 0.0 +h 13 10 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.3 0.38 0.44 0.48 0.5 0.0 +g 13 11 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.4 0.41 0.42 0.46 0.5 0.0 +h 13 11 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.2 -0.22 -0.25 -0.30 -0.4 0.0 +g 13 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0 -0.10 -0.26 -0.35 -0.5 0.0 +h 13 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.5 -0.57 -0.53 -0.43 -0.4 0.0 +g 13 13 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.1 -0.18 -0.26 -0.36 -0.4 0.0 +h 13 13 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.9 -0.82 -0.79 -0.71 -0.6 0.0 diff --git a/src/astroedu/planetmagfields/data/ganymede.dat b/src/astroedu/planetmagfields/data/ganymede.dat new file mode 100644 index 0000000..883e679 --- /dev/null +++ b/src/astroedu/planetmagfields/data/ganymede.dat @@ -0,0 +1,10 @@ +g 1 0 -711.0 +g 1 1 46.8 +h 1 0 0.0 +h 1 1 22.3 +g 2 0 0.9 +g 2 1 27.0 +g 2 2 -0.4 +h 2 0 0.0 +h 2 1 1.8 +h 2 2 -11.0 diff --git a/src/astroedu/planetmagfields/data/jupiter.dat b/src/astroedu/planetmagfields/data/jupiter.dat new file mode 100644 index 0000000..ec8c732 --- /dev/null +++ b/src/astroedu/planetmagfields/data/jupiter.dat @@ -0,0 +1,120 @@ + 1 410244.7 1.0 1 0 + 2 -71498.3 1.0 1 1 + 3 21330.5 1.0 1 1 + 4 11670.4 1.0 2 0 + 5 -56835.8 1.0 2 1 + 6 48689.5 1.0 2 2 + 7 -42027.3 1.0 2 1 + 8 19353.2 1.0 2 2 + 9 4018.6 1.0 3 0 + 10 -37791.1 1.0 3 1 + 11 15926.3 1.0 3 2 + 12 -2710.5 .99 3 3 + 13 -32957.3 1.0 3 1 + 14 42084.5 1.0 3 2 + 15 -27544.2 .99 3 3 + 16 -34645.4 1.0 4 0 + 17 -8247.6 1.0 4 1 + 18 -2406.1 1.0 4 2 + 19 -11083.8 1.0 4 3 + 20 -17837.2 .99 4 4 + 21 31994.5 1.0 4 1 + 22 27811.2 1.0 4 2 + 23 -926.1 1.0 4 3 + 24 367.1 .97 4 4 + 25 -18023.6 1.0 5 0 + 26 4683.9 1.0 5 1 + 27 16160.0 1.0 5 2 + 28 -16402.0 .99 5 3 + 29 -2600.7 .99 5 4 + 30 -3660.7 .96 5 5 + 31 45347.9 1.0 5 1 + 32 -749.0 1.0 5 2 + 33 6268.5 .99 5 3 + 34 10859.6 .98 5 4 + 35 9608.4 .96 5 5 + 36 -20819.6 .99 6 0 + 37 9992.9 .99 6 1 + 38 11791.8 .99 6 2 + 39 -12574.7 .99 6 3 + 40 2669.7 .99 6 4 + 41 1113.2 .98 6 5 + 42 7584.9 .95 6 6 + 43 14533.1 .99 6 1 + 44 -10592.9 .99 6 2 + 45 568.6 .99 6 3 + 46 12871.7 .97 6 4 + 47 -4147.8 .98 6 5 + 48 3604.4 .95 6 6 + 49 598.4 .99 7 0 + 50 4665.9 .99 7 1 + 51 -6495.7 .98 7 2 + 52 -2516.5 .98 7 3 + 53 -6448.5 .99 7 4 + 54 1855.3 .97 7 5 + 55 -2892.9 .97 7 6 + 56 2968.0 .95 7 7 + 57 -7626.3 .99 7 1 + 58 -10948.4 .98 7 2 + 59 2633.3 .98 7 3 + 60 5394.2 .96 7 4 + 61 -6050.8 .97 7 5 + 62 -1526.0 .97 7 6 + 63 -5684.2 .95 7 7 + 64 10059.2 .98 8 0 + 65 1934.4 .98 8 1 + 66 -6702.9 .97 8 2 + 67 153.7 .97 8 3 + 68 -4124.2 .98 8 4 + 69 -867.2 .96 8 5 + 70 -3740.6 .96 8 6 + 71 -732.4 .96 8 7 + 72 -2433.2 .95 8 8 + 73 -2409.7 .97 8 1 + 74 -11614.6 .97 8 2 + 75 9287.0 .97 8 3 + 76 -911.9 .95 8 4 + 77 2754.5 .96 8 5 + 78 -2446.1 .96 8 6 + 79 1207.3 .96 8 7 + 80 -2887.3 .95 8 8 + 81 9671.8 .95 9 0 + 82 -3046.2 .96 9 1 + 83 260.9 .95 9 2 + 84 2071.3 .95 9 3 + 85 3329.6 .97 9 4 + 86 -2523.1 .95 9 5 + 87 1787.1 .95 9 6 + 88 -1148.2 .95 9 7 + 89 1276.5 .96 9 8 + 90 -1976.8 .95 9 9 + 91 -8467.4 .95 9 1 + 92 -1383.8 .95 9 2 + 93 5697.7 .95 9 3 + 94 -2056.3 .93 9 4 + 95 3081.5 .95 9 5 + 96 -721.2 .95 9 6 + 97 1352.5 .95 9 7 + 98 -210.1 .95 9 8 + 99 1567.6 .95 9 9 +100 -2299.5 .90 10 0 +101 2009.7 .90 10 1 +102 2127.8 .91 10 2 +103 3498.3 .91 10 3 +104 2967.6 .95 10 4 +105 16.3 .94 10 5 +106 1806.5 .94 10 6 +107 -46.5 .94 10 7 +108 2897.8 .95 10 8 +109 574.5 .95 10 9 +110 1298.9 .94 10 10 +111 -4692.6 .88 10 1 +112 4445.8 .89 10 2 +113 -2378.6 .92 10 3 +114 -2204.3 .92 10 4 +115 164.1 .94 10 5 +116 -1361.6 .94 10 6 +117 -2031.5 .94 10 7 +118 1411.8 .94 10 8 +119 -714.3 .95 10 9 +120 1676.5 .94 10 10 diff --git a/src/astroedu/planetmagfields/data/mercury.dat b/src/astroedu/planetmagfields/data/mercury.dat new file mode 100644 index 0000000..fe1b6d6 --- /dev/null +++ b/src/astroedu/planetmagfields/data/mercury.dat @@ -0,0 +1,4 @@ +g 1 0 -190 +g 2 0 -74.6 +g 3 0 -22.0 +g 4 0 -5.7 diff --git a/src/astroedu/planetmagfields/data/neptune.dat b/src/astroedu/planetmagfields/data/neptune.dat new file mode 100644 index 0000000..63bbb99 --- /dev/null +++ b/src/astroedu/planetmagfields/data/neptune.dat @@ -0,0 +1,18 @@ +g 1 0 9732 +g 1 1 3220 +h 1 0 0.0 +h 1 1 -9889 +g 2 0 7448 +g 2 1 664 +g 2 2 4499 +h 2 0 0.0 +h 2 1 11230 +h 2 2 -70 +g 3 0 -6592 +g 3 1 4098 +g 3 2 -3581 +g 3 3 484 +h 3 0 0.0 +h 3 1 -3669 +h 3 2 1791 +h 3 3 -770 diff --git a/src/astroedu/planetmagfields/data/saturn.dat b/src/astroedu/planetmagfields/data/saturn.dat new file mode 100644 index 0000000..dc4631f --- /dev/null +++ b/src/astroedu/planetmagfields/data/saturn.dat @@ -0,0 +1,14 @@ +g 1 0 21141 +g 2 0 1583 +g 3 0 2262 +g 4 0 95 +g 5 0 10.3 +g 6 0 17.4 +g 7 0 -68.8 +g 8 0 -15.5 +g 9 0 -24.2 +g 10 0 9.0 +g 11 0 11.3 +g 12 0 -2.8 +g 13 0 -2.4 +g 14 0 -0.8 diff --git a/src/astroedu/planetmagfields/data/uranus.dat b/src/astroedu/planetmagfields/data/uranus.dat new file mode 100644 index 0000000..44a7bc3 --- /dev/null +++ b/src/astroedu/planetmagfields/data/uranus.dat @@ -0,0 +1,18 @@ +g 1 0 11893 +g 1 1 11579 +h 1 0 0 +h 1 1 -15684 +g 2 0 -6030 +g 2 1 -12587 +g 2 2 196 +h 2 0 0 +h 2 1 6116 +h 2 2 4759 +g 3 0 2705 +g 3 1 1188 +g 3 2 -4808 +g 3 3 -2412 +h 3 0 0 +h 3 1 -7095 +h 3 2 -1616 +h 3 3 -2608 diff --git a/src/astroedu/planetmagfields/libbfield.py b/src/astroedu/planetmagfields/libbfield.py new file mode 100644 index 0000000..b40c279 --- /dev/null +++ b/src/astroedu/planetmagfields/libbfield.py @@ -0,0 +1,93 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import numpy as np +import matplotlib.pyplot as plt +from .libgauss import get_data,gen_idx,get_grid,getB,getBm0 +from .plotlib import * +from importlib_resources import files +import sys + +planetlist = ["mercury","earth","jupiter","saturn","uranus","neptune","ganymede"] + +datDir = files('astroedu.planetmagfields').joinpath('data/') + +def getBr(datDir=datDir,planet="earth",r=1,info=True): + + if planet not in planetlist: + print("Planet must be one of the following!") + print(planetlist) + sys.exit() + + if planet == '--help': + print("Usage: ./magField.py ") + print("Example: ./magField.py earth") + sys.exit() + + g,h,lmax,idx = get_data(datDir,planet) + + p2D,th2D = get_grid() + + if planet in ["mercury","saturn"]: + Br = getBm0(lmax,g,r,p2D,th2D) * 1e-3 + dipTheta = 0. + dipPhi = 0. + else: + Br = getB(lmax,g,h,idx,r,p2D,th2D,planet=planet) * 1e-3 + + dipTheta = np.arctan(np.sqrt(g[idx[1,1]]**2 + h[idx[1,1]]**2)/g[idx[1,0]]) * 180./np.pi + dipPhi = np.arctan(h[idx[1,1]]/g[idx[1,1]]) * 180./np.pi + + if info: + print(("Planet: %s" %planet.capitalize())) + #print(("Depth (fraction of surface radius) = %.2f" %r)) + print(("l_max = %d" %lmax)) + print(("Dipole tilt (degrees) = %f" %dipTheta)) + + return p2D, th2D, Br, dipTheta, dipPhi + +def plotAllFields(datDir=datDir,r=1.0,levels=30,cmap='RdBu_r',proj='Mollweide'): + + print("") + print('|=========|======|=======|') + print(('|%-8s | %-2s| %-5s |' %('Planet','Theta','Phi'))) + print('|=========|======|=======|') + + plt.figure(figsize=(12,12)) + for k, planet in enumerate(planetlist): + p2D,th2D,Br,dipTheta,dipPhi = getBr(datDir=datDir,planet=planet,r=r,info=False) + + if planet == "ganymede": + nplot = 8 + else: + nplot = k+1 + + if proj.lower() == 'hammer': + ax = plt.subplot(3,3,nplot) + else: + import cartopy.crs as ccrs + projection = eval('ccrs.'+proj+'()') + ax = plt.subplot(3,3,nplot,projection=projection) + + plotB_subplot(p2D,th2D,Br,ax,planet=planet,levels=levels,cmap=cmap,proj=proj) + + if planet in ["mercury","saturn"]: + print(('|%-8s | %-4.1f | %-5.1f |' %(planet.capitalize(),dipTheta, dipPhi))) + else: + print(('|%-8s | %-3.1f | %-5.1f |' %(planet.capitalize(),dipTheta, dipPhi))) + + print('|---------|------|-------|') + + if r == 1: + plt.suptitle(r'Radial magnetic field ($\mu$T) at surface', fontsize=20) + else: + plt.suptitle(r'Radial magnetic field ($\mu$T) at $r/r_{\rm surface} = %.2f$' %r, fontsize=20) + + + + +def plotMagField(datDir=datDir,planet="earth",r=1,levels=30,proj='moll',cmap='RdBu_r'): + + p2D, th2D, Br, dum1,dum2 = getBr(datDir=datDir,planet=planet,r=r) + plt.figure(figsize=(12,6.75)) + plotB(p2D,th2D,Br,planet=planet,levels=levels,proj=proj,cmap=cmap,r=r) diff --git a/src/astroedu/planetmagfields/libgauss.py b/src/astroedu/planetmagfields/libgauss.py new file mode 100644 index 0000000..053f9bd --- /dev/null +++ b/src/astroedu/planetmagfields/libgauss.py @@ -0,0 +1,349 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +import numpy as np +import scipy.special as sp +from copy import deepcopy +from importlib_resources import files + +def gen_idx(lmax): + + ''' + Generate index to convert from (l,m) to [idx] + ''' + + idx = np.zeros([lmax+1,lmax+1]) + count = 0 + + for l in range(lmax+1): + for m in range(l+1): + + idx[l,m] = count + count += 1 + + return np.int32(idx) + +def get_data(datDir,planet="earth"): + + datfile = datDir.joinpath(planet+'.dat') + + if planet == "earth": + datfile = datDir.joinpath('IGRF13.dat') + dat = np.loadtxt(datfile,usecols=[1,2,-2]) + gh = np.genfromtxt(datfile,usecols=[0],dtype='str') + + mask = gh == 'g' + gDat = dat[mask,:] + g = gDat[:,-1] + gl = gDat[:,0] + + mask = gh == 'h' + hDat = dat[mask,:] + h = hDat[:,-1] + hm = hDat[:,1] + hIdx = np.where(hm == 1.)[0] + h = np.insert(h,hIdx,0.) + + lmax = np.int32(gl.max()) + + elif planet in ["mercury","saturn"]: + dat = np.loadtxt(datfile,usecols=[3]) + g = dat.flatten() + lmax = len(g) + h = np.zeros_like(g) + + elif planet == "jupiter": + dat = np.loadtxt(datfile) + l_dat = dat[:,-2] + ghlm = dat[:,1] + + lmax = np.int32(l_dat.max()) + + g = [] + h = [] + + ######################## + # Separate glm and hlm + ######################## + + for i in range(1,lmax+1): + mask = l_dat == i + n = len(l_dat[mask]) + half = int(n/2) + + g.append(ghlm[mask][:half+1]) + h.append(np.concatenate([[0.],ghlm[mask][half+1:]])) + + g = np.concatenate(g) + h = np.concatenate(h) + + elif planet in ['uranus','neptune','ganymede']: + dat = np.loadtxt(datfile,usecols=[1,2,3]) + gh = np.genfromtxt(datfile,usecols=[0],dtype='str') + + mask = gh == 'g' + gDat = dat[mask,:] + g = gDat[:,-1] + gl = gDat[:,0] + + mask = gh == 'h' + hDat = dat[mask,:] + h = hDat[:,-1] + + gl = np.int32(gl) + + lmax = np.int32(gl.max()) + + # Insert (0,0) -> 0 for less confusion + + g = np.insert(g,0,0.) + h = np.insert(h,0,0.) + + if planet in ['mercury','saturn']: + idx = np.arange(lmax+1) + else: + idx = gen_idx(lmax) + + return g,h,lmax,idx + +def get_grid(nphi=256,ntheta=128): + + phi = np.linspace(0.,2*np.pi,nphi) + x,w = sp.roots_legendre(ntheta) + theta = np.sort(np.arccos(x)) + + p2D = np.zeros([nphi, ntheta]) + th2D = np.zeros([nphi, ntheta]) + + for i in range(nphi): + p2D[i,:] = phi[i] + + for j in range(ntheta): + th2D[:,j] = theta[j] + + return p2D, th2D + +def gen_arr(lmax, l1,m1,mode='g'): + + ''' + Generate Gauss coefficient array for testing + purposes + ''' + + idx = np.zeros([lmax+1,lmax+1]) + lArr = [] + mArr = [] + + count = 0 + + glm = [] + hlm = [] + + for l in range(lmax+1): + for m in range(l+1): + + if l in l1 and m in m1: + if mode == 'g' or mode == 'gh': + glm.append(1.) + else: + glm.append(0.) + if mode == 'h' or mode == 'gh': + hlm.append(1.) + else: + hlm.append(0.) + else: + glm.append(0.) + hlm.append(0.) + + idx[l,m] = count + lArr.append(l) + mArr.append(m) + + count += 1 + + glm = np.array(glm) + hlm = np.array(hlm) + lArr = np.array(lArr) + mArr = np.array(mArr) + idx = np.int32(idx) + + return glm, hlm, lArr, mArr, idx + +def getB(lmax,glm,hlm,idx,r,p2D,th2D,planet="earth"): + + Br = np.zeros_like(p2D) + + for l in range(1,lmax+1): + for m in range(l+1): + ylm = sp.sph_harm(m, l, p2D, th2D) + + # Include Condon-Shortley Phase for Earth but not other planets + # Scipy sph_harm has the phase included by default + + if planet in ["earth"]: + fac_m = 1. + else: + fac_m = (-1)**m + + if m != 0: + fac_m *= np.sqrt(2) + + fac = fac_m * (l+1) * r**(-l-2) * np.sqrt((4.*np.pi)/(2*l+1)) + + G = glm[idx[l,m]] * np.real(ylm) + H = hlm[idx[l,m]] * np.imag(ylm) + + Br += np.real(fac * (G + H)) + + return Br + +def getBm0(lmax,g,r,p2D,th2D): + + Br = np.zeros_like(p2D) + + for l in range(1,lmax+1): + ylm = sp.sph_harm(0, l, p2D, th2D) + fac = (l + 1) * r**(-l-2) * np.sqrt((4.*np.pi)/(2*l+1)) + + Br += fac * g[l] * np.real(ylm) + + return Br + +def get_spec(glm,hlm,idx,lmax,planet='earth',r=1): + E = np.zeros(lmax+1) + + if planet in ['mercury','saturn']: + for l in range(1,lmax+1): + E[l] = (l+1) * r**(-2*l-4) *(np.abs(glm[l])**2 + np.abs(hlm[l])**2) + emag_10 = E[1] + else: + for l in range(1,lmax+1): + for m in range(l+1): + E[l] += (l+1) * r**(-2*l-4) *(np.abs(glm[idx[l,m]])**2 + np.abs(hlm[idx[l,m]])**2) + + emag_10 = 2 * r**(-2*l-4)* np.abs(glm[idx[1,0]])**2 + return E, emag_10 + +def filt_Gauss(glm,hlm,lmax,idx,larr=None,marr=None,lCutMin=0,lCutMax=None,mmin=0,mmax=None): + + glm_filt = deepcopy(glm) + hlm_filt = deepcopy(hlm) + + if lCutMax is None: + lCutMax = lmax + if mmax is None: + mmax = lmax + + if larr is not None: + if max(larr) > lmax: + print("Error! Values in filter array must be <= lmax=%d" %lmax) + else: + for ell in range(lmax+1): + if ell not in larr: + glm_filt[idx[ell,:]] = 0. + hlm_filt[idx[ell,:]] = 0. + else: + if lCutMax > lmax or lCutMin > lmax: + print("Error! lCutMin/lCutMax must be <= lmax = %d" %lmax) + else: + for ell in range(lCutMin): + glm_filt[idx[ell,:]] = 0. + hlm_filt[idx[ell,:]] = 0. + for ell in range(lCutMax,lmax+1): + glm_filt[idx[ell,:]] = 0. + hlm_filt[idx[ell,:]] = 0. + + if marr is not None: + if max(marr) > lmax: + print("Error! Values in filter array must be <= lmax=%d" %lmax) + else: + for m in range(lmax+1): + if m not in marr: + glm_filt[idx[:,m]] = 0. + hlm_filt[idx[:,m]] = 0. + else: + if mmin > lmax or mmax > lmax: + print("Error! mmin/mmax must be <= lmax = %d" %lmax) + else: + for m in range(mmin): + glm_filt[idx[:,m]] = 0. + hlm_filt[idx[:,m]] = 0. + for m in range(mmax+1,lmax+1): + glm_filt[idx[:,m]] = 0. + hlm_filt[idx[:,m]] = 0. + + return glm_filt,hlm_filt + +def filt_Gaussm0(glm,hlm,lmax,larr=None,lCutMin=0,lCutMax=None): + + glm_filt = deepcopy(glm) + hlm_filt = deepcopy(hlm) + + if lCutMax is None: + lCutMax = lmax + + if larr is not None: + if max(larr) > lmax: + print("Error! Values in filter array must be <= lmax=%d" %lmax) + else: + for ell in range(lmax+1): + if ell not in larr: + glm_filt[ell] = 0. + hlm_filt[ell] = 0. + else: + if lCutMax > lmax or lCutMin > lmax: + print("Error! lCutMin/lCutMax must be <= lmax = %d" %lmax) + else: + for ell in range(lCutMin): + glm_filt[ell] = 0. + hlm_filt[ell] = 0. + for ell in range(lCutMax,lmax+1): + glm_filt[ell] = 0. + hlm_filt[ell] = 0. + + return glm_filt,hlm_filt + + +def sphInt(f,g,phi,th2D,theta): + + thetaInt = np.trapz(f * g * np.sin(th2D),theta,axis=1) + phiInt = np.trapz(thetaInt,phi) + + return phiInt + + +def getGauss(lmax,Br,r,phi,theta,th2D,p2D): + + ''' + Get Gauss coefficients from a surface field + ''' + glm = [] + hlm = [] + + comp = complex(0,1) + + for l in range(0,lmax+1): + for m in range(0,l+1): + + ylm = (-1)**m * sp.sph_harm(m, l, p2D, th2D) + + ylm_conj = np.conjugate(ylm) + + fac = r**(l+2)/(l+1) + + if m==0: + fac *= 0.5 + + I1 = sphInt(Br,ylm,phi,th2D,theta) + I2 = sphInt(Br,ylm_conj,phi,th2D,theta) + + g = fac * (I2 + I1) + h = comp * fac * (I2 - I1) + + glm.append(g) + hlm.append(h) + + glm = np.array(glm) + hlm = np.array(hlm) + + return glm, hlm diff --git a/src/astroedu/planetmagfields/planet.py b/src/astroedu/planetmagfields/planet.py new file mode 100644 index 0000000..a9f1a51 --- /dev/null +++ b/src/astroedu/planetmagfields/planet.py @@ -0,0 +1,160 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import numpy as np +import matplotlib.pyplot as plt +from .libgauss import get_data, filt_Gauss, filt_Gaussm0,getB, getBm0, get_spec +from .libbfield import getBr +from .plotlib import plotB, plotSurf, plot_spec +from importlib_resources import files +import sys + +planetlist = ["mercury","earth","jupiter","saturn","uranus","neptune","ganymede"] + +datDir = files('astroedu.planetmagfields').joinpath('data/') + +class planet: + + def __init__(self,name='earth',datDir=datDir): + + self.name = name.lower() + if self.name not in planetlist: + print("Planet must be one of the following!") + print(planetlist) + + self.datDir = datDir + self.glm, self.hlm, self.lmax, self.idx = \ + get_data(self.datDir,planet=self.name) + + self.p2D, self.th2D, self.Br, self.dipTheta, self.dipPhi = \ + getBr(datDir=self.datDir,planet=self.name,r=1,info=True) + + self.phi = self.p2D[:,0] + self.theta = self.th2D[0,:] + self.r = 1 + + def plot(self,r=1,levels=30,cmap='RdBu_r',proj='Mollweide'): + plt.figure(figsize=(12,6.75)) + + if r == 1: + ax,cbar = plotSurf(self.p2D,self.th2D,self.Br,levels=levels,cmap=cmap,proj=proj) + else: + self.p2D, self.th2D, self.Br, self.dipTheta, self.dipPhi = \ + getBr(datDir=self.datDir,planet=self.name,r=r,info=False) + self.r = r + ax,cbar = plotSurf(self.p2D,self.th2D,self.Br,levels=levels,cmap=cmap,proj=proj) + + cbar.ax.set_xlabel(r'Radial magnetic field ($\mu$T)',fontsize=25) + cbar.ax.tick_params(labelsize=20) + + if r==1: + radLabel = ' Surface' + else: + radLabel = r' $r/r_{\rm surface}=%.2f$' %r + + if proj.lower() != 'hammer' and self.name == 'earth': + ax.coastlines() + + if r==1: + radLabel = ' Surface' + else: + radLabel = r' $r/r_{\rm surface}=%.2f$' %r + + ax.set_title(self.name.capitalize() + radLabel,fontsize=25,pad=20) + plt.tight_layout() + + def writeVtsFile(self,potExtra=False,ratio_out=2,nrout=32): + from .potextra import extrapot, writeVts + + rout = np.linspace(1,ratio_out,nrout) + if potExtra: + brout, btout, bpout = extrapot(self.lmax,1.,self.Br,rout) + else: + brout = self.Br + btout = np.zeros_like(self.Br) + bpout = np.zeros_like(self.Br) + + writeVts(self.name,brout,btout,bpout,rout,self.theta,self.phi) + + ## Filtered plots + + def plot_filt(self,r=1,larr=None,marr=None,lCutMin=0,lCutMax=None,mmin=0,mmax=None,levels=30,cmap='RdBu_r',proj='Mollweide'): + + self.larr_filt = larr + self.marr_filt = marr + self.lCutMin = lCutMin + self.lCutMax = lCutMax + self.mmin_filt = mmin + self.mmax_filt = mmax + + if self.lCutMax is None: + self.lCutMax = self.lmax + if self.mmax_filt is None: + self.mmax_filt = self.lmax + + self.r_filt = r + + if self.name in ['mercury','saturn']: + self.glm_filt,self.hlm_filt =\ + filt_Gaussm0(self.glm,self.hlm,self.lmax,larr=self.larr_filt,\ + lCutMin=self.lCutMin,lCutMax=self.lCutMax) + self.Br_filt = 1e-3*getBm0(self.lmax,self.glm_filt,self.r_filt,self.p2D,self.th2D) + else: + self.glm_filt,self.hlm_filt =\ + filt_Gauss(self.glm,self.hlm,self.lmax,self.idx,larr=self.larr_filt,\ + marr=self.marr_filt,lCutMin=self.lCutMin,lCutMax=self.lCutMax,mmin=self.mmin_filt,mmax=self.mmax_filt) + + self.Br_filt = 1e-3*getB(self.lmax,self.glm_filt,self.hlm_filt,self.idx,self.r_filt,self.p2D,self.th2D,planet=self.name) + + plt.figure(figsize=(12,6.75)) + + ax,cbar = plotSurf(self.p2D,self.th2D,self.Br_filt,levels=levels,cmap=cmap,proj=proj) + + if r==1: + radLabel = ' Surface' + else: + radLabel = r' $r/r_{\rm surface}=%.2f$' %r + + if self.larr_filt is not None: + elllabel = r', $l = %s$' %np.str(self.larr_filt) + else: + if self.lCutMin > 0: + if self.lCutMax < self.lmax: + elllabel = r', $ %d \leq l \leq %d$' %(self.lCutMin,self.lCutMax) + else: + elllabel = r', $l \geq %d$' %self.lCutMin + + elif self.lCutMax < self.lmax: + elllabel = r', $l \leq %d$' %self.lCutMax + + if self.marr_filt is not None: + elllabel += r', $m = %s$' %np.str(self.marr_filt) + else: + if self.mmin_filt > 0: + if self.mmax_filt < self.lmax: + elllabel += r', $ %d \leq m \leq %d$' %(self.mmin_filt,self.mmax_filt) + else: + elllabel += r', $m \geq %d$' %self.mmin_filt + elif self.mmax_filt < self.lmax: + elllabel += r', $m \leq %d$' %self.mmax_filt + + cbar.ax.set_xlabel(r'Radial magnetic field ($\mu$T)',fontsize=25) + cbar.ax.tick_params(labelsize=20) + + if proj.lower() != 'hammer' and self.name == 'earth': + ax.coastlines() + ax.set_title(self.name.capitalize() + radLabel + elllabel,fontsize=25,pad=20) + plt.tight_layout() + + + def spec(self,r=1,iplot=True): + self.emag_spec, emag_10 = get_spec(self.glm,self.hlm,self.idx,self.lmax,planet=self.name,r=r) + l = np.arange(self.lmax+1) + + self.dip_tot = self.emag_spec[1]/sum(self.emag_spec) + self.dipolarity = emag_10/sum(self.emag_spec) + if iplot: + plt.figure(figsize=(7,7)) + plot_spec(l,self.emag_spec,r,self.name) + plt.tight_layout() + plt.show() \ No newline at end of file diff --git a/src/astroedu/planetmagfields/plotlib.py b/src/astroedu/planetmagfields/plotlib.py new file mode 100644 index 0000000..929d621 --- /dev/null +++ b/src/astroedu/planetmagfields/plotlib.py @@ -0,0 +1,146 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.colors as colors + +def hammer2cart(ttheta, pphi, colat=False): + """ + This function is used to define the Hammer projection for + default plotting. Copied from MagIC python plotting script: + + https://github.com/magic-sph/magic/blob/master/python/magic/plotlib.py + """ + + if not colat: # for lat and phi \in [-pi, pi] + xx = 2.*np.sqrt(2.) * np.cos(ttheta)*np.sin(pphi/2.)\ + /np.sqrt(1.+np.cos(ttheta)*np.cos(pphi/2.)) + yy = np.sqrt(2.) * np.sin(ttheta)\ + /np.sqrt(1.+np.cos(ttheta)*np.cos(pphi/2.)) + else: # for colat and phi \in [0, 2pi] + xx = -2.*np.sqrt(2.) * np.sin(ttheta)*np.cos(pphi/2.)\ + /np.sqrt(1.+np.sin(ttheta)*np.sin(pphi/2.)) + yy = np.sqrt(2.) * np.cos(ttheta)\ + /np.sqrt(1.+np.sin(ttheta)*np.sin(pphi/2.)) + return xx, yy + +def plotSurf(p2D,th2D,B,levels=60,cmap='RdBu_r',proj='Mollweide'): + + bmax = np.abs(B).max() + digits = int(np.log10(bmax)) + 1 + + if digits > 1: + bmax = np.round(bmax) + else: + bmax = np.round(bmax,decimals=1) + + divnorm = colors.TwoSlopeNorm(vmin=-bmax, vcenter=0, vmax=bmax) + cs = np.linspace(-bmax,bmax,levels) + + + lon2D = p2D - np.pi + lat2D = np.pi/2 - th2D + + try: + import cartopy.crs as ccrs + except: + print("cartopy library not available, using Hammer projection") + proj = 'hammer' + + if proj.lower() == 'hammer': + ax = plt.axes() + xx,yy = hammer2cart(lat2D,lon2D) + cont = ax.contourf(xx,yy,B,cs,cmap=cmap,norm=divnorm,extend='both') + else: + projection = eval('ccrs.'+proj+'()') + + ax = plt.axes(projection=projection) + + cont = ax.contourf(lon2D*180/np.pi,lat2D*180/np.pi,B,cs, \ + transform=ccrs.PlateCarree(),cmap=cmap,norm=divnorm,extend='both') + + cbar = plt.colorbar(cont,orientation='horizontal',fraction=0.06, pad=0.04,ticks=[-bmax,0,bmax]) + + ax.axis('equal') + ax.axis('off') + + return ax, cbar + +def plotB(p2D,th2D,B,r=1,planet="earth",levels=60,cmap='RdBu_r',proj='Mollweide'): + + planet = planet.lower() + + ax,cbar = plotSurf(p2D,th2D,B,levels=levels,cmap=cmap,proj=proj) + cbar.ax.set_xlabel(r'Radial magnetic field ($\mu$T)',fontsize=25) + cbar.ax.tick_params(labelsize=20) + + if proj.lower() != 'hammer' and planet == 'earth': + ax.coastlines() + + if r==1: + radLabel = ' Surface' + else: + radLabel = r' $r/r_{\rm surface}=%.2f$' %r + + ax.set_title(planet.capitalize() + radLabel,fontsize=25,pad=20) + +def plotB_subplot(p2D,th2D,B,ax,planet="earth",levels=60,cmap='RdBu_r',proj='Mollweide'): + planet = planet.lower() + + bmax = np.abs(B).max() + digits = int(np.log10(bmax)) + 1 + + if digits > 1: + bmax = np.round(bmax) + else: + bmax = np.round(bmax,decimals=1) + + p2D -= np.pi + th2D -= np.pi/2 + th2D = -th2D + + cs = np.linspace(-bmax,bmax,levels) + divnorm = colors.TwoSlopeNorm(vmin=-bmax, vcenter=0, vmax=bmax) + + try: + import cartopy.crs as ccrs + except: + print("cartopy library not available, using Hammer projection") + proj = 'hammer' + + if proj.lower() == 'hammer': + xx,yy = hammer2cart(th2D,p2D) + cont = ax.contourf(xx,yy,B,cs,cmap=cmap,norm=divnorm,extend='both') + else: + if planet == "earth": + ax.coastlines() + + cont = ax.contourf(p2D*180/np.pi,th2D*180/np.pi,B,cs, \ + transform=ccrs.PlateCarree(),cmap=cmap,norm=divnorm,extend='both') + + cbar = plt.colorbar(cont,orientation='horizontal',fraction=0.06, pad=0.04,ticks=[-bmax,0,bmax]) + cbar.ax.tick_params(labelsize=15) + + ax.set_title(planet.capitalize(),fontsize=20) + ax.axis('equal') + ax.axis('off') + +def plot_spec(l,E,r,planet): + + plt.semilogy(l[1:],E[1:],'-o',lw=1.2,color='#449c99',mfc='#347b79',mec='#347b79',ms=8) + plt.xlabel(r'$l$',fontsize=30) + plt.ylabel(r'$R_l$ (nT$^2$)',fontsize=30) + plt.xticks(l) + + plt.tick_params(which='major',labelsize=15,length=10) + plt.tick_params(which='minor',length=5) + + plt.grid(True,alpha=0.5) + + if r==1: + radLabel = ' Surface' + else: + radLabel = r' $r/r_{\rm surface}=%.2f$' %r + + plt.title('Lowes spectrum, '+ planet.capitalize() + radLabel,fontsize=20,pad=10) diff --git a/src/astroedu/planetmagfields/potextra.py b/src/astroedu/planetmagfields/potextra.py new file mode 100644 index 0000000..533337b --- /dev/null +++ b/src/astroedu/planetmagfields/potextra.py @@ -0,0 +1,120 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import numpy as np + +def extrapot(lmax,rcmb,brcmb,rout): + + nphi, ntheta = brcmb.shape + nrout = len(rout) + + polar_opt = 1e-10 + + lmax = int(nphi/3) + mmax = lmax + + try: + import shtns + except ImportError: + print("Potential extrapolation requires the SHTns library") + print("It can be obtained here: https://bitbucket.org/nschaeff/shtns") + + norm=shtns.sht_orthonormal | shtns.SHT_NO_CS_PHASE + + sh = shtns.sht(lmax,mmax=mmax,norm=norm) + ntheta, nphi = sh.set_grid(ntheta, nphi, polar_opt=polar_opt) + + + L = sh.l * (sh.l + 1) + + brlm = sh.analys(brcmb.T) + bpolcmb = np.zeros_like(brlm) + bpolcmb[1:] = rcmb**2 * brlm[1:]/L[1:] + btor = np.zeros_like(brlm) + + brout = np.zeros([ntheta,nphi,nrout]) + btout = np.zeros([ntheta,nphi,nrout]) + bpout = np.zeros([ntheta,nphi,nrout]) + + for k,radius in enumerate(rout): + print(("%d/%d" %(k,nrout))) + + radratio = rcmb/radius + bpol = bpolcmb * radratio**(sh.l) + brlm = bpol * L/radius**2 + brout[...,k] = sh.synth(brlm) + + dbpoldr = -sh.l/radius * bpol + slm = dbpoldr + + btout[...,k], bpout[...,k] = sh.synth(slm,btor) + + brout = np.transpose(brout,(1,0,2)) + btout = np.transpose(btout,(1,0,2)) + bpout = np.transpose(bpout,(1,0,2)) + + return brout, btout, bpout + +def get_grid(r,theta,phi): + + nr,ntheta,nphi = len(r),len(theta),len(phi) + + r3D = np.zeros([nphi,ntheta,nr]) + th3D = np.zeros([nphi,ntheta,nr]) + p3D = np.zeros([nphi,ntheta,nr]) + + for i in range(nr): + r3D[...,i] = r[i] + for j in range(ntheta): + th3D[:,j,:] = theta[j] + for k in range(nphi): + p3D[k,...] = phi[k] + + s = r3D * np.sin(th3D) + x = s * np.cos(p3D) + y = s * np.sin(p3D) + z = r3D * np.cos(th3D) + + return r3D,th3D,p3D, x,y,z, s + +def get_cart(vr,vt,vp,r3D,th3D,p3D): + + vs = vr * np.sin(th3D) + vt * np.cos(th3D) + vz = vr * np.cos(th3D) - vt * np.sin(th3D) + + vx = vs * np.cos(p3D) - vp * np.sin(p3D) + vy = vs * np.sin(p3D) + vp * np.cos(p3D) + + return vx,vy,vz + +def writeVts(name,br,btheta,bphi,r,theta,phi): + + r3D,th3D,p3D, x,y,z, s = get_grid(r,theta,phi) + + print("grid shape=",th3D.shape) + + bx,by,bz = get_cart(br, btheta, bphi,r3D,th3D,p3D) + + try: + try: # Version 2 changed naming convention of functions + #import evtk + from evtk.hl import structuredToVTK + #gridToVTK = evtk.hl.structuredToVTK + gridToVTK = structuredToVTK + except: + import evtk + gridToVTK = evtk.hl.gridToVTK + except: + print("movie2vtk requires the use of evtk library!") + print("You can get it from https://github.com/paulo-herrera/PyEVTK") + + + br = np.asfortranarray(br) + bx = np.asfortranarray(bx) + by = np.asfortranarray(by) + bz = np.asfortranarray(bz) + + gridToVTK("%s"%name,x,y,z,pointData= {"radius":r3D, + "Radial mag field":br, + "Mag Field":(bx, by,bz)}) +