diff --git a/.github/ISSUE_TEMPLATE/bug_report.md b/.github/ISSUE_TEMPLATE/bug_report.md index 4395f95..903de26 100644 --- a/.github/ISSUE_TEMPLATE/bug_report.md +++ b/.github/ISSUE_TEMPLATE/bug_report.md @@ -31,9 +31,9 @@ assignees: '' -- **OS:** [e.g., macOS 13.0, Ubuntu 22.04, Windows 11] +- **OS:** [e.g., macOS 26.0, Ubuntu 22.04, Windows 11] - **Python Version:** [e.g., 3.9.7, 3.10.5] -- **mrfmsim Version:** [e.g., 0.1.0, or commit hash if installed from source] +- **mrfmsim Version:** [e.g., 0.1.0] - **Installation Method:** [e.g., pip, conda, from source] ## Error Messages/Stack Traces @@ -54,4 +54,5 @@ Paste error message here ## Possible Solution - + diff --git a/.github/ISSUE_TEMPLATE/documentation.md b/.github/ISSUE_TEMPLATE/documentation.md index 0510d7b..2ac2224 100644 --- a/.github/ISSUE_TEMPLATE/documentation.md +++ b/.github/ISSUE_TEMPLATE/documentation.md @@ -24,6 +24,5 @@ assignees: '' ## Additional Context - - - + diff --git a/.github/ISSUE_TEMPLATE/feature_request.md b/.github/ISSUE_TEMPLATE/feature_request.md index db3f68c..c3d019e 100644 --- a/.github/ISSUE_TEMPLATE/feature_request.md +++ b/.github/ISSUE_TEMPLATE/feature_request.md @@ -28,11 +28,10 @@ assignees: '' ## Additional Context - + ## Potential Implementation - - - - + diff --git a/.github/pull_request_template.md b/.github/pull_request_template.md index eeaf923..2aeaa48 100644 --- a/.github/pull_request_template.md +++ b/.github/pull_request_template.md @@ -1,3 +1,15 @@ +## Description + + + +## Changes Made + + + +- +- +- + ## Type of Pull Request @@ -11,10 +23,6 @@ - [ ] Test coverage improvement - [ ] Other (please describe): -## Description - - - ## Related Issue(s) @@ -22,25 +30,17 @@ -## Changes Made - - - -- -- -- - ## Breaking Changes -## Testing Performed +## Testing -- [ ] All existing tests pass (`tox` or `pytest`) +- [ ] All existing tests pass (`tox`) - [ ] Added new tests for new functionality - [ ] Manual testing performed (describe below) @@ -58,31 +58,18 @@ - [ ] Added usage examples (if applicable) - [ ] Documentation builds without errors (`cd docs && make html`) -## Code Quality +## Checklist -- [ ] Code follows PEP 8 style guidelines +- [ ] Code follows PEP 8 style guidelines (use ``Black`` formatter's default settings) - [ ] Used meaningful variable and function names - [ ] Code is modular and readable - -## Checklist - - - -- [ ] My code follows the project's coding standards -- [ ] I have performed a self-review of my own code -- [ ] I have commented my code, particularly in hard-to-understand areas -- [ ] My changes generate no new warnings or errors -- [ ] I have checked that this PR focuses on a single feature or bug fix +- [ ] Commented code, particularly in hard-to-understand areas +- [ ] Changes generate no new warnings or errors +- [ ] PR focuses on a single feature or bug fix ## Additional Notes - - - -## Screenshots (if applicable) - - - - + diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 7a61f84..39e7e6f 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -10,13 +10,21 @@ and this project adheres to [unreleased] ------------ -Minor fixes to documentation. +Added Cylinder Magnet and minor fixes to docstrings and documentation. Changed ^^^^^^^ - Change the CERMITD experiment graph to avoid output duplication. +- Change the API to use the latest version of the mmodel package (0.9.0). +- Modify the Nan treatment of ``polarization.rel_dpol_sat_td`` (#4) +- Reformat the docstrings and documentations +Added +^^^^^ + +- Add ``CylinderMagnetApprox`` class to approximate the magnetic field of a cylinder magnet (#4) +- Add pull request and issue templates. [0.3.0] - 2025-02-25 ---------------------- diff --git a/README.rst b/README.rst index 73ca6fc..01f5eef 100644 --- a/README.rst +++ b/README.rst @@ -65,6 +65,11 @@ To test in different environments:: tox +Contributing +^^^^^^^^^^^^ + +We welcome contributions! Please see our `Contributing Guide `_ for details on how to follow our development guidelines and submit pull requests. + .. rubric:: References .. [#Sidles1995jan] Sidles, J. A.; Garbini, J. J.; Bruland, K. J.; Rugar, D.; diff --git a/docs/basic/geometry.rst b/docs/basic/geometry.rst index c7d0e5e..757a4c4 100644 --- a/docs/basic/geometry.rst +++ b/docs/basic/geometry.rst @@ -1,9 +1,9 @@ Experiment geometry ============================= -For experiments in the *mrfmsim*, two experimental geometry +For experiments in the *mrfmsim*, two experimental geometries are used: "hangdown" and SPAM (springiness preservation by aligning magnetization). -The two geometry are shown in the figure below. +The two geometries are shown in the figure below. .. image:: ../_static/mrfm-setups.pdf :width: 800px diff --git a/docs/basic/inheritance.rst b/docs/basic/inheritance.rst index 8768778..e2514b2 100644 --- a/docs/basic/inheritance.rst +++ b/docs/basic/inheritance.rst @@ -67,5 +67,5 @@ ExperimentGroup ---------------- The ``mrfmsim.group.ExperimentGroup`` class adds additional properties -``epxeriments`` and ``epxeriment_defaults`` alongside the parent class's +``experiments`` and ``experiment_defaults`` alongside the parent class's ``models`` and ``model_defaults`` properties. diff --git a/docs/basic/plugin.rst b/docs/basic/plugin.rst index c5882a9..9d04ed7 100644 --- a/docs/basic/plugin.rst +++ b/docs/basic/plugin.rst @@ -5,13 +5,13 @@ Ancillary features are added to the mrfmsim ecosystem through plugin packages. The separation of the core functionality and the plugins allows for an easy maintenance of the mrfmsim package. The core functionalities do not depend on the plugins. -mrfmsim-ymal: YAML configuration plugin +mrfmsim-yaml: YAML configuration plugin --------------------------------------- mrfmsim classes can be fully defined using YAML configuration files. The plugin creates a ``MrfmSimLoader`` and a ``MrfmSimDumper`` that can be used to load and dump the mrfmsim classes to and from YAML files. For detailed usage of the package -see the `mrfmsim-ymal `_. +see the `mrfmsim-yaml `_. To define an experiment or experiment group with YAML files, the available tags are: "!func", "!import", "!nodes", "!Graph", "!Experiment", and "!ExperimentGroup". diff --git a/docs/contribute.rst b/docs/contribute.rst new file mode 100644 index 0000000..9daa7b4 --- /dev/null +++ b/docs/contribute.rst @@ -0,0 +1,158 @@ +Contribute +======================== + +We welcome contributions to mrfmsim! This guide will help you get started. + +Getting Started +--------------- + +1. **Fork the Repository** + + Fork the mrfmsim repository on GitHub and clone your fork locally: + + .. code-block:: bash + + git clone https://github.com/Marohn-Group/mrfmsim.git + cd mrfmsim + + +2. **Create a Branch** + + Create a new branch for your feature or bug fix: + + .. code-block:: bash + + git checkout -b feature/your-branch-name + +Development Guidelines +---------------------- + +Code Style +~~~~~~~~~~ + +- Follow PEP 8 style guidelines for Python code (use ``Black`` formatter's default settings) +- Use meaningful variable and function names +- Keep functions focused and modular +- Add type hints where appropriate + +Documentation +~~~~~~~~~~~~~ + +- Write clear docstrings for all public classes, methods, and functions +- Use reStructuredText (reST) format for docstrings +- Refer to the `Sphinx docstring guide `__ for detailed formatting information +- Include: + + - Brief description of what the function/class does + - Parameter descriptions with types and units in **Sphinx** docstring format + - Return value description + - Usage examples where helpful + - Mathematical formulas using Sphinx math directive when relevant + +Example docstring format: + +.. code-block:: python + + def example_method(x, y, z): + r"""Calculate some property. + + Brief description of what this method does and any important + details about the calculation or algorithm. + + :param float x: :math:`x` coordinate [nm] + :param float y: :math:`y` coordinate [nm] + :param float z: :math:`z` coordinate [nm] + :return: calculated value + :rtype: float + """ + pass + + +Building Documentation +---------------------- + +To build the documentation locally: + +.. code-block:: bash + + cd docs + make html + +The built documentation will be in ``docs/_build/html/``. + +Testing +~~~~~~~ + +- Write tests for all new features and bug fixes +- Tests should be fast and independent from one another +- Ensure all tests pass before submitting a pull request +- Aim for good test coverage of your code + +Install test dependencies: + +.. code-block:: bash + + pip install ".[test]" + +Run tests with tox: + +.. code-block:: bash + + tox + +This will run tests across multiple Python environments and check code coverage. + + +Submitting a Pull Request +-------------------------- + +**Important**: Each pull request should focus on a single feature or bug fix. This makes the code easier to review, test, and maintain. + +1. **Push to Your Fork** + + .. code-block:: bash + + git push origin feature/your-feature-name + +2. **Create Pull Request** + + - Go to the mrfmsim repository on GitHub + - Click "New Pull Request" + - Select your fork and the ``develop`` branch as the base branch + - Fill out the pull request template with: + + - Type of pull request (bug fix, feature, documentation, etc.) + - Description of changes + - Related issue numbers (if applicable) + - Any breaking changes + - Testing performed + +3. **Review Process** + + - Maintainers will review your pull request + - Address any feedback or requested changes + - Once approved, your changes will be merged + + +Reporting Issues +---------------- + +When reporting issues, please include: + +- A clear description of the problem +- Steps to reproduce the issue +- Expected vs. actual behavior +- Your environment (Python version, OS, etc.) +- Any relevant error messages or stack traces + +Questions? +---------- + +If you have questions about contributing, feel free to: + +- Open an issue on GitHub +- Check existing documentation +- Reach out to the maintainers + +Thank you for contributing to mrfmsim! + diff --git a/docs/experiment/CermitSat.rst b/docs/experiment/CermitSat.rst index df868c5..db2e26d 100644 --- a/docs/experiment/CermitSat.rst +++ b/docs/experiment/CermitSat.rst @@ -101,7 +101,7 @@ in the Table below. ======================================== ================================================= Variables appearing in eqs. :eq:`eq:x`, :eq:`eq:x-vars`, and subsequent equations. -he magnetic field :math:`B_0` is applied along the :math:`z` axis, +The magnetic field :math:`B_0` is applied along the :math:`z` axis, and the equilibrium magnetization :math:`M_0` also lies along the :math:`z` axis. Steady-state solution diff --git a/docs/experiment/CermitSingleSpin.rst b/docs/experiment/CermitSingleSpin.rst index 44109e1..7748f12 100644 --- a/docs/experiment/CermitSingleSpin.rst +++ b/docs/experiment/CermitSingleSpin.rst @@ -5,7 +5,7 @@ Overview -------- The experiments calculate the effective force on a cantilever from a single electron spin located directly below a magnet-tipped cantilever in the "hangdown" and SPAM geometries. -There are two experiments in the group: cantilever spin constant shift calculated using the Trapezoid rule ("CermitSingleSpinApprox"), and cantilever spring constant shift calculated directly using elliptic integrals ("CermitSingleSpin"). +There are two experiments in the group: cantilever spring constant shift calculated using the Trapezoid rule ("CermitSingleSpinApprox"), and cantilever spring constant shift calculated directly using elliptic integrals ("CermitSingleSpin"). The exact solution only works for spin directly under a **spherical** magnet. The experiments are adopted from Section 3.5 of Eric Moore's dissertation. [#Moore2011Sep]_ diff --git a/docs/experiment/CermitTD.rst b/docs/experiment/CermitTD.rst index 1bcf269..0d19972 100644 --- a/docs/experiment/CermitTD.rst +++ b/docs/experiment/CermitTD.rst @@ -26,7 +26,7 @@ At low-:math:`B_1` region, we can determine numerically the final magnetization .. math:: :label: eq:Mz(tau0) - M_z^\mathrm{finial} + M_z^\mathrm{final} \approx e^{-R(\tau_i, \tau_f)} M_z^\mathrm{initial} diff --git a/docs/experiment/IBMCyclic.rst b/docs/experiment/IBMCyclic.rst index 2f2bc54..8c85896 100644 --- a/docs/experiment/IBMCyclic.rst +++ b/docs/experiment/IBMCyclic.rst @@ -34,7 +34,7 @@ as the MRFM signal. For :math:`n` independent configurations of the spin ensembl the sample variance :math:`s^2_{\Delta N}` is .. math:: - s^2_{\Delta N} = \frac{1}{n-1} \sum^n_{j=1} (\Delta N_j - \overline{\Delta N}) + s^2_{\Delta N} = \frac{1}{n-1} \sum^n_{j=1} (\Delta N_j - \overline{\Delta N})^2 where :math:`\overline{\Delta N}` is the Boltzmann polarization. diff --git a/docs/index.rst b/docs/index.rst index bd48edc..8610f65 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -15,6 +15,7 @@ License overview tests changelog + contribute .. toctree:: diff --git a/docs/overview.rst b/docs/overview.rst index 3403598..70f97d7 100644 --- a/docs/overview.rst +++ b/docs/overview.rst @@ -93,7 +93,7 @@ To access an experiment model from a experiment group: .. code:: python - # print c summary + # print a summary print(CermitESRGroup) # list experiments print(list(CermitESRGroup.experiments.keys())) @@ -184,7 +184,7 @@ can be very computationally intensive. If we want to simulate the change of the signal over a wide range of external field (:math:`B_0`) and microwave frequency -(:math:`f_\mathrm{rf}``), we would want to avoid repeated +(:math:`f_\mathrm{rf}`), we would want to avoid repeated calculations of unnecessary components that are independent of the two parameters. @@ -263,7 +263,7 @@ This is equivalent to the result from the following loops: .. note:: - Note that for individual parameters, the loop shortcut can achieve + For individual parameters, the loop shortcut can achieve optimal looping. However, for multiple parameters, users must decide which parameter to loop first. Since all nodes that are dependent on "f_rf" also depend on "B_0", we loop "f_rf" first. diff --git a/docs/reference/ref_magnetization.rst b/docs/reference/ref_magnetization.rst index 5f0e851..2b449e9 100644 --- a/docs/reference/ref_magnetization.rst +++ b/docs/reference/ref_magnetization.rst @@ -43,7 +43,7 @@ the equilibrium magnetization tends towards the Curie-Weiss law, .. math:: - mu_z^{\text{eq}} + \mu_z^{\text{eq}} \approx \dfrac{\hbar^2 \gamma^2 \: J (J + 1)}{3 \: k_b T} B_0 diff --git a/docs/setup/cantilever.rst b/docs/setup/cantilever.rst index 9c78698..f81836d 100644 --- a/docs/setup/cantilever.rst +++ b/docs/setup/cantilever.rst @@ -9,11 +9,11 @@ Example Usage ^^^^^^^^^^^^^ Create a cantilever object with a spring constant of 780 aN/nm and -frequency of 4795 Hz: +frequency of 4975 Hz: .. code-block:: python - cantilever = Cantilever(k_c=780, f_c=4795) + cantilever = Cantilever(k_c=780, f_c=4975) To print out the cantilever summary: diff --git a/docs/setup/coordinates.rst b/docs/setup/coordinates.rst index 7b1f79b..ab740a2 100644 --- a/docs/setup/coordinates.rst +++ b/docs/setup/coordinates.rst @@ -39,8 +39,8 @@ With the right figure setup, magnet = RectangularMagnet([70, 70, 1500], mu0_Ms=1800.0, magnet_origin=[0.0, 750.0, 0.0]) The grid has a width of 220 nm in :math:`y`-direction, and the magnet has a height of -1500 nm. The grid origin is set to be 110 nm below left of the :math:`xz`-plane, -and the magnet origin is set to be 750 nm right to the :math:`xz`-plane. The magnet and the +1500 nm. The grid origin is set to be 110 nm to the left of the :math:`xz`-plane, +and the magnet origin is set to be 750 nm to the right of the :math:`xz`-plane. The magnet and the grid are flush, and the tip-sample separation is 0 nm. .. note:: diff --git a/docs/setup/grid.rst b/docs/setup/grid.rst index 2ca0483..66f3693 100644 --- a/docs/setup/grid.rst +++ b/docs/setup/grid.rst @@ -14,8 +14,8 @@ To create a grid object: grid = Grid(grid_shape=[21, 11, 101], grid_step=[20.0, 4.0, 20.0], grid_origin=[0.0, 0.0, 0.0]) The resulting grid is rectangular with the shape of (21, 11, 101), the distance between -two edge grid points in each direction are (200, 40, 2000) nm, and the full grid size -is (220 :math:`\times` 44 :math:`\times` 2020) nm. +two edge grid points in each direction are (400, 40, 2000) nm, and the full grid size +is (420 :math:`\times` 44 :math:`\times` 2020) nm. The returned grid points are numpy "ogrid" points: `np.ogrid `__. diff --git a/docs/setup/magnet.rst b/docs/setup/magnet.rst index 7dbcfba..1cd5b0d 100644 --- a/docs/setup/magnet.rst +++ b/docs/setup/magnet.rst @@ -12,6 +12,7 @@ Currently, the two types of magnets supported are mrfmsim.component.magnet.SphereMagnet mrfmsim.component.magnet.RectangularMagnet + mrfmsim.component.cylindermagnet.CylinderMagnetApprox Example Usage ------------- @@ -41,3 +42,12 @@ To display the magnet information: :members: :undoc-members: :show-inheritance: + + +:mod:`cylindermagnet` module +---------------------------- + +.. automodule:: mrfmsim.component.cylindermagnet + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/tests.rst b/docs/tests.rst index 641996b..f553cad 100644 --- a/docs/tests.rst +++ b/docs/tests.rst @@ -36,8 +36,8 @@ Magnet Test The RectMagnet is tested through its symmetry factors: it is symmetric in the :math:`z` direction when x and y are 0. The gradient functions are tested against a numerical derivative calculation from the field. :math:`B_z` is an even function in the :math:`x` direction, -:math:`B_{zx}` is an odd function in the x-direction and -:math: `B_{zxx}` is an even function in the :math:`x` direction. +:math:`B_{zx}` is an odd function in the :math:`x` direction, and +:math:`B_{zxx}` is an even function in the :math:`x` direction. Polarization Test @@ -68,7 +68,7 @@ Note to see the notebook and examples go to the the effect of amplitude is also small (There is a comparison in "test-micrometertip_smallampapprox_vs_exactsol.ipynb" on small amplitude, however, the result is unchecked). -- CmeritSingleSpinESR_Approx is compared with CermitARP_SmallTip. +- CermitSingleSpinESR_Approx is compared with CermitARP_SmallTip. - SingleSpinESR - The approximation is tested against both SPAM and hangdown geometry of the analytical solution diff --git a/mrfmsim/component/__init__.py b/mrfmsim/component/__init__.py index 34b9c7d..defd808 100644 --- a/mrfmsim/component/__init__.py +++ b/mrfmsim/component/__init__.py @@ -1,5 +1,7 @@ from .base import ComponentBase from .magnet import SphereMagnet, RectangularMagnet +from .cylindermagnet import CylinderMagnetApprox from .cantilever import Cantilever from .grid import Grid from .sample import Sample + diff --git a/mrfmsim/component/cantilever.py b/mrfmsim/component/cantilever.py index 280925b..2d8c02c 100644 --- a/mrfmsim/component/cantilever.py +++ b/mrfmsim/component/cantilever.py @@ -10,8 +10,13 @@ class Cantilever(ComponentBase): """Cantilever object. - :param float k: spring constant [aN/nm] - :param float f: mechanical resonance frequency [Hz] + :param float k_c: spring constant [aN/nm] + :param float f_c: mechanical resonance frequency [Hz] + + :ivar float k2f_modulated: spring constant to frequency ratio for modulated + cantilever (lockin) [Hz.nm/aN] + :ivar float k2f: spring constant to frequency ratio for non-modulated + cantilever [Hz.nm/aN] """ k_c: float = field(metadata={"unit": "aN/nm", "format": ".3e"}) diff --git a/mrfmsim/component/cylindermagnet.py b/mrfmsim/component/cylindermagnet.py new file mode 100644 index 0000000..caeeda8 --- /dev/null +++ b/mrfmsim/component/cylindermagnet.py @@ -0,0 +1,565 @@ +import numpy as np +import numba as nb +from dataclasses import dataclass, field +from mrfmsim.component import ComponentBase + + +@dataclass +class CylinderMagnetApprox(ComponentBase): + """Cylinder magnet object approximated by Rectangular Magnets. + + :param float magnet_radius: cylinder magnet radius [nm] + :param float magnet_length: cylinder magnet length [nm] + :param tuple magnet_origin: the position of the magnet origin + :math:`(x, y, z)` [nm] + :param float mu0_Ms: saturation magnetization [mT] + """ + + magnet_radius: float = field(metadata={"unit": "nm", "format": ".1f"}) + magnet_length: float = field(metadata={"unit": "nm", "format": ".1f"}) + magnet_origin: tuple[float, float, float] = field( + metadata={"unit": "nm", "format": ".1f"} + ) + mu0_Ms: float = field(metadata={"unit": "mT"}) + + def __post_init__(self): + d = self.magnet_radius / 10 + + self._range = np.array( + [ + [ + self.magnet_origin[0] - 3 * d, + self.magnet_origin[0] + 3 * d, + self.magnet_origin[1] - 10 * d, + self.magnet_origin[1] + 10 * d, + self.magnet_origin[2] - self.magnet_length / 2, + self.magnet_origin[2] + self.magnet_length / 2, + ], + [ + self.magnet_origin[0] - 5 * d, + self.magnet_origin[0] - 3 * d, + self.magnet_origin[1] - 9 * d, + self.magnet_origin[1] + 9 * d, + self.magnet_origin[2] - self.magnet_length / 2, + self.magnet_origin[2] + self.magnet_length / 2, + ], + [ + self.magnet_origin[0] + 3 * d, + self.magnet_origin[0] + 5 * d, + self.magnet_origin[1] - 9 * d, + self.magnet_origin[1] + 9 * d, + self.magnet_origin[2] - self.magnet_length / 2, + self.magnet_origin[2] + self.magnet_length / 2, + ], + [ + self.magnet_origin[0] - 7 * d, + self.magnet_origin[0] - 5 * d, + self.magnet_origin[1] - 8 * d, + self.magnet_origin[1] + 8 * d, + self.magnet_origin[2] - self.magnet_length / 2, + self.magnet_origin[2] + self.magnet_length / 2, + ], + [ + self.magnet_origin[0] + 5 * d, + self.magnet_origin[0] + 7 * d, + self.magnet_origin[1] - 8 * d, + self.magnet_origin[1] + 8 * d, + self.magnet_origin[2] - self.magnet_length / 2, + self.magnet_origin[2] + self.magnet_length / 2, + ], + [ + self.magnet_origin[0] - 8 * d, + self.magnet_origin[0] - 7 * d, + self.magnet_origin[1] - 7 * d, + self.magnet_origin[1] + 7 * d, + self.magnet_origin[2] - self.magnet_length / 2, + self.magnet_origin[2] + self.magnet_length / 2, + ], + [ + self.magnet_origin[0] + 7 * d, + self.magnet_origin[0] + 8 * d, + self.magnet_origin[1] - 7 * d, + self.magnet_origin[1] + 7 * d, + self.magnet_origin[2] - self.magnet_length / 2, + self.magnet_origin[2] + self.magnet_length / 2, + ], + [ + self.magnet_origin[0] - 9 * d, + self.magnet_origin[0] - 8 * d, + self.magnet_origin[1] - 5 * d, + self.magnet_origin[1] + 5 * d, + self.magnet_origin[2] - self.magnet_length / 2, + self.magnet_origin[2] + self.magnet_length / 2, + ], + [ + self.magnet_origin[0] + 8 * d, + self.magnet_origin[0] + 9 * d, + self.magnet_origin[1] - 5 * d, + self.magnet_origin[1] + 5 * d, + self.magnet_origin[2] - self.magnet_length / 2, + self.magnet_origin[2] + self.magnet_length / 2, + ], + [ + self.magnet_origin[0] - 10 * d, + self.magnet_origin[0] - 9 * d, + self.magnet_origin[1] - 3 * d, + self.magnet_origin[1] + 3 * d, + self.magnet_origin[2] - self.magnet_length / 2, + self.magnet_origin[2] + self.magnet_length / 2, + ], + [ + self.magnet_origin[0] + 9 * d, + self.magnet_origin[0] + 10 * d, + self.magnet_origin[1] - 3 * d, + self.magnet_origin[1] + 3 * d, + self.magnet_origin[2] - self.magnet_length / 2, + self.magnet_origin[2] + self.magnet_length / 2, + ], + ] + ) + self._pre_term = self.mu0_Ms / (4 * np.pi) + + def Bz_method(self, x, y, z): + r"""Calculate magnetic field :math:`B_z` [mT]. + + Approximating Cylinder Magnet by 11 Rectangular Magnets. When viewed from the + vertical direction, we are using a row of rectangles to approximate a circle, + these rectangular blocks are arranged side by side. + + The magnetic field of each rectangular magnet is calculated following the + method described in Ravaud2009 [#]_. The magnet is set up so that the + :math:`x` and :math:`y` dimensions are centered about the zero point. + The translation in :math:`z` shifts the tip of the magnet in the + :math:`z`-direction to be the given distance from the surface. + Using the Coulombian model, assuming a uniform magnetization throughout + the volume of the magnet and modeling each face of the magnet as a + layer of continuous current density. The total field is found by + summing over the faces. + The magnetic field is given by: + + .. math:: + B_z = \dfrac{\mu_0 M_s}{4\pi} \sum_{i=1}^{2} + \sum_{j=1}^2 \sum_{k=1}^2(-1)^{i+j+k} + \arctan \left( \dfrac{(x - x_i)(y - y_i))}{(z - z_k)R} \right) + + Here :math:`(x,y,z)` are the coordinates for the location at which we + want to know the field; + The magnet spans from :math:`x_1` to :math:`x_2` in the :math:`x`-direction, + :math:`y_1` to :math:`y_2` in the :math:`y`-direction, and :math:`z_1` to + :math:`z_2` in the :math:`z`-direction; + + .. math:: + R = \sqrt{(x - x_i)^2 + (y - y_j)^2 + (z - z_k)^2} + + where :math:`\mu_0 M_s` is the magnet's saturation magnetization in mT. + + .. [#] Ravaud, R. and Lemarquand, G. "Magnetic field produced by a + parallelepipedic magnet of various and uniform polarization" , + *PIER*, **2009**, *98*, 207-219 + [`10.2528/PIER09091704 `__]. + + + :param float x: :math:`x` coordinate of sample grid [nm] + :param float y: :math:`y` coordinate of sample grid [nm] + :param float z: :math:`z` coordinate of sample grid [nm] + """ + dx11, dx12 = x - self._range[0][0], x - self._range[0][1] + dy11, dy12 = y - self._range[0][2], y - self._range[0][3] + dz11, dz12 = z - self._range[0][4], z - self._range[0][5] + + dx21, dx22 = x - self._range[1][0], x - self._range[1][1] + dy21, dy22 = y - self._range[1][2], y - self._range[1][3] + dz21, dz22 = z - self._range[1][4], z - self._range[1][5] + + dx31, dx32 = x - self._range[2][0], x - self._range[2][1] + dy31, dy32 = y - self._range[2][2], y - self._range[2][3] + dz31, dz32 = z - self._range[2][4], z - self._range[2][5] + + dx41, dx42 = x - self._range[3][0], x - self._range[3][1] + dy41, dy42 = y - self._range[3][2], y - self._range[3][3] + dz41, dz42 = z - self._range[3][4], z - self._range[3][5] + + dx51, dx52 = x - self._range[4][0], x - self._range[4][1] + dy51, dy52 = y - self._range[4][2], y - self._range[4][3] + dz51, dz52 = z - self._range[4][4], z - self._range[4][5] + + dx61, dx62 = x - self._range[5][0], x - self._range[5][1] + dy61, dy62 = y - self._range[5][2], y - self._range[5][3] + dz61, dz62 = z - self._range[5][4], z - self._range[5][5] + + dx71, dx72 = x - self._range[6][0], x - self._range[6][1] + dy71, dy72 = y - self._range[6][2], y - self._range[6][3] + dz71, dz72 = z - self._range[6][4], z - self._range[6][5] + + dx81, dx82 = x - self._range[7][0], x - self._range[7][1] + dy81, dy82 = y - self._range[7][2], y - self._range[7][3] + dz81, dz82 = z - self._range[7][4], z - self._range[7][5] + + dx91, dx92 = x - self._range[8][0], x - self._range[8][1] + dy91, dy92 = y - self._range[8][2], y - self._range[8][3] + dz91, dz92 = z - self._range[8][4], z - self._range[8][5] + + dx101, dx102 = x - self._range[9][0], x - self._range[9][1] + dy101, dy102 = y - self._range[9][2], y - self._range[9][3] + dz101, dz102 = z - self._range[9][4], z - self._range[9][5] + + dx111, dx112 = x - self._range[10][0], x - self._range[10][1] + dy111, dy112 = y - self._range[10][2], y - self._range[10][3] + dz111, dz112 = z - self._range[10][4], z - self._range[10][5] + + return self._pre_term * ( + self._bz(dx11, dx12, dy11, dy12, dz11, dz12) + + self._bz(dx21, dx22, dy21, dy22, dz21, dz22) + + self._bz(dx31, dx32, dy31, dy32, dz31, dz32) + + self._bz(dx41, dx42, dy41, dy42, dz41, dz42) + + self._bz(dx51, dx52, dy51, dy52, dz51, dz52) + + self._bz(dx61, dx62, dy61, dy62, dz61, dz62) + + self._bz(dx71, dx72, dy71, dy72, dz71, dz72) + + self._bz(dx81, dx82, dy81, dy82, dz81, dz82) + + self._bz(dx91, dx92, dy91, dy92, dz91, dz92) + + self._bz(dx101, dx102, dy101, dy102, dz101, dz102) + + self._bz(dx111, dx112, dy111, dy112, dz111, dz112) + ) + + @staticmethod + @nb.vectorize( + [ + nb.float64( + nb.float64, + nb.float64, + nb.float64, + nb.float64, + nb.float64, + nb.float64, + ) + ], + nopython=True, + target="parallel", + ) + def _bz(dx1, dx2, dy1, dy2, dz1, dz2): + """Calculate the summation term for magnetic field optimized by numba. + + See method Bz_method for the explanation. + + :param float dx1: distance between grid and one end of magnet + in :math:`x` direction [nm] + :param float dx2: distance between grid and other end of magnet + in :math:`x` direction [nm] + :param float dy1: distance between grid and one end of magnet + in :math:`y` direction [nm] + :param float dy2: distance between grid and other end of magnet + in :math:`y` direction [nm] + :param float dz1: distance between grid and one end of magnet + in :math:`z` direction [nm] + :param float dz2: distance between grid and other end of magnet + in :math:`z` direction [nm] + """ + + return ( + -np.arctan2(dx1 * dy1, (np.sqrt(dx1**2 + dy1**2 + dz1**2) * dz1)) + + np.arctan2(dx2 * dy1, (np.sqrt(dx2**2 + dy1**2 + dz1**2) * dz1)) + + np.arctan2(dx1 * dy2, (np.sqrt(dx1**2 + dy2**2 + dz1**2) * dz1)) + - np.arctan2(dx2 * dy2, (np.sqrt(dx2**2 + dy2**2 + dz1**2) * dz1)) + + np.arctan2(dx1 * dy1, (np.sqrt(dx1**2 + dy1**2 + dz2**2) * dz2)) + - np.arctan2(dx2 * dy1, (np.sqrt(dx2**2 + dy1**2 + dz2**2) * dz2)) + - np.arctan2(dx1 * dy2, (np.sqrt(dx1**2 + dy2**2 + dz2**2) * dz2)) + + np.arctan2(dx2 * dy2, (np.sqrt(dx2**2 + dy2**2 + dz2**2) * dz2)) + ) + + def Bzx_method(self, x, y, z): + r"""Calculate magnetic field gradient :math:`B_{zx}`. + + Approximating Cylinder Magnet by 11 Rectangular Magnets. When viewed from the + vertical direction, we are using a row of rectangles to approximate a circle, + these rectangular blocks are arranged side by side. + + The magnetic field gradient for RectangularMagnet is: + :math:`B_{zx} = \dfrac{\partial{B_z}}{\partial x}` is + given by the following: + + .. math:: + B_{zx} = \dfrac{\mu_0 M_s}{4 \pi} \sum_{i=1}^2 \sum_{j=1}^2 + \sum_{k=1}^2(-1)^{i+j+k} + \left( \dfrac{(y-y_j)(z-z_k)}{ R((x-x_i)^2 + (z-z_k)^2))} + \right) + + As described above, :math:`(x,y,z)` are coordinates for the location + at which we want to know the field gradient; the magnet spans from + x1 to x2 in the :math:`x`-direction, y1 to y2 in the :math:`y`-direction, and + from z1 to z2 in the :math:`z`-direction; + + .. math:: + R = \sqrt{(x - x_i)^2 + (y - y_j)^2 + (z - z_k)^2} + + :math:`\mu_0 M_s` is the magnet's saturation magnetization in mT. + + :param float x: :math:`x` coordinate [nm] + :param float y: :math:`y` coordinate [nm] + :param float z: :math:`z` coordinate [nm] + """ + + dx11, dx12 = x - self._range[0][0], x - self._range[0][1] + dy11, dy12 = y - self._range[0][2], y - self._range[0][3] + dz11, dz12 = z - self._range[0][4], z - self._range[0][5] + + dx21, dx22 = x - self._range[1][0], x - self._range[1][1] + dy21, dy22 = y - self._range[1][2], y - self._range[1][3] + dz21, dz22 = z - self._range[1][4], z - self._range[1][5] + + dx31, dx32 = x - self._range[2][0], x - self._range[2][1] + dy31, dy32 = y - self._range[2][2], y - self._range[2][3] + dz31, dz32 = z - self._range[2][4], z - self._range[2][5] + + dx41, dx42 = x - self._range[3][0], x - self._range[3][1] + dy41, dy42 = y - self._range[3][2], y - self._range[3][3] + dz41, dz42 = z - self._range[3][4], z - self._range[3][5] + + dx51, dx52 = x - self._range[4][0], x - self._range[4][1] + dy51, dy52 = y - self._range[4][2], y - self._range[4][3] + dz51, dz52 = z - self._range[4][4], z - self._range[4][5] + + dx61, dx62 = x - self._range[5][0], x - self._range[5][1] + dy61, dy62 = y - self._range[5][2], y - self._range[5][3] + dz61, dz62 = z - self._range[5][4], z - self._range[5][5] + + dx71, dx72 = x - self._range[6][0], x - self._range[6][1] + dy71, dy72 = y - self._range[6][2], y - self._range[6][3] + dz71, dz72 = z - self._range[6][4], z - self._range[6][5] + + dx81, dx82 = x - self._range[7][0], x - self._range[7][1] + dy81, dy82 = y - self._range[7][2], y - self._range[7][3] + dz81, dz82 = z - self._range[7][4], z - self._range[7][5] + + dx91, dx92 = x - self._range[8][0], x - self._range[8][1] + dy91, dy92 = y - self._range[8][2], y - self._range[8][3] + dz91, dz92 = z - self._range[8][4], z - self._range[8][5] + + dx101, dx102 = x - self._range[9][0], x - self._range[9][1] + dy101, dy102 = y - self._range[9][2], y - self._range[9][3] + dz101, dz102 = z - self._range[9][4], z - self._range[9][5] + + dx111, dx112 = x - self._range[10][0], x - self._range[10][1] + dy111, dy112 = y - self._range[10][2], y - self._range[10][3] + dz111, dz112 = z - self._range[10][4], z - self._range[10][5] + + return self._pre_term * ( + self._bzx(dx11, dx12, dy11, dy12, dz11, dz12) + + self._bzx(dx21, dx22, dy21, dy22, dz21, dz22) + + self._bzx(dx31, dx32, dy31, dy32, dz31, dz32) + + self._bzx(dx41, dx42, dy41, dy42, dz41, dz42) + + self._bzx(dx51, dx52, dy51, dy52, dz51, dz52) + + self._bzx(dx61, dx62, dy61, dy62, dz61, dz62) + + self._bzx(dx71, dx72, dy71, dy72, dz71, dz72) + + self._bzx(dx81, dx82, dy81, dy82, dz81, dz82) + + self._bzx(dx91, dx92, dy91, dy92, dz91, dz92) + + self._bzx(dx101, dx102, dy101, dy102, dz101, dz102) + + self._bzx(dx111, dx112, dy111, dy112, dz111, dz112) + ) + + @staticmethod + @nb.vectorize( + [ + nb.float64( + nb.float64, + nb.float64, + nb.float64, + nb.float64, + nb.float64, + nb.float64, + ) + ], + nopython=True, + target="parallel", + ) + def _bzx(dx1, dx2, dy1, dy2, dz1, dz2): + """Calculate the summation term for magnetic field gradient. + + Optimized with numba. See method ``Bzx_method`` for the explanation. + + :param float dx1: distance between grid and one end of magnet + in :math:`x` direction [nm] + :param float dx2: distance between grid and other end of magnet + in :math:`x` direction [nm] + :param float dy1: distance between grid and one end of magnet + in :math:`y` direction [nm] + :param float dy2: distance between grid and other end of magnet + in :math:`y` direction [nm] + :param float dz1: distance between grid and one end of magnet + in :math:`z` direction [nm] + :param float dz2: distance between grid and other end of magnet + in :math:`z` direction [nm] + """ + + return ( + -dy1 * dz1 / (np.sqrt(dx1**2 + dy1**2 + dz1**2) * (dx1**2 + dz1**2)) + + dy1 * dz1 / (np.sqrt(dx2**2 + dy1**2 + dz1**2) * (dx2**2 + dz1**2)) + + dy2 * dz1 / (np.sqrt(dx1**2 + dy2**2 + dz1**2) * (dx1**2 + dz1**2)) + - dy2 * dz1 / (np.sqrt(dx2**2 + dy2**2 + dz1**2) * (dx2**2 + dz1**2)) + + dy1 * dz2 / (np.sqrt(dx1**2 + dy1**2 + dz2**2) * (dx1**2 + dz2**2)) + - dy1 * dz2 / (np.sqrt(dx2**2 + dy1**2 + dz2**2) * (dx2**2 + dz2**2)) + - dy2 * dz2 / (np.sqrt(dx1**2 + dy2**2 + dz2**2) * (dx1**2 + dz2**2)) + + dy2 * dz2 / (np.sqrt(dx2**2 + dy2**2 + dz2**2) * (dx2**2 + dz2**2)) + ) + + def Bzxx_method(self, x, y, z): + r"""Calculate magnetic field second derivative :math:`B_{zxx}`. + + Approximating Cylinder Magnet by 11 Rectangular Magnets. When viewed from the + vertical direction, we are using a row of rectangles to approximate a circle, + these rectangular blocks are arranged side by side. + + The magnetic field second derivative for RectangularMagnet is: + :math:`B_{zxx} \equiv \partial^2 B_z / \partial x^2` + [ :math:`\mathrm{mT} \; \mathrm{nm}^{-2}`] + The magnetic field's second derivative is given by the following: + + .. math:: + B_{zxx} = \dfrac{\partial^2 B_z}{\partial x^2} + = \dfrac{\mu_0 M_s}{4 \pi} \sum_{i=1}^2 + \sum_{j=1}^2 \sum_{k=1}^2(-1)^{i+j+k} + \left( \dfrac{-(x-x_i)(y-y_j)(z-z_k) + (3(x-x_i)^2 +2(y-y_j)^2 + 3(z-z_k)^2)} + {((x-x_i)^2 + (y-y_j)^2 + (z-z_k)^2)^{3/2} + ((x-x_i)^2 + (z-z_k)^2)^2} \right) + + with the variables defined above. + + :param float x: :math:`x` coordinate [nm] + :param float y: :math:`y` coordinate [nm] + :param float z: :math:`z` coordinate [nm] + """ + + dx11, dx12 = x - self._range[0][0], x - self._range[0][1] + dy11, dy12 = y - self._range[0][2], y - self._range[0][3] + dz11, dz12 = z - self._range[0][4], z - self._range[0][5] + + dx21, dx22 = x - self._range[1][0], x - self._range[1][1] + dy21, dy22 = y - self._range[1][2], y - self._range[1][3] + dz21, dz22 = z - self._range[1][4], z - self._range[1][5] + + dx31, dx32 = x - self._range[2][0], x - self._range[2][1] + dy31, dy32 = y - self._range[2][2], y - self._range[2][3] + dz31, dz32 = z - self._range[2][4], z - self._range[2][5] + + dx41, dx42 = x - self._range[3][0], x - self._range[3][1] + dy41, dy42 = y - self._range[3][2], y - self._range[3][3] + dz41, dz42 = z - self._range[3][4], z - self._range[3][5] + + dx51, dx52 = x - self._range[4][0], x - self._range[4][1] + dy51, dy52 = y - self._range[4][2], y - self._range[4][3] + dz51, dz52 = z - self._range[4][4], z - self._range[4][5] + + dx61, dx62 = x - self._range[5][0], x - self._range[5][1] + dy61, dy62 = y - self._range[5][2], y - self._range[5][3] + dz61, dz62 = z - self._range[5][4], z - self._range[5][5] + + dx71, dx72 = x - self._range[6][0], x - self._range[6][1] + dy71, dy72 = y - self._range[6][2], y - self._range[6][3] + dz71, dz72 = z - self._range[6][4], z - self._range[6][5] + + dx81, dx82 = x - self._range[7][0], x - self._range[7][1] + dy81, dy82 = y - self._range[7][2], y - self._range[7][3] + dz81, dz82 = z - self._range[7][4], z - self._range[7][5] + + dx91, dx92 = x - self._range[8][0], x - self._range[8][1] + dy91, dy92 = y - self._range[8][2], y - self._range[8][3] + dz91, dz92 = z - self._range[8][4], z - self._range[8][5] + + dx101, dx102 = x - self._range[9][0], x - self._range[9][1] + dy101, dy102 = y - self._range[9][2], y - self._range[9][3] + dz101, dz102 = z - self._range[9][4], z - self._range[9][5] + + dx111, dx112 = x - self._range[10][0], x - self._range[10][1] + dy111, dy112 = y - self._range[10][2], y - self._range[10][3] + dz111, dz112 = z - self._range[10][4], z - self._range[10][5] + + return self._pre_term * ( + self._bzxx(dx11, dx12, dy11, dy12, dz11, dz12) + + self._bzxx(dx21, dx22, dy21, dy22, dz21, dz22) + + self._bzxx(dx31, dx32, dy31, dy32, dz31, dz32) + + self._bzxx(dx41, dx42, dy41, dy42, dz41, dz42) + + self._bzxx(dx51, dx52, dy51, dy52, dz51, dz52) + + self._bzxx(dx61, dx62, dy61, dy62, dz61, dz62) + + self._bzxx(dx71, dx72, dy71, dy72, dz71, dz72) + + self._bzxx(dx81, dx82, dy81, dy82, dz81, dz82) + + self._bzxx(dx91, dx92, dy91, dy92, dz91, dz92) + + self._bzxx(dx101, dx102, dy101, dy102, dz101, dz102) + + self._bzxx(dx111, dx112, dy111, dy112, dz111, dz112) + ) + + @staticmethod + @nb.vectorize( + [ + nb.float64( + nb.float64, + nb.float64, + nb.float64, + nb.float64, + nb.float64, + nb.float64, + ) + ], + nopython=True, + target="parallel", + ) + def _bzxx(dx1, dx2, dy1, dy2, dz1, dz2): + """The summation term for the second derivative of magnetic field. + + Optimized by numba. See Bzxx_method for the explanation. + + :param float dx1: distance between grid and one end of magnet + in :math:`x` direction [nm] + :param float dx2: distance between grid and other end of magnet + in :math:`x` direction [nm] + :param float dy1: distance between grid and one end of magnet + in :math:`y` direction [nm] + :param float dy2: distance between grid and other end of magnet + in :math:`y` direction [nm] + :param float dz1: distance between grid and one end of magnet + in :math:`z` direction [nm] + :param float dz2: distance between grid and other end of magnet + in :math:`z` direction [nm] + """ + + return ( + +dx1 + * dy1 + * dz1 + * (3.0 * dx1**2 + 2.0 * dy1**2 + 3.0 * dz1**2) + / ((dx1**2 + dy1**2 + dz1**2) ** 1.5 * (dx1**2 + dz1**2) ** 2) + - dx2 + * dy1 + * dz1 + * (3.0 * dx2**2 + 2.0 * dy1**2 + 3.0 * dz1**2) + / ((dx2**2 + dy1**2 + dz1**2) ** 1.5 * (dx2**2 + dz1**2) ** 2) + - dx1 + * dy2 + * dz1 + * (3.0 * dx1**2 + 2.0 * dy2**2 + 3.0 * dz1**2) + / ((dx1**2 + dy2**2 + dz1**2) ** 1.5 * (dx1**2 + dz1**2) ** 2) + + dx2 + * dy2 + * dz1 + * (3.0 * dx2**2 + 2.0 * dy2**2 + 3.0 * dz1**2) + / ((dx2**2 + dy2**2 + dz1**2) ** 1.5 * (dx2**2 + dz1**2) ** 2) + - dx1 + * dy1 + * dz2 + * (3.0 * dx1**2 + 2.0 * dy1**2 + 3.0 * dz2**2) + / ((dx1**2 + dy1**2 + dz2**2) ** 1.5 * (dx1**2 + dz2**2) ** 2) + + dx2 + * dy1 + * dz2 + * (3.0 * dx2**2 + 2.0 * dy1**2 + 3.0 * dz2**2) + / ((dx2**2 + dy1**2 + dz2**2) ** 1.5 * (dx2**2 + dz2**2) ** 2) + + dx1 + * dy2 + * dz2 + * (3.0 * dx1**2 + 2.0 * dy2**2 + 3.0 * dz2**2) + / ((dx1**2 + dy2**2 + dz2**2) ** 1.5 * (dx1**2 + dz2**2) ** 2) + - dx2 + * dy2 + * dz2 + * (3.0 * dx2**2 + 2.0 * dy2**2 + 3.0 * dz2**2) + / ((dx2**2 + dy2**2 + dz2**2) ** 1.5 * (dx2**2 + dz2**2) ** 2) + ) diff --git a/mrfmsim/component/grid.py b/mrfmsim/component/grid.py index e4ed68f..eb18d17 100644 --- a/mrfmsim/component/grid.py +++ b/mrfmsim/component/grid.py @@ -14,18 +14,18 @@ class Grid(ComponentBase): The grid array uses numpy's open mesh-grid, which has speed and storage benefits. - :param np.array length: array of lengths along (x, y, z) - :param np.array step: a list of step sizes - :param np.array origin: the grid origin - - :param tuple[int] shape: grid dimension + :param tuple[int, int, int] grid_shape: grid dimension (number of points in x, y, z direction) - :param list step: grid setup size in x, y, z direction - :param list origin: grid origin - :param float voxel: the volume of each grid voxel - :param np.array range: range in (x, y, z direction), shape (3, 2) - :param np.array length: actual lengths of the grid. This is recalculated + :param list[float] grid_step: grid step size in x, y, z direction [nm] + :param list[float] grid_origin: the grid origin [nm] + + :ivar ndarray grid_length: array of lengths along (x, y, z) + :ivar float grid_voxel: the volume of each grid voxel + :ivar ndarray grid_range: range in (x, y, z direction), shape (3, 2) + :ivar ndarray grid_length: actual lengths of the grid. This is recalculated based on the rounded value of grid shape and step size. + :ivar ndarray grid_extents: the grid extents in (x, y, z direction), shape (3, 2) + :ivar ndarray grid_array: the grid array in (x, y, z direction), shape (3, n) """ grid_shape: tuple[int] @@ -91,7 +91,7 @@ def extend_grid_by_length(self, ext_length): This is used to extend the grid by the cantilever motion. The length needs to be more than the step size to count. - :param int ext_pts: distance (one side) to extend along x direction + :param int ext_length: distance (one side) to extend along x direction (cantilever motion direction). The length should be a list of three dimensions. """ diff --git a/mrfmsim/component/magnet.py b/mrfmsim/component/magnet.py index 2b4e060..f227a96 100644 --- a/mrfmsim/component/magnet.py +++ b/mrfmsim/component/magnet.py @@ -8,8 +8,9 @@ class SphereMagnet(ComponentBase): """Spherical magnet object with its Bz, Bzx, Bzxx calculations. - :param float radius: sphere magnet radius [nm] - :param list origin: the position of the magnet origin (x, y, z) + :param float magnet_radius: sphere magnet radius [nm] + :param list magnet_origin: the position of the magnet origin + :math:`(x, y, z)` [nm] :param float mu0_Ms: saturation magnetization [mT] """ @@ -20,9 +21,9 @@ class SphereMagnet(ComponentBase): def Bz_method(self, x, y, z): r"""Calculate magnetic field :math:`B_z` [mT]. - :param float x: x coordinate of sample grid [nm] - :param float y: y coordinate of sample grid [nm] - :param float z: z coordinate of sample grid [nm] + :param float x: :math:`x` coordinate of sample grid [nm] + :param float y: :math:`y` coordinate of sample grid [nm] + :param float z: :math:`z` coordinate of sample grid [nm] The magnetic field is calculated as @@ -57,12 +58,15 @@ def Bz_method(self, x, y, z): def _bz(dx, dy, dz): """Internal calculation for bz, optimized with numba. - :param dx: normalized distances to the center of the magnet in x - :param dx: normalized distances to the center of the magnet in y - :param dx: normalized distances to the center of the magnet in z - :type: np.array + :param float dx: normalized distances to the center of the magnet in + :math:`x` direction [nm] + :param float dy: normalized distances to the center of the magnet in + :math:`y` direction [nm] + :param float dz: normalized distances to the center of the magnet in + :math:`z` direction [nm] + :return: bz without pre-term - :rtype: np.array + :rtype: float """ return ( @@ -80,14 +84,14 @@ def Bzx_method(self, x, y, z): gradient is calculated as .. math:: - B_{zx} = \dfrac{\partial B_z}{\partial z} + B_{zx} = \dfrac{\partial B_z}{\partial x} = \dfrac{\mu_0 M_s}{r} X \: \left( \dfrac{1}{R^5} - 5 \dfrac{Z^2}{R^7} \right) R = \sqrt{X^2+Y^2+Z^2} - :param float x: x coordinate of sample grid [nm] - :param float y: y coordinate of sample grid [nm] - :param float z: z coordinate of sample grid [nm] + :param float x: :math:`x` coordinate of sample grid [nm] + :param float y: :math:`y` coordinate of sample grid [nm] + :param float z: :math:`z` coordinate of sample grid [nm] :return: magnetic field gradient :rtype: np.array """ @@ -108,12 +112,15 @@ def Bzx_method(self, x, y, z): def _bzx(dx, dy, dz): """Internal calculation for bzx, optimized with numba. - :param dx: normalized distances to the center of the magnet in x - :param dy: normalized distances to the center of the magnet in y - :param dz: normalized distances to the center of the magnet in z - :type: np.array + :param float dx: normalized distances to the center of the magnet in + :math:`x` direction [nm] + :param float dy: normalized distances to the center of the magnet in + :math:`y` direction [nm] + :param float dz: normalized distances to the center of the magnet in + :math:`z` direction [nm] + :return: bzx without pre-term - :rtype: np.array + :rtype: float """ return dx * ( @@ -131,15 +138,15 @@ def Bzxx_method(self, x, y, z): second derivative is calculated as .. math:: - B_{zxx} = \dfrac{\partial B_z}{\partial z} + B_{zxx} = \dfrac{\partial^2 B_z}{\partial x^2} = \dfrac{\mu_0 M_s}{r^2} \: \left( \dfrac{1}{R^5} - 5 \dfrac{X^2}{R^7} - 5 \dfrac{Z^2}{R^7} + 35 \dfrac{X^2 Z^2}{R^9} \right) R = \sqrt{X^2+Y^2+Z^2} - :param float x: x coordinate of sample grid [nm] - :param float y: y coordinate of sample grid [nm] - :param float z: z coordinate of sample grid [nm] + :param float x: :math:`x` coordinate of sample grid [nm] + :param float y: :math:`y` coordinate of sample grid [nm] + :param float z: :math:`z` coordinate of sample grid [nm] :return: magnetic field second derivative :rtype: np.array """ @@ -160,12 +167,15 @@ def Bzxx_method(self, x, y, z): def _bzxx(dx, dy, dz): """Internal calculation for bzxx, optimized with numba. - :param dx: normalized distances to the center of the magnet in x - :param dy: normalized distances to the center of the magnet in y - :param dz: normalized distances to the center of the magnet in z - :type: np.array + :param float dx: normalized distances to the center of the magnet in + :math:`x` direction [nm] + :param float dy: normalized distances to the center of the magnet in + :math:`y` direction [nm] + :param float dz: normalized distances to the center of the magnet in + :math:`z` direction [nm] + :return: bzxx without pre-term - :rtype: np.array + :rtype: float """ return ( @@ -180,12 +190,11 @@ def _bzxx(dx, dy, dz): class RectangularMagnet(ComponentBase): """Rectangular magnet object with the bz, bzx, bzxx calculations. - :param list length: length of rectangular magnet in (x, y, z) direction [nm] - :param list origin: the position of the magnet origin (x, y, z) + :param list magnet_length: length of rectangular magnet in + :math:`(x, y, z)` direction [nm] + :param list magnet_origin: the position of the magnet origin + :math:`(x, y, z)` [nm] :param float mu0_Ms: saturation magnetization [mT] - :ivar np.array _range: - [-xrange, +xrange, -yrange, +yrange, -zrange, zrange] - range in the x, y and z direction after the center point to the origin. """ magnet_length: list[float] = field(metadata={"unit": "nm", "format": ".1f"}) @@ -204,11 +213,18 @@ def __post_init__(self): def Bz_method(self, x, y, z): r"""Calculate magnetic field :math:`B_z` [mT]. - The magnetic field is calculated following Ravaud2009. + The magnetic field is calculated following Ravaud2009 [#]_. + + The magnet is set up so that the + :math:`x` and :math:`y` dimensions are centered about the zero point. + The translation in :math:`z` shifts the tip of the magnet in the + :math:`z`-direction to be the given distance from the surface. + Using the Coulombian model, assuming a uniform magnetization throughout the volume of the magnet and modeling each face of the magnet as a layer of continuous current density. The total field is found by summing over the faces. + The magnetic field is given by: .. math:: @@ -218,29 +234,24 @@ def Bz_method(self, x, y, z): Here :math:`(x,y,z)` are the coordinates for the location at which we want to know the field; - The magnet spans from x1 to x2 in the ``x``-direction, - y1 to y2 in the ``y``-direction, and z1 to z2 in - the ``z``-direction; + The magnet spans from x1 to x2 in the :math:`x`-direction, + y1 to y2 in the :math:`y`-direction, and z1 to z2 in + the :math:`z`-direction; .. math:: R = \sqrt{(x - x_i)^2 + (y - y_j)^2 + (z - z_k)^2} where :math:`\mu_0 M_s` is the magnet's saturation magnetization in mT. - **Reference**: - Ravaud, R. and Lemarquand, G. "Magnetic field produced by a - parallelepipedic magnet of various and uniform polarization" , - *PIER*, **2009**, *98*, 207-219 - [`10.2528/PIER09091704 `__]. + .. [#] Ravaud, R. and Lemarquand, G. "Magnetic field produced by a + parallelepipedic magnet of various and uniform polarization" , + *PIER*, **2009**, *98*, 207-219 + [`10.2528/PIER09091704 `__]. - - set the magnet up so that the x and y dimensions are centered about - the zero point. - - The translation in z should shift the tip of the magnet in the - z-direction to be the given distance from the surface. - :param float x: x coordinate of sample grid [nm] - :param float y: y coordinate of sample grid [nm] - :param float z: z coordinate of sample grid [nm] + :param float x: :math:`x` coordinate of sample grid [nm] + :param float y: :math:`y` coordinate of sample grid [nm] + :param float z: :math:`z` coordinate of sample grid [nm] """ dx1, dx2 = x - self._range[0], x - self._range[1] @@ -267,13 +278,20 @@ def Bz_method(self, x, y, z): def _bz(dx1, dx2, dy1, dy2, dz1, dz2): """Calculate the summation term for magnetic field optimized by numba. - See method bz for the explanation. - :param float dx1, dx2: distance between grid and the one end of magnet - in x direction [nm] - :param float dy1, dy2: distance between grid and the one end of magnet - in y direction [nm] - :param float dz1, dz2: distance between grid and the one end of magnet - in z direction [nm] + See method Bz_method for the explanation. + + :param float dx1: distance between grid and one end of magnet in + :math:`x` direction [nm] + :param float dx2: distance between grid and other end of magnet in + :math:`x` direction [nm] + :param float dy1: distance between grid and one end of magnet in + :math:`y` direction [nm] + :param float dy2: distance between grid and other end of magnet in + :math:`y` direction [nm] + :param float dz1: distance between grid and one end of magnet in + :math:`z` direction [nm] + :param float dz2: distance between grid and other end of magnet in + :math:`z` direction [nm] """ return ( @@ -308,13 +326,13 @@ def Bzx_method(self, x, y, z): from z1 to z2 in the ``z``-direction; .. math:: - R = \sqrt{(x - x_i) + (y - y_j) + (z - z_k)} + R = \sqrt{(x - x_i)^2 + (y - y_j)^2 + (z - z_k)^2} :math:`\mu_0 M_s` is the magnet's saturation magnetization in mT. - :param float x: ``x`` coordinate [nm] - :param float y: ``y`` coordinate [nm] - :param float z: ``z`` coordinate [nm] + :param float x: :math:`x` coordinate [nm] + :param float y: :math:`y` coordinate [nm] + :param float z: :math:`z` coordinate [nm] """ dx1, dx2 = x - self._range[0], x - self._range[1] @@ -341,14 +359,20 @@ def Bzx_method(self, x, y, z): def _bzx(dx1, dx2, dy1, dy2, dz1, dz2): """Calculate the summation term for magnetic field gradient. - Optimized with numba. See method bzx for the explanation. - - :param np.array dx1, dx2: distance between grid and the 2 ends of - magnet in x direction [nm] - :param np.array dy1, dy2: distance between grid and the 2 ends of - magnet in y direction [nm] - :param np.array dz1, dz2: distance between grid and the 2 ends of - magnet in z direction [nm] + Optimized with numba. See method Bzx_method for the explanation. + + :param float dx1: distance between grid and one end of magnet in + :math:`x` direction [nm] + :param float dx2: distance between grid and other end of magnet in + :math:`x` direction [nm] + :param float dy1: distance between grid and one end of magnet in + :math:`y` direction [nm] + :param float dy2: distance between grid and other end of magnet in + :math:`y` direction [nm] + :param float dz1: distance between grid and one end of magnet in + :math:`z` direction [nm] + :param float dz2: distance between grid and other end of magnet in + :math:`z` direction [nm] """ return ( @@ -363,14 +387,14 @@ def _bzx(dx1, dx2, dy1, dy2, dz1, dz2): ) def Bzxx_method(self, x, y, z): - r"""Calculate magnetic field second derivative :math:`B_{zxx}`/ + r"""Calculate magnetic field second derivative :math:`B_{zxx}`. :math:`B_{zxx} \equiv \partial^2 B_z / \partial x^2` [ :math:`\mathrm{mT} \; \mathrm{nm}^{-2}`] The magnetic field's second derivative is given by the following: .. math:: - B_{zxx} = \dfrac{\partial B_z}{\partial z} + B_{zxx} = \dfrac{\partial^2 B_z}{\partial x^2} = \dfrac{\mu_0 M_s}{4 \pi} \sum_{i=1}^2 \sum_{j=1}^2 \sum_{k=1}^2(-1)^{i+j+k} \left( \dfrac{-(x-x_i)(y-y_j)(z-z_k) @@ -379,9 +403,10 @@ def Bzxx_method(self, x, y, z): ((x-x_i)^2 + (z-z_k)^2)^2} \right) with the variables defined above. - :param float x: ``x`` coordinate [nm] - :param float y: ``y`` coordinate [nm] - :param float z: ``z`` coordinate [nm] + + :param float x: :math:`x` coordinate [nm] + :param float y: :math:`y` coordinate [nm] + :param float z: :math:`z` coordinate [nm] """ dx1, dx2 = x - self._range[0], x - self._range[1] @@ -408,13 +433,20 @@ def Bzxx_method(self, x, y, z): def _bzxx(dx1, dx2, dy1, dy2, dz1, dz2): """The summation term for the second derivative of magnetic field. - Optimized by numba. See bzxx method for the explanation. - :param float dx1, dx2: distance between grid and the one end of magnet - in x direction [nm] - :param float dy1, dy2: distance between grid and the one end of magnet - in y direction [nm] - :param float dz1, dz2: distance between grid and the one end of magnet - in z direction [nm] + Optimized by numba. See Bzxx_method for the explanation. + + :param float dx1: distance between grid and one end of magnet in + :math:`x` direction [nm] + :param float dx2: distance between grid and other end of magnet in + :math:`x` direction [nm] + :param float dy1: distance between grid and one end of magnet in + :math:`y` direction [nm] + :param float dy2: distance between grid and other end of magnet in + :math:`y` direction [nm] + :param float dz1: distance between grid and one end of magnet in + :math:`z` direction [nm] + :param float dz2: distance between grid and other end of magnet in + :math:`z` direction [nm] """ return ( diff --git a/mrfmsim/component/sample.py b/mrfmsim/component/sample.py index 7063158..056695c 100644 --- a/mrfmsim/component/sample.py +++ b/mrfmsim/component/sample.py @@ -28,9 +28,9 @@ class Sample(ComponentBase): :param float temperature: the sample temperature [K] :param float spin_density: the sample spin density :math:`\rho` [1/nm^3] :param float Gamma: spin gyromagnetic ratio [rad/s.mT] - defaults to None if spin_type is one of the preset - :param float L: spin angular momentum [unitless] - default to None if spin_type is one of the preset + defaults to None if spin is one of the preset types + :param float J: spin angular momentum [unitless] + defaults to None if spin is one of the preset types :param float dB_hom: homogeneous linewidth [mT] :param float dB_sat: saturation linewidth [mT] """ diff --git a/mrfmsim/formula/field.py b/mrfmsim/formula/field.py index 5f8a966..ed38d65 100644 --- a/mrfmsim/formula/field.py +++ b/mrfmsim/formula/field.py @@ -1,4 +1,4 @@ -"""Cacluations related to the magnetic field.""" +"""Calculations related to the magnetic field.""" import numpy as np import numba as nb diff --git a/mrfmsim/formula/magnetization.py b/mrfmsim/formula/magnetization.py index 63b1b88..9062f00 100644 --- a/mrfmsim/formula/magnetization.py +++ b/mrfmsim/formula/magnetization.py @@ -9,13 +9,10 @@ def mz_eq(B_tot, Gamma, J, temperature): r"""Magnetization per spin at the thermal equilibrium using the Brillouin function. - :param float gamma: the gyromagnetic ratio. [rad/s.mT] - :param float j: total spin angular momentum + :param float B_tot: total magnetic field [mT] + :param float Gamma: the gyromagnetic ratio [rad/s.mT] + :param float J: total spin angular momentum :param float temperature: the spin temperature [K] - :param float spin_density: the sample spin - density :math:`\rho` [1/nm^3] - :param float b0: the external magnetic field [mT] - :param float bz: tip magnetic field in z [mT] :return: equilibrium per-spin magnetization [aN.nm/mT] The outputs are calculated from the sample properties @@ -82,9 +79,10 @@ def mz_eq(B_tot, Gamma, J, temperature): def mz2_eq(Gamma, J): r"""Compute the magnetization variance per spin. - :param float gamma: the gyromagnetic ratio. [rad/s.mT] + :param float Gamma: the gyromagnetic ratio [rad/s.mT] :param float J: total spin angular momentum - :return: Mz2_eq, rhoMz_eq + :return: magnetization variance per spin [aN^2.nm^2/mT^2] + :rtype: float """ mu = HBAR * Gamma * np.sqrt(J * (J + 1) / 3.0) # aN.nm/mT diff --git a/mrfmsim/formula/misc.py b/mrfmsim/formula/misc.py index 48b4a36..81d704f 100644 --- a/mrfmsim/formula/misc.py +++ b/mrfmsim/formula/misc.py @@ -9,7 +9,10 @@ def convert_grid_pts(distance, grid_step): """Convert distance to ext points. - :param int distance: distance in the x direction. + :param float distance: distance in the :math:`x` direction [nm] + :param list[float] grid_step: grid step size [nm] + :return: number of grid points + :rtype: int """ return math.floor(distance / grid_step[0]) @@ -76,16 +79,18 @@ def singlespin_analytical( \frac{4\hat{y}^3}{3\pi(1+\hat{y}^2)^2}[(1+\hat{y}^2)E(-1/\hat{y}^2) - (\hat{y}^2 -1)K(-1/\hat{y}^2) - :param float Gamma: the gyromagnetic ratio + :param float Gamma: the gyromagnetic ratio [rad/s.mT] + :param str geometry: experiment geometry ('spam' or 'hangdown') :param float J: the spin angular momentum - :param array magnet_origin: the origin of the magnet - :param float magnet_radius: the radius of the magnet - :param float mu0_Ms: the saturation magnetization of the magnet - :param float mu_z: the spin magnetic moment - :param float magnet_spin_distance: magnet-spin distance + :param float magnet_spin_dist: magnet-spin distance [nm] + :param list magnet_origin: the origin of the magnet [nm] + :param float magnet_radius: the radius of the magnet [nm] + :param float mu0_Ms: the saturation magnetization of the magnet [mT] + :param float x_0p: zero-to-peak amplitude of the cantilever [nm] - :return: The analytical solution for a single spin (effective force). + :return: The analytical solution for a single spin (effective force) [aN]. The spin constant shift is the effective force divided by the 0 to peak amplitude. + :rtype: float """ mu_z = HBAR * Gamma * J diff --git a/mrfmsim/formula/polarization.py b/mrfmsim/formula/polarization.py index bb8ea30..fe162fc 100644 --- a/mrfmsim/formula/polarization.py +++ b/mrfmsim/formula/polarization.py @@ -2,8 +2,7 @@ # -*- coding: utf-8 -*- # Peter Sun -"""Collection of calculations of relative changes in polarization. -""" +"""Collection of calculations of relative changes in polarization.""" import numba import numpy as np @@ -233,6 +232,10 @@ def rel_dpol_sat_td(Bzx, B1, ext_B_offset, ext_pts, Gamma, T2, tip_v): """Relative change in polarization for time-dependent saturation. The result is not a steady-state solution because it ignores T1 relaxation. + In the case where Bzx is 0, and B_offset is symmetric, the division will be nan. + Here we try to adjust the nan values to the average of the surrounding values. + However, if the resulting value is still nan or there are nan values at the boundary, + an ValueError is raised. """ # ignore division error the Exp takes care of the inf, and nan np.seterr(divide="ignore", invalid="ignore") @@ -244,10 +247,19 @@ def rel_dpol_sat_td(Bzx, B1, ext_B_offset, ext_pts, Gamma, T2, tip_v): div = np.divide(atan_omega_f - atan_omega_i, Bzx) - # adjust for the center slice of the discontinuous issue - center_index = div.shape[0] // 2 # if the grid is even it should not be a problem - if np.all(np.isnan(div[center_index])): - div[center_index] = (div[center_index + 1] + div[center_index - 1]) / 2 + # adjust the nan values to the average of the surrounding values + for idx in np.where(np.isnan(div))[0]: + if idx == 0 or idx == len(div) - 1: + raise ValueError( + "Nan values at the boundary, check the Bzx and B_offset values." + ) + value = (div[idx + 1] + div[idx - 1]) / 2 + + if np.isnan(value): + raise ValueError( + "Nan value from division, check the Bzx and B_offset values." + ) + div[idx] = value rt = Gamma * B1**2 * np.abs(div) / tip_v dpol = np.exp(-rt) diff --git a/pyproject.toml b/pyproject.toml index fa54b2d..3d7f4b0 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -22,7 +22,7 @@ build-backend = "poetry.core.masonry.api" [tool.poetry.dependencies] python = ">=3.10" numpy = ">=1.26.4" -mmodel = { git = "https://www.github.com/Marohn-Group/mmodel.git", branch = "develop-0.8.1" } +mmodel = ">=0.9.0" numba = ">=0.60.0" scipy = ">=1.14.1" tox = { version = ">=3.25.0", optional = true } diff --git a/tests/test_component/test_cylindermagnet.py b/tests/test_component/test_cylindermagnet.py new file mode 100644 index 0000000..2874528 --- /dev/null +++ b/tests/test_component/test_cylindermagnet.py @@ -0,0 +1,140 @@ +import numpy as np +from mrfmsim.component import CylinderMagnetApprox + + +"""Test CylinderMagnet module in mrfmsim.component. + + +CylinderMagnetApprox +-------------------- + +All cylinder magnets are tested using the parameters: +radius = 0.5 nm +length = 10 nm +mu0_Ms = 1 mT +origin = [0, 0, 0] nm + +Bz_method +^^^^^^^^^^ + +1. Test Bz output has the same shape as the input grid +2. Test Bz output is symmetric along x-axis +3. Test Bz in the near field is close to the exact value +4. Test Bz in the far field is close to the exact value + +Bzx_method +^^^^^^^^^^ + +1. Test Bzx output has the same shape as the input grid +2. Test Bzx output is symmetric along x-axis +3. Test Bzx in the near field is close to the exact value +4. Test Bzx in the far field is close to the exact value +5. Test Bzx when x is 0, Bzx is 0 + +Bzxx_method +^^^^^^^^^^^ + +1. Test Bzxx output has the same shape as the input grid +2. Test Bzxx output is symmetric along x-axis +3. Test Bzxx in the near field is close to the exact value +4. Test Bzxx in the far field is close to the exact value +""" + + +class TestCylinderMagnet: + def setup_method(self): + self.magnet = CylinderMagnetApprox( + magnet_radius=0.5, magnet_length=10, magnet_origin=[0, 0, 0], mu0_Ms=1 + ) + + def test_Bz_shape(self): + """Test the shape of the magnetic field Bz to be the same as the input grid.""" + grid = np.ogrid[-1:1:2j, -1:1:5j, -1:1:10j] + Bz = self.magnet.Bz_method(grid[0], grid[1], grid[2]) + assert Bz.shape == (2, 5, 10) + + def test_Bzx_shape(self): + """Test the shape of the magnetic field Bzx to be the same as the input grid.""" + grid = np.ogrid[-1:1:2j, -1:1:5j, -1:1:10j] + Bzx = self.magnet.Bzx_method(grid[0], grid[1], grid[2]) + assert Bzx.shape == (2, 5, 10) + + def test_Bzxx_shape(self): + """Test the shape of the magnetic field Bzxx to be the same as the input grid.""" + grid = np.ogrid[-1:1:2j, -1:1:5j, -1:1:10j] + Bzxx = self.magnet.Bzxx_method(grid[0], grid[1], grid[2]) + assert Bzxx.shape == (2, 5, 10) + + def test_Bz_symmetry(self): + """Test the symmetry of the magnetic field Bz along x-axis.""" + Bz1 = self.magnet.Bz_method(1, 2, 6) + Bz2 = self.magnet.Bz_method(-1, 2, 6) + assert np.allclose(Bz1, Bz2) + + def test_Bzx_symmetry(self): + """Test the symmetry of the magnetic field Bzx along x-axis.""" + Bzx1 = self.magnet.Bzx_method(1, 2, 6) + Bzx2 = self.magnet.Bzx_method(-1, 2, 6) + assert np.allclose(Bzx1, -Bzx2) # Bzx is antisymmetric + + def test_Bzxx_symmetry(self): + """Test the symmetry of the magnetic field Bzxx along x-axis.""" + Bzxx1 = self.magnet.Bzxx_method(1, 2, 6) + Bzxx2 = self.magnet.Bzxx_method(-1, 2, 6) + assert np.allclose(Bzxx1, Bzxx2) # Bzxx is symmetric + + def test_Bz_near_field(self): + """Test the magnetic field Bz at near field. + + Test at poles (1.0, 2.0, 6.0). + """ + Bz = self.magnet.Bz_method(1, 2, 6) + assert np.allclose(Bz, 0.003938311762612047, atol=1e-5) + + def test_Bzx_near_field(self): + """Test the magnetic field Bzx at near field. + + Test at poles (1.0, 2.0, 6.0). + """ + Bzx = self.magnet.Bzx_method(1, 2, 6) + assert np.allclose(Bzx, -0.0022319789019074046, atol=1e-5) + + def test_Bzxx_near_field(self): + """Test the magnetic field Bzxx at near field. + + Test at poles (1.0, 2.0, 6.0). + """ + Bzxx = self.magnet.Bzxx_method(1, 2, 6) + assert np.allclose(Bzxx, -0.00035144076191428254, atol=1e-5) + + def test_Bz_far_field(self): + """Test the magnetic field Bz at far field. + + Test at (10.0, 0.0, 0.0). + """ + Bz = self.magnet.Bz_method(10, 0, 0) + assert np.allclose(Bz, (-0.00045051510382442363), atol=1e-4) + + def test_Bzx_far_field(self): + """Test the magnetic field Bzx at far field. + + Test at (10.0, 0.0, 0.0). + """ + Bzx = self.magnet.Bzx_method(10, 0, 0) + assert np.allclose(Bzx, 0.00010817803386829378, atol=1e-4) + + def test_Bzxx_far_field(self): + """Test the magnetic field Bzxx at far field. + + Test at (10.0, 0.0, 0.0). + """ + Bzxx = self.magnet.Bzxx_method(10, 0, 0) + assert np.allclose(Bzxx, -3.245766094831442e-05, atol=1e-4) + + def test_Bzx_x_is_0(self): + """Test the magnetic field Bzx when x is 0. + + Test at (0.0, 10.0, 0.0). + """ + Bzx = self.magnet.Bzx_method(0, 10, 0) + assert np.allclose(Bzx, 0, atol=1e-10) diff --git a/tests/test_formula/test_polarization.py b/tests/test_formula/test_polarization.py index e382463..7307a1b 100644 --- a/tests/test_formula/test_polarization.py +++ b/tests/test_formula/test_polarization.py @@ -178,6 +178,47 @@ def test_rel_dpol_sat_td_symmetry(sample_e): assert np.array_equal(rpol_a, rpol_b) +def test_rel_dpol_sat_td_nan(sample_e): + """Test rel_dpol_sat_td does not return nan. + + When Bzx is 0 and B_offset is symmetric, the result should not be nan. + """ + Bzx = np.array([1, 0, -1]) + ext_B_offset_a = np.array([2, 0, 0, 0, 2]) + ext_B_offset_b = np.array([0, 2, 2, 2, 0]) + + rpol_a = pol.rel_dpol_sat_td( + Bzx, 1.0, ext_B_offset_a, 1, sample_e.Gamma, sample_e.T2, 2000 + ) + rpol_b = pol.rel_dpol_sat_td( + Bzx, 1.0, ext_B_offset_b, 1, sample_e.Gamma, sample_e.T2, 2000 + ) + assert not np.any(np.isnan(rpol_a)) + assert not np.any(np.isnan(rpol_b)) + + +def test_rel_dpol_sat_td_nan_boundary(sample_e): + """Test rel_dpol_sat_td raises an error if nan values are at the boundary.""" + Bzx = np.array([0, 0, -1]) + ext_B_offset = np.array([0, 0, 0, 0, 2]) + + with pytest.raises(ValueError, match="Nan values at the boundary"): + pol.rel_dpol_sat_td( + Bzx, 1.0, ext_B_offset, 1, sample_e.Gamma, sample_e.T2, 2000 + ) + + +def test_rel_dpol_sat_td_nan_division(sample_e): + """Test rel_dpol_sat_td raises an error if nan values are from division.""" + Bzx = np.array([2, 0, 0, -1]) + ext_B_offset = np.array([1, 0, 0, 0, 0, 1]) + + with pytest.raises(ValueError, match="Nan value from division"): + pol.rel_dpol_sat_td( + Bzx, 1.0, ext_B_offset, 1, sample_e.Gamma, sample_e.T2, 2000 + ) + + def test_rel_dpol_sat_td_without_td(sample_e): """Test rel_dpol_sat_td completely saturate spins if no td component.