From 549e133e558af13a1378530b0f7281f6efd40733 Mon Sep 17 00:00:00 2001 From: johannes-moegerle Date: Tue, 10 Mar 2026 09:38:50 +0100 Subject: [PATCH 01/10] bugfix calc nu angular --- src/rydstate/species/sqdt/species_object_sqdt.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/rydstate/species/sqdt/species_object_sqdt.py b/src/rydstate/species/sqdt/species_object_sqdt.py index d41d04c..75532e9 100644 --- a/src/rydstate/species/sqdt/species_object_sqdt.py +++ b/src/rydstate/species/sqdt/species_object_sqdt.py @@ -230,7 +230,7 @@ def calc_nu( if not isinstance(angular_ket, AngularKetLS): raise NotImplementedError("calc_nu is only implemented for AngularKetLS.") - l, j_tot, s_tot = angular_ket.l_r, angular_ket.j_tot, angular_ket.s_r + l, j_tot, s_tot = angular_ket.l_r, angular_ket.j_tot, angular_ket.s_tot if n <= nist_n_max and use_nist_data: # try to use NIST data if (n, l, j_tot, s_tot) in self._nist_energy_levels: From 1690816b6d08eee3b5aa3448198445d2dc50ad81 Mon Sep 17 00:00:00 2001 From: johannes-moegerle Date: Mon, 9 Mar 2026 14:25:17 +0100 Subject: [PATCH 02/10] add NotSet for m quantum number --- src/rydstate/angular/__init__.py | 2 ++ src/rydstate/angular/angular_ket.py | 26 ++++++++++++++------------ src/rydstate/angular/angular_state.py | 5 +++-- src/rydstate/angular/utils.py | 16 +++++++++++++++- src/rydstate/rydberg/rydberg_sqdt.py | 11 ++++++----- 5 files changed, 40 insertions(+), 20 deletions(-) diff --git a/src/rydstate/angular/__init__.py b/src/rydstate/angular/__init__.py index 66f0142..6d9c0ed 100644 --- a/src/rydstate/angular/__init__.py +++ b/src/rydstate/angular/__init__.py @@ -1,11 +1,13 @@ from rydstate.angular import utils from rydstate.angular.angular_ket import AngularKetFJ, AngularKetJJ, AngularKetLS from rydstate.angular.angular_state import AngularState +from rydstate.angular.utils import NotSet __all__ = [ "AngularKetFJ", "AngularKetJJ", "AngularKetLS", "AngularState", + "NotSet", "utils", ] diff --git a/src/rydstate/angular/angular_ket.py b/src/rydstate/angular/angular_ket.py index ccc0b71..fbedbc8 100644 --- a/src/rydstate/angular/angular_ket.py +++ b/src/rydstate/angular/angular_ket.py @@ -14,6 +14,7 @@ ) from rydstate.angular.utils import ( InvalidQuantumNumbersError, + NotSet, check_spin_addition_rule, get_possible_quantum_number_values, minus_one_pow, @@ -67,9 +68,9 @@ class AngularKetBase(ABC): f_tot: float """Total atom angular quantum number (including nuclear, core electron and rydberg electron contributions).""" - m: float | None + m: float | NotSet """Magnetic quantum number, which is the projection of `f_tot` onto the quantization axis. - If None, only reduced matrix elements can be calculated. + If NotSet, only reduced matrix elements can be calculated. """ def __init__( @@ -80,7 +81,7 @@ def __init__( s_r: float = 0.5, l_r: int | None = None, f_tot: float | None = None, # noqa: ARG002 - m: float | None = None, + m: float | NotSet = NotSet, species: str | SpeciesObject | None = None, ) -> None: """Initialize the Spin ket. @@ -114,7 +115,7 @@ def __init__( self.l_r = int(l_r) # f_tot will be set in the subclasses - self.m = None if m is None else float(m) + self.m = NotSet if isinstance(m, NotSet) else float(m) def _post_init(self) -> None: self.quantum_numbers = tuple(getattr(self, qn) for qn in self.quantum_number_names) @@ -132,7 +133,7 @@ def sanity_check(self, msgs: list[str] | None = None) -> None: if self.s_r != 0.5: msgs.append(f"Rydberg electron spin s_r must be 1/2, but {self.s_r=}") - if self.m is not None and not -self.f_tot <= self.m <= self.f_tot: + if not isinstance(self.m, NotSet) and not -self.f_tot <= self.m <= self.f_tot: msgs.append(f"m must be between -f_tot and f_tot, but {self.f_tot=}, {self.m=}") if msgs: @@ -149,7 +150,7 @@ def __setattr__(self, key: str, value: object) -> None: def __repr__(self) -> str: args = ", ".join(f"{qn}={val}" for qn, val in zip(self.quantum_number_names, self.quantum_numbers, strict=True)) - if self.m is not None: + if not isinstance(self.m, NotSet): args += f", m={self.m}" return f"{self.__class__.__name__}({args})" @@ -466,15 +467,16 @@ def calc_matrix_element(self, other: AngularKetBase, operator: AngularOperatorTy The dimensionless angular matrix element. """ - if self.m is None or other.m is None: - raise ValueError("m must be set to calculate the matrix element.") + if isinstance(self.m, NotSet) or isinstance(other.m, NotSet): + raise RuntimeError("m must be set to calculate the matrix element.") # noqa: TRY004 prefactor = self._calc_wigner_eckart_prefactor(other, kappa, q) reduced_matrix_element = self.calc_reduced_matrix_element(other, operator, kappa) return prefactor * reduced_matrix_element def _calc_wigner_eckart_prefactor(self, other: AngularKetBase, kappa: int, q: int) -> float: - assert self.m is not None and other.m is not None, "m must be set to calculate the Wigner-Eckart prefactor." # noqa: PT018 + if isinstance(self.m, NotSet) or isinstance(other.m, NotSet): + raise RuntimeError("m must be set to calculate the Wigner-Eckart prefactor.") # noqa: TRY004 return minus_one_pow(self.f_tot - self.m) * calc_wigner_3j(self.f_tot, kappa, other.f_tot, -self.m, q, other.m) def _kronecker_delta_non_involved_spins(self, other: AngularKetBase, qn: AngularMomentumQuantumNumbers) -> int: @@ -572,7 +574,7 @@ def __init__( l_tot: int | None = None, j_tot: float | None = None, f_tot: float | None = None, - m: float | None = None, + m: float | NotSet = NotSet, species: str | SpeciesObject | None = None, ) -> None: """Initialize the Spin ket.""" @@ -635,7 +637,7 @@ def __init__( j_r: float | None = None, j_tot: float | None = None, f_tot: float | None = None, - m: float | None = None, + m: float | NotSet = NotSet, species: str | SpeciesObject | None = None, ) -> None: """Initialize the Spin ket.""" @@ -698,7 +700,7 @@ def __init__( f_c: float | None = None, j_r: float | None = None, f_tot: float | None = None, - m: float | None = None, + m: float | NotSet = NotSet, species: str | SpeciesObject | None = None, ) -> None: """Initialize the Spin ket.""" diff --git a/src/rydstate/angular/angular_state.py b/src/rydstate/angular/angular_state.py index 5ec8cfc..9d996e0 100644 --- a/src/rydstate/angular/angular_state.py +++ b/src/rydstate/angular/angular_state.py @@ -13,6 +13,7 @@ AngularKetLS, ) from rydstate.angular.angular_matrix_element import is_angular_momentum_quantum_number +from rydstate.angular.utils import NotSet if TYPE_CHECKING: from collections.abc import Iterator, Sequence @@ -209,8 +210,8 @@ def calc_matrix_element( "Different m values are not supported yet for AngularState.calc_matrix_element." ) - if self.kets[0].m is None or other.kets[0].m is None: - raise ValueError("m must be set for all kets to calculate the matrix element.") + if isinstance(self.kets[0].m, NotSet) or isinstance(other.kets[0].m, NotSet): + raise RuntimeError("m must be set for all kets to calculate the matrix element.") # noqa: TRY004 prefactor = self.kets[0]._calc_wigner_eckart_prefactor(other.kets[0], kappa, q) # noqa: SLF001 reduced_matrix_element = self.calc_reduced_matrix_element(other, operator, kappa) diff --git a/src/rydstate/angular/utils.py b/src/rydstate/angular/utils.py index fa41889..5a2721c 100644 --- a/src/rydstate/angular/utils.py +++ b/src/rydstate/angular/utils.py @@ -1,6 +1,7 @@ from __future__ import annotations import contextlib +import typing as t from typing import TYPE_CHECKING, Literal import numpy as np @@ -12,6 +13,19 @@ CouplingScheme = Literal["LS", "JJ", "FJ"] +@t.runtime_checkable +class NotSet(t.Protocol): + """Singleton for a not set value and type at the same time. + + See Also: + https://stackoverflow.com/questions/77571796/how-to-create-singleton-object-which-could-be-used-both-as-type-and-value-simi + + """ + + @staticmethod + def __not_set() -> None: ... + + class InvalidQuantumNumbersError(ValueError): def __init__(self, ket: AngularKetBase, msg: str = "") -> None: _msg = f"Invalid quantum numbers for {ket!r}" @@ -74,7 +88,7 @@ def quantum_numbers_to_angular_ket( l_tot: int | None = None, j_tot: float | None = None, f_tot: float | None = None, - m: float | None = None, + m: float | NotSet = NotSet, ) -> AngularKetBase: r"""Return an AngularKet object in the corresponding coupling scheme from the given quantum numbers. diff --git a/src/rydstate/rydberg/rydberg_sqdt.py b/src/rydstate/rydberg/rydberg_sqdt.py index 13ce75a..f99c1ee 100644 --- a/src/rydstate/rydberg/rydberg_sqdt.py +++ b/src/rydstate/rydberg/rydberg_sqdt.py @@ -7,6 +7,7 @@ import numpy as np +from rydstate.angular import NotSet from rydstate.angular.utils import quantum_numbers_to_angular_ket from rydstate.radial import RadialKet from rydstate.rydberg.rydberg_base import RydbergStateBase @@ -47,7 +48,7 @@ def __init__( l_tot: int | None = None, j_tot: float | None = None, f_tot: float | None = None, - m: float | None = None, + m: float | NotSet = NotSet, ) -> None: r"""Initialize the Rydberg state. @@ -331,7 +332,7 @@ def __init__( l: int, j: float | None = None, f: float | None = None, - m: float | None = None, + m: float | NotSet = NotSet, nu: float | None = None, ) -> None: r"""Initialize the Rydberg state. @@ -379,7 +380,7 @@ def __init__( s_tot: int, j_tot: int | None = None, f_tot: float | None = None, - m: float | None = None, + m: float | NotSet = NotSet, nu: float | None = None, ) -> None: r"""Initialize the Rydberg state. @@ -428,7 +429,7 @@ def __init__( j_r: float, j_tot: int | None = None, f_tot: float | None = None, - m: float | None = None, + m: float | NotSet = NotSet, nu: float | None = None, ) -> None: r"""Initialize the Rydberg state. @@ -477,7 +478,7 @@ def __init__( j_r: float, f_c: float | None = None, f_tot: float | None = None, - m: float | None = None, + m: float | NotSet = NotSet, nu: float | None = None, ) -> None: r"""Initialize the Rydberg state. From 0021bbf0b4dd6f54e460761904efdce6bd673e09 Mon Sep 17 00:00:00 2001 From: johannes-moegerle Date: Mon, 9 Mar 2026 18:49:05 +0100 Subject: [PATCH 03/10] update basis to accept n as tuple and also m --- src/rydstate/basis/basis_sqdt.py | 163 +++++++++++++++++-------------- tests/test_basis.py | 6 +- 2 files changed, 94 insertions(+), 75 deletions(-) diff --git a/src/rydstate/basis/basis_sqdt.py b/src/rydstate/basis/basis_sqdt.py index 8892db1..2248d40 100644 --- a/src/rydstate/basis/basis_sqdt.py +++ b/src/rydstate/basis/basis_sqdt.py @@ -1,11 +1,14 @@ from __future__ import annotations import logging +from typing import TypeVar import numpy as np +from rydstate.angular.utils import NotSet from rydstate.basis.basis_base import BasisBase from rydstate.rydberg import ( + RydbergStateSQDT, RydbergStateSQDTAlkali, RydbergStateSQDTAlkalineFJ, RydbergStateSQDTAlkalineJJ, @@ -15,114 +18,130 @@ logger = logging.getLogger(__name__) +_RydbergStateSQDT = TypeVar("_RydbergStateSQDT", bound=RydbergStateSQDT) -class BasisSQDTAlkali(BasisBase[RydbergStateSQDTAlkali]): - def __init__(self, species: str | SpeciesObjectSQDT, n_min: int = 1, n_max: int | None = None) -> None: + +class BasisSQDT(BasisBase[_RydbergStateSQDT]): + def _get_m_range(self, m: tuple[float, float] | None | NotSet, f_tot: float | np.floating) -> list[NotSet | float]: + if isinstance(m, NotSet): + return [NotSet] + if m is None: + m = (-np.inf, np.inf) + + m_min = max(-f_tot, m[0]) + m_max = min(f_tot, m[1]) + return [float(_m) for _m in np.arange(m_min, m_max + 1)] + + +class BasisSQDTAlkali(BasisSQDT[RydbergStateSQDTAlkali]): + def __init__( + self, species: str | SpeciesObjectSQDT, n: tuple[int, int], m: tuple[float, float] | None | NotSet = NotSet + ) -> None: if isinstance(species, str): species = SpeciesObjectSQDT.from_name(species) self.species = species - if n_max is None: - raise ValueError("n_max must be given") - - s = 1 / 2 i_c = self.species.i_c if self.species.i_c is not None else 0 - self.states = [] - for n in range(n_min, n_max + 1): - for l in range(n): - if not self.species.is_allowed_shell(n, l, s): + + _s = 1 / 2 + for _n in range(n[0], n[1] + 1): + for _l in range(_n): + if not self.species.is_allowed_shell(_n, _l, _s): continue - for j in np.arange(abs(l - s), l + s + 1): - for f in np.arange(abs(j - i_c), j + i_c + 1): - state = RydbergStateSQDTAlkali(species, n=n, l=l, j=float(j), f=float(f)) - self.states.append(state) + for _j in np.arange(abs(_l - _s), _l + _s + 1): + for _f in np.arange(abs(_j - i_c), _j + i_c + 1): + for _m in self._get_m_range(m, _f): + state = RydbergStateSQDTAlkali(species, n=_n, l=_l, j=float(_j), f=float(_f), m=_m) + self.states.append(state) -class BasisSQDTAlkalineLS(BasisBase[RydbergStateSQDTAlkalineLS]): - def __init__(self, species: str | SpeciesObjectSQDT, n_min: int = 1, n_max: int | None = None) -> None: +class BasisSQDTAlkalineLS(BasisSQDT[RydbergStateSQDTAlkalineLS]): + def __init__( + self, species: str | SpeciesObjectSQDT, n: tuple[int, int], m: tuple[float, float] | None | NotSet = NotSet + ) -> None: if isinstance(species, str): species = SpeciesObjectSQDT.from_name(species) self.species = species - if n_max is None: - raise ValueError("n_max must be given") - i_c = self.species.i_c if self.species.i_c is not None else 0 - self.states = [] - for n in range(n_min, n_max + 1): - for l in range(n): - for s_tot in [0, 1]: - if not self.species.is_allowed_shell(n, l, s_tot): - continue - for j_tot in range(abs(l - s_tot), l + s_tot + 1): - for f_tot in np.arange(abs(j_tot - i_c), j_tot + i_c + 1): - state = RydbergStateSQDTAlkalineLS( - species, n=n, l=l, s_tot=s_tot, j_tot=j_tot, f_tot=float(f_tot) - ) - self.states.append(state) - -class BasisSQDTAlkalineJJ(BasisBase[RydbergStateSQDTAlkalineJJ]): - def __init__(self, species: str | SpeciesObjectSQDT, n_min: int = 0, n_max: int | None = None) -> None: + for _n in range(n[0], n[1] + 1): + for _l in range(_n): + for _s_tot in [0, 1]: + if not self.species.is_allowed_shell(_n, _l, _s_tot): + continue + for _j_tot in range(abs(_l - _s_tot), _l + _s_tot + 1): + for _f_tot in np.arange(abs(_j_tot - i_c), _j_tot + i_c + 1): + for _m in self._get_m_range(m, _f_tot): + state = RydbergStateSQDTAlkalineLS( + species, n=_n, l=_l, s_tot=_s_tot, j_tot=_j_tot, f_tot=float(_f_tot), m=_m + ) + self.states.append(state) + + +class BasisSQDTAlkalineJJ(BasisSQDT[RydbergStateSQDTAlkalineJJ]): + def __init__( + self, species: str | SpeciesObjectSQDT, n: tuple[int, int], m: tuple[float, float] | None | NotSet = NotSet + ) -> None: if isinstance(species, str): species = SpeciesObjectSQDT.from_name(species) self.species = species - if n_max is None: - raise ValueError("n_max must be given") - i_c = self.species.i_c if self.species.i_c is not None else 0 - j_c = 0.5 - s_r = 0.5 self.states = [] - for n in range(n_min, n_max + 1): - for l_r in range(n): - if self.species.is_allowed_shell(n, l_r, 0) != self.species.is_allowed_shell(n, l_r, 1): + + _j_c = 0.5 + _s_r = 0.5 + for _n in range(n[0], n[1] + 1): + for _l_r in range(_n): + if self.species.is_allowed_shell(_n, _l_r, 0) != self.species.is_allowed_shell(_n, _l_r, 1): logger.warning( "For l=%d, n=%d one of the singlet/triplet states is not allowed. " "In JJ coupling the state does not exist, thus skipping this shell", - *(l_r, n), + *(_l_r, _n), ) - if not all(self.species.is_allowed_shell(n, l_r, s_tot) for s_tot in [0, 1]): + if not all(self.species.is_allowed_shell(_n, _l_r, s_tot) for s_tot in [0, 1]): continue - for j_r in np.arange(abs(l_r - s_r), l_r + s_r + 1): - for j_tot in range(int(abs(j_r - j_c)), int(j_r + j_c + 1)): - for f_tot in np.arange(abs(j_tot - i_c), j_tot + i_c + 1): - state = RydbergStateSQDTAlkalineJJ( - species, n=n, l=l_r, j_r=float(j_r), j_tot=j_tot, f_tot=float(f_tot) - ) - self.states.append(state) - - -class BasisSQDTAlkalineFJ(BasisBase[RydbergStateSQDTAlkalineFJ]): - def __init__(self, species: str | SpeciesObjectSQDT, n_min: int = 0, n_max: int | None = None) -> None: + for _j_r in np.arange(abs(_l_r - _s_r), _l_r + _s_r + 1): + for _j_tot in range(int(abs(_j_r - _j_c)), int(_j_r + _j_c + 1)): + for _f_tot in np.arange(abs(_j_tot - i_c), _j_tot + i_c + 1): + for _m in self._get_m_range(m, _f_tot): + state = RydbergStateSQDTAlkalineJJ( + species, n=_n, l=_l_r, j_r=float(_j_r), j_tot=_j_tot, f_tot=float(_f_tot), m=_m + ) + self.states.append(state) + + +class BasisSQDTAlkalineFJ(BasisSQDT[RydbergStateSQDTAlkalineFJ]): + def __init__( + self, species: str | SpeciesObjectSQDT, n: tuple[int, int], m: tuple[float, float] | None | NotSet = NotSet + ) -> None: if isinstance(species, str): species = SpeciesObjectSQDT.from_name(species) self.species = species - if n_max is None: - raise ValueError("n_max must be given") - i_c = self.species.i_c if self.species.i_c is not None else 0 - j_c = 0.5 - s_r = 0.5 self.states = [] - for n in range(n_min, n_max + 1): - for l_r in range(n): - if self.species.is_allowed_shell(n, l_r, 0) != self.species.is_allowed_shell(n, l_r, 1): + + _j_c = 0.5 + _s_r = 0.5 + for _n in range(n[0], n[1] + 1): + for _l_r in range(_n): + if self.species.is_allowed_shell(_n, _l_r, 0) != self.species.is_allowed_shell(_n, _l_r, 1): logger.warning( "For l=%d, n=%d one of the singlet/triplet states is not allowed. " "In FJ coupling the state does not exist, thus skipping this shell", - *(l_r, n), + *(_l_r, _n), ) - if not all(self.species.is_allowed_shell(n, l_r, s_tot) for s_tot in [0, 1]): + if not all(self.species.is_allowed_shell(_n, _l_r, s_tot) for s_tot in [0, 1]): continue - for j_r in np.arange(abs(l_r - s_r), l_r + s_r + 1): - for f_c in np.arange(abs(j_c - i_c), j_c + i_c + 1): - for f_tot in np.arange(abs(f_c - j_r), f_c + j_r + 1): - state = RydbergStateSQDTAlkalineFJ( - species, n=n, l=l_r, j_r=float(j_r), f_c=float(f_c), f_tot=float(f_tot) - ) - self.states.append(state) + for _j_r in np.arange(abs(_l_r - _s_r), _l_r + _s_r + 1): + for _f_c in np.arange(abs(_j_c - i_c), _j_c + i_c + 1): + for _f_tot in np.arange(abs(_f_c - _j_r), _f_c + _j_r + 1): + for _m in self._get_m_range(m, _f_tot): + state = RydbergStateSQDTAlkalineFJ( + species, n=_n, l=_l_r, j_r=float(_j_r), f_c=float(_f_c), f_tot=float(_f_tot), m=_m + ) + self.states.append(state) diff --git a/tests/test_basis.py b/tests/test_basis.py index dd49e48..5753666 100644 --- a/tests/test_basis.py +++ b/tests/test_basis.py @@ -6,7 +6,7 @@ @pytest.mark.parametrize("species_name", ["Rb", "Na", "H"]) def test_alkali_basis(species_name: str) -> None: """Test alkali basis creation.""" - basis = BasisSQDTAlkali(species_name, n_min=1, n_max=20) + basis = BasisSQDTAlkali(species_name, n=(1, 20)) basis.sort_states("n", "l_r") lowest_n_state = {"Rb": (4, 2), "Na": (3, 0), "H": (1, 0)}[species_name] assert (basis.states[0].n, basis.states[0].l) == lowest_n_state @@ -35,7 +35,7 @@ def test_alkali_basis(species_name: str) -> None: @pytest.mark.parametrize("species_name", ["Sr88", "Sr87", "Yb174", "Yb171"]) def test_alkaline_basis(species_name: str) -> None: """Test alkaline basis creation.""" - basis = BasisSQDTAlkalineLS(species_name, n_min=30, n_max=35) + basis = BasisSQDTAlkalineLS(species_name, n=(30, 35)) basis.sort_states("n", "l_r") assert (basis.states[0].n, basis.states[0].l) == (30, 0) assert (basis.states[-1].n, basis.states[-1].l) == (35, 34) @@ -62,5 +62,5 @@ def test_alkaline_basis(species_name: str) -> None: assert np.shape(me_matrix) == (len(basis.states), len(basis.states)) assert np.count_nonzero(me_matrix) > 0 - basis = BasisSQDTAlkalineLS(species_name, n_min=30, n_max=35) + basis = BasisSQDTAlkalineLS(species_name, n=(30, 35)) basis.filter_states("l_r", (6, 10)) From 42d1895d599606f20f92156bd5f5a07bb36be04c Mon Sep 17 00:00:00 2001 From: johannes-moegerle Date: Mon, 9 Mar 2026 12:30:47 +0100 Subject: [PATCH 04/10] implement lifetimes again --- docs/examples.rst | 1 + docs/examples/lifetimes.ipynb | 244 +++++++++++++++++++++++++++ src/rydstate/rydberg/rydberg_sqdt.py | 207 ++++++++++++++++++++++- 3 files changed, 451 insertions(+), 1 deletion(-) create mode 100644 docs/examples/lifetimes.ipynb diff --git a/docs/examples.rst b/docs/examples.rst index a6e3f4c..687821a 100644 --- a/docs/examples.rst +++ b/docs/examples.rst @@ -26,6 +26,7 @@ Some examples demonstrating the usage of the RydbergState class. .. nbgallery:: examples/dipole_matrix_elements + examples/lifetimes Comparisons diff --git a/docs/examples/lifetimes.ipynb b/docs/examples/lifetimes.ipynb new file mode 100644 index 0000000..3381acd --- /dev/null +++ b/docs/examples/lifetimes.ipynb @@ -0,0 +1,244 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Lifetimes" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We calculate the lifetime of Rydberg states and analyze which transitions contribute to it. The calculations can be performed with a few lines of code; most of the code within this notebook is only needed for plotting." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We import the libraries that we will use within the notebook." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "from typing import Any\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from scipy.optimize import curve_fit\n", + "\n", + "import rydstate" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As a quick example, we show how to calculate the lifetime of the Rubidium $|60S,m=1/2\\rangle$ state via `state.get_lifetime`. At zero temperature, the lifetime is determined by the spontaneous decay. If the temperature is non-zero, black body radiation can drive transitions to neighboring Rydberg states, reducing the lifetime." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Lifetime at T=0: 235.77 microsecond\n", + "Lifetime at T=300K: 100.50 microsecond\n" + ] + } + ], + "source": [ + "state = rydstate.RydbergStateSQDTAlkali(\"Rb\", n=60, l=0, j=0.5, m=0.5)\n", + "\n", + "temperature = 300 # Kelvin\n", + "lifetime_0 = state.get_lifetime()\n", + "lifetime = state.get_lifetime(temperature, temperature_unit=\"K\")\n", + "\n", + "print(f\"Lifetime at T=0: {lifetime_0.to('mus'):.2f}\")\n", + "print(f\"Lifetime at T={temperature}K: {lifetime.to('mus'):.2f}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Transition rates contributing to the Rydberg lifetime" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To analyze which transitions contribute to the lifetime, we can obtain the transition rates from spontaneous decay (`state.get_spontaneous_transition_rates`) and black body radiation (`state.get_black_body_transition_rates`)." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Number of possible spontaneous decay transitions: 275\n", + "Number of considered BBR transitions: 435\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj8AAAHcCAYAAAA5lMuGAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAeg9JREFUeJzt3X1czff/P/BHF06XKpROWepYzUUuIkouhmlysU2YYUZytQtMa2ZsijFrmLn+aLYRGxM2ZmwmYRckFwkpxhYZnUKqiYrO6/eHb++ft04X5zgpzuN+u50b5/V+vl/v53md865n7/N6v98mQggBIiIiIiNhWtMJEBERET1KLH6IiIjIqLD4ISIiIqPC4oeIiIiMCosfIiIiMiosfoiIiMiosPghIiIio8Lih4iIiIwKix8iIiIyKix+iIiIyKiw+CEiIiKjwuKHiGqNWbNmwcTEpKbTMFoxMTEwMTHBhQsXDBpLVNuw+KFqd+rUKbz88stwd3eHpaUlGjVqhOeffx7Lli2r6dQAAAcPHsSsWbOQm5tb06lQFZX+4r3/0bBhQ/To0QO//PKL3rHa4s3NzdGoUSOMGjUKly9frjS3B7dV3mP//v2GGo5qVRv3j4d9j8pTG18rVQ/zmk6AnmwHDx5Ejx490LhxY4wbNw5KpRKXLl3CoUOHsGTJEkyaNKmmU8TBgwfx0UcfYdSoUXBwcKjpdEgHs2fPhkqlghACWVlZiImJQd++ffHTTz/hhRde0Dv2/vjCwkIcOnQIMTEx+PPPP5GSkgJLS8tyc/rmm29kz9etW4e4uLgy7c2bN3+IV149RowYgaFDh8LCwkJqK2//0Bb7qOn7HpWHPwuMB4sfqlZz586Fvb09jhw5UuaHSXZ2ds0kRU+MPn36oH379tLzMWPGwNnZGd99912ZgkaX2Afjx44dC0dHR8ybNw/bt2/HK6+8Um5Or732muz5oUOHEBcXV6a9PAUFBbCxsalSrKGZmZnBzMzM4LHVRd/3iIhfe1G1+vvvv+Ht7a31r6iGDRtK/y+d63HmzBm88sorsLOzQ4MGDTB58mQUFhaWWff48ePo06cP7OzsYGtri549e+LQoUOymNI+z58/L/0lZ29vj9DQUNy6dUuKee+99wAAKpVKOoxeOo/h4sWLeOutt9C0aVNYWVmhQYMGGDx4cJl5DlXZ1v0uX76M0aNHw9nZGRYWFvD29sbq1av1ep2jRo2Ch4dHmXUfnD/z33//ISwsDB4eHrCwsEDDhg3x/PPPIykpqcy696vqGOg6Dn/++Sc6dOgAS0tLPP300/jiiy8qzKMqHBwcYGVlBXPzyv+u0yUWALp27Qrg3mfaUErHKzU1Fa+++irq1auHLl26ADD8Z68q7/+D83gq2j/Km/NjqH1TH9reI13GsaKfBVXZZ/XZx7Zs2QITExP89ttvZZZ98cUXMDExQUpKit79l+rRoweeffZZJCUloU+fPqhbty4aNWqEJUuWVLruk4hHfqhaubu7IyEhASkpKWjZsmWl8a+88go8PDwQFRWFQ4cOYenSpbhx4wbWrVsnxZw+fRpdu3aFnZ0dpk6dijp16uCLL75A9+7d8dtvv8Hf379MnyqVClFRUUhKSsJXX32Fhg0bYt68eRg4cCD++usvfPfdd1i0aBEcHR0BAE5OTgCAI0eO4ODBgxg6dCieeuopXLhwAStXrkT37t2RmpoKa2vrKm+rVFZWFjp27AgTExNMnDgRTk5O+OWXXzBmzBjk5+cjLCxMr9dZmTfeeANbtmzBxIkT0aJFC1y/fh1//vkn0tLS0K5du3LX03UMqjIOp06dQq9eveDk5IRZs2bh7t27mDlzJpydnXV6TXl5ebh27RqEEMjOzsayZctw8+ZNrUdZdInVpvSXYL169XTKsSoGDx4MLy8vfPLJJxBCADD8Z0+f97+y/eNBhtw39aHtParqOFb0Wqu6z+ozxv369YOtrS02bdqEbt26yZbFxsbC29tb+tmp7z4M3NvnXF1d8eKLLyI0NBTBwcH48ssv8c477+C5555Dq1atdB7vx5ogqka7d+8WZmZmwszMTAQEBIipU6eKX3/9VRQXF8viZs6cKQCIl156Sdb+1ltvCQDixIkTUltwcLBQKBTi77//ltquXLki6tatK5599tkyfY4ePVrW54ABA0SDBg2k5wsWLBAARHp6epn8b926VaYtISFBABDr1q3TeVtCCDFmzBjh4uIirl27JmsfOnSosLe3l7ZZ1dcZEhIi3N3dy+RZmlMpe3t7MWHChDJxlanqGNy/zcrGITg4WFhaWoqLFy9KbampqcLMzExU5cfSmjVrBIAyDwsLCxETE6N37P3xe/bsEVevXhWXLl0SW7ZsEU5OTsLCwkJcunSp0vzuN2HChHJfU+l4DRs2rMwyQ3/2qvL+l772+/eF8vYPbbGG3jcry7Mq75Eun9/yXmtV91l997Fhw4aJhg0birt370ptmZmZwtTUVMyePVtq07f/K1euCADCyclJNjapqakCgFi7dq3OfT7u+LUXVavnn38eCQkJeOmll3DixAnMnz8fQUFBaNSoEbZv314mfsKECbLnpROif/75ZwBASUkJdu/ejeDgYDRp0kSKc3Fxwauvvoo///wT+fn5sj7eeOMN2fOuXbvi+vXrZeK0sbKykv5/584dXL9+HZ6ennBwcNB6qLmybQkh8P333+PFF1+EEALXrl2THkFBQcjLy0NSUpJer7MyDg4OSExMxJUrV3RaT9cxACoeh5KSEvz6668IDg5G48aNpZjmzZsjKChIp9xWrFiBuLg4xMXF4dtvv0WPHj0wduxY/PDDDw8VCwCBgYFwcnKCm5sbXn75ZdjY2GD79u146qmndMqxKh4cL8Dwnz193/+qetT7JlC190ifz+/9qrrPAvqP8ZAhQ5CdnS07A3DLli3QaDQYMmSI1KZv/6dOnQIAzJw5UzY2derUAQAoFAqd+nsSsPihatehQwf88MMPuHHjBg4fPozp06fjv//+w8svv4zU1FRZrJeXl+z5008/DVNTU+lw9tWrV3Hr1i00bdq0zHaaN28OjUaDS5cuydrv/wUL/P9D4jdu3Kg099u3byMyMhJubm6wsLCAo6MjnJyckJubi7y8vDLxlW3r6tWryM3NxapVq+Dk5CR7hIaGArg3EVyf11mZ+fPnIyUlBW5ubvDz88OsWbPwzz//GHwMKhuHq1ev4vbt22XeawBaX29F/Pz8EBgYiMDAQAwfPhw7d+5EixYtMHHiRBQXF+sdC/z/YmnLli3o27cvrl27Vm1nNqlUqjJthv7s6fv+V9Wj3jeBqr1H+nx+H3xdVdlnAf3HuHfv3rC3t0dsbKzUFhsbCx8fHzzzzDNSm779lxY/wcHBsvYzZ84A0H2/exKw+KFHRqFQoEOHDvjkk0+wcuVK3LlzB5s3b65wHUNc8K68M1LE/82tqMikSZMwd+5cvPLKK9i0aRN2796NuLg4NGjQABqNRudtla7z2muvSUchHnx07ty5qi8NQPljVFJSInv+yiuv4J9//sGyZcvg6uqKBQsWwNvbW+u1bu6n6xgADzfmD8PU1BQ9evRAZmYmzp0791CxpcXSoEGDsH37drRs2RKvvvoqbt68afC87z86UcrQnz193//q9LCfk6q8R/p8fu+nyz6r7xhbWFggODgYW7duxd27d3H58mUcOHBAdtTnYfo/efIklEolGjVqJGs/ceIEzM3N0aJFi0rH4UnDCc9UI0pPT83MzJS1nzt3TvZX8Pnz56HRaKSzmZycnGBtbY2zZ8+W6fPMmTMwNTWFm5ubTrlUVGBt2bIFISEhWLhwodRWWFio90XQnJycULduXZSUlCAwMLDcuJKSkiq/znr16mnN5+LFi2XaXFxc8NZbb+Gtt95CdnY22rVrh7lz56JPnz7l5lIdY2BlZaW14ND2enV19+5dAKhSkVLVWDMzM0RFRaFHjx5Yvnw5pk2b9tB5VsbQ4w7o9/5X9Q+Q6tg3dVHee6TLOGp7rVXdZ0vpM8bAva++1q5di/j4eKSlpUEIUab40bf/U6dOoU2bNmXaT548iWeeeaZGr9VUU3jkh6rVvn37tP4VVzqH58HDrStWrJA9L70KdOmObWZmhl69euHHH3+UnaqalZWFDRs2oEuXLrCzs9Mpx9Jrqmj7YWhmZlYm/2XLlpU5qlJVZmZmGDRoEL7//nvp9NX7Xb16VYqr6ut8+umnkZeXh5MnT0pxmZmZ2Lp1q/S8pKSkzCH+hg0bwtXVFUVFRZXmbOgxCAoKwrZt25CRkSG1p6Wl4ddff9Wrz1J37tzB7t27oVAoKr2IoC6xANC9e3f4+flh8eLFWi+/YGiGHPeHef8r2j8ezNfQ+6autL1Huoyjttda1X32YcYYuDd/qX79+oiNjUVsbCz8/Pxkfwjq239JSQnS0tK0Fj8nTpxA69atK83tScQjP1StJk2ahFu3bmHAgAFo1qwZiouLcfDgQcTGxsLDw0P6zrxUeno6XnrpJfTu3RsJCQn49ttv8eqrr8p23I8//hhxcXHo0qUL3nrrLZibm+OLL75AUVER5s+fr3OOvr6+AIAPP/wQQ4cORZ06dfDiiy/CxsYGL7zwAr755hvY29ujRYsWSEhIwJ49e9CgQQO9x+TTTz/Fvn374O/vj3HjxqFFixbIyclBUlIS9uzZg5ycHJ1e59ChQ/H+++9jwIABePvtt3Hr1i2sXLkSzzzzjDQR87///sNTTz2Fl19+GW3atIGtrS327NmDI0eOyP4i1qY6xuCjjz7Crl270LVrV7z11lu4e/culi1bBm9vb1kRV5lffvlFmreQnZ2NDRs24Ny5c5g2bVqZX7S6xJbnvffew+DBgxETE6N1krIhGXLcH+b9L2//0MbQ+6Y+HnyPdBnH8l5rVfbZhxlj4N7k44EDB2Ljxo0oKCjAZ599Jluub//nzp1DYWFhmeLn9u3bOH/+PEJCQqoyrE+eGjjDjIzIL7/8IkaPHi2aNWsmbG1thUKhEJ6enmLSpEkiKytLiis99TU1NVW8/PLLom7duqJevXpi4sSJ4vbt22X6TUpKEkFBQcLW1lZYW1uLHj16iIMHD8piSvu8evWqrF3bKbpz5swRjRo1EqamprJlN27cEKGhocLR0VHY2tqKoKAgcebMGeHu7i5CQkL02pYQQmRlZYkJEyYINzc3UadOHaFUKkXPnj3FqlWrdH6dQty7pEDLli2FQqEQTZs2Fd9++63sVPeioiLx3nvviTZt2oi6desKGxsb0aZNG/G///2vTF8PquoY6DoOv/32m/D19RUKhUI0adJEREdHlzk9vzzaTl+3tLQUPj4+YuXKlUKj0egVe3/8kSNHymy3pKREPP300+Lpp5+WnZZckaqc6v7geAlh2M9eVd//8j6v2vaP8mINvW9qo8t7pMvnt7zXKkTl++zD7GOl4uLiBABhYmJS5pIK+va/adMmAUCkpKTI2g8fPiwAiB07dlQ5vyeJiRDVPAORqApmzZqFjz76CFevXpUuLkZERFQdOOeHiIiIjAqLHyIiIjIqLH6IiIjIqHDODxERERkVHvkhIiIio8Lih4iIiIwKL3JYi2g0Gly5cgV169Y1yD2tiIiIjIUQAv/99x9cXV1halrxsR0WP7XIlStXqvXeN0RERE+6S5cu4amnnqowhsVPLVK3bl0A99646r4HDhER0ZMkPz8fbm5u0u/SirD4qUVKv+qys7Nj8UNERKSHqkwb4YRnIiIiMiosfoiIiMio1MriZ8WKFfDw8IClpSX8/f1x+PDhCuM3b96MZs2awdLSEq1atcLPP/8sWy6EQGRkJFxcXGBlZYXAwECcO3dOWn7hwgWMGTMGKpUKVlZWePrppzFz5kwUFxfL+jl58iS6du0KS0tLuLm5Yf78+TrnQkRERDWr1s35iY2NRXh4OKKjo+Hv74/FixcjKCgIZ8+eRcOGDcvEHzx4EMOGDUNUVBReeOEFbNiwAcHBwUhKSkLLli0BAPPnz8fSpUuxdu1aqFQqREREICgoCKmpqbC0tMSZM2eg0WjwxRdfwNPTEykpKRg3bhwKCgrw2WefAbg3kapXr14IDAxEdHQ0Tp06hdGjR8PBwQHjx4+vci4PSwiBu3fvoqSkxCD9EdGjYWZmBnNzc17GgqgWqHW3t/D390eHDh2wfPlyAPeufePm5oZJkyZh2rRpZeKHDBmCgoIC7NixQ2rr2LEjfHx8EB0dDSEEXF1d8e6772LKlCkAgLy8PDg7OyMmJgZDhw7VmseCBQuwcuVK/PPPPwCAlStX4sMPP4RarYZCoQAATJs2Ddu2bcOZM2eqlEtl8vPzYW9vj7y8PK0TnouLi5GZmYlbt25V2hcR1T7W1tZwcXGRfoYQkeFU9jv0frXqyE9xcTGOHTuG6dOnS22mpqYIDAxEQkKC1nUSEhIQHh4uawsKCsK2bdsAAOnp6VCr1QgMDJSW29vbw9/fHwkJCeUWP3l5eahfv75sO88++6zsh1ZQUBDmzZuHGzduoF69epXm8qCioiIUFRVJz/Pz87XGAfeKwPT0dJiZmcHV1RUKhYJ/QRI9JoQQKC4uxtWrV5Geng4vL69KL8JGRNWnVhU/165dQ0lJCZydnWXtzs7O0tGVB6nVaq3xarVaWl7aVl7Mg86fP49ly5ZJX3mV9qNSqcr0UbqsXr16lebyoKioKHz00Udalz2ouLhYOgpmbW1dpXWIqPawsrJCnTp1cPHiRRQXF8PS0rKmUyIyWvzT4wGXL19G7969MXjwYIwbN65atzV9+nTk5eVJj0uXLlW6Dv9aJHp8cf8lqh1q1Z7o6OgIMzMzZGVlydqzsrKgVCq1rqNUKiuML/23Kn1euXIFPXr0QKdOnbBq1aoqbef+bVSWy4MsLCykCxrywoZERESPRq0qfhQKBXx9fREfHy+1aTQaxMfHIyAgQOs6AQEBsngAiIuLk+JVKhWUSqUsJj8/H4mJibI+L1++jO7du8PX1xdr1qwp8xdaQEAAfv/9d9y5c0e2naZNm6JevXpVyoWIiIhqAVHLbNy4UVhYWIiYmBiRmpoqxo8fLxwcHIRarRZCCDFixAgxbdo0Kf7AgQPC3NxcfPbZZyItLU3MnDlT1KlTR5w6dUqK+fTTT4WDg4P48ccfxcmTJ0X//v2FSqUSt2/fFkII8e+//wpPT0/Rs2dP8e+//4rMzEzpUSo3N1c4OzuLESNGiJSUFLFx40ZhbW0tvvjiC51yqUheXp4AIPLy8sosu337tkhNTZVyJqLHD/djoupT0e/QB9WqCc/AvdPFr169isjISKjVavj4+GDXrl3SROKMjAzZUZlOnTphw4YNmDFjBj744AN4eXlh27ZtsuvqTJ06FQUFBRg/fjxyc3PRpUsX7Nq1S5pwGBcXh/Pnz+P8+fNl7gQr/u9KAPb29ti9ezcmTJgAX19fODo6IjIyUrrGT1VzqQ7t21dr92UcPapbfOn7uXPnTmRlZaFevXpo06YNIiMj0blz5+pJshweHh4ICwtDWFjYI91uTRk1ahTWrl0LADA3N0f9+vXRunVrDBs2DKNGjeIcFCIySrXuOj/GrKJrFBQWFiI9PR0qlarMWSK1vfh59tlnUVxcjKioKDRp0gRZWVmIj4+Ht7c3XnrppepJshzGWPxkZWVhzZo1KCkpQVZWFnbt2oWoqCh07doV27dvh7l5rfsb6IlV0X5MRA9Hl+v88M8+qla5ubn4448/MG/ePPTo0QPu7u7w8/PD9OnTZYWPiYkJVq5ciT59+sDKygpNmjTBli1bZH2dOnUKzz33HKysrNCgQQOMHz8eN2/elJaPGjUKwcHB+Oyzz+Di4oIGDRpgwoQJ0jyt7t274+LFi3jnnXdgYmIiXSfp+vXrGDZsGBo1agRra2u0atUK3333nWzb3bt3x9tvv42pU6eifv36UCqVmDVrVpnXOnbsWDg5OcHOzg7PPfccTpw4IYtZuXIlnn76aSgUCjRt2hTffPONtOzChQswMTFBcnKyrE8TExPs378fAHDjxg0MHz4cTk5OsLKygpeXF9asWVPhe2BhYQGlUolGjRqhXbt2+OCDD/Djjz/il19+QUxMjE75//TTT+jQoQMsLS3h6OiIAQMGSMu++eYbtG/fHnXr1oVSqcSrr76K7OxsAPeOoHp6esouHwEAycnJMDExwfnz5yt8DUREhsTih6qVra0tbG1tsW3bNtkFHbWJiIjAoEGDcOLECQwfPhxDhw5FWloaAKCgoABBQUGoV68ejhw5gs2bN2PPnj2YOHGirI99+/bh77//xr59+7B27VrExMRIv+B/+OEHPPXUU5g9ezYyMzORmZkJ4N5f476+vti5cydSUlIwfvx4jBgxosw95dauXQsbGxskJiZi/vz5mD17NuLi4qTlgwcPRnZ2Nn755RccO3YM7dq1Q8+ePZGTkwMA2Lp1KyZPnox3330XKSkpeP311xEaGop9+/ZVeTwjIiKQmpqKX375BWlpaVi5ciUcHR2rvH6p5557Dm3atMEPP/xQ5fx37tyJAQMGoG/fvjh+/Dji4+Ph5+cnrX/nzh3MmTMHJ06cwLZt23DhwgWMGjUKwL3idvTo0WUKtTVr1uDZZ5+Fp6enzq+BiEhv1Tz/iHSg74RnX99H+7jf6dPyhzZbtmwR9erVE5aWlqJTp05i+vTp4sSJE7IYAOKNN96Qtfn7+4s333xTCCHEqlWrRL169cTNmzel5Tt37hSmpqbSZPiQkBDh7u4u7t69K8UMHjxYDBkyRHru7u4uFi1apD3R+/Tr10+8++670vNu3bqJLl26yGI6dOgg3n//fSGEEH/88Yews7MThYWFspinn35amhTfqVMnMW7cONnywYMHi759+wohhEhPTxcAxPHjx6XlN27cEADEvn37hBBCvPjiiyI0NLTS/EuFhISI/v37a102ZMgQ0bx58yrnHxAQIIYPH17lbR85ckQAEP/9958QQojLly8LMzMzkZiYKIQQori4WDg6OoqYmJgq9/m444Rnouqjy4RnHvmhajdo0CBcuXIF27dvR+/evbF//360a9dO9pULgDKXBAgICJCO/KSlpaFNmzawsbGRlnfu3BkajQZnz56V2ry9vWFmZiY9d3Fxkb56KU9JSQnmzJmDVq1aoX79+rC1tcWvv/6KjIwMWVzr1q1lz+/v+8SJE7h58yYaNGggHe2ytbVFeno6/v77b+k1PDjBu3PnztJrrIo333wTGzduhI+PD6ZOnYqDBw9Wed0HCSGkr/6qkn9ycjJ69uxZbn/Hjh3Diy++iMaNG6Nu3bro1q0bAEjj6Orqin79+mH16tUA7n2FVlRUhMGDB+v9GugRaN9e/iB6AnCmIz0SlpaWeP755/H8888jIiICY8eOxcyZM6WvRQylTp06sucmJibQaDQVrrNgwQIsWbIEixcvRqtWrWBjY4OwsDAUFxdXue+bN2/CxcVFmptzPwcHhyrlXnrmlbjvHIT7rysFAH369MHFixfx888/Iy4uDj179sSECRPKzKWpirS0NOmWLVXJ38rKqty+Sr+WDAoKwvr16+Hk5ISMjAwEBQXJxnHs2LEYMWIEFi1ahDVr1mDIkCG8XQsRPXI88kM1okWLFigoKJC1HTp0qMzz5s2bAwCaN2+OEydOyNY5cOAATE1N0bRp0ypvV6FQoKSkRNZ24MAB9O/fH6+99hratGmDJk2a4K+//tLp9bRr1w5qtRrm5ubw9PSUPUrn5DRv3hwHDhwos+0WLVoAAJycnABAmosEQDb5uZSTkxNCQkLw7bffYvHixWWuRl4Ve/fuxalTpzBo0KAq59+6desyF/EsdebMGVy/fh2ffvopunbtimbNmmk94ta3b1/Y2Nhg5cqV2LVrF0aPHq1z7kRED4vFD1Wr69ev47nnnsO3336LkydPIj09HZs3b8b8+fPRv39/WezmzZuxevVq/PXXX5g5cyYOHz4sTWgePnw4LC0tERISgpSUFOzbtw+TJk3CiBEjytxMtiIeHh74/fffcfnyZVy7dg0A4OXlhbi4OBw8eBBpaWl4/fXXy9ympDKBgYEICAhAcHAwdu/ejQsXLuDgwYP48MMPcfT/rg3w3nvvISYmBitXrsS5c+fw+eef44cffsCUKVMA3Duy0rFjR3z66adIS0vDb7/9hhkzZsi2ExkZiR9//BHnz5/H6dOnsWPHDqlALE9RURHUajUuX76MpKQkfPLJJ+jfvz9eeOEFjBw5ssr5z5w5E9999x1mzpyJtLQ0nDp1CvPmzQMANG7cGAqFAsuWLcM///yD7du3Y86cOWVyMTMzw6hRozB9+nR4eXnx6udEVDOqfQYSVdnjeIXnyiY8FxYWimnTpol27doJe3t7YW1tLZo2bSpmzJghbt26JcUBECtWrBDPP/+8sLCwEB4eHiI2NlbW18mTJ0WPHj2EpaWlqF+/vhg3bpw0mVYI7ZN7J0+eLLp16yY9T0hIEK1btxYWFhai9ON//fp10b9/f2FraysaNmwoZsyYIUaOHCnrq1u3bmLy5Mmyvvv37y9CQkKk5/n5+WLSpEnC1dVV1KlTR7i5uYnhw4eLjIwMKeZ///ufaNKkiahTp4545plnxLp162R9pqamioCAAGFlZSV8fHzE7t27ZROe58yZI5o3by6srKxE/fr1Rf/+/cU///xTduDvGxMAAoAwNzcXTk5OIjAwUKxevVqUlJTIYquS//fffy98fHyEQqEQjo6OYuDAgdKyDRs2CA8PD2FhYSECAgLE9u3by0zgFkKIv//+WwAQ8+fPLzfvJ1Vt3Y8rVNEZD0S1iC4TnnmRw1pE34sc1qTUVPnz//sGR2cmJibYunUrgoODHzonqt3++OMP9OzZE5cuXdLpqN2ToLbuxxV6cJKzrlc5JXpEdLnIISc8E9EjUVRUhKtXr2LWrFkYPHiw0RU+RFR7cM4PET0S3333Hdzd3ZGbm4v58+fXdDpEZMR45IdqBX77+uQbNWqUwS9tQESkDx75ISIiIqPCIz9ERKSb+ydBcwI0PYZ45IeIiIiMCosfIiIiMiosfoiIiMiosPghIiIio8IJz0+CB6/AWt2qYYLjhQsXoFKpcPz4cfj4+Bikz9p01eiYmBiEhYUhNze3plPRysPDA2FhYQgLCwNguLF7lO/Bs88+izfeeAOvvvpqtW9Lm2vXrqFFixZISkrCU089VSM5EFHV8MgPVbtRo0bBxMREejRo0AC9e/fGyZMnazo1mZiYGFmetra28PX1xQ8//FDTqQEAunfvLuVmaWmJZ555BlFRUdVyjaTMzEz06dOnyvGzZs3SWrTq2o++tm/fjqysLAwdOlRq8/DwkMbLzMwMrq6uGDNmDG7cuCHF7N+/X/aeW1lZwdvbG6tWrZL1f/9nuE6dOlCpVJg6dSoKCwulGEdHR4wcORIzZ86s9tdLRA+HxQ89Er1790ZmZiYyMzMRHx8Pc3NzvPDCCzWdVhl2dnZSnsePH0dQUBBeeeUVnD17tqZTAwCMGzcOmZmZOHv2LKZPn47IyEhER0cbfDtKpRIWFha1pp/KLF26FKGhoTA1lf9Imz17NjIzM5GRkYH169fj999/x9tvv11m/bNnzyIzMxOpqal4/fXX8eabbyI+Pl4WU/oZ/ueff7Bo0SJ88cUXZQqd0NBQrF+/Hjk5OYZ/kURkMCx+6JGwsLCAUqmEUqmEj48Ppk2bhkuXLuHq1ata40tKSjBmzBioVCpYWVmhadOmWLJkSZm41atXw9vbGxYWFnBxccHEiRPLzWHmzJlwcXGp8IiTiYmJlKeXlxc+/vhjmJqayta5ceMGRo4ciXr16sHa2hp9+vTBuXPnZP3ExMSgcePGsLa2xoABA3D9+nVp2YULF2BqaoqjD3x9uHjxYri7u0Oj0ZSbn7W1NZRKJdzd3REaGorWrVsjLi5OWv7333+jf//+cHZ2hq2tLTp06IA9e/bI+sjOzsaLL74IKysrqFQqrF+/Xus4bNu2TXr+/vvv45lnnoG1tTWaNGmCiIgI3LlzR3qtH330EU6cOCEdHYmJidHaz6lTp/Dcc8/BysoKDRo0wPjx43Hz5k1p+ahRoxAcHIzPPvsMLi4uaNCgASZMmCBtS5urV69i7969ePHFF8ssq1u3LpRKJRo1aoQePXogJCQESUlJZeIaNmwIpVIJlUqFt99+GyqVqkxc6WfYzc0NwcHBCAwMlI09AHh7e8PV1RVbt24tN18iqnksfuiRu3nzJr799lt4enqiQYMGWmM0Gg2eeuopbN68GampqYiMjMQHH3yATZs2STErV67EhAkTMH78eJw6dQrbt2+Hp6dnmb6EEJg0aRLWrVuHP/74A61bt65SniUlJVi7di0AoF27dlL7qFGjcPToUWzfvh0JCQkQQqBv377SL+jExESMGTMGEydORHJyMnr06IGPP/5YWt/DwwOBgYFYs2aNbHtr1qzBqFGjyhy90EYIgT/++ANnzpyBQqGQ2m/evIm+ffsiPj4ex48fR+/evfHiiy8iIyNDlv+lS5ewb98+bNmyBf/73/+QnZ1d4fbq1q2LmJgYpKamYsmSJfjyyy+xaNEiAMCQIUPw7rvvwtvbWzpqNmTIkDJ9FBQUICgoCPXq1cORI0ewefNm7Nmzp0zBum/fPvz999/Yt28f1q5di5iYGKmY0ubPP/+EtbU1mjdvXuFruHz5Mn766Sf4+/uXGyOEwK5du5CRkVFhXEpKCg4ePCgb+1J+fn74448/KsyFiGqYoFojLy9PABB5eXlllt2+fVukpqaK27dvl13R1/fRPu5z+rT8oU1ISIgwMzMTNjY2wsbGRgAQLi4u4tixY1JMenq6ACCOHz9e7vhMmDBBDBo0SHru6uoqPvzww3LjAYjNmzeLV199VTRv3lz8+++/5cYKIcSaNWsEAClPU1NTYWFhIdasWSPF/PXXXwKAOHDggNR27do1YWVlJTZt2iSEEGLYsGGib9++sr6HDBki7O3tpeexsbGiXr16orCwUAghxLFjx4SJiYlIT08vN79u3bqJOnXqCBsbG1GnTh0BQFhaWspy0cbb21ssW7ZMCCHE2bNnBQBx+PBhaXlaWpoAIBYtWiS1ARBbt24tt88FCxYI3/s+CzNnzhRt2rQpE3d/P6tWrRL16tUTN2/elJbv3LlTmJqaCrVaLYS491lxd3cXd+/elWIGDx4shgwZUm4uixYtEk2aNCnT7u7uLhQKhbCxsRGWlpYCgPD39xc3btyQYvbt2yd7z83NzYWpqan4+OOPZX3d/xm2sLAQAISpqanYsmVLme2+8847onv37lpzrXA/rq207f/l/DwgqkkV/Q59EI/80CPRo0cPJCcnIzk5GYcPH0ZQUBD69OmDixcvlrvOihUr4OvrCycnJ9ja2mLVqlXSEYzs7GxcuXIFPXv2rHC777zzDhITE/H777+jUaNGleZZt25dKc/jx4/jk08+wRtvvIGffvoJAJCWlgZzc3PZUYEGDRqgadOmSEtLk2IePGoQEBAgex4cHAwzMzPp65GYmBj06NEDHh4eFeY3fPhwJCcn48CBA+jTpw8+/PBDdOrUSVp+8+ZNTJkyBc2bN4eDgwNsbW2RlpYmjVtp/r6+vtI6zZo1g4ODQ4XbjY2NRefOnaFUKmFra4sZM2bIjiZVRVpaGtq0aQMbGxuprXPnztBoNLI5Vd7e3jAzM5Oeu7i4VHhk6vbt27C0tNS67L333kNycjJOnjwpzeHp168fSkpKZHF//PGH9L5/9dVX+OSTT7By5UpZTOlnODExESEhIQgNDcWgQYPKbNPKygq3bt2qYCSIqKax+KFHwsbGBp6envD09ESHDh3w1VdfoaCgAF9++aXW+I0bN2LKlCkYM2YMdu/ejeTkZISGhqK4uBjAvV8wVfH888/j8uXL+PXXX6sUb2pqKuXZunVrhIeHo3v37pg3b17VXmgVKRQKjBw5EmvWrEFxcTE2bNiA0aNHV7qevb29NIabNm3C8uXLZXN6pkyZgq1bt+KTTz6RfqG3atVKGjd9JCQkYPjw4ejbty927NiB48eP48MPP3yoPitSp04d2XMTE5MK50E5OjrKzuB6cJmnpye8vLzw3HPPYfHixTh48CD27dsni1OpVPD09IS3tzdCQ0MxYsQIzJ07VxZT+hlu06YNVq9ejcTERHz99ddltpmTkwMnJ6eqvlwiqgEsfqhGmJiYwNTUFLdv39a6/MCBA+jUqRPeeusttG3bFp6envj777+l5XXr1oWHh0eZM3Ie9NJLL2HDhg0YO3YsNm7cqFeuZmZmUp7NmzfH3bt3kZiYKC2/fv06zp49ixYtWkgx9y8HgEOHDpXpd+zYsdizZw/+97//4e7duxg4cKBOedna2mLy5MmYMmWKdLr7gQMHMGrUKAwYMACtWrWCUqnEhQsXpHWaNWuGu3fv4tixY1Lb2bNnK7z+0MGDB+Hu7o4PP/wQ7du3h5eXV5kjdgqFoszRlAc1b94cJ06cQEFBgdR24MABmJqaomnTpjq8crm2bdtCrVaXWwDdr/SIUnmfu/vjKooxNTXFBx98gBkzZpSJS0lJQdu2bauQORHVFBY/9EgUFRVBrVZDrVYjLS0NkyZNws2bN7WeoQMAXl5eOHr0KH799Vf89ddfiIiIwJEjR2Qxs2bNwsKFC7F06VKcO3cOSUlJWLZsWZm+BgwYgG+++QahoaHYsmVLhXkKIaQ809PTsWrVKvz666/o37+/lFf//v0xbtw4/Pnnnzhx4gRee+01NGrUSIp5++23sWvXLnz22Wc4d+4cli9fjl27dpXZVvPmzdGxY0e8//77GDZsWJWPZt3v9ddfx19//YXvv/9eyu+HH35AcnIyTpw4gVdffVV21KRp06bo3bs3Xn/9dSQmJuLYsWMYO3Zshdv28vJCRkYGNm7ciL///htLly4tczaTh4cH0tPTkZycjGvXrqGoqKhMP8OHD4elpSVCQkKQkpKCffv2YdKkSRgxYgScnZ11fu2l2rZtC0dHRxw4cKDMsv/++w9qtRqZmZk4fPgw3nvvPTg5Ocm+KgTufY2qVqtx8eJFbN68Gd988430fpZn8ODBMDMzw4oVK6S2W7du4dixY+jVq5fer4eIHoFqn4FEVab3hOcaVNUJzwCkR926dUWHDh1kk0UfnPBcWFgoRo0aJezt7YWDg4N48803xbRp08pMqo2OjhZNmzYVderUES4uLmLSpEnSMjwwaTc2NlZYWlqK77//XmuepROeSx8WFhbimWeeEXPnzpVNwM3JyREjRowQ9vb2wsrKSgQFBYm//vpL1tfXX38tnnrqKWFlZSVefPFF8dlnn8kmPN8fhwcmIJenW7duYvLkyWXaX3/9deHt7S1KSkpEenq66NGjh7CyshJubm5i+fLlZdbLzMwU/fr1ExYWFqJx48Zi3bp1wt3dvcIJz++9955o0KCBsLW1FUOGDBGLFi2SvZ7CwkIxaNAg4eDgIABIk8Qf7OfkyZOiR48ewtLSUtSvX1+MGzdO/Pfff9LykJAQ0b9/f9nrmzx5sujWrVuFYzN16lQxdOhQWZu7u7vs/XRychJ9+/aVTaovnfBc+jA3NxcqlUpMmTJFNjFbW15CCBEVFSWcnJyk2A0bNoimTZuWm2dt3Y8rxAnP9JjQZcKziRDVcHlY0kt+fj7s7e2Rl5cHOzs72bLCwkKkp6dDpVKVO7mzJqSmyp//3zc/VEVz5szB5s2ba93Vrh83arUa3t7eSEpKgru7e43l0bFjR7z99tvl3mKjtu7HFXrw9jlHj8rbquF2N0T6qOh36IP4tRdRDbh58yZSUlKwfPlyTJo0qabTeewplUp8/fXXOp+BZkjXrl3DwIEDMWzYsBrLgYiqhjc2JaoBEydOxHfffYfg4OAqneVFlavpG9g6Ojpi6tSpNZoDEVUNix+iGlDZVYuJiKj68GsvIiIiMiosfh4znJ9O9Pji/ktUO7D4eUyUXvWWl80nenyV7r8PXsWaiB6tWjfnZ8WKFViwYAHUajXatGmDZcuWwc/Pr9z4zZs3IyIiAhcuXICXlxfmzZuHvn37SsuFEJg5cya+/PJL5ObmonPnzli5ciW8vLykmLlz52Lnzp1ITk6GQqEoc7XbmJgYhIaGat1+VlYWGjZsiP3796NHjx5llmdmZkKpVOo4CmWZmZnBwcFBuseRtbU1TExMHrrfh/XgXQcKC2smD6LaTAiBW7duITs7Gw4ODrJ7lxHRo1erip/Y2FiEh4cjOjoa/v7+WLx4MYKCgnD27Fk0bNiwTPzBgwcxbNgwREVF4YUXXsCGDRsQHByMpKQktGzZEgAwf/58LF26FGvXroVKpUJERASCgoKQmpoqXWejuLgYgwcPRkBAgNZ79QwZMgS9e/eWtY0aNQqFhYVl8jp79qzs+gLa8tZXaRFV0U0eH7UHU+HPdKLyOTg4GOSPISJ6OLXqIof+/v7o0KEDli9fDgDQaDRwc3PDpEmTMG3atDLxQ4YMQUFBAXbs2CG1dezYET4+PoiOjoYQAq6urnj33XcxZcoUAEBeXh6cnZ0RExODoUOHyvqLiYlBWFhYhfc5AoCrV6+iUaNG+PrrrzFixAgAkI783Lhxo9I7ZJenqhdoKikpwZ07d/TahqE9eFPr/7vLAhE9oE6dOo/nER9e5JAeE7pc5LDWHPkpLi7GsWPHMH36dKnN1NQUgYGBSEhI0LpOQkICwsPDZW1BQUHYtm0bACA9PR1qtRqBgYHScnt7e/j7+yMhIaFM8VNV69atg7W1NV5++eUyy3x8fFBUVISWLVti1qxZ6Ny5c7n9FBUVye6BlJ+fX6Xtm5mZ1ZofollZ8uePy0VriYjIeNWaCc/Xrl1DSUlJmRscOjs7Q61Wa11HrVZXGF/6ry59VsXXX3+NV199VXYzSBcXF0RHR+P777/H999/Dzc3N3Tv3h1JSUnl9hMVFQV7e3vp4ebmpndOREREVDW15sjP4yIhIQFpaWn45ptvZO1NmzZF06ZNpeedOnXC33//jUWLFpWJLTV9+nTZkav8/HwWQERERNWs1hz5cXR0hJmZGbIe+B4lKyur3AmCSqWywvjSf3XpszJfffUVfHx84OvrW2msn58fzp8/X+5yCwsL2NnZyR5ERERUvWpN8aNQKODr64v4+HipTaPRID4+HgEBAVrXCQgIkMUDQFxcnBSvUqmgVCplMfn5+UhMTCy3z4rcvHkTmzZtwpgxY6oUn5ycDBcXF523Q0RERNWnVn3tFR4ejpCQELRv3x5+fn5YvHgxCgoKpGvsjBw5Eo0aNUJUVBQAYPLkyejWrRsWLlyIfv36YePGjTh69ChWrVoFADAxMUFYWBg+/vhjeHl5Sae6u7q6ym6CmJGRgZycHGRkZKCkpATJyckAAE9PT9ja2kpxsbGxuHv3Ll577bUyuS9evBgqlQre3t4oLCzEV199hb1792L37t3VNFpERESkj1pV/AwZMgRXr15FZGQk1Go1fHx8sGvXLmnCckZGBkxN///Bqk6dOmHDhg2YMWMGPvjgA3h5eWHbtm3SNX4AYOrUqSgoKMD48eORm5uLLl26YNeuXdI1fgAgMjISa9eulZ63bdsWALBv3z50795dav/6668xcOBAraeyFxcX491338Xly5dhbW2N1q1bY8+ePVovfEhEREQ1p1Zd58fY6XKNgtpC2yVAiOgJwuv80GNCl9+htWbODxEREdGjwOKHiIiIjAqLHyIiIjIqLH6IiIjIqLD4ISIiIqPC4oeIiIiMCosfIiIiMiosfoiIiMiosPghIiIio8Lih4iIiIwKix8iIiIyKix+iIiIyKiw+CEiIiKjwuKHiIiIjAqLHyIiIjIqLH6IiIjIqLD4ISIiIqPC4oeIiIiMCosfIiIiMiosfoiIiMiosPghIiIio8Lih4iIiIwKix8iIiIyKix+iIiIyKiw+CEiIiKjwuKHiIiIjAqLHyIiIjIqLH6IiIjIqLD4ISIiIqPC4oeIiIiMCosfIiIiMiosfoiIiMiosPghIiIio8Lih4iIiIwKix8iIiIyKix+iIiIyKjUuuJnxYoV8PDwgKWlJfz9/XH48OEK4zdv3oxmzZrB0tISrVq1ws8//yxbLoRAZGQkXFxcYGVlhcDAQJw7d04WM3fuXHTq1AnW1tZwcHDQuh0TE5Myj40bN8pi9u/fj3bt2sHCwgKenp6IiYnR+fUTERFR9apVxU9sbCzCw8Mxc+ZMJCUloU2bNggKCkJ2drbW+IMHD2LYsGEYM2YMjh8/juDgYAQHByMlJUWKmT9/PpYuXYro6GgkJibCxsYGQUFBKCwslGKKi4sxePBgvPnmmxXmt2bNGmRmZkqP4OBgaVl6ejr69euHHj16IDk5GWFhYRg7dix+/fXXhxsUIiIiMigTIYSo6SRK+fv7o0OHDli+fDkAQKPRwM3NDZMmTcK0adPKxA8ZMgQFBQXYsWOH1NaxY0f4+PggOjoaQgi4urri3XffxZQpUwAAeXl5cHZ2RkxMDIYOHSrrLyYmBmFhYcjNzS2zLRMTE2zdulVW8Nzv/fffx86dO2WF19ChQ5Gbm4tdu3ZV6fXn5+fD3t4eeXl5sLOzq9I6Na19e/nzo0drJg8iqibadvL727jTUy2hy+/QWnPkp7i4GMeOHUNgYKDUZmpqisDAQCQkJGhdJyEhQRYPAEFBQVJ8eno61Gq1LMbe3h7+/v7l9lmRCRMmwNHREX5+fli9ejXurxsry0WboqIi5Ofnyx5ERERUvcxrOoFS165dQ0lJCZydnWXtzs7OOHPmjNZ11Gq11ni1Wi0tL20rL6aqZs+ejeeeew7W1tbYvXs33nrrLdy8eRNvv/12hbnk5+fj9u3bsLKyKtNnVFQUPvroI53yICIioodTa4qf2i4iIkL6f9u2bVFQUIAFCxZIxY8+pk+fjvDwcOl5fn4+3NzcHipPIiIiqlit+drL0dERZmZmyMrKkrVnZWVBqVRqXUepVFYYX/qvLn1Wlb+/P/79918UFRVVmIudnZ3Woz4AYGFhATs7O9mDiIiIqletKX4UCgV8fX0RHx8vtWk0GsTHxyMgIEDrOgEBAbJ4AIiLi5PiVSoVlEqlLCY/Px+JiYnl9llVycnJqFevHiwsLKqUCxEREdUOteprr/DwcISEhKB9+/bw8/PD4sWLUVBQgNDQUADAyJEj0ahRI0RFRQEAJk+ejG7dumHhwoXo168fNm7ciKNHj2LVqlUA7p2hFRYWho8//hheXl5QqVSIiIiAq6ur7KytjIwM5OTkICMjAyUlJUhOTgYAeHp6wtbWFj/99BOysrLQsWNHWFpaIi4uDp988ol0BhkAvPHGG1i+fDmmTp2K0aNHY+/evdi0aRN27tz5aAaPiIiIqqRWFT9DhgzB1atXERkZCbVaDR8fH+zatUuaSJyRkQFT0/9/sKpTp07YsGEDZsyYgQ8++ABeXl7Ytm0bWrZsKcVMnToVBQUFGD9+PHJzc9GlSxfs2rULlpaWUkxkZCTWrl0rPW/bti0AYN++fejevTvq1KmDFStW4J133oEQAp6envj8888xbtw4aR2VSoWdO3finXfewZIlS/DUU0/hq6++QlBQULWNFxEREemuVl3nx9jxOj9EVOvwOj/0mNDld+hDHfm5c+cO1Go1bt26BScnJ9SvX/9huiMiIiKqdjpPeP7vv/+wcuVKdOvWDXZ2dvDw8EDz5s3h5OQEd3d3jBs3DkeOHKmOXImIiIgemk7Fz+effw4PDw+sWbMGgYGB2LZtG5KTk/HXX38hISEBM2fOxN27d9GrVy/07t27zA1EiYiIiGqaTl97HTlyBL///ju8vb21Lvfz88Po0aMRHR2NNWvW4I8//oCXl5dBEiUiIiIyBJ2Kn++++65KcRYWFnjjjTf0SoiIiIioOtWaixwSERERPQp6Fz+3b9/GrVu3pOcXL17E4sWLsXv3boMkRkRERFQd9C5++vfvj3Xr1gEAcnNz4e/vj4ULF6J///5YuXKlwRIkIiIiMiS9i5+kpCR07doVALBlyxY4Ozvj4sWLWLduHZYuXWqwBImIiIgMSe/i59atW6hbty4AYPfu3Rg4cCBMTU3RsWNHXLx40WAJEhERERmS3sWPp6cntm3bhkuXLuHXX39Fr169AADZ2dmPza0ZiIiIyPjoXfxERkZiypQp8PDwgJ+fHwICAgDcOwpUemNQIiIiotpG73t7+fv7IyMjA5mZmfDx8ZHae/bsiYEDBxoiNyIiIiKD07v4cXd3R/369dGmTRv4+PhIDzMzM3zyySdYu3atIfMkIiIiMgi9i5/09HQcP34cycnJOH78ODZt2oQrV64AAOf8EBERUa31UEd+3N3dERwcLLUlJCQgJCQEs2fPNkRuRERERAZn0NtbBAQEYMmSJfjss88M2S0RERGRwehd/BQXF2tt9/LywunTp/VOiIiIiKg66f21l62tLVq0aIG2bdvCx8cHbdu2haurK5YtW4bAwEBD5khERERkMHoXP3v37sWJEydw4sQJrF+/HtOnT0dhYSEAoHfv3oiMjESrVq3QqlUrNGvWzGAJExERET0MvYufLl26oEuXLtJzjUaDs2fPIjk5GcnJyTh8+DC+/PJLZGdno6SkxCDJEhERET0svYufB5mamqJ58+Zo3rw5hg0bJrVnZWUZahNERERED82gZ3tp4+zsXN2bICIiIqqyai9+iIiIiGoTFj9ERERkVFj8EBERkVHRufiJjIzEsWPHqiMXIiIiomqnc/Hz77//ok+fPnjqqafw5ptv4pdffin3as9EREREtY3Oxc/q1auhVqvx3XffoW7duggLC4OjoyMGDRqEdevWIScnpzryJCIiIjIIveb8mJqaomvXrpg/fz7Onj2LxMRE+Pv744svvoCrqyueffZZfPbZZ7h8+bKh8yUiIiJ6KAaZ8Ny8eXNMnToVBw4cwKVLlxASEoI//vgD3333nSG6JyIiIjIYg13huZSTkxPGjBmDMWPGGLprIiIioofGU92JiIjIqLD4ISIiIqPC4oeIiIiMCosfIiIiMioGL34uXbqE0aNH673+ihUr4OHhAUtLS/j7++Pw4cMVxm/evBnNmjWDpaUlWrVqhZ9//lm2XAiByMhIuLi4wMrKCoGBgTh37pwsZu7cuejUqROsra3h4OBQZhsnTpzAsGHD4ObmBisrKzRv3hxLliyRxezfvx8mJiZlHmq1Wr+BICIiomph8OInJycHa9eu1Wvd2NhYhIeHY+bMmUhKSkKbNm0QFBSE7OxsrfEHDx7EsGHDMGbMGBw/fhzBwcEIDg5GSkqKFDN//nwsXboU0dHRSExMhI2NDYKCglBYWCjFFBcXY/DgwXjzzTe1bufYsWNo2LAhvv32W5w+fRoffvghpk+fjuXLl5eJPXv2LDIzM6VHw4YN9RoLIiIiqh4mQgihywrbt2+vcPk///yDd999FyUlJTon4+/vjw4dOkhFhUajgZubGyZNmoRp06aViR8yZAgKCgqwY8cOqa1jx47w8fFBdHQ0hBBwdXXFu+++iylTpgAA8vLy4OzsjJiYGAwdOlTWX0xMDMLCwpCbm1tprhMmTEBaWhr27t0L4N6Rnx49euDGjRtajx5VRX5+Puzt7ZGXlwc7Ozu9+njU2reXPz96tGbyIKJqom0nv7+NOz3VErr8DtX5Oj/BwcEwMTFBRTWTiYmJrt2iuLgYx44dw/Tp06U2U1NTBAYGIiEhQes6CQkJCA8Pl7UFBQVh27ZtAID09HSo1WoEBgZKy+3t7eHv74+EhIQyxY8u8vLyUL9+/TLtPj4+KCoqQsuWLTFr1ix07ty53D6KiopQVFQkPc/Pz9c7HyIiIqoanb/2cnFxwQ8//ACNRqP1kZSUpFci165dQ0lJCZydnWXtzs7O5c6bUavVFcaX/qtLn1Vx8OBBxMbGYvz48VKbi4sLoqOj8f333+P777+Hm5sbunfvXuF4REVFwd7eXnq4ubnpnRMRERFVjc7Fj6+vL44dO1bu8sqOCj3uUlJS0L9/f8ycORO9evWS2ps2bYrXX38dvr6+6NSpE1avXo1OnTph0aJF5fY1ffp05OXlSY9Lly49ipdARERk1HT+2uu9995DQUFBucs9PT2xb98+nRNxdHSEmZkZsrKyZO1ZWVlQKpVa11EqlRXGl/6blZUFFxcXWYyPj4/OOaampqJnz54YP348ZsyYUWm8n58f/vzzz3KXW1hYwMLCQuc8iIiISH86H/np2rUrevfuXe5yGxsbdOvWTedEFAoFfH19ER8fL7VpNBrEx8cjICBA6zoBAQGyeACIi4uT4lUqFZRKpSwmPz8fiYmJ5fZZntOnT6NHjx4ICQnB3Llzq7ROcnKyrOgiIiKimqfTkZ+MjAw0bty4yvGXL19Go0aNqhwfHh6OkJAQtG/fHn5+fli8eDEKCgoQGhoKABg5ciQaNWqEqKgoAMDkyZPRrVs3LFy4EP369cPGjRtx9OhRrFq1CsC9r+DCwsLw8ccfw8vLCyqVChEREXB1dUVwcLDsdeXk5CAjIwMlJSVITk4GcO8olq2tLVJSUvDcc88hKCgI4eHh0nwhMzMzODk5AQAWL14MlUoFb29vFBYW4quvvsLevXuxe/fuKr9+IiIiegSEDho2bCjGjx8vDh8+XG5Mbm6uWLVqlfD29hZLlizRpXshhBDLli0TjRs3FgqFQvj5+YlDhw5Jy7p16yZCQkJk8Zs2bRLPPPOMUCgUwtvbW+zcuVO2XKPRiIiICOHs7CwsLCxEz549xdmzZ2UxISEhAkCZx759+4QQQsycOVPrcnd3d6mPefPmiaefflpYWlqK+vXri+7du4u9e/fq9Nrz8vIEAJGXl6fTejXJ11f+IKInjLadnDs91UK6/A7V6To/169fx9y5c7F69WpYWlrC19cXrq6usLS0xI0bN5CamorTp0+jXbt2iIiIQN++fQ1dqz3ReJ0fIqp1eJ0fekzo8jtUpzk/DRo0wOeff47MzEwsX74cXl5euHbtmnS7iOHDh+PYsWNISEhg4UNERES1ks5newGAlZUVXn75Zbz88suGzoeIiIioWvGu7kRERGRUWPwQERGRUWHxQ0REREaFxQ8REREZFRY/REREZFT0Ln5u376NW7duSc8vXryIxYsX84rGREREVKvpXfz0798f69atAwDk5ubC398fCxcuRP/+/bFy5UqDJUhERERkSHoXP0lJSejatSsAYMuWLXB2dsbFixexbt06LF261GAJEhERERmS3sXPrVu3ULduXQDA7t27MXDgQJiamqJjx464ePGiwRIkIiIiMiS9ix9PT09s27YNly5dwq+//opevXoBALKzsx+b+1IRERGR8dG7+ImMjMSUKVPg4eEBPz8/BAQEALh3FKht27YGS5CIiIjIkPS6txcAvPzyy+jSpQsyMzPh4+Mjtffs2RMDBgwwRG5EREREBvdQ1/k5d+4cFi5ciM6dO+Py5csAgLNnz+LatWsGSY6IiIjI0PQufr7//nsEBQXBysoKSUlJKCoqAgDk5eXhk08+MViCRERERIakd/Hz8ccfIzo6Gl9++SXq1KkjtXfu3BlJSUkGSY6IiIjI0PQufs6ePYtnn322TLu9vT1yc3MfJiciIiKiaqN38aNUKnH+/Pky7X/++SeaNGnyUEkRERERVRe9i59x48Zh8uTJSExMhImJCa5cuYL169djypQpePPNNw2ZIxEREZHB6H2q+7Rp06DRaNCzZ0/cunULzz77LCwsLDBlyhRMmjTJkDkSERERGYzexY+JiQk+/PBDvPfeezh//jxu3ryJFi1awNbW1pD5ERERERmU3sVPRkYG3NzcoFAo0KJFizLLGjdu/NDJERERERma3nN+VCoVrl69Wqb9+vXrUKlUD5UUERERUXXRu/gRQsDExKRM+82bN2FpaflQSRERERFVF52/9goPDwdwb85PREQErK2tpWUlJSVITEyU3euLiIiIqDbRufg5fvw4gHtHfk6dOgWFQiEtUygUaNOmDaZMmWK4DImIiIgMSOfiZ9++fQCA0NBQLFmyBHZ2dgZPioiIiKi66H2215o1awAAqampyMjIQHFxsWz5Sy+99HCZEREREVUDvYuf9PR0BAcH49SpUzAxMYEQAgCkSdAlJSWGyZCIiIjIgPQ+2+vtt9+GSqVCdnY2rK2tcfr0afz+++9o37499u/fb8AUiYiIiAxH7yM/CQkJ2Lt3LxwdHWFqagpTU1N06dIFUVFRePvtt6WJ0URERES1id5HfkpKSlC3bl0AgKOjI65cuQIAcHd3x9mzZw2THREREZGB6X3kp2XLljhx4gRUKhX8/f0xf/58KBQKrFq1Ck2aNDFkjkREREQGo3fxM2PGDBQUFAAAZs+ejRdeeAFdu3ZFgwYNEBsba7AEiYiIiAxJ7+InKChI+r+npyfOnDmDnJwc1KtXT+ttL4iIiIhqA73n/GhTv359mJiY4Pbt23r3sWLFCnh4eMDS0hL+/v44fPhwhfGbN29Gs2bNYGlpiVatWuHnn3+WLRdCIDIyEi4uLrCyskJgYCDOnTsni5k7dy46deoEa2trODg4aN1ORkYG+vXrB2trazRs2BDvvfce7t69K4vZv38/2rVrBwsLC3h6eiImJkbn109ERETVy6DFT1FRERYuXKj3Xd1jY2MRHh6OmTNnIikpCW3atEFQUBCys7O1xh88eBDDhg3DmDFjcPz4cQQHByM4OBgpKSlSzPz587F06VJER0cjMTERNjY2CAoKQmFhoRRTXFyMwYMH480339S6nZKSEvTr1w/FxcU4ePAg1q5di5iYGERGRkox6enp6NevH3r06IHk5GSEhYVh7Nix+PXXX/UaCyIiIqomQkeFhYVi2rRpwtfXVwQEBIitW7cKIYRYvXq1cHFxEU899ZT49NNPde1WCCGEn5+fmDBhgvS8pKREuLq6iqioKK3xr7zyiujXr5+szd/fX7z++utCCCE0Go1QKpViwYIF0vLc3FxhYWEhvvvuuzL9rVmzRtjb25dp//nnn4WpqalQq9VS28qVK4WdnZ0oKioSQggxdepU4e3tLVtvyJAhIigoqJJX/f/l5eUJACIvL6/K69Q0X1/5g4ieMNp2cu70VAvp8jtU5yM/kZGRWLlyJTw8PHDhwgUMHjwY48ePx6JFi/D555/jwoULeP/993UuwoqLi3Hs2DEEBgZKbaampggMDERCQoLWdRISEmTxwL25SKXx6enpUKvVshh7e3v4+/uX22d522nVqhWcnZ1l28nPz8fp06erlIs2RUVFyM/Plz2IiIioeuk84Xnz5s1Yt24dXnrpJaSkpKB169a4e/cuTpw48VATna9du4aSkhJZgQEAzs7OOHPmjNZ11Gq11ni1Wi0tL20rL6YqytvO/dsoLyY/Px+3b9+GlZVVmX6joqLw0UcfVTkPIiIieng6H/n5999/4evrC+DetX4sLCzwzjvv8AwvPUyfPh15eXnS49KlSzWdEhER0RNP5+KnpKQECoVCem5ubg5bW9uHTsTR0RFmZmbIysqStWdlZUGpVGpdR6lUVhhf+q8ufeqynfu3UV6MnZ2d1qM+AGBhYQE7OzvZg4iIiKqXzsWPEAKjRo3CwIEDMXDgQBQWFuKNN96Qnpc+dKVQKODr64v4+HipTaPRID4+HgEBAVrXCQgIkMUDQFxcnBSvUqmgVCplMfn5+UhMTCy3z/K2c+rUKdlZZ3FxcbCzs0OLFi2qlAsRERHVDjrP+QkJCZE9f+211wyWTHh4OEJCQtC+fXv4+flh8eLFKCgoQGhoKABg5MiRaNSoEaKiogAAkydPRrdu3bBw4UL069cPGzduxNGjR7Fq1SoAgImJCcLCwvDxxx/Dy8sLKpUKERERcHV1RXBwsLTdjIwM5OTkICMjAyUlJUhOTgZw7+KNtra26NWrF1q0aIERI0Zg/vz5UKvVmDFjBiZMmAALCwsAwBtvvIHly5dj6tSpGD16NPbu3YtNmzZh586dBhsfIiIiMoDqP/lMN8uWLRONGzcWCoVC+Pn5iUOHDknLunXrJkJCQmTxmzZtEs8884xQKBTC29tb7Ny5U7Zco9GIiIgI4ezsLCwsLETPnj3F2bNnZTEhISECQJnHvn37pJgLFy6IPn36CCsrK+Ho6CjeffddcefOHVk/+/btEz4+PkKhUIgmTZqINWvW6PTaeao7EdU6PNWdHhO6/A41EUKIGqy96D75+fmwt7dHXl7eYzP/p317+fOjR2smDyKqJtp28vvbuNNTLaHL71CDXuGZiIiIqLZj8UNERERGhcUPERERGRUWP0RERGRUdD7V/X7x8fGIj49HdnY2NBqNbNnq1asfKjEiIiKi6qB38fPRRx9h9uzZaN++PVxcXHh7CyIiInos6F38REdHIyYmBiNGjDBkPkRERETVSu85P8XFxejUqZMhcyEiIiKqdnoXP2PHjsWGDRsMmQsRERFRtdP7a6/CwkKsWrUKe/bsQevWrVGnTh3Z8s8///yhkyMiIiIyNL2Ln5MnT8LHxwcAkJKSIlvGyc9ERERUW+ld/Ozbt8+QeRARERE9ErzIIRERERmVh7rIYW5uLr7++mukpaUBAFq0aIExY8bA3t7eIMkRERERGZreR36OHj2Kp59+GosWLUJOTg5ycnKwaNEiPP3000hKSjJkjkREREQGo/eRn3feeQcvvfQSvvzyS5ib3+vm7t27GDt2LMLCwvD7778bLEkiIiIiQ9G7+Dl69Kis8AEAc3NzTJ06Fe3btzdIckRERESGpvfXXnZ2dsjIyCjTfunSJdStW/ehkiIiIiKqLnoXP0OGDMGYMWMQGxuLS5cu4dKlS9i4cSPGjh2LYcOGGTJHIiIiIoPR+2uvzz77DCYmJhg5ciTu3r0LIQQUCgXefPNNfPrpp4bMkYiIiMhg9C5+FAoFlixZgqioKPz9998AgKeffhrW1tYGS46IiIjI0HQqfsLDwzFnzhzY2NggPDy8wlje24uIiIhqI52Kn+PHj+POnTvS/8vDe3sRERFRbaVT8XP//bx4by8iIiJ6HOl9tldGRgaEEOUuIyIiIqqN9C5+VCoVrl69Wqb9+vXrUKlUD5UUERERUXXRu/gRQmid23Pz5k1YWlo+VFJERERE1UXnU91Lz/IyMTFBRESE7NT2kpISJCYmwsfHx2AJEhERERmSzsVP6VleQgicOnUKCoVCWqZQKNCmTRtMmTLFcBkSERERGZDOxU/pWV6hoaFYsmQJ7OzsDJ4UERERUXXR+wrPa9asMWQeRERERI8Er/BMRERERoVXeCYiIiKjwis8ExERkVHRe87P7du3IYSQTnW/ePEitm7dihYtWqBXr14GS5AeP+3by58fPVozeRAREWmj90UO+/fvj3Xr1gEAcnNz4efnh4ULF6J///5YuXKlwRIkIiIiMiS9i5+kpCR07doVALBlyxYolUpcvHgR69atw9KlSx8qqRUrVsDDwwOWlpbw9/fH4cOHK4zfvHkzmjVrBktLS7Rq1Qo///yzbLkQApGRkXBxcYGVlRUCAwNx7tw5WUxOTg6GDx8OOzs7ODg4YMyYMbh586a0fNasWTAxMSnzsLGxkWJiYmLKLOfVromIiGoXvYufW7duoW7dugCA3bt3Y+DAgTA1NUXHjh1x8eJFvROKjY1FeHg4Zs6ciaSkJLRp0wZBQUHIzs7WGn/w4EEMGzYMY8aMwfHjxxEcHIzg4GCkpKRIMfPnz8fSpUsRHR2NxMRE2NjYICgoCIWFhVLM8OHDcfr0acTFxWHHjh34/fffMX78eGn5lClTkJmZKXu0aNECgwcPluVjZ2cni3mYsSAiIqJqIPTUqlUrsWTJEpGRkSHs7OzEwYMHhRBCHD16VDg7O+vbrfDz8xMTJkyQnpeUlAhXV1cRFRWlNf6VV14R/fr1k7X5+/uL119/XQghhEajEUqlUixYsEBanpubKywsLMR3330nhBAiNTVVABBHjhyRYn755RdhYmIiLl++rHW7ycnJAoD4/fffpbY1a9YIe3t73V7wffLy8gQAkZeXp3cfj5qvr/xRXhsRPaYq28mJagldfofqfeQnMjISU6ZMgYeHB/z9/REQEADg3lGgtm3b6tVncXExjh07hsDAQKnN1NQUgYGBSEhI0LpOQkKCLB4AgoKCpPj09HSo1WpZjL29Pfz9/aWYhIQEODg4oP19M3UDAwNhamqKxMRErdv96quv8Mwzz0hf/ZW6efMm3N3d4ebmhv79++P06dPlvt6ioiLk5+fLHkRERFS99C5+Xn75ZWRkZODo0aPYtWuX1N6zZ08sWrRIrz6vXbuGkpISODs7y9qdnZ2hVqu1rqNWqyuML/23spiGDRvKlpubm6N+/fpat1tYWIj169djzJgxsvamTZti9erV+PHHH/Htt99Co9GgU6dO+Pfff7XmHhUVBXt7e+nh5uamNY6IiIgMR+9T3QFAqVRCqVTK2vz8/B4qocfB1q1b8d9//yEkJETWHhAQIB0BA4BOnTqhefPm+OKLLzBnzpwy/UyfPl12pez8/HwWQERERNXsoYqf+Ph4xMfHIzs7GxqNRrZs9erVOvfn6OgIMzMzZGVlydqzsrLKFFmllEplhfGl/2ZlZcHFxUUW4+PjI8U8OKH67t27yMnJ0brdr776Ci+88EKZo0kPqlOnDtq2bYvz589rXW5hYQELC4sK+yAiIiLD0vtrr48++gi9evVCfHw8rl27hhs3bsge+lAoFPD19UV8fLzUptFoEB8fLzuicr+AgABZPADExcVJ8SqVCkqlUhaTn5+PxMREKSYgIAC5ubk4duyYFLN3715oNBr4+/vL+k5PT8e+ffvKfOWlTUlJCU6dOiUruoiIiKhm6X3kJzo6GjExMRgxYoQh80F4eDhCQkLQvn17+Pn5YfHixSgoKEBoaCgAYOTIkWjUqBGioqIAAJMnT0a3bt2wcOFC9OvXDxs3bsTRo0exatUqAPfuMxYWFoaPP/4YXl5eUKlUiIiIgKurK4KDgwEAzZs3R+/evTFu3DhER0fjzp07mDhxIoYOHQpXV1dZfqtXr4aLiwv69OlTJvfZs2ejY8eO8PT0RG5uLhYsWICLFy9i7NixBh0jIiIi0p/exU9xcTE6depkyFwAAEOGDMHVq1cRGRkJtVoNHx8f7Nq1S/qKKSMjA6am//+AVadOnbBhwwbMmDEDH3zwAby8vLBt2za0bNlSipk6dSoKCgowfvx45ObmokuXLti1a5fsAoTr16/HxIkT0bNnT5iammLQoEFlLtao0WgQExODUaNGwczMrEzuN27cwLhx46BWq1GvXj34+vri4MGDaNGihaGHiYiIiPRkIoQQ+qz4/vvvw9bWFhEREYbOyWjl5+fD3t4eeXl5sLOzq+l0qkTbfbx4by+iJ0hlOzl3cKoldPkdqveRn8LCQqxatQp79uxB69atUadOHdnyzz//XN+uiYiIiKqN3sXPyZMnpbOl7r+VBHBvng0RERFRbaR38bNv3z5D5kFERET0SDzUdX4AIDU1FRkZGSguLpbaTExM8OKLLz5s10REREQGp3fx888//2DAgAE4deoUTExMUDpvuvQrr5KSEsNkSERERGRAel/kcPLkyVCpVMjOzoa1tTVOnz6N33//He3bt8f+/fsNmCIRERGR4eh95CchIQF79+6Fo6MjTE1NYWpqii5duiAqKgpvv/02jh8/bsg8iYiIiAxC7yM/JSUlqFu3LoB79+S6cuUKAMDd3R1nz541THZEREREBqb3kZ+WLVvixIkTUKlU8Pf3x/z586FQKLBq1So0adLEkDkSERERGYzexc+MGTNQUFAA4N49rV544QV07doVDRo0QGxsrMESJCIiIjIkvYufoKAg6f+enp44c+YMcnJyUK9ePV7kkIiIiGotveb83LlzBz179sS5c+dk7fXr12fhQ0RERLWaXsVPnTp1cPLkSUPnQkRERFTt9D7b67XXXsPXX39tyFyIiIiIqp3ec37u3r2L1atXY8+ePfD19YWNjY1sOe/qTkRERLWR3sVPSkoK2rVrBwD466+/ZMs474eIiIhqK52Ln9mzZ2PKlCm8qzsRERE9lnSe8/PRRx/h5s2b1ZELERERUbXTufgpvXs7ERER0eNIr7O9OKeHiIiIHld6TXh+5plnKi2AcnJy9EqIiIiIqDrpVfx89NFHsLe3N3QuRERERNVOr+Jn6NChaNiwoaFzISIiIqp2Os/54XwfIiIiepzxbC8iIiIyKjp/7aXRaKojDyIiIqJHQu8bmxIRERE9jlj8EBERkVFh8UNERERGhcUPERERGRUWP0RERGRUWPwQERGRUWHxQ0REREaFxQ8REREZFRY/REREZFRY/BAREZFRqZXFz4oVK+Dh4QFLS0v4+/vj8OHDFcZv3rwZzZo1g6WlJVq1aoWff/5ZtlwIgcjISLi4uMDKygqBgYE4d+6cLCYnJwfDhw+HnZ0dHBwcMGbMGNy8eVNafuHCBZiYmJR5HDp0SKdciIiIqGbVuuInNjYW4eHhmDlzJpKSktCmTRsEBQUhOztba/zBgwcxbNgwjBkzBsePH0dwcDCCg4ORkpIixcyfPx9Lly5FdHQ0EhMTYWNjg6CgIBQWFkoxw4cPx+nTpxEXF4cdO3bg999/x/jx48tsb8+ePcjMzJQevr6+OuVCRERENctE1LLbtPv7+6NDhw5Yvnw5gHs3UnVzc8OkSZMwbdq0MvFDhgxBQUEBduzYIbV17NgRPj4+iI6OhhACrq6uePfddzFlyhQAQF5eHpydnRETE4OhQ4ciLS0NLVq0wJEjR9C+fXsAwK5du9C3b1/8+++/cHV1xYULF6BSqXD8+HH4+Phozb2yXCqTn58Pe3t75OXlwc7OrspjVpP+b7gkR49qbyOix1RlOzl3cKoldPkdWquO/BQXF+PYsWMIDAyU2kxNTREYGIiEhASt6yQkJMjiASAoKEiKT09Ph1qtlsXY29vD399fiklISICDg4NU+ABAYGAgTE1NkZiYKOv7pZdeQsOGDdGlSxds375dp1weVFRUhPz8fNmDiIiIqletKn6uXbuGkpISODs7y9qdnZ2hVqu1rqNWqyuML/23spiGDRvKlpubm6N+/fpSjK2tLRYuXIjNmzdj586d6NKlC4KDg2UFUGW5PCgqKgr29vbSw83NTWscERERGY55TSfwuHB0dER4eLj0vEOHDrhy5QoWLFiAl156Sa8+p0+fLuszPz+fBRAREVE1q1VHfhwdHWFmZoasrCxZe1ZWFpRKpdZ1lEplhfGl/1YW8+CE6rt37yInJ6fc7QL35iedP3++yrk8yMLCAnZ2drIHERERVa9aVfwoFAr4+voiPj5eatNoNIiPj0dAQIDWdQICAmTxABAXFyfFq1QqKJVKWUx+fj4SExOlmICAAOTm5uLYsWNSzN69e6HRaODv719uvsnJyXBxcalyLkRERFTzat3XXuHh4QgJCUH79u3h5+eHxYsXo6CgAKGhoQCAkSNHolGjRoiKigIATJ48Gd26dcPChQvRr18/bNy4EUePHsWqVasAACYmJggLC8PHH38MLy8vqFQqREREwNXVFcHBwQCA5s2bo3fv3hg3bhyio6Nx584dTJw4EUOHDoWrqysAYO3atVAoFGjbti0A4IcffsDq1avx1VdfSblXlgsRERHVvFpX/AwZMgRXr15FZGQk1Go1fHx8sGvXLmkicUZGBkxN//8Bq06dOmHDhg2YMWMGPvjgA3h5eWHbtm1o2bKlFDN16lQUFBRg/PjxyM3NRZcuXbBr1y5YWlpKMevXr8fEiRPRs2dPmJqaYtCgQVi6dKkstzlz5uDixYswNzdHs2bNEBsbi5dfflmnXIiIiKhm1brr/BgzXueHiGodXueHHhOP7XV+iIiIiKobix8iIiIyKix+iIiIyKiw+CEiIiKjwuKHiIiIjAqLHyIiIjIqLH6IiIjIqLD4ISIiIqPC4oeIiIiMSq27vQU9mXhBWCIiqi145IeIiIiMCosfIiIiMiosfoiIiMiosPghIiIio8Lih4iIiIwKix8iIiIyKix+iIiIyKiw+CEiIiKjwuKHiIiIjAqLHyIiIjIqLH6IiIjIqLD4ISIiIqPC4oeIiIiMCosfIiIiMirmNZ0AERHVEu3by58fPVozeRBVMx75ISIiIqPC4oeIiIiMCosfIiIiMiosfoiIiMiosPghIiIio8Lih4iIiIwKix8iIiIyKix+iIiIyKiw+CEiIiKjwuKHiIiIjAqLHyIiIjIqtbL4WbFiBTw8PGBpaQl/f38cPny4wvjNmzejWbNmsLS0RKtWrfDzzz/LlgshEBkZCRcXF1hZWSEwMBDnzp2TxeTk5GD48OGws7ODg4MDxowZg5s3b0rL9+/fj/79+8PFxQU2Njbw8fHB+vXrZX3ExMTAxMRE9rC0tHzI0SAiIiJDqnXFT2xsLMLDwzFz5kwkJSWhTZs2CAoKQnZ2ttb4gwcPYtiwYRgzZgyOHz+O4OBgBAcHIyUlRYqZP38+li5diujoaCQmJsLGxgZBQUEoLCyUYoYPH47Tp08jLi4OO3bswO+//47x48fLttO6dWt8//33OHnyJEJDQzFy5Ejs2LFDlo+dnR0yMzOlx8WLFw08QkRERPQwTIQQoqaTuJ+/vz86dOiA5cuXAwA0Gg3c3NwwadIkTJs2rUz8kCFDUFBQICtCOnbsCB8fH0RHR0MIAVdXV7z77ruYMmUKACAvLw/Ozs6IiYnB0KFDkZaWhhYtWuDIkSNo/393Nd61axf69u2Lf//9F66urlpz7devH5ydnbF69WoA9478hIWFITc3V6/Xnp+fD3t7e+Tl5cHOzk6vPh41bTeBrqyNN4omqqX03aF5N3iqBXT5HVqrjvwUFxfj2LFjCAwMlNpMTU0RGBiIhIQEreskJCTI4gEgKChIik9PT4darZbF2Nvbw9/fX4pJSEiAg4ODVPgAQGBgIExNTZGYmFhuvnl5eahfv76s7ebNm3B3d4ebmxv69++P06dPl7t+UVER8vPzZQ9j0b69/EFERPSo1Kri59q1aygpKYGzs7Os3dnZGWq1Wus6arW6wvjSfyuLadiwoWy5ubk56tevX+52N23ahCNHjiA0NFRqa9q0KVavXo0ff/wR3377LTQaDTp16oR///1Xax9RUVGwt7eXHm5ublrjiIiIyHBqVfHzuNi3bx9CQ0Px5ZdfwtvbW2oPCAjAyJEj4ePjg27duuGHH36Ak5MTvvjiC639TJ8+HXl5edLj0qVLj+olEBERGa1aVfw4OjrCzMwMWVlZsvasrCwolUqt6yiVygrjS/+tLObBCdV3795FTk5Ome3+9ttvePHFF7Fo0SKMHDmywtdTp04dtG3bFufPn9e63MLCAnZ2drIHERERVa9aVfwoFAr4+voiPj5eatNoNIiPj0dAQIDWdQICAmTxABAXFyfFq1QqKJVKWUx+fj4SExOlmICAAOTm5uLYsWNSzN69e6HRaODv7y+17d+/H/369cO8efNkZ4KVp6SkBKdOnYKLi0sVXj0RERE9CuY1ncCDwsPDERISgvbt28PPzw+LFy9GQUGBNLdm5MiRaNSoEaKiogAAkydPRrdu3bBw4UL069cPGzduxNGjR7Fq1SoAgImJCcLCwvDxxx/Dy8sLKpUKERERcHV1RXBwMACgefPm6N27N8aNG4fo6GjcuXMHEydOxNChQ6Uzvfbt24cXXngBkydPxqBBg6S5QAqFQpr0PHv2bHTs2BGenp7Izc3FggULcPHiRYwdO/ZRDiERERFVoNYVP0OGDMHVq1cRGRkJtVoNHx8f7Nq1S5qwnJGRAVPT/3/AqlOnTtiwYQNmzJiBDz74AF5eXti2bRtatmwpxUydOhUFBQUYP348cnNz0aVLF+zatUt2AcL169dj4sSJ6NmzJ0xNTTFo0CAsXbpUWr527VrcunULUVFRUuEFAN26dcP+/fsBADdu3MC4ceOgVqtRr149+Pr64uDBg2jRokV1DRcRERHpqNZd58eYGdN1fnhZEKJaiDs0PcYe2+v8EBEREVU3Fj9ERERkVFj8EBERkVFh8UNERERGhcUPERERGRUWP0RERGRUWPwQERGRUWHxQ0REREal1l3hmYwXr5NGRESPAo/8EBERkVHhkR8iImP14OFWIiPBIz9ERERkVFj8EBERkVFh8UNERERGhcUPERERGRUWP0RERGRUWPwQERGRUWHxQ0REREaF1/mhWotXfCYiourA4oeIiAyPf71QLcavvYiIiMio8MgPEZEx4JEYIgmP/BAREZFR4ZEfeqzwj1ciInpYPPJDRERERoXFDxERERkVFj9ERERkVDjnhx57988D4hwgov9TG3eM2pgTGSUe+SEiIiKjwuKHiIiIjAq/9qInDk+HJ6PDDz2RTlj8EBFRzWDRRjWExQ8ZBf6MpScGP8xED43FDxkl/v6gx8aDH1Yiemgsfoj+DwsiqnH8EPJ0eHokWPwQVYA/h8lgtBU2/IBVjgUhVYNaear7ihUr4OHhAUtLS/j7++Pw4cMVxm/evBnNmjWDpaUlWrVqhZ9//lm2XAiByMhIuLi4wMrKCoGBgTh37pwsJicnB8OHD4ednR0cHBwwZswY3Lx5UxZz8uRJdO3aFZaWlnBzc8P8+fN1zoUeb+3byx/ltdETrCofAn4oqhfHlx5SrTvyExsbi/DwcERHR8Pf3x+LFy9GUFAQzp49i4YNG5aJP3jwIIYNG4aoqCi88MIL2LBhA4KDg5GUlISWLVsCAObPn4+lS5di7dq1UKlUiIiIQFBQEFJTU2FpaQkAGD58ODIzMxEXF4c7d+4gNDQU48ePx4YNGwAA+fn56NWrFwIDAxEdHY1Tp05h9OjRcHBwwPjx46ucCxmHB/+g5x/91UzXAdeGb9TjrSqfAaL/U+uKn88//xzjxo1DaGgoACA6Oho7d+7E6tWrMW3atDLxS5YsQe/evfHee+8BAObMmYO4uDgsX74c0dHREEJg8eLFmDFjBvr37w8AWLduHZydnbFt2zYMHToUaWlp2LVrF44cOYL2/7fDLFu2DH379sVnn30GV1dXrF+/HsXFxVi9ejUUCgW8vb2RnJyMzz//XCp+KsuFqCJV+X2t7XextphUm//f2KKgisVAVTvXp9Co7r6JqqI6P7v0WKlVxU9xcTGOHTuG6dOnS22mpqYIDAxEQkKC1nUSEhIQHh4uawsKCsK2bdsAAOnp6VCr1QgMDJSW29vbw9/fHwkJCRg6dCgSEhLg4OAgFT4AEBgYCFNTUyQmJmLAgAFISEjAs88+C4VCIdvOvHnzcOPGDdSrV6/SXB5UVFSEoqIi6XleXh6Ae0eZHhclJfLn+fmVt1UlRhtD9v3F2W4PxPyG/6VW3FaVGH3XM3TfN0XJfc/zcSZVPgDNqtBmqJja2Lc27Lvqfeu10+mz3uPUd7f79sPffpM/r2qboWIeZr0HVXffBlT6u1MIUWlsrSp+rl27hpKSEjg7O8vanZ2dcebMGa3rqNVqrfFqtVpaXtpWUcyDX6mZm5ujfv36shiVSlWmj9Jl9erVqzSXB0VFReGjjz4q0+7m5qY1/nFgb195W1Viqrtv/+rsXJ/12Pej7Vsb9l29fde2zwD7rp2fLwP477//YF9J37Wq+DE206dPlx0p0mg0yMnJQYMGDWBiYmLQbeXn58PNzQ2XLl2CnZ2dQfumsjjejxbH+9HjmD9aHO/KCSHw33//wdXVtdLYWlX8ODo6wszMDFlZWbL2rKwsKJVKresolcoK40v/zcrKgouLiyzGx8dHisnOzpb1cffuXeTk5Mj60bad+7dRWS4PsrCwgIWFhazNwcFBa6yh2NnZccd5hDjejxbH+9HjmD9aHO+KVXbEp1StOtVdoVDA19cX8fHxUptGo0F8fDwCAgK0rhMQECCLB4C4uDgpXqVSQalUymLy8/ORmJgoxQQEBCA3NxfHjh2TYvbu3QuNRgN/f38p5vfff8edO3dk22natCnq1atXpVyIiIioFhC1zMaNG4WFhYWIiYkRqampYvz48cLBwUGo1WohhBAjRowQ06ZNk+IPHDggzM3NxWeffSbS0tLEzJkzRZ06dcSpU6ekmE8//VQ4ODiIH3/8UZw8eVL0799fqFQqcfv2bSmmd+/eom3btiIxMVH8+eefwsvLSwwbNkxanpubK5ydncWIESNESkqK2Lhxo7C2thZffPGFTrnUlLy8PAFA5OXl1XQqRoHj/WhxvB89jvmjxfE2rFpX/AghxLJly0Tjxo2FQqEQfn5+4tChQ9Kybt26iZCQEFn8pk2bxDPPPCMUCoXw9vYWO3fulC3XaDQiIiJCODs7CwsLC9GzZ09x9uxZWcz169fFsGHDhK2trbCzsxOhoaHiv//+k8WcOHFCdOnSRVhYWIhGjRqJTz/9tEzuleVSUwoLC8XMmTNFYWFhTadiFDjejxbH+9HjmD9aHG/DMhGiCueEERERET0hatWcHyIiIqLqxuKHiIiIjAqLHyIiIjIqLH6IiIjIqLD4MRIrVqyAh4cHLC0t4e/vj8OHD9d0Sk+EqKgodOjQAXXr1kXDhg0RHByMs2fPymIKCwsxYcIENGjQALa2thg0aFCZi2GS7j799FOYmJggLCxMauNYG97ly5fx2muvoUGDBrCyskKrVq1w9L4beQohEBkZCRcXF1hZWSEwMBDnzp2rwYwfXyUlJYiIiIBKpYKVlRWefvppzJkzR3avKo63YbD4MQKxsbEIDw/HzJkzkZSUhDZt2iAoKKjMVa1Jd7/99hsmTJiAQ4cOIS4uDnfu3EGvXr1QUFAgxbzzzjv46aefsHnzZvz222+4cuUKBg4cWINZP/6OHDmCL774Aq1bt5a1c6wN68aNG+jcuTPq1KmDX375BampqVi4cKF0YVcAmD9/PpYuXYro6GgkJibCxsYGQUFBKCwsrMHMH0/z5s3DypUrsXz5cqSlpWHevHmYP38+li1bJsVwvA2kRk+0p0fCz89PTJgwQXpeUlIiXF1dRVRUVA1m9WTKzs4WAMRvv/0mhLh3ccw6deqIzZs3SzFpaWkCgEhISKipNB9r//33n/Dy8hJxcXGiW7duYvLkyUIIjnV1eP/990WXLl3KXa7RaIRSqRQLFiyQ2nJzc4WFhYX47rvvHkWKT5R+/fqJ0aNHy9oGDhwohg8fLoTgeBsSj/w84YqLi3Hs2DEEBgZKbaampggMDERCQkINZvZkysvLAwDUr18fAHDs2DHcuXNHNv7NmjVD48aNOf56mjBhAvr16ycbU4BjXR22b9+O9u3bY/DgwWjYsCHatm2LL7/8Ulqenp4OtVotG3N7e3v4+/tzzPXQqVMnxMfH46+//gIAnDhxAn/++Sf69OkDgONtSLXqxqZkeNeuXUNJSQmcnZ1l7c7Ozjhz5kwNZfVk0mg0CAsLQ+fOndGyZUsAgFqthkKhKHPDWmdnZ6jV6hrI8vG2ceNGJCUl4ciRI2WWcawN759//sHKlSsRHh6ODz74AEeOHMHbb78NhUKBkJAQaVy1/XzhmOtu2rRpyM/PR7NmzWBmZoaSkhLMnTsXw4cPBwCOtwGx+CEykAkTJiAlJQV//vlnTafyRLp06RImT56MuLg4WFpa1nQ6RkGj0aB9+/b45JNPAABt27ZFSkoKoqOjERISUsPZPXk2bdqE9evXY8OGDfD29kZycjLCwsLg6urK8TYwfu31hHN0dISZmVmZM16ysrKgVCprKKsnz8SJE7Fjxw7s27cPTz31lNSuVCpRXFyM3NxcWTzHX3fHjh1DdnY22rVrB3Nzc5ibm+O3337D0qVLYW5uDmdnZ461gbm4uKBFixaytubNmyMjIwMApHHlzxfDeO+99zBt2jQMHToUrVq1wogRI/DOO+8gKioKAMfbkFj8POEUCgV8fX0RHx8vtWk0GsTHxyMgIKAGM3syCCEwceJEbN26FXv37oVKpZIt9/X1RZ06dWTjf/bsWWRkZHD8ddSzZ0+cOnUKycnJ0qN9+/YYPny49H+OtWF17ty5zKUb/vrrL7i7uwMAVCoVlEqlbMzz8/ORmJjIMdfDrVu3YGoq/7VsZmYGjUYDgONtUDU945qq38aNG4WFhYWIiYkRqampYvz48cLBwUGo1eqaTu2x9+abbwp7e3uxf/9+kZmZKT1u3bolxbzxxhuicePGYu/eveLo0aMiICBABAQE1GDWT477z/YSgmNtaIcPHxbm5uZi7ty54ty5c2L9+vXC2tpafPvtt1LMp59+KhwcHMSPP/4oTp48Kfr37y9UKpW4fft2DWb+eAoJCRGNGjUSO3bsEOnp6eKHH34Qjo6OYurUqVIMx9swWPwYiWXLlonGjRsLhUIh/Pz8xKFDh2o6pScCAK2PNWvWSDG3b98Wb731lqhXr56wtrYWAwYMEJmZmTWX9BPkweKHY214P/30k2jZsqWwsLAQzZo1E6tWrZIt12g0IiIiQjg7OwsLCwvRs2dPcfbs2RrK9vGWn58vJk+eLBo3biwsLS1FkyZNxIcffiiKioqkGI63YZgIcd+lI4mIiIiecJzzQ0REREaFxQ8REREZFRY/REREZFRY/BAREZFRYfFDRERERoXFDxERERkVFj9ERERkVFj8EBERkVFh8UNERERGhcUPERERGRUWP0RPqO7duyMsLKzW9lddfVL14/tGjzsWP0S13KhRo2BiYgITExMoFAp4enpi9uzZuHv3boXr/fDDD5gzZ47B8jB0f1V16dIljB49Gq6urlAoFHB3d8fkyZNx/fr1R54LUPO/+Es/D59++qmsfdu2bTAxMamhrIgeLyx+iB4DvXv3RmZmJs6dO4d3330Xs2bNwoIFC7TGFhcXAwDq16+PunXrGiwHQ/dXFf/88w/at2+Pc+fO4bvvvsP58+cRHR2N+Ph4BAQEICcn55HmU1tYWlpi3rx5uHHjRk2nYjCln1uiR4HFD9FjwMLCAkqlEu7u7njzzTcRGBiI7du3A7h3JGLixIkICwuDo6MjgoKCpPbSIxTdu3fH22+/jalTp6J+/fpQKpWYNWuWbBsajQbz58+Hp6cnLCws0LhxY8ydO1da/uARj9LtTpw4Efb29nB0dERERASEEFLMrl270KVLFzg4OKBBgwZ44YUX8Pfff1f5dU+YMAEKhQK7d+9Gt27d0LhxY/Tp0wd79uzB5cuX8eGHH0qxHh4eWLx4sWx9Hx8f2eusLJ/KxmnUqFH47bffsGTJEulo3IULF6q8/e7du2PSpEkICwtDvXr14OzsjC+//BIFBQUIDQ1F3bp14enpiV9++aXCcQkMDIRSqURUVFSFcdWZ0927dyt87zUaDaKioqBSqWBlZYU2bdpgy5Ytsu1q+9xq8/fff8PExAQ7duxAz549YW1tjaZNmyIxMbHC109UHhY/RI8hKysr2V/Ka9euhUKhwIEDBxAdHa11nbVr18LGxgaJiYmYP38+Zs+ejbi4OGn59OnT8emnnyIiIgKpqanYsGEDnJ2dK8xj7dq1MDc3x+HDh7FkyRJ8/vnn+Oqrr6TlBQUFCA8Px9GjRxEfHw9TU1MMGDAAGo2m0teYk5ODX3/9FW+99RasrKxky5RKJYYPH47Y2FjZL9zKVCWfisZpyZIlCAgIwLhx45CZmYnMzEy4ublVeful/Ts6OuLw4cOYNGkS3nzzTQwePBidOnVCUlISevXqhREjRuDWrVvl9mFmZoZPPvkEy5Ytw7///qvT9g2VU2XvfVRUFNatW4fo6GicPn0a77zzDl577TX89ttvsj4q+9wCwIkTJ2BiYoLPP/8cEREROHHiBBo3boxp06Y99GsnIyWIqFYLCQkR/fv3F0IIodFoRFxcnLCwsBBTpkwRQgjRrVs30bZt2zLrdevWTUyePFn6f5cuXWTLO3ToIN5//30hhBD5+fnCwsJCfPnll+XmcX9/pc+bN28uNBqN1Pb++++L5s2bl9vH1atXBQBx6tQprX3e79ChQwKA2Lp1q9bln3/+uQAgsrKyhBBCuLu7i0WLFsli2rRpI2bOnKlTPhWNU0U5V2X7D/Z/9+5dYWNjI0aMGCG1ZWZmCgAiISFBa873fx46duwoRo8eLYQQYuvWreLBH+nVlVNl731hYaGwtrYWBw8elG17zJgxYtiwYVIf2j632kRGRop69eqJ7OxsqW3p0qXC29u7SusTPYhHfogeAzt27ICtrS0sLS3Rp08fDBkyRPbVha+vb6V9tG7dWvbcxcUF2dnZAIC0tDQUFRWhZ8+eOuXVsWNH2STbgIAAnDt3DiUlJQCAc+fOYdiwYWjSpAns7Ozg4eEBAMjIyKjyNkQlR3YUCkWV+6pKPhWNkyHc37+ZmRkaNGiAVq1aSW2lR9uqss158+Zh7dq1SEtLe+Q5VfTenz9/Hrdu3cLzzz8PW1tb6bFu3TrZ14xV+dwC94789O/fH05OTlJbeno6PD09dX+xRADMazoBIqpcjx49sHLlSigUCri6usLcXL7r2tjYVNpHnTp1ZM9NTEykr3se/FrJUF588UW4u7vjyy+/hKurKzQaDVq2bFmlya2enp4wMTFBWloaBgwYUGZ5WloanJyc4ODgAAAwNTUtUyjduXNH53wqGqeKVGX75fV/f1tpQVGVbT777LMICgrC9OnTMWrUqFqREwDcvHkTALBz5040atRItszCwkL6f1U+t8C94mf69OmytuTkZDz77LNVWp/oQTzyQ/QYsLGxgaenJxo3blym8DEELy8vWFlZIT4+Xqf1HpxweujQIXh5ecHMzAzXr1/H2bNnMWPGDPTs2RPNmzfX6eykBg0a4Pnnn8f//vc/3L59W7ZMrVZj/fr1sl/4Tk5OyMzMlJ7n5+cjPT1dev6w+ZRSKBTSka37Vbb96vLpp5/ip59+QkJCwiPNqaL3vkWLFrCwsEBGRgY8PT1lD13nSOXl5eHChQto27atrD05ORk+Pj4P+zLISPHIDxHB0tIS77//PqZOnQqFQoHOnTvj6tWrOH36NMaMGVPuehkZGQgPD8frr7+OpKQkLFu2DAsXLgQA1KtXDw0aNMCqVavg4uKCjIwMnSeoLl++HJ06dUJQUBA+/vhjqFQqnD59Gu+99x6eeeYZREZGSrHPPfccYmJi8OKLL8LBwQGRkZEwMzOTlhsiH+DeGVSJiYm4cOECbG1tUb9+fZiamla6/erSqlUrDB8+HEuXLi2zrDpzqui9r1u3LqZMmYJ33nkHGo0GXbp0QV5eHg4cOAA7OzuEhIRUeTsnT56Eubm57Gu4ixcv4saNGyx+SG8sfogIABAREQFzc3NERkbiypUrcHFxwRtvvFHhOiNHjsTt27fh5+cHMzMzTJ48GePHjwdw7yuXjRs34u2330bLli3RtGlTLF26FN27d69yTl5eXjhy5AhmzZqFV155BdnZ2RBCYODAgfjmm29gbW0txU6fPh3p6el44YUXYG9vjzlz5siOchgiHwCYMmUKQkJC0KJFC9y+fRvp6enw8PCodPvVafbs2YiNjS3TXp05VfTeA8CcOXPg5OSEqKgo/PPPP3BwcEC7du3wwQcf6LSdEydOoGnTprC0tJTajh8/DgcHB2nOFpGuTERlswmJiLTo3r07fHx8ylxHprrNnDkTn3/+OeLi4tCxY8dHum0iejLwyA8RPVY++ugjeHh44NChQ/Dz84OpKacuEpFuWPwQ0WMnNDS0plMgoscYv/YiIiIio8LjxURERGRUWPwQERGRUWHxQ0REREaFxQ8REREZFRY/REREZFRY/BAREZFRYfFDRERERoXFDxERERkVFj9ERERkVFj8EBERkVH5fzCkrdrKniBFAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Calculate the transition rates\n", + "state_sp, transition_rates_sp = state.get_spontaneous_transition_rates(unit=\"1/mus\")\n", + "print(f\"Number of possible spontaneous decay transitions: {len(transition_rates_sp)}\")\n", + "\n", + "state_bbr, transition_rates_bbr = state.get_black_body_transition_rates(temperature, \"kelvin\", unit=\"1/mus\")\n", + "print(f\"Number of considered BBR transitions: {len(transition_rates_bbr)}\")\n", + "\n", + "# Plot the transition rates\n", + "fig, ax = plt.subplots(figsize=(6, 5))\n", + "\n", + "n_list = np.arange(0, np.max([s.n for s in state_bbr]) + 1)\n", + "rates_summed = {}\n", + "for key, state, rates in [\n", + " (\"BBR\", state_bbr, transition_rates_bbr),\n", + " (\"SP\", state_sp, transition_rates_sp),\n", + "]:\n", + " rates_summed[key] = np.zeros(len(n_list))\n", + " for i, s in enumerate(state):\n", + " rates_summed[key][s.n] += rates[i]\n", + "\n", + "ax.bar(n_list, rates_summed[\"SP\"], label=\"Spontaneous Decay\", color=\"blue\", alpha=0.8)\n", + "ax.bar(n_list, rates_summed[\"BBR\"], label=\"Black Body Radiation (BBR)\", color=\"red\", alpha=0.8)\n", + "ax.legend()\n", + "\n", + "ax.set_xlabel(\"Principal Quantum Number $n$\")\n", + "ax.set_ylabel(r\"Transition Rates (1 / $\\mu$s)\")\n", + "ax.set_title(\"Spontaneous and BBR Transition Rates vs $n$\")\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Lifetime scaling with the principal quantum number" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As a more sophisticated example, we study how the lifetime scales with the effective principal quantum number $\\nu$. Our numerics reproduce the $\\nu^3$ scaling which one expects for states with a small angular quantum number $l$. For circular states, the lifetime scales as $\\nu^5$." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAHECAYAAAATY9HhAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAdPJJREFUeJzt3Xd8VFX+//HXzKR3kkAqEDqG3kIRiCgKqNjLuggoigKxIt9VdlfQXV11/elaNoCigl3ctaGuWFATEAQBoyJdIpCQQhLS+8z9/REYGAiQhCQzSd7PxyMPmHvuvfOZuQPzybnnfI7JMAwDERERERdidnYAIiIiIidSgiIiIiIuRwmKiIiIuBwlKCIiIuJylKCIiIiIy1GCIiIiIi5HCYqIiIi4HCUoIiIi4nLcnB2Aq7HZbBw8eBB/f39MJpOzwxEREWkxDMOgqKiIyMhIzOaz6wNRgnKCgwcP0rFjR2eHISIi0mIdOHCA6OjoszqHEpQT+Pv7AzVvbkBAgJOjERERaTkKCwvp2LGj/bv0bChBOcHR2zoBAQFKUERERBqgMYZIaJDsEYmJicTGxjJs2DBnhyIiItLmmbSasaPCwkICAwMpKChQD4qIiEg9NOZ3qHpQRERExOVoDEoDWK1WqqqqnB2GtAIWiwU3NzdNaRcROYESlHoqLi4mLS0N3RmTxuLj40NERAQeHh7ODkVExGUoQakHq9VKWloaPj4+tG/fXr/1ylkxDIPKykoOHTpEamoqPXr0OOvCRiIirYUSlHqoqqrCMAzat2+Pt7e3s8ORVsDb2xt3d3f27dtHZWUlXl5ezg5JRMQl6Ne1I+ozzVg9J9KY1GsiInIy/c94REJCAtu2beOHH35wdigiIiJtnm7xiIiItCU2K+xbB8VZ4BcGnUeB2eLsqE6iBMUJrDaDjal5ZBeV08Hfi7guwVjMum0kIiJNbNtKWHU/FB48ti0gEiY+AbGXOS+uWugWTzNbtTWD0U98zQ1Lv+fud1K4Yen3jH7ia1ZtzWiy5zx06BCzZ8+mU6dOeHp6Eh4ezoQJE/juu+8adL6HHnqIgQMH1vu45cuXExQU1KDnFBGRs7RtJbw7zTE5ASjMqNm+baVz4joFJSjNaNXWDGa/sYWMgnKH7ZkF5cx+Y0uTJSlXX301P/74I6+++iq7du1i5cqVnHfeeeTm5jbJ84mIiIuxWWt6TjhWw8tWDTUlvY5sW/VAzX4uQglKM7HaDB7+eBu1lXc7uu3hj7dhtTVuAbj8/HzWrFnDE088wbhx4+jcuTNxcXHMnz+fyy47dXfet99+S1xcHL6+vgQFBXHuueeyb98+li9fzsMPP8xPP/2EyWTCZDKxfPlyAJ5++mn69euHr68vHTt2ZM6cORQXF9vPd/PNN1NQUGA/7qGHHgKgoqKCefPmERUVha+vL8OHD+fbb7+1x7Jv3z4mT55Mu3bt8PX1pU+fPvzvf/9r1PdJRKRV27fO3nNiqzaRu8OXPZ+EUXzQ88gOBhSm1+znIjQGpZlsTM07qefkeAaQUVDOxtQ8RnYLabTn9fPzw8/Pjw8//JARI0bg6el5xmOqq6u54oormDlzJm+//TaVlZVs3LgRk8nE9ddfz9atW1m1ahVfffUVAIGBgUDNdNnnnnuOLl26sHfvXubMmcOf/vQnFi1axKhRo3jmmWdYsGABO3futMcGcMcdd7Bt2zbeeecdIiMj+eCDD5g4cSK//PILPXr0ICEhgcrKSpKTk/H19WXbtm32Y0VEpA6Ks7BVmcjb40veDl+sFTWDYvP3+uAfVeGwn6tQgnJEYmIiiYmJWK1N072VXXTq5KQh+9WVm5sby5cvZ+bMmSxZsoTBgwcTHx/PH/7wB/r371/rMYWFhRQUFHDppZfSrVs3AM455xx7u5+fH25uboSHhzscd88999j/HhMTwyOPPMKsWbNYtGgRHh4eBAYGYjKZHI7bv38/y5YtY//+/URGRgIwb948Vq1axbJly/jHP/7B/v37ufrqq+nXrx8AXbt2bZT3RkSkLbAWF3P4k83kfdwBa2VNYuLuV01obBGBMWWOO/uFOSHC2ukWzxFNXQelg3/dKoTWdb/6uPrqqzl48CArV65k4sSJfPvttwwePNh+a+ZEwcHB3HTTTUyYMIHJkyfz7LPPkpFx5vExX331FRdccAFRUVH4+/szdepUcnNzKS0tPeUxv/zyC1arlZ49e9p7e/z8/EhKSuK3334D4K677uKRRx7h3HPPZeHChfz8888Neh9ERNoSa2EhhxIT2XP+BRx69UOslRY8/KuJGH6YbhdnE9S1DJM9CzBBQFTNlGMXoQSlmcR1CSYi0ItTTSY2ARGBNVOOm4KXlxcXXnghDz74IOvWreOmm25i4cKFp9x/2bJlrF+/nlGjRrFixQp69uzJ999/f8r9f//9dy699FL69+/Pe++9x+bNm0lMTASgsrLylMcVFxdjsVjYvHkzKSkp9p/t27fz7LPPAnDrrbeyd+9epk6dyi+//MLQoUN5/vnnG/hOiIi0btb8fA499zx7LhhPzvP/xlZYiEfXrkTe8we6TjpEUJfy4xIT4Og308THXaoeihKUZmIxm1g4ORbgpCTl6OOFk2ObrR5KbGwsJSUlp91n0KBBzJ8/n3Xr1tG3b1/eeustADw8PE66FbZ582ZsNhtPPfUUI0aMoGfPnhw86DiVrbbjBg0ahNVqJTs7m+7duzv8HH8rqGPHjsyaNYv333+f++67j6VLl57NyxcRaXWqDx8m+1/P1CQmixZhKyrCs0d3ov71NF0/XkngrIWY/vAaBEQ4HhgQCde95nJ1UDQGpRlN7BvB4hsH8/DH2xwGzIYHerFwciwT+0ac5uiGyc3N5dprr2XGjBn0798ff39/Nm3axD//+U8uv/zyWo9JTU3lxRdf5LLLLiMyMpKdO3eye/dupk2bBtSML0lNTSUlJYXo6Gj8/f3p3r07VVVVPP/880yePJnvvvuOJUuWOJw3JiaG4uJiVq9ezYABA/Dx8aFnz55MmTKFadOm8dRTTzFo0CAOHTrE6tWr6d+/P5dccgn33HMPkyZNomfPnhw+fJhvvvnGYUyMiEhbVp2bS96yZeS99TbGkVvqnr16ETpnDv4Xjsd0/HpfsZdB70taRCVZDHFQUFBgAEZBQcFJbWVlZca2bduMsrKys3qOaqvNWLcnx/jwxzRj3Z4co9pqO6vznU55ebnxwAMPGIMHDzYCAwMNHx8fo1evXsZf//pXo7S0tNZjMjMzjSuuuMKIiIgwPDw8jM6dOxsLFiwwrFar/ZxXX321ERQUZADGsmXLDMMwjKefftqIiIgwvL29jQkTJhivvfaaARiHDx+2n3vWrFlGSEiIARgLFy40DMMwKisrjQULFhgxMTGGu7u7ERERYVx55ZXGzz//bBiGYdxxxx1Gt27dDE9PT6N9+/bG1KlTjZycnCZ7z5pbY32uRKSVsVYbxt5kw/j5PzV/Wqsdmquys43Mx58wtg8cZGzr1dvY1qu3sffKq4zCr74ybEf+v25up/sOrS+TYRiNW3ijhSssLCQwMJCCggICAgIc2srLy0lNTaVLly54eTX+YFZpm/S5EpGTnKYkfVXICHJffon8Fe9iVNRMEfbq14/QObPxO+88TCbnLZ1yuu/Q+tItHhEREVdytCT9CaU9qzKzyL3/DvJ/D8SoqgbAe8AAQu9IwHf0aKcmJk1BCYqIiIirqKUkfVWJhZztfhTs9cGwmYBqvAcPpv0dCfiMHNnqEpOjlKCIiIi4iuNK0lcWW8jd5kd+qg8YNUmIT4cKQvsU4fN/t2Pq6jo1S5qCEpQjmrqSrIiIyBkVZ1FZZCFnmz8Fv3vbExPfsCOJSYcjdaVKsp0YZPNQgnJEQkICCQkJ9gE+IiIizalibyq5Sz6l4OsOxxKT8HJC+xTj0/6EgpcuVJK+qShBERERcaKKPXvIWfIChf/7H9hsgAm/yHJC+xThHVJ1wt6mmtk8LlSSvqkoQREREXGC8p27yFmymKJVn8ORih9+559P6KQ+eG/+cy1HuGZJ+qaiBEVERKQZlW/fTs6ixRR9+aV9m/+FFxI6exZesTVLotAt+hR1UB53uZL0TUUJitiZTCY++OADrrjiiiZ7jm+//ZZx48Zx+PBhgoKCmux5RERcTdnWX8lZvJji1atrNphM+E+cQOis2Xj16um4c0sqSd9EtFigM9iskLoGfvlvzZ+2pp85lJmZyZ133knXrl3x9PSkY8eOTJ48mdVH/6EAGRkZTJo0qcljcYaHHnqIgQMH1vu45cuXK5ESkbNS9vPPHLh9Fr9fc01NcmIyEXDJJXRd+RHR//rXycnJUWYLdBkD/a6p+bMNJSegHpTmd5ryxU3Vbff7779z7rnnEhQUxJNPPkm/fv2oqqri888/JyEhgR07dgA4rB5cm6qqKtzd3ZskxvqorKzEw8PD2WGIiNT8gnmKXo7SH38kZ9FiStasqdnXbCZw8qWE3H47nl27OjHolkE9KM3paPni45MTgMKMmu3bVjbJ086ZMweTycTGjRu5+uqr6dmzJ3369GHu3Ll8//339v1MJhMffvghUJPUmEwmVqxYQXx8PF5eXrz55psAvPLKK/Tp0wdPT08iIiK44447HI5JSUmxnzM/Px+TycS3335ba2y5ubnccMMNREVF4ePjQ79+/Xj77bcd9jnvvPO44447uOeeewgNDWXChAm1nuvbb78lLi4OX19fgoKCOPfcc9m3bx/Lly/n4Ycf5qeffsJkMmEymVi+fDkATz/9NP369cPX15eOHTsyZ84ciouL7ee7+eabKSgosB/30EMPAVBRUcG8efOIiorC19eX4cOHO7zGffv2MXnyZNq1a4evry99+vThf//7X10ul4i0FNtWwjN94dVL4b1bav58pi+l/32W/TNmsO+GP9YkJxYLgVdeSbf/fUrkE08oOakj9aA0l1rKFx9jACZY9UDNPcdG7MbLy8tj1apVPProo/j6+p7UfqbbFw888ABPPfUUgwYNwsvLi8WLFzN37lwef/xxJk2aREFBAd99912D4ysvL2fIkCHcf//9BAQE8OmnnzJ16lS6detGXFycfb9XX32V2bNnn/K5qqurueKKK5g5cyZvv/02lZWVbNy4EZPJxPXXX8/WrVtZtWoVX331FYC91o3ZbOa5556jS5cu7N27lzlz5vCnP/2JRYsWMWrUKJ555hkWLFjAzp07AfDz8wPgjjvuYNu2bbzzzjtERkbywQcfMHHiRH755Rd69OhBQkIClZWVJCcn4+vry7Zt2+zHikgrcMJ6OYYBpdke5HxdQWn2kpp93NwIvOJyQm+/HY+OHZ0XawulBKW5HFe+uHYGFKbX7NdlTKM97Z49ezAMg969ezfo+HvuuYerrrrK/viRRx7hvvvu4+6777ZvGzZsWIPji4qKYt68efbHd955J59//jnvvvuuQ4LSo0cP/vnPf57yPIWFhRQUFHDppZfSrVs3AM455xx7u5+fH25ubifdxrrnnnvsf4+JieGRRx5h1qxZLFq0CA8PDwIDAzGZTA7H7d+/n2XLlrF//34iIyMBmDdvHqtWrWLZsmX84x//YP/+/Vx99dX069cPgK76jUmk9TjuF07DgNIsDw796k/ZIc+adrNBUC8TIc99ikfHTk4NtSVTgtJcirMad786MozaemzqbujQofa/Z2dnc/DgQS644IKzDcvOarXyj3/8g3fffZf09HQqKyupqKjAx8fHYb8hQ4ac9jzBwcHcdNNNTJgwgQsvvJDx48dz3XXXERERcdrjvvrqKx577DF27NhBYWEh1dXVlJeXU1paelIMR/3yyy9YrVZ69nQc2FZRUUFISAgAd911F7Nnz+aLL75g/PjxXH311fTv3/9Mb4eItAT71mEUHKQk05OcX/0py6kZE2cyGwR1KyWkdxHuvjao3gcoQWkojUFpLnUtS9zI5Yt79OiByWSyD4Str+NvC3l7e592X7O55uN0fFJUVXViFURHTz75JM8++yz3338/33zzDSkpKUyYMIHKSseyzrXdnjrRsmXLWL9+PaNGjWLFihX07NnTYYzNiX7//XcuvfRS+vfvz3vvvcfmzZtJTEwEOOn5j1dcXIzFYmHz5s2kpKTYf7Zv386zzz4LwK233srevXuZOnUqv/zyC0OHDuX5558/42sQEddmGAZFScn8/mUoB5JCKMvxwGQxaNejmG6XZhE+pKAmOYFG/4WzrVGC0lw6j6qZrcOplsU2QUBUo5cvDg4OZsKECSQmJlJSUnJSe35+fp3P5e/vT0xMjMPU5OO1b98eqJmufNTxA2Zr891333H55Zdz4403MmDAALp27cquXbvqHNOJBg0axPz581m3bh19+/blrbfeAsDDw+OkhSA3b96MzWbjqaeeYsSIEfTs2ZODBx1vw9V23KBBg7BarWRnZ9O9e3eHn+NvBXXs2JFZs2bx/vvvc99997F06dIGvy4RcS7DMChavZrfr76GtMffoDzPA5PFRnCvo4lJIe4+NseD2sB6OU1JCUpzMVtqphIDJycpTVu++OgqzXFxcbz33nvs3r2b7du389xzzzFy5Mh6neuhhx7iqaee4rnnnmP37t1s2bLF3jPg7e3NiBEjePzxx9m+fTtJSUn89a9/Pe35evTowZdffsm6devYvn07t99+O1lZ9f+tIzU1lfnz57N+/Xr27dvHF198we7du+3jUGJiYkhNTSUlJYWcnBwqKiro3r07VVVVPP/88+zdu5fXX3+dJUuWOJw3JiaG4uJiVq9eTU5ODqWlpfTs2ZMpU6Ywbdo03n//fVJTU9m4cSOPPfYYn376KVAztuXzzz8nNTWVLVu28M033ziMiRGRlsGw2Sj84gtSr7yKtIQ7KN+2DZO3NyEDoPvkQ4QNKsTd+4TEpIl+4WxrlKAckZiYSGxs7FkN+Dyj2Mvgutcg4IRxEQGRNdubqA5K165d2bJlC+PGjeO+++6jb9++XHjhhaxevZrFixfX61zTp0/nmWeeYdGiRfTp04dLL72U3bt329tfeeUVqqurGTJkCPfccw+PPPLIac/317/+lcGDBzNhwgTOO+88wsPDG1TJ1sfHhx07dtinUd92220kJCRw++23A3D11VczceJExo0bR/v27Xn77bcZMGAATz/9NE888QR9+/blzTff5LHHHnM476hRo5g1axbXX3897du3tw/UXbZsGdOmTeO+++6jV69eXHHFFfzwww906lRzv9lqtZKQkMA555zDxIkT6dmzJ4sWLar36xIR5zBsNgpXrSL1iitJv+tuKnbswOzjQ8htt9H969V0WPgEbl41C/s5alvr5TQlk3G2oyhbmcLCQgIDAykoKCAgIMChrby8nNTUVLp06YKXl1fDn+Q0hX2k7Wm0z5WInDXDaqXws1XkLFlM5Z7fADD7+RE8bSrB06ZhOb40Q62FN6Pa1Ho5Jzrdd2h9aRaPMxwtXywiIi7BqK6m8NNPyVm8hMrffwfAHBBA8LRpBE+9EcuR2kkOtF5Ok1KCIiIibZZRVUXBx5+Q88ISqvbtB8ASGEjwzTfRbsoULP7+pz+BfuFsMkpQRESk9TnDrXSjspL8jz4i94UXqUpLA8DSrh3BM26m3Q1/xOJ35tIG0rSUoIiISOtymkVZbd0nUvD+++S8+CLVB2tKIlhCQgiZMYN2f7gecx1qLknzUIIiIiKtxwlr5BxlO5xB/qOzyP29E9W5BQBY2ocSeuutBF13HeYzFKKU5qcEpQE08Ukakz5PIo2klkVZbdWQ/5svudv9qC63AAW4dehAyMyZBF17DWbNnHNZSlDqwWKpuX9ZWVl5xrLvInVVWloKgLu7u5MjEWnhjluU1VZt4vAeH3J3+GEtr/m/282nmtBziglckIi51/nOjFTqQAlKPbi5ueHj48OhQ4dwd3e3rz0j0hCGYVBaWkp2djZBQUH2BFhEGqg4C1uVicN7fMnd4Yu1oubflLtPNSGxxQR1KcVkASrznBun1IkSlHowmUxERESQmprKvn37nB2OtBJBQUEOa/iISP1Zi4s5/Mlm8j7ugLXySGLiW01obDGBXUoxHf/7pNbIaRGUoNSTh4cHPXr0OO1qtyJ15e7urp4TkbNgLSwk7403yHv1NWwFBYAFD/9qQmKLCOxc5piYYKqZzaM1cloEJSgNYDabVZJcRMSJrPn55L32Onmvv46tqAgAj65dCZ0cR0DGv05ITEBr5LQ8SlBERKTFqD58mLxXX+Xw629gKykBwLNHd0Jnz8Z/wgRMFgtsG3SKOihtd42clkgJioiIuLzqvDzyli0j7823MI7MfPPs1YvQOXPwv3A8puMnLWiNnFZBCYqIiLis6pwccl9+hcPvvINRVgaAZ+w5tJ8zB7/zz3dMTI6nNXJaPCUoIiLicqqyssl75WUOv7MCo6ICAK++fQlNmIPfeedhMpmcHKE0NSUoIiLiMqoyM8ld+hL5//kPxpHZkl4D+tM+IQHfMWOUmLQhSlBERKR5nGaF4ar0dHKWLqXgvfcxqqoA8B48mNCEOfiOGqXEpA1qlQlKTEwMAQEBmM1m2rVrxzfffOPskERE2rZTrDBcOeh+cr/+jfwPPoTqagB8hg0jNGEOPsOHKzFpw1plggKwbt06/Pz8nB2GiIjUssJwZZGFnA2lFCx9HIyaJMRn5AhCZ8/GNy7OSYGKK2m1CYqIiLiAE1YYrii0kLvNn4J93vbExDcaQh9/DZ+hw5wYqLgal1vtLjk5mcmTJxMZGYnJZOLDDz88aZ/ExERiYmLw8vJi+PDhbNy40aHdZDIRHx/PsGHDePPNN5spchEROcmRFYYrCtxIXx/E3s86UPC7DxgmfCPKiRl/iE6jD+ITUu7sSMXFuFwPSklJCQMGDGDGjBlcddVVJ7WvWLGCuXPnsmTJEoYPH84zzzzDhAkT2LlzJx06dABg7dq1REVFkZGRwfjx4+nXrx/9+/dv7pciItLmlf/6EznftaPogBdHy837RZYT2rcI7+CqYzsWZzknQHFZLpegTJo0iUmTJp2y/emnn2bmzJncfPPNACxZsoRPP/2UV155hQceeACAqKgoACIiIrj44ovZsmXLKROUiooKKo7MsQcoLCxsrJciItJmlW/fTs7iJRR98QXgDYB/dBmhfYrwald98gFaYVhO4HK3eE6nsrKSzZs3M378ePs2s9nM+PHjWb9+PVDTA1N0ZOGo4uJivv76a/r06XPKcz722GMEBgbafzp27Ni0L0JEpBUr2/orBxLuIPXKq2qSE5MJ/67QZeIhokcfriU5MUFAlFYYlpO4XA/K6eTk5GC1WgkLc8y0w8LC2LFjBwBZWVlceeWVAFitVmbOnMmwYaceeDV//nzmzp1rf1xYWKgkRUSknsp+/pmcxEUUJyXVbDCZCJg0idDZs/Cs2n5kFo+J42fyaIVhOZ0WlaDURdeuXfnpp5/qvL+npyeenp5NGJGISOtV+uOP5CxaTMmaNTUbzGYCLr2E0Fmz8Oza9chePeC617TCsNRLi0pQQkNDsVgsZGU5DqbKysoiPDzcSVGJiLQ9pZs2kbNoESXram6vY7EQeNllhN5+Gx4xMScfoBWGpZ5aVILi4eHBkCFDWL16NVdccQUANpuN1atXc8cdd5zVuRMTE0lMTMRqtTZCpCIirY9hGJRu/IGcxERKj5Z3cHMj8IrLCb3tNjw6dTr9CbTCsNSDyyUoxcXF7Nmzx/44NTWVlJQUgoOD6dSpE3PnzmX69OkMHTqUuLg4nnnmGUpKSuyzehoqISGBhIQECgsLCQwMPNuXISLSahiGQen69RxatIiyTZtrNrq7E3TllYTcdhse0VHODVBaJZdLUDZt2sS4cePsj48OYJ0+fTrLly/n+uuv59ChQyxYsIDMzEwGDhzIqlWrTho4KyIiZ8cwDErWriUncRFlKSkAmNzdCbr2GkJuvRX3yEjnBiitmskwDOPMu7UdR3tQCgoKCAgIcHY4IiJN4zQrCxuGQXFSEjmLFlP+888AmDw9CbruOkJuvQV3/UIop9CY36Eu14PiLBqDIiJtxilWFjYmPE5xph85iYso37YNAJOXF+2uv57gW2bgfqRat0hzUA/KCdSDIiKtWi0rCxsGFKV5k/OrHxX57gCYvL1p98cbCLn5ZtxCQ50UrLQ06kEREZH6O2FlYcMGRWle5PzqT0VBTWJidod2N91K8M034xYc7MRgpa1TgiIi0lYcWVnYsEHhgZoek8rCo4mJjXY9SgjuVYzbVXGg5EScTAmKiEgbYRQcpDDVm5xt/lQW1fz3b3a3EdyrmOCeJVg8jtz20crC4gKUoByhQbIi0loZVVUUrPyYnH8nUpXRDgCzh42QXsW063FcYnKUVhYWF6BBsifQIFkRaS2MykryP/qI3BdepCotDQCLJwT3KqxJTNxP/O/fVLM+zj2/qAS9NIgGyYqIyCnZKispeP99cl58keqDGQBYgoMJuWUG7QYGYP54Zi1HaWVhcS1KUEREWglbRQX5//0vuUtfojozEwBL+1BCbrmFdtdfj9nbu2ZHb0+tLCwuTwmKiEgLZysvJ//dd2sSk0OHAHDr0IGQmTMJuvYazF5ejgdoZWFpAZSgHKFBsiLS0thKSzm84l1yX34Za04OAG7h4YTcNpOgq6/G7Ol56oO1srC4OA2SPYEGyYqIq7OVlHD47bfJfWUZ1rw8ANwjIwm5/XYCr7wCs4eHkyOUtkqDZEVE2iBrcTGH33yLvGXLsObnA+DesSOhs24n8LLLMLm7OzdAkUakBEVExMVZCwvJe+MN8l59DVtBAQDunTsROms2gZdeosREWiUlKCIiLspaUEDea6+T99pr2IqKAPDo2pXQ2bMImDQJk5v+C5fWS59uERFnsFlPOYum+vBh8l59lcOvv4GtpAQAzx7dCZ09G/8JEzBZNNtGWj8lKCIizW3bylrrkFSPfJC8dRkcfvMtbKWlAHj26lWTmFx0ISaz2UkBizQ/JShHaJqxiDSLbSvh3WnAsQmU1eVmclOKOfzywxjWmiTEM/Yc2s+Zg9/55ysxkTZJ04xPoGnGItJkbFZ4pq+956SqzEzedj8O/+ZjT0y82kPoQ/+uSUxMJmdGK1JvmmYsItIS7VsHhQepKjWTu92P/N98MWw1SYhXSCXt+xThG1GBqasHKDmRNk4JiohIM6lK3UnOpkAK9vrYExPv0ApC+xbjG1ZxLCcpznJekCIuQgmKiEgTq0xLJ/eFF8j/4H2o9gXAp30FoX2L8OlQeXJniV9Y8wcp4mKUoIiINJHK/fvJeeEFCj5aCdXVAPhEQWivXHw7VNRyhKlmVeHOo5o3UBEXpARFRKSRVaSmkrvkBQo++QSOzAz0PfdcQhPm4OOVdmQWj4njZ/LUPAYmPq5VhUVQgiIi0mgqfvuNnCUvUPjpp2CzAeAbP5b2s2fjPXDgkb0Gw3Wv1VoHhYmPQ+xlzR63iCtSgnKE6qCISEOV79pF7pIlFH62Co5UbvAbN47QOXPw7tf35ANiL4Pel5yykqyIqA7KSVQHRUTqqnzHDnIWL6Ho88/t2/wvHE/o7Nl4xcY6MTIR51AdFBERJyr79VdyFi+m+KvVNRtMJvwnTCB09iy8evVybnAirYQSFBGROir75RdyEhdR/O23NRtMJgImTSJ09iw8e/RwamwirY0SFBGRMyhLSeHQokWUJK+p2WA2E3DpJYTOmoVn167ODU6klVKCIiJyCqWbN5OTuIiSdetqNlgsBE6eTOis2/GIiXFqbCKtnRIUEZETlGzcSE7iIko3bKjZ4OZG4BWXE3rbbXh06uTc4ETaCCUoItJ22KynnNprGAal339fk5hs2lSzv7s7QVdeSchtt+ERHeXEwEXaHiUoItI2bFtZa3E0Y8LjlBwOJWfRIsq2bAHA5O5O0LXXEHLrrbhHRjopYJG2TQmKiLR+21YeKS9/rOyTYUDJjjwOvTeX8lwPAEyengRddx0ht96Ce5gW7BNxJiUoR6iSrEgrZbPW9JwcSU4MA4oPepKz1Z/yw0cSEwu0u3EawbfcgnuHDk4MVkSOUiXZE6iSrEgrk7oGXr0Uw4CidC9yfvWjwp6Y2GjXo5SQ3sW43b4SuoxxcrAiLZsqyYqI1JFRmEHRfi9yfvWnosAdALObjXY9SgjuVYKbV82ifhRnOTFKETmREhQRaZUMq5XCVavIeXYJlfuDATC7H01MinHzPKHz2E9jTkRciRIUEWlVjOpqCv/3P3IWL6EyNRUAswcE9ywiuGcxFo8T72qbICCyZsqxiLgMJSgi0ioYVVUUfPwJOS8soWrffgDMgYGE3DSddkODsXx6ey1HmWr+mPi4vR6KiLgGJSgi0qIZlZUUrFxJzgsvUnXgAACWoCCCZ8yg3R9vwOLnV7Ojr3etdVCY+DjEXuaEyEXkdJSgiEiLZKuspOD9D8h98UWqDtYkHZbgYEJumUG7P/wBs6+v4wGxl0HvS05ZSVZEXIsSFBFpUWwVFeT/97/kLn2J6sxMACztQwm55RbaXX89Zm/vUx9stmgqsUgLoQRFRFoEW3k5+e/+h9yXXqI6OxsAtw4dCLn1VoKuuxazl5eTIxSRxqQERURcmq20lMMr3iX35Zex5uQA4BYeTshtMwm6+mrMnp5OjlBEmoISFBFxSbaSEg6//Ta5ryzDmpcHgHtkJCG3307glVdg9vBwcoQi0pSUoIiIS7EWF3P4zbfIW7YMa34+AO7R0YTOup3Ayy7DpMREpE1QgiIiLsFaWEjeG2+Q9+pr2AoKAHDv3InQWbMJvPQSTO7uTo5QRJqTEhQRcSprQQF5r71O3muvYSsqAsCjSxdCZ88i4OKLMbnpvymRtkj/8o9ITEwkMTERq9Xq7FBE2oTqw4fJe/VVDr/+BraSEgA8uncjdPZsAiZOxGRRfRKRtsxkGMaJC1O0aY25VLRIm2SznrYYWnVeHnnLlnP4zTexlZYC4NmzJ6Fz5uB/0YWYzGZnRS4iZ6kxv0PVgyIijWfbylOUk3+C6g6jyH1lGYfffhujrAwAz3POIXTObPwvuECJiYg4UIIiIo1j20p4dxrg2ClblZVF3gMJHP69HUZlFQBeffsSOmcOfuPOw2QyNX+sIuLylKCIyNmzWWt6To5LTqpKzeRu9yP/N18Mmwmowqt/f9onzMF37FglJiJyWkpQROTs7Vtnv61TVWIhZ7sfBXt9jiQm4B1aSWifInzvT8DUdawzIxWRFkIJioicveIsKost5G7zI/93HziSmPi0ryC0bxE+HSoxmYCSbOfGKSIthhIUETkrlfv3k7N0FQVfdgDjSGLSoaKmxySs0nFnvzAnRCgiLZESFBFpkIrUVHKXvEDBJ5+A1QqY8A070mPS/oTEBFPNbJ7Oo5wRqoi0QEpQRKReKn77jZwlL1D46adgswHgO3YMoZP64fPTg7UccWQw7MTHHeqhiIicjhIUEamT8l27yF2yhMLPVsGR+o5+48YROmc23v361ezUq/Mp6qA8DrGXOSFqEWmplKCIyGmV79hBzqLFFH3xhX2b3/gLCJ09G+8+fRx3jr0Mel9y2kqyIiJ1oQRFRGpV9uuv5CxeTPFXq+3b/CdMIHT2LLx69z71gWYLdBnTDBGKSGumBEVEHJT9/DM5ixZT/O23NRtMJgImTSJk1u149ezp1NhEpO1QgiIiAJT++CM5ixZTsmZNzQazmYBLLiF01u14duvm3OBEpM1RgiLSxpVu3kxO4iJK1q2r2WCxEDh5MiG334Znly7ODU5E2iwlKCJtVMmGjeQsWkTphg01G9zcCLz8MkJvuw2Pzp2dG5yItHlKUETaEMMwKP3+e3ISF1G6aVPNRnd3gq68kpDbZuIRHe3cAEVEjlCCItIGGIZBydrvyFm0iLIffwTA5O5O4DVXEzpzJu6RkU6OUETEUatNUEpLSznnnHO49tpr+X//7/85OxwRpzAMg5LkZA4tWkT5Tz8DYPLwIOi66wi59Rbcw8OdHKGISO1abYLy6KOPMmLECGeHIeIUhmFQ/M035CQuovzXXwEweXnR7vrrCb5lBu4dOjg5QhGR02uVCcru3bvZsWMHkydPZuvWrc4OR6TZGDYbRV99Rc7iJVRs3w6AydubdjfcQMiMm3ELDXVyhCIidWN2dgAnSk5OZvLkyURGRmIymfjwww9P2icxMZGYmBi8vLwYPnw4GzdudGifN28ejz32WDNFLOJ8hs1G4apVpF5xJel33U3F9u2YfXwImTmT7qu/IuxP/6fkRERaFJdLUEpKShgwYACJiYm1tq9YsYK5c+eycOFCtmzZwoABA5gwYQLZ2dkAfPTRR/Ts2ZOeqngprYnNCqlr4Jf/1vxpswJgWK0UfPIpey+7jPR77qVi1y7Mvr6EzLqdbqu/osN9c3ELDnZy8CIi9WcyjCPLkjZAVVUVmZmZlJaW0r59e4Ib+T9Ck8nEBx98wBVXXGHfNnz4cIYNG8a///1vAGw2Gx07duTOO+/kgQceYP78+bzxxhtYLBaKi4upqqrivvvuY8GCBbU+R0VFBRUVFfbHhYWFdOzYkYKCAgICAhr19Yg0yLaVJ60QbPhFUuh7PTkrN1KZmgqA2d+f4GnTCJ42FUtgoLOiFZE2rLCwkMDAwEb5Dq33GJSioiLeeOMN3nnnHTZu3EhlZSWGYWAymYiOjuaiiy7itttuY9iwYWcVWG0qKyvZvHkz8+fPt28zm82MHz+e9evXA/DYY4/Zb+8sX76crVu3njI5Obr/ww8/3OixijSKbSvh3WlAze8Rhg0KfvcmZ1s1VcUrADAHBhJy03Ta3XgjFn9/JwYrItJ46nWL5+mnnyYmJoZly5Yxfvx4PvzwQ1JSUti1axfr169n4cKFVFdXc9FFFzFx4kR2797dqMHm5ORgtVoJCwtz2B4WFkZmZmaDzjl//nwKCgrsPwcOHGiMUEXOns1a03OCgWGFw7/58NunHcjY2I6qYjcsHjbaD4PuX35O6OzZSk5EpFWpVw/KDz/8QHJyMn369Km1PS4ujhkzZrBkyRKWLVvGmjVr6NGjR6ME2hA33XTTGffx9PTE09Oz6YMRqa9967AdPkhBqg+52/yoKq3552rxtBLSu4R23UswuxuQ+zMEjHFysCIijateCcrbb79dp/08PT2ZNWtWgwI6ndDQUCwWC1lZWQ7bs7KyCFfBKWlFbBUV5P93JbmfhFFdZgHA4mUlpHcx7bqXYnY7buhYcdYpziIi0nI1eBZPWVkZpaWl9sf79u3jmWee4fPPP2+UwGrj4eHBkCFDWL16tX2bzWZj9erVjBw58qzOnZiYSGxsbJOMnRGpK1t5OXmvvc5vF15E1sufUF1mwc3bStjgArpfmkVI7xLH5ATAL6z2kzVQTEwM/fv3Z+DAgYwbN67WfT755BN69epFjx49eOmll+zbDxw4wHnnnUdsbCz9+/fnP//5T6PGJiJtR4Nn8Vx00UVcddVVzJo1i/z8fHr37o27uzs5OTk8/fTTzJ49u0EBFRcXs2fPHgAGDRrE008/zbhx4wgODqZTp06sWLGC6dOn88ILLxAXF8czzzzDu+++y44dO04am9IQjTkCWaSubKWlHF7xLrkvv4w1JwcAt/AwQrplERSZgdlS2z9TEwREwj2/gNnSaLHExMSwdetW/Pz8am2vrq4mNjaWb775hsDAQIYMGcK6desICQkhIyODrKwsBg4cSGZmJkOGDGHXrl34+vo2Wnwi4roa8zu0wT0oW7ZsYcyYmvve//3vfwkLC2Pfvn289tprPPfccw0OaNOmTQwaNIhBgwYBMHfuXAYNGmSfiXP99dfz//7f/2PBggUMHDiQlJQUVq1a1SjJiUhzs5WUkPvSS+wZfyHZTzyBNScHt8gIwh96iG5ffEHwfU8cyT1MJxx55PHExxs1OamLjRs30qdPH6KiovDz82PSpEl88cUXAERERDBw4EAAwsPDCQ0NJS8vr1njE5HWocGl7ktLS/E/Mmvgiy++4KqrrsJsNjNixAj27dvX4IDOO+88ztSpc8cdd3DHHXc0+DlEnM1aXMzhN98ib9kyrPn5ALhHRxM663YCL7sMk4dHzY6xl8F1r51UB4WAyJrkJPayRo/NZDIRHx+P2WzmnnvuYcqUKQ7tBw8eJCoqyv44KiqK9PT0k86zefNmrFYrHTt2bPQYRaT1a3CC0r17dz788EOuvPJKPv/8c+69914AsrOzW+StkcTERBITE7Farc4ORVoxa2EheW+8Qd6rr2ErKADAvXMnQm+fReDkSzG5u598UOxl0PsS2LeuZkCsXxh0HtVkPSdr164lKiqKjIwMxo8fT79+/ejfv3+9zpGXl8e0adNYunRpk8QoIq1fg2/xLFiwgHnz5hETE8Pw4cPtg1S/+OIL++2ZliQhIYFt27bxww8/ODsUaYWsBQUcev7f7LlgPDnPPY+toACPLl2I/OcTdPv0U4KuurL25OQoswW6jIF+19T8eRbJycUXX8z06dPtj7/55htCQ0PtyfnR3pGIiAguvvhitmzZ4nB8ZGSkQ49Jeno6kZGR9scVFRVcccUVPPDAA4waNarBcYpI29bgBOWaa65h//79bNq0iVWrVtm3X3DBBfzrX/9qlOBEWrrqw4fJfuYZ9px/ATmJidiKivDo3o3Ip/4fXT/5uOZ2jtvZLSq+cOFCevbsyR/+8AcKCgr4+OOPGTBgAC+//HKt+594SyY+Pp6ysjK+//57SkpKKCoqAmoGrH/99dcn1T2Ki4tj69atpKenU1xczGeffcaECRMAMAyDm266ifPPP5+pU6ee1esSkbatwf8zHjhwgI4dO55UfyQuLu6sgxJp6arz8shbtpzDb76J7ch0fM8ePQhNmIP/RRdhMjfOOp1fffUV2dnZbN68mUWLFnHFFVeQn5/Pf/7zn1MumBkVFcWaNWvsj81mM97e3mRnZ5OVlcWVV14JgNVqZebMmfap90cHpbu5ufHUU08xbtw4bDYbf/rTnwgJCQHgu+++Y8WKFfTv39++Evnrr79Ov379GuX1ikjb0eBpxmazmeDgYAYMGMDAgQPtP5WVlTz33HO8+uqrjR1rs9A0Yzkb1Tk55L6yjMNvv41RVgaA5znnEDpnNv4XXNBoiclRTz75JKNHj7bfYh0+fDh33XXXSQNbj7d06VLmzp1r7ylJSUlh6NChpKWlqeChiJwVpy4WeFRqaio//vgjKSkp/Pjjj7z77rscPFgzy6AlfrFrkKycjarsbPJefoXDK1ZglJcD4NWnD6EJc/AbNw6T6cRpwo2jd+/erFy5kpEjR7Jq1SpMJhNPPPEE559/PhEREbUeExUVRXFxMYWFhfj5+XHvvfcyZcoUJSci4lIanKB07tyZzp07c8UVV9i3rV+/nunTp/O3v/2tMWJrVgkJCSQkJNizP5G6qMrMJPell8l/912MykoAvPr3p33CHHzHjm2yxOSoyZMns3r1ajp16kRMTAzvv/8+X3/9NSNGjOCvf/0rM2fOPOmYo4Ng09LSWLZsGZmZmXz00UdNGqeISH01+BbPqXz22Wc8+OCDbNq0qTFP22x0i0fqourgQXKWLqXgv+9hVFUB4D1oEKEJCfieO6rJE5OzkZOTQ/v27Zk4cSK7du0iOTnZoa6JiEhDucQtnsrKSjyOFpM6To8ePfj111/PKigRV1WZlk7uiy+S/8EHcDQxGTqE9gkJ+IwY4dKJyVGhoaF4enqyb98+kpKSlJyIiEtqcILi5+dHbGwsgwYNYuDAgQwaNIjIyEief/55xo8f35gxijhd5f795LzwAgUfrYTqagB8hg8nNGEOvi1w5lr5kXEyIiKuqsG3eNauXctPP/3ETz/9REpKClu3brX/pzdx4kSGDh1Kv3796NevH717927UoJvC8YNkd+3apVs8AkBFaiq5L7xIwccfw5EB1L6jRhE6ZzY+Q4c6OToREdfSmLd4Gm0Mis1mY+fOnaSkpJCSkmJPXrKzs1vUzBiNQRGAit9+I2fJCxR++inYbAD4jh1D6OzZ+LTASskiIs3BJROUU8nKympRKw0rQWnbKnbvJmfxEgo/+wyO/NPwO+88QhPm4K1iYyIip+W0QbL79++nU6dOdd4/PT1dA/CkRSjfuZOcRYsp+vxz+za/8RcQOns23ieUehcRkaZXr7KWw4YN4/bbbz/tgnoFBQUsXbqUvn378t577511gCJnzWaF1DXwy39r/rQdu+VYvm0bB+64g9TLr7AnJ/4TJtDlww/o+O9/KzkREXGSevWgbNu2jUcffZQLL7wQLy8vhgwZQmRkJF5eXhw+fJht27bx66+/MnjwYP75z39y8cUXN1XcInWzbSWsuh8KDx7bFhBJWbc7yFn1K8XffFOzzWQiYNJEQmbNwusUa9iIiEjzadAYlLKyMj799FPWrl3Lvn37KCsrIzQ0lEGDBjFhwgT69u3bFLE2C41BaUW2rYR3pwHHPuJlOe4c+tWfkgyvmg1mMwGXXELorNvx7NbNOXGKiLQSLWqQbEuhacatjM0Kz/S195yUHvIg51c/SjKPJCYmg8AeJkKe+QTPrkpMREQagxKUJqQelFYidQ28eikl2R7kbPWnNNuzZrvJIDCmlNDYYjz8rTD9E+gyxrmxioi0Ei5R6l7EVRmGQen69eSsDqH00JHExGwQFFNKSGwxHn7H1eUpznJOkCIiclpKUKTVMAyDku/WkbNoEWVbtgCemMwGgV1LCT2nGHffWgoG+rWcGj0iIm2JEhRp8QzDoCQ5mUOLFlH+088AmDw8COpRSUi3bNx9aqtkbIKASOg8qnmDFRGROlGCIi2WYRgUf/MtOYsWUb51KwAmT0/a/eF6gm+5Bfec74/M4jFx/EyemsfAxMfBbGnusEVEpA7qVajtRGvWrOHGG29k5MiRpKenA/D666+zdu3aRglOpDaGzUbhl1+SevXVpM2ZQ/nWrZi8vQmeMYPuX31J2Pz5uHfoALGXwXWvQUCE4wkCImu2x17mnBcgIiJn1OAelPfee4+pU6cyZcoUfvzxRyoqKoCaSrL/+Mc/+N///tdoQYpATWJS9MWX5CxeTMXOnQCYfXxoN2UKwTffhFtw8MkHxV4GvS+BfetqBsT6hdXc1lHPiYiIS2vwNONBgwZx7733Mm3aNPz9/fnpp5/o2rUrP/74I5MmTSIzM7OxY21SqoPiugyrlcJVq8hZvJjKPb8BYPb1pd3UGwmePh23du2cHKGIiICLTDPeuXMnY8eOPWl7YGAg+fn5ZxOTUyQkJJCQkGB/c8X5jOpqCv/3P3IWL6EyNRUAs78/wdOmETxtKhZdJxGRVqvBCUp4eDh79uwhJibGYfvatWvp2rXr2cYlbZhRXU3Bx5+Qu2QJlfv2AWAODCR4+jSCb7wRi3q2RERavQYnKDNnzuTuu+/mlVdewWQycfDgQdavX8+8efN48MEHGzNGaSOMqioKPvqInBdepOrAAQAsQUEE33wz7ab8EYufn5MjFBGR5tLgBOWBBx7AZrNxwQUXUFpaytixY/H09GTevHnceeedjRmjtHK2ykoKPviQ3BdeoOpgzdo5luBgQmbcTLsbbsDs6+vkCEVEpLmd9Vo8lZWV7Nmzh+LiYmJjY/Fr4b/lai2e5mOrqCD/vffIXfoS1RkZAFhCQwm55RbaXX8dZh8fJ0coIiL14RKDZAHKy8v5+eefyc7OxmazOczcuewy1ZiQ2tnKy8l/9z/kvvQS1dnZALh16EDIrbcSdN21mL28nByhiIg4W4MTlFWrVjF16lRyc3NPajOZTFittZUXl7bMVlbG4RUryH35ZayHcgBwCw8nZOatBF1zDWZPTydHKCIirqLBCcqdd97Jddddx4IFCwgL04Jrcmq2khIOv/MOua8sw3okoXWLjCD0ttsJvOpKzB4eTo5QRERcTYMTlKysLObOnavkRE7JWlzC4bfeIm/ZMqyHDwPgHh1N6KzbCbzsMkxKTERE5BQanKBcc801fPvtt3Tr1q0x43Ga4yvJytmxFhVx+I03yFv+KtaCAgDcO3ci9PZZBE6+FJO7u5MjFBERV9fgWTylpaVce+21tG/fnn79+uF+wpfOXXfd1SgBNjfN4mk4a0EBea+9Tt5rr2ErKgLAo0sXQmfPIuDiizG5afFsEZHWzCVm8bz99tt88cUXeHl58e2332IymextJpOpxSYoUn/Vhw+T99prHH79DWzFxQB4dOtG6OzZBEyaiMmihflERKR+Gpyg/OUvf+Hhhx/mgQcewGw2N2ZM0kJU5+WRt2w5h998E1tpKQCePXoQmjAH/4suwqTPhYiINFCDE5TKykquv/56JSdtUHVODrmvLOPw229jlJUB4Nm7N6FzZuM/frwSExEROWsNTlCmT5/OihUr+POf/9yY8YgLq8rOJu/lVzi8YgVGeTkAXn36EJowB79x4xxu84mIiJyNBicoVquVf/7zn3z++ef079//pEGyTz/99FkHJ66hKjOT3JdeJv/ddzEqKwHw6t+f0Dmz8YuPV2IiIiKNrsEJyi+//MKgQYMA2Lp1q0ObvrBaEJsV9q2D4izwC4POo8BcM6i16uBBcpYupeC/72FUVQHgPXAgoQkJ+I4+V9dZRESaTIMTlG+++aYx4xBn2LYSVt0PhQePbQuIpHLw/eR+nUr+Bx/A0cRk6BDaJyTgM2KEEhMREWlyKkzRVm1bCe9OA46VwaksspCzoZSCpY+DUZOE+AwfTuicOfgOj3NSoCIi0hbVK0GZO3cuf//73/H19WXu3Lmn3VdjUFyYzVrTc3IkOakotJC7zZ+Cfd72xMQ3CkIfexWfOCUmIiLS/OqVoPz4449UHeny//HHH0+5n24BuLh96+y3dWzV8PuX7bFV1UwN9o0oJ7RPET6hVdC+wplRiohIG1avBOX4cSevvvoq0dHRJ9VBMQyDAwcONE500jSKs+x/NbtBux4lVOS7E9qnCO+Qqlr3ExERaU4NHoPSpUsXMjIy6NChg8P2vLw8unTpokX3XJmf4wrU7fsVUWunl59WqhYREedocMnPU60xWFxcjJeXV4MDkmbQeRQERAI1WcnJyYkJAqJq9hMREXGCevegHB0cazKZWLBgAT4+PvY2q9XKhg0bGDhwYKMF2FwSExNJTExsGz0/ZgtMfOLILB4Tx8/kOZq0MPFxez0UERGR5mYyTtUVcgrjxo0DICkpiZEjR+Lh4WFv8/DwICYmhnnz5tGjR4/GjbSZNOZS0S6v1jooUTXJSexlzotLRERapMb8Dq13gnLUzTffzLPPPtvqvsTbVIICp60kKyIiUh8ukaC0Vm0uQREREWkkjfkd2uBBsgBr1qzhxhtvZOTIkaSnpwPw+uuvs3bt2rMKSkRERNq2Bico7733HhMmTMDb25sff/yRioqaol4FBQX84x//aLQARUREpO1pcILyyCOPsGTJEpYuXYq7u7t9+7nnnsuWLVsaJTgRERFpmxqcoOzcuZOxY8eetD0wMJD8/PyziUlERETauAYnKOHh4ezZs+ek7WvXrqVr165nFZSIiIi0bQ1OUGbOnMndd9/Nhg0bMJlMHDx4kDfffJN58+Yxe/bsxoxRRERE2pgGr8XzwAMPYLPZuOCCCygtLWXs2LF4enoyb9487rzzzsaMUURERNqYetVB+fnnn+nbt6/DCsaVlZXs2bOH4uJiYmNj8fPza5JAm4vqoIiIiDSM0+qgDBo0iJycHAC6du1Kbm4uHh4exMbGEhcX1+KTExEREXEN9UpQgoKCSE1NBeD333/HZrM1SVAiIiLSttVrDMrVV19NfHw8ERERmEwmhg4disVS+7ote/fubZQARUREpO2pV4Ly4osvctVVV7Fnzx7uuusuZs6cib+/f1PFJiIiIm1UvWfxTJw4EYDNmzdz9913K0ERERGRRtfgacbLli1rzDhERERE7OqVoMydO5e///3v+Pr6Mnfu3NPu+/TTT59VYCIiItJ21StB+fHHH6mqqrL//VRMJtPZRXUW8vPzGT9+PNXV1VRXV3P33Xczc+ZMp8UjIiIi9VevQm11kZaWxt/+9jdefPHFxjxtnVmtVioqKvDx8aGkpIS+ffuyadMmQkJC6nS8CrWJiIg0jNMKtdVFbm4uL7/8cmOfts4sFgs+Pj4AVFRUYBgGjZyDiYiISBNr9ATlbCUnJzN58mQiIyMxmUx8+OGHJ+2TmJhITEwMXl5eDB8+nI0bNzq05+fnM2DAAKKjo/m///s/QkNDmyl6ERERaQwul6CUlJQwYMAAEhMTa21fsWIFc+fOZeHChWzZsoUBAwYwYcIEsrOz7fsEBQXx008/kZqayltvvUVWVlZzhS8iIiKNwOUSlEmTJvHII49w5ZVX1tr+9NNPM3PmTG6++WZiY2NZsmQJPj4+vPLKKyftGxYWxoABA1izZs0pn6+iooLCwkKHHxEREXGuetdBueqqq07bnp+f39BYzqiyspLNmzczf/58+zaz2cz48eNZv349AFlZWfj4+ODv709BQQHJycnMnj37lOd87LHHePjhh5ssZhEREam/eicogYGBZ2yfNm1agwM6nZycHKxWK2FhYQ7bw8LC2LFjBwD79u3jtttusw+OvfPOO+nXr98pzzl//nyHmi6FhYV07NixSeIXERGRuql3guLqFWTj4uJISUmp8/6enp54eno2XUAiIiJSby43BuV0QkNDsVgsJw16zcrKIjw83ElRiYiISGNrUQmKh4cHQ4YMYfXq1fZtNpuN1atXM3LkyLM6d2JiIrGxsQwbNuxswxQREZGz1ODFAptKcXExe/bssT9OTU0lJSWF4OBgOnXqxNy5c5k+fTpDhw4lLi6OZ555hpKSEm6++eazet6EhAQSEhLsVfBERETEeVwuQdm0aRPjxo2zPz46gHX69OksX76c66+/nkOHDrFgwQIyMzMZOHAgq1atOmngrIiIiLRcjb4WT0untXhEREQapjG/Q12uB8VZEhMTSUxMxGq1OjsUERGROrHaDDam5pFdVE4Hfy/iugRjMZucHVajUA/KCdSDIiIiLcGqrRk8/PE2MgrK7dsiAr1YODmWiX0jnBKTS69mLCIiIk1r1dYMZr+xxSE5AcgsKGf2G1tYtTXDSZE1HiUoIiIiLYjVZvDwx9swAExVWPx2YPGpmf169JbIwx9vw2pr2TdINAZFRESkBfls+3ZyTN/iHb0Di+8eTOZqqku6Uba/O1CTpGQUlLMxNY+R3UKcG+xZUIJyhAbJioiIK7LarPyc8zNJB5JITk9m9+HdeB03xMRWFYStIoya1OTYANnsovKTztWSaJDsCTRIVkREnK2gooDv0r8jOT2ZtelrKagosLeZMFNV2glrcW+qi3sfSU5Onrnz9swRzd6DomnGIiIiLURdpgIbhsFv+b+RlJZEcloyPx36CatxrEc/wCOAc6POJT46nhERo7j0mS1kFpRTWw+DCQgPrHmelkwJioiISBM53VTgcecEszFjI8lpySSnJXOw5KDDsd2DujM2eizx0fH0b98fN/Oxr+yFk2OZ/cYWTOCQpJiOa2/p9VB0i+cEusUjIiKN4ehUYIcEwq0AN78duPntwDtwL1W2Cnubh9mDuIg44qPjGRs9lki/yDOevzXXQVEPyhEaJCsiIo3l2FRgG2bvA/akxOJ1rD5JlQ3CfMLsvSRxEXF4u3nX+Tkm9o3gwthwVZJtK9SDIiIiZ6OwspDlWz5n0YaPsfjtwuxWYm8zDBO2sk5UF/emurgXb0y9klHdQ50YbeNSD4qIiIiLMAyD1IJUktOSSUpL4sfsH7EaVtyDjrRbvagu7kl18TlYS3piWH3txx4qrqj9pKIERUREpL4qrBVsytxkH+CaVpzm0B7h3Zl9aZ2pLu6NtbQzYKn1PB38vZoh2pZJCYqIiLRZ9VkNOLs0mzVpa0hKS+L7jO8pqy6zt7mb3YkLj2NM9JiaAa6+0Yx+4msyS2svltZapgI3JSUoIiLSJp1pFozNsLE1Z6u9l2R73naH4zt4d7AnJCMiRuDj7uPQ3hamAjclDZI9gQbJioi0frVNAQYwmcux+O5m9IBsUks2k1eed6wNE/1C+zEmegzx0fH0Du6NyXT6BMMVpwI3JQ2SbQKaZiwi0jY4rAaMgckj59g0YJ9UTCYbm3Nr9vVz92Nk5Ejio+MZHTWaEO/6lY5v7VOBm5J6UE6gHhQRkdYteXcGN7+zwp6UmD1yHdptFaFUF/fmgfiruHHgebhb3J0UacujHhQREZF6yCnLsQ9wXZO2Dp9Oxwa4GoYFa0mXI7VJemNU1dQlCbbEKjlxIiUoIiLS6tgMG9tzt9sX3/s191fH9mq/minAxb2pLukBNs+TzqEpwM6lBEVERFqFkqoS1h9cf6SXZA255Y63bvqE9CE+Op5zo8Zw20sZZBVUturVgFs6JSgiIuJS6lObZH/hfnsvyaasTVTbqu1tPm4+jIocxdjosYyJHkOo97GS8g9NDtEUYBenBEVERFzGmablVlmr2JK9xd5L8nvh7w7Hd/LvxNjosYyNHsuQsCF4WDxqfZ6JfSNYfOPgk54rvBVPAW5pNIvnBJrFIyLiHKeqTWK2FGHx20lcbBa/lWyhpOrY4ntuJjeGhA2x1yaJCYyp13PWp7dGzkyzeJqA6qCIiDjPibVJzF4HcfPbjpvfTsxeaZhMBj/n1+wb7BXMmKiaCq4jI0fi7+Hf4Oe1mE2M7Fa/2ibSPNSDcgL1oIiINL9vdh1g5rtvYzlam8S9yKHdWhZFdXFvHhp/DX8YMAqzyeykSOV01IMiIiIt3oHCAySn16xzsyHjB7w7VtnbDJsH1SXda6YBF/fGqK75svOli5KTNkIJioiINIsqWxUp2SkkpyWTlJZEakGqQ7utMtheLM1a2hWMk7+iVJuk7VCCIiIiTSavPI+16WtJTktmXfo6iqqO3bqxmCwMDhvM2KixjI4ay41L9pJVUKHaJAIoQRERkTOoz0wXwzDYeXinvZfkl0O/YByXcrTzbMfoqNGM7TiWUZGjCPA4Nk7hocneqk0idkpQRETklM5UlwSgtKqUDRkb7ONJskuzHc7RO7g3Y6LGEN8xnr4hfbGYLbU+l2qTyPE0i+cEmsUjIlLjVHVJTIDJPY8p40o5ZPuRHzJ+oNJWaW/3sngxImIEYzuOZUzUGMJ9w+v1vKpN0nJpFo+IiDQpx7okAFYs3vvt04AtXll8eODY/lF+UfbaJMPCh+Hl1vDBrKpNIqAERUREarExNY+M4lzcAnbh5rcDN79dmCxl9nbDMGMt68R1sROY2n8i3YK6YTKpl0MajxKUI1RJVkTaOsMw2HV4F2vS1/DBji/x67Edk+lYH4pR7UN1Sc8jU4F7gs2HwSMH0r1dlBOjltZKY1BOoDEoItKWlFWX8UPmDyQdSCI5PZnMkkyHdmt5eE1dkuLeWMs6Ao4DXN+eOUK3Y8ROY1BERKTBMoozSE5LJjk9mQ0ZG6iwVtjbPC2eDI8YzujIMTzzsRvZed6qSyJOoQRFRKSVs9qs/Jzzs72XZPfh3Q7t4b7hxEfH2we4ert5A9DOmqG6JOI0SlBERFqA+k69Lago4Lv070hOT2Zt+loKKgrsbWaTmQHtBzA2eixjo8fSI6hHrQNcVZdEnEkJioiIi6tLsTTDMPgt/zeS0pJITkvmp0M/YTWODfr39/CvqeAaPZbRkaMJ8gqq03NP7BvBhbHhqksizU6DZE+gQbIi4kpOVywNUxV3X2qixLKV5LRkDpYcdNine1B3ey/JgPYDcDPrd1JpWhokKyLSBpxcLA1MbgVH6pLswOK7h5f3VNnbPMwexEXE2ZOSKD9N/5WWSwmKiIiL2piaR0ZBKWbvA8eSEq8Mh31sVQGc1zGea2MvIi48Dh93HydFK9K4lKCIiLiYwspC1qWv481fV+Hb43vMbiX2NsMwYSvreKRYWm9sFRFcNHQQ53VUb4m0LkpQRESczDAMUgtSSU5LJiktiR+zf7QPcDW7gWH1orq4poKrtaQnhtXP4fgO/g1f90bEVSlBERFxggprBZsyN9UUTEtLJq04zaG9a2BXRkeN4T/J/hw6FIFxQgVXULE0ad2UoIiINJPs0mzWpK0hKS2J7zO+p6z62OJ77mZ34sLjGBNdsyJwR/+OAPTzVrE0aZuUoIiINBGbYWNrzlZ7L8n2vO0O7R28O9gTkhERI2od4KpiadJWKUERETmN+lZwLaosYv3B9SSlJbE2fS155Xn2NhMm+oX2Y0z0GOKj4+kd3LvWCq4nUrE0aYuUoByRmJhIYmIiVqv1zDuLSJtQ1wquvxf+bu8l2ZK1hWqj2r6/n7sfoyJH1VRwjRpNiHfDVv61mE1aNVjaFFWSPYEqyYoInKmCazX3TrZQ5vYryWnJ7C/a77BPTEAMY6PHEh8dz6CwQbib3ZsrbBGnUiVZEZEmVGsFV0sRliPF0tx8d7N0d6W9zd3sztCwofYKrp0COjV/0CKtjBIUEZET2Cu4eh3EzW87bn47sXg7TgO2VfszNmoM15xzESMiR+Dr7uukaEVaJyUoIiJHlFSVsP7gel7/9TN8e6zD7Fbs0G4tiz5WwbU8kolDBnNBZ1VwFWkKSlBEpE3bX7ifpLQkktOS2ZS1iWpbzQDXmgqunlSX9Kip4FrcC8Pq73CsKriKNB0lKCLSplRZq9iSvYWktCTWpK3h98LfHdo7B3RmTNRY/rvGn5xDkRi1/DepCq4iTU8Jioi0ejllOaxNX0tyWjLrDq6jpOrY4ntuJjeGhA9hbFTNANeYwBgABviogquIMylBEZFWxzAMtudtt/eSbM3ZinFcmhHsFcyYqDHEd4xnZMRI/Dz8TjqHKriKOJcSFBFpFUqrSlmfsZ41aWtITkvmUNkhh/bYkFh7bZLYkFjMJvMZz6kKriLOowRFRFxGfcvKHyg8QHJ6TQXXHzJ/oMpWZW/zdvNmZMRI4jvGMzpqNB18OjQoJlVwFXEOJSgi4hLqUla+ylZFSnYKSQeSSE5PJrUg1eEc0X7RxHeMZ2zUWIaGD8XD4tGsr0FEGo8SFBFxulOVlc8sKGf228ncclEF+fzMuvR1FFUV2dvdTG4MChtEfHQ8Y6LH0CWgS50W3xMR16cERUSc6uSy8gZmzwx7BVez9wFW/H6stZ1nO8ZEj2FM9BhGRY4iwENrZom0RkpQRMSpNqbmkVFYiMVvT806N347MLsXOuxjLY/gsh7juaHfBPqG9MVitjgpWhFpLkpQRMQp0orSSE5L5r/bv8SvZwomc7W9zbC5U13SHeuRsvJGdSAjhw9kQHuVlRdpK5SgiEizqLZVk5KdUjPr5kAyvxX8Zm8zmcFW2c6+zo21tCsY7g7Hq6y8SNuiBEVEmszh8sOsTV/LmrQ1rD24lqLKYwNcLSYLAzsMZHTUGF5c5cmh3CAMTh7gqrLyIm2TEhQRaTSGYbDr8C6S02pqk/yc8zM2w2ZvD/QMZHTUaOKj4xkVOYpAz0AAok0qKy8ijlpdgnLgwAGmTp1KdnY2bm5uPPjgg1x77bXODkuk1SqrLmNjxsaapCQ9mcySTIf2nu162iu49gvtV+sAV5WVF5ETmQzDOLH0QIuWkZFBVlYWAwcOJDMzkyFDhrBr1y58fX3rdHxhYSGBgYEUFBQQEKDpiyK1OVh80N5LsjFzIxXWCnubp8WTEREjGBs9ljFRY4jwq3tyUd9KsiLiWhrzO7TV9aBEREQQEVHzH2J4eDihoaHk5eXVOUERkZNV26r5+dDPJKclk5SWxJ78PQ7tEb4RjI2uWQ04LjwOL7eGDWhVWXkROcrlEpTk5GSefPJJNm/eTEZGBh988AFXXHGFwz6JiYk8+eSTZGZmMmDAAJ5//nni4uJOOtfmzZuxWq107NixmaIXcV317Z0oqChgbfpaktOSWZu+lsLKY7VJzCYzA9sPZEz0GMZGj6VHUA9VcBWRRuVyCUpJSQkDBgxgxowZXHXVVSe1r1ixgrlz57JkyRKGDx/OM888w4QJE9i5cycdOhxbDCwvL49p06axdOnS5gxfxCXVZZ0bwzDYk7/Hfusm5VCKwwDXAI8Azo06l/joeM6NPJcgr6Dmfhki0oa49BgUk8l0Ug/K8OHDGTZsGP/+978BsNlsdOzYkTvvvJMHHngAgIqKCi688EJmzpzJ1KlTT/scFRUVVFQcu39eWFhIx44dNQZFWo1TrXNjAjBVcdelJkosv7AmbQ0HSw467NM9qLt9gGv/9v1xM7vc7zQi4kLa7BiUyspKNm/ezPz58+3bzGYz48ePZ/369UDNb4E33XQT559//hmTE4DHHnuMhx9+uMliFnGmk9e5AZNbwZGS8tux+P7GK3uq7G0eZg/iIuKIj45nbPRYIv0imz9oERFaWIKSk5OD1WolLCzMYXtYWBg7duwA4LvvvmPFihX079+fDz/8EIDXX3+dfv361XrO+fPnM3fuXPvjoz0oIq3BxtQ8MgpKMXsfsK9zY/HKcNjHVhXAeR3juTb2IuLC4/Bx93FStCIix7SoBKUuRo8ejc1mO/OOR3h6euLp6dmEEYk0v8LKQtalr+ONXz/Dt8f3mN1K7W2GYcJW1pHq4nOoLu6FrSKCi4YO4ryOWudGRFxHi0pQQkNDsVgsZGVlOWzPysoiPDzcSVGJOJ9hGOwt2GufBpySnYLVsAJgdgPD6kV1cc+adW5KemJY/RyO1zo3IuJqWlSC4uHhwZAhQ1i9erV94KzNZmP16tXccccdZ3XuxMREEhMTsVqtjRCpSNOrsFawKXMTSWlJJKclk16c7tDeNbArY6LG8m6yH4cORWBwcgVXrXMjIq7K5RKU4uJi9uw5VgQqNTWVlJQUgoOD6dSpE3PnzmX69OkMHTqUuLg4nnnmGUpKSrj55pvP6nkTEhJISEiwj0AWcUVZJVmsSV9DUloSGzI2UFZdZm9zN7sTFx5nL5gW7R8NQF9vrXMjIi2Py00z/vbbbxk3btxJ26dPn87y5csB+Pe//20v1DZw4ECee+45hg8f3ijPr1L34kpsho2tOVtJSktiTdoatudtd2jv4N3BXixtRMSIUw5wrUsdFBGRs9WY36Eul6A4mxIUcbaiyiLWHVxnr+CaV55nbzNhol9oP3svSe/g3nWu4Kp1bkSkqbXZOihNSWNQxFkMw+D3wt/tFVy3ZG2h2qi2t/u5+zEqchRjo8cyOmo0Id4NW6tG69yISEuiHpQTqAdFmkOltZJNWZvsScmBogMO7TEBMfYKroPCBuFudndSpCIidaceFJEW6FDpIdakryE5LZn1B9dTWn2sNomb2Y1hYcPst246BXRyYqQiIs6nBEWkHuozjsNm2NiWu80+DXhb7jaH9lDvUMZEjSE+Op4RkSPwdfdtjpcgItIiKEERqaO6zIQprixmfcZ6ktOSWZO2htzyXIdz9A3pW9NL0nEs5wSfg9lkbtbXICLSUihBOUKDZOV0TrUicGZBOXNWfMENWcVkWX9kc9Zmqm3HBrj6uPlwbtS5jIkaw5joMYR6hzZv4CIiLZQGyZ5Ag2TlRFabwegnvj6u56Qai8/v9sX3zJ45Dvt38u9UM8C1YzxDOgzB3aIBriLSNmiQrEgz2piaR2bxIdwCd9YkJb67MVkq7O2GYcZa2oUb+k5k2oCJxATGOC9YEZFWQgmKSC1sho3tedtJTkvmo51f4ddzl2N7tR/W4l5UF/emuqQH2LwYMGogMYFaEVhEpDEoQRE5oqSqhO8Pfk9yes0A10NlhxzarWVRNQlJcW9s5VGA4wBXrQgsItJ4lKAcoUGybdOBwgMkpyeTdCCJTVmbqLJV2du83bwZGTGSMVFjefJDE9mHPU8aJAtNsyJwfn4+48ePp7q6murqau6++25mzpx50n6ffPIJ9913Hzabjfvvv59bb721Tm0iIq5Og2RPoEGyrVuVrYqU7BSSDiSRnJ5MakGqQ3u0XzTxHeMZGzWWoeFD8bB4AMdm8UDtKwIvvnFwoy66Z7VaqaiowMfHh5KSEvr27cumTZsICTlWqr66uprY2Fi++eYbAgMDGTJkCOvWrSMkJOS0bSIiTUWDZEXqIa88j7Xpa0k6kMT6g+spqiqyt7mZ3BgUNoj46HjGRI+hS0CXWhffm9g3gsU3Dj6pDkp4E60IbLFY8PGpWZm4oqICwzA48XeJjRs30qdPH6Kiasa9TJo0iS+++IIbbrjhtG0iIi2BEhRpdQzDYEfejpp1btKT+eXQLxjH9Xu082zHmOiauiSjIkcR4FG3LH9i3wgujA1vthWB8/PziY+PZ/fu3Tz55JOEhjrWUDl48KA9AQGIiooiPT39jG0iIi2BEhRpFUqrStmQsYGktCTWpK8huzTbob13cG/7Ojd9Q/piMVsa9DzNuSJwUFAQP/30E1lZWVx11VVcc801hIWFNctzi4g4mxIUabHSitLsvSQ/ZPxApa3S3ubt5s3wiOGMjR7LmKgxhPuGOzHS2l188cW0b9+eV199FYBvvvmGa6+9lqysLCyWYwlUWFgYAwYMYM2aNVxzzTX27ZGRkQ69Iunp6cTFxZ2xTUSkJVCCIi1Gta2alOwUktOTST6QzG8Fvzm0R/lF2XtJhoUPw9Pi2azxLVy4kLfffpvBgwfzwgsvkJyczF//+lfuuusubrnllpP2j4qKIjX12CDd+Ph4ysrK+P777+nevTs+Pj74+/tTUFBAcnIys2fPdjg+Li6OrVu3kp6eTmBgIJ999hkPPvjgGdtERFoCJShHaJqxazpcfpi16WtZk7aGtQfXUlR5bICrxWRhYIeBNWXlo+PpGti11gGuzeGrr74iOzubzZs3s2jRIq644gry8/P5z3/+Q8+ePWs9JioqijVr1tgfm81mvL29yc7Oxt3dndtuu80+OPbOO++kX79+AAwcOJCUlBTc3Nx46qmnGDduHDabjT/96U/2WTqnaxMRaQk0zfgEmmbsXIZhsOvwrppbN2nJ/JzzMzbDZm8P9AxkTNQYxkaPZVTkKAI9A50Y7TFPPvkko0ePZuTIkQAMHz6cu+66iylTppzymKVLlzJ37lyKimqSrpSUFIYOHUpaWhrh4a53S0pE5Ew0zVhcntVm1Hm2S1l1GRszNtrHk2SWZDq092zX095L0i+0X4MHuDal3r17s3LlSkaOHMmqVaswmUw88cQTnH/++URE1D4FOSoqiuLiYgoLC/Hz8+Pee+9lypQpSk5ERFCCIk1g1daMk+qFRJxQL+Rg8UF7L8nGzI1UWI8tvudl8XIY4Brh17g1RprC5MmTWb16NZ06dSImJob333+fr7/+mhEjRvDXv/611iqwR6cBp6WlsWzZMjIzM/noo4+aO3QREZekWzwn0C2es3O04uqJHyoTVszeB7hkxGEOlG9mT/4eh/YI3wj7ANe48Di83Fr/ujY5OTm0b9+eiRMnsmvXLpKTkx1ql4iItDS6xSMuyWozePjjbceSE3Mpbn67cPPbgZvvLkxupXxz5O6N2WRmYPuBjIkeQ3x0PN2DujttgKuzhIaG4unpyb59+0hKSlJyIiJyHCUo0mg27M0lq/x3PEJ2YPHbgcV7HybTsb4Uw+pNdXFPZg2bzE2DJhDkFeS8YF1EeXn5mXcSEWmDlKDIWSmvLmdjZs0A18/3foNvV8cKrtbyMKqLe2MtPgdrWUfAQhfvgUpORETktJSgHKE6KHWXWZJpH+C6IWMD5dZjvQCGzQ1raTeqi3tTXdQbo7rdScd38G/940tEROTsaJDsCTRI9mRWm5Vfcn6xJyU7D+90aA/zCTsy42Ys898qIyvfdtIgWQATNav/rr3//CZbYE9ERJxHg2SlyRVWFrIufR1JaUmsTV9LfkW+vc2EiQHtB9hn3fRs19M+wPWhS2tm8ZjAIUk5mo4snByr5ERERM5ICYoANRVc9xbsJTktmaS0JFKyU7Aax253+Xv4c27kuYyNHsvoqNG08zr51g3AxL4RLL5x8El1UMJPqIMiIiJyOkpQ2rAKawWbMjeRlJZEcloy6cXpDu3dArvV3LqJHsPADgNxN7vX6bwT+0ZwYWx4nSvJioiInEgJShuTVZLFmvQ1JKUlsSFjA2XVZfY2d7M7ceFx9ls30f7RDX4ei9nEyG5anE5ERBpGCUorZ7VZ2Zq7leS0ZNakrWF73naH9g7eHezF0oZHDMfH3cdJkYqIiByjBKUVKqosYt3BdSSnJbM2fS155Xn2NhMm+oX2s/eS9A7u3eYquIqIiOtTgtIKGIZBamEqa9Jqbt38mPUj1Ua1vd3P3Y9RkaPsA1xDvHXrRUREXJsSlBaq0lrJpqxN9tokB4oOOLTHBMQQHx3P2OixDAobVOcBriIiIq5ACcoRLaGS7KHSQzUDXA8ksT5j/UkDXIeGDbXfuukU0MmJkYqIiJwdVZI9gStVkrUZNn7N+ZXk9Jpekm252xzaQ71DaxKSqLGMiByBr7uvkyIVERFRJdkWxWoz6lUPpLiymPUZ60k6UFPBNbc816G9b0hfxnas6SU5J/gczCZzU78EERGRZqcEpQmt2ppxUkXViFoqqv5e8HvNWJL0ZDZnbabadmyAq6+7L6MiRzEmagxjoscQ6h3arK9BRETEGZSgNJFVW2vWpDnx/llmQTmz39jI3MvcKXX7hTXpa9hXuM9hn84Bne1jSYZ0GIK7RQNcRUSkbVGC0gSsNoOHP97muFiepQiL307c/Hbg5rubF3dX2NvczG4MCRvC2KiapCQmMKbZYxYREXElSlCawMbUPPttHYvPHjw7rMLineawj63aj9GRY7gm9kJGRozEz8PPGaGKiIi4JCUoTSC76NiYEww3e3JiLYuiurg31cW9sZVHcfGQwVzYOcpJUYqIiLguJShNoIO/l/3v1rKOlB28BmtJT4zqgFPuJyIiIsdojmoTiOsSTESgFzWTiS1UFwx1SE5M1MzmiesS7KQIRUREXJsSlCZgMZtYODkWgBMrnhx9vHBy7GnroYiIiLRlSlCayMS+ESy+cTDhgY63ccIDvVh842CHOigiIiLiSGNQmtDEvhFcGBter0qyIiIiogSlyVnMJkZ2C3F2GCIiIi2KbvGIiIiIy1GCckRiYiKxsbEMGzbM2aGIiIi0eSbDME5cLqZNa8ylokVERNqSxvwOVQ+KiIiIuBwlKCIiIuJylKCIiIiIy1GCIiIiIi5HCYqIiIi4HCUoIiIi4nKUoIiIiIjLUan7ExwtC1NYWOjkSERERFqWo9+djVFiTQnKCYqKigDo2LGjkyMRERFpmYqKiggMDDyrc6iS7AlsNhsHDx7E398fk6ltrDpcWFhIx44dOXDggKrnujBdp8ah9/GY1v5etIbX15Jew9FYt23bRq9evTCbz24UiXpQTmA2m4mOjnZ2GE4REBDg8v8ARNepseh9PKa1vxet4fW1pNcQFRV11skJaJCsiIiIuCAlKCIiIuJylKAInp6eLFy4EE9PT2eHIqeh69Q49D4e09rfi9bw+lrSa2jsWDVIVkRERFyOelBERETE5ShBEREREZejBEVERERcjhIUERERcTlKUNqQxYsX079/f3vBn5EjR/LZZ5/Z28vLy0lISCAkJAQ/Pz+uvvpqsrKynBixPP7445hMJu655x77Nl2nM3vooYcwmUwOP71797a3t7X3MD09nRtvvJGQkBC8vb3p168fmzZtsrcbhsGCBQuIiIjA29ub8ePHs3v3bidGXD8xMTEnXW+TyURCQgLg+tfbarXy4IMP0qVLF7y9venWrRt///vfHdazcbVrVFRUxD333EPnzp3x9vZm1KhR/PDDD40bryFtxsqVK41PP/3U2LVrl7Fz507jz3/+s+Hu7m5s3brVMAzDmDVrltGxY0dj9erVxqZNm4wRI0YYo0aNcnLUbdfGjRuNmJgYo3///sbdd99t367rdGYLFy40+vTpY2RkZNh/Dh06ZG9vS+9hXl6e0blzZ+Omm24yNmzYYOzdu9f4/PPPjT179tj3efzxx43AwEDjww8/NH766SfjsssuM7p06WKUlZU5MfK6y87OdrjWX375pQEY33zzjWEYrn+9H330USMkJMT45JNPjNTUVOM///mP4efnZzz77LP2fVztGl133XVGbGyskZSUZOzevdtYuHChERAQYKSlpTVavEpQ2rh27doZL730kpGfn2+4u7sb//nPf+xt27dvNwBj/fr1ToywbSoqKjJ69OhhfPnll0Z8fLw9QdF1qpuFCxcaAwYMqLWtrb2H999/vzF69OhTtttsNiM8PNx48skn7dvy8/MNT09P4+23326OEBvd3XffbXTr1s2w2Wwt4npfcsklxowZMxy2XXXVVcaUKVMMw3C9a1RaWmpYLBbjk08+cdg+ePBg4y9/+UujxatbPG2U1WrlnXfeoaSkhJEjR7J582aqqqoYP368fZ/evXvTqVMn1q9f78RI26aEhAQuueQSh+sB6DrVw+7du4mMjKRr165MmTKF/fv3A23vPVy5ciVDhw7l2muvpUOHDgwaNIilS5fa21NTU8nMzHR4PwIDAxk+fHiLfD8qKyt54403mDFjBiaTqUVc71GjRrF69Wp27doFwE8//cTatWuZNGkS4HrXqLq6GqvVipeXl8N2b29v1q5d22jxarHANuaXX35h5MiRlJeX4+fnxwcffEBsbCwpKSl4eHgQFBTksH9YWBiZmZnOCbaNeuedd9iyZYvD/dyjMjMzdZ3qYPjw4SxfvpxevXqRkZHBww8/zJgxY9i6dWubew/37t3L4sWLmTt3Ln/+85/54YcfuOuuu/Dw8GD69On21xwWFuZwXEt9Pz788EPy8/O56aabgJbxb+aBBx6gsLCQ3r17Y7FYsFqtPProo0yZMgXA5a6Rv78/I0eO5O9//zvnnHMOYWFhvP3226xfv57u3bs3WrxKUNqYXr16kZKSQkFBAf/973+ZPn06SUlJzg5Ljjhw4AB33303X3755Um/nUjdHf3NE6B///4MHz6czp078+677+Lt7e3EyJqfzWZj6NCh/OMf/wBg0KBBbN26lSVLljB9+nQnR9f4Xn75ZSZNmkRkZKSzQ6mzd999lzfffJO33nqLPn36kJKSwj333ENkZKTLXqPXX3+dGTNmEBUVhcViYfDgwdxwww1s3ry50Z5Dt3jaGA8PD7p3786QIUN47LHHGDBgAM8++yzh4eFUVlaSn5/vsH9WVhbh4eHOCbYN2rx5M9nZ2QwePBg3Nzfc3NxISkriueeew83NjbCwMF2nBggKCqJnz57s2bOnzX3WIyIiiI2Nddh2zjnn2G95HX3NJ85qaYnvx759+/jqq6+49dZb7dtawvX+v//7Px544AH+8Ic/0K9fP6ZOncq9997LY489BrjmNerWrRtJSUkUFxdz4MABNm7cSFVVFV27dm20eJWgtHE2m42KigqGDBmCu7s7q1evtrft3LmT/fv3M3LkSCdG2LZccMEF/PLLL6SkpNh/hg4dypQpU+x/13Wqv+LiYn777TciIiLa3Gf93HPPZefOnQ7bdu3aRefOnQHo0qUL4eHhDu9HYWEhGzZsaHHvx7Jly+jQoQOXXHKJfVtLuN6lpaWYzY5fxxaLBZvNBrj2NfL19SUiIoLDhw/z+eefc/nllzdevI02rFdc3gMPPGAkJSUZqampxs8//2w88MADhslkMr744gvDMGqm4nXq1Mn4+uuvjU2bNhkjR440Ro4c6eSo5fhZPIah61QX9913n/Htt98aqampxnfffWeMHz/eCA0NNbKzsw3DaFvv4caNGw03Nzfj0UcfNXbv3m28+eabho+Pj/HGG2/Y93n88ceNoKAg46OPPjJ+/vln4/LLL29R04wNwzCsVqvRqVMn4/777z+pzdWv9/Tp042oqCj7NOP333/fCA0NNf70pz/Z93G1a7Rq1Srjs88+M/bu3Wt88cUXxoABA4zhw4cblZWVjRavEpQ2ZMaMGUbnzp0NDw8Po3379sYFF1xgT04MwzDKysqMOXPmGO3atTN8fHyMK6+80sjIyHBixGIYJycouk5ndv311xsRERGGh4eHERUVZVx//fUOdT/a2nv48ccfG3379jU8PT2N3r17Gy+++KJDu81mMx588EEjLCzM8PT0NC644AJj586dToq2YT7//HMDqDVuV7/ehYWFxt1332106tTJ8PLyMrp27Wr85S9/MSoqKuz7uNo1WrFihdG1a1fDw8PDCA8PNxISEoz8/PxGjddkGMeVqhMRERFxARqDIiIiIi5HCYqIiIi4HCUoIiIi4nKUoIiIiIjLUYIiIiIiLkcJioiIiLgcJSgiIiLicpSgiIiIiMtRgiIiIiIuRwmKtHiGYXDbbbcRHByMyWQiJSWl1m1N5bzzzuOee+5psvOfjcaOrSleqyu/f3Jqum7S1JSgiMu76aabMJlMJ/1MnDgRgFWrVrF8+XI++eQTMjIy6Nu3b63bztap/kN+//33+fvf/37W5z+d498DDw8Punfvzt/+9jeqq6tPe1xjx9Ycr7U2Bw4cYMaMGURGRuLh4UHnzp25++67yc3NbfZYwDW+nI9+Jh5//HGH7R9++CEmk8lJUYk0HjdnByBSFxMnTmTZsmUO2zw9PQH47bffiIiIYNSoUfa22rY1leDg4CZ/Djj2HlRUVPC///2PhIQE3N3dmT9//kn7VlZW4uHh0eixNddrPd7evXsZOXIkPXv25O2336ZLly78+uuv/N///R+fffYZ33//vVPicgVeXl488cQT3H777bRr187Z4TSKo59dEfWgSIvg6elJeHi4w0+7du246aabuPPOO9m/fz8mk4mYmJhatwHYbDYee+wxunTpgre3NwMGDOC///2vw/PYbDb++c9/0r17dzw9PenUqROPPvooN910E0lJSTz77LP2nozff/8dOPbb9IsvvkhkZCQ2m83hnJdffjkzZsxweI4zxXG696Bz587Mnj2b8ePHs3LlSnsMd9xxB/fccw+hoaFMmDDBIbajzjvvPO666y7+9Kc/ERwcTHh4OA899NAZX//xx594vjvuuIM77riDwMBAQkNDefDBBzl+DdJVq1YxevRogoKCCAkJ4dJLL+W333474+s9KiEhAQ8PD7744gvi4+Pp1KkTkyZN4quvviI9PZ2//OUv9n1jYmJ45plnHI4fOHCgw2usSzyne59O91moy/Ofd9553Hnnndxzzz20a9eOsLAwli5dSklJCTfffDP+/v50796dzz777Izvzfjx4wkPD+exxx475T5NGVN1dfVpr31dPuun+uye6OKLL2b69On2x9988w2hoaFYrdYzvU3SQilBkRbt2Wef5W9/+xvR0dFkZGTwww8/1LoN4LHHHuO1115jyZIl/Prrr9x7773ceOONJCUl2c83f/58Hn/8cR588EG2bdvGW2+9RVhYGM8++ywjR45k5syZZGRkkJGRQceOHR1iufbaa8nNzeWbb76xb8vLy2PVqlVMmTLFvq0ucdSFt7c3lZWV9sevvvoqHh4efPfddyxZsuSUx7366qv4+vqyYcMG/vnPf/K3v/2NL7/88rSv/3ReffVV3Nzc2LhxI88++yxPP/00L730kr29pKSEuXPnsmnTJlavXo3ZbObKK688KZGrTV5eHp9//jlz5szB29vboS08PJwpU6awYsUK6rMoe13jOdX7VJfPwpm8+uqrhIaGsnHjRu68805mz57Ntddey6hRo9iyZQsXXXQRU6dOpbS09LTnsVgs/OMf/+D5558nLS2tXjE0RkxnuvZ1/azX5bMbFRVFenq6/XF8fDxlZWV8//33Z/W6xYUZIi5u+vTphsViMXx9fR1+Hn30UcMwDONf//qX0blzZ4djTtxWXl5u+Pj4GOvWrXPY75ZbbjFuuOEGwzAMo7Cw0PD09DSWLl1aaxzx8fHG3Xfffdrtl19+uTFjxgx72wsvvGBERkYaVqu1znGc6j24/PLLDcMwDJvNZnz55ZeGp6enMW/ePHsMgwYNOmPM8fHxxujRox32GTZsmHH//fef8fWf6nznnHOOYbPZ7Nvuv/9+45xzzjnlOQ4dOmQAxi+//FLrOY/3/fffG4DxwQcf1Nr+9NNPG4CRlZVlGIZhdO7c2fjXv/7lsM+AAQOMhQsX1jmeozGd6n06Xcx1ef4Tz11dXW34+voaU6dOtW/LyMgwAGP9+vWnjPv4z8SIESPsn7sPPvjAOP6/9qaK6UzXvq6f9VN9dk+0cOFCo1evXg7bQkJCjPfff/+Mx0rLpDEo0iKMGzeOxYsXO2yrz7iDPXv2UFpayoUXXuiwvbKykkGDBgGwfft2KioquOCCCxoc55QpU5g5cyaLFi3C09OTN998kz/84Q+YzeY6x3Eqn3zyCX5+flRVVWGz2fjjH//o0E0/ZMiQOsXYv39/h8cRERFkZ2c3+PWPGDHCYVDmyJEjeeqpp7BarVgsFnbv3s2CBQvYsGEDOTk59p6K/fv313nwsnGGHpL6jFmoazynep8aw/HntlgshISE0K9fP/u2o71WdX2+J554gvPPP5958+Y1a0ynu/b1+azX5bN7Yg9KSkoK+fn5jBw5so6vUFoaJSjSIvj6+tK9e/cGH19cXAzAp59+SlRUlEPb0cG2J95CaIjJkydjGAaffvopw4YNY82aNfzrX/+qVxyncjRJ8/DwIDIyEjc3x3++vr6+dYrR3d3d4bHJZMJmszXK66/N5MmT6dy5M0uXLrWP0enbt6/D7alT6d69OyaTie3bt3PllVee1L59+3bat29PUFAQAGaz+aRkpqqqqkHxnOp9Op26PP+pzn38tqNf+nW5DQYwduxYJkyYwPz587nppptcIqb6fNbr8tmNioqiuLiYwsJC/Pz8uPfee5kyZQrh4eF1ikdaHiUo0ibExsbi6enJ/v37iY+Pr3WfHj164O3tzerVq7n11ltPavfw8DjjgDwvLy+uuuoq3nzzTfbs2UOvXr0YPHhwveI4lbNN0s7kTK//VDZs2ODw+Pvvv6dHjx5YLBZyc3PZuXMnS5cuZcyYMQCsXbu2zucOCQnhwgsvZNGiRdx7770OSVRmZiZvvvkmCQkJ9m3t27cnIyPD/riwsJDU1FT747ON56hTfRbO9PxN6fHHH2fgwIH06tWr2WI63bU/m896bY4mOWlpaSxbtozMzEw++uijsz6vuC4lKNIiVFRUkJmZ6bDNzc2N0NDQOh3v7+/PvHnzuPfee7HZbIwePZqCggK+++47AgICmD59Ol5eXtx///386U9/wsPDg3PPPZdDhw7x66+/cssttxATE8OGDRv4/fff8fPzIzg42H7r5nhTpkzh0ksv5ddff+XGG2+sdxzOcqbXfyr79+9n7ty53H777WzZsoXnn3+ep556CoB27doREhLCiy++SEREBPv37+eBBx6oV1z//ve/GTVqFBMmTOCRRx5xmGbcs2dPFixYYN/3/PPPZ/ny5UyePJmgoCAWLFiAxWKxtzdGPMApPwtnev6m1K9fP6ZMmcJzzz3nsL0pYzrdtW/sz/rRBOW+++5j165dJCcnExAQ0CivQ1yTEhRpEVatWkVERITDtl69erFjx446n+Pvf/877du357HHHmPv3r0EBQUxePBg/vznP9v3efDBB3Fzc2PBggUcPHiQiIgIZs2aBcC8efOYPn06sbGxlJWVkZqaap/CfLzzzz+f4OBgdu7cyR//+McGxeEsp3v9pzJt2jTKysqIi4vDYrFw9913c9tttwE1txfeeecd7rrrLvr27UuvXr147rnnOO+88+ocU48ePfjhhx946KGHuO6668jOzsYwDK666ipef/11fHx87PvOnz+f1NRULr30UgIDA/n73//u0FvQGPHAqT8LZ3r+pva3v/2NFStWOGxryphOd+2hcT/roaGheHp6sm/fPpKSkk66bSStj8k40+gzEZFTOO+88xg4cOBJdTaa2sKFC3n66af58ssvGTFiRLM+t4g0D/WgiEiL8/DDDxMTE8P3339PXFxcrbfaRKRlU4IiIi3SzTff7OwQRKQJ6RaPiIiIuBz1i4qIiIjLUYIiIiIiLkcJioiIiLgcJSgiIiLicpSgiIiIiMtRgiIiIiIuRwmKiIiIuBwlKCIiIuJylKCIiIiIy1GCIiIiIi5HCYqIiIi4nP8PefvJ4rcPYuoAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "n_list = list(range(30, 90, 5))\n", + "\n", + "\n", + "def fit_function(x: np.ndarray, a: float, b: float) -> np.typing.NDArray[Any]:\n", + " return a * x + b\n", + "\n", + "\n", + "# Calculate lifetimes for S states\n", + "kets_s = [rydstate.RydbergStateSQDTAlkali(\"Rb\", n=n, l=0, j=0.5, m=0.5) for n in n_list]\n", + "nu_s = [ket.nu for ket in kets_s]\n", + "lifetimes_s = [ket.get_lifetime(unit=\"mus\") for ket in kets_s]\n", + "popt_s, _ = curve_fit(fit_function, np.log(nu_s), np.log(lifetimes_s))\n", + "\n", + "# Calculate lifetimes for circular states\n", + "kets_circular = [rydstate.RydbergStateSQDTAlkali(\"Rb\", n=n, l=n - 1, j=n - 0.5, m=n - 0.5) for n in n_list]\n", + "nu_circular = [ket.nu for ket in kets_circular]\n", + "lifetimes_circular = [ket.get_lifetime(unit=\"mus\") for ket in kets_circular]\n", + "popt_circular, _ = curve_fit(fit_function, np.log(nu_circular), np.log(lifetimes_circular))\n", + "\n", + "# Plot the scaling of the lifetimes\n", + "fig, ax = plt.subplots(figsize=(6, 5))\n", + "\n", + "ax.plot(nu_s, lifetimes_s, \"o\", label=\"S states\")\n", + "ax.plot(nu_circular, lifetimes_circular, \"o\", label=\"Circular states\")\n", + "\n", + "fit_s = np.exp(fit_function(np.log(nu_s), *popt_s))\n", + "fit_circular = np.exp(fit_function(np.log(nu_circular), *popt_circular))\n", + "ax.plot(nu_s, fit_s)\n", + "ax.plot(nu_circular, fit_circular)\n", + "\n", + "ax.text(nu_s[2], fit_s[2], rf\"$\\propto \\nu^{{{popt_s[0]:.2f}}}$\", verticalalignment=\"top\")\n", + "ax.text(\n", + " nu_circular[2],\n", + " fit_circular[2],\n", + " rf\"$\\propto \\nu^{{{popt_circular[0]:.2f}}}$\",\n", + " verticalalignment=\"top\",\n", + ")\n", + "\n", + "ax.legend()\n", + "ax.set_yscale(\"log\")\n", + "ax.set_xscale(\"log\")\n", + "ax.set_xlabel(r\"Effective Principal Quantum Number $\\nu$\")\n", + "ax.set_ylabel(r\"Lifetime ($\\mu$s)\")\n", + "\n", + "ax.set_xticks([30, 40, 50, 60, 70, 80, 90])\n", + "ax.get_xaxis().set_major_formatter(plt.ScalarFormatter())\n", + "\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "rydstate", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.2" + } + }, + "nbformat": 4, + "nbformat_minor": 1 +} diff --git a/src/rydstate/rydberg/rydberg_sqdt.py b/src/rydstate/rydberg/rydberg_sqdt.py index f99c1ee..5f7e5ce 100644 --- a/src/rydstate/rydberg/rydberg_sqdt.py +++ b/src/rydstate/rydberg/rydberg_sqdt.py @@ -6,6 +6,7 @@ from typing import TYPE_CHECKING, overload import numpy as np +from scipy.special import exprel from rydstate.angular import NotSet from rydstate.angular.utils import quantum_numbers_to_angular_ket @@ -19,7 +20,7 @@ from typing_extensions import Self from rydstate.angular.angular_ket import AngularKetBase, AngularKetFJ, AngularKetJJ, AngularKetLS - from rydstate.units import MatrixElementOperator, PintFloat + from rydstate.units import MatrixElementOperator, NDArray, PintArray, PintFloat logger = logging.getLogger(__name__) @@ -319,6 +320,210 @@ def calc_matrix_element( reduced_matrix_element = self.calc_reduced_matrix_element(other, operator, unit) return prefactor * reduced_matrix_element + @overload + def get_spontaneous_transition_rates(self, unit: None = None) -> tuple[list[RydbergStateSQDT], PintArray]: ... + + @overload + def get_spontaneous_transition_rates(self, unit: str) -> tuple[list[RydbergStateSQDT], NDArray]: ... + + def get_spontaneous_transition_rates( + self, unit: str | None = None + ) -> tuple[list[RydbergStateSQDT], NDArray | PintArray]: + """Calculate the spontaneous transition rates for the Rydberg state. + + The spontaneous transition rates are given by the Einstein A coefficients. + + Args: + unit: The unit to which to convert the result. + Default None will return a `pint.Quantity`. + + Returns: + The relevant states and the transition rates. + + """ + relevant_states_masked, transition_rates_au = self._get_transition_rates_au(only_spontaneous=True) + + if unit == "a.u.": + return relevant_states_masked, transition_rates_au + transition_rates = ureg.Quantity(transition_rates_au, "1/atomic_unit_of_time") + if unit is None: + return relevant_states_masked, transition_rates + return relevant_states_masked, transition_rates.to(unit).magnitude + + @overload + def get_black_body_transition_rates( + self, temperature: float | PintFloat, temperature_unit: str | None = None, unit: None = None + ) -> tuple[list[RydbergStateSQDT], PintArray]: ... + + @overload + def get_black_body_transition_rates( + self, temperature: PintFloat, *, unit: str + ) -> tuple[list[RydbergStateSQDT], NDArray]: ... + + @overload + def get_black_body_transition_rates( + self, temperature: float, temperature_unit: str, unit: str + ) -> tuple[list[RydbergStateSQDT], NDArray]: ... + + def get_black_body_transition_rates( + self, temperature: float | PintFloat, temperature_unit: str | None = None, unit: str | None = None + ) -> tuple[list[RydbergStateSQDT], NDArray | PintArray]: + """Calculate the black body transition rates of the Rydberg state. + + The black body transition rates are given by the Einstein B coefficients, + with a weight factor given by Planck's law. + + Args: + temperature: The temperature, for which to calculate the black body transition rates. + temperature_unit: The unit of the temperature. + Default None will assume the temperature is given as `pint.Quantity`. + unit: The unit to which to convert the result. + Default None will return a `pint.Quantity`. + + Returns: + The relevant states and the transition rates. + + """ + temperature_au = ureg.Quantity(temperature, temperature_unit).to_base_units().magnitude + relevant_states_masked, transition_rates_au = self._get_transition_rates_au( + temperature_au, only_spontaneous=False + ) + + if unit == "a.u.": + return relevant_states_masked, transition_rates_au + transition_rates = ureg.Quantity(transition_rates_au, "1/atomic_unit_of_time") + if unit is None: + return relevant_states_masked, transition_rates + return relevant_states_masked, transition_rates.to(unit).magnitude + + def _get_transition_rates_au( + self, + temperature_au: float | None = None, + *, + only_spontaneous: bool = False, + ) -> tuple[list[RydbergStateSQDT], NDArray]: + r"""Calculate the transition rates in atomic units. + + The transition rates are given by the Einstein A coefficients for spontaneous transitions, + and by the Einstein B coefficients with a weight factor given by Planck's law for black body transitions. + + Concretely the transition rates are calculated as + + .. math:: + \Gamma^{spontaneous}_{self \to other} = \frac{4}{3} \frac{\alpha}{c^2} \omega^3 + |\langle self || r^k_radial \hat{O}_{k_angular} || other \rangle|^2 + + and + + .. math:: + \Gamma^{blackbody}_{self \to other} = \Gamma^{spontaneous}_{self \to other} \frac{1}{\exp(\omega / T) - 1} + + """ + if self.species.number_valence_electrons == 2 and self.angular.coupling_scheme != "LS": + raise NotImplementedError( + "For alkaline earth atoms transition rates are only implemented for LS coupling scheme." + ) + from rydstate.basis import BasisSQDTAlkali, BasisSQDTAlkalineLS # noqa: PLC0415 + + basis_class = BasisSQDTAlkali if self.species.number_valence_electrons == 1 else BasisSQDTAlkalineLS + + m = self.angular.m + if isinstance(m, NotSet): + raise RuntimeError("m quantum number must be defined to calculate transition rates.") # noqa: TRY004 + + basis = basis_class(self.species, n=(1, int(self.nu + 35)), m=(m - 1, m + 1)) + basis.filter_states("l_r", (self.angular.l_r - 1, self.angular.l_r + 1)) + + if only_spontaneous: + basis.filter_states("nu", (0, self.nu)) + + relevant_states = basis.states + energy_differences_au = self.get_energy("hartree") - np.array( + [s.get_energy("hartree") for s in relevant_states] + ) + electric_dipole_moments_au = np.zeros(len(relevant_states)) + for q in [-1, 0, 1]: + # the different entries are only at most once nonzero -> we can just add the arrays + el_di_m = np.array( + [s.calc_matrix_element(self, "electric_dipole", q, unit="a.u.") for s in relevant_states] + ) + electric_dipole_moments_au += el_di_m + + transition_rates_au = ( + (4 / 3) + * np.abs(electric_dipole_moments_au) ** 2 + * energy_differences_au**2 + / ureg.Quantity(1, "speed_of_light").to_base_units().magnitude ** 3 + ) + + if only_spontaneous: + transition_rates_au *= energy_differences_au + else: + assert temperature_au is not None, "Temperature must be given for black body transitions." + if temperature_au == 0: + transition_rates_au *= 0 + else: # for numerical stability we use 1 / exprel(x) = x / (exp(x) - 1) + transition_rates_au *= temperature_au / exprel(energy_differences_au / temperature_au) + + mask = transition_rates_au != 0 + relevant_states_masked = [ket for ket, is_relevant in zip(relevant_states, mask, strict=True) if is_relevant] + transition_rates_au = transition_rates_au[mask] + return relevant_states_masked, transition_rates_au + + @overload + def get_lifetime( + self, + temperature: float | PintFloat | None = None, + temperature_unit: str | None = None, + unit: None = None, + ) -> PintFloat: ... + + @overload + def get_lifetime(self, *, unit: str) -> float: ... + + @overload + def get_lifetime(self, temperature: PintFloat, *, unit: str) -> float: ... + + @overload + def get_lifetime(self, temperature: float, temperature_unit: str, unit: str) -> float: ... + + def get_lifetime( + self, + temperature: float | PintFloat | None = None, + temperature_unit: str | None = None, + unit: str | None = None, + ) -> float | PintFloat: + """Calculate the lifetime of the Rydberg state. + + The lifetime is the inverse of the sum of all transition rates. + + Args: + temperature: The temperature, for which to calculate the black body transition rates. + Default None will not include black body transitions. + temperature_unit: The unit of the temperature. + Default None will assume the temperature is given as `pint.Quantity`. + unit: The unit to which to convert the result. + Default None will return a `pint.Quantity`. + + Returns: + The lifetime of the state. + + """ + _, transition_rates = self.get_spontaneous_transition_rates() + transition_rates_au = transition_rates.to_base_units().magnitude + if temperature is not None: + _, black_body_transition_rates = self.get_black_body_transition_rates(temperature, temperature_unit) + transition_rates_au = np.append(transition_rates_au, black_body_transition_rates.to_base_units().magnitude) + + lifetime_au: float = 1 / np.sum(transition_rates_au) + + if unit == "a.u.": + return lifetime_au + lifetime = ureg.Quantity(lifetime_au, "atomic_unit_of_time") + if unit is None: + return lifetime + return lifetime.to(unit).magnitude + class RydbergStateSQDTAlkali(RydbergStateSQDT): """Create an Alkali Rydberg state, including the radial and angular states.""" From b6ecabae786d9cc58c269b7f7db72654b1968604 Mon Sep 17 00:00:00 2001 From: johannes-moegerle Date: Mon, 9 Mar 2026 15:55:16 +0100 Subject: [PATCH 05/10] speed up check is not set --- src/rydstate/angular/angular_ket.py | 15 ++++++++------- src/rydstate/angular/angular_state.py | 6 +++--- src/rydstate/angular/utils.py | 9 ++++++++- src/rydstate/basis/basis_sqdt.py | 5 +++-- src/rydstate/rydberg/rydberg_sqdt.py | 6 +++--- 5 files changed, 25 insertions(+), 16 deletions(-) diff --git a/src/rydstate/angular/angular_ket.py b/src/rydstate/angular/angular_ket.py index fbedbc8..cf9380e 100644 --- a/src/rydstate/angular/angular_ket.py +++ b/src/rydstate/angular/angular_ket.py @@ -17,6 +17,7 @@ NotSet, check_spin_addition_rule, get_possible_quantum_number_values, + is_not_set, minus_one_pow, try_trivial_spin_addition, ) @@ -115,7 +116,7 @@ def __init__( self.l_r = int(l_r) # f_tot will be set in the subclasses - self.m = NotSet if isinstance(m, NotSet) else float(m) + self.m = NotSet if is_not_set(m) else float(m) def _post_init(self) -> None: self.quantum_numbers = tuple(getattr(self, qn) for qn in self.quantum_number_names) @@ -133,7 +134,7 @@ def sanity_check(self, msgs: list[str] | None = None) -> None: if self.s_r != 0.5: msgs.append(f"Rydberg electron spin s_r must be 1/2, but {self.s_r=}") - if not isinstance(self.m, NotSet) and not -self.f_tot <= self.m <= self.f_tot: + if not is_not_set(self.m) and not -self.f_tot <= self.m <= self.f_tot: msgs.append(f"m must be between -f_tot and f_tot, but {self.f_tot=}, {self.m=}") if msgs: @@ -150,7 +151,7 @@ def __setattr__(self, key: str, value: object) -> None: def __repr__(self) -> str: args = ", ".join(f"{qn}={val}" for qn, val in zip(self.quantum_number_names, self.quantum_numbers, strict=True)) - if not isinstance(self.m, NotSet): + if not is_not_set(self.m): args += f", m={self.m}" return f"{self.__class__.__name__}({args})" @@ -467,16 +468,16 @@ def calc_matrix_element(self, other: AngularKetBase, operator: AngularOperatorTy The dimensionless angular matrix element. """ - if isinstance(self.m, NotSet) or isinstance(other.m, NotSet): - raise RuntimeError("m must be set to calculate the matrix element.") # noqa: TRY004 + if is_not_set(self.m) or is_not_set(other.m): + raise RuntimeError("m must be set to calculate the matrix element.") prefactor = self._calc_wigner_eckart_prefactor(other, kappa, q) reduced_matrix_element = self.calc_reduced_matrix_element(other, operator, kappa) return prefactor * reduced_matrix_element def _calc_wigner_eckart_prefactor(self, other: AngularKetBase, kappa: int, q: int) -> float: - if isinstance(self.m, NotSet) or isinstance(other.m, NotSet): - raise RuntimeError("m must be set to calculate the Wigner-Eckart prefactor.") # noqa: TRY004 + if is_not_set(self.m) or is_not_set(other.m): + raise RuntimeError("m must be set to calculate the Wigner-Eckart prefactor.") return minus_one_pow(self.f_tot - self.m) * calc_wigner_3j(self.f_tot, kappa, other.f_tot, -self.m, q, other.m) def _kronecker_delta_non_involved_spins(self, other: AngularKetBase, qn: AngularMomentumQuantumNumbers) -> int: diff --git a/src/rydstate/angular/angular_state.py b/src/rydstate/angular/angular_state.py index 9d996e0..7bd5ff8 100644 --- a/src/rydstate/angular/angular_state.py +++ b/src/rydstate/angular/angular_state.py @@ -13,7 +13,7 @@ AngularKetLS, ) from rydstate.angular.angular_matrix_element import is_angular_momentum_quantum_number -from rydstate.angular.utils import NotSet +from rydstate.angular.utils import is_not_set if TYPE_CHECKING: from collections.abc import Iterator, Sequence @@ -210,8 +210,8 @@ def calc_matrix_element( "Different m values are not supported yet for AngularState.calc_matrix_element." ) - if isinstance(self.kets[0].m, NotSet) or isinstance(other.kets[0].m, NotSet): - raise RuntimeError("m must be set for all kets to calculate the matrix element.") # noqa: TRY004 + if is_not_set(self.kets[0].m) or is_not_set(other.kets[0].m): + raise RuntimeError("m must be set for all kets to calculate the matrix element.") prefactor = self.kets[0]._calc_wigner_eckart_prefactor(other.kets[0], kappa, q) # noqa: SLF001 reduced_matrix_element = self.calc_reduced_matrix_element(other, operator, kappa) diff --git a/src/rydstate/angular/utils.py b/src/rydstate/angular/utils.py index 5a2721c..b577eed 100644 --- a/src/rydstate/angular/utils.py +++ b/src/rydstate/angular/utils.py @@ -2,11 +2,13 @@ import contextlib import typing as t -from typing import TYPE_CHECKING, Literal +from typing import TYPE_CHECKING, Any, Literal import numpy as np if TYPE_CHECKING: + from typing_extensions import TypeIs + from rydstate.angular.angular_ket import AngularKetBase from rydstate.species.species_object import SpeciesObject @@ -26,6 +28,11 @@ class NotSet(t.Protocol): def __not_set() -> None: ... +def is_not_set(obj: Any) -> TypeIs[NotSet]: # noqa: ANN401 + """Check if the obj is the NotSet singleton.""" + return obj is NotSet + + class InvalidQuantumNumbersError(ValueError): def __init__(self, ket: AngularKetBase, msg: str = "") -> None: _msg = f"Invalid quantum numbers for {ket!r}" diff --git a/src/rydstate/basis/basis_sqdt.py b/src/rydstate/basis/basis_sqdt.py index 2248d40..ca5c735 100644 --- a/src/rydstate/basis/basis_sqdt.py +++ b/src/rydstate/basis/basis_sqdt.py @@ -5,7 +5,8 @@ import numpy as np -from rydstate.angular.utils import NotSet +from rydstate.angular import NotSet +from rydstate.angular.utils import is_not_set from rydstate.basis.basis_base import BasisBase from rydstate.rydberg import ( RydbergStateSQDT, @@ -23,7 +24,7 @@ class BasisSQDT(BasisBase[_RydbergStateSQDT]): def _get_m_range(self, m: tuple[float, float] | None | NotSet, f_tot: float | np.floating) -> list[NotSet | float]: - if isinstance(m, NotSet): + if is_not_set(m): return [NotSet] if m is None: m = (-np.inf, np.inf) diff --git a/src/rydstate/rydberg/rydberg_sqdt.py b/src/rydstate/rydberg/rydberg_sqdt.py index 5f7e5ce..39e89e7 100644 --- a/src/rydstate/rydberg/rydberg_sqdt.py +++ b/src/rydstate/rydberg/rydberg_sqdt.py @@ -9,7 +9,7 @@ from scipy.special import exprel from rydstate.angular import NotSet -from rydstate.angular.utils import quantum_numbers_to_angular_ket +from rydstate.angular.utils import is_not_set, quantum_numbers_to_angular_ket from rydstate.radial import RadialKet from rydstate.rydberg.rydberg_base import RydbergStateBase from rydstate.species import SpeciesObjectSQDT @@ -428,8 +428,8 @@ def _get_transition_rates_au( basis_class = BasisSQDTAlkali if self.species.number_valence_electrons == 1 else BasisSQDTAlkalineLS m = self.angular.m - if isinstance(m, NotSet): - raise RuntimeError("m quantum number must be defined to calculate transition rates.") # noqa: TRY004 + if is_not_set(m): + raise RuntimeError("m quantum number must be defined to calculate transition rates.") basis = basis_class(self.species, n=(1, int(self.nu + 35)), m=(m - 1, m + 1)) basis.filter_states("l_r", (self.angular.l_r - 1, self.angular.l_r + 1)) From 23815f74624a689e1261ac608dc30aff4fc2d13a Mon Sep 17 00:00:00 2001 From: johannes-moegerle Date: Mon, 9 Mar 2026 15:57:36 +0100 Subject: [PATCH 06/10] cache quantum_numbers_to_angular_ket to speed up lifetime calculations --- src/rydstate/angular/utils.py | 2 ++ src/rydstate/rydberg/rydberg_sqdt.py | 6 +++++- 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/src/rydstate/angular/utils.py b/src/rydstate/angular/utils.py index b577eed..10382c9 100644 --- a/src/rydstate/angular/utils.py +++ b/src/rydstate/angular/utils.py @@ -2,6 +2,7 @@ import contextlib import typing as t +from functools import lru_cache from typing import TYPE_CHECKING, Any, Literal import numpy as np @@ -82,6 +83,7 @@ def get_possible_quantum_number_values(s_1: float, s_2: float, s_tot: float | No return [float(s) for s in np.arange(abs(s_1 - s_2), s_1 + s_2 + 1, 1)] +@lru_cache(maxsize=1_000) def quantum_numbers_to_angular_ket( species: str | SpeciesObject, s_c: float | None = None, diff --git a/src/rydstate/rydberg/rydberg_sqdt.py b/src/rydstate/rydberg/rydberg_sqdt.py index 39e89e7..a23937b 100644 --- a/src/rydstate/rydberg/rydberg_sqdt.py +++ b/src/rydstate/rydberg/rydberg_sqdt.py @@ -90,7 +90,7 @@ def __init__( l_tot=l_tot, j_tot=j_tot, f_tot=f_tot, - m=m, + m=m, # type: ignore [arg-type] ) self.n = n @@ -468,6 +468,10 @@ def _get_transition_rates_au( mask = transition_rates_au != 0 relevant_states_masked = [ket for ket, is_relevant in zip(relevant_states, mask, strict=True) if is_relevant] transition_rates_au = transition_rates_au[mask] + + if np.any(transition_rates_au < 0): + raise RuntimeError("Got negative transition rates, which should not happen.") + return relevant_states_masked, transition_rates_au @overload From e97783bb0553762d2228cdb4c7d36a3bd9f38e11 Mon Sep 17 00:00:00 2001 From: johannes-moegerle Date: Tue, 10 Mar 2026 10:35:44 +0100 Subject: [PATCH 07/10] add property i_c_number --- src/rydstate/angular/angular_ket.py | 5 ++--- src/rydstate/basis/basis_sqdt.py | 8 ++++---- src/rydstate/rydberg/rydberg_sqdt.py | 2 +- src/rydstate/species/species_object.py | 5 +++++ tests/test_all_elements.py | 2 +- 5 files changed, 13 insertions(+), 9 deletions(-) diff --git a/src/rydstate/angular/angular_ket.py b/src/rydstate/angular/angular_ket.py index cf9380e..30f2b7f 100644 --- a/src/rydstate/angular/angular_ket.py +++ b/src/rydstate/angular/angular_ket.py @@ -96,10 +96,9 @@ def __init__( species = SpeciesObjectSQDT.from_name(species) # use i_c = 0 for species without defined nuclear spin (-> ignore hyperfine) - species_i_c = species.i_c if species.i_c is not None else 0 - if i_c is not None and i_c != species_i_c: + if i_c is not None and i_c != species.i_c_number: raise ValueError(f"Nuclear spin i_c={i_c} does not match the species {species} with i_c={species.i_c}.") - i_c = species_i_c + i_c = species.i_c_number s_c = 0.5 * (species.number_valence_electrons - 1) if i_c is None: raise ValueError("Nuclear spin i_c must be set or a species must be given.") diff --git a/src/rydstate/basis/basis_sqdt.py b/src/rydstate/basis/basis_sqdt.py index ca5c735..03fafd0 100644 --- a/src/rydstate/basis/basis_sqdt.py +++ b/src/rydstate/basis/basis_sqdt.py @@ -42,7 +42,7 @@ def __init__( species = SpeciesObjectSQDT.from_name(species) self.species = species - i_c = self.species.i_c if self.species.i_c is not None else 0 + i_c = self.species.i_c_number self.states = [] _s = 1 / 2 @@ -65,7 +65,7 @@ def __init__( species = SpeciesObjectSQDT.from_name(species) self.species = species - i_c = self.species.i_c if self.species.i_c is not None else 0 + i_c = self.species.i_c_number self.states = [] for _n in range(n[0], n[1] + 1): @@ -90,7 +90,7 @@ def __init__( species = SpeciesObjectSQDT.from_name(species) self.species = species - i_c = self.species.i_c if self.species.i_c is not None else 0 + i_c = self.species.i_c_number self.states = [] _j_c = 0.5 @@ -123,7 +123,7 @@ def __init__( species = SpeciesObjectSQDT.from_name(species) self.species = species - i_c = self.species.i_c if self.species.i_c is not None else 0 + i_c = self.species.i_c_number self.states = [] _j_c = 0.5 diff --git a/src/rydstate/rydberg/rydberg_sqdt.py b/src/rydstate/rydberg/rydberg_sqdt.py index a23937b..bc7b6ba 100644 --- a/src/rydstate/rydberg/rydberg_sqdt.py +++ b/src/rydstate/rydberg/rydberg_sqdt.py @@ -572,7 +572,7 @@ def __repr__(self) -> str: species, n, nu = self.species.name, self.n, self.nu l, j, f, m = self.l, self.j, self.f, self.m n_str = f", {n=}" if n is not None else "" - f_string = f", {f=}" if self.species.i_c not in (None, 0) else "" + f_string = f", {f=}" if self.species.i_c_number != 0 else "" return f"{self.__class__.__name__}({species}{n_str}, {nu=}, {l=}, {j=}{f_string}, {m=})" diff --git a/src/rydstate/species/species_object.py b/src/rydstate/species/species_object.py index 82568be..a01b84f 100644 --- a/src/rydstate/species/species_object.py +++ b/src/rydstate/species/species_object.py @@ -62,6 +62,11 @@ class SpeciesObject(ABC): defined in: Y. Fei et al., Chin. Phys. B 18, 4349 (2009), https://iopscience.iop.org/article/10.1088/1674-1056/18/10/025 """ + @property + def i_c_number(self) -> float: + """Return a numerical value for i_c, i.e. either i_c or 0 if i_c is None.""" + return self.i_c if self.i_c is not None else 0 + def __repr__(self) -> str: return f"{self.__class__.__name__}()" diff --git a/tests/test_all_elements.py b/tests/test_all_elements.py index ca4d7a0..a8031a3 100644 --- a/tests/test_all_elements.py +++ b/tests/test_all_elements.py @@ -11,7 +11,7 @@ @pytest.mark.parametrize("species_name", SpeciesObjectSQDT.get_available_species()) def test_sqdt_species(species_name: str) -> None: species = SpeciesObjectSQDT.from_name(species_name) - i_c = species.i_c if species.i_c is not None else 0 + i_c = species.i_c_number state: RydbergStateSQDT if species.number_valence_electrons == 1: From aaa93bb4e34861d5edb6ffcb7c108e8c3c2be52e Mon Sep 17 00:00:00 2001 From: johannes-moegerle Date: Tue, 10 Mar 2026 10:36:26 +0100 Subject: [PATCH 08/10] various species object improvements --- src/rydstate/species/species_object.py | 12 ++++++++++-- src/rydstate/species/sqdt/species_object_sqdt.py | 4 ---- 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/src/rydstate/species/species_object.py b/src/rydstate/species/species_object.py index a01b84f..e44e05a 100644 --- a/src/rydstate/species/species_object.py +++ b/src/rydstate/species/species_object.py @@ -4,16 +4,24 @@ import logging from abc import ABC from functools import cache, cached_property -from typing import TYPE_CHECKING, ClassVar, overload +from typing import TYPE_CHECKING, ClassVar, TypeVar, overload from rydstate.units import rydberg_constant, ureg if TYPE_CHECKING: - from typing_extensions import Self + from collections.abc import Callable + + from typing_extensions import ParamSpec, Self from rydstate.radial.model import PotentialType from rydstate.units import PintFloat + P = ParamSpec("P") + R = TypeVar("R") + + def cache(func: Callable[P, R]) -> Callable[P, R]: ... # type: ignore [misc] + + logger = logging.getLogger(__name__) diff --git a/src/rydstate/species/sqdt/species_object_sqdt.py b/src/rydstate/species/sqdt/species_object_sqdt.py index 75532e9..584c6af 100644 --- a/src/rydstate/species/sqdt/species_object_sqdt.py +++ b/src/rydstate/species/sqdt/species_object_sqdt.py @@ -71,10 +71,6 @@ def _setup_nist_energy_levels(self, file: Path) -> None: # noqa: C901, PLR0912 Args: file: Path to the NIST energy levels file. - n_max: Maximum principal quantum number for which to load the NIST energy levels. - For large quantum numbers, the NIST data is not accurate enough - (it does not even show fine structure splitting), - so we limit the maximum principal quantum number to 15 by default. """ if not file.exists(): From 0b47c9abd7e22a50bda90b0d3469a85ac891018c Mon Sep 17 00:00:00 2001 From: johannes-moegerle Date: Tue, 10 Mar 2026 10:36:42 +0100 Subject: [PATCH 09/10] add test lifetime --- tests/test_lifetimes.py | 58 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 58 insertions(+) create mode 100644 tests/test_lifetimes.py diff --git a/tests/test_lifetimes.py b/tests/test_lifetimes.py new file mode 100644 index 0000000..3331ba3 --- /dev/null +++ b/tests/test_lifetimes.py @@ -0,0 +1,58 @@ +import numpy as np +import pytest +from rydstate import RydbergStateSQDTAlkali +from rydstate.angular.angular_ket import AngularKetLS +from rydstate.rydberg.rydberg_sqdt import RydbergStateSQDT +from rydstate.species.sqdt.species_object_sqdt import SpeciesObjectSQDT + + +# Reference values from NIST Atomic Spectra Database (ASD), Einstein A coefficients: +# H 2p -> 1s: A = 6.2648e8 s^-1 +# H 3p -> all lower: A_tot = 1.8971e8 s^-1 (A(3p->1s)=1.6725e8, A(3p->2s)=2.245e7) +# H 3d -> 2p: A = 6.4651e7 s^-1 +# Source: https://physics.nist.gov/PhysRefData/ASD/lines_form.html +@pytest.mark.parametrize( + ("n", "l", "j", "expected_gamma"), + [ + (2, 1, 1.5, 6.2648e8), + (3, 1, 1.5, 1.8971e8), + (3, 2, 2.5, 6.4651e7), + ], +) +def test_hydrogen_textbook_lifetimes(n: int, l: int, j: float, expected_gamma: float) -> None: + """Test that calculated H lifetimes match NIST textbook values within 5%.""" + state = RydbergStateSQDTAlkali("H_textbook", n=n, l=l, j=j, m=0.5) + tau = state.get_lifetime(unit="s") + np.testing.assert_allclose(tau, 1 / expected_gamma, rtol=0.05) + + +def test_bbr_shortens_lifetime() -> None: + """Test that black body radiation at 300 K shortens the lifetime relative to T=0.""" + state = RydbergStateSQDTAlkali("Rb", n=30, l=0, j=0.5, m=0.5) + tau_0 = state.get_lifetime(unit="mus") + tau_300 = state.get_lifetime(300, temperature_unit="K", unit="mus") + assert tau_300 < tau_0 + + +@pytest.mark.parametrize("species_name", SpeciesObjectSQDT.get_available_species()) +def test_lifetime_n_scaling(species_name: str) -> None: + """Test that Rydberg state lifetimes scale as nu^3 (effective quantum number).""" + if species_name in ["Sr87", "Yb171", "Yb173"]: + pytest.skip("No quantum defect data available") + if species_name == "Yb174": + pytest.skip("Quantum defects not correct for low n states") + + n1, n2 = 30, 60 + species = SpeciesObjectSQDT.from_name(species_name) + s_tot = (species.number_valence_electrons / 2) % 1 + f = abs(s_tot - species.i_c_number) + angular = AngularKetLS(l_r=0, s_tot=s_tot, f_tot=f, m=f, species=species) + + state1 = RydbergStateSQDT.from_angular_ket(species, angular, n=n1) + state2 = RydbergStateSQDT.from_angular_ket(species, angular, n=n2) + tau1 = state1.get_lifetime(unit="mus") + tau2 = state2.get_lifetime(unit="mus") + nu1 = state1.nu + nu2 = state2.nu + expected_ratio = (nu2 / nu1) ** 3 + np.testing.assert_allclose(tau2 / tau1, expected_ratio, rtol=0.05) From ef105b4b8eb855c48558f9a95c0a2ea95928b46d Mon Sep 17 00:00:00 2001 From: johannes-moegerle Date: Tue, 10 Mar 2026 11:36:46 +0100 Subject: [PATCH 10/10] various small improvements --- src/rydstate/angular/angular_ket.py | 2 +- src/rydstate/basis/basis_sqdt.py | 30 ++++++++++--------- src/rydstate/rydberg/rydberg_sqdt.py | 16 +++++----- .../species/sqdt/species_object_sqdt.py | 4 +-- tests/test_lifetimes.py | 7 ++--- 5 files changed, 30 insertions(+), 29 deletions(-) diff --git a/src/rydstate/angular/angular_ket.py b/src/rydstate/angular/angular_ket.py index 30f2b7f..a0ca7ad 100644 --- a/src/rydstate/angular/angular_ket.py +++ b/src/rydstate/angular/angular_ket.py @@ -159,7 +159,7 @@ def __str__(self) -> str: def __eq__(self, other: object) -> bool: if not isinstance(other, AngularKetBase): - raise NotImplementedError(f"Cannot compare {self!r} with {other!r}.") + return NotImplemented if type(self) is not type(other): return False if self.m != other.m: diff --git a/src/rydstate/basis/basis_sqdt.py b/src/rydstate/basis/basis_sqdt.py index 03fafd0..77fbe85 100644 --- a/src/rydstate/basis/basis_sqdt.py +++ b/src/rydstate/basis/basis_sqdt.py @@ -97,13 +97,14 @@ def __init__( _s_r = 0.5 for _n in range(n[0], n[1] + 1): for _l_r in range(_n): - if self.species.is_allowed_shell(_n, _l_r, 0) != self.species.is_allowed_shell(_n, _l_r, 1): - logger.warning( - "For l=%d, n=%d one of the singlet/triplet states is not allowed. " - "In JJ coupling the state does not exist, thus skipping this shell", - *(_l_r, _n), - ) - if not all(self.species.is_allowed_shell(_n, _l_r, s_tot) for s_tot in [0, 1]): + allowed = [self.species.is_allowed_shell(_n, _l_r, s) for s in [0, 1]] + if not all(allowed): + if any(allowed): + logger.warning( + "For l=%d, n=%d one of the singlet/triplet states is not allowed. " + "In JJ coupling the state does not exist, thus skipping this shell", + *(_l_r, _n), + ) continue for _j_r in np.arange(abs(_l_r - _s_r), _l_r + _s_r + 1): for _j_tot in range(int(abs(_j_r - _j_c)), int(_j_r + _j_c + 1)): @@ -130,13 +131,14 @@ def __init__( _s_r = 0.5 for _n in range(n[0], n[1] + 1): for _l_r in range(_n): - if self.species.is_allowed_shell(_n, _l_r, 0) != self.species.is_allowed_shell(_n, _l_r, 1): - logger.warning( - "For l=%d, n=%d one of the singlet/triplet states is not allowed. " - "In FJ coupling the state does not exist, thus skipping this shell", - *(_l_r, _n), - ) - if not all(self.species.is_allowed_shell(_n, _l_r, s_tot) for s_tot in [0, 1]): + allowed = [self.species.is_allowed_shell(_n, _l_r, s) for s in [0, 1]] + if not all(allowed): + if any(allowed): + logger.warning( + "For l=%d, n=%d one of the singlet/triplet states is not allowed. " + "In FJ coupling the state does not exist, thus skipping this shell", + *(_l_r, _n), + ) continue for _j_r in np.arange(abs(_l_r - _s_r), _l_r + _s_r + 1): for _f_c in np.arange(abs(_j_c - i_c), _j_c + i_c + 1): diff --git a/src/rydstate/rydberg/rydberg_sqdt.py b/src/rydstate/rydberg/rydberg_sqdt.py index bc7b6ba..a6934ac 100644 --- a/src/rydstate/rydberg/rydberg_sqdt.py +++ b/src/rydstate/rydberg/rydberg_sqdt.py @@ -413,6 +413,8 @@ def _get_transition_rates_au( \Gamma^{spontaneous}_{self \to other} = \frac{4}{3} \frac{\alpha}{c^2} \omega^3 |\langle self || r^k_radial \hat{O}_{k_angular} || other \rangle|^2 + where :math:`\alpha = 1/c` in atomic units. + and .. math:: @@ -443,27 +445,25 @@ def _get_transition_rates_au( ) electric_dipole_moments_au = np.zeros(len(relevant_states)) for q in [-1, 0, 1]: - # the different entries are only at most once nonzero -> we can just add the arrays el_di_m = np.array( [s.calc_matrix_element(self, "electric_dipole", q, unit="a.u.") for s in relevant_states] ) - electric_dipole_moments_au += el_di_m + electric_dipole_moments_au += np.abs(el_di_m) ** 2 transition_rates_au = ( - (4 / 3) - * np.abs(electric_dipole_moments_au) ** 2 - * energy_differences_au**2 - / ureg.Quantity(1, "speed_of_light").to_base_units().magnitude ** 3 + (4 / 3) * electric_dipole_moments_au / ureg.Quantity(1, "speed_of_light").to_base_units().magnitude ** 3 ) if only_spontaneous: - transition_rates_au *= energy_differences_au + transition_rates_au *= energy_differences_au**3 else: assert temperature_au is not None, "Temperature must be given for black body transitions." if temperature_au == 0: transition_rates_au *= 0 else: # for numerical stability we use 1 / exprel(x) = x / (exp(x) - 1) - transition_rates_au *= temperature_au / exprel(energy_differences_au / temperature_au) + transition_rates_au *= ( + energy_differences_au**2 * temperature_au / exprel(energy_differences_au / temperature_au) + ) mask = transition_rates_au != 0 relevant_states_masked = [ket for ket, is_relevant in zip(relevant_states, mask, strict=True) if is_relevant] diff --git a/src/rydstate/species/sqdt/species_object_sqdt.py b/src/rydstate/species/sqdt/species_object_sqdt.py index 584c6af..2e9bc9a 100644 --- a/src/rydstate/species/sqdt/species_object_sqdt.py +++ b/src/rydstate/species/sqdt/species_object_sqdt.py @@ -164,11 +164,11 @@ def get_ionization_energy(self, unit: None = None) -> PintFloat: ... @overload def get_ionization_energy(self, unit: str) -> float: ... - def get_ionization_energy(self, unit: str | None = "hartree") -> PintFloat | float: + def get_ionization_energy(self, unit: str | None = None) -> PintFloat | float: """Return the ionization energy in the desired unit. Args: - unit: Desired unit for the ionization energy. Default is atomic units "hartree". + unit: Desired unit for the ionization energy. Default is None (returns a Pint quantity). Returns: Ionization energy in the desired unit. diff --git a/tests/test_lifetimes.py b/tests/test_lifetimes.py index 3331ba3..915db18 100644 --- a/tests/test_lifetimes.py +++ b/tests/test_lifetimes.py @@ -1,9 +1,8 @@ import numpy as np import pytest -from rydstate import RydbergStateSQDTAlkali -from rydstate.angular.angular_ket import AngularKetLS -from rydstate.rydberg.rydberg_sqdt import RydbergStateSQDT -from rydstate.species.sqdt.species_object_sqdt import SpeciesObjectSQDT +from rydstate import RydbergStateSQDT, RydbergStateSQDTAlkali +from rydstate.angular import AngularKetLS +from rydstate.species import SpeciesObjectSQDT # Reference values from NIST Atomic Spectra Database (ASD), Einstein A coefficients: