diff --git a/.github/actions/setup_env/action.yml b/.github/actions/setup_env/action.yml index 8c64f24dff4..a0e4abb5e35 100644 --- a/.github/actions/setup_env/action.yml +++ b/.github/actions/setup_env/action.yml @@ -18,7 +18,7 @@ runs: - name: Generate Cache Key run: | file_hash=$(cat conda-${{ inputs.os-label }}.lock | shasum -a 256 | cut -d' ' -f1) - echo "file_hash=$file_hash" >> "${GITHUB_OUTPUT}" + echo "file_hash=tardis-conda-env-${{ inputs.os-label }}-${file_hash}-v1" >> "${GITHUB_OUTPUT}" id: cache-environment-key shell: bash diff --git a/.github/actions/setup_lfs/action.yml b/.github/actions/setup_lfs/action.yml index 86a3b0464d4..d937ddfc316 100644 --- a/.github/actions/setup_lfs/action.yml +++ b/.github/actions/setup_lfs/action.yml @@ -1,12 +1,16 @@ name: "Setup LFS" -description: "Pull LFS repositories and caches them" +description: "Sets up Git LFS, retrieves LFS cache and fails if cache is not available" inputs: regression-data-repo: - description: "tardis regression data repository" + description: "Repository containing regression data (format: owner/repo)" required: false default: "tardis-sn/tardis-regression-data" + atom-data-sparse: + description: "If true, only downloads atom_data/kurucz_cd23_chianti_H_He.h5 instead of full regression data" + required: false + default: 'false' runs: using: "composite" @@ -16,37 +20,31 @@ runs: with: repository: ${{ inputs.regression-data-repo }} path: tardis-regression-data + sparse-checkout: ${{ inputs.atom-data-sparse == 'true' && 'atom_data/kurucz_cd23_chianti_H_He.h5' || '' }} + lfs: false - name: Create LFS file list - run: git lfs ls-files -l | cut -d' ' -f1 | sort > .lfs-assets-id + run: | + if [ "${{ inputs.atom-data-sparse }}" == "true" ]; then + echo "Using atom data sparse checkout" + echo "atom_data/kurucz_cd23_chianti_H_He.h5" > .lfs-files-list + else + echo "Using full repository checkout" + git lfs ls-files -l | cut -d' ' -f1 | sort > .lfs-files-list + fi working-directory: tardis-regression-data shell: bash - + - name: Restore LFS cache uses: actions/cache/restore@v4 id: lfs-cache-regression-data with: path: tardis-regression-data/.git/lfs - key: ${{ runner.os }}-lfs-${{ hashFiles('tardis-regression-data/.lfs-assets-id') }}-v1 - - - name: Git LFS Pull - run: git lfs pull - working-directory: tardis-regression-data - if: steps.lfs-cache-regression-data.outputs.cache-hit != 'true' - shell: bash + key: tardis-regression-${{ inputs.atom-data-sparse == 'true' && 'atom-data-sparse' || 'full-data' }}-${{ hashFiles('tardis-regression-data/.lfs-files-list') }}-${{ inputs.regression-data-repo }}-v1 + fail-on-cache-miss: true - name: Git LFS Checkout run: git lfs checkout working-directory: tardis-regression-data if: steps.lfs-cache-regression-data.outputs.cache-hit == 'true' shell: bash - - - name: Save LFS cache if not found - # uses fake ternary - # for reference: https://github.com/orgs/community/discussions/26738#discussioncomment-3253176 - if: ${{ steps.lfs-cache-regression-data.outputs.cache-hit != 'true' && !contains(github.ref, 'merge') && always() || false }} - uses: actions/cache/save@v4 - id: lfs-cache-regression-data-save - with: - path: tardis-regression-data/.git/lfs - key: ${{ runner.os }}-lfs-${{ hashFiles('tardis-regression-data/.lfs-assets-id') }}-v1 diff --git a/.github/workflows/benchmarks.yml b/.github/workflows/benchmarks.yml index e1e84afb6ca..db9f730debb 100644 --- a/.github/workflows/benchmarks.yml +++ b/.github/workflows/benchmarks.yml @@ -29,6 +29,12 @@ defaults: shell: bash -l {0} jobs: + test-cache: + uses: ./.github/workflows/lfs-cache.yml + with: + atom-data-sparse: false + regression-data-repo: tardis-sn/tardis-regression-data + build: if: github.repository_owner == 'tardis-sn' && (github.event_name == 'push' || @@ -37,6 +43,7 @@ jobs: (github.event_name == 'pull_request_target' && contains(github.event.pull_request.labels.*.name, 'benchmarks'))) runs-on: ubuntu-latest + needs: [test-cache] steps: - uses: actions/checkout@v4 if: github.event_name != 'pull_request_target' @@ -54,13 +61,10 @@ jobs: run: git fetch origin master:master if: github.event_name == 'pull_request_target' - - uses: actions/checkout@v4 + - name: Setup LFS + uses: ./.github/actions/setup_lfs with: - repository: tardis-sn/tardis-regression-data - path: tardis-regression-data - lfs: true - sparse-checkout: | - atom_data/kurucz_cd23_chianti_H_He.h5 + atom-data-sparse: true - name: Setup Mamba uses: mamba-org/setup-micromamba@v1 diff --git a/.github/workflows/build-docs.yml b/.github/workflows/build-docs.yml index 22a9d1369c6..b9a928d4aef 100644 --- a/.github/workflows/build-docs.yml +++ b/.github/workflows/build-docs.yml @@ -36,6 +36,12 @@ defaults: shell: bash -l {0} jobs: + test-cache: + uses: ./.github/workflows/lfs-cache.yml + with: + atom-data-sparse: true + regression-data-repo: tardis-sn/tardis-regression-data + check-for-changes: runs-on: ubuntu-latest if: ${{ !github.event.pull_request.draft }} @@ -77,7 +83,7 @@ jobs: build-docs: runs-on: ubuntu-latest - needs: check-for-changes + needs: [test-cache, check-for-changes] if: needs.check-for-changes.outputs.trigger-check-outcome == 'success' || needs.check-for-changes.outputs.docs-check-outcome == 'success' steps: - uses: actions/checkout@v4 @@ -90,13 +96,10 @@ jobs: ref: ${{ github.event.pull_request.head.sha }} if: github.event_name == 'pull_request_target' - - uses: actions/checkout@v4 + - name: Setup LFS + uses: ./.github/actions/setup_lfs with: - repository: tardis-sn/tardis-regression-data - path: tardis-regression-data - lfs: true - sparse-checkout: | - atom_data/kurucz_cd23_chianti_H_He.h5 + atom-data-sparse: true - name: Setup environment uses: ./.github/actions/setup_env diff --git a/.github/workflows/lfs-cache.yml b/.github/workflows/lfs-cache.yml new file mode 100644 index 00000000000..500fc14237d --- /dev/null +++ b/.github/workflows/lfs-cache.yml @@ -0,0 +1,77 @@ +name: Save LFS Cache + +on: + workflow_call: + inputs: + atom-data-sparse: + description: "If true, only downloads atom_data/kurucz_cd23_chianti_H_He.h5" + required: false + default: false + type: boolean + regression-data-repo: + description: "Repository containing regression data (format: owner/repo)" + required: false + default: "tardis-sn/tardis-regression-data" + type: string + +defaults: + run: + shell: bash -l {0} + +concurrency: + # Only one workflow can run at a time + # the workflow group is a unique identifier and contains the workflow name, pull request number, atom data sparse, and regression data repo + group: ${{ github.workflow }}-${{ github.event.pull_request.number || github.ref }}-${{ inputs.atom-data-sparse == 'true' && 'atom-data-sparse' || 'full-data' }}-${{ inputs.regression-data-repo }} + cancel-in-progress: true + + +jobs: + lfs-cache: + if: github.repository_owner == 'tardis-sn' + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + with: + repository: ${{ inputs.regression-data-repo }} + path: tardis-regression-data + sparse-checkout: ${{ inputs.atom-data-sparse == 'true' && 'atom_data/kurucz_cd23_chianti_H_He.h5' || '' }} + + - name: Create LFS file list + run: | + if [ "${{ inputs.atom-data-sparse }}" == "true" ]; then + echo "Using atom data sparse checkout" + echo "atom_data/kurucz_cd23_chianti_H_He.h5" > .lfs-files-list + else + echo "Using full repository checkout" + git lfs ls-files -l | cut -d' ' -f1 | sort > .lfs-files-list + fi + working-directory: tardis-regression-data + + + - name: Test cache availability + uses: actions/cache/restore@v4 + id: test-lfs-cache-regression-data + with: + path: tardis-regression-data/.git/lfs + key: tardis-regression-${{ inputs.atom-data-sparse == 'true' && 'atom-data-sparse' || 'full-data' }}-${{ hashFiles('tardis-regression-data/.lfs-files-list') }}-${{ inputs.regression-data-repo }}-v1 + lookup-only: true + + - name: Git LFS Pull Atom Data + run: git lfs pull --include=atom_data/kurucz_cd23_chianti_H_He.h5 + if: ${{ inputs.atom-data-sparse == true && steps.test-lfs-cache-regression-data.outputs.cache-hit != 'true' }} + working-directory: tardis-regression-data + + - name: Git LFS Pull Full Data + run: git lfs pull + if: ${{ inputs.atom-data-sparse == false && steps.test-lfs-cache-regression-data.outputs.cache-hit != 'true' }} + working-directory: tardis-regression-data + + - name: Git LFS Checkout + run: git lfs checkout + working-directory: tardis-regression-data + + - name: Save LFS cache if not found + uses: actions/cache/save@v4 + with: + path: tardis-regression-data/.git/lfs + key: tardis-regression-${{ inputs.atom-data-sparse == true && 'atom-data-sparse' || 'full-data' }}-${{ hashFiles('tardis-regression-data/.lfs-files-list') }}-${{ inputs.regression-data-repo }}-v1 diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 7513c56e9d5..b0b4353a78f 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -38,9 +38,16 @@ concurrency: cancel-in-progress: true jobs: + test-cache: + uses: ./.github/workflows/lfs-cache.yml + with: + atom-data-sparse: false + regression-data-repo: tardis-sn/tardis-regression-data + tests: name: ${{ matrix.continuum }} continuum ${{ matrix.os }} ${{ inputs.pip_git && 'pip tests enabled' || '' }} if: github.repository_owner == 'tardis-sn' + needs: [test-cache] runs-on: ${{ matrix.os }} strategy: fail-fast: false diff --git a/CHANGELOG.md b/CHANGELOG.md index 710ccec0cc1..84ff205c398 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,11 @@ ## Changelog -### release-2025.01.19 (2025/01/18 20:07) +### release-2025.01.26 (2025/01/25 20:04) +- [2932](https://github.com/tardis-sn/tardis/pull/2932) opacity_state_to_numba -> to_numba and moved to as a method of class OpacityState (2932) (@Sonu0305) +- [2897](https://github.com/tardis-sn/tardis/pull/2897) Ionization rates (2897) (@andrewfullard) +- [2938](https://github.com/tardis-sn/tardis/pull/2938) Fixes spelling errors in the codebase (2938) (@Sonu0305) +- [2952](https://github.com/tardis-sn/tardis/pull/2952) Post-release 2025.01.19 (2952) (@tardis-bot) +### release-2025.01.19 (2025/01/14 16:54) - [2800](https://github.com/tardis-sn/tardis/pull/2800) V inner formal integral (2800) (@Rodot-) - [2946](https://github.com/tardis-sn/tardis/pull/2946) Fixes indentation for tracking key to be nested under montecarlo key in YAML (2946) (@Sonu0305) - [2907](https://github.com/tardis-sn/tardis/pull/2907) Add missing init.py file in opacitites module (2907) (@KasukabeDefenceForce) diff --git a/CITATION.cff b/CITATION.cff index 2a2b4e1bab6..10aeeb89116 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -3,8 +3,8 @@ cff-version: 1.0.3 message: If you use this software, please cite it using these metadata. # FIXME title as repository name might not be the best name, please make human readable -title: 'tardis-sn/tardis: TARDIS v2025.01.12' -doi: 10.5281/zenodo.14633332 +title: 'tardis-sn/tardis: TARDIS v2025.01.26' +doi: 10.5281/zenodo.14740815 # FIXME splitting of full names is error prone, please check if given/family name are correct authors: - given-names: Wolfgang @@ -347,7 +347,7 @@ authors: - given-names: Atharwa family-names: Kharkar affiliation: -version: release-2025.01.12 -date-released: 2025-01-12 +version: release-2025.01.26 +date-released: 2025-01-26 repository-code: https://github.com/tardis-sn/tardis license: cc-by-4.0 diff --git a/README.rst b/README.rst index 18e8c49bc27..6ea4e58dbd6 100644 --- a/README.rst +++ b/README.rst @@ -110,14 +110,14 @@ The following BibTeX entries are needed for the references: adsnote = {Provided by the SAO/NASA Astrophysics Data System} } -.. |CITATION| replace:: kerzendorf_2025_14633332 +.. |CITATION| replace:: kerzendorf_2025_14740815 -.. |DOI_BADGE| image:: https://img.shields.io/badge/DOI-10.5281/zenodo.14633332-blue - :target: https://doi.org/10.5281/zenodo.14633332 +.. |DOI_BADGE| image:: https://img.shields.io/badge/DOI-10.5281/zenodo.14740815-blue + :target: https://doi.org/10.5281/zenodo.14740815 .. code-block:: bibtex - @software{kerzendorf_2025_14633332, + @software{kerzendorf_2025_14740815, author = {Kerzendorf, Wolfgang and Sim, Stuart and Vogl, Christian and @@ -225,13 +225,13 @@ The following BibTeX entries are needed for the references: Nayak U, Ashwin and Kumar, Atul and Kharkar, Atharwa}, - title = {tardis-sn/tardis: TARDIS v2025.01.12}, + title = {tardis-sn/tardis: TARDIS v2025.01.26}, month = jan, year = 2025, publisher = {Zenodo}, - version = {release-2025.01.12}, - doi = {10.5281/zenodo.14633332}, - url = {https://doi.org/10.5281/zenodo.14633332}, + version = {release-2025.01.26}, + doi = {10.5281/zenodo.14740815}, + url = {https://doi.org/10.5281/zenodo.14740815}, } ******* diff --git a/docs/contributing/development/continuous_integration.rst b/docs/contributing/development/continuous_integration.rst index 45db365ccf1..82aec6f8964 100644 --- a/docs/contributing/development/continuous_integration.rst +++ b/docs/contributing/development/continuous_integration.rst @@ -27,6 +27,30 @@ TARDIS Pipelines Brief description of pipelines already implemented on TARDIS +Cache Keys in TARDIS CI +----------------------- + +TARDIS uses specific cache key formats to efficiently store and retrieve data during CI runs: + +1. **Regression Data Cache Keys** + - Format: ``tardis-regression---v1`` + - Examples: + - ``tardis-regression-atom-data-sparse--v1`` - For atomic data cache + - ``tardis-regression-full-data--v1`` - For full TARDIS regression data cache + - Used in: ``setup_lfs`` action + +2. **Environment Cache Keys** + - Format: ``tardis-conda-env---v1`` + - Examples: + - ``tardis-conda-env-linux--v1`` - For Linux conda environment + - ``tardis-conda-env-macos--v1`` - For macOS conda environment + - Used in: ``setup_env`` action + +.. warning:: + - The version suffix (-v1) allows for future cache invalidation if needed. + - Sometimes the cache might not be saved due to race conditions between parallel jobs. Please check workflow runs when testing new regression data for cache misses to avoid consuming LFS quota. + + Streamlined Steps for TARDIS Pipelines ======================================== diff --git a/docs/contributing/development/running_tests.rst b/docs/contributing/development/running_tests.rst index a11ac2eca7e..c1150e51340 100644 --- a/docs/contributing/development/running_tests.rst +++ b/docs/contributing/development/running_tests.rst @@ -62,7 +62,7 @@ Or, to run tests for a particular file or directory To prevent leaking LFS quota, tests have been disabled on forks. If, by any chance, you need to run tests on your fork, make sure to run the tests workflow on master branch first. The LFS cache generated in the master branch should be available in all child branches. - You can check if cache was generated by looking in the ``Restore LFS Cache`` step of the workflow run. + You can check if cache was generated by looking in the ``Setup LFS`` step of the workflow run. Cache can also be found under the "Management" Section under "Actions" tab. Generating Plasma Reference diff --git a/docs/physics/plasma/equilibrium/tardis_solver_cmfgen.ipynb b/docs/physics/plasma/equilibrium/tardis_solver_cmfgen.ipynb new file mode 100644 index 00000000000..e990cb0a9a4 --- /dev/null +++ b/docs/physics/plasma/equilibrium/tardis_solver_cmfgen.ipynb @@ -0,0 +1,381 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Rates using CMFGEN data" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/afullard/tardis/tardis/__init__.py:20: UserWarning: Astropy is already imported externally. Astropy should be imported after TARDIS.\n", + " warnings.warn(\n" + ] + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "a9004b0713d345fc84b58a8cafbfde62", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Iterations: 0/? [00:00" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiUAAAHACAYAAAB9DBhHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABaB0lEQVR4nO3deZhU9Z3v8fepU0tX7wv0xq4isogiGMWIoiQYTBycmIlJTDRmmTBRoyFeE0ye7DM4GeMQx4hxRI0xLjcX9ZpIFOZGwA0jCAEVCWpDN9BN002vtS/n/lHdbTc00NVU1anu/ryepx66Tp2q861jzdQnvzq/39ewLMtCRERExGYOuwsQERERAYUSERERyRIKJSIiIpIVFEpEREQkKyiUiIiISFZQKBEREZGsoFAiIiIiWUGhRERERLKCQomIiIhkBYUSERERyQpDKpRs3LiRK664gurqagzD4Jlnnkn6NSzL4s477+T000/H4/Ewbtw4/u3f/i31xYqIiEhSnHYXkAyfz8dZZ53F9ddfz1VXXTWo17j55ptZu3Ytd955J2eeeSZtbW00NTWluFIRERFJljFUG/IZhsHTTz/NlVde2bMtHA7zgx/8gN///ve0trYyY8YM/v3f/5358+cDsHPnTmbOnMlbb73FlClT7ClcRERE+jWkfr45keuvv55XXnmFJ554gu3bt/NP//RPfOITn2D37t0A/PGPf+SUU07hT3/6E5MmTWLixIl87Wtf4/DhwzZXLiIiIsMmlLz//vs8/vjj/OEPf2DevHmceuqp3HrrrVx44YU89NBDAHzwwQfs3buXP/zhDzzyyCM8/PDDbNmyhc985jM2Vy8iIiJD6pqS43nzzTexLIvTTz+9z/ZQKERZWRkA8XicUCjEI4880rPfqlWrmD17Nrt27dJPOiIiIjYaNqEkHo9jmiZbtmzBNM0+j+Xn5wNQVVWF0+nsE1ymTp0KQG1trUKJiIiIjYZNKJk1axaxWIzGxkbmzZvX7z4f/ehHiUajvP/++5x66qkA/P3vfwdgwoQJGatVREREjjakZt90dnby3nvvAYkQctddd3HJJZdQWlrK+PHj+eIXv8grr7zCL3/5S2bNmkVTUxN/+ctfOPPMM7n88suJx+Oce+655Ofns2LFCuLxODfccAOFhYWsXbvW5ncnIiIysg2pULJ+/XouueSSo7Zfd911PPzww0QiEX7+85/zyCOPsH//fsrKypg7dy4/+clPOPPMMwE4cOAAN910E2vXriUvL49Fixbxy1/+ktLS0ky/HREREellSIUSERERGb6GzZRgERERGdoUSkRERCQrDInZN/F4nAMHDlBQUIBhGHaXIyIiIgNgWRYdHR1UV1fjcJx4HGRIhJIDBw4wbtw4u8sQERGRQairq2Ps2LEn3G9IhJKCggIg8aYKCwttrkZEREQGor29nXHjxvV8j5/IkAgl3T/ZFBYWKpSIiIgMMQO99EIXuoqIiEhWUCgRERGRrKBQIiIiIllhSFxTIiIiMlTFYjEikYjdZaSFy+XCNM2UvZ5CiYiISBpYlkVDQwOtra12l5JWxcXFVFZWpmQdMYUSERGRNOgOJOXl5eTm5g67xT8ty8Lv99PY2AhAVVXVSb+mQomIiEiKxWKxnkBSVlZmdzlp4/V6AWhsbKS8vPykf8rRha4iIiIp1n0NSW5urs2VpF/3e0zFdTMKJSIiImky3H6y6U8q36NCiYiIiGQFhRIRERHJCgolIiIi0se9997LpEmTyMnJYfbs2bz00ksZOa5CiU3icYtdDR3E45bdpYiIiPR48sknueWWW/j+97/P1q1bmTdvHosWLaK2tjbtx1YoscmDr9Rw2YqN/G7TXrtLERER6XHXXXfx1a9+la997WtMnTqVFStWMG7cOFauXJn2Y2udEpv8bV8bAG/WtnDdBRPtLUZERNLOsiwCkVjGj+t1mQOeIRMOh9myZQvf+973+mxfuHAhr776ajrK60OhxCYHDncy2djHnkMFdpciIiIZEIjEmPbDFzJ+3Hd+ehm57oF93Tc1NRGLxaioqOizvaKigoaGhnSU14d+vrHJhc3/h3We2ziv+WksS9eViIhI9jhyZMWyrIysuaKREhsEIzFOjewCE6bH3uWwL0xZvsfuskREJI28LpN3fnqZLccdqFGjRmGa5lGjIo2NjUeNnqSDQokNDrQGqDQO8bbbxfjgAWqafAolIiLDnGEYA/4ZxS5ut5vZs2ezbt06/vEf/7Fn+7p161i8eHHaj6+fb2ywryXA68XtfG5MFVuK2qk51Gl3SSIiIgAsXbqUBx54gAcffJCdO3fy7W9/m9raWpYsWZL2Y2d3ZBum6ptbqfMkrsDelWMwsaEOGG9vUSIiIsDVV19Nc3MzP/3pT6mvr2fGjBmsWbOGCRMmpP3YCiU2aD9YQwNOpu+Js2+0i8qDfwc+andZIiIiAHzzm9/km9/8ZsaPq59vbBBp3su0bSY/ejzO6dtNjMPv212SiIiI7RRKbBBrrWF0U2Jq1YQGcAR2a1qwiIiMeAolNgiH9lDW7qGhfDaVLW6crgMcbA/ZXZaIiIitFEoyLBSNYdFALHcB70z7CrHcS7BczdQ0+ewuTURExFYKJRlW3xrEYRwm5q4CwJ87jlgoyN6mdpsrExERsZdCSYbtawngDAUIekoB8OdWEOw0aT7wgc2ViYiI2EuhJMMONLdh+aMEcxKhJOAdTbzNReTg322uTERExF4KJRnWdnAPsU4vEXeiO7DlcJLbUY7R+p7NlYmIiNhLoSTDwk17cARG9dlW6ivHCr9HLK5pwSIiMnIplGRYvG0P3iNCSX6kAtPdwIHWgE1ViYiI2E+hJMMiwRrywiV9N5oVxI02TQsWERHbbdy4kSuuuILq6moMw+CZZ57J2LEVSjIoHI0Tox5PLHGRa443saqrP7eCuC9M7aEWO8sTERHB5/Nx1llncc8992T82GrIl0H1bQFMRws4ygAYN6WY3dta8OdWEOpw0rb/78Dp9hYpIiIj2qJFi1i0aJEtx1YoyaB9LQHMkL9nOvD4syvZva2FiCsfs62Q2KHdNlcoIiJpY1kQ8Wf+uK5cMIzMH3cQFEoy6EBzO0Yg0hNKSqvy8LojBMIuijorCbTtsrlCERFJm4gf/q0688e9/QC48zJ/3EFI6pqSlStXMnPmTAoLCyksLGTu3Ln8+c9/Pu5zNmzYwOzZs8nJyeGUU07hvvvuO6mCh7K2g3uJdnoIe4oBKCjLobjEBUBhsJxYfA+RWNzGCkVEROyT1EjJ2LFjueOOOzjttNMA+O1vf8vixYvZunUr06dPP2r/mpoaLr/8cr7+9a/z6KOP8sorr/DNb36T0aNHc9VVV6XmHQwhoaY9GL5R4AEHUXLyXJSOLaT+YDseqxKHaxt1h/2cMjrf7lJFRCTVXLmJUQs7jjtEJBVKrrjiij73//Vf/5WVK1eyadOmfkPJfffdx/jx41mxYgUAU6dOZfPmzdx5550jMpTEW/eQEywj7oFcTxTDMCg7ZTRsaSfkqcCIdLCn2adQIiIyHBnGkPkZxS6DvqYkFovxhz/8AZ/Px9y5c/vd57XXXmPhwoV9tl122WWsWrWKSCSCy+Xq93mhUIhQKNRzv719eHTQjQT3UBwuoR0oLEyc+pIxieXm/bkVRDvj7Gs4CGdU2FiliIiMZJ2dnbz33oetT2pqati2bRulpaWMHz8+rcdOep2SHTt2kJ+fj8fjYcmSJTz99NNMmzat330bGhqoqOj7BVtRUUE0GqWpqemYx1i+fDlFRUU9t3HjxiVbZtaJxOJY1gE80cRFrgWjEsNpJZWJ1Bz0lhFry6HzgC52FRER+2zevJlZs2Yxa9YsAJYuXcqsWbP44Q9/mPZjJz1SMmXKFLZt20ZrayurV6/muuuuY8OGDccMJsYR05Asy+p3e2/Lli1j6dKlPffb29uHfDBpaAviMA/jMBJrlBSNKQYgt8iNaUSJ4cTTWU68SdOCRUTEPvPnz+/5rs60pEOJ2+3uudB1zpw5vPHGG/zqV7/iN7/5zVH7VlZW0tDQ0GdbY2MjTqeTsrKyYx7D4/Hg8XiSLS2r1bX4McM+Ql3TgYvGJ96/YRgU5Vsc7oBiXyUtnRopERGRkemkl5m3LKvP9R+9zZ07l3Xr1vXZtnbtWubMmXPM60mGq/3NHdBrjZLCUR9eDV1Skfg7L1pOjFqCkZgtNYqIiNgpqVBy++2389JLL7Fnzx527NjB97//fdavX88111wDJH52ufbaa3v2X7JkCXv37mXp0qXs3LmTBx98kFWrVnHrrbem9l0MAa2NdUR8LkLda5SU5vQ8Vjqpa9TIUYHD0UTtYRtW/BMREbFZUj/fHDx4kC996UvU19dTVFTEzJkzef755/n4xz8OQH19PbW1tT37T5o0iTVr1vDtb3+bX//611RXV3P33XePyOnAoUN7MHylWC4Tgxi5he6ex0onlgENBHIrcfgDfNDYyekVBfYVKyIiYoOkQsmqVauO+/jDDz981LaLL76YN998M6mihiOrdQ+eQBm4wOuKYDg+vNC3ewaOL7eCSIfBwYY6OLPKrlJFRERscdLXlMjARIM15IcSP9MUFJp9Hisq9wJxYk4vRnspgXpd7CoiIiOPQkkGRGNxYvED5ERLACgs67vkr9NlkuuOAJDrqyTarFAiIiIjj0JJBtS3BXGYzZhG13TgqsKj9inpasxXFCwnHPx7RusTERHJBgolGbCvJYAZ8RH2dC2cNmHUUfuUjE0EFU+8gjj78YWiGa1RRETEbgolGbC/xYcRDBP0dK1RMvrojo1lkysBCHsqMGMt1DT5MlqjiIiI3RRKMqDlYC3RTpNgTuKakvxea5R0Kx2TGCnx51bg6Ayzp2l4NCEUEZGhY/ny5Zx77rkUFBRQXl7OlVdeya5dmbvOUaEkA0JNe7B8JVgOFxAnv/joJfSLu1Z1DeaUEmvPoXnf+xmuUkRERroNGzZwww03sGnTJtatW0c0GmXhwoX4fJkZvU+6940kL96yN7FGSQHkOMM4zKOzoLfAhdMIE8WNq7OC4MG/A/MzXquIiIxczz//fJ/7Dz30EOXl5WzZsoWLLroo7cdXKMmAqL+G0lAZvgIoyDf73ccwDAq7GvPl+yvxtezMcJUiIpJOlmURiAYyflyv04thGCfesR9tbW0AlJaWprKkY1IoSbNoLI4V3483OhkfUFh29PUk3Uor8jjcESUvUsHhyHuZK1JERNIuEA1w3mPnZfy4r3/hdXJdR0+wOBHLsli6dCkXXnghM2bMSENlR9M1JWl2sCOEw9mE00pMBy6sPHqNkm7djfkMRwWWVU+bP5KRGkVERI504403sn37dh5//PGMHVMjJWm277AfR6yTsLtr4bTxZcfct+zU0bDuIIHcClyBDmqafZydW5yhSkVEJJ28Ti+vf+F1W46brJtuuolnn32WjRs3Mnbs2DRU1T+FkjTbd9iH4Q8TzOlao6Qi/5j79m7MR3uM2oMtnD2uOBNliohImhmGMaifUTLJsixuuukmnn76adavX8+kSZMyenz9fJNmhxv3E/E7ekJJQT9rlHQrHO3FIE7c9GB1lnF4v3rgiIhI5txwww08+uijPPbYYxQUFNDQ0EBDQwOBQGYu0FUoSbNQUw34CombibVJCkqOHUpM00GuOwRAjq+S8EGFEhERyZyVK1fS1tbG/Pnzqaqq6rk9+eSTGTm+fr5Js3jLXlz+UZAPHjOE6Tp+DiwpdeNrgIJQJU3t72SoShERkcTPN3bSSEmaxf17KAglLm7NzzvxPPHSsUUAuOMVBCM1tn9AREREMkWhJI1icYtobB/eSNdFrsf56aZb2eQKACLucpzRQzT7wmmtUUREJFsolKTRwfYgTrMJl9UVSioLTvic0q7ZNv7cSlx+n7oFi4jIiKFQkkb7WgIYVjtR14nXKOnW3ZgvlFOC1e5mX31DWmsUERHJFgolabS/xYcR6LVGyQBGSnLyXLiMxAwcs7OaDk0LFhGREUKhJI2aDx4g6jcGtEZJb4X5cQBy/RUED6kxn4iIjAwKJWkUOFRD3JdP1Jn4SWagoaS0MrHqa260Al+nQomIiIwMCiVpFGvtWqMEcDnCuDzmgJ5XNinxHMNRSShaSzyuacEiIjL8KZSkkdW5h/xA1xoluQMPFqWnlQMQ9JbjCRzmYEcwLfWJiIhkE4WSNInHLeKxOnK71igpKPYM+LmlVYnGfH5vOWZnmJpDnWmpUUREJJsolKRJY0cIh+swrviJuwMfqaDMi2HFiJtujI4yGg7UpatMERGRPlauXMnMmTMpLCyksLCQuXPn8uc//zkjx1YoSZN9LX4Mq5V49xol4068Rkk3h8Mg1534ycblq6Jjvy52FRGRzBg7dix33HEHmzdvZvPmzVx66aUsXryYt99+O+3HVkO+NNl32I8ZDPVMBy6qKkzq+SVlicZ8eaEKOg6rMZ+IyFBnWRZWIJDx4xpeL4Zx4t5r3a644oo+9//1X/+VlStXsmnTJqZPn57q8vpQKEmTpkP1hP0GQU/XNSVl3qSeXza2mH0NPtxWBb7AjnSUKCIiGWQFAuw6Z3bGjzvlzS0YubmDem4sFuMPf/gDPp+PuXPnpriyo+nnmzQJHNqD1ekl4k6s4lpQNrA1SrqNOr0SgKi7Aiu4n5imBYuISIbs2LGD/Px8PB4PS5Ys4emnn2batGlpP65GStIk3rIX0zcKcsFpRPB4kzvVJeNLAPB7K/B0tnGgNcC40sElXRERsZ/h9TLlzS22HDdZU6ZMYdu2bbS2trJ69Wquu+46NmzYkPZgolCSJlbnHvKDZYRzIc8bT/r5JV2N+cKeIpwdTj5obFcoEREZwgzDGPTPKJnmdrs57bTTAJgzZw5vvPEGv/rVr/jNb36T1uPq55s0iMctYtE68sKJ0Y6CInfSr+H2OnEbiQuiDF81TfveS2mNIiIiA2VZFqFQKO3H0UhJGhzqDGE6m3DHpwJQWJE3qNcpyI/R3AGeQCWd9TuBC1NYpYiIyNFuv/12Fi1axLhx4+jo6OCJJ55g/fr1PP/882k/tkJJGuxrCWAYrVhm13TgMaWDep2yygKaOyy80QqaW99KZYkiIiL9OnjwIF/60peor6+nqKiImTNn8vzzz/Pxj3887cdWKEmDfYd9GMFgzxolhdVFg3qdskmjYXcjDqMCf+D/pbJEERGRfq1atcq2Y+uakjQ4dKiRaMD4MJSMSv7KZ4CyrmnBAW8Fpu8g4WjyF8yKiIgMFQolaeA/VEOs00PYUwwkv0ZJt5KuxnwB72g8HQHqWvypKlFERCTrKJSkQfzwXpz+RK8bB1Fy8lyDep2CkhwcVgTL4cTRUcLegy2pLFNERCSrJBVKli9fzrnnnktBQQHl5eVceeWV7Nq167jPWb9+fWJu9hG3d99996QKz2qde8ntCiV5ObGkeg70ZjgMvF2N+Zz+Kg7XDeNzJiIiI15SoWTDhg3ccMMNbNq0iXXr1hGNRlm4cCE+n++Ez921axf19fU9t8mTJw+66GxmWRbxaG2vNUoGN0rSraQ08fyccCXtB9PfoVFERMQuSc2+OXKO8kMPPUR5eTlbtmzhoosuOu5zy8vLKS4uTrrAoeZQZwiH8xA5sdNoAwpHndzqfaPGl7LvoB9PvJyGdoUSEREZvk7qmpK2tjYASktPvA7HrFmzqKqqYsGCBbz44osnc9istq8lAI5WcHStUTK25KReb9TpVQBE3ZWEfDUnW56IiEjWGvQ6JZZlsXTpUi688EJmzJhxzP2qqqq4//77mT17NqFQiN/97ncsWLCA9evXH3N0JRQK9VnOtr29fbBlZty+lgBmKECoezrwmJMLJaUTS4H38edW4G5rIhiJkeMyU1CpiIhIdhl0KLnxxhvZvn07L7/88nH3mzJlClOmTOm5P3fuXOrq6rjzzjuPGUqWL1/OT37yk8GWZqvGQweJ+CGQk7jQdbBrlHQr7mrMF3Hl4+50sLfZz5TKgpOuU0REJNsM6uebm266iWeffZYXX3yRsWPHJv38888/n927dx/z8WXLltHW1tZzq6urG0yZtvA37iHW6SbsSaziOtg1Srq53CZuOgEwfRXU1TecdI0iIiIDsXz5cgzD4JZbbsnI8ZIaKbEsi5tuuomnn36a9evXM2nSpEEddOvWrVRVVR3zcY/Hg8fjGdRr2y1+eC8OXymWx8QgRm5B8h2Cj5SfH+NwJzgDVbTU7oRZw3PmkoiIZI833niD+++/n5kzZ2bsmEmNlNxwww08+uijPPbYYxQUFNDQ0EBDQwOBQKBnn2XLlnHttdf23F+xYgXPPPMMu3fv5u2332bZsmWsXr2aG2+8MXXvIpt07CE3kPjpJtcdxXAMbo2S3kZVJX6uyYlW0nZo+0m/noiIyPF0dnZyzTXX8N///d+UlJzctZHJSGqkZOXKlQDMnz+/z/aHHnqIL3/5ywDU19dTW1vb81g4HObWW29l//79eL1epk+fznPPPcfll19+cpVnIcuyiEXrKAiV0AYUFKam3+HoUyr4++5DmEY5rR1/TclriohIZlmWRTSc+R5mTrcj6UU8b7jhBj75yU/ysY99jJ///OdpquxoSf98cyIPP/xwn/u33XYbt912W1JFDVVNnWGcZiOe2HgACk7yItduZVOq4IVDBLwV0FZ74ieIiEjWiYbj3H/zhowf959/dTEuz8BnbT7xxBO8+eabvPHGG2msqn+p+Z/yAsD+1gCG2YrDSPx8U1RdnJLXLa1O/HwT9I7C3d5BZyhKvkf/6UREJLXq6uq4+eabWbt2LTk5JzdRYzD0zZZC+1r8EP5wjZKicSdeVG4gcovcOKwQccODq7OYPYc6mTG2OCWvLSIimeF0O/jnX11sy3EHasuWLTQ2NjJ79uyebbFYjI0bN3LPPfcQCoUwzfStlaVQkkIHDzURC1gEU7RGSTfDMPC6/PiiHkx/Jfv31yqUiIgMMYZhJPUzih0WLFjAjh07+my7/vrrOeOMM/jud7+b1kACCiUp5WusIdLpJOgpBqCgLDWhBKC41IWvEdzhStr2vQPnZW6KloiIjAwFBQVHrdKel5dHWVnZcVdvT5WT6n0jfcUO1+LwlWA5XECcvKKTX6OkW/n4UQC4rQqam7el7HVFRESyhUZKUsjo2Is3UAYu8LoiOMzUZb5RU6pgcw0xVwWdbZm/IlpEREam9evXZ+xYGilJEcuyiIX3UhBMXNxakJ/a393KThkNgD+3AvPwvpS+toiISDZQKEmRw74wTmcj3mhi5btUrVHSrajcC1acqDMXT4eDVn84pa8vIiJiN4WSFNnXEgBnCyZda5RUFqb09Z0uE4/RAYDLN5qaxvaUvr6IiIjdFEpSZF9LAEc4QMjTtUbJ+LKUHyMvLwaAGayicd97KX99EREROymUpEhDUzPRQIxg18JpBaNzU36Msq7GfO5oBc371JhPRESGF4WSFPEd2pNYo6R74bSy1C/PW3FqJQAORwVNLQolIiLZbiA944a6VL5HhZIUiTXvxfAXETfdgEV+SepDyagzxgAQzKkg0rw75a8vIiKp4XK5APD7/TZXkn7d77H7PZ8MrVOSKu17yfGXQQF4zAimM/V5r2RM4uLZYE4prsOHsSwr6XbUIiKSfqZpUlxcTGNjIwC5ubnD7v9fW5aF3++nsbGR4uLilCxBr1CSApZlEY/UUhgspbMACvLTMwDlLXBhxgPEHF5cvkKaOsOMLvCk5VgiInJyKisTP7l3B5Phqri4uOe9niyFkhRo8UdwOhrwRmbQCRSUpicoJBrz+eiMeXH6y9lz8DCjC6rSciwRETk5hmFQVVVFeXk5kUjE7nLSwuVypbRJn0JJCuxvCeBwtuAkMfOmsCq1a5T0VljqpPMQOCOVNO19F05TKBERyWamaaa9u+5woQtdU2Bfix8r6ifiToSS4nGlaTtWxYTEcvOueCWNDVvTdhwREZFMUyhJgfqmFuLBXmuUlOen7VjlZ4wFIOaqoKVJ04JFRGT4UChJgc7GGsI+J4GuNUoKSlM/Hbhb6anlAARyK7AO1aTtOCIiIpmmUJIC0cO10JlHzJlowleQhoXTuiUa88WImR5cbXHi8eG/MI+IiIwMCiUpYLTX4vGPAsDtiOByp++CJtN04KEtcSxfKQfaAmk7loiISCYplJwky7KwQnsoDCSuJ8nPS/8x8/ISU8vMYAU1++rTf0AREZEMUCg5SW2BCKZ5kNxo10WuJelfzKykMtGYzxWr5NCebWk/noiISCYolJykfV1rlLisrjVKKtI386Zb9eREDxyHUcHB+tfTfjwREZFMUCg5Sfta/BD3EXUlQknR+LK0H3PUtHEABHMq8R96O+3HExERyQSFkpN0oLmNeCDas0ZJYUVB2o9ZOrYIgFBOCa7mphHRGltERIY/hZKT1HlwLyGfSdDTdU1JGtco6ZaT58IZ7wTA3VHMwfZQ2o8pIiKSbgolJylyeA/4vETciRGSdK5R0ltujh8AV6CKD/YfyMgxRURE0kmh5CQZbbW4fIk1SpxGBI83Mz0OSysTc4+dkWoaP/hbRo4pIiKSTgolJ8GyLAjtpah7jZLczF3bMeb0xAwc06imvn5Txo4rIiKSLgolJ6E9EMV0NJAb6bqepNidsWOXzxgPQNBbRaBxR8aOKyIiki4KJSdhX6sfh6sFd7xr5k0auwMfqWxCMQAhTwnOphbNwBERkSFPoeQk7GsJYFmdxLrXKBlXmrFje7xOXPGuHjgdpRzq1AwcEREZ2hRKTsKB5nbigUjPdODCqsKMHj/XGwTAFajkg30NGT22iIhIqimUnIT2xj2E/WbPwmmZWKOkt7KqxM9FzugYDn7wZkaPLSIikmoKJSch0ryXeKeHsKcYyNwaJd3GnJFYbt40qjigHjgiIjLEKZScBEdbLR5fYpTENKLk5LkyevyKMxMzcPy5VQQb3srosUVERFJNoeQkGMG95AcTDfjycuIYhpHR43f3wIm4C3E2tWf02CIiIqmmUDJI7cEIpqOB/FDX9SRFmVujpJvLY+KOtwDg7iihWTNwRERkCFMoGaSmjhCm63DPGiUF5bm21JGXmwgizmAV7+/XDBwRERm6kgoly5cv59xzz6WgoIDy8nKuvPJKdu3adcLnbdiwgdmzZ5OTk8Mpp5zCfffdN+iCs0WLP4JFB5hda5SMKbGljrIxiWnIrmg1DR9staUGERGRVEgqlGzYsIEbbriBTZs2sW7dOqLRKAsXLsTn8x3zOTU1NVx++eXMmzePrVu3cvvtt/Otb32L1atXn3Txdmrr9BMLRAjkJK4pKRpTbEsdY6dNAMBwVLH/gHrgiIjI0JVUS9vnn3++z/2HHnqI8vJytmzZwkUXXdTvc+677z7Gjx/PihUrAJg6dSqbN2/mzjvv5Kqrrhpc1Vkg0FJPyO+wbY2SbuUzxsMfD+HPrSK0/xlbahAREUmFk7qmpK0tscx5aemxl1d/7bXXWLhwYZ9tl112GZs3byYSifT7nFAoRHt7e59btgm1HyLucxH2JGbAZHqNkm4l1flgxYi68nA2+22pQUREJBUGHUosy2Lp0qVceOGFzJgx45j7NTQ0UFFR0WdbRUUF0WiUpqamfp+zfPlyioqKem7jxo0bbJlpE+44hMtfimWYOIiRW5D52TcATpeJx2oFEjNwWnxhW+oQERE5WYMOJTfeeCPbt2/n8ccfP+G+R67f0d3R9ljreixbtoy2traeW11d3WDLTBu/v54CX+J6klxPDMOR2TVKesvLT4w4OYNVfHDgoG11iIiInIykrinpdtNNN/Hss8+yceNGxo4de9x9KysraWjoO1W1sbERp9NJWVlZv8/xeDx4PJ7BlJYx/kADJeFSWoH8wkGdxpQZNa6Yw7vAGaum/v2tMDn7RpZEREROJKmREsuyuPHGG3nqqaf4y1/+wqRJk074nLlz57Ju3bo+29auXcucOXNwuTK7LHsqhSLN5MS6ugOPzrO1lrHTJib+MKuo2/+arbWIiIgMVlKh5IYbbuDRRx/lscceo6CggIaGBhoaGggEAj37LFu2jGuvvbbn/pIlS9i7dy9Lly5l586dPPjgg6xatYpbb701de/CBtHIYQxH9xolxbbWUj4jMTLiz60ivH+7rbWIiIgMVlKhZOXKlbS1tTF//nyqqqp6bk8++WTPPvX19dTW1vbcnzRpEmvWrGH9+vWcffbZ/OxnP+Puu+8e0tOBAczONoJda5QU2hxKiivzMKwoMWcOzqagrbWIiIgMVlIXQ3RfoHo8Dz/88FHbLr74Yt58881kDpX1HCF/zxolhWVeW2sxTQduq4WQMRpXRwltgQhF3qH705iIiIxM6n0zCMFIDGcwQtCTWFrerjVKessviALgDFXxwf5Gm6sRERFJnkLJILT4w5hhL5bDCVacPBs6BB9p9ITEqI0zXs3+97fYXI2IiEjyFEoGodUfwYwkZtw4CeAw7T+N46YnZkJZZhX7971qczUiIiLJs//bdAhq6fTjiOQD4HL0v1R+pn04A6eS4L4dNlcjIiKSPIWSQfC3NkF3KHHHba4moXB0LkY8Qtx042zKjqAkIiKSDIWSQQi0H8KMFQCQ482OU+hwGHhoBsDdWUJHUMFERESGluz4Rh1iQm2HMGOJkZLcguxZDr+gKDFqY4aqqDlwyOZqREREkqNQMgi+zv04rcSFrvlF9i4x31v5xNFAYgbOPs3AERGRIUahZBD8wYMYJEZK8krzba7mQ2PPPAWAuLOKutqXba5GREQkOQolgxAIHQJHYoQkd1SBzdV8qHxaomOzP7eCUO1bNlcjIiKSHIWSQYgHmog6u64pGV1oczUfKijLwREPYjmcmM3ZMStIRERkoBRKBsEItBFxd4WS0uy5psQwDDzGYQBcnaX4w1GbKxIRERk4hZJBMP0+Iq5EKMnJz67GdwVFiaaJrnAlHxxosrkaERGRgVMoSVI8buEIWcQdiTDiLbC/701vFadUAOCIV1P3nmbgiIjI0KFQkqSOYBQznOgKbFgRXG7T5or6GjfzNABi7mrq9r5kczUiIiIDp1CSpO4OwZBoxpdtRk+rBiDgHU2o9m2bqxERERk4hZIktfhCGNHENGCXmX1LuecVeXDE/GA4MA9pBo6IiAwdCiVJ6mw7jCPa3YzPsrmaoxmGQY4jMQPH7SsjGInZXJGIiMjAKJQkyd/WiKOr701ObnZdT9KtoCTxrxmp4n31wBERkSFCoSRJwbZDmPHsa8bXW+VpietKHFYV+zQDR0REhgiFkiR1dh7AtBKhpKA4e/re9Db+7MkARN3V1O3ZaHM1IiIiA6NQkiSfvx4HXR2Cy7Kn701vo6dUARD0jiK4Z6fN1YiIiAyMQkmS/MFGLEcijOSOyp6+N715C9yYsQ4AzCbD5mpEREQGRqEkSXFfExFXV4fg8iKbqzk2j6MFAJe/lFBUM3BERCT7KZQkyfK39PS9yS3JtbmaYysqS/ynNcNV1NSrB46IiGQ/hZIkmX5/z0hJtjXj661y8lgAHFSzVzNwRERkCFAoSZIZMsBInLZsDiUTzjkdgIinin0fbLC5GhERkRNTKElCKBrDDCXWJnHEg5hm9p6+UZMT3YJDOSUEa3bZXI2IiMiJZe+3ahZq80dwRBLXkTiNoM3VHJ8n14Uz2gqA2aT/zCIikv30bZWEFn8EsyuUZGMzviN5nK0AOP2lhKNqziciItlNoSQJbe1tGLHEGiVuT/Y14ztS0ahEbx5ntJq9B5ttrkZEROT4FEqS4G/9sBmfJ89pczUnVj1lPAAGVezZvdnmakRERI5PoSQJ/rZDOLqa8eUV5NhczYmNP+cMAEI5VdS+t97eYkRERE5AoSQJnR0HMEmEkvyS7GzG19uo08oBiLgLCdbstrkaERGR41MoSUJn5wGMrmZ8BWXZu8R8N5fHxBVNXEtiNps2VyMiInJ8CiVJ8AUOEje7fr4ZnZ3N+I7kdrYB4PKXEYlpBo6IiGQvhZIkxDobiXb3vakotreYASouT6w6a0arqD142OZqREREjk2hJBm+wx92CC722lzMwIyZOjHxh1HFnvc0A0dERLKXQkkSDH+AqDOxeFo2973pbfzsqUBiBs6ed/9iczUiIiLHlv2LbWQRM9CV4aw4Hu/QOHVlk8rAihN15RGr+cDuckRERI4p6ZGSjRs3csUVV1BdXY1hGDzzzDPH3X/9+vUYhnHU7d133x1szbawLAtHyA2AafkxHIbNFQ2M02XiijYB4DisGTgiIpK9kg4lPp+Ps846i3vuuSep5+3atYv6+vqe2+TJk5M9tK06QlHMSOI6EifZ3YzvSB53O5CYgROLZ//y+CIiMjIl/RvEokWLWLRoUdIHKi8vp7i4OOnnZYtWXwRHNA9c4HJmfzO+3koqPHQ2gBmrpvZQC5MqSu0uSURE5CgZu9B11qxZVFVVsWDBAl588cVMHTZlWjr9OKKJ6cDunKHx0023sdNPTfzhqOKDv79ubzEiIiLHkPZQUlVVxf3338/q1at56qmnmDJlCgsWLGDjxo3HfE4oFKK9vb3PzW6+tiYc8USHYE/e0Jh50238nMQMnGBOFXvfHXqBUERERoa0TyGZMmUKU6ZM6bk/d+5c6urquPPOO7nooov6fc7y5cv5yU9+ku7SkuJvPYjD6lrNtXBorFHSrWR8CUY8SsyZQ+j9vXaXIyIi0i9b1ik5//zz2b372A3ili1bRltbW8+trq4ug9X1r7PtAA6rq+9NaYHN1STHNB24Yt0zcIbGVGYRERl5bPmG2rp1K1VVVcd83OPx4PF4MljRibV17sdhJEZKCkdlfzO+I3nc7YStSpyBUuJxC8cQmdIsIiIjR9KhpLOzk/fee6/nfk1NDdu2baO0tJTx48ezbNky9u/fzyOPPALAihUrmDhxItOnTyccDvPoo4+yevVqVq9enbp3kQF+XwO55gwAcsuHXigprsqh40BiBs6+plbGl5fYXZKIiEgfSYeSzZs3c8kll/TcX7p0KQDXXXcdDz/8MPX19dTW1vY8Hg6HufXWW9m/fz9er5fp06fz3HPPcfnll6eg/MwJtzfgcs0FIK+82N5iBmH8madTdyCIZVbx3q7XGV/+CbtLEhER6SPpUDJ//nws69gLcD388MN97t92223cdtttSReWdTo+bMbnLcqxuZjkjZ8zjVdeeJNgTiV73/l/ME+hREREsosa8g2QEQgQNxPLzA+VZny9FY8pwoiHiZtuQh/ss7scERGRoyiUDFB3Mz4jHsHlGXo9ZBwOA1fsEACGZuCIiEgWUigZoD7N+IyhOXPF4+kAwBkcRVw9cEREJMsolAxAJBbHDOcC4DSGVjO+3krHJN6DGaviwOE2m6sRERHpS6FkAFr9ERyxxEWuLmfU5moGb8JZZwAQdyZm4IiIiGQThZIBaPWFMLpCiXvoTbzpMX7ONACCORXs3fE/NlcjIiLSl0LJALS3HcYRS6zmOhRn3nQrrCjAEQtiOZwEPthvdzkiIiJ9KJQMgK+1sacZX25xns3VDJ5hGLhijQA4DrttrkZERKQvhZIB6Gjdj0FX35vSQpurOTlubycAztCo4y6CJyIikmkKJQPQ1rEfHF2hZPTQ7hlTNjbxPox4FQ0t7TZXIyIi8iGFkgHo6DxA3Ex8medVFNtbzEmaePZ0AOJmFbt3brK5GhERkQ8plAxAtO0A4a6+N0OxQ3BvE86dCkDQW86e7f/P5mpEREQ+pFAyAPG2ZqLdoWQINuPrLa8sDzPqA8NB4P0DdpcjIiLSQ6FkIHxhLCPR7yYnb+hOCYauGThWVw+cFpcudhURkayhUDIADn+i140jFsB0Dv1TllsQBsAZKueD+iabqxEREUkY+t+wGWBGPIl/rYDNlaTGxLNOB8CMn8obr622uRoREZEEhZITsCwLRzhxHYnpGLrN+HqbtuhcAHwFkzj05nM2VyMiIpKgUHICvnAMM9rdjC9iczWpUVRZgDt2CMswyTlQQDASs7skERERhZITafWHMeKJUOLxGjZXkzqjKhPdjl3Bqbz59nabqxEREVEoOaG29g6MnmZ8w6dfzJkL5wAQyZnG9lcfsrkaERERhZIT6mw92NP3Jrc43+ZqUmfi+adhxMOEPSXE3/nA7nJEREQUSk6kvaX+w2Z8o4b2aq69OV0mec56AHJaJ9HQ0mlzRSIiMtIplJxAa3sdOBLXlBSXl9lcTWpNOqsKADM+nb9u+qPN1YiIyEinUHICHa37iA2TZnxHOvNTcwHw5Z9K3etar0REROylUHICodYDRFxd15QMs1BSMqYIV7QJy2HiqnMRi2vJeRERsY9CyQnE2pqJunIByC3w2FxN6pWN9gPg9k/lrfd0wauIiNhHoeQEjM5EnxisOO5cp73FpMGMjyWmBkc903jzZU0NFhER+yiUnIDR1YzPjPtxOIbP4mndTpk3FSMeIZRTiv9vO+wuR0RERjCFkhMwwi4ATMtvcyXp4XKb5Dr2A5BzeBztgbDNFYmIyEilUHICzogXANMI2VxJ+oyfnpjq7IxN4/U3/sfmakREZKRSKDmOaCyOEU2EEqdreDTj689Z/3ARkJgavPvl39tcjYiIjFQKJcfRFojg6Op7M5ya8R2pdHwxrmgzlsOF+UEcy9LUYBERyTyFkuNo6QxgWIlQ4h2G04G7GYZBUXEbAG7/FN7f32hzRSIiMhIplBxHZ0tjTyjJLS2wuZr0mrHgHABi7mm8vvG3NlcjIiIjkULJcbS31mMYib43RaNKbK4mvSZfMrNravAo2jZvsrscEREZgRRKjuNwax2WIzFSUlwxvJrxHcmd4ySHfQB4mioJRmI2VyQiIiONQslxtB/eS9SZCCX5laU2V5N+Y6ckfqIyo9PY8rc3bK5GRERGGoWS4wi19GrGV15kczXpN2vxJQAE8ifz1voHbK5GRERGGoWS44i0HCZuugHwFrhtrib9Rp06Cmf0MHGHC/7us7scEREZYRRKjqerGZ8Rj+LymDYXk36GYVCY3wyA2zeZgy3tNlckIiIjiULJ8XS1uzHjnRjG8F08rbdpl8wEIO6czisbn7C5GhERGUmSDiUbN27kiiuuoLq6GsMweOaZZ074nA0bNjB79mxycnI45ZRTuO+++wZTa8aZYScAjmHajK8/ZyycgxGPEvSO5uDL6oMjIiKZk3Qo8fl8nHXWWdxzzz0D2r+mpobLL7+cefPmsXXrVm6//Xa+9a1vsXr16qSLzTRHdzM+x/Btxnckj9eFJ14HgLuxjFhcS86LiEhmOJN9wqJFi1i0aNGA97/vvvsYP348K1asAGDq1Kls3ryZO++8k6uuuirZw2eUEc0FNzhdUbtLyaiqUz3U7E1MDX7r7+9y1hlT7S5JRERGgLRfU/Laa6+xcOHCPtsuu+wyNm/eTCTSf+fdUChEe3t7n1umBcIxHLHEaq7DuRlff2b/4wIAArmT2bx2aPzUJiIiQ1/aQ0lDQwMVFRV9tlVUVBCNRmlqaur3OcuXL6eoqKjnNm7cuHSXeZQWXwiHlQgl3sLhPx24t/IplTijrcRNN5F3+v9vJCIikmoZmX1z5MwVy7L63d5t2bJltLW19dzq6urSXuOR2lqbgMTCaXmlw3/htN4MwyDP2wCAp2MS7YGRc02NiIjYJ+2hpLKykoaGhj7bGhsbcTqdlJX130/G4/FQWFjY55ZpHS0NYCRCSVH58F9i/khnzJsGJKYGv/ryszZXIyIiI0HaQ8ncuXNZt25dn21r165lzpw5uFyudB9+0Jpaaol3NeMrqRxtczWZN+OTF2DEYwS9FexZ/3/tLkdEREaApENJZ2cn27ZtY9u2bUBiyu+2bduora0FEj+9XHvttT37L1myhL1797J06VJ27tzJgw8+yKpVq7j11ltT8w7SpLW5hlhX35v8qlE2V5N5OXluPLHEf1Pn/vyen9xERETSJelQsnnzZmbNmsWsWbMAWLp0KbNmzeKHP/whAPX19T0BBWDSpEmsWbOG9evXc/bZZ/Ozn/2Mu+++O+unAwcP7SPiSlzomjs68z8fZYPR4xP/OiNTeX/fPnuLERGRYS/pdUrmz59/3P/V/PDDDx+17eKLL+bNN99M9lC2irS2YBmJfje5+SNr9k232Z9eQN2vPyCQezqvvfAbTvvaz+0uSUREhjH1vjmGeHtixokjFsR0jczTVD1jAs5oG3HTg//NPXaXIyIiw9zI/LYdACOQmK5sxn02V2IfwzDwuvcD4GqfQCgysla2FRGRzFIoOQYj2NWMj5HTjK8/k+eeAoDlmMZfN79oczUiIjKcKZQcgyOSk/h3BDXj68/Ziy8BK04wt4p3n3/c7nJERGQYUyg5Bkcs0SHY6eq/P89I4S304okkZlMZtSPzgl8REckMhZJ+xOIWRlczPvcIa8bXn5KqxGiRGZ7CwabDNlcjIiLDlUJJP9oDEQwrsXBabqHH5mrsN+fTHwMg5J3CxhdW2VyNiIgMVwol/WjraO/pe1NQNrKa8fVn/KzTMKPtxJw5tGx6y+5yRERkmFIo6Ufb4Q+b8RVUlttcjf0Mh0GOs2vJ+dZq4nEtOS8iIqmnUNKPpuYa4mbimpLSqpHXjK8/k2aPBcAypvG3t7fYXI2IiAxHCiX9aGncQ9SZGCkpUigB4NyrFnZNDR7Dm398yO5yRERkGFIo6YfvUC3RrmZ83tEFNleTHXJL8/FE6gCIvx+zuRoRERmOFEr6EWluTfxhxcnJ09oc3QpGdQDgDE2mwx+wuRoRERluFEr6EetqxmfG/DgcWqek25zFlwAQyjmD9Wt/Z3M1IiIy3CiU9MefmF3iGMHN+PpzytwZmNEOYk4vDRs22V2OiIgMMwol/XCEzMS/I7wZ35EMh4HH2AOAo7kcy9LUYBERSR2Fkn4YasZ3TGNnlgFgMJW3dr9rczUiIjKcKJT0w+hpxhe1uZLsM/fzi8GKE8gdx0sP/cjuckREZBhRKOmHI56YDuzK1UWuR8ofVYA3vheA3L3T2FVTY3NFIiIyXCiUHCEYiYGa8R3XvC99BAB//gX8z/3ftbkaEREZLhRKjtDaGQASoaRwlJrx9WfyRWfije3BcjjJ2zuD3Xs0WiIiIidPoeQIrYcbsByJn2+KqypsriZ7XXRt92jJXNbd/z2bqxERkeFAoeQIh5r2EDcTIyWl1dU2V5O9Tps3A2+sBsthkrtnBu/t3Wt3SSIiMsQplBzh8MH3ibi6mvGNUTO+47nouvMB8Oefz7r7NFoiIiInR6HkCJ31tcTNxAWuuaPUjO94TrtwOt7YB1gOk5y90/mgttbukkREZAhTKDlCqLkFACMexZVj2lxN9rv4y3MBCBScxwsrNVoiIiKDp1ByhGhPM75ODEPrlJzIqR+djjf2PpaRGC2p2Vdnd0kiIjJEKZQcyadmfMma/5ULAAgUfIQX7tVoiYiIDI5CyZG6m/EZAZsLGTpOmTud3Nh7PaMle/fvt7skEREZghRKjuCIJC5yNRxBmysZWi7+ykcB8OV/hOfv0SqvIiKSPIWSI3Q34zPVjC8pidGS3WA4cO+dQe2BA3aXJCIiQ4xCyREMNeMbtPlfmweAv2AOf9ZoiYiIJEmhpJd43MJQM75Bm3TeNHJjf+8ZLanTaImIiCRBoaSXjmAEjMRIiZrxDc4lX78IAH/+bP78X8tsrkZERIYShZJeWluaiDsSIyUlY6psrmZomviRaeTGdoHhwFU7nX0N9XaXJCIiQ4RCSS+HDn3YjK9s3Dibqxm6Lvn6xQD488/hubs1WiIiIgOjUNJL44FdPc34iseU21zN0JUYLXm369qS6Rw42GB3SSIiMgQolPTSsb8Oy5FYPC23LN/maoa2S74xHwB//iz++CuNloiIyIkplPQSbEo043PEgjhdasZ3MibOmUZubGfPaEl9o0ZLRETk+BRKeom0JZaWN2Pqe5MKl35jPlhx/AXn8Kf//L7d5YiISJYbVCi59957mTRpEjk5OcyePZuXXnrpmPuuX78ewzCOur377ruDLjpdLF8cAIelUJIKE+ZMJze+EwBn7TTqDx20uSIREclmSYeSJ598kltuuYXvf//7bN26lXnz5rFo0SJqa2uP+7xdu3ZRX1/fc5s8efKgi06brmZ8BmrGlyqXLlnQNVoyi+fu0miJiIgcW9Kh5K677uKrX/0qX/va15g6dSorVqxg3LhxrFy58rjPKy8vp7Kysudmmtl3zYYRyUn8a6oZX6pMmD2tZ7TErJvGwaZGmysSEZFslVQoCYfDbNmyhYULF/bZvnDhQl599dXjPnfWrFlUVVWxYMECXnzxxePuGwqFaG9v73PLBCOWCCWmK5KR440UC/6la7Qk/2z+9EuNloiISP+SCiVNTU3EYjEqKir6bK+oqKChof/ZFVVVVdx///2sXr2ap556iilTprBgwQI2btx4zOMsX76coqKintu4DC1k1tOMz6tmfKk0/pxp5MbfAcBRO53Gw002VyQiItloUBe6GkbfL23Lso7a1m3KlCl8/etf55xzzmHu3Lnce++9fPKTn+TOO+885usvW7aMtra2nltdXd1gyhyExNokHjXjS7kFN3ys69qSmfzxF7fbXY6IiGShpELJqFGjME3zqFGRxsbGo0ZPjuf8889n9+7dx3zc4/FQWFjY55Zu4Wi8VzO+9B9vpBl/9jRy428D4Ng3naaWZpsrEhGRbJNUKHG73cyePZt169b12b5u3TouuOCCAb/O1q1bqarKroZ3La2tPc34Ro2ttrma4eljN3y869qSM3n637XKq4iI9OVM9glLly7lS1/6EnPmzGHu3Lncf//91NbWsmTJEiDx08v+/ft55JFHAFixYgUTJ05k+vTphMNhHn30UVavXs3q1atT+05O0qHGPcSciZGSURMn2VzN8DTu7GnkxtfiN2di7p/L6sd/xVWfv9nuskREJEskHUquvvpqmpub+elPf0p9fT0zZsxgzZo1TJgwAYD6+vo+a5aEw2FuvfVW9u/fj9frZfr06Tz33HNcfvnlqXsXKXCwbhdR52gASsdppCRdFt16Jc/8xw4CeRMw/rSdDVV/5OL5V9hdloiIZAHDsizL7iJOpL29naKiItra2tJ2fcnqVXfQ8MZHwIrzL/deisPUCvzp8vbaV9nwf9qxHG5yO17ivB/+I9OmzLS7LBERSbFkv7/1zdsl0HgYADPmVyBJs+kLL+DMuYGu2TjzeO3nD2oJehERUSjpFm5NLC3vUDO+jJh3/T8yfsJeAIJ5/8Cfvnc7/oBW0hURGckUSrrEA2rGl2lX3P5VSnMTi6pFzM/xu//1dWKxuM1ViYiIXRRKugW7mvEZfpsLGVk+d+cN5Fk7sRwmVvBzrPrBN+wuSUREbKJQ0i3iBsAwQzYXMrIYDoMv3v0NcsIfEHN64eAiHrnru3aXJSIiNlAo6WLEcgEwXVGbKxl5nB4nn7vzs7hD9UTcxUT+dibPPHq33WWJiEiGKZR0MayuUOLN+hnSw1JeaSFXLPsozkgrQW81h5938/KGP9pdloiIZJBCSY/Eaq45BWrGZ5fK0ycy/7qxOGIBAvmns/u/d7Jr1w67yxIRkQxRKCHR5dgyEn1vCkYX2VzNyDbloo9w1qURjHgMf/4cXvn54xw6fMjuskREJAMUSoAOf4C4IzFSUqpmfLa74PNXMuH0xBomgbyP8X+/8xOCIV2ALCIy3CmUAI0H9xJzJkZKxpxyms3VCMAnb/0ao4q2AxByX8nvvn0D8biu9xERGc4USoD6mneIOXMAKJkw1uZqpNtn77iZfLaB4SAa+QwPLvsXu0sSEZE0UigBDtcmfiow4jFyCr02VyPdDMPgi/fcjDe8k7jpJtZ0Ob//xe12lyUiImmiUAL4DjYDYMY6MQzD5mqkN9Np8oX//DKe4F6irnz875zNHx+9x+6yREQkDRRK6N2Mr9PmSqQ/OUV5/ONPP44rdIhwziga1+bxv3/zb1iWrjERERlOFEqAmC/RBM5AzfiyVdn4sVy6ZDLOSAfB3Am0/PVMHvjGEvY37LO7NBERSRGFEoBQ4jQ4CNhciBzPaeedwyVfqcYdqiPqyiPsuJq1t/4fVj/wC7tLExGRFFAoAYh0reJqBu2tQ07o9I+ey/X//UVK87dhxKP482fS9NpUHvjGEhoa6+0uT0REToJCCWDEEtOBHWrGNyQ4nSafv3Mp87+Q3zNqEjI+y5qlT/LMQ/9pd3kiIjJICiUAVmI1VzNXF04OJdPmX8D1919DSe6bGPEogfyZHHz5NB5Y8k0ONTfaXZ6IiCRJoQT4sBmf2+Y6JFlOl5Mv3HUrF3/O++GoCZ/hjzc/xrO/09RhEZGhRKEEsLr63uSOKrS5Ehms6ZfO48u/+QLF3i09oyb16yfwwJJvcri12e7yRERkAEZ8KIlEY8TMRN+bMjXjG9JcbhfX/Of/Yt5nPX1GTZ658Xesefw+u8sTEZETGPGh5ODBOqJdzfiqT5tsczWSCmd+7GK+/JsvUJTzRs+oSd3/jOGBf7mR1rbDdpcnIiLHMOJDyb6/v4XlcAJQecqpNlcjqeJyu/jiiu9ywWeceIJdoybWp3nqxkd4ZtUvicfidpcoIiJHGPGhpLl2DwCOWAh3Xo69xUjKnb3wUq67/wsUeV5PjJrkzeTAX2fy8JdX8NBt32L/gb12lygiIl1GfCjpbPiwGZ8MTy63iy/+ahkXXGXgCe7GMkwCeWfjb7+SNbdvYtVXb2Xd6ofUS0dExGZOuwuwW6jFD4Ajrr43w93Zl32csy+Dt158iTcee5mQdSbhnNHA5ex+IUrdM/+BObGRK779XUpLR9tdrojIiDPiQ0nMHwPAsBRKRooZl8xjxiXzCLR18sc7V9JeV0YoZyKBvDlwCJ5aug7DtZkzPnsRH/34lXaXKyIyYoz4UBIPmokfsQy/3aVIhnmL8vnsz/4XAG/83+fY8exbhBwzCeVUAp9i+/8Os+uRfyPnDB+Ll95OnjfP3oJFRIa5ER9KiLjAA4aa8Y1o5y7+JOcu/iTtjU386c778TWOI5wzhkDe+QTq4Ikl/we82zn3+iuYed58u8sVERmWRnwoMeLexL/uiM2VSDYoLB/FF35xO5ZlsfH3T/DeX/YRNs8k6B0HjOOVBwJsvuceLE8t+ZNzmffFL1M9ZqLdZYuIDAsjPpRg5QJgajaw9GIYBhd/8fNc/EU4tGcvz6/4LcG20wh7KgnkTQOmEayBZ3/0Nu7Q/8Xy7qNoWgmXfvGrlI6qsLt8EZEhacSHEovEaq7uApfNlUi2Gj1xAl9a8UPi8ThvrlnDO/+zhUhrMRHnqcScuQScZwJnEnwX/vDd13GFd0NePWWzqllwzdfIz1dPJRGRgVAo6WrGlzeqyOZKJNs5HA7mfOpTzPnUpwCIRmO8/tRq3nvpHaLtowi7TyHqyifqmgXMYt9WeOz1/4czuhsKGqn6yGnMv/o6vDlee9+IiEiWGtGhxLKsnmZ8xdWVNlcjQ43TafLRz36Wj342cT8cDPLS409Q+8Ye4r4Kwu5JRNxFRNxzIAYfvAa1G9fhjNSD0YiR20FOVQ4TPnIOsy/5BDke/YYoIiPbiA4lzU0HiToT15SMmXy6zdXIUOfOyWHB9V+G6xP3gx0d/OV3v+fg9gZigWrCnoldIymTgclgQeAAtDwDb/1hLc5IPYbRCF1hZdJHZjP70kW43W4b35WISOaM6FBSt3MHGCYAY08/w+ZqZLjJKSjg8m8u6bnf2dLC5ufXcGD7+wSbLKxwMXFHJWF32THDyvY/PI8z0gDGQYzcTlzFDvKrS6meOpVp582joEA/O4rI8DGiQ0ljzR7gVMyoD7dXQ+eSXvklJcz//DXw+b7bO1ta2Pzn59i//T1CzWCFS4g7Kgi7R3WFldOA0xJhpQXaW+DA27Dlf2/CFWnFEWsFow3MDoycMK4iJ4VjyqieNpVpH5lHXn6+HW9XRCRpIzqUtB9sAk5VMz6xVX5JCfO/8EX4Qt/tnS2H+eua56jf8X5XWCnCooiYWUzUVYjlcBH2jAZ69emJfxhc9r0FbzzxyofBhQ5wBMAZwnBHMb0G7kI33lHFlFZXMWbKGYw9bSoetyeD715E5EMjOpSEWhJhRM34JBvll5Ry6TVf6vexkK+TXW/8ldq33qb9wGFCrVGsoBsrln/i4AIQB3xdt3rYvwN2vHAYI74eZ9SHI+br6gflA0cAwwyBM4bDZWF4DJxeJ+78HHKK88kvK6OkopLyCRMYXTUBl1vT60VkcAYVSu69917+4z/+g/r6eqZPn86KFSuYN2/eMfffsGEDS5cu5e2336a6uprbbruNJUuWHHP/TIn64l1/KZTI0OLJy2fm/EuZOf/SY+4T7Oxk5183se+dnbTvaybSGSUeckDUjRX3ArlYRh5xRx4xM4+46cZyuIi4i4Hi/l802nXzAU1HPliLEa/BEQ/iiAVxWEEMKwhWCAiBEcUwIlhGDMOMgiOO4YxjOMFwgukxMT1OXLke3LlePPl5eIsKKSgrI6+giLySIgpKRpOfX4jDNE/6HIpI9kk6lDz55JPccsst3HvvvXz0ox/lN7/5DYsWLeKdd95h/PjxR+1fU1PD5Zdfzte//nUeffRRXnnlFb75zW8yevRorrrqqpS8icGKB41EMz7UjE+Gn5z8fGZd+jFmXfqxAe3f0tDA3nfeonHvHtoPNhNq9RH1xYgFu4KM5QLLA5YHy8hJ3Bwe4o4cYs7E2iuWwyTmyCPmHGDzwjgQ7rod9/8MO7pu+8CK4YhHMawoRjyKw4qAlbgP3f/G6E5QBjEwYljEMYwYGHEsI45BHBwWGBaGwwJHHBxgOMBwGImw5DAwnA4M08AwHZgOBw6XA8M0MZ0mDpcT0+nCdDkx3S5cbhem243p8uDO8eDK8WC63LjcnsQtx4Pb48blycHtycPjzcGtqeAiPQzLsqxknnDeeedxzjnnsHLlyp5tU6dO5corr2T58uVH7f/d736XZ599lp07d/ZsW7JkCX/729947bXXBnTM9vZ2ioqKaGtro7AwdatjrvrKjwm6L8Ib2cBXVv0kZa8rMtJEQ2GaGvZxsHYPrQ0H6Ww+TKC1g1BnkFggSjxsYUXBijog7oC4E8tyguUEywW4wHBh4cYy3FiGC8vhJu5wE3e4EklhuLLiGJYFxDGsxA2srn/jYFmJAIXV9bfVsx36Pkb330dt732jz34fbut+fuJxo2uLYST2sbC6tnW/fq+/DavnFYxer2kZif16Xt+g53V6b+vz3A8f7PP4h3vRte3Dry7DAMvo9ZDR92utz2O9Gb03Whj0t2PvGoxEiMVIvA+j145Hvnifeuh6d93PsXp26L7Xs79xdA1GP3/h6Ofx3jUYfe4cvbn363dtnDL/As65dCGplOz3d1IjJeFwmC1btvC9732vz/aFCxfy6quv9vuc1157jYUL+77Jyy67jFWrVhGJRHC5jv79ORQKEQqFeu63t7cnU+bAxdSMTyQVnB43lRNOoXLCKSl/bcuyCHR20NbcRGdrC76ONoLt7QR9nYR8AcKBEJFAkGgwRCwSIx6OEY/GsCIW8ShYcQviBsQNrHjXt5dlgOXAshxgOcAyMXBgYZL4//YmWGbXkgEOLBxd27uGUjAT24yu5xgOLKNrP8NM/G04sAzzxIHKcHR9aZpYx99TBurIE6kTOyBvv/CXlIeSZCUVSpqamojFYlRU9G04VlFRQUNDQ7/PaWho6Hf/aDRKU1MTVVVVRz1n+fLl/OQn6R+5MJwBPMF6XBX6fVokWxmGQW5BIbkFQ7OHUCwWIxIKEQ4GCAUDRMIhosEQ4VCISCREJBQiHg4TjkSIRcLEIpGuW5RYLEY8GiEeixOPx4lHYsTjMax4HOIWsVgMK2ZhxeO9biT+jXWNcFiAZZEYgLGwLAviiTGODwdcrJ4vbqt7gKX7DRw5mAJgGUc8bvR6Qq//CW71/rvXY73+tnqNGPTev2d7r2EOy+q1b9dfFh8e2zCMxPvps0evWnvGQvpu6zVw0XX/iNEKq/9Rh6PvH2e4xDrO/aNGaI54ncSQ1dHb6Gf7MWvpf3vvcSvPOPtbYAzqQlfjiJNgWdZR2060f3/buy1btoylS5f23G9vb2fcuHGDKfW4vrLqxyl/TRGR3kzTxMzNJSc31+5SRLJeUqFk1KhRmKZ51KhIY2PjUaMh3SorK/vd3+l0UlZW1u9zPB4PHo/WShARERlJkrp6zO12M3v2bNatW9dn+7p167jgggv6fc7cuXOP2n/t2rXMmTOn3+tJREREZGRK+pL2pUuX8sADD/Dggw+yc+dOvv3tb1NbW9uz7siyZcu49tpre/ZfsmQJe/fuZenSpezcuZMHH3yQVatWceutt6buXYiIiMiQl/Q1JVdffTXNzc389Kc/pb6+nhkzZrBmzRomTJgAQH19PbW1tT37T5o0iTVr1vDtb3+bX//611RXV3P33XfbvkaJiIiIZJek1ymxQ7rWKREREZH0Sfb7exivSCQiIiJDiUKJiIiIZAWFEhEREckKCiUiIiKSFRRKREREJCsolIiIiEhWUCgRERGRrKBQIiIiIllBoURERESyQtLLzNuhe9HZ9vZ2mysRERGRger+3h7o4vFDIpR0dHQAMG7cOJsrERERkWR1dHRQVFR0wv2GRO+beDzOgQMHKCgowDCMlL1ue3s748aNo66uTj11kqDzNjg6b4Oj85Y8nbPB0XkbnOOdN8uy6OjooLq6GofjxFeMDImREofDwdixY9P2+oWFhfoADoLO2+DovA2OzlvydM4GR+dtcI513gYyQtJNF7qKiIhIVlAoERERkawwokOJx+PhRz/6ER6Px+5ShhSdt8HReRscnbfk6ZwNjs7b4KTyvA2JC11FRERk+BvRIyUiIiKSPRRKREREJCsolIiIiEhWGNGh5N5772XSpEnk5OQwe/ZsXnrpJbtLymo//vGPMQyjz62ystLusrLOxo0bueKKK6iursYwDJ555pk+j1uWxY9//GOqq6vxer3Mnz+ft99+255is8SJztmXv/zloz57559/vj3FZonly5dz7rnnUlBQQHl5OVdeeSW7du3qs48+a0cbyHnT5+1oK1euZObMmT1rkcydO5c///nPPY+n6rM2YkPJk08+yS233ML3v/99tm7dyrx581i0aBG1tbV2l5bVpk+fTn19fc9tx44ddpeUdXw+H2eddRb33HNPv4//4he/4K677uKee+7hjTfeoLKyko9//OM97RRGohOdM4BPfOITfT57a9asyWCF2WfDhg3ccMMNbNq0iXXr1hGNRlm4cCE+n69nH33WjjaQ8wb6vB1p7Nix3HHHHWzevJnNmzdz6aWXsnjx4p7gkbLPmjVCfeQjH7GWLFnSZ9sZZ5xhfe9737Opouz3ox/9yDrrrLPsLmNIAaynn3665348HrcqKyutO+64o2dbMBi0ioqKrPvuu8+GCrPPkefMsizruuuusxYvXmxLPUNFY2OjBVgbNmywLEuftYE68rxZlj5vA1VSUmI98MADKf2sjciRknA4zJYtW1i4cGGf7QsXLuTVV1+1qaqhYffu3VRXVzNp0iQ+97nP8cEHH9hd0pBSU1NDQ0NDn8+ex+Ph4osv1mfvBNavX095eTmnn346X//612lsbLS7pKzS1tYGQGlpKaDP2kAded666fN2bLFYjCeeeAKfz8fcuXNT+lkbkaGkqamJWCxGRUVFn+0VFRU0NDTYVFX2O++883jkkUd44YUX+O///m8aGhq44IILaG5utru0IaP786XPXnIWLVrE73//e/7yl7/wy1/+kjfeeINLL72UUChkd2lZwbIsli5dyoUXXsiMGTMAfdYGor/zBvq8HcuOHTvIz8/H4/GwZMkSnn76aaZNm5bSz9qQaMiXLkd2HLYsK6VdiIebRYsW9fx95plnMnfuXE499VR++9vfsnTpUhsrG3r02UvO1Vdf3fP3jBkzmDNnDhMmTOC5557j05/+tI2VZYcbb7yR7du38/LLLx/1mD5rx3as86bPW/+mTJnCtm3baG1tZfXq1Vx33XVs2LCh5/FUfNZG5EjJqFGjME3zqATX2Nh4VNKTY8vLy+PMM89k9+7ddpcyZHTPVtJn7+RUVVUxYcIEffaAm266iWeffZYXX3yxTzd1fdaO71jnrT/6vCW43W5OO+005syZw/LlyznrrLP41a9+ldLP2ogMJW63m9mzZ7Nu3bo+29etW8cFF1xgU1VDTygUYufOnVRVVdldypAxadIkKisr+3z2wuEwGzZs0GcvCc3NzdTV1Y3oz55lWdx444089dRT/OUvf2HSpEl9HtdnrX8nOm/90eetf5ZlEQqFUvtZS9FFuEPOE088YblcLmvVqlXWO++8Y91yyy1WXl6etWfPHrtLy1rf+c53rPXr11sffPCBtWnTJutTn/qUVVBQoHN2hI6ODmvr1q3W1q1bLcC66667rK1bt1p79+61LMuy7rjjDquoqMh66qmnrB07dlif//znraqqKqu9vd3myu1zvHPW0dFhfec737FeffVVq6amxnrxxRetuXPnWmPGjBnR5+xf/uVfrKKiImv9+vVWfX19z83v9/fso8/a0U503vR569+yZcusjRs3WjU1Ndb27dut22+/3XI4HNbatWsty0rdZ23EhhLLsqxf//rX1oQJEyy3222dc845faaEydGuvvpqq6qqynK5XFZ1dbX16U9/2nr77bftLivrvPjiixZw1O26666zLCsxVfNHP/qRVVlZaXk8Huuiiy6yduzYYW/RNjveOfP7/dbChQut0aNHWy6Xyxo/frx13XXXWbW1tXaXbav+zhdgPfTQQz376LN2tBOdN33e+veVr3yl5/ty9OjR1oIFC3oCiWWl7rOmLsEiIiKSFUbkNSUiIiKSfRRKREREJCsolIiIiEhWUCgRERGRrKBQIiIiIllBoURERESygkKJiIiIZAWFEhEREckKCiUiw8j8+fO55ZZbMna8L3/5y1x55ZUZO95gZfq8iMjgKJSIiIhIVlAoEREZhFgsRjwet7sMkWFFoURkmAqHw9x2222MGTOGvLw8zjvvPNavXw9AW1sbXq+X559/vs9znnrqKfLy8ujs7ARg//79XH311ZSUlFBWVsbixYvZs2fPoOqZP38+3/rWt7jtttsoLS2lsrKSH//4xz2P79mzB8Mw2LZtW8+21tZWDMPoqXv9+vUYhsELL7zArFmz8Hq9XHrppTQ2NvLnP/+ZqVOnUlhYyOc//3n8fn+f40ejUW688UaKi4spKyvjBz/4Ab1bfx3vfAE8/PDDFBcX86c//Ylp06bh8XjYu3fvoM6FiPRPoURkmLr++ut55ZVXeOKJJ9i+fTv/9E//xCc+8Ql2795NUVERn/zkJ/n973/f5zmPPfYYixcvJj8/H7/fzyWXXEJ+fj4bN27k5ZdfJj8/n0984hOEw+FB1fTb3/6WvLw8Xn/9dX7xi1/w05/+lHXr1iX9Oj/+8Y+55557ePXVV6mrq+Ozn/0sK1as4LHHHuO5555j3bp1/Nd//ddRx3Y6nbz++uvcfffd/Od//icPPPDAgM5XN7/fz/Lly3nggQd4++23KS8vH9R5EJFjSFlfYxGx3cUXX2zdfPPN1nvvvWcZhmHt37+/z+MLFiywli1bZlmWZT311FNWfn6+5fP5LMuyrLa2NisnJ8d67rnnLMuyrFWrVllTpkyx4vF4z/NDoZDl9XqtF154wbIsy7ruuuusxYsXD7i2Cy+8sM+2c8891/rud79rWZZl1dTUWIC1devWnsdbWloswHrxxRcty7KsF1980QKs//mf/+nZZ/ny5RZgvf/++z3bvvGNb1iXXXZZn2NPnTq1z3v57ne/a02dOtWyLGtA5+uhhx6yAGvbtm0Der8ikjynrYlIRNLizTffxLIsTj/99D7bQ6EQZWVlAHzyk5/E6XTy7LPP8rnPfY7Vq1dTUFDAwoULAdiyZQvvvfceBQUFfV4jGAzy/vvvD6qumTNn9rlfVVVFY2PjSb1ORUUFubm5nHLKKX22/fWvf+3znPPPPx/DMHruz507l1/+8pfEYrEBnS8At9t91HsQkdRRKBEZhuLxOKZpsmXLFkzT7PNYfn4+kPiC/cxnPsNjjz3G5z73OR577DGuvvpqnE5nz2vMnj37qJ94AEaPHj2oulwuV5/7hmH0XCzqcCR+TbZ6XecRiURO+DqGYRz3dQdiIOcLwOv19gk2IpJaCiUiw9CsWbOIxWI0NjYyb968Y+53zTXXsHDhQt5++21efPFFfvazn/U8ds455/Dkk09SXl5OYWFh2mvuDjr19fXMmjULoM9Frydr06ZNR92fPHkypmkO+HyJSHrpQleRYej000/nmmuu4dprr+Wpp56ipqaGN954g3//939nzZo1PftdfPHFVFRUcM011zBx4kTOP//8nseuueYaRo0axeLFi3nppZeoqalhw4YN3Hzzzezbty/lNXu9Xs4//3zuuOMO3nnnHTZu3MgPfvCDlL1+XV0dS5cuZdeuXTz++OP813/9FzfffDMw8PMlIumlUCIyTD300ENce+21fOc732HKlCn8wz/8A6+//jrjxo3r2ccwDD7/+c/zt7/9jWuuuabP83Nzc9m4cSPjx4/n05/+NFOnTuUrX/kKgUAgbSMnDz74IJFIhDlz5nDzzTfz85//PGWvfe211xIIBPjIRz7CDTfcwE033cQ///M/9zw+kPMlIullWL1/wBURERGxiUZKREREJCsolIjISautrSU/P/+Yt9raWrtLFJEhQD/fiMhJi0ajx11+fuLEiT1TjUVEjkWhRERERLKCfr4RERGRrKBQIiIiIllBoURERESygkKJiIiIZAWFEhEREckKCiUiIiKSFRRKREREJCsolIiIiEhW+P9t8cYR4KajJAAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "spontaneous_recombination_rate_old.loc[1,0,:].plot()" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiUAAAHACAYAAAB9DBhHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABaB0lEQVR4nO3deZhU9Z3v8fepU0tX7wv0xq4isogiGMWIoiQYTBycmIlJTDRmmTBRoyFeE0ye7DM4GeMQx4hxRI0xLjcX9ZpIFOZGwA0jCAEVCWpDN9BN002vtS/n/lHdbTc00NVU1anu/ryepx66Tp2q861jzdQnvzq/39ewLMtCRERExGYOuwsQERERAYUSERERyRIKJSIiIpIVFEpEREQkKyiUiIiISFZQKBEREZGsoFAiIiIiWUGhRERERLKCQomIiIhkBYUSERERyQpDKpRs3LiRK664gurqagzD4Jlnnkn6NSzL4s477+T000/H4/Ewbtw4/u3f/i31xYqIiEhSnHYXkAyfz8dZZ53F9ddfz1VXXTWo17j55ptZu3Ytd955J2eeeSZtbW00NTWluFIRERFJljFUG/IZhsHTTz/NlVde2bMtHA7zgx/8gN///ve0trYyY8YM/v3f/5358+cDsHPnTmbOnMlbb73FlClT7ClcRERE+jWkfr45keuvv55XXnmFJ554gu3bt/NP//RPfOITn2D37t0A/PGPf+SUU07hT3/6E5MmTWLixIl87Wtf4/DhwzZXLiIiIsMmlLz//vs8/vjj/OEPf2DevHmceuqp3HrrrVx44YU89NBDAHzwwQfs3buXP/zhDzzyyCM8/PDDbNmyhc985jM2Vy8iIiJD6pqS43nzzTexLIvTTz+9z/ZQKERZWRkA8XicUCjEI4880rPfqlWrmD17Nrt27dJPOiIiIjYaNqEkHo9jmiZbtmzBNM0+j+Xn5wNQVVWF0+nsE1ymTp0KQG1trUKJiIiIjYZNKJk1axaxWIzGxkbmzZvX7z4f/ehHiUajvP/++5x66qkA/P3vfwdgwoQJGatVREREjjakZt90dnby3nvvAYkQctddd3HJJZdQWlrK+PHj+eIXv8grr7zCL3/5S2bNmkVTUxN/+ctfOPPMM7n88suJx+Oce+655Ofns2LFCuLxODfccAOFhYWsXbvW5ncnIiIysg2pULJ+/XouueSSo7Zfd911PPzww0QiEX7+85/zyCOPsH//fsrKypg7dy4/+clPOPPMMwE4cOAAN910E2vXriUvL49Fixbxy1/+ktLS0ky/HREREellSIUSERERGb6GzZRgERERGdoUSkRERCQrDInZN/F4nAMHDlBQUIBhGHaXIyIiIgNgWRYdHR1UV1fjcJx4HGRIhJIDBw4wbtw4u8sQERGRQairq2Ps2LEn3G9IhJKCggIg8aYKCwttrkZEREQGor29nXHjxvV8j5/IkAgl3T/ZFBYWKpSIiIgMMQO99EIXuoqIiEhWUCgRERGRrKBQIiIiIllhSFxTIiIiMlTFYjEikYjdZaSFy+XCNM2UvZ5CiYiISBpYlkVDQwOtra12l5JWxcXFVFZWpmQdMYUSERGRNOgOJOXl5eTm5g67xT8ty8Lv99PY2AhAVVXVSb+mQomIiEiKxWKxnkBSVlZmdzlp4/V6AWhsbKS8vPykf8rRha4iIiIp1n0NSW5urs2VpF/3e0zFdTMKJSIiImky3H6y6U8q36NCiYiIiGQFhRIRERHJCgolIiIi0se9997LpEmTyMnJYfbs2bz00ksZOa5CiU3icYtdDR3E45bdpYiIiPR48sknueWWW/j+97/P1q1bmTdvHosWLaK2tjbtx1YoscmDr9Rw2YqN/G7TXrtLERER6XHXXXfx1a9+la997WtMnTqVFStWMG7cOFauXJn2Y2udEpv8bV8bAG/WtnDdBRPtLUZERNLOsiwCkVjGj+t1mQOeIRMOh9myZQvf+973+mxfuHAhr776ajrK60OhxCYHDncy2djHnkMFdpciIiIZEIjEmPbDFzJ+3Hd+ehm57oF93Tc1NRGLxaioqOizvaKigoaGhnSU14d+vrHJhc3/h3We2ziv+WksS9eViIhI9jhyZMWyrIysuaKREhsEIzFOjewCE6bH3uWwL0xZvsfuskREJI28LpN3fnqZLccdqFGjRmGa5lGjIo2NjUeNnqSDQokNDrQGqDQO8bbbxfjgAWqafAolIiLDnGEYA/4ZxS5ut5vZs2ezbt06/vEf/7Fn+7p161i8eHHaj6+fb2ywryXA68XtfG5MFVuK2qk51Gl3SSIiIgAsXbqUBx54gAcffJCdO3fy7W9/m9raWpYsWZL2Y2d3ZBum6ptbqfMkrsDelWMwsaEOGG9vUSIiIsDVV19Nc3MzP/3pT6mvr2fGjBmsWbOGCRMmpP3YCiU2aD9YQwNOpu+Js2+0i8qDfwc+andZIiIiAHzzm9/km9/8ZsaPq59vbBBp3su0bSY/ejzO6dtNjMPv212SiIiI7RRKbBBrrWF0U2Jq1YQGcAR2a1qwiIiMeAolNgiH9lDW7qGhfDaVLW6crgMcbA/ZXZaIiIitFEoyLBSNYdFALHcB70z7CrHcS7BczdQ0+ewuTURExFYKJRlW3xrEYRwm5q4CwJ87jlgoyN6mdpsrExERsZdCSYbtawngDAUIekoB8OdWEOw0aT7wgc2ViYiI2EuhJMMONLdh+aMEcxKhJOAdTbzNReTg322uTERExF4KJRnWdnAPsU4vEXeiO7DlcJLbUY7R+p7NlYmIiNhLoSTDwk17cARG9dlW6ivHCr9HLK5pwSIiMnIplGRYvG0P3iNCSX6kAtPdwIHWgE1ViYiI2E+hJMMiwRrywiV9N5oVxI02TQsWERHbbdy4kSuuuILq6moMw+CZZ57J2LEVSjIoHI0Tox5PLHGRa443saqrP7eCuC9M7aEWO8sTERHB5/Nx1llncc8992T82GrIl0H1bQFMRws4ygAYN6WY3dta8OdWEOpw0rb/78Dp9hYpIiIj2qJFi1i0aJEtx1YoyaB9LQHMkL9nOvD4syvZva2FiCsfs62Q2KHdNlcoIiJpY1kQ8Wf+uK5cMIzMH3cQFEoy6EBzO0Yg0hNKSqvy8LojBMIuijorCbTtsrlCERFJm4gf/q0688e9/QC48zJ/3EFI6pqSlStXMnPmTAoLCyksLGTu3Ln8+c9/Pu5zNmzYwOzZs8nJyeGUU07hvvvuO6mCh7K2g3uJdnoIe4oBKCjLobjEBUBhsJxYfA+RWNzGCkVEROyT1EjJ2LFjueOOOzjttNMA+O1vf8vixYvZunUr06dPP2r/mpoaLr/8cr7+9a/z6KOP8sorr/DNb36T0aNHc9VVV6XmHQwhoaY9GL5R4AEHUXLyXJSOLaT+YDseqxKHaxt1h/2cMjrf7lJFRCTVXLmJUQs7jjtEJBVKrrjiij73//Vf/5WVK1eyadOmfkPJfffdx/jx41mxYgUAU6dOZfPmzdx5550jMpTEW/eQEywj7oFcTxTDMCg7ZTRsaSfkqcCIdLCn2adQIiIyHBnGkPkZxS6DvqYkFovxhz/8AZ/Px9y5c/vd57XXXmPhwoV9tl122WWsWrWKSCSCy+Xq93mhUIhQKNRzv719eHTQjQT3UBwuoR0oLEyc+pIxieXm/bkVRDvj7Gs4CGdU2FiliIiMZJ2dnbz33oetT2pqati2bRulpaWMHz8+rcdOep2SHTt2kJ+fj8fjYcmSJTz99NNMmzat330bGhqoqOj7BVtRUUE0GqWpqemYx1i+fDlFRUU9t3HjxiVbZtaJxOJY1gE80cRFrgWjEsNpJZWJ1Bz0lhFry6HzgC52FRER+2zevJlZs2Yxa9YsAJYuXcqsWbP44Q9/mPZjJz1SMmXKFLZt20ZrayurV6/muuuuY8OGDccMJsYR05Asy+p3e2/Lli1j6dKlPffb29uHfDBpaAviMA/jMBJrlBSNKQYgt8iNaUSJ4cTTWU68SdOCRUTEPvPnz+/5rs60pEOJ2+3uudB1zpw5vPHGG/zqV7/iN7/5zVH7VlZW0tDQ0GdbY2MjTqeTsrKyYx7D4/Hg8XiSLS2r1bX4McM+Ql3TgYvGJ96/YRgU5Vsc7oBiXyUtnRopERGRkemkl5m3LKvP9R+9zZ07l3Xr1vXZtnbtWubMmXPM60mGq/3NHdBrjZLCUR9eDV1Skfg7L1pOjFqCkZgtNYqIiNgpqVBy++2389JLL7Fnzx527NjB97//fdavX88111wDJH52ufbaa3v2X7JkCXv37mXp0qXs3LmTBx98kFWrVnHrrbem9l0MAa2NdUR8LkLda5SU5vQ8Vjqpa9TIUYHD0UTtYRtW/BMREbFZUj/fHDx4kC996UvU19dTVFTEzJkzef755/n4xz8OQH19PbW1tT37T5o0iTVr1vDtb3+bX//611RXV3P33XePyOnAoUN7MHylWC4Tgxi5he6ex0onlgENBHIrcfgDfNDYyekVBfYVKyIiYoOkQsmqVauO+/jDDz981LaLL76YN998M6mihiOrdQ+eQBm4wOuKYDg+vNC3ewaOL7eCSIfBwYY6OLPKrlJFRERscdLXlMjARIM15IcSP9MUFJp9Hisq9wJxYk4vRnspgXpd7CoiIiOPQkkGRGNxYvED5ERLACgs67vkr9NlkuuOAJDrqyTarFAiIiIjj0JJBtS3BXGYzZhG13TgqsKj9inpasxXFCwnHPx7RusTERHJBgolGbCvJYAZ8RH2dC2cNmHUUfuUjE0EFU+8gjj78YWiGa1RRETEbgolGbC/xYcRDBP0dK1RMvrojo1lkysBCHsqMGMt1DT5MlqjiIiI3RRKMqDlYC3RTpNgTuKakvxea5R0Kx2TGCnx51bg6Ayzp2l4NCEUEZGhY/ny5Zx77rkUFBRQXl7OlVdeya5dmbvOUaEkA0JNe7B8JVgOFxAnv/joJfSLu1Z1DeaUEmvPoXnf+xmuUkRERroNGzZwww03sGnTJtatW0c0GmXhwoX4fJkZvU+6940kL96yN7FGSQHkOMM4zKOzoLfAhdMIE8WNq7OC4MG/A/MzXquIiIxczz//fJ/7Dz30EOXl5WzZsoWLLroo7cdXKMmAqL+G0lAZvgIoyDf73ccwDAq7GvPl+yvxtezMcJUiIpJOlmURiAYyflyv04thGCfesR9tbW0AlJaWprKkY1IoSbNoLI4V3483OhkfUFh29PUk3Uor8jjcESUvUsHhyHuZK1JERNIuEA1w3mPnZfy4r3/hdXJdR0+wOBHLsli6dCkXXnghM2bMSENlR9M1JWl2sCOEw9mE00pMBy6sPHqNkm7djfkMRwWWVU+bP5KRGkVERI504403sn37dh5//PGMHVMjJWm277AfR6yTsLtr4bTxZcfct+zU0bDuIIHcClyBDmqafZydW5yhSkVEJJ28Ti+vf+F1W46brJtuuolnn32WjRs3Mnbs2DRU1T+FkjTbd9iH4Q8TzOlao6Qi/5j79m7MR3uM2oMtnD2uOBNliohImhmGMaifUTLJsixuuukmnn76adavX8+kSZMyenz9fJNmhxv3E/E7ekJJQT9rlHQrHO3FIE7c9GB1lnF4v3rgiIhI5txwww08+uijPPbYYxQUFNDQ0EBDQwOBQGYu0FUoSbNQUw34CombibVJCkqOHUpM00GuOwRAjq+S8EGFEhERyZyVK1fS1tbG/Pnzqaqq6rk9+eSTGTm+fr5Js3jLXlz+UZAPHjOE6Tp+DiwpdeNrgIJQJU3t72SoShERkcTPN3bSSEmaxf17KAglLm7NzzvxPPHSsUUAuOMVBCM1tn9AREREMkWhJI1icYtobB/eSNdFrsf56aZb2eQKACLucpzRQzT7wmmtUUREJFsolKTRwfYgTrMJl9UVSioLTvic0q7ZNv7cSlx+n7oFi4jIiKFQkkb7WgIYVjtR14nXKOnW3ZgvlFOC1e5mX31DWmsUERHJFgolabS/xYcR6LVGyQBGSnLyXLiMxAwcs7OaDk0LFhGREUKhJI2aDx4g6jcGtEZJb4X5cQBy/RUED6kxn4iIjAwKJWkUOFRD3JdP1Jn4SWagoaS0MrHqa260Al+nQomIiIwMCiVpFGvtWqMEcDnCuDzmgJ5XNinxHMNRSShaSzyuacEiIjL8KZSkkdW5h/xA1xoluQMPFqWnlQMQ9JbjCRzmYEcwLfWJiIhkE4WSNInHLeKxOnK71igpKPYM+LmlVYnGfH5vOWZnmJpDnWmpUUREJJsolKRJY0cIh+swrviJuwMfqaDMi2HFiJtujI4yGg7UpatMERGRPlauXMnMmTMpLCyksLCQuXPn8uc//zkjx1YoSZN9LX4Mq5V49xol4068Rkk3h8Mg1534ycblq6Jjvy52FRGRzBg7dix33HEHmzdvZvPmzVx66aUsXryYt99+O+3HVkO+NNl32I8ZDPVMBy6qKkzq+SVlicZ8eaEKOg6rMZ+IyFBnWRZWIJDx4xpeL4Zx4t5r3a644oo+9//1X/+VlStXsmnTJqZPn57q8vpQKEmTpkP1hP0GQU/XNSVl3qSeXza2mH0NPtxWBb7AjnSUKCIiGWQFAuw6Z3bGjzvlzS0YubmDem4sFuMPf/gDPp+PuXPnpriyo+nnmzQJHNqD1ekl4k6s4lpQNrA1SrqNOr0SgKi7Aiu4n5imBYuISIbs2LGD/Px8PB4PS5Ys4emnn2batGlpP65GStIk3rIX0zcKcsFpRPB4kzvVJeNLAPB7K/B0tnGgNcC40sElXRERsZ/h9TLlzS22HDdZU6ZMYdu2bbS2trJ69Wquu+46NmzYkPZgolCSJlbnHvKDZYRzIc8bT/r5JV2N+cKeIpwdTj5obFcoEREZwgzDGPTPKJnmdrs57bTTAJgzZw5vvPEGv/rVr/jNb36T1uPq55s0iMctYtE68sKJ0Y6CInfSr+H2OnEbiQuiDF81TfveS2mNIiIiA2VZFqFQKO3H0UhJGhzqDGE6m3DHpwJQWJE3qNcpyI/R3AGeQCWd9TuBC1NYpYiIyNFuv/12Fi1axLhx4+jo6OCJJ55g/fr1PP/882k/tkJJGuxrCWAYrVhm13TgMaWDep2yygKaOyy80QqaW99KZYkiIiL9OnjwIF/60peor6+nqKiImTNn8vzzz/Pxj3887cdWKEmDfYd9GMFgzxolhdVFg3qdskmjYXcjDqMCf+D/pbJEERGRfq1atcq2Y+uakjQ4dKiRaMD4MJSMSv7KZ4CyrmnBAW8Fpu8g4WjyF8yKiIgMFQolaeA/VEOs00PYUwwkv0ZJt5KuxnwB72g8HQHqWvypKlFERCTrKJSkQfzwXpz+RK8bB1Fy8lyDep2CkhwcVgTL4cTRUcLegy2pLFNERCSrJBVKli9fzrnnnktBQQHl5eVceeWV7Nq167jPWb9+fWJu9hG3d99996QKz2qde8ntCiV5ObGkeg70ZjgMvF2N+Zz+Kg7XDeNzJiIiI15SoWTDhg3ccMMNbNq0iXXr1hGNRlm4cCE+n++Ez921axf19fU9t8mTJw+66GxmWRbxaG2vNUoGN0rSraQ08fyccCXtB9PfoVFERMQuSc2+OXKO8kMPPUR5eTlbtmzhoosuOu5zy8vLKS4uTrrAoeZQZwiH8xA5sdNoAwpHndzqfaPGl7LvoB9PvJyGdoUSEREZvk7qmpK2tjYASktPvA7HrFmzqKqqYsGCBbz44osnc9istq8lAI5WcHStUTK25KReb9TpVQBE3ZWEfDUnW56IiEjWGvQ6JZZlsXTpUi688EJmzJhxzP2qqqq4//77mT17NqFQiN/97ncsWLCA9evXH3N0JRQK9VnOtr29fbBlZty+lgBmKECoezrwmJMLJaUTS4H38edW4G5rIhiJkeMyU1CpiIhIdhl0KLnxxhvZvn07L7/88nH3mzJlClOmTOm5P3fuXOrq6rjzzjuPGUqWL1/OT37yk8GWZqvGQweJ+CGQk7jQdbBrlHQr7mrMF3Hl4+50sLfZz5TKgpOuU0REJNsM6uebm266iWeffZYXX3yRsWPHJv38888/n927dx/z8WXLltHW1tZzq6urG0yZtvA37iHW6SbsSaziOtg1Srq53CZuOgEwfRXU1TecdI0iIiIDsXz5cgzD4JZbbsnI8ZIaKbEsi5tuuomnn36a9evXM2nSpEEddOvWrVRVVR3zcY/Hg8fjGdRr2y1+eC8OXymWx8QgRm5B8h2Cj5SfH+NwJzgDVbTU7oRZw3PmkoiIZI833niD+++/n5kzZ2bsmEmNlNxwww08+uijPPbYYxQUFNDQ0EBDQwOBQKBnn2XLlnHttdf23F+xYgXPPPMMu3fv5u2332bZsmWsXr2aG2+8MXXvIpt07CE3kPjpJtcdxXAMbo2S3kZVJX6uyYlW0nZo+0m/noiIyPF0dnZyzTXX8N///d+UlJzctZHJSGqkZOXKlQDMnz+/z/aHHnqIL3/5ywDU19dTW1vb81g4HObWW29l//79eL1epk+fznPPPcfll19+cpVnIcuyiEXrKAiV0AYUFKam3+HoUyr4++5DmEY5rR1/TclriohIZlmWRTSc+R5mTrcj6UU8b7jhBj75yU/ysY99jJ///OdpquxoSf98cyIPP/xwn/u33XYbt912W1JFDVVNnWGcZiOe2HgACk7yItduZVOq4IVDBLwV0FZ74ieIiEjWiYbj3H/zhowf959/dTEuz8BnbT7xxBO8+eabvPHGG2msqn+p+Z/yAsD+1gCG2YrDSPx8U1RdnJLXLa1O/HwT9I7C3d5BZyhKvkf/6UREJLXq6uq4+eabWbt2LTk5JzdRYzD0zZZC+1r8EP5wjZKicSdeVG4gcovcOKwQccODq7OYPYc6mTG2OCWvLSIimeF0O/jnX11sy3EHasuWLTQ2NjJ79uyebbFYjI0bN3LPPfcQCoUwzfStlaVQkkIHDzURC1gEU7RGSTfDMPC6/PiiHkx/Jfv31yqUiIgMMYZhJPUzih0WLFjAjh07+my7/vrrOeOMM/jud7+b1kACCiUp5WusIdLpJOgpBqCgLDWhBKC41IWvEdzhStr2vQPnZW6KloiIjAwFBQVHrdKel5dHWVnZcVdvT5WT6n0jfcUO1+LwlWA5XECcvKKTX6OkW/n4UQC4rQqam7el7HVFRESyhUZKUsjo2Is3UAYu8LoiOMzUZb5RU6pgcw0xVwWdbZm/IlpEREam9evXZ+xYGilJEcuyiIX3UhBMXNxakJ/a393KThkNgD+3AvPwvpS+toiISDZQKEmRw74wTmcj3mhi5btUrVHSrajcC1acqDMXT4eDVn84pa8vIiJiN4WSFNnXEgBnCyZda5RUFqb09Z0uE4/RAYDLN5qaxvaUvr6IiIjdFEpSZF9LAEc4QMjTtUbJ+LKUHyMvLwaAGayicd97KX99EREROymUpEhDUzPRQIxg18JpBaNzU36Msq7GfO5oBc371JhPRESGF4WSFPEd2pNYo6R74bSy1C/PW3FqJQAORwVNLQolIiLZbiA944a6VL5HhZIUiTXvxfAXETfdgEV+SepDyagzxgAQzKkg0rw75a8vIiKp4XK5APD7/TZXkn7d77H7PZ8MrVOSKu17yfGXQQF4zAimM/V5r2RM4uLZYE4prsOHsSwr6XbUIiKSfqZpUlxcTGNjIwC5ubnD7v9fW5aF3++nsbGR4uLilCxBr1CSApZlEY/UUhgspbMACvLTMwDlLXBhxgPEHF5cvkKaOsOMLvCk5VgiInJyKisTP7l3B5Phqri4uOe9niyFkhRo8UdwOhrwRmbQCRSUpicoJBrz+eiMeXH6y9lz8DCjC6rSciwRETk5hmFQVVVFeXk5kUjE7nLSwuVypbRJn0JJCuxvCeBwtuAkMfOmsCq1a5T0VljqpPMQOCOVNO19F05TKBERyWamaaa9u+5woQtdU2Bfix8r6ifiToSS4nGlaTtWxYTEcvOueCWNDVvTdhwREZFMUyhJgfqmFuLBXmuUlOen7VjlZ4wFIOaqoKVJ04JFRGT4UChJgc7GGsI+J4GuNUoKSlM/Hbhb6anlAARyK7AO1aTtOCIiIpmmUJIC0cO10JlHzJlowleQhoXTuiUa88WImR5cbXHi8eG/MI+IiIwMCiUpYLTX4vGPAsDtiOByp++CJtN04KEtcSxfKQfaAmk7loiISCYplJwky7KwQnsoDCSuJ8nPS/8x8/ISU8vMYAU1++rTf0AREZEMUCg5SW2BCKZ5kNxo10WuJelfzKykMtGYzxWr5NCebWk/noiISCYolJykfV1rlLisrjVKKtI386Zb9eREDxyHUcHB+tfTfjwREZFMUCg5Sfta/BD3EXUlQknR+LK0H3PUtHEABHMq8R96O+3HExERyQSFkpN0oLmNeCDas0ZJYUVB2o9ZOrYIgFBOCa7mphHRGltERIY/hZKT1HlwLyGfSdDTdU1JGtco6ZaT58IZ7wTA3VHMwfZQ2o8pIiKSbgolJylyeA/4vETciRGSdK5R0ltujh8AV6CKD/YfyMgxRURE0kmh5CQZbbW4fIk1SpxGBI83Mz0OSysTc4+dkWoaP/hbRo4pIiKSTgolJ8GyLAjtpah7jZLczF3bMeb0xAwc06imvn5Txo4rIiKSLgolJ6E9EMV0NJAb6bqepNidsWOXzxgPQNBbRaBxR8aOKyIiki4KJSdhX6sfh6sFd7xr5k0auwMfqWxCMQAhTwnOphbNwBERkSFPoeQk7GsJYFmdxLrXKBlXmrFje7xOXPGuHjgdpRzq1AwcEREZ2hRKTsKB5nbigUjPdODCqsKMHj/XGwTAFajkg30NGT22iIhIqimUnIT2xj2E/WbPwmmZWKOkt7KqxM9FzugYDn7wZkaPLSIikmoKJSch0ryXeKeHsKcYyNwaJd3GnJFYbt40qjigHjgiIjLEKZScBEdbLR5fYpTENKLk5LkyevyKMxMzcPy5VQQb3srosUVERFJNoeQkGMG95AcTDfjycuIYhpHR43f3wIm4C3E2tWf02CIiIqmmUDJI7cEIpqOB/FDX9SRFmVujpJvLY+KOtwDg7iihWTNwRERkCFMoGaSmjhCm63DPGiUF5bm21JGXmwgizmAV7+/XDBwRERm6kgoly5cv59xzz6WgoIDy8nKuvPJKdu3adcLnbdiwgdmzZ5OTk8Mpp5zCfffdN+iCs0WLP4JFB5hda5SMKbGljrIxiWnIrmg1DR9staUGERGRVEgqlGzYsIEbbriBTZs2sW7dOqLRKAsXLsTn8x3zOTU1NVx++eXMmzePrVu3cvvtt/Otb32L1atXn3Txdmrr9BMLRAjkJK4pKRpTbEsdY6dNAMBwVLH/gHrgiIjI0JVUS9vnn3++z/2HHnqI8vJytmzZwkUXXdTvc+677z7Gjx/PihUrAJg6dSqbN2/mzjvv5Kqrrhpc1Vkg0FJPyO+wbY2SbuUzxsMfD+HPrSK0/xlbahAREUmFk7qmpK0tscx5aemxl1d/7bXXWLhwYZ9tl112GZs3byYSifT7nFAoRHt7e59btgm1HyLucxH2JGbAZHqNkm4l1flgxYi68nA2+22pQUREJBUGHUosy2Lp0qVceOGFzJgx45j7NTQ0UFFR0WdbRUUF0WiUpqamfp+zfPlyioqKem7jxo0bbJlpE+44hMtfimWYOIiRW5D52TcATpeJx2oFEjNwWnxhW+oQERE5WYMOJTfeeCPbt2/n8ccfP+G+R67f0d3R9ljreixbtoy2traeW11d3WDLTBu/v54CX+J6klxPDMOR2TVKesvLT4w4OYNVfHDgoG11iIiInIykrinpdtNNN/Hss8+yceNGxo4de9x9KysraWjoO1W1sbERp9NJWVlZv8/xeDx4PJ7BlJYx/kADJeFSWoH8wkGdxpQZNa6Yw7vAGaum/v2tMDn7RpZEREROJKmREsuyuPHGG3nqqaf4y1/+wqRJk074nLlz57Ju3bo+29auXcucOXNwuTK7LHsqhSLN5MS6ugOPzrO1lrHTJib+MKuo2/+arbWIiIgMVlKh5IYbbuDRRx/lscceo6CggIaGBhoaGggEAj37LFu2jGuvvbbn/pIlS9i7dy9Lly5l586dPPjgg6xatYpbb701de/CBtHIYQxH9xolxbbWUj4jMTLiz60ivH+7rbWIiIgMVlKhZOXKlbS1tTF//nyqqqp6bk8++WTPPvX19dTW1vbcnzRpEmvWrGH9+vWcffbZ/OxnP+Puu+8e0tOBAczONoJda5QU2hxKiivzMKwoMWcOzqagrbWIiIgMVlIXQ3RfoHo8Dz/88FHbLr74Yt58881kDpX1HCF/zxolhWVeW2sxTQduq4WQMRpXRwltgQhF3qH705iIiIxM6n0zCMFIDGcwQtCTWFrerjVKessviALgDFXxwf5Gm6sRERFJnkLJILT4w5hhL5bDCVacPBs6BB9p9ITEqI0zXs3+97fYXI2IiEjyFEoGodUfwYwkZtw4CeAw7T+N46YnZkJZZhX7971qczUiIiLJs//bdAhq6fTjiOQD4HL0v1R+pn04A6eS4L4dNlcjIiKSPIWSQfC3NkF3KHHHba4moXB0LkY8Qtx042zKjqAkIiKSDIWSQQi0H8KMFQCQ482OU+hwGHhoBsDdWUJHUMFERESGluz4Rh1iQm2HMGOJkZLcguxZDr+gKDFqY4aqqDlwyOZqREREkqNQMgi+zv04rcSFrvlF9i4x31v5xNFAYgbOPs3AERGRIUahZBD8wYMYJEZK8krzba7mQ2PPPAWAuLOKutqXba5GREQkOQolgxAIHQJHYoQkd1SBzdV8qHxaomOzP7eCUO1bNlcjIiKSHIWSQYgHmog6u64pGV1oczUfKijLwREPYjmcmM3ZMStIRERkoBRKBsEItBFxd4WS0uy5psQwDDzGYQBcnaX4w1GbKxIRERk4hZJBMP0+Iq5EKMnJz67GdwVFiaaJrnAlHxxosrkaERGRgVMoSVI8buEIWcQdiTDiLbC/701vFadUAOCIV1P3nmbgiIjI0KFQkqSOYBQznOgKbFgRXG7T5or6GjfzNABi7mrq9r5kczUiIiIDp1CSpO4OwZBoxpdtRk+rBiDgHU2o9m2bqxERERk4hZIktfhCGNHENGCXmX1LuecVeXDE/GA4MA9pBo6IiAwdCiVJ6mw7jCPa3YzPsrmaoxmGQY4jMQPH7SsjGInZXJGIiMjAKJQkyd/WiKOr701ObnZdT9KtoCTxrxmp4n31wBERkSFCoSRJwbZDmPHsa8bXW+VpietKHFYV+zQDR0REhgiFkiR1dh7AtBKhpKA4e/re9Db+7MkARN3V1O3ZaHM1IiIiA6NQkiSfvx4HXR2Cy7Kn701vo6dUARD0jiK4Z6fN1YiIiAyMQkmS/MFGLEcijOSOyp6+N715C9yYsQ4AzCbD5mpEREQGRqEkSXFfExFXV4fg8iKbqzk2j6MFAJe/lFBUM3BERCT7KZQkyfK39PS9yS3JtbmaYysqS/ynNcNV1NSrB46IiGQ/hZIkmX5/z0hJtjXj661y8lgAHFSzVzNwRERkCFAoSZIZMsBInLZsDiUTzjkdgIinin0fbLC5GhERkRNTKElCKBrDDCXWJnHEg5hm9p6+UZMT3YJDOSUEa3bZXI2IiMiJZe+3ahZq80dwRBLXkTiNoM3VHJ8n14Uz2gqA2aT/zCIikv30bZWEFn8EsyuUZGMzviN5nK0AOP2lhKNqziciItlNoSQJbe1tGLHEGiVuT/Y14ztS0ahEbx5ntJq9B5ttrkZEROT4FEqS4G/9sBmfJ89pczUnVj1lPAAGVezZvdnmakRERI5PoSQJ/rZDOLqa8eUV5NhczYmNP+cMAEI5VdS+t97eYkRERE5AoSQJnR0HMEmEkvyS7GzG19uo08oBiLgLCdbstrkaERGR41MoSUJn5wGMrmZ8BWXZu8R8N5fHxBVNXEtiNps2VyMiInJ8CiVJ8AUOEje7fr4ZnZ3N+I7kdrYB4PKXEYlpBo6IiGQvhZIkxDobiXb3vakotreYASouT6w6a0arqD142OZqREREjk2hJBm+wx92CC722lzMwIyZOjHxh1HFnvc0A0dERLKXQkkSDH+AqDOxeFo2973pbfzsqUBiBs6ed/9iczUiIiLHlv2LbWQRM9CV4aw4Hu/QOHVlk8rAihN15RGr+cDuckRERI4p6ZGSjRs3csUVV1BdXY1hGDzzzDPH3X/9+vUYhnHU7d133x1szbawLAtHyA2AafkxHIbNFQ2M02XiijYB4DisGTgiIpK9kg4lPp+Ps846i3vuuSep5+3atYv6+vqe2+TJk5M9tK06QlHMSOI6EifZ3YzvSB53O5CYgROLZ//y+CIiMjIl/RvEokWLWLRoUdIHKi8vp7i4OOnnZYtWXwRHNA9c4HJmfzO+3koqPHQ2gBmrpvZQC5MqSu0uSURE5CgZu9B11qxZVFVVsWDBAl588cVMHTZlWjr9OKKJ6cDunKHx0023sdNPTfzhqOKDv79ubzEiIiLHkPZQUlVVxf3338/q1at56qmnmDJlCgsWLGDjxo3HfE4oFKK9vb3PzW6+tiYc8USHYE/e0Jh50238nMQMnGBOFXvfHXqBUERERoa0TyGZMmUKU6ZM6bk/d+5c6urquPPOO7nooov6fc7y5cv5yU9+ku7SkuJvPYjD6lrNtXBorFHSrWR8CUY8SsyZQ+j9vXaXIyIi0i9b1ik5//zz2b372A3ili1bRltbW8+trq4ug9X1r7PtAA6rq+9NaYHN1STHNB24Yt0zcIbGVGYRERl5bPmG2rp1K1VVVcd83OPx4PF4MljRibV17sdhJEZKCkdlfzO+I3nc7YStSpyBUuJxC8cQmdIsIiIjR9KhpLOzk/fee6/nfk1NDdu2baO0tJTx48ezbNky9u/fzyOPPALAihUrmDhxItOnTyccDvPoo4+yevVqVq9enbp3kQF+XwO55gwAcsuHXigprsqh40BiBs6+plbGl5fYXZKIiEgfSYeSzZs3c8kll/TcX7p0KQDXXXcdDz/8MPX19dTW1vY8Hg6HufXWW9m/fz9er5fp06fz3HPPcfnll6eg/MwJtzfgcs0FIK+82N5iBmH8madTdyCIZVbx3q7XGV/+CbtLEhER6SPpUDJ//nws69gLcD388MN97t92223cdtttSReWdTo+bMbnLcqxuZjkjZ8zjVdeeJNgTiV73/l/ME+hREREsosa8g2QEQgQNxPLzA+VZny9FY8pwoiHiZtuQh/ss7scERGRoyiUDFB3Mz4jHsHlGXo9ZBwOA1fsEACGZuCIiEgWUigZoD7N+IyhOXPF4+kAwBkcRVw9cEREJMsolAxAJBbHDOcC4DSGVjO+3krHJN6DGaviwOE2m6sRERHpS6FkAFr9ERyxxEWuLmfU5moGb8JZZwAQdyZm4IiIiGQThZIBaPWFMLpCiXvoTbzpMX7ONACCORXs3fE/NlcjIiLSl0LJALS3HcYRS6zmOhRn3nQrrCjAEQtiOZwEPthvdzkiIiJ9KJQMgK+1sacZX25xns3VDJ5hGLhijQA4DrttrkZERKQvhZIB6Gjdj0FX35vSQpurOTlubycAztCo4y6CJyIikmkKJQPQ1rEfHF2hZPTQ7hlTNjbxPox4FQ0t7TZXIyIi8iGFkgHo6DxA3Ex8medVFNtbzEmaePZ0AOJmFbt3brK5GhERkQ8plAxAtO0A4a6+N0OxQ3BvE86dCkDQW86e7f/P5mpEREQ+pFAyAPG2ZqLdoWQINuPrLa8sDzPqA8NB4P0DdpcjIiLSQ6FkIHxhLCPR7yYnb+hOCYauGThWVw+cFpcudhURkayhUDIADn+i140jFsB0Dv1TllsQBsAZKueD+iabqxEREUkY+t+wGWBGPIl/rYDNlaTGxLNOB8CMn8obr622uRoREZEEhZITsCwLRzhxHYnpGLrN+HqbtuhcAHwFkzj05nM2VyMiIpKgUHICvnAMM9rdjC9iczWpUVRZgDt2CMswyTlQQDASs7skERERhZITafWHMeKJUOLxGjZXkzqjKhPdjl3Bqbz59nabqxEREVEoOaG29g6MnmZ8w6dfzJkL5wAQyZnG9lcfsrkaERERhZIT6mw92NP3Jrc43+ZqUmfi+adhxMOEPSXE3/nA7nJEREQUSk6kvaX+w2Z8o4b2aq69OV0mec56AHJaJ9HQ0mlzRSIiMtIplJxAa3sdOBLXlBSXl9lcTWpNOqsKADM+nb9u+qPN1YiIyEinUHICHa37iA2TZnxHOvNTcwHw5Z9K3etar0REROylUHICodYDRFxd15QMs1BSMqYIV7QJy2HiqnMRi2vJeRERsY9CyQnE2pqJunIByC3w2FxN6pWN9gPg9k/lrfd0wauIiNhHoeQEjM5EnxisOO5cp73FpMGMjyWmBkc903jzZU0NFhER+yiUnIDR1YzPjPtxOIbP4mndTpk3FSMeIZRTiv9vO+wuR0RERjCFkhMwwi4ATMtvcyXp4XKb5Dr2A5BzeBztgbDNFYmIyEilUHICzogXANMI2VxJ+oyfnpjq7IxN4/U3/sfmakREZKRSKDmOaCyOEU2EEqdreDTj689Z/3ARkJgavPvl39tcjYiIjFQKJcfRFojg6Op7M5ya8R2pdHwxrmgzlsOF+UEcy9LUYBERyTyFkuNo6QxgWIlQ4h2G04G7GYZBUXEbAG7/FN7f32hzRSIiMhIplBxHZ0tjTyjJLS2wuZr0mrHgHABi7mm8vvG3NlcjIiIjkULJcbS31mMYib43RaNKbK4mvSZfMrNravAo2jZvsrscEREZgRRKjuNwax2WIzFSUlwxvJrxHcmd4ySHfQB4mioJRmI2VyQiIiONQslxtB/eS9SZCCX5laU2V5N+Y6ckfqIyo9PY8rc3bK5GRERGGoWS4wi19GrGV15kczXpN2vxJQAE8ifz1voHbK5GRERGGoWS44i0HCZuugHwFrhtrib9Rp06Cmf0MHGHC/7us7scEREZYRRKjqerGZ8Rj+LymDYXk36GYVCY3wyA2zeZgy3tNlckIiIjiULJ8XS1uzHjnRjG8F08rbdpl8wEIO6czisbn7C5GhERGUmSDiUbN27kiiuuoLq6GsMweOaZZ074nA0bNjB79mxycnI45ZRTuO+++wZTa8aZYScAjmHajK8/ZyycgxGPEvSO5uDL6oMjIiKZk3Qo8fl8nHXWWdxzzz0D2r+mpobLL7+cefPmsXXrVm6//Xa+9a1vsXr16qSLzTRHdzM+x/Btxnckj9eFJ14HgLuxjFhcS86LiEhmOJN9wqJFi1i0aNGA97/vvvsYP348K1asAGDq1Kls3ryZO++8k6uuuirZw2eUEc0FNzhdUbtLyaiqUz3U7E1MDX7r7+9y1hlT7S5JRERGgLRfU/Laa6+xcOHCPtsuu+wyNm/eTCTSf+fdUChEe3t7n1umBcIxHLHEaq7DuRlff2b/4wIAArmT2bx2aPzUJiIiQ1/aQ0lDQwMVFRV9tlVUVBCNRmlqaur3OcuXL6eoqKjnNm7cuHSXeZQWXwiHlQgl3sLhPx24t/IplTijrcRNN5F3+v9vJCIikmoZmX1z5MwVy7L63d5t2bJltLW19dzq6urSXuOR2lqbgMTCaXmlw3/htN4MwyDP2wCAp2MS7YGRc02NiIjYJ+2hpLKykoaGhj7bGhsbcTqdlJX130/G4/FQWFjY55ZpHS0NYCRCSVH58F9i/khnzJsGJKYGv/ryszZXIyIiI0HaQ8ncuXNZt25dn21r165lzpw5uFyudB9+0Jpaaol3NeMrqRxtczWZN+OTF2DEYwS9FexZ/3/tLkdEREaApENJZ2cn27ZtY9u2bUBiyu+2bduora0FEj+9XHvttT37L1myhL1797J06VJ27tzJgw8+yKpVq7j11ltT8w7SpLW5hlhX35v8qlE2V5N5OXluPLHEf1Pn/vyen9xERETSJelQsnnzZmbNmsWsWbMAWLp0KbNmzeKHP/whAPX19T0BBWDSpEmsWbOG9evXc/bZZ/Ozn/2Mu+++O+unAwcP7SPiSlzomjs68z8fZYPR4xP/OiNTeX/fPnuLERGRYS/pdUrmz59/3P/V/PDDDx+17eKLL+bNN99M9lC2irS2YBmJfje5+SNr9k232Z9eQN2vPyCQezqvvfAbTvvaz+0uSUREhjH1vjmGeHtixokjFsR0jczTVD1jAs5oG3HTg//NPXaXIyIiw9zI/LYdACOQmK5sxn02V2IfwzDwuvcD4GqfQCgysla2FRGRzFIoOQYj2NWMj5HTjK8/k+eeAoDlmMZfN79oczUiIjKcKZQcgyOSk/h3BDXj68/Ziy8BK04wt4p3n3/c7nJERGQYUyg5Bkcs0SHY6eq/P89I4S304okkZlMZtSPzgl8REckMhZJ+xOIWRlczPvcIa8bXn5KqxGiRGZ7CwabDNlcjIiLDlUJJP9oDEQwrsXBabqHH5mrsN+fTHwMg5J3CxhdW2VyNiIgMVwol/WjraO/pe1NQNrKa8fVn/KzTMKPtxJw5tGx6y+5yRERkmFIo6Ufb4Q+b8RVUlttcjf0Mh0GOs2vJ+dZq4nEtOS8iIqmnUNKPpuYa4mbimpLSqpHXjK8/k2aPBcAypvG3t7fYXI2IiAxHCiX9aGncQ9SZGCkpUigB4NyrFnZNDR7Dm398yO5yRERkGFIo6YfvUC3RrmZ83tEFNleTHXJL8/FE6gCIvx+zuRoRERmOFEr6EWluTfxhxcnJ09oc3QpGdQDgDE2mwx+wuRoRERluFEr6EetqxmfG/DgcWqek25zFlwAQyjmD9Wt/Z3M1IiIy3CiU9MefmF3iGMHN+PpzytwZmNEOYk4vDRs22V2OiIgMMwol/XCEzMS/I7wZ35EMh4HH2AOAo7kcy9LUYBERSR2Fkn4YasZ3TGNnlgFgMJW3dr9rczUiIjKcKJT0w+hpxhe1uZLsM/fzi8GKE8gdx0sP/cjuckREZBhRKOmHI56YDuzK1UWuR8ofVYA3vheA3L3T2FVTY3NFIiIyXCiUHCEYiYGa8R3XvC99BAB//gX8z/3ftbkaEREZLhRKjtDaGQASoaRwlJrx9WfyRWfije3BcjjJ2zuD3Xs0WiIiIidPoeQIrYcbsByJn2+KqypsriZ7XXRt92jJXNbd/z2bqxERkeFAoeQIh5r2EDcTIyWl1dU2V5O9Tps3A2+sBsthkrtnBu/t3Wt3SSIiMsQplBzh8MH3ibi6mvGNUTO+47nouvMB8Oefz7r7NFoiIiInR6HkCJ31tcTNxAWuuaPUjO94TrtwOt7YB1gOk5y90/mgttbukkREZAhTKDlCqLkFACMexZVj2lxN9rv4y3MBCBScxwsrNVoiIiKDp1ByhGhPM75ODEPrlJzIqR+djjf2PpaRGC2p2Vdnd0kiIjJEKZQcyadmfMma/5ULAAgUfIQX7tVoiYiIDI5CyZG6m/EZAZsLGTpOmTud3Nh7PaMle/fvt7skEREZghRKjuCIJC5yNRxBmysZWi7+ykcB8OV/hOfv0SqvIiKSPIWSI3Q34zPVjC8pidGS3WA4cO+dQe2BA3aXJCIiQ4xCyREMNeMbtPlfmweAv2AOf9ZoiYiIJEmhpJd43MJQM75Bm3TeNHJjf+8ZLanTaImIiCRBoaSXjmAEjMRIiZrxDc4lX78IAH/+bP78X8tsrkZERIYShZJeWluaiDsSIyUlY6psrmZomviRaeTGdoHhwFU7nX0N9XaXJCIiQ4RCSS+HDn3YjK9s3Dibqxm6Lvn6xQD488/hubs1WiIiIgOjUNJL44FdPc34iseU21zN0JUYLXm369qS6Rw42GB3SSIiMgQolPTSsb8Oy5FYPC23LN/maoa2S74xHwB//iz++CuNloiIyIkplPQSbEo043PEgjhdasZ3MibOmUZubGfPaEl9o0ZLRETk+BRKeom0JZaWN2Pqe5MKl35jPlhx/AXn8Kf//L7d5YiISJYbVCi59957mTRpEjk5OcyePZuXXnrpmPuuX78ewzCOur377ruDLjpdLF8cAIelUJIKE+ZMJze+EwBn7TTqDx20uSIREclmSYeSJ598kltuuYXvf//7bN26lXnz5rFo0SJqa2uP+7xdu3ZRX1/fc5s8efKgi06brmZ8BmrGlyqXLlnQNVoyi+fu0miJiIgcW9Kh5K677uKrX/0qX/va15g6dSorVqxg3LhxrFy58rjPKy8vp7Kysudmmtl3zYYRyUn8a6oZX6pMmD2tZ7TErJvGwaZGmysSEZFslVQoCYfDbNmyhYULF/bZvnDhQl599dXjPnfWrFlUVVWxYMECXnzxxePuGwqFaG9v73PLBCOWCCWmK5KR440UC/6la7Qk/2z+9EuNloiISP+SCiVNTU3EYjEqKir6bK+oqKChof/ZFVVVVdx///2sXr2ap556iilTprBgwQI2btx4zOMsX76coqKintu4DC1k1tOMz6tmfKk0/pxp5MbfAcBRO53Gw002VyQiItloUBe6GkbfL23Lso7a1m3KlCl8/etf55xzzmHu3Lnce++9fPKTn+TOO+885usvW7aMtra2nltdXd1gyhyExNokHjXjS7kFN3ys69qSmfzxF7fbXY6IiGShpELJqFGjME3zqFGRxsbGo0ZPjuf8889n9+7dx3zc4/FQWFjY55Zu4Wi8VzO+9B9vpBl/9jRy428D4Ng3naaWZpsrEhGRbJNUKHG73cyePZt169b12b5u3TouuOCCAb/O1q1bqarKroZ3La2tPc34Ro2ttrma4eljN3y869qSM3n637XKq4iI9OVM9glLly7lS1/6EnPmzGHu3Lncf//91NbWsmTJEiDx08v+/ft55JFHAFixYgUTJ05k+vTphMNhHn30UVavXs3q1atT+05O0qHGPcSciZGSURMn2VzN8DTu7GnkxtfiN2di7p/L6sd/xVWfv9nuskREJEskHUquvvpqmpub+elPf0p9fT0zZsxgzZo1TJgwAYD6+vo+a5aEw2FuvfVW9u/fj9frZfr06Tz33HNcfvnlqXsXKXCwbhdR52gASsdppCRdFt16Jc/8xw4CeRMw/rSdDVV/5OL5V9hdloiIZAHDsizL7iJOpL29naKiItra2tJ2fcnqVXfQ8MZHwIrzL/deisPUCvzp8vbaV9nwf9qxHG5yO17ivB/+I9OmzLS7LBERSbFkv7/1zdsl0HgYADPmVyBJs+kLL+DMuYGu2TjzeO3nD2oJehERUSjpFm5NLC3vUDO+jJh3/T8yfsJeAIJ5/8Cfvnc7/oBW0hURGckUSrrEA2rGl2lX3P5VSnMTi6pFzM/xu//1dWKxuM1ViYiIXRRKugW7mvEZfpsLGVk+d+cN5Fk7sRwmVvBzrPrBN+wuSUREbKJQ0i3iBsAwQzYXMrIYDoMv3v0NcsIfEHN64eAiHrnru3aXJSIiNlAo6WLEcgEwXVGbKxl5nB4nn7vzs7hD9UTcxUT+dibPPHq33WWJiEiGKZR0MayuUOLN+hnSw1JeaSFXLPsozkgrQW81h5938/KGP9pdloiIZJBCSY/Eaq45BWrGZ5fK0ycy/7qxOGIBAvmns/u/d7Jr1w67yxIRkQxRKCHR5dgyEn1vCkYX2VzNyDbloo9w1qURjHgMf/4cXvn54xw6fMjuskREJAMUSoAOf4C4IzFSUqpmfLa74PNXMuH0xBomgbyP8X+/8xOCIV2ALCIy3CmUAI0H9xJzJkZKxpxyms3VCMAnb/0ao4q2AxByX8nvvn0D8biu9xERGc4USoD6mneIOXMAKJkw1uZqpNtn77iZfLaB4SAa+QwPLvsXu0sSEZE0UigBDtcmfiow4jFyCr02VyPdDMPgi/fcjDe8k7jpJtZ0Ob//xe12lyUiImmiUAL4DjYDYMY6MQzD5mqkN9Np8oX//DKe4F6irnz875zNHx+9x+6yREQkDRRK6N2Mr9PmSqQ/OUV5/ONPP44rdIhwziga1+bxv3/zb1iWrjERERlOFEqAmC/RBM5AzfiyVdn4sVy6ZDLOSAfB3Am0/PVMHvjGEvY37LO7NBERSRGFEoBQ4jQ4CNhciBzPaeedwyVfqcYdqiPqyiPsuJq1t/4fVj/wC7tLExGRFFAoAYh0reJqBu2tQ07o9I+ey/X//UVK87dhxKP482fS9NpUHvjGEhoa6+0uT0REToJCCWDEEtOBHWrGNyQ4nSafv3Mp87+Q3zNqEjI+y5qlT/LMQ/9pd3kiIjJICiUAVmI1VzNXF04OJdPmX8D1919DSe6bGPEogfyZHHz5NB5Y8k0ONTfaXZ6IiCRJoQT4sBmf2+Y6JFlOl5Mv3HUrF3/O++GoCZ/hjzc/xrO/09RhEZGhRKEEsLr63uSOKrS5Ehms6ZfO48u/+QLF3i09oyb16yfwwJJvcri12e7yRERkAEZ8KIlEY8TMRN+bMjXjG9JcbhfX/Of/Yt5nPX1GTZ658Xesefw+u8sTEZETGPGh5ODBOqJdzfiqT5tsczWSCmd+7GK+/JsvUJTzRs+oSd3/jOGBf7mR1rbDdpcnIiLHMOJDyb6/v4XlcAJQecqpNlcjqeJyu/jiiu9ywWeceIJdoybWp3nqxkd4ZtUvicfidpcoIiJHGPGhpLl2DwCOWAh3Xo69xUjKnb3wUq67/wsUeV5PjJrkzeTAX2fy8JdX8NBt32L/gb12lygiIl1GfCjpbPiwGZ8MTy63iy/+ahkXXGXgCe7GMkwCeWfjb7+SNbdvYtVXb2Xd6ofUS0dExGZOuwuwW6jFD4Ajrr43w93Zl32csy+Dt158iTcee5mQdSbhnNHA5ex+IUrdM/+BObGRK779XUpLR9tdrojIiDPiQ0nMHwPAsBRKRooZl8xjxiXzCLR18sc7V9JeV0YoZyKBvDlwCJ5aug7DtZkzPnsRH/34lXaXKyIyYoz4UBIPmokfsQy/3aVIhnmL8vnsz/4XAG/83+fY8exbhBwzCeVUAp9i+/8Os+uRfyPnDB+Ll95OnjfP3oJFRIa5ER9KiLjAA4aa8Y1o5y7+JOcu/iTtjU386c778TWOI5wzhkDe+QTq4Ikl/we82zn3+iuYed58u8sVERmWRnwoMeLexL/uiM2VSDYoLB/FF35xO5ZlsfH3T/DeX/YRNs8k6B0HjOOVBwJsvuceLE8t+ZNzmffFL1M9ZqLdZYuIDAsjPpRg5QJgajaw9GIYBhd/8fNc/EU4tGcvz6/4LcG20wh7KgnkTQOmEayBZ3/0Nu7Q/8Xy7qNoWgmXfvGrlI6qsLt8EZEhacSHEovEaq7uApfNlUi2Gj1xAl9a8UPi8ThvrlnDO/+zhUhrMRHnqcScuQScZwJnEnwX/vDd13GFd0NePWWzqllwzdfIz1dPJRGRgVAo6WrGlzeqyOZKJNs5HA7mfOpTzPnUpwCIRmO8/tRq3nvpHaLtowi7TyHqyifqmgXMYt9WeOz1/4czuhsKGqn6yGnMv/o6vDlee9+IiEiWGtGhxLKsnmZ8xdWVNlcjQ43TafLRz36Wj342cT8cDPLS409Q+8Ye4r4Kwu5JRNxFRNxzIAYfvAa1G9fhjNSD0YiR20FOVQ4TPnIOsy/5BDke/YYoIiPbiA4lzU0HiToT15SMmXy6zdXIUOfOyWHB9V+G6xP3gx0d/OV3v+fg9gZigWrCnoldIymTgclgQeAAtDwDb/1hLc5IPYbRCF1hZdJHZjP70kW43W4b35WISOaM6FBSt3MHGCYAY08/w+ZqZLjJKSjg8m8u6bnf2dLC5ufXcGD7+wSbLKxwMXFHJWF32THDyvY/PI8z0gDGQYzcTlzFDvKrS6meOpVp582joEA/O4rI8DGiQ0ljzR7gVMyoD7dXQ+eSXvklJcz//DXw+b7bO1ta2Pzn59i//T1CzWCFS4g7Kgi7R3WFldOA0xJhpQXaW+DA27Dlf2/CFWnFEWsFow3MDoycMK4iJ4VjyqieNpVpH5lHXn6+HW9XRCRpIzqUtB9sAk5VMz6xVX5JCfO/8EX4Qt/tnS2H+eua56jf8X5XWCnCooiYWUzUVYjlcBH2jAZ69emJfxhc9r0FbzzxyofBhQ5wBMAZwnBHMb0G7kI33lHFlFZXMWbKGYw9bSoetyeD715E5EMjOpSEWhJhRM34JBvll5Ry6TVf6vexkK+TXW/8ldq33qb9wGFCrVGsoBsrln/i4AIQB3xdt3rYvwN2vHAYI74eZ9SHI+br6gflA0cAwwyBM4bDZWF4DJxeJ+78HHKK88kvK6OkopLyCRMYXTUBl1vT60VkcAYVSu69917+4z/+g/r6eqZPn86KFSuYN2/eMfffsGEDS5cu5e2336a6uprbbruNJUuWHHP/TIn64l1/KZTI0OLJy2fm/EuZOf/SY+4T7Oxk5183se+dnbTvaybSGSUeckDUjRX3ArlYRh5xRx4xM4+46cZyuIi4i4Hi/l802nXzAU1HPliLEa/BEQ/iiAVxWEEMKwhWCAiBEcUwIlhGDMOMgiOO4YxjOMFwgukxMT1OXLke3LlePPl5eIsKKSgrI6+giLySIgpKRpOfX4jDNE/6HIpI9kk6lDz55JPccsst3HvvvXz0ox/lN7/5DYsWLeKdd95h/PjxR+1fU1PD5Zdfzte//nUeffRRXnnlFb75zW8yevRorrrqqpS8icGKB41EMz7UjE+Gn5z8fGZd+jFmXfqxAe3f0tDA3nfeonHvHtoPNhNq9RH1xYgFu4KM5QLLA5YHy8hJ3Bwe4o4cYs7E2iuWwyTmyCPmHGDzwjgQ7rod9/8MO7pu+8CK4YhHMawoRjyKw4qAlbgP3f/G6E5QBjEwYljEMYwYGHEsI45BHBwWGBaGwwJHHBxgOMBwGImw5DAwnA4M08AwHZgOBw6XA8M0MZ0mDpcT0+nCdDkx3S5cbhem243p8uDO8eDK8WC63LjcnsQtx4Pb48blycHtycPjzcGtqeAiPQzLsqxknnDeeedxzjnnsHLlyp5tU6dO5corr2T58uVH7f/d736XZ599lp07d/ZsW7JkCX/729947bXXBnTM9vZ2ioqKaGtro7AwdatjrvrKjwm6L8Ib2cBXVv0kZa8rMtJEQ2GaGvZxsHYPrQ0H6Ww+TKC1g1BnkFggSjxsYUXBijog7oC4E8tyguUEywW4wHBh4cYy3FiGC8vhJu5wE3e4EklhuLLiGJYFxDGsxA2srn/jYFmJAIXV9bfVsx36Pkb330dt732jz34fbut+fuJxo2uLYST2sbC6tnW/fq+/DavnFYxer2kZif16Xt+g53V6b+vz3A8f7PP4h3vRte3Dry7DAMvo9ZDR92utz2O9Gb03Whj0t2PvGoxEiMVIvA+j145Hvnifeuh6d93PsXp26L7Xs79xdA1GP3/h6Ofx3jUYfe4cvbn363dtnDL/As65dCGplOz3d1IjJeFwmC1btvC9732vz/aFCxfy6quv9vuc1157jYUL+77Jyy67jFWrVhGJRHC5jv79ORQKEQqFeu63t7cnU+bAxdSMTyQVnB43lRNOoXLCKSl/bcuyCHR20NbcRGdrC76ONoLt7QR9nYR8AcKBEJFAkGgwRCwSIx6OEY/GsCIW8ShYcQviBsQNrHjXt5dlgOXAshxgOcAyMXBgYZL4//YmWGbXkgEOLBxd27uGUjAT24yu5xgOLKNrP8NM/G04sAzzxIHKcHR9aZpYx99TBurIE6kTOyBvv/CXlIeSZCUVSpqamojFYlRU9G04VlFRQUNDQ7/PaWho6Hf/aDRKU1MTVVVVRz1n+fLl/OQn6R+5MJwBPMF6XBX6fVokWxmGQW5BIbkFQ7OHUCwWIxIKEQ4GCAUDRMIhosEQ4VCISCREJBQiHg4TjkSIRcLEIpGuW5RYLEY8GiEeixOPx4lHYsTjMax4HOIWsVgMK2ZhxeO9biT+jXWNcFiAZZEYgLGwLAviiTGODwdcrJ4vbqt7gKX7DRw5mAJgGUc8bvR6Qq//CW71/rvXY73+tnqNGPTev2d7r2EOy+q1b9dfFh8e2zCMxPvps0evWnvGQvpu6zVw0XX/iNEKq/9Rh6PvH2e4xDrO/aNGaI54ncSQ1dHb6Gf7MWvpf3vvcSvPOPtbYAzqQlfjiJNgWdZR2060f3/buy1btoylS5f23G9vb2fcuHGDKfW4vrLqxyl/TRGR3kzTxMzNJSc31+5SRLJeUqFk1KhRmKZ51KhIY2PjUaMh3SorK/vd3+l0UlZW1u9zPB4PHo/WShARERlJkrp6zO12M3v2bNatW9dn+7p167jgggv6fc7cuXOP2n/t2rXMmTOn3+tJREREZGRK+pL2pUuX8sADD/Dggw+yc+dOvv3tb1NbW9uz7siyZcu49tpre/ZfsmQJe/fuZenSpezcuZMHH3yQVatWceutt6buXYiIiMiQl/Q1JVdffTXNzc389Kc/pb6+nhkzZrBmzRomTJgAQH19PbW1tT37T5o0iTVr1vDtb3+bX//611RXV3P33XfbvkaJiIiIZJek1ymxQ7rWKREREZH0Sfb7exivSCQiIiJDiUKJiIiIZAWFEhEREckKCiUiIiKSFRRKREREJCsolIiIiEhWUCgRERGRrKBQIiIiIllBoURERESyQtLLzNuhe9HZ9vZ2mysRERGRger+3h7o4vFDIpR0dHQAMG7cOJsrERERkWR1dHRQVFR0wv2GRO+beDzOgQMHKCgowDCMlL1ue3s748aNo66uTj11kqDzNjg6b4Oj85Y8nbPB0XkbnOOdN8uy6OjooLq6GofjxFeMDImREofDwdixY9P2+oWFhfoADoLO2+DovA2OzlvydM4GR+dtcI513gYyQtJNF7qKiIhIVlAoERERkawwokOJx+PhRz/6ER6Px+5ShhSdt8HReRscnbfk6ZwNjs7b4KTyvA2JC11FRERk+BvRIyUiIiKSPRRKREREJCsolIiIiEhWGNGh5N5772XSpEnk5OQwe/ZsXnrpJbtLymo//vGPMQyjz62ystLusrLOxo0bueKKK6iursYwDJ555pk+j1uWxY9//GOqq6vxer3Mnz+ft99+255is8SJztmXv/zloz57559/vj3FZonly5dz7rnnUlBQQHl5OVdeeSW7du3qs48+a0cbyHnT5+1oK1euZObMmT1rkcydO5c///nPPY+n6rM2YkPJk08+yS233ML3v/99tm7dyrx581i0aBG1tbV2l5bVpk+fTn19fc9tx44ddpeUdXw+H2eddRb33HNPv4//4he/4K677uKee+7hjTfeoLKyko9//OM97RRGohOdM4BPfOITfT57a9asyWCF2WfDhg3ccMMNbNq0iXXr1hGNRlm4cCE+n69nH33WjjaQ8wb6vB1p7Nix3HHHHWzevJnNmzdz6aWXsnjx4p7gkbLPmjVCfeQjH7GWLFnSZ9sZZ5xhfe9737Opouz3ox/9yDrrrLPsLmNIAaynn3665348HrcqKyutO+64o2dbMBi0ioqKrPvuu8+GCrPPkefMsizruuuusxYvXmxLPUNFY2OjBVgbNmywLEuftYE68rxZlj5vA1VSUmI98MADKf2sjciRknA4zJYtW1i4cGGf7QsXLuTVV1+1qaqhYffu3VRXVzNp0iQ+97nP8cEHH9hd0pBSU1NDQ0NDn8+ex+Ph4osv1mfvBNavX095eTmnn346X//612lsbLS7pKzS1tYGQGlpKaDP2kAded666fN2bLFYjCeeeAKfz8fcuXNT+lkbkaGkqamJWCxGRUVFn+0VFRU0NDTYVFX2O++883jkkUd44YUX+O///m8aGhq44IILaG5utru0IaP786XPXnIWLVrE73//e/7yl7/wy1/+kjfeeINLL72UUChkd2lZwbIsli5dyoUXXsiMGTMAfdYGor/zBvq8HcuOHTvIz8/H4/GwZMkSnn76aaZNm5bSz9qQaMiXLkd2HLYsK6VdiIebRYsW9fx95plnMnfuXE499VR++9vfsnTpUhsrG3r02UvO1Vdf3fP3jBkzmDNnDhMmTOC5557j05/+tI2VZYcbb7yR7du38/LLLx/1mD5rx3as86bPW/+mTJnCtm3baG1tZfXq1Vx33XVs2LCh5/FUfNZG5EjJqFGjME3zqATX2Nh4VNKTY8vLy+PMM89k9+7ddpcyZHTPVtJn7+RUVVUxYcIEffaAm266iWeffZYXX3yxTzd1fdaO71jnrT/6vCW43W5OO+005syZw/LlyznrrLP41a9+ldLP2ogMJW63m9mzZ7Nu3bo+29etW8cFF1xgU1VDTygUYufOnVRVVdldypAxadIkKisr+3z2wuEwGzZs0GcvCc3NzdTV1Y3oz55lWdx444089dRT/OUvf2HSpEl9HtdnrX8nOm/90eetf5ZlEQqFUvtZS9FFuEPOE088YblcLmvVqlXWO++8Y91yyy1WXl6etWfPHrtLy1rf+c53rPXr11sffPCBtWnTJutTn/qUVVBQoHN2hI6ODmvr1q3W1q1bLcC66667rK1bt1p79+61LMuy7rjjDquoqMh66qmnrB07dlif//znraqqKqu9vd3myu1zvHPW0dFhfec737FeffVVq6amxnrxxRetuXPnWmPGjBnR5+xf/uVfrKKiImv9+vVWfX19z83v9/fso8/a0U503vR569+yZcusjRs3WjU1Ndb27dut22+/3XI4HNbatWsty0rdZ23EhhLLsqxf//rX1oQJEyy3222dc845faaEydGuvvpqq6qqynK5XFZ1dbX16U9/2nr77bftLivrvPjiixZw1O26666zLCsxVfNHP/qRVVlZaXk8Huuiiy6yduzYYW/RNjveOfP7/dbChQut0aNHWy6Xyxo/frx13XXXWbW1tXaXbav+zhdgPfTQQz376LN2tBOdN33e+veVr3yl5/ty9OjR1oIFC3oCiWWl7rOmLsEiIiKSFUbkNSUiIiKSfRRKREREJCsolIiIiEhWUCgRERGRrKBQIiIiIllBoURERESygkKJiIiIZAWFEhEREckKCiUiw8j8+fO55ZZbMna8L3/5y1x55ZUZO95gZfq8iMjgKJSIiIhIVlAoEREZhFgsRjwet7sMkWFFoURkmAqHw9x2222MGTOGvLw8zjvvPNavXw9AW1sbXq+X559/vs9znnrqKfLy8ujs7ARg//79XH311ZSUlFBWVsbixYvZs2fPoOqZP38+3/rWt7jtttsoLS2lsrKSH//4xz2P79mzB8Mw2LZtW8+21tZWDMPoqXv9+vUYhsELL7zArFmz8Hq9XHrppTQ2NvLnP/+ZqVOnUlhYyOc//3n8fn+f40ejUW688UaKi4spKyvjBz/4Ab1bfx3vfAE8/PDDFBcX86c//Ylp06bh8XjYu3fvoM6FiPRPoURkmLr++ut55ZVXeOKJJ9i+fTv/9E//xCc+8Ql2795NUVERn/zkJ/n973/f5zmPPfYYixcvJj8/H7/fzyWXXEJ+fj4bN27k5ZdfJj8/n0984hOEw+FB1fTb3/6WvLw8Xn/9dX7xi1/w05/+lHXr1iX9Oj/+8Y+55557ePXVV6mrq+Ozn/0sK1as4LHHHuO5555j3bp1/Nd//ddRx3Y6nbz++uvcfffd/Od//icPPPDAgM5XN7/fz/Lly3nggQd4++23KS8vH9R5EJFjSFlfYxGx3cUXX2zdfPPN1nvvvWcZhmHt37+/z+MLFiywli1bZlmWZT311FNWfn6+5fP5LMuyrLa2NisnJ8d67rnnLMuyrFWrVllTpkyx4vF4z/NDoZDl9XqtF154wbIsy7ruuuusxYsXD7i2Cy+8sM+2c8891/rud79rWZZl1dTUWIC1devWnsdbWloswHrxxRcty7KsF1980QKs//mf/+nZZ/ny5RZgvf/++z3bvvGNb1iXXXZZn2NPnTq1z3v57ne/a02dOtWyLGtA5+uhhx6yAGvbtm0Der8ikjynrYlIRNLizTffxLIsTj/99D7bQ6EQZWVlAHzyk5/E6XTy7LPP8rnPfY7Vq1dTUFDAwoULAdiyZQvvvfceBQUFfV4jGAzy/vvvD6qumTNn9rlfVVVFY2PjSb1ORUUFubm5nHLKKX22/fWvf+3znPPPPx/DMHruz507l1/+8pfEYrEBnS8At9t91HsQkdRRKBEZhuLxOKZpsmXLFkzT7PNYfn4+kPiC/cxnPsNjjz3G5z73OR577DGuvvpqnE5nz2vMnj37qJ94AEaPHj2oulwuV5/7hmH0XCzqcCR+TbZ6XecRiURO+DqGYRz3dQdiIOcLwOv19gk2IpJaCiUiw9CsWbOIxWI0NjYyb968Y+53zTXXsHDhQt5++21efPFFfvazn/U8ds455/Dkk09SXl5OYWFh2mvuDjr19fXMmjULoM9Frydr06ZNR92fPHkypmkO+HyJSHrpQleRYej000/nmmuu4dprr+Wpp56ipqaGN954g3//939nzZo1PftdfPHFVFRUcM011zBx4kTOP//8nseuueYaRo0axeLFi3nppZeoqalhw4YN3Hzzzezbty/lNXu9Xs4//3zuuOMO3nnnHTZu3MgPfvCDlL1+XV0dS5cuZdeuXTz++OP813/9FzfffDMw8PMlIumlUCIyTD300ENce+21fOc732HKlCn8wz/8A6+//jrjxo3r2ccwDD7/+c/zt7/9jWuuuabP83Nzc9m4cSPjx4/n05/+NFOnTuUrX/kKgUAgbSMnDz74IJFIhDlz5nDzzTfz85//PGWvfe211xIIBPjIRz7CDTfcwE033cQ///M/9zw+kPMlIullWL1/wBURERGxiUZKREREJCsolIjISautrSU/P/+Yt9raWrtLFJEhQD/fiMhJi0ajx11+fuLEiT1TjUVEjkWhRERERLKCfr4RERGRrKBQIiIiIllBoURERESygkKJiIiIZAWFEhEREckKCiUiIiKSFRRKREREJCsolIiIiEhW+P9t8cYR4KajJAAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "spontaneous_recombination_rate.loc[1,0,:].plot()" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [], + "source": [ + "from tardis.plasma.equilibrium.level_populations import LevelPopulationSolver\n", + "from tardis.plasma.equilibrium.rate_matrix import RateMatrix\n", + "\n", + "rate_matrix_solver = RateMatrix(rate_solvers, cmfgen_atom_data.levels)\n", + "\n", + "rate_matrix = rate_matrix_solver.solve(rad_field, electron_dist)\n", + "\n", + "lte_rate_matrix = RateMatrix(lte_rate_solvers, cmfgen_atom_data.levels).solve(rad_field, electron_dist)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "solver = LevelPopulationSolver(rate_matrix, cmfgen_atom_data.levels)\n", + "\n", + "level_pops = solver.solve()\n", + "\n", + "lte_level_pops = LevelPopulationSolver(lte_rate_matrix, cmfgen_atom_data.levels).solve()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.scatter(cmfgen_atom_data.levels.loc[1,0].energy * u.erg.to('eV'), level_pops.loc[1,0,:][0], marker='x', label='TARDIS')\n", + "plt.scatter(cmfgen_atom_data.levels.loc[1,0].energy * u.erg.to('eV'), lte_level_pops.loc[1,0,:][0], marker='x', label='TARDIS col only')\n", + "plt.xlabel(\"Energy (eV)\")\n", + "plt.ylabel(\"Population\")\n", + "plt.semilogy()\n", + "plt.legend()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "tardis", + "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.12.4" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/docs/resources/credits.rst b/docs/resources/credits.rst index 22ef542db45..b3429c6cc67 100644 --- a/docs/resources/credits.rst +++ b/docs/resources/credits.rst @@ -74,14 +74,14 @@ The following BibTeX entries are needed for the references: adsnote = {Provided by the SAO/NASA Astrophysics Data System} } -.. |CITATION| replace:: kerzendorf_2025_14633332 +.. |CITATION| replace:: kerzendorf_2025_14740815 -.. |DOI_BADGE| image:: https://img.shields.io/badge/DOI-10.5281/zenodo.14633332-blue - :target: https://doi.org/10.5281/zenodo.14633332 +.. |DOI_BADGE| image:: https://img.shields.io/badge/DOI-10.5281/zenodo.14740815-blue + :target: https://doi.org/10.5281/zenodo.14740815 .. code-block:: bibtex - @software{kerzendorf_2025_14633332, + @software{kerzendorf_2025_14740815, author = {Kerzendorf, Wolfgang and Sim, Stuart and Vogl, Christian and @@ -189,12 +189,12 @@ The following BibTeX entries are needed for the references: Nayak U, Ashwin and Kumar, Atul and Kharkar, Atharwa}, - title = {tardis-sn/tardis: TARDIS v2025.01.12}, + title = {tardis-sn/tardis: TARDIS v2025.01.26}, month = jan, year = 2025, publisher = {Zenodo}, - version = {release-2025.01.12}, - doi = {10.5281/zenodo.14633332}, - url = {https://doi.org/10.5281/zenodo.14633332}, + version = {release-2025.01.26}, + doi = {10.5281/zenodo.14740815}, + url = {https://doi.org/10.5281/zenodo.14740815}, } diff --git a/tardis/opacities/opacity_state.py b/tardis/opacities/opacity_state.py index 3efe183a21a..6db3b5791b9 100644 --- a/tardis/opacities/opacity_state.py +++ b/tardis/opacities/opacity_state.py @@ -7,105 +7,6 @@ from tardis.opacities.tau_sobolev import calculate_sobolev_line_opacity from tardis.transport.montecarlo.configuration import montecarlo_globals - -class OpacityState: - def __init__( - self, - electron_density, - t_electrons, - line_list_nu, - tau_sobolev, - beta_sobolev, - continuum_state, - ): - """ - Opacity State in Python - - Parameters - ---------- - electron_density : pd.DataFrame - t_electrons : numpy.ndarray - line_list_nu : pd.DataFrame - tau_sobolev : pd.DataFrame - beta_sobolev : pd.DataFrame - continuum_state: tardis.opacities.continuum.continuum_state.ContinuumState - """ - self.electron_density = electron_density - self.t_electrons = t_electrons - self.line_list_nu = line_list_nu - - self.tau_sobolev = tau_sobolev - - self.beta_sobolev = beta_sobolev - - # Continuum Opacity Data - self.continuum_state = continuum_state - - @classmethod - def from_legacy_plasma(cls, plasma, tau_sobolev): - """ - Generates an OpacityStatePython object from a tardis BasePlasma - - Parameters - ---------- - plasma : tardis.plasma.BasePlasma - legacy base plasma - tau_sobolev : pd.DataFrame - Expansion Optical Depths - - Returns - ------- - OpacityStatePython - """ - if hasattr(plasma, "photo_ion_cross_sections"): - continuum_state = ContinuumState.from_legacy_plasma(plasma) - else: - continuum_state = None - - return cls( - plasma.electron_densities, - plasma.t_electrons, - plasma.atomic_data.lines.nu, - tau_sobolev, - plasma.beta_sobolev, - continuum_state, - ) - - @classmethod - def from_plasma(cls, plasma, tau_sobolev, beta_sobolev): - """ - Generates an OpacityStatePython object from a tardis BasePlasma - - Parameters - ---------- - plasma : tardis.plasma.BasePlasma - legacy base plasma - tau_sobolev : pd.DataFrame - Expansion Optical Depths - beta_sobolev : pd.DataFrame - Modified expansion Optical Depths - - Returns - ------- - OpacityStatePython - """ - if hasattr(plasma, "photo_ion_cross_sections"): - continuum_state = ContinuumState.from_legacy_plasma(plasma) - else: - continuum_state = None - - atomic_data = plasma.atomic_data - - return cls( - plasma.electron_densities, - plasma.t_electrons, - atomic_data.lines.nu, - tau_sobolev, - beta_sobolev, - continuum_state, - ) - - opacity_state_spec = [ ("electron_density", float64[:]), ("t_electrons", float64[:]), @@ -131,7 +32,6 @@ def from_plasma(cls, plasma, tau_sobolev, beta_sobolev): ("k_packet_idx", int64), ] - @jitclass(opacity_state_spec) class OpacityStateNumba: def __init__( @@ -242,130 +142,227 @@ def __getitem__(self, i: slice): self.k_packet_idx, ) +class OpacityState: + def __init__( + self, + electron_density, + t_electrons, + line_list_nu, + tau_sobolev, + beta_sobolev, + continuum_state, + ): + """ + Opacity State in Python -def opacity_state_to_numba( - opacity_state: OpacityState, - macro_atom_state: MacroAtomState, - line_interaction_type, -) -> OpacityStateNumba: - """ - Initialize the OpacityStateNumba object and copy over the data over from OpacityState class + Parameters + ---------- + electron_density : pd.DataFrame + t_electrons : numpy.ndarray + line_list_nu : pd.DataFrame + tau_sobolev : pd.DataFrame + beta_sobolev : pd.DataFrame + continuum_state: tardis.opacities.continuum.continuum_state.ContinuumState + """ + self.electron_density = electron_density + self.t_electrons = t_electrons + self.line_list_nu = line_list_nu - Parameters - ---------- - opacity_state : tardis.opacities.opacity_state.OpacityState - line_interaction_type : enum - """ + self.tau_sobolev = tau_sobolev - electron_densities = opacity_state.electron_density.values - t_electrons = opacity_state.t_electrons - line_list_nu = opacity_state.line_list_nu.values + self.beta_sobolev = beta_sobolev - # NOTE: Disabled line scattering is handled by the opacitystate solver - tau_sobolev = np.ascontiguousarray( - opacity_state.tau_sobolev, dtype=np.float64 - ) + # Continuum Opacity Data + self.continuum_state = continuum_state - if line_interaction_type == "scatter": - # to adhere to data types, we must have an array of minimum size 1 - array_size = 1 - transition_probabilities = np.zeros( - (array_size, array_size), dtype=np.float64 - ) # to adhere to data types - line2macro_level_upper = np.zeros(array_size, dtype=np.int64) - macro_block_references = np.zeros(array_size, dtype=np.int64) - transition_type = np.zeros(array_size, dtype=np.int64) - destination_level_id = np.zeros(array_size, dtype=np.int64) - transition_line_id = np.zeros(array_size, dtype=np.int64) - else: - transition_probabilities = np.ascontiguousarray( - macro_atom_state.transition_probabilities.values.copy(), - dtype=np.float64, - ) - line2macro_level_upper = macro_atom_state.line2macro_level_upper - # TODO: Fix setting of block references for non-continuum mode + @classmethod + def from_legacy_plasma(cls, plasma, tau_sobolev): + """ + Generates an OpacityStatePython object from a tardis BasePlasma - macro_block_references = np.asarray( - macro_atom_state.macro_block_references + Parameters + ---------- + plasma : tardis.plasma.BasePlasma + legacy base plasma + tau_sobolev : pd.DataFrame + Expansion Optical Depths + + Returns + ------- + OpacityStatePython + """ + if hasattr(plasma, "photo_ion_cross_sections"): + continuum_state = ContinuumState.from_legacy_plasma(plasma) + else: + continuum_state = None + + return cls( + plasma.electron_densities, + plasma.t_electrons, + plasma.atomic_data.lines.nu, + tau_sobolev, + plasma.beta_sobolev, + continuum_state, ) - transition_type = macro_atom_state.transition_type.values + @classmethod + def from_plasma(cls, plasma, tau_sobolev, beta_sobolev): + """ + Generates an OpacityStatePython object from a tardis BasePlasma - # Destination level is not needed and/or generated for downbranch - destination_level_id = macro_atom_state.destination_level_id.values - transition_line_id = macro_atom_state.transition_line_id.values + Parameters + ---------- + plasma : tarids.plasma.BasePlasma + legacy base plasma + tau_sobolev : pd.DataFrame + Expansion Optical Depths + beta_sobolev : pd.DataFrame + Modified expansion Optical Depths - if montecarlo_globals.CONTINUUM_PROCESSES_ENABLED: - bf_threshold_list_nu = ( - opacity_state.continuum_state.bf_threshold_list_nu.values - ) - p_fb_deactivation = np.ascontiguousarray( - opacity_state.continuum_state.p_fb_deactivation.values.copy(), - dtype=np.float64, - ) + Returns + ------- + OpacityStatePython + """ + if hasattr(plasma, "photo_ion_cross_sections"): + continuum_state = ContinuumState.from_legacy_plasma(plasma) + else: + continuum_state = None - phot_nus = opacity_state.continuum_state.phot_nus - photo_ion_block_references = ( - opacity_state.continuum_state.photo_ion_block_references - ) - photo_ion_nu_threshold_mins = ( - opacity_state.continuum_state.photo_ion_nu_threshold_mins.values - ) - photo_ion_nu_threshold_maxs = ( - opacity_state.continuum_state.photo_ion_nu_threshold_maxs.values + atomic_data = plasma.atomic_data + + return cls( + plasma.electron_densities, + plasma.t_electrons, + atomic_data.lines.nu, + tau_sobolev, + beta_sobolev, + continuum_state, ) - chi_bf = opacity_state.continuum_state.chi_bf.values - x_sect = opacity_state.continuum_state.x_sect.values + def to_numba( + self, + macro_atom_state: MacroAtomState, + line_interaction_type, + ) -> OpacityStateNumba: + """ + Initialize the OpacityStateNumba object and copy over the data over from OpacityState class - phot_nus = phot_nus.values - ff_opacity_factor = ( - opacity_state.continuum_state.ff_cooling_factor - / np.sqrt(t_electrons) - ).astype(np.float64) - emissivities = opacity_state.continuum_state.emissivities.values - photo_ion_activation_idx = ( - opacity_state.continuum_state.photo_ion_activation_idx.values + Parameters + ---------- + macro_atom_state : tardis.opacities.macro_atom.macroatom_state.MacroAtomState + line_interaction_type : enum + """ + + electron_densities = self.electron_density.values + t_electrons = self.t_electrons + line_list_nu = self.line_list_nu.values + + # NOTE: Disabled line scattering is handled by the opacitystate solver + tau_sobolev = np.ascontiguousarray( + self.tau_sobolev, dtype=np.float64 ) - k_packet_idx = np.int64(opacity_state.continuum_state.k_packet_idx) - else: - bf_threshold_list_nu = np.zeros(0, dtype=np.float64) - p_fb_deactivation = np.zeros((0, 0), dtype=np.float64) - photo_ion_nu_threshold_mins = np.zeros(0, dtype=np.float64) - photo_ion_nu_threshold_maxs = np.zeros(0, dtype=np.float64) - photo_ion_block_references = np.zeros(0, dtype=np.int64) - chi_bf = np.zeros((0, 0), dtype=np.float64) - x_sect = np.zeros(0, dtype=np.float64) - phot_nus = np.zeros(0, dtype=np.float64) - ff_opacity_factor = np.zeros(0, dtype=np.float64) - emissivities = np.zeros((0, 0), dtype=np.float64) - photo_ion_activation_idx = np.zeros(0, dtype=np.int64) - k_packet_idx = np.int64(-1) - return OpacityStateNumba( - electron_densities, - t_electrons, - line_list_nu, - tau_sobolev, - transition_probabilities, - line2macro_level_upper, - macro_block_references, - transition_type, - destination_level_id, - transition_line_id, - bf_threshold_list_nu, - p_fb_deactivation, - photo_ion_nu_threshold_mins, - photo_ion_nu_threshold_maxs, - photo_ion_block_references, - chi_bf, - x_sect, - phot_nus, - ff_opacity_factor, - emissivities, - photo_ion_activation_idx, - k_packet_idx, - ) + if line_interaction_type == "scatter": + # to adhere to data types, we must have an array of minimum size 1 + array_size = 1 + transition_probabilities = np.zeros( + (array_size, array_size), dtype=np.float64 + ) # to adhere to data types + line2macro_level_upper = np.zeros(array_size, dtype=np.int64) + macro_block_references = np.zeros(array_size, dtype=np.int64) + transition_type = np.zeros(array_size, dtype=np.int64) + destination_level_id = np.zeros(array_size, dtype=np.int64) + transition_line_id = np.zeros(array_size, dtype=np.int64) + else: + transition_probabilities = np.ascontiguousarray( + macro_atom_state.transition_probabilities.values.copy(), + dtype=np.float64, + ) + line2macro_level_upper = macro_atom_state.line2macro_level_upper + + # TODO: Fix setting of block references for non-continuum mode + + macro_block_references = np.asarray( + macro_atom_state.macro_block_references + ) + + transition_type = macro_atom_state.transition_type.values + + # Destination level is not needed and/or generated for downbranch + destination_level_id = macro_atom_state.destination_level_id.values + transition_line_id = macro_atom_state.transition_line_id.values + + if montecarlo_globals.CONTINUUM_PROCESSES_ENABLED: + bf_threshold_list_nu = ( + self.continuum_state.bf_threshold_list_nu.values + ) + p_fb_deactivation = np.ascontiguousarray( + self.continuum_state.p_fb_deactivation.values.copy(), + dtype=np.float64, + ) + + phot_nus = self.continuum_state.phot_nus + photo_ion_block_references = ( + self.continuum_state.photo_ion_block_references + ) + photo_ion_nu_threshold_mins = ( + self.continuum_state.photo_ion_nu_threshold_mins.values + ) + photo_ion_nu_threshold_maxs = ( + self.continuum_state.photo_ion_nu_threshold_maxs.values + ) + + chi_bf = self.continuum_state.chi_bf.values + x_sect = self.continuum_state.x_sect.values + + phot_nus = phot_nus.values + ff_opacity_factor = ( + self.continuum_state.ff_cooling_factor + / np.sqrt(t_electrons) + ).astype(np.float64) + emissivities = self.continuum_state.emissivities.values + photo_ion_activation_idx = ( + self.continuum_state.photo_ion_activation_idx.values + ) + k_packet_idx = np.int64(self.continuum_state.k_packet_idx) + else: + bf_threshold_list_nu = np.zeros(0, dtype=np.float64) + p_fb_deactivation = np.zeros((0, 0), dtype=np.float64) + photo_ion_nu_threshold_mins = np.zeros(0, dtype=np.float64) + photo_ion_nu_threshold_maxs = np.zeros(0, dtype=np.float64) + photo_ion_block_references = np.zeros(0, dtype=np.int64) + chi_bf = np.zeros((0, 0), dtype=np.float64) + x_sect = np.zeros(0, dtype=np.float64) + phot_nus = np.zeros(0, dtype=np.float64) + ff_opacity_factor = np.zeros(0, dtype=np.float64) + emissivities = np.zeros((0, 0), dtype=np.float64) + photo_ion_activation_idx = np.zeros(0, dtype=np.int64) + k_packet_idx = np.int64(-1) + + return OpacityStateNumba( + electron_densities, + t_electrons, + line_list_nu, + tau_sobolev, + transition_probabilities, + line2macro_level_upper, + macro_block_references, + transition_type, + destination_level_id, + transition_line_id, + bf_threshold_list_nu, + p_fb_deactivation, + photo_ion_nu_threshold_mins, + photo_ion_nu_threshold_maxs, + photo_ion_block_references, + chi_bf, + x_sect, + phot_nus, + ff_opacity_factor, + emissivities, + photo_ion_activation_idx, + k_packet_idx, + ) def opacity_state_initialize( diff --git a/tardis/opacities/tests/test_opacity_state_numba.py b/tardis/opacities/tests/test_opacity_state_numba.py index 167fd2673b6..96111f6d6d6 100644 --- a/tardis/opacities/tests/test_opacity_state_numba.py +++ b/tardis/opacities/tests/test_opacity_state_numba.py @@ -1,5 +1,5 @@ import pytest -from tardis.opacities.opacity_state import opacity_state_to_numba +from tardis.opacities.opacity_state import OpacityState from tardis.opacities.opacity_solver import OpacitySolver from tardis.opacities.macro_atom.macroatom_solver import MacroAtomSolver import numpy.testing as npt @@ -35,8 +35,8 @@ def test_opacity_state_to_numba( ) else: macro_atom_state = None - actual = opacity_state_to_numba( - opacity_state, macro_atom_state, line_interaction_type + actual = opacity_state.to_numba( + macro_atom_state, line_interaction_type ) if sliced: diff --git a/tardis/plasma/electron_energy_distribution/base.py b/tardis/plasma/electron_energy_distribution/base.py index bbd9cc48542..892e3376a74 100644 --- a/tardis/plasma/electron_energy_distribution/base.py +++ b/tardis/plasma/electron_energy_distribution/base.py @@ -21,7 +21,7 @@ class ThermalElectronEnergyDistribution(ElectronEnergyDistribution): temperature : Quantity Electron temperatures in Kelvin. number_density : Quantity - Electron number densities in g/cm^3. + Electron number densities in cm^-3. """ temperature: u.Quantity diff --git a/tardis/plasma/equilibrium/rates/__init__.py b/tardis/plasma/equilibrium/rates/__init__.py index 2d7e4b79bc2..c34e7f5448d 100644 --- a/tardis/plasma/equilibrium/rates/__init__.py +++ b/tardis/plasma/equilibrium/rates/__init__.py @@ -2,9 +2,24 @@ UpsilonCMFGENSolver, UpsilonRegemorterSolver, ) +from tardis.plasma.equilibrium.rates.collisional_ionization_rates import ( + CollisionalIonizationRateSolver, +) +from tardis.plasma.equilibrium.rates.collisional_ionization_strengths import ( + CollisionalIonizationSeaton, +) from tardis.plasma.equilibrium.rates.collisional_rates import ( ThermalCollisionalRateSolver, ) +from tardis.plasma.equilibrium.rates.photoionization_rates import ( + AnalyticPhotoionizationRateSolver, + EstimatedPhotoionizationRateSolver, +) +from tardis.plasma.equilibrium.rates.photoionization_strengths import ( + AnalyticPhotoionizationCoeffSolver, + EstimatedPhotoionizationCoeffSolver, + SpontaneousRecombinationCoeffSolver, +) from tardis.plasma.equilibrium.rates.radiative_rates import ( RadiativeRatesSolver, ) diff --git a/tardis/plasma/equilibrium/rates/collisional_ionization_rates.py b/tardis/plasma/equilibrium/rates/collisional_ionization_rates.py new file mode 100644 index 00000000000..d35ee40dc79 --- /dev/null +++ b/tardis/plasma/equilibrium/rates/collisional_ionization_rates.py @@ -0,0 +1,50 @@ +from tardis.plasma.equilibrium.rates.collisional_ionization_strengths import ( + CollisionalIonizationSeaton, +) + + +class CollisionalIonizationRateSolver: + """Solver for collisional ionization and recombination rates.""" + + def __init__(self, photoionization_cross_sections): + self.photoionization_cross_sections = photoionization_cross_sections + + def solve(self, electron_temperature, saha_factor, approximation="seaton"): + """Solve the collisional ionization and recombination rates. + + Parameters + ---------- + electron_temperature : u.Quantity + Electron temperatures per cell + saha_factor : pandas.DataFrame, dtype float + The Saha factor for each cell. Indexed by atom number, ion number, level number. + approximation : str, optional + The rate approximation to use, by default "seaton" + + Returns + ------- + pd.DataFrame + Collisional ionization rates + pd.DataFrame + Collisional recombination rates + + Raises + ------ + ValueError + If an unsupported approximation is requested. + """ + if approximation == "seaton": + strength_solver = CollisionalIonizationSeaton( + self.photoionization_cross_sections + ) + else: + raise ValueError(f"approximation {approximation} not supported") + + collision_ionization_rates = strength_solver.solve(electron_temperature) + + # Inverse of the ionization rate for equilibrium + collision_recombination_rates = collision_ionization_rates.multiply( + saha_factor.loc[collision_ionization_rates.index] + ) + + return collision_ionization_rates, collision_recombination_rates diff --git a/tardis/plasma/equilibrium/rates/collisional_ionization_strengths.py b/tardis/plasma/equilibrium/rates/collisional_ionization_strengths.py new file mode 100644 index 00000000000..3b6db851057 --- /dev/null +++ b/tardis/plasma/equilibrium/rates/collisional_ionization_strengths.py @@ -0,0 +1,59 @@ +import astropy.units as u +import numpy as np +import pandas as pd + +from tardis import constants as const + +H = const.h.cgs +K_B = const.k_B.cgs + + +class CollisionalIonizationSeaton: + """Solver for collisional ionization rate coefficients in the Seaton approximation.""" + + def __init__(self, photoionization_cross_sections): + self.photoionization_cross_sections = photoionization_cross_sections + + def solve(self, electron_temperature): + """ + Parameters + ---------- + electron_temperature : u.Quantity + The electron temperature in K. + + Returns + ------- + pandas.DataFrame, dtype float + The rate coefficient for collisional ionization in the Seaton + approximation. Multiply with the electron density and the + level number density to obtain the total rate. + + Notes + ----- + The rate coefficient for collisional ionization in the Seaton approximation + is calculated according to Eq. 9.60 in [1]. + + References + ---------- + .. [1] Hubeny, I. and Mihalas, D., "Theory of Stellar Atmospheres". 2014. + """ + photo_ion_cross_sections_threshold = ( + self.photoionization_cross_sections.groupby(level=[0, 1, 2]).first() + ) + nu_i = photo_ion_cross_sections_threshold["nu"] + u0s = ( + nu_i.values[np.newaxis].T * u.Hz / electron_temperature * (H / K_B) + ) + factor = np.exp(-u0s) / u0s + factor = pd.DataFrame(factor, index=nu_i.index) + coll_ion_coeff = 1.55e13 * photo_ion_cross_sections_threshold["x_sect"] + coll_ion_coeff = factor.multiply(coll_ion_coeff, axis=0) + coll_ion_coeff = coll_ion_coeff.divide( + np.sqrt(electron_temperature), axis=1 + ) + + ion_number = coll_ion_coeff.index.get_level_values("ion_number").values + coll_ion_coeff[ion_number == 0] *= 0.1 + coll_ion_coeff[ion_number == 1] *= 0.2 + coll_ion_coeff[ion_number >= 2] *= 0.3 + return coll_ion_coeff diff --git a/tardis/plasma/equilibrium/rates/photoionization_rates.py b/tardis/plasma/equilibrium/rates/photoionization_rates.py new file mode 100644 index 00000000000..720ecea1945 --- /dev/null +++ b/tardis/plasma/equilibrium/rates/photoionization_rates.py @@ -0,0 +1,214 @@ +from tardis.plasma.equilibrium.rates.photoionization_strengths import ( + AnalyticPhotoionizationCoeffSolver, + EstimatedPhotoionizationCoeffSolver, + SpontaneousRecombinationCoeffSolver, +) + + +class AnalyticPhotoionizationRateSolver: + """Solve the photoionization and spontaneous recombination rates in the + case where the radiation field is computed analytically. + """ + + def __init__(self, photoionization_cross_sections): + self.photoionization_cross_sections = photoionization_cross_sections + + self.spontaneous_recombination_rate_coeff_solver = ( + SpontaneousRecombinationCoeffSolver( + self.photoionization_cross_sections + ) + ) + + def compute_rates( + self, + photoionization_rate_coeff, + stimulated_recombination_rate_coeff, + spontaneous_recombination_rate_coeff, + level_number_density, + ion_number_density, + electron_number_density, + saha_factor, + ): + """Compute the photoionization and spontaneous recombination rates + + Parameters + ---------- + photoionization_rate_coeff : pd.DataFrame + The photoionization rate coefficients for each transition. + Columns are cells. + stimulated_recombination_rate_coeff : pd.DataFrame + The stimulated recombination rate coefficients for each transition. + Columns are cells. + spontaneous_recombination_rate_coeff : pd.DataFrame + The spontaneous recombination rate coefficients for each transition. + Columns are cells. + level_number_density : pd.DataFrame + The electron energy level number density. Columns are cells. + ion_number_density : pd.DataFrame + The ion number density. Columns are cells. + electron_number_density : u.Quantity + The free electron number density per cell. + saha_factor : pd.DataFrame + The LTE population factor. Columns are cells. + + Returns + ------- + pd.DataFrame + Photoionization rate for each electron energy level. Columns are cells + pd.DataFrame + Spontaneous recombination rate for each electron energy level. Columns are cells + """ + photoionization_rate = ( + photoionization_rate_coeff * level_number_density + - saha_factor + * stimulated_recombination_rate_coeff + * ion_number_density + * electron_number_density + ) + spontaneous_recombination_rate = ( + saha_factor + * spontaneous_recombination_rate_coeff + * ion_number_density + * electron_number_density + ) + + return photoionization_rate, spontaneous_recombination_rate + + def solve( + self, + dilute_blackbody_radiationfield_state, + electron_energy_distribution, + level_number_density, + ion_number_density, + saha_factor, + ): + """Solve the photoionization and spontaneous recombination rates in the + case where the radiation field is not estimated. + + Parameters + ---------- + dilute_blackbody_radiationfield_state : DiluteBlackBodyRadiationFieldState + A dilute black body radiation field state. + electron_energy_distribution : ThermalElectronEnergyDistribution + Electron properties. + level_number_density : pd.DataFrame + Electron energy level number density. Columns are cells. + ion_number_density : pd.DataFrame + Ion number density. Columns are cells. + saha_factor : pd.DataFrame + Saha factor: the LTE level number density divided by the LTE ion + number density and the electron number density. + + Returns + ------- + pd.DataFrame + Photoionization rate. Columns are cells. + pd.DataFrame + Spontaneous recombination rate. Columns are cells. + """ + photoionization_rate_coeff_solver = AnalyticPhotoionizationCoeffSolver( + self.photoionization_cross_sections + ) + + photoionization_rate_coeff, stimulated_recombination_rate_coeff = ( + photoionization_rate_coeff_solver.solve( + dilute_blackbody_radiationfield_state, + electron_energy_distribution.temperature, + ) + ) + + spontaneous_recombination_rate_coeff = ( + self.spontaneous_recombination_rate_coeff_solver.solve( + electron_energy_distribution.temperature + ) + ) + + return self.compute_rates( + photoionization_rate_coeff, + stimulated_recombination_rate_coeff, + spontaneous_recombination_rate_coeff, + level_number_density, + ion_number_density, + electron_energy_distribution.number_density, + saha_factor, + ) + + +class EstimatedPhotoionizationRateSolver(AnalyticPhotoionizationRateSolver): + """Solve the photoionization and spontaneous recombination rates in the + case where the radiation field is estimated by Monte Carlo processes. + """ + + def __init__( + self, photoionization_cross_sections, level2continuum_edge_idx + ): + super().__init__( + photoionization_cross_sections, + ) + self.level2continuum_edge_idx = level2continuum_edge_idx + + def solve( + self, + electron_energy_distribution, + radfield_mc_estimators, + time_simulation, + volume, + level_number_density, + ion_number_density, + saha_factor, + ): + """Solve the photoionization and spontaneous recombination rates in the + case where the radiation field is estimated by Monte Carlo processes. + + Parameters + ---------- + electron_energy_distribution : ThermalElectronEnergyDistribution + Electron properties. + radfield_mc_estimators : RadiationFieldMCEstimators + Estimators of the radiation field properties. + time_simulation : u.Quantity + Time of simulation. + volume : u.Quantity + Volume per cell. + level_number_density : pd.DataFrame + Electron energy level number density. Columns are cells. + ion_number_density : pd.DataFrame + Ion number density. Columns are cells. + saha_factor : pd.DataFrame + Saha factor: the LTE level number density divided by the LTE ion + number density and the electron number density. + + Returns + ------- + pd.DataFrame + Photoionization rate. Columns are cells. + pd.DataFrame + Spontaneous recombination rate. Columns are cells. + """ + photoionization_rate_coeff_solver = EstimatedPhotoionizationCoeffSolver( + self.level2continuum_edge_idx + ) + + photoionization_rate_coeff, stimulated_recombination_rate_coeff = ( + photoionization_rate_coeff_solver.solve( + radfield_mc_estimators, + time_simulation, + volume, + ) + ) + + spontaneous_recombination_rate_coeff = ( + self.spontaneous_recombination_rate_coeff_solver.solve( + electron_energy_distribution.temperature + ) + ) + + return self.compute_rates( + photoionization_rate_coeff, + stimulated_recombination_rate_coeff, + spontaneous_recombination_rate_coeff, + level_number_density, + ion_number_density, + electron_energy_distribution.number_density, + saha_factor, + ) diff --git a/tardis/plasma/equilibrium/rates/photoionization_strengths.py b/tardis/plasma/equilibrium/rates/photoionization_strengths.py new file mode 100644 index 00000000000..f974dbdbea6 --- /dev/null +++ b/tardis/plasma/equilibrium/rates/photoionization_strengths.py @@ -0,0 +1,328 @@ +import numpy as np +import pandas as pd + +from tardis import constants as const +from tardis.transport.montecarlo.estimators.util import ( + bound_free_estimator_array2frame, + integrate_array_by_blocks, +) + +C = const.c.cgs.value +H = const.h.cgs.value +K_B = const.k_B.cgs.value + + +class SpontaneousRecombinationCoeffSolver: + def __init__( + self, + photoionization_cross_sections, + ): + self.photoionization_cross_sections = photoionization_cross_sections + self.nu = self.photoionization_cross_sections.nu.values + + self.photoionization_block_references = np.pad( + self.photoionization_cross_sections.nu.groupby(level=[0, 1, 2]) + .count() + .values.cumsum(), + [1, 0], + ) + + self.photoionization_index = ( + self.photoionization_cross_sections.index.unique() + ) + + @property + def common_prefactor(self): + """Used to multiply with both spontaneous recombination and + photoionization coefficients. + + Returns + ------- + pd.DataFrame + A dataframe of the prefactor. + """ + return ( + 4.0 + * np.pi + * self.photoionization_cross_sections.x_sect + / (H * self.nu) + ) + + def calculate_photoionization_boltzmann_factor(self, electron_temperature): + """Calculate the Boltzmann factor at each photoionization frequency + + Parameters + ---------- + electron_temperature : Quantity + Electron temperature in each shell. + + Returns + ------- + numpy.ndarray + The Boltzmann factor per shell per photoionization frequency. + """ + return np.exp(-self.nu[np.newaxis].T / electron_temperature * (H / K_B)) + + def solve(self, electron_temperature): + """ + Calculate the spontaneous recombination rate coefficient. + + Parameters + ---------- + electron_temperature : u.Quantity + Electron temperature in each cell. + + Returns + ------- + pd.DataFrame + The calculated spontaneous recombination rate coefficient. + + Notes + ----- + Equation 13 in Lucy 2003. + """ + prefactor = self.common_prefactor * (2 * H * self.nu**3.0) / (C**2.0) + photoionization_boltzmann_factor = pd.DataFrame( + self.calculate_photoionization_boltzmann_factor( + electron_temperature + ), + index=prefactor.index, + ) + spontaneous_recombination_rate_coeff = ( + photoionization_boltzmann_factor.multiply( + prefactor, + axis=0, + ) + ) + spontaneous_recombination_rate_coeff_integrated = ( + integrate_array_by_blocks( + spontaneous_recombination_rate_coeff.to_numpy(), + self.nu, + self.photoionization_block_references, + ) + ) + + return pd.DataFrame( + spontaneous_recombination_rate_coeff_integrated, + index=self.photoionization_index, + ) + + +class AnalyticPhotoionizationCoeffSolver(SpontaneousRecombinationCoeffSolver): + def __init__( + self, + photoionization_cross_sections, + ): + super().__init__(photoionization_cross_sections) + + def calculate_mean_intensity_photoionization_df( + self, + dilute_blackbody_radiationfield_state, + ): + """Calculates the mean intensity of the radiation field at each photoionization frequency. + + Parameters + ---------- + dilute_blackbody_radiationfield_state : DilutePlanckianRadiationField + The radiation field. + + Returns + ------- + pd.DataFrame + DataFrame of mean intensities indexed by photoionization levels and + columns of cells. + """ + mean_intensity = ( + dilute_blackbody_radiationfield_state.calculate_mean_intensity( + self.nu + ) + ) + return pd.DataFrame( + mean_intensity, + index=self.photoionization_cross_sections.index, + columns=np.arange( + len(dilute_blackbody_radiationfield_state.temperature) + ), + ) + + def calculate_photoionization_rate_coeff( + self, + mean_intensity_photoionization_df, + ): + """ + Calculate the photoionization rate coefficient. + + Parameters + ---------- + dilute_blackbody_radiationfield_state : DiluteBlackBodyRadiationFieldState + A dilute black body radiation field state. + + Returns + ------- + pd.DataFrame + The calculated photoionization rate coefficient. + + Notes + ----- + Equation 16 in Lucy 2003. + """ + photoionization_rate_coeff = mean_intensity_photoionization_df.multiply( + self.common_prefactor, + axis=0, + ) + photoionization_rate_coeff = integrate_array_by_blocks( + photoionization_rate_coeff.values, + self.nu, + self.photoionization_block_references, + ) + photoionization_rate_coeff = pd.DataFrame( + photoionization_rate_coeff, + index=self.photoionization_index, + ) + return photoionization_rate_coeff + + def calculate_stimulated_recombination_rate_coeff( + self, + mean_intensity_photoionization_df, + photoionization_boltzmann_factor, + ): + """ + Calculate the photoionization rate coefficient. + + Parameters + ---------- + mean_intensity_photoionization_df : pd.DataFrame + Mean intensity at each photoionization frequency. + photoionization_boltzmann_factor : np.ndarray + Boltzmann factor for each photoionization frequency. + + Returns + ------- + pd.DataFrame + The stimulated recombination rate coefficient. + + Notes + ----- + Equation 15 in Lucy 2003. + """ + stimulated_recombination_rate_coeff = ( + mean_intensity_photoionization_df * photoionization_boltzmann_factor + ) + + stimulated_recombination_rate_coeff = ( + mean_intensity_photoionization_df.multiply( + self.common_prefactor, + axis=0, + ) + ) + stimulated_recombination_rate_coeff = integrate_array_by_blocks( + stimulated_recombination_rate_coeff.values, + self.nu, + self.photoionization_block_references, + ) + stimulated_recombination_rate_coeff = pd.DataFrame( + stimulated_recombination_rate_coeff, + index=self.photoionization_index, + ) + return stimulated_recombination_rate_coeff + + def solve( + self, + dilute_blackbody_radiationfield_state, + electron_temperature, + ): + """ + Prepares the ionization and recombination coefficients by grouping them for + ion numbers. + + Parameters + ---------- + dilute_blackbody_radiationfield_state : DiluteBlackBodyRadiationFieldState + The dilute black body radiation field state. + electron_temperature : u.Quantity + Electron temperature in each shell. + + Returns + ------- + photoionization_rate_coeff + Photoionization rate coefficient grouped by atomic number and ion number. + recombination_rate_coeff + Radiative recombination rate coefficient grouped by atomic number and ion number. + """ + photoionization_boltzmann_factor = ( + self.calculate_photoionization_boltzmann_factor( + electron_temperature + ) + ) + + mean_intensity_photoionization_df = ( + self.calculate_mean_intensity_photoionization_df( + dilute_blackbody_radiationfield_state + ) + ) + # Equation 16 Lucy 2003 + photoionization_rate_coeff = self.calculate_photoionization_rate_coeff( + mean_intensity_photoionization_df, + ) + # Equation 15 Lucy 2003. Must be multiplied by Saha LTE factor Phi_ik + stimulated_recombination_rate_coeff = ( + self.calculate_stimulated_recombination_rate_coeff( + mean_intensity_photoionization_df, + photoionization_boltzmann_factor, + ) + ) + + return ( + photoionization_rate_coeff, + stimulated_recombination_rate_coeff, + ) + + +class EstimatedPhotoionizationCoeffSolver: + def __init__( + self, + level2continuum_edge_idx, + ): + self.level2continuum_edge_idx = level2continuum_edge_idx + + def solve( + self, + radfield_mc_estimators, + time_simulation, + volume, + ): + """ + Solve for the continuum properties. + + Parameters + ---------- + radfield_mc_estimators : RadiationFieldMCEstimators + The Monte Carlo estimators for the radiation field. + time_simulation : float + The simulation time. + volume : float + The volume of the cells. + + Returns + ------- + ContinuumProperties + The calculated continuum properties. + """ + # TODO: the estimators are computed in the form epsilon_nu * distance * xsection / comoving_nu + # with the stimulated recombination multiplied by a Boltzmann factor exp(-h * comoving_nu / k * electron_temp) + # This is why this method does not match the one in AnalyticPhotoionizationCoeffSolver + photoionization_normalization = (time_simulation * volume * H) ** -1 + + photoionization_rate_coeff = bound_free_estimator_array2frame( + radfield_mc_estimators.photo_ion_estimator, + self.level2continuum_edge_idx, + ) + photoionization_rate_coeff *= photoionization_normalization + + stimulated_recombination_rate_coeff = bound_free_estimator_array2frame( + radfield_mc_estimators.stim_recomb_estimator, + self.level2continuum_edge_idx, + ) + stimulated_recombination_rate_coeff *= photoionization_normalization + + return photoionization_rate_coeff, stimulated_recombination_rate_coeff diff --git a/tardis/spectrum/formal_integral.py b/tardis/spectrum/formal_integral.py index 12c0ca56313..d754d492d43 100644 --- a/tardis/spectrum/formal_integral.py +++ b/tardis/spectrum/formal_integral.py @@ -9,8 +9,9 @@ from scipy.interpolate import interp1d from tardis import constants as const +from tardis.opacities.continuum.continuum_state import ContinuumState from tardis.opacities.opacity_state import ( - opacity_state_to_numba, + OpacityState, opacity_state_initialize, ) from tardis.spectrum.formal_integral_cuda import ( @@ -287,8 +288,7 @@ def __init__( self.transport.montecarlo_configuration ) if plasma and opacity_state and macro_atom_state: - self.opacity_state = opacity_state_to_numba( - opacity_state, + self.opacity_state = opacity_state.to_numba( macro_atom_state, transport.line_interaction_type, ) diff --git a/tardis/transport/montecarlo/base.py b/tardis/transport/montecarlo/base.py index 5fac3521335..9e708c96f06 100644 --- a/tardis/transport/montecarlo/base.py +++ b/tardis/transport/montecarlo/base.py @@ -7,8 +7,9 @@ from tardis import constants as const from tardis.io.logger import montecarlo_tracking as mc_tracker from tardis.io.util import HDFWriterMixin +from tardis.opacities.continuum.continuum_state import ContinuumState from tardis.opacities.opacity_state import ( - opacity_state_to_numba, + OpacityState, ) from tardis.transport.montecarlo.configuration.base import ( MonteCarloConfiguration, @@ -114,9 +115,9 @@ def initialize_transport_state( ) geometry_state = simulation_state.geometry.to_numba() - - opacity_state_numba = opacity_state_to_numba( - opacity_state, macro_atom_state, self.line_interaction_type + opacity_state_numba = opacity_state.to_numba( + macro_atom_state, + self.line_interaction_type, ) opacity_state_numba = opacity_state_numba[ simulation_state.geometry.v_inner_boundary_index : simulation_state.geometry.v_outer_boundary_index