diff --git a/docs/introduction/fitting_example.rst b/docs/introduction/fitting_example.rst index ff98cc870..401b9da45 100644 --- a/docs/introduction/fitting_example.rst +++ b/docs/introduction/fitting_example.rst @@ -339,7 +339,7 @@ to be modeled with the Method class. Al27 = Isotope(symbol="27Al") origin_offset = spectral_dims[0]["origin_offset"] - B0 = Al27.ref_freq_to_B0(origin_offset / 1e6) + B0 = Al27.ref_freq_to_B0(origin_offset) MAS = BlochDecayCTSpectrum( channels=["27Al"], diff --git a/docs/introduction/getting_started.rst b/docs/introduction/getting_started.rst index 16f53d757..8219a03a5 100644 --- a/docs/introduction/getting_started.rst +++ b/docs/introduction/getting_started.rst @@ -137,8 +137,8 @@ familiar to an NMR spectroscopist. from mrsimulator.spin_system.isotope import Isotope # Set the magnetic flux density in T from the proton - # frequency of TMS, a primary reference, in MHz - B0 = Isotope(symbol="1H").ref_freq_to_B0(400) + # frequency of TMS, a primary reference, in Hz + B0 = Isotope(symbol="1H").ref_freq_to_B0(400e6) # Create a BlochDecaySpectrum instance method = BlochDecaySpectrum( diff --git a/src/mrsimulator/simulator/__init__.py b/src/mrsimulator/simulator/__init__.py index f349761d1..7aec3b563 100644 --- a/src/mrsimulator/simulator/__init__.py +++ b/src/mrsimulator/simulator/__init__.py @@ -420,7 +420,7 @@ def run( )(jobs) B0 = method.spectral_dimensions[0].events[0].magnetic_flux_density - w_ref = method.channels[0].B0_to_ref_freq(B0) * 1e6 + w_ref = method.channels[0].B0_to_ref_freq(B0) for seq in method.spectral_dimensions: seq.origin_offset = w_ref diff --git a/src/mrsimulator/spin_system/isotope.py b/src/mrsimulator/spin_system/isotope.py index 740b2c86e..4528aeeae 100644 --- a/src/mrsimulator/spin_system/isotope.py +++ b/src/mrsimulator/spin_system/isotope.py @@ -218,7 +218,7 @@ def larmor_freq(self, B0=9.4): float B0: magnetic field strength in T Returns: - float: Larmor frequency in MHz + float: Larmor frequency in Hz Example ------- @@ -226,13 +226,13 @@ def larmor_freq(self, B0=9.4): >>> silicon = Isotope(symbol="29Si") >>> freq = silicon.larmor_freq(B0 = 9.4) """ - return -self.gyromagnetic_ratio * B0 + return -self.gyromagnetic_ratio * B0 * 1.0e6 - def ref_freq_to_B0(self, ref_freq=400): + def ref_freq_to_B0(self, ref_freq=400e6): """Return the magnetic field strength B0 given the primary reference frequency. Args: - float ref_freq: primary reference frequency in MHz + float ref_freq: primary reference frequency in Hz Returns: float: magnetic flux density in T @@ -241,10 +241,10 @@ def ref_freq_to_B0(self, ref_freq=400): ------- >>> H1 = Isotope(symbol="1H") - >>> B0 = H1.ref_freq_to_B0(ref_freq = 400) + >>> B0 = H1.ref_freq_to_B0(ref_freq = 400e6) """ ref_ratio = self.reference.ratio / 100 # normalize reference ratio to 1 - return 0.02348731439404777 * ref_freq / ref_ratio + return 1.0e-6 * 0.02348731439404777 * ref_freq / ref_ratio def B0_to_ref_freq(self, B0=9.4): """Return the primary reference frequency given the magnetic field strength B0. @@ -253,7 +253,7 @@ def B0_to_ref_freq(self, B0=9.4): float B0: magnetic flux density in T Returns: - float: primary reference frequency in MHz + float: primary reference frequency in Hz Example ------- @@ -262,7 +262,7 @@ def B0_to_ref_freq(self, B0=9.4): >>> B0 = H1.B0_to_ref_freq(B0 = 9.4) """ ref_ratio = self.reference.ratio / 100 # normalize reference ratio to 1 - return B0 * ref_ratio / 0.02348731439404777 + return 1.0e6 * B0 * ref_ratio / 0.02348731439404777 @property def ref_larmor_ratio(self): diff --git a/src/mrsimulator/spin_system/tests/test_isotope.py b/src/mrsimulator/spin_system/tests/test_isotope.py index 4e4970745..bbd832245 100644 --- a/src/mrsimulator/spin_system/tests/test_isotope.py +++ b/src/mrsimulator/spin_system/tests/test_isotope.py @@ -1,3 +1,4 @@ +import numpy as np import pytest from mrsimulator.spin_system.isotope import get_isotope_dict from mrsimulator.spin_system.isotope import Isotope @@ -16,9 +17,10 @@ def test_isotope(): assert silicon.quadrupole_moment == 0.0 assert silicon.spin == 0.5 assert silicon.efg_to_Cq == 0 - assert silicon.larmor_freq(B0=11.75) == 99.46962016339306 - assert silicon.B0_to_ref_freq(B0=11.75) == 99.3895867929281 - assert silicon.ref_freq_to_B0(ref_freq=99.3895867929281) == 11.75 + assert np.allclose(silicon.larmor_freq(B0=11.75), 99.46962016339306e6) + silicon_ref_freq = 99.3895867929281e6 + assert np.allclose(silicon.B0_to_ref_freq(B0=11.75), silicon_ref_freq) + assert np.allclose(silicon.ref_freq_to_B0(ref_freq=silicon_ref_freq), 11.75) proton = Isotope(symbol="1H") assert proton.atomic_number == 1 @@ -27,14 +29,15 @@ def test_isotope(): assert proton.quadrupole_moment == 0.0 assert proton.spin == 0.5 assert silicon.efg_to_Cq == 0 - assert proton.larmor_freq(B0=9.40) == -400.2283045725638 - assert proton.B0_to_ref_freq(B0=9.40) == 400.21604182989006 - assert proton.ref_freq_to_B0(ref_freq=400.21604182989006) == 9.40 + assert np.allclose(proton.larmor_freq(B0=9.40), -400.2283045725638e6) + proton_ref_freq = 400.21604182989006e6 + assert np.allclose(proton.B0_to_ref_freq(B0=9.40), proton_ref_freq) + assert np.allclose(proton.ref_freq_to_B0(ref_freq=proton_ref_freq), 9.40) americium = Isotope(symbol="243Am") assert americium.atomic_number == 95 assert americium.reference.ratio == 10.956895712592464 - assert americium.B0_to_ref_freq(B0=11.75) == 54.81406791045811 + assert np.allclose(americium.B0_to_ref_freq(B0=11.75), 54.81406791045811e6) nitrogen = Isotope(symbol="14N") assert nitrogen.atomic_number == 7 @@ -43,9 +46,10 @@ def test_isotope(): assert nitrogen.quadrupole_moment == 0.0193 assert nitrogen.spin == 1 assert nitrogen.efg_to_Cq == 4534820.208433115 - assert nitrogen.larmor_freq(B0=18.79) == -57.830093199115574 - assert nitrogen.B0_to_ref_freq(B0=18.79) == 57.81099284148487 - assert nitrogen.ref_freq_to_B0(ref_freq=57.81099284148487) == 18.79 + assert np.allclose(nitrogen.larmor_freq(B0=18.79), -57.830093199115574e6) + nitrogen_ref_freq = 57.81099284148487e6 + assert np.allclose(nitrogen.B0_to_ref_freq(B0=18.79), nitrogen_ref_freq) + assert np.allclose(nitrogen.ref_freq_to_B0(ref_freq=nitrogen_ref_freq), 18.79) error = "Isotope symbol `x` not recognized." with pytest.raises(Exception, match=error): diff --git a/src/mrsimulator/utils/spectral_fitting.py b/src/mrsimulator/utils/spectral_fitting.py index d23923997..63d136a1d 100644 --- a/src/mrsimulator/utils/spectral_fitting.py +++ b/src/mrsimulator/utils/spectral_fitting.py @@ -614,7 +614,7 @@ def _generate_distribution_spectrum( """ method = kernel.method isotope = method.channels[0] - larmor_freq = isotope.larmor_freq(B0=method.magnetic_flux_density) + larmor_freq = isotope.larmor_freq(B0=method.magnetic_flux_density) / 1.0e6 # MHz larmor_freq *= isotope.ref_larmor_ratio exp_spectrum = method.experiment diff --git a/tests/2D_spectrum_tests/test_das.py b/tests/2D_spectrum_tests/test_das.py index 1fc0ed71b..0b63a215f 100644 --- a/tests/2D_spectrum_tests/test_das.py +++ b/tests/2D_spectrum_tests/test_das.py @@ -137,7 +137,7 @@ def test_DAS(): eta = site.quadrupolar.eta iso = site.isotropic_chemical_shift iso_obs = quad_iso_shift( - 0.5, Cq, eta, spin, iso, larmor_freq, O17.B0_to_ref_freq(B0) + 0.5, Cq, eta, spin, iso, larmor_freq, O17.B0_to_ref_freq(B0) / 1e6 ) # get the index where there is a signal id1 = np.real(dataset_das[i] / dataset_das[i].max()) diff --git a/tests/2D_spectrum_tests/test_mqmas.py b/tests/2D_spectrum_tests/test_mqmas.py index a957ab5e8..b8ac3ff29 100644 --- a/tests/2D_spectrum_tests/test_mqmas.py +++ b/tests/2D_spectrum_tests/test_mqmas.py @@ -131,7 +131,7 @@ def test_ThreeQ_VAS_spin_3halves(): # ref: D. Massiot et al. / Solid State Nuclear Magnetic Resonance 6 (1996) 73-83 spin = method.channels[0].spin v0 = method.channels[0].gyromagnetic_ratio * B0 * 1e6 - vref = method.channels[0].B0_to_ref_freq(B0) * 1e6 + vref = method.channels[0].B0_to_ref_freq(B0) vq = (3 * 3.5e6) / (2 * spin * (2 * spin - 1)) v_iso = -9 * 17 / 8 + 1e6 / 8 * (vq**2 / abs(v0 * vref)) * ((0.36**2) / 3 + 1) @@ -196,7 +196,7 @@ def test_MQMAS_spin_5halves(): # ref: D. Massiot et al. / Solid State Nuclear Magnetic Resonance 6 (1996) 73-83 spin = method.channels[0].spin v0 = method.channels[0].gyromagnetic_ratio * 7 * 1e6 - vref = method.channels[0].B0_to_ref_freq(7) * 1e6 + vref = method.channels[0].B0_to_ref_freq(7) vq = 3 * 3.22e6 / (2 * spin * (2 * spin - 1)) v_iso = -(17 / 31) * 64.5 - (8e6 / 93) * (vq**2 / abs(v0 * vref)) * ( (0.66**2) / 3 + 1