{ "cells": [ { "cell_type": "raw", "metadata": {}, "source": [ "%This notebook demonstrates the use of the workpackage template, replace with your own.\n", "\n", "\\documentclass[english]{workpackage}[1996/06/02]\n", "\n", "% input the common preamble content (required by the ipnb2latex converter)\n", "\\input{header.tex}\n", "\n", "\n", "% then follows the rest of the preamble to be placed before the begin document\n", "% this preamble content is special to the documentclass you defined above.\n", "\\WPproject{Computational Radiometry} % project name\n", "\\WPequipment{} % equipment name\n", "\\WPsubject{06 Diverse pyradi utilities} % main heading \n", "\\WPconclusions{} \n", "\\WPclassification{} \n", "\\WPdocauthor{CJ Willers}\n", "\\WPcurrentpackdate{\\today}\n", "\\WPcurrentpacknumber{} % work package number\n", "\\WPdocnumber{} % this doc number hosts all the work packages\n", "\\WPprevpackdate{} % work package which this one supersedes\n", "\\WPprevpacknumber{} % work package which this one supersedes\n", "\\WPsuperpackdate{} % work package which comes after this one\n", "\\WPsuperpacknumber{} % work package which comes after this one\n", "\\WPdocontractdetails{false}\n", "\\WPcontractname{} % contract name \n", "\\WPorderno{} % contract order number\n", "\\WPmilestonenumber{} % contract milestone number\n", "\\WPmilestonetitle{} % contract milestone title\n", "\\WPcontractline{} % contract milestone line number \n", "\\WPdocECPnumber{} % ecp\\ecr number\n", "\\WPdistribution{}\n", "\n", " \n", "\n", "% this is entered just before the end{document}\n", "\\newcommand{\\atendofdoc}{\n", "\\bibliographystyle{IEEEtran}\n", "\\bibliography{references}\n", "}\n", "\n", "%and finally the document begin.\n", "\\begin{document}\n", "\\WPlayout\n" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "# 6 Diverse pyradi utilities" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "This notebook forms part of a series on [computational optical radiometry](https://github.com/NelisW/ComputationalRadiometry#computational-optical-radiometry-with-pyradi). The notebooks can be downloaded from [Github](https://github.com/NelisW/ComputationalRadiometry#computational-optical-radiometry-with-pyradi). These notebooks are constantly revised and updated, please revisit from time to time. \n", "\n", "\n", "[](http://dx.doi.org/10.5281/zenodo.9910)\n" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "The date of this document and module versions used in this document are given at the end of the file. \n", "Feedback is appreciated: neliswillers at gmail dot com." ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "## Overview" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "The pyradi library has a reasonably complete collection of Planck radiator models, both spectral and wide band. A comprehensive collection of physical constants, pertinent to optical radiation is also included. This notebook introduces these functions in the [`pyradi.ryplanck`](http://nelisw.github.io/pyradi-docs/_build/html/ryplanck.html#pyradi.ryplanck) library." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "run_control": { "frozen": false, "read_only": false } }, "outputs": [], "source": [ "from IPython.display import display\n", "from IPython.display import Image\n", "from IPython.display import HTML" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "## File and directory services" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "### List files in a directory" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "The [`ryfiles.listFiles`](http://nelisw.github.io/pyradi-docs/_build/html/ryfiles.html#pyradi.ryfiles.listFiles) function \n", "returns a list of file paths to files in a file system, searching a directory structure along the specified path, looking for files that matches the glob pattern. If specified, the search will continue into sub-directories. A list of matching names is returned.\n", "The function supports a local or network reachable filesystem, but not URLs.\n", "\n", "The function signature is: \n", "\n", "`listFiles(root, patterns='*', recurse=1, return_folders=0, useRegex=False)`\n", "\n", "- `root (string)` directory root from where the search must take place.\n", "- `patterns (string)` glob/regex pattern for filename matching.\n", "- `recurse (unt)` flag to indicate if subdirectories must also be searched (optional).\n", "- `return_folders (int)` flag to indicate if folder names must also be returned (optional).\n", "- `useRegex (bool)` flag to indicate if patterns areregular expression strings (optional)." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "run_control": { "frozen": false, "read_only": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['05a-PlottingWithPyradi-GeneralAndCartesian.ipynb', '05b-PlottingWithPyradi-Polar-and-3D.ipynb', '.ipynb_checkpoints\\\\05a-PlottingWithPyradi-GeneralAndCartesian-checkpoint.ipynb', '.ipynb_checkpoints\\\\05b-PlottingWithPyradi-Polar-and-3D-checkpoint.ipynb']\n" ] } ], "source": [ "import pyradi.ryfiles as ryfiles\n", "print(ryfiles.listFiles('.', '05*.ipynb'))" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "run_control": { "frozen": false, "read_only": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['05a-PlottingWithPyradi-GeneralAndCartesian.bbl', '05a-PlottingWithPyradi-GeneralAndCartesian.bib', '05a-PlottingWithPyradi-GeneralAndCartesian.blg', '05a-PlottingWithPyradi-GeneralAndCartesian.ipynb', '05a-PlottingWithPyradi-GeneralAndCartesian.log', '05a-PlottingWithPyradi-GeneralAndCartesian.pdf', '05a-PlottingWithPyradi-GeneralAndCartesian.tex', '05b-PlottingWithPyradi-Polar-and-3D.bbl', '05b-PlottingWithPyradi-Polar-and-3D.bib', '05b-PlottingWithPyradi-Polar-and-3D.blg', '05b-PlottingWithPyradi-Polar-and-3D.ipynb', '05b-PlottingWithPyradi-Polar-and-3D.log', '05b-PlottingWithPyradi-Polar-and-3D.pdf', '05b-PlottingWithPyradi-Polar-and-3D.tex', '05', '.git\\\\objects\\\\05', '.git\\\\objects\\\\06\\\\05114a0c65030f3ab63c80822e07a03700ef35', '.git\\\\objects\\\\51\\\\059c45e9c9d5cfd46021370481d725b6eecccd', '.git\\\\objects\\\\7d\\\\05a03633f8985bca427403a7ff84ba35e199e6', '.git\\\\objects\\\\99\\\\05b062b719b832d586eb36cbd3ec21fe74d328', '.git\\\\objects\\\\d8\\\\052cd1780d21a06306ec3090e7657f3ea3e17e', '.git\\\\objects\\\\d9\\\\05835262f64f8ca948425e72c5f4ee07305246', '.ipynb_checkpoints\\\\05a-PlottingWithPyradi-GeneralAndCartesian-checkpoint.ipynb', '.ipynb_checkpoints\\\\05a-PlottingWithPyradi-GeneralAndCartesian-checkpoint.pdf', '.ipynb_checkpoints\\\\05b-PlottingWithPyradi-Polar-and-3D-checkpoint.ipynb', '05\\\\05-testfile.txt', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_13_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_14_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_15_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_17_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_19_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_22_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_22_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_22_2.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_25_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_25_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_25_2.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_29_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_31_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_32_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_34_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_34_2.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_34_3.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_34_4.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_34_5.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_35_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_36_2.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_38_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_38_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_38_2.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_38_3.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_38_4.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_38_5.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_40_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_40_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_42_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_43_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_44_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_46_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_47_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_48_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_50_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_51_2.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_52_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_53_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_55_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_56_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_57_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_60_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_62_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_64_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_66_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_68_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_70_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_72_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_73_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_74_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_75_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_76_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_77_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_78_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_79_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_81_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_83_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_84_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_85_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_87_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_90_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_91_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_91_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_93_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_96_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_97_1.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_10_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_11_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_12_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_13_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_14_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_15_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_16_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_17_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_17_1.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_18_1.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_22_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_23_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_24_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_25_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_26_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_27_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_28_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_29_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_32_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_33_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_35_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_36_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_38_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_39_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_39_1.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_40_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_41_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_46_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_47_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_49_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_50_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_54_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_55_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_58_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_59_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_59_1.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_61_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_62_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_64_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_65_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_72_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_73_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_74_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_75_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_75_1.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_83_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_84_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_85_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_86_0.png']\n" ] } ], "source": [ "import pyradi.ryfiles as ryfiles\n", "print(ryfiles.listFiles('.', '05*', return_folders=1))" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "run_control": { "frozen": false, "read_only": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[]\n" ] } ], "source": [ "print(ryfiles.listFiles('.', '05*.ipynb', useRegex=1))" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "run_control": { "frozen": false, "read_only": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['05a-PlottingWithPyradi-GeneralAndCartesian.bbl', '05a-PlottingWithPyradi-GeneralAndCartesian.bib', '05a-PlottingWithPyradi-GeneralAndCartesian.blg', '05a-PlottingWithPyradi-GeneralAndCartesian.ipynb', '05a-PlottingWithPyradi-GeneralAndCartesian.log', '05a-PlottingWithPyradi-GeneralAndCartesian.pdf', '05a-PlottingWithPyradi-GeneralAndCartesian.tex', '05b-PlottingWithPyradi-Polar-and-3D.bbl', '05b-PlottingWithPyradi-Polar-and-3D.bib', '05b-PlottingWithPyradi-Polar-and-3D.blg', '05b-PlottingWithPyradi-Polar-and-3D.ipynb', '05b-PlottingWithPyradi-Polar-and-3D.log', '05b-PlottingWithPyradi-Polar-and-3D.pdf', '05b-PlottingWithPyradi-Polar-and-3D.tex', 'tape7-05b', '.git\\\\objects\\\\08\\\\2acc4bfa2505b8d6fe0f85446c972a05a42888', '.git\\\\objects\\\\0d\\\\7f4c205a5f765006ae18faf50ffcbe4a430767', '.git\\\\objects\\\\2f\\\\f2c19f11277ca6a4895b1fbb5fcaf805abafbc', '.git\\\\objects\\\\38\\\\4b27fe04aa105acd84ebdbb16611729d4342f4', '.git\\\\objects\\\\40\\\\dcbc05a06258a6fc781a76ef07a6ee877b65b4', '.git\\\\objects\\\\6b\\\\42dbd56bd6e9ba05a54d592cfba05733d3aa20', '.git\\\\objects\\\\78\\\\656e36efc86279d39fe23f5b7b73dac4f005b9', '.git\\\\objects\\\\7b\\\\dc9c6407a05b826590e66950ca58484a07c69b', '.git\\\\objects\\\\7d\\\\05a03633f8985bca427403a7ff84ba35e199e6', '.git\\\\objects\\\\99\\\\05b062b719b832d586eb36cbd3ec21fe74d328', '.git\\\\objects\\\\9a\\\\e05ae9a9ff4ed6b6d749d024585c86433296da', '.git\\\\objects\\\\b2\\\\bd7a8d1a1ba253d4b5cdc12a05a1b76c8a573c', '.git\\\\objects\\\\c0\\\\91f1f7728fb9da43f08705ad8ebe3968cefea6', '.git\\\\objects\\\\e8\\\\8705a4a9b8a5682b325f4fe715b42065c08d43', '.ipynb_checkpoints\\\\05a-PlottingWithPyradi-GeneralAndCartesian-checkpoint.ipynb', '.ipynb_checkpoints\\\\05a-PlottingWithPyradi-GeneralAndCartesian-checkpoint.pdf', '.ipynb_checkpoints\\\\05b-PlottingWithPyradi-Polar-and-3D-checkpoint.ipynb', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_13_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_14_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_15_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_17_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_19_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_22_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_22_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_22_2.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_25_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_25_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_25_2.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_29_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_31_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_32_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_34_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_34_2.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_34_3.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_34_4.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_34_5.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_35_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_36_2.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_38_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_38_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_38_2.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_38_3.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_38_4.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_38_5.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_40_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_40_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_42_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_43_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_44_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_46_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_47_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_48_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_50_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_51_2.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_52_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_53_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_55_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_56_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_57_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_60_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_62_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_64_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_66_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_68_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_70_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_72_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_73_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_74_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_75_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_76_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_77_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_78_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_79_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_81_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_83_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_84_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_85_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_87_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_90_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_91_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_91_1.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_93_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_96_0.png', 'pic\\\\05a-PlottingWithPyradi-GeneralAndCartesian_97_1.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_10_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_11_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_12_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_13_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_14_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_15_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_16_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_17_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_17_1.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_18_1.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_22_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_23_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_24_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_25_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_26_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_27_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_28_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_29_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_32_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_33_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_35_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_36_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_38_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_39_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_39_1.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_40_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_41_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_46_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_47_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_49_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_50_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_54_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_55_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_58_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_59_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_59_1.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_61_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_62_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_64_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_65_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_72_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_73_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_74_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_75_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_75_1.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_83_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_84_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_85_0.png', 'pic\\\\05b-PlottingWithPyradi-Polar-and-3D_86_0.png']\n" ] } ], "source": [ "print(ryfiles.listFiles('.', '05[ab].*', useRegex=1))" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "run_control": { "frozen": false, "read_only": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['05a-PlottingWithPyradi-GeneralAndCartesian.bbl', '05a-PlottingWithPyradi-GeneralAndCartesian.bib', '05a-PlottingWithPyradi-GeneralAndCartesian.blg', '05a-PlottingWithPyradi-GeneralAndCartesian.ipynb', '05a-PlottingWithPyradi-GeneralAndCartesian.log', '05a-PlottingWithPyradi-GeneralAndCartesian.pdf', '05a-PlottingWithPyradi-GeneralAndCartesian.tex', '05b-PlottingWithPyradi-Polar-and-3D.bbl', '05b-PlottingWithPyradi-Polar-and-3D.bib', '05b-PlottingWithPyradi-Polar-and-3D.blg', '05b-PlottingWithPyradi-Polar-and-3D.ipynb', '05b-PlottingWithPyradi-Polar-and-3D.log', '05b-PlottingWithPyradi-Polar-and-3D.pdf', '05b-PlottingWithPyradi-Polar-and-3D.tex']\n" ] } ], "source": [ "print(ryfiles.listFiles('.', '^05[ab].*\\.[^h].*', useRegex=1))" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "### Working with filenames and directories" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "The information in this section does not describe any pyradi functionality, but may be useful in this context.\n", "\n", "To split a qualified path into the folder and filename use [`os.path.split`](https://docs.python.org/2/library/os.path.html), yielding a tuple with the path as the first element and the filename as the second element." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "run_control": { "frozen": false, "read_only": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "('05', '05-testfile.txt')\n", "05-testfile.txt\n", "05\n" ] } ], "source": [ "import os.path\n", "print(type(os.path.split('05\\\\05-testfile.txt')))\n", "print(os.path.split('05\\\\05-testfile.txt'))\n", "print(os.path.basename('05\\\\05-testfile.txt'))\n", "print(os.path.dirname('05\\\\05-testfile.txt'))" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "To join a folder and a filename use [`os.path.join`](https://docs.python.org/2/library/os.path.html), but there is a catch here; `os.path.join` does not expect the tuple, but simply a number of arguments. Study the documentation carefully, because [`os.path.join`](https://docs.python.org/2/library/os.path.html) may discard some arguments.\n", "\n", "Alternatively, use `os.sep` or `os.path.sep` to get the pathname separator ('/' for POSIX and '\\\\' for Windows) and then use this character to join the list." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "run_control": { "frozen": false, "read_only": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "05\\05-testfile.txt\n", "05\\05-testfile.txt\n", "05\\05-testfile.txt\n" ] } ], "source": [ "print(os.path.join(* os.path.split('05\\\\05-testfile.txt')))\n", "print(os.sep.join(os.path.split('05\\\\05-testfile.txt')))\n", "print(os.path.sep.join(os.path.split('05\\\\05-testfile.txt')))" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "To test for the existance of a file use `os.path.lexists` or `os.path.exists`, see the documentation to confirm the difference between the two." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "run_control": { "frozen": false, "read_only": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "True\n", "False\n" ] } ], "source": [ "print(os.path.lexists('05\\\\05-testfile.txt'))\n", "print(os.path.lexists('06\\\\05-testfile.txt'))" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "### Downloading compressed and uncompressed files from the internet" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "The [`ryfiles.downloadFileUrl`](http://nelisw.github.io/pyradi-docs/_build/html/ryfiles.html#pyradi.ryfiles.downloadFileUrl) function downloads a file from a URL. The URL is used to download a file, to the saveFilename specified. If no saveFilename is given, the basename of the URL is used. The name of the downloaded file is returned, or None if the download file.\n", "\n", "The function signature is: \n", "\n", " `downloadFileUrl(url, saveFilename=None)`\n", "\n", "- `url (string)` the url to be accessed.\n", "- `saveFilename (string)` path to where the file must be saved (optional).\n", "\n", "The [`ryfiles.downloadUntar`](http://nelisw.github.io/pyradi-docs/_build/html/ryfiles.html#pyradi.ryfiles.downloadUntar) downloads and untars a compressed tar archive, and saves all files to the specified directory. The tarfilename is used to open the tar file, extracting to the destinationDir specified. If no destinationDir is given, the local directory '.' is used. Returns a list of the intarred filenames, or None of not succesful. Before downloading, a check is done to determine if the file was already downloaded and exists in the local file system.\n", "\n", "The file signature is: \n", "\n", " `downloadUntar(tgzFilename, url, destinationDir=None, tarFilename=None)`\n", "\n", "- `tgzFilename (string)` the name of the tar archive file.\n", "- `url (string)` url where to look for the file (not including the filename).\n", "- `destinationDir (string)` to where the files must be extracted (optional).\n", "- `tarFilename (string)` downloaded tar filename (optional).\n", "\n", "The [`ryfiles.untarTarfile`](http://nelisw.github.io/pyradi-docs/_build/html/ryfiles.html#pyradi.ryfiles.untarTarfile) untars a tar archive, and save all files to the specified directory. The tarfilename is used to open a file, extraxting to the saveDirname specified. If no saveDirname is given, the local directory '.' is used. Returns a list of filenames saved, or None if failed.\n", "\n", "The file signature is: \n", "\n", " `untarTarfile(tarfilename, saveDirname=None)`\n", "\n", "- `tarfilename (string)` the name of the tar archive.\n", "- `saveDirname (string)` to where the files must be extracted\n", "\n", "\n", "The [`ryfiles.unzipGZipfile`](http://nelisw.github.io/pyradi-docs/_build/html/ryfiles.html#pyradi.ryfiles.unzipGZipfile)\n", "unzips a file that was compressed using the gzip format. The zipfilename is used to open a file, to the saveFilename specified. If no saveFilename is given, the basename of the zipfilename is used, but with the file extension removed. Returns the filename of the saved file, or None if failed.\n", "\n", "The file signature is: \n", " \n", " `unzipGZipfile(zipfilename, saveFilename=None)`\n", "\n", "- `zipfilename (string)` the zipfilename to be decompressed.\n", "- `saveFilename (string)` to where the file must be saved (optional)." ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "The following example downloads a file, uncompress and untar the file, saving all the files in the tar archive to the destination directory (the current working directory in this case)." ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "run_control": { "frozen": false, "read_only": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "filesAvailable are ['arrayplotdemo.txt']\n" ] } ], "source": [ "tgzFilename = 'arrayplotdemo.tgz'\n", "destinationDir = '.'\n", "tarFilename = 'arrayplotdemo.tar'\n", "url = 'https://raw.githubusercontent.com/NelisW/pyradi/master/pyradi/data/'\n", "dlNames = ryfiles.downloadUntar(tgzFilename, url, destinationDir, tarFilename)\n", "\n", "if dlNames:\n", " print('filesAvailable are {}'.format(dlNames))" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "## Pulse detection: probability of detection and false alarm rate" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "The [`ryutils.detectSignalToNoise`](http://nelisw.github.io/pyradi-docs/_build/html/ryutils.html#pyradi.ryutils.detectSignalToNoise)\n", "function solves for the signal to noise ratio, given the threshold to noise ratio and probability of detection. Using the theory of matched filter design, calculate the signal to noise ratio, to achieve a required probability of detection. The function returns the signal to noise ratio required to achieve the probability of detection.\n", "Reference: \"Electro-optics handbook,\" Tech. Rep. EOH-11, RCA, 1974. RCA Technical Series Publication. \n", "\n", "When there is a signal present, the probability of detection (signal plus noise exceeds the threshold) is given by\n", "\n", "$$\n", "P_d =\n", "\\frac{1}{2}\\left[\n", "1+{\\rm erf}\\left(\n", "\\frac{i_s-i_t}{\\sqrt{2}i_n}\n", "\\right)\n", "\\right]\n", "$$\n", "\n", "where ${\\rm erf}$ is the error function:\n", "\n", "$$\n", "{\\rm erf}(z)=\\frac{2}{\\sqrt{\\pi}}\\int_0^z e^{-t^2}dt.\n", "$$\n", "\n", "\n", "\n", "The function signature is:\n", "\n", " `detectSignalToNoise(ThresholdToNoise, pD)` \n", "\n", "- `ThresholdToNoise (float)` the threshold to noise ratio (unitless).\n", "- `pD (float)` the probability of detection (unitless).\n", "\n", "The [`ryutils.detectThresholdToNoise`](http://nelisw.github.io/pyradi-docs/_build/html/ryutils.html#pyradi.ryutils.detectThresholdToNoise)\n", "function solve for threshold to noise ratio, given pulse width and FAR, for matched filter. Using the theory of matched filter design, calculate the threshold to noise ratio, to achieve a required false alarm rate. The function returns the threshold to noise ratio. Reference R. D. Hippenstiel, Detection Theory: Applications and Digital Signal Pro-cessing, CRC Press, 2002.\n", "\n", "The average false alarm rate is given by\n", "$$\n", "{FAR}=\n", "\\frac{1}{2 t_p \\sqrt{3}}\\exp^{-i_t^2/(2i_n^2)},\n", "$$\n", "where\n", "$t_p$ is the pulse width,\n", "$i_t$ is the threshold value, and\n", "$i_n$ is the rms noise value at the input to the threshold detector.\n", "\n", "\n", "The function signature is:\n", " `detectThresholdToNoise(pulseWidth, FAR)`\n", " \n", "- `pulseWidth (float)` the signal pulse width in [s].\n", "- `FAR (float)` the false alarm rate in [alarms/s].\n", "\n" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "run_control": { "frozen": false, "read_only": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "For a laser pulse with width=1e-07, a FAR=15 and Pd=0.999,\n", "the Threshold to Noise ratio must be 4.933070468777455\n", "and the Signal to Noise ratio must be 8.023302774945268\n", " \n" ] } ], "source": [ "import pyradi.ryutils as ryutils\n", "\n", "pulsewidth = 100e-9\n", "FAR = 15\n", "probDetection = 0.999\n", "ThresholdToNoise = ryutils.detectThresholdToNoiseTpFAR(pulsewidth,FAR)\n", "SignalToNoise = ryutils.detectSignalToNoiseThresholdToNoisePd(ThresholdToNoise, probDetection)\n", "print('For a laser pulse with width={0}, a FAR={1} and Pd={2},'.format(pulsewidth,FAR,probDetection))\n", "print('the Threshold to Noise ratio must be {0}'.format(ThresholdToNoise))\n", "print('and the Signal to Noise ratio must be {0}'.format(SignalToNoise))\n", "print(' ')" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "## Solving the range equation" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "It is frequently necessary to determine the operational detection distance of a source and sensor combination. The problem is usually stated as follows: \"What operating detection range can be achieved with a given source intensity, atmospheric attenuation, and sensor sensitivity?\" \n", "The inband irradiance is given by (Sec. 7.5)\n", "\n", "$$\n", "E_{\\cal S}=\\frac{v_{\\cal S} }{k \\,\\widehat{{\\cal R}}\\,Z_t\\,A_{1}}\n", "=\n", "\\int_{A_0} \\left(\n", "\\frac{1}{R_{01}^2}\n", "\\int_{0}^{\\infty}\n", " \\epsilon_{0\\lambda} L_{0\\lambda}\n", "\\,\\tau_{a\\lambda}(R) {\\cal S}_\\lambda\\,\n", "d\\lambda\\right) dA_{0}\\,\\cos\\theta_0\n", "$$\n", "\n", "The objective is to solve for $R$ where the irradiance $E_{{\\cal S}}=E_{{\\cal S}\\theta}$ is at the threshold at which the range must be determined. Normally, $E_{{\\cal S}\\theta}=SNR\\times NEE$ where the SNR is selected to yield a given probability of detection. \n", "\n", "Assuming a target smaller than the sensor field of view, and using the effective transmittance the above equation can be simplified to the form\n", "\n", "$$\n", "E_{{\\cal S}\\theta}=\\frac{I_{{\\cal S}} \\tau_{\\rm eff}(R)}{R^2}\n", "$$\n", "\n", "where $I_{{\\cal S}}$ is the source intensity. \n", "This equation can be solved is $\\tau_{\\rm eff}(R)$ is available in lookup table form.\n", "\n", "The [`ryutils.rangeEquation`](http://nelisw.github.io/pyradi-docs/_build/html/ryutils.html#pyradi.ryutils.rangeEquation) function solves the range equation for arbitrary transmittance vs range, with the equation given by\n", "\n", "$$\n", "E=\\frac{I_{{\\cal S}} \\tau_{\\rm eff}(R)}{R^n}\n", "$$\n", "\n", "where $E$ is the threshold irradiance in [W/m2], and $I$ is the intensity in [W/sr]. This range equation holds for\n", "the case where the target is smaller than the field of view.\n", "\n", "The range $R$ must be in [m], and $\\tau_a(R)$\n", "is calculated from a lookup table of atmospheric transmittance vs. range.\n", "The transmittance lookup table can be calculated from the simple Bouguer law,\n", "or it can have any abritrary shape, provided it decreases with increasing range.\n", "The user supplies the lookup table in the form of an array of range values and\n", "an associated array of transmittance values. The range values need not be on\n", "constant linear range increment.\n", "\n", "The parameter $n$\n", "\n", "* $n$=2 (default value) the general case of a radiating source smaller than the field of view.\n", "\n", "* $n$=4 the special case of a laser rangefinder illuminating a splot smaller than the sensor field of view, viewed against the sky. In this case there is an $R^2$ attenuation from the laser to the source and another $R^2$ attenuation from the source to the receiver, hence $R^4$ overall.\n", "\n", "The function signature is: \n", " `rangeEquation(Intensity, Irradiance, rangeTab, tauTab, rangeGuess = 1, n = 2)`\n", " \n", "- `Intensity (float or np.array[N,] or [N,1])` in [W/sr].\n", "- `Irradiance (float or np.array[N,] or [N,1])` in [W/m2].\n", "- `rangeTab (np.array[N,] or [N,1])` range vector for tauTab lookup in [m].\n", "- `tauTab (np.array[N,] or [N,1])` transmittance vector for lookup range in [m].\n", "- `rangeGuess (float)` starting value range estimate in [m] (optional).\n", "- `n (float)` range power (2 or 4) (optional).\n", "\n", "If the range solution is doubtful (e.g. not a trustworthy solution) the returned value is made negative. The following example attempts to solve the equation for three cases, only one of which provides a stable solution." ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "run_control": { "frozen": false, "read_only": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Range equation solver with irradiance threshold of 1e-99 W/m2:\n", " intensity of 200 W/sr at range -11381.128909552699 m the irradiance is 4.462603202968597e-07 W/m2, error is 4.46e+92. \n", "Check maximum range in lookup table\n", "\n", "Range equation solver with irradiance threshold of 1e-05 W/m2:\n", " intensity of 200 W/sr at range 3452.03409008682 m the irradiance is 1.000000000000012e-05 W/m2, error is 1.19e-14. \n", "\n", "\n", "Range equation solver with irradiance threshold of 1.0 W/m2:\n", " intensity of 200 W/sr at range -14.12716136943599 m the irradiance is 0.49750424774943763 W/m2, error is -5.02e-01. \n", "Check range resolution in lookup table\n", "\n" ] } ], "source": [ "import numpy as np\n", "from scipy.interpolate import interp1d\n", "\n", "rangeTab = np.linspace(0, 10000, 1000)\n", "tauTab = np.exp(- 0.00015 * rangeTab)\n", "Intensity=200\n", "Irradiancetab=[10e-100, 10e-6, 10e-1]\n", "for Irradiance in Irradiancetab:\n", " r = ryutils.rangeEquation(Intensity = Intensity, Irradiance = Irradiance, rangeTab = rangeTab,\n", " tauTab = tauTab, rangeGuess = 1, n = 2)\n", "\n", " #test the solution by calculating the irradiance at this range.\n", " tauTable = interp1d(rangeTab, tauTab, kind = 'linear')\n", "\n", " if np.abs(r[0]) < rangeTab[2]:\n", " rr = rangeTab[2]\n", " strError = \"Check range resolution in lookup table\"\n", " elif np.abs(r[0]) > rangeTab[-1]:\n", " rr = rangeTab[-1]\n", " strError = \"Check maximum range in lookup table\"\n", " else:\n", " rr = r[0]\n", " strError = \"\"\n", "\n", " irrad = Intensity * tauTable(rr) / rr ** 2\n", "\n", " print('Range equation solver with irradiance threshold of {} W/m2:'.format(Irradiance))\n", " print(' intensity of {4} W/sr at range {0} m the irradiance is {1} W/m2, error is {2:.2e}. \\n{3}\\n'.format(\n", " r[0],irrad, (irrad - Irradiance) / Irradiance, strError, Intensity))" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "The next example plots the detection range for a missile given the threshold irradiance, viewing an aircraft with signature of 200W/sr, through an atmosphere with an attenuation coefficient of 0.15 km-1. The graphs shows that when viewing the target through the atmosphere, an irradiance threshold (sensitivity) of 1 $\\mu$W/m$^2$ will provide a detection distance of 8 km. A missile with sensitivity of 0.1 $\\mu$W/m$^2$ provides a detection distance of 15 km. The example shows that you need ten times better sensitivity to double the detection distance. Increased detection distance comes at a tremendous cost!" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "run_control": { "frozen": false, "read_only": false } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwsAAAFUCAYAAABr4aCxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABPcUlEQVR4nO3dd5wV1f3/8dcbUGyoiNgABUuwYIuKWIMN0ajYxVgw+vtiTIzfxBY1MZaIUb8aE7+WfFHsXWxorFHXkqjYYyWiqBBQxI4FFT+/P86sDte7u3f37u5seT8fj3lw98ycmc89O3u5n5lz5igiMDMzMzMzK9Wl6ADMzMzMzKxtcrJgZmZmZmZlOVkwMzMzM7OynCyYmZmZmVlZThbMzMzMzKwsJwtmZmZmZlaWkwUzK5Sk2ZJWbIH9DpT0jKRPJB3W3Pu3tkPSUEnTio4DQNIBkh4pOo6GtJc42ypJJ0q6sug4zFqDkwWzZiTpDUmfZ1+A35Z0qaRFio6rrZBUI+n/5csiYpGIeL0FDnc0UBMRPSLinGp3JmmUpKckfSxpmqQzJHXLrV9C0s2SPpX0pqSflNTfStIrkj6T9ICkFeo4zt6SXiopu7eOsmNyPy/XHF+Ys3P2lGr305LHlxSSVm6tmFpSQ+dNybaDJN0taZak702SlP19fZF9/syWNKllo68zzvkljc8+D0PS0Cr21WYSQbPOysmCWfPbMSIWAdYB1gWOLTacTmsF4MWmVMwnATkLAb8ClgQ2BLYCjsytPw/4Elga2Ae4QNIa2f6WBG4CjgeWAJ4Erqvj8A8Cq0nqnYtlbWChkrKNgIdy9bYH7mqG91mVlthnc5PUtegYcuo8b8r4CrgeOKie/R2aJeCLRMTA5g21UR4B9gXeLjCGdqE9/M1Y5+ZkwayFRMTbwN2kpAEAScdIei3rGvOSpF1y6w6Q9IikMyV9IGmKpO1y6wdIeiir+3dJ5+Vvg0saIumfkj6U9Fx9V/Oyq9A3Sno3O85huXULZld3P8hiPCp/Za/0qm7+SrCknpJuz/b7Qfa6b7ZuDLAZcG521fPc0v1JWkzS5Vn9NyX9TlKXStqn5P3dD2yRO9YPKtj3PySdLel94MQyv88LIuLhiPgyIv4DXAVsktVfGNgNOD4iZkfEI8AEYL+s+q7AixFxQ0R8ke1/bUmrljnOdOB1YPOs6IekpOfBkrIupKSj1vbAHVk8v5H0n+xcmSRpq6z8xOyK75WSPgYOKGm30aQvrEdn7XZbVt7QeTtP20nqJek2pbswT0g6RbkuL5JWVboz8n4W3571Hb8kxtoE6blsm71y646QNFPSDEk/zZVfKukCSXdI+hTYQtJqSlfiP5T0oqSdctvPcwdMJV12JA3L4v5I0vmSHlTJHbMKz9OGzpt5RMSkiBhHE5Pg+kj6n+zva7GS3+mHkl6XtHFWPjVr41F17Sv7G/lz9n7mVnDsn0p6OTu/Xpd0cFa+MHAnsJy+u1uyXHYe35Cdx59Iej77Gz82i22qpGG5/S8naUJ2vk2W9F+5dYMlPZmdq+9I+lNW3l/ps2m0pOnZOXVESejzK32mfJKdQ+uXHLOuz9jv/R1m7T4uO85/sr+ZtpTUWifmZMGshSh9Sd4OmJwrfo30hXkx4CTgSknL5tZvCEwiXb0+AxgnSdm6q4GJQC/Sl81vv1BI6gP8DTiFdOX6SOBGZVeiS+LqAtwGPAf0IV0h/5WkbbNNTgBWypZtgTq/FJTRBbiEdFV/eeBz4FyAiPgt8DDfXfk8tEz9/yW1zYrAj4D9gZ/m1tfXPt+KiC1LjvXvCvf9OrAUMKaC97o5331p+wEwNztOreeA2ivEa2Q/18b3KelcqOsK8kN8lxhsnr2XR0rKHouILwEkzZeV3StpIHAosEFE9CD9Dt/I7XsEMB5YnJTwfCsixmZlZ2TttmO2qpLzNt925wGfAsuQzp9vz6HsC+C9pPN5KWBv4HxJa9Rz/HyMtW2wdrZN7R2aZbL4+pCuvJ8nqWeu6k+y2HoAj5P+Bu7JYvglcFXWdvVSuks0nnTHsBfpfNy4ZLOKzlMaPm+a4o9K3ZT+oQq6/0jqIulCYC1gWER8lHsP/yK9x6uBa4ENgJVJdwzOVfN1sZwJ7AAsSvqbPFvSD7O/k+2A6bm7JdOzOjsCVwA9gWdIF2a6kH7/JwP/l9v/NcA0YDlgd+BUZQk08BfgLxGxKOkz7/qS2LYAVgGGAcdI2jq3bidSuyxOSvJqL4A09BkL3/87vAz4mtS+62bHmycBNSuKkwWz5neLpE+AqaT/BE+oXZFdWZ4eEd9kX3JeBQbn6r4ZERdGxFzSfx7LAktLWp70H/Xvs6t2tVcga+0L3BERd2T7vpd01Xn7MvFtAPSOiJOzfb0OXAiMzNbvCYyJiPcjYipQcX//iHgvIm6MiM8i4hPSl7MfVVI3u4q2F3BsRHwSEW8AZzHvVday7dNM+54eEf8bEV9HxOcN7O+nwPrAmVnRIsBHJZt9RPpiWsn6Uvm7CJuRkoWHS8oezG2/OfBc1uZzge7A6pLmi4g3IuK13LaPRsQt2XlS7/usVcF5+23bkbrU7AackJ0HL5F+V7V2AN6IiEuytn4auJH0Ja4aXwEnR8RXEXEHMBvIf/m/NSL+ERHfkO72LQKclv0N3A/cTkpcGrI96S7RTdn7PYfvd7Wp9Dxt7HnRkN+QkuE+wFjgNkkr1bP9fKQv0kuQuk9+lls3JfsdzSV1metHat85EXEP6ffcLONGIuJvEfFaJA+SkrjNGqj2cETcnf0ObgB6k36fX5G+wPeXtLikfsCmwG8i4ouIeBa4iO/+9r8CVpa0ZHZ357GS45wUEZ9GxPOkCyH5c+SR7DN3LilxWTsrb+gzFnJ/h6QkaTvgV9mxZgJnl2xvVhgnC2bNb+fsiu5QYFXS1UUAJO0v6dns1v6HwKD8enJfOnL/cS9CuiL2fsl/5lNzr1cA9qjdb7bvTUlfUkqtQLqtn9/2OL77MrNcyb7frOhdp/e3kKT/U+rm8zHpCvniFd5OXxKYv+R4b5K++NSqq32aY99TqYCknYHTgO0iYlZWPJv0H37eosAnFa4v9RCwVnZlfAjpi8UrwLJZ2aZ8f7zCHQARMZk0tuJEYKakayUtl9u2oveZV8F5m99nb6BbSVnpubphyfm3D+nOQDXey7441vqMec+NfAzLAVOzL2q1Ss+Huszz9xERQbpqnVfpedrY86JeEfF4lgzPiYjLgH9Q/oJBrZVJV7hPqr1LlfNO7vXn2f5LyxaRtHyui9DspsQtaTtJj2XdhD7MYl6ygWqlsczKvrR/Gy/zfnbm2zT/uz6IdIfnFaUuczuUHKf0szD/t5RPEj8DFlAaf9DQZ2zpflcgJW4zctv/H+mul1nhnCyYtZDsCtmlZFeflZ5+cyGpi0iviFgceAEo1z2h1AxgCUkL5cr65V5PBa6IiMVzy8IRcVqZfU0lXTXMb9sjImq/VMwo2ffyJfU/Iw32rZX/kncE6Wruhtlt/dor4bXv8XtPcMmZRbrKl39K0PLAf+qpU6lK9l1fbABIGk76He6YXWms9W+gm6RVcmVr8103pRf57qpjbVeclaij73l2JXI6MBp4KyJqv4Q9mpUtAuSvgG5P6oZWW//qiNiU9H4DOL0R73Oe9RWet/k675K6U/TNlZWeqw+WnH+LRMQhFcbXVPn9Tgf6Zd1FauXPh0+p+xyfQe69Zd2L8u+1MRo6b6oV1P/58jKp28+dlXTBKnuAiLfiuy5Cje6WJKk76c7SmcDS2fl1B5V9ZlRiOumzM3+35tvfdUS8GhF7k76Ynw6Mz/4+a5V+Fk6nYQ19xsK872sqMAdYMrf9ohFRTXc0s2bjZMGsZf0Z2EbSOsDCpP8g3oVvu7IMqmQnEfEmqVvRiUqPJdyI1Ge31pXAjpK2ldRV0gJKjxws9yVmIvCx0iDYBbPtB0naIFt/PXCs0mDlvqT+3HnPAj/J6g1n3m5GPUhX9T6UtAS5LliZd0jdJMq9x7nZscdI6pF9ST08e29VaY59S9qS1Ld4t4iYWLL/T0lPOzpZ0sKSNiFdsb0i2+RmYJCk3SQtAPwe+Fd2t6AuD2cxPpwreyQre7K2C5GkAUD32n0pzS+xZfYl7AvS76PBQaY5pb+jRp23WVvfRDpXF1IaxL1/bpPbgR9I2k/SfNmygaTV6jh+JTE21uOkhODo7PhDSX9P12brnwV2zeJfmXmfPvQ3YE1JO2dXkX9BE++KVHDezEPJAqS7ZGR/592z14tnf/8LSOomaR9Ssn53AzFcQ7rq/fcGuiw1iqTuWayQBgIvkCVWpeYndZt7F/haaTD4sNz6d4BekhZrShyRulL+kzSWYwFJa5F+n1dlce4rqXd2l+nDrFr+7+X47DxYg5RY1fUUs7yGPmNLY5xB6np1lqRFlcaRrCSpoi6cZi3NyYJZC4qId4HLSU87eYnUT/5R0n+Aa5K6CVRqH9LjMt8jDWS+jnQ1qvY/xBGk//TfJV2pOooyf+PZl7kdSf22p5Cuul9EGhwKaQDrm9m6e/j+F5f/zup/mMV0S27dn4EFs30+xvcf5fkXYHelp8SUGwvxS9KXuNdJX4yvBi4us11TVLvv40ltdEeu28WdufU/J733maR+4IdExIvw7XmwG2kMxwekwaMN9Ud+kHS1Mz9x1sNZWb4L0o/JuiBlupO6Sc0idZNYinReVGocabzDh5JuaeJ5eyiprd4mnT/X8N25+gnpy+BI0lXat0lXdLuXO34d+z8RuCzbZs9GvDeyGL4kDU7djtRO5wP755K3s0l98t8hjTm4Kld3FrAHaeDye8DqpER+TmPjyNR53uS6+NTe3VuBlPzV3nn4nDSQGlI3llNIf/+zSOf7zhHR4FwLWZelk4H7JfVv4vsoNSmLrw8pYfmcee/s1R77E+AwUjL/AWkg+oTc+ldI7fJ69vternQfFdgb6E86324mjae5N1s3HHhRqQvVX4CRkZ5YVutB0kMq7gPOzMZr1KuCz9hy9iclTi+R2mE85buRmrU6pe6WZtbeSLoOeCUiSq/eN/dxhgJXRkRTu1pYC5J0B3BupEG9bZKk04FlIqIxT9ZqF7KuTNOAfSLigaLjseaRJU1TgPlKxsKYdTq+s2DWTmRdNVbKblEPJ91JuKXgsKx4NUCb+pKqNI/CWlm3mcGkbh83Fx1Xc8m6+yyedQE6jtS/vvQpOmZmHYJnDTRrP5Yh9W/uRbqSeUhEPFNsSFa0iDij6BjK6EHqOrIcqXvNWcCthUbUvDYidWOr7Tayc1T4GFozs/bG3ZDMzMzMzKysVumGJOlipSnYX8iVraP0XOVnlaZaH5xbd6zSlOyTlJvxUNJ6StO6T5Z0Tu2TFbKnLlyXlT9e1wCtauubmZmZmXUmrTVm4VLSEwfyziBNBLMO6TGCZwBIWp30lIw1sjrn67sJnS4gPWN8lWyp3edBwAcRsTLpKRb5Z4rnVVvfzMzMzKzTaJUxCxHxUJmr9cF3M1cuxncTnYwAro2IOcAUSZOBwZLeABaNiEcBJF0O7AzcmdU5Mas/HjhXkiLXx0rSstXUL2fxxRePlVdultnurQk+/fRTFl544YY3tGbnti+W279Ybv/iuO2L5fYv1lNPPTUrInq39nGLHOD8K+BuSWeS7nBsnJX3Yd6nSkzLyr7KXpeW19aZChARX0v6iDQIdFZu+z5V1gdA0mjS3Ql69+7NmWeeWfEbtuY1e/ZsFlmk0ROGWjNw2xfL7V8st39x3PbFcvsXa4sttniziOMWmSwcAvw6Im7MJtUZB2xN+anp65qyvvbKf33rqGCbSuqnwoixwFiAgQMHxtChQ8ttZq2gpqYGt38x3PbFcvsXy+1fHLd9sdz+nVOR8yyMIj0GEuAGoHaA8zSgX267vqQuStOy16Xl89SR1I3Uren9kuNVW9/MzMzMrFMpMlmYDvwoe70l8Gr2egIwMntC0QDSQOSJETED+ETSkOwpRvvz3XO7J5CSD4DdgftLxxtUW9/MzMzMrLNplW5Ikq4BhgJLSpoGnAD8F/CX7Er+F2TjACLiRUnXkya6+Rr4RUTMzXZ1COnJSguSBibfmZWPA67IBkO/T3qaUu2xn82euNSk+mZmZmZmnVVrPQ1p7zpWrVfH9mOAMWXKnwQGlSn/Atijjn2tU019MzMzM7POqshuSGZmZmZm1oY5WTAzMzMzs7KcLJiZmZmZWVlOFszMzMzMrCwnC2bt2Ndfw1lnwezZRUdiZmZmHZGTBbN27Omn4eij4cADwTODmJmZWXNzsmDWjg0eDKeeCjfcAH/6U9HRmJmZWUfjZMGsnTv6aNhtt/Tv/fcXHY2ZmZl1JE4WzNo5CS65BH7wA9hrL5g6teiIzMzMrKNwsmDWAfToATffDHPmpLsMX3xRdERmZmbWEThZMOsgVl0VLrsMnngCDjus6GjMzMysI3CyYNaB7LILHHssXHhhWszMzMyq4WTBrIP5wx9gm23g0ENh4sSiozEzM7P2zMmCWQfTtStccw0su2wavzBzZtERmZmZWXvlZMGsA+rVC266CWbNgpEj00zPZmZmZo3lZMGsg/rhD+Gvf4UHHkjjGMzMzMway8mCWQc2ahT8/Odw5plplmczMzOzxnCyYNbBnX02bLQR/PSn8OKLRUdjZmZm7YmTBbMObv75012FRRZJj1b96KOiIzIzM7P2wsmCWSfQp09KGKZMgf33h2++KToiMzMzaw+cLJh1EpttBmedBRMmwB//WHQ0ZmZm1h60SrIg6WJJMyW9UFL+S0mTJL0o6Yxc+bGSJmfrts2Vryfp+WzdOZKUlXeXdF1W/rik/nXEUVV9s/bul7+En/wEjj8e7rqr6GjMzMysrWutOwuXAsPzBZK2AEYAa0XEGsCZWfnqwEhgjazO+ZK6ZtUuAEYDq2RL7T4PAj6IiJWBs4HT64ij2vpm7ZoEY8fCmmumpOH114uOyMzMzNqyVkkWIuIh4P2S4kOA0yJiTrZN7TyzI4BrI2JOREwBJgODJS0LLBoRj0ZEAJcDO+fqXJa9Hg9sVXvXoFa19c06ioUXThO2RaQZnj/7rOiIzMzMrK0qcszCD4DNsm4/D0raICvvA0zNbTctK+uTvS4tn6dORHwNfAT0KjletfXNOoyVVoKrr4bnnoOf/SwlDmZmZmaluhV87J7AEGAD4HpJKwLlruhHPeU0sI4KtqmkftpQGk3qykTv3r2pqakpt5m1gtmzZ7v9q7DggjBq1ApceukAevZ8lV12+U/Fdd32xXL7F8vtXxy3fbHc/p1TkcnCNOCmrEvQREnfAEtm5f1y2/UFpmflfcuUk6szTVI3YDG+3+2p2voARMRYYCzAwIEDY+jQoRW+XWtuNTU1uP2rs/nmMGsWnH/+Kuyxxypsumll9dz2xXL7F8vtXxy3fbHc/p1Tkd2QbgG2BJD0A2B+YBYwARiZPaFoAGkg8sSImAF8ImlINp5gf+DWbF8TgFHZ692B+7Mk5FvV1jfriLp0gSuugP79YY89YMaMoiMyMzOztqS1Hp16DfAoMFDSNEkHARcDK2aPU70WGBXJi8D1wEvAXcAvImJutqtDgItIg55fA+7MyscBvSRNBg4Hjskd+9lcKI2ub9bRLb443HwzfPxxShi+/LLoiMzMzKytaJVuSBGxdx2r9q1j+zHAmDLlTwKDypR/AexRx77Wqaa+WWcwaBCMGwd77w1HHAH/+79FR2RmZmZtQZFjFsysDRk5Ep54Av70Jxg8GPbbr+iIzMzMrGhFjlkwszbm9NNh6FAYPRqeeaboaMzMzKxoThbM7FvdusF110GvXrDrrvB+2WeCmZmZWWfhZMHM5rHUUnDjjTB9OvzkJzB3bsN1zMzMrGNysmBm37PhhmmQ8913w4knFh2NmZmZFcXJgpmV9V//BQceCKecArfe2vD2ZmZm1vE4WTCzsiQ47zxYf/30ZKRJk4qOyMzMzFqbkwUzq9MCC6TxC927pwHPs2cXHZGZmZm1JicLZlav5ZeHa6+FV15J3ZIiio7IzMzMWouTBTNr0FZbwR//CDfcAGedVXQ0ZmZm1lqcLJhZRY46CnbbDX7zG3j66cWLDsfMzMxagZMFM6uIBJdcAgMHwsknr87UqUVHZGZmZi3NyYKZVaxHD7j5Zvjqqy7stht88UXREZmZmVlLcrJgZo0ycCAce+wrPPEE/PKXRUdjZmZmLcnJgpk12qabzuK44+Cii+DCC4uOxszMzFqKkwUza5KTT4Zhw+DQQ2HixKKjMTMzs5bgZMHMmqRrV7j6alhuufSUpJkzi47IzMzMmpuTBTNrsl690gzPs2bBXnvB118XHZGZmZk1JycLZlaVH/4Q/vpXqKmBY48tOhozMzNrTt2KDsDM2r9Ro+CJJ+DMM2GDDWDPPYuOyMzMzJqD7yyYWbP4059g443hwAPhhReKjsbMzMyag5MFM2sW888PN9wAiywCu+4KH31UdERmZmZWrVZJFiRdLGmmpO9db5R0pKSQtGSu7FhJkyVNkrRtrnw9Sc9n686RpKy8u6TrsvLHJfWvI46q6ptZ/ZZbLiUMU6bA/vvDN98UHZGZmZlVo7XuLFwKDC8tlNQP2AZ4K1e2OjASWCOrc76krtnqC4DRwCrZUrvPg4APImJl4Gzg9DriqLa+mTVgs81Sl6QJE+DUU4uOxszMzKrRKslCRDwEvF9m1dnA0UDkykYA10bEnIiYAkwGBktaFlg0Ih6NiAAuB3bO1bksez0e2Kr2rkGtauubWeUOPRT22Qd+/3u4886iozEzM7OmKmzMgqSdgP9ExHMlq/oAU3M/T8vK+mSvS8vnqRMRXwMfAb3K7Lea+mZWIQnGjoW11oKf/ARef73oiMzMzKwpCnl0qqSFgN8Cw8qtLlMW9ZTXV6eS/VZaP20ojSZ1ZaJ3797U1NSU28xawezZs93+Bam07Y8+egEOPng9hg2bw7nnPs0CC3gQQ3PwuV8st39x3PbFcvt3TkXNs7ASMAB4Luvt0xd4WtJg0hX/frlt+wLTs/K+ZcrJ1ZkmqRuwGN/v9lRtfQAiYiwwFmDgwIExdOjQit6wNb+amhrc/sVoTNv37Ak//vF8XHnl5lxxRbrrYNXxuV8st39x3PbFcvt3ToV0Q4qI5yNiqYjoHxH9SV/WfxgRbwMTgJHZE4oGkAYiT4yIGcAnkoZk4wn2B27NdjkBGJW93h24PxuXkD9mVfXNrGm22w5OOgmuugrOPbfoaMzMzKwxWuvRqdcAjwIDJU2TdFBd20bEi8D1wEvAXcAvImJutvoQ4CLSoOfXgNqhk+OAXpImA4cDx+SO/Wxu942ub2bV++1vYccd4fDD4ZFHio7GzMzMKtUq3ZAiYu8G1vcv+XkMMKbMdk8Cg8qUfwHsUce+16mmvplVr0sXuOIK2GAD2GMPeOqpNCeDmZmZtW2ewdnMWsVii8FNN8HHH6eE4csvi47IzMzMGuJkwcxazaBBcPHF8M9/pi5JZmZm1rY5WTCzVrXXXnDEEXDeeXD55UVHY2ZmZvVxsmBmre6002DoUDj4YHjmmaKjMTMzs7o4WTCzVtetG1x3HSy5JOy6K7z3XtERmZmZWTlOFsysEEstBePHw/TpsM8+MHduw3XMzMysdTlZMLPCbLhhmqjt7rvhhBOKjsbMzMxKOVkws0L913/BQQfBmDFwyy1FR2NmZmZ5dU7KJunACvfxdUT4mSZm1mTnngvPPQf77w9PPAEDBxYdkZmZmUH9MziPBR6uYB8bAE4WzKzJFlgAbrwR1lsvDXh+7DHo0aPoqMzMzKy+ZOHziNiioR1I+qAZ4zGzTmr55dMTkrbZBg48EK6/HqSiozIzM+vc6huz8MMK97FBcwRiZrbllmkOhvHj4cwzi47GzMzM6kwWIuLVSnYQEZObLxwz6+yOPBJ23x2OOQbuv7/oaMzMzDq3+rohfUtSN2BvYF1gkfy6iBjdAnGZWSclwcUXw0svwV57wVNPpS5KZmZm1voqfXTqlcAxwDfAOyWLmVmz6tEDbroJvvwSdtsNvvii6IjMzMw6p4ruLADDgX4R8UlLBmNmVmvgQLj8cth5Zzj0ULjooqIjMjMz63wqvbPwErBESwZiZlZqxAj47W9h3Di48MKiozEzM+t8Kr2zsC9wkaR7KOl65AnZzKwlnXQSPPlkuruw1lqw4YZFR2RmZtZ5VJosHABsBvQEPs+VB56QzcxaUNeucPXVacK23XZLA56XXrroqMzMzDqHSpOF/wbWjYiXWzIYM7NyllgiDXjeeGMYORLuvRe6VfrpZWZmZk1W6ZiFd4C3WjIQM7P6rLsujB0LNTVpDgYzMzNreZVemzsbuErSacDM/IqIeL3ZozIzK2O//WDiRDjrLNhggzQPg5mZmbWcSu8snAfsBPwTmJxbKprlWdLFkmZKeiFX9j+SXpH0L0k3S1o8t+5YSZMlTZK0ba58PUnPZ+vOkaSsvLuk67LyxyX1ryOOquqbWfHOOit1RzrwQHjhhYa3NzMzs6arKFmIiC51LF0rPM6lpLka8u4FBkXEWsC/gWMBJK0OjATWyOqcL6n2OBcAo4FVsqV2nwcBH0TEyqS7IKfXEUe19c2sYPPPDzfcAIsuCrvuCh9+WHREZmZmHVeldxaqEhEPAe+XlN0TEV9nPz4G9M1ejwCujYg5ETGFdAdjsKRlgUUj4tGIqH0K0865Opdlr8cDW9XeNahVbX0zazuWWy4lDFOmwP77wzffFB2RmZlZx1RnsiDpsrrWlWx3STPEcSBwZ/a6DzA1t25aVtYne11aPk+dLAH5COhVcoxq65tZG7LppvCnP8Ftt8GYMUVHY2Zm1jHVN8B5d0mXAg1dYd8F+GlTA5D0W+Br4KraojKbRT3l9dWZ51BV1k8bSqNJXZno3bs3NTU15TazVjB79my3f0HaStsPGgTbbLMqJ5ywNPPP/zwbbvh+w5U6gLbS/p2V2784bvtiuf07p/qShXeBiyvYx9tNPbikUcAOwFZZ1yBIV/z75TbrC0zPyvuWKc/XmSapG7AYJd2emqE+ABExFhgLMHDgwBg6dGglb9VaQE1NDW7/YrSltt9wwzTg+bTT1uLJJ2GllYqOqOW1pfbvjNz+xXHbF8vt3znV2Q0pIvpHxIAKllWbcmBJw4HfADtFxGe5VROAkdkTigaQBiJPjIgZwCeShmTjCfYHbs3VGZW93h24P5d81L6fquqbWdu00EJpwjYpDXj+7LOG65iZmVllWmWAs6RrgEeBgZKmSToIOBfoAdwr6VlJfwWIiBeB64GXgLuAX0TE3GxXhwAXkQY9v8Z34xzGAb0kTQYOB76dsknSs7lQGl3fzNq+FVeEq6+G55+H0aPBqb6ZmVnzqHRStqpExN5lisfVs/0Y4HtDFiPiSWBQmfIvgD3q2Nc61dQ3s/Zh+HA4+WQ4/vjUNemXvyw6IjMzs/avVe4smJm1huOOg512gsMPh4cfLjoaMzOz9s/Jgpl1GF26wOWXw4ABsMceMH16w3XMzMysbhUnC5JWk3S8pPOyn1eVtFbLhWZm1niLLQY33wyzZ8Puu8OXXxYdkZmZWftVUbIgaQ/gQdLkZftlxYsAf2qhuMzMmmyNNeDii+HRR1OXJDMzM2uaSu8snAwMi4ifAbVPJnoOWLtFojIzq9Kee8KRR8J558FlFc1Hb2ZmZqUqTRaWIiUH8N3MxkEdsxybmbUFf/wjbLEF/Oxn8MwzRUdjZmbW/lSaLDzFd92Pao0EJjZvOGZmzadbN7j2WlhyyTRh23vvFR2RmZlZ+1JpsnAYcIqkB4GFJd0N/AH4dYtFZmbWDJZaCm68MT0Zaa+9YM6coiMyMzNrPypKFiLiFWBV4Dzgd8AlwJoR8WoLxmZm1iwGD4axY+G++2C//WDu3IbrmJmZWSNmcI6Iz4DrWzAWM7MWM2oUzJqVBj337Al//StIRUdlZmbWtlWULEh6mPKDmecA04CbIuK25gzMzKy5HXEEvP8+nHoqLLFEGgBtZmZmdat0zEIN0J8018KV2b8rAE8C7wAXSzq6BeIzM2tWp5wCBx8Mp50G//M/RUdjZmbWtlXaDWkYsG1EvFxbIOkq4LKI2FDSTcC1wBktEKOZWbOR0twLH34IRx+duiT9v/9XdFRmZmZtU6XJwqrA6yVlbwIDASJioqSlmjMwM7OW0rUrXH45fPRRusvQsyfstlvRUZmZmbU9lXZDegi4RNLKkhaQtDJwIfAIgKQ1gRktFKOZWbObf34YPx6GDIGf/ATuvbfoiMzMzNqeSpOFUdm2LwGfAi8CXYEDsvVfAns3d3BmZi1p4YXh9tth1VVhl13gsceKjsjMzKxtqXSehfcjYiSwALAcsGBE7B0Rs7L1kyLiyRaM08ysRfTsCXffDcssA9tvDy+8UHREZmZmbUeldxaQtBiwPrAGMFTSlpK2bLHIzMxayTLLpG5ICywAw4bBlClFR2RmZtY2VDrPwgGk2ZtnA5/lVgWwYvOHZWbWugYMgHvugc03h222gUceSUmEmZlZZ1bpnYUxwO4RsXREDMgtThTMrMMYNAjuuAPefhu23RY++KDoiMzMzIpVabLQDbinJQMxM2sLhgyBW26BV16BHXaATz8tOiIzM7PiVJosnA78TlLFYxzMzNqrrbeGq69OT0fafXf48suiIzIzMytGpV/+fw38DvhE0lv5pZLKki6WNFPSC7myJSTdK+nV7N+euXXHSposaZKkbXPl60l6Plt3jiRl5d0lXZeVPy6pfx1xVFXfzDqP3XaDsWPhrrtg//1h7tyiIzIzM2t9lSYL+wJbA9sD+5UslbgUGF5SdgxwX0SsAtyX/Yyk1YGRpKcuDQfOl9Q1q3MBMBpYJVtq93kQ8EFErAycTboTUk619c2sEznoIDjjDLjuOvjFLyCi6IjMzMxaV0VPQ4qIB6s5SEQ8VOZq/QhgaPb6MqAG+E1Wfm1EzAGmSJoMDJb0BrBoRDwKIOlyYGfgzqzOidm+xgPnSlLEd/+1S1q2mvpm1jkddRS8/z6cdhr06gVjxhQdkZmZWeupKFkAkLQOsBmwJKDa8oj4fROPvXREzMj2MUPSUll5HyA/j+q0rOyr7HVpeW2dqdm+vpb0EdALmJXbvk+V9c2skzr11JQwnHoqLLEEHHFE0RGZmZm1jkrnWRhN6p5zD7Ad6Wr8MODWFohJZcqinvL66lSy30rrpw1TW4wG6N27NzU1NeU2s1Ywe/Zst39BOmPb77kn/Pvfq3PkkUvxzjuvsP32bxcWS2ds/7bE7V8ct32x3P6dU6V3Fo4GhkfEw5I+iIhdJG1HGlvQVO9IWja7q7AsMDMrnwb0y23XF5ielfctU56vM01SN2Ax4P2S41VbH4CIGAuMBRg4cGAMHTq0ojdrza+mpga3fzE6a9tvthnsuCOcddaqDBmyKrvuWkwcnbX92wq3f3Hc9sVy+3dOlQ5wXioiHs5efyOpS0TcCexYxbEnAKOy16P47i7FBGBk9oSiAaSByBOzLkufSBqSPcVo/5I6tfvaHbi/dLxBtfXNzOafH266CTbcEPbeG/7+96IjMjMza1mVJgvTcgOU/w2MkLQZUNHTxyVdAzwKDJQ0TdJBwGnANpJeBbbJfiYiXgSuB14C7gJ+ERG1Dy08BLgImAy8RuoOBTAO6JUNhj6c7MlK2bGfzYXS6PpmZnkLLwx/+xsMHAg77wyPP150RGZmZi2n0m5IZwCrAW8AJ5OeGDQ/cFgllSNi7zpWbVXH9mOA7z1zJCKeBAaVKf8C2KOOfa1TTX0zs1I9e8Ldd8Omm8L228NDD8EaaxQdlZmZWfNr8M5C1mXnIeBegKz7UU+gZ0Rc0LLhmZm1TcsuC/feC927w7Bh8MYbRUdkZmbW/BpMFrK++88D3+TKvoyI2S0ZmJlZW7fiinDPPfD557DNNvB2cQ9IMjMzaxGVjll4BvhBSwZiZtYeDRqUxjBMnw7Dh8OHHxYdkZmZWfOpdMxCDXCXpEtJk5d9+6SgiLi4+cMyM2s/NtoIbr4ZdtghLXfdBYssUnRUZmZm1as0WdgEmAL8qKQ8ACcLZtbpDRsGV18NI0fCdtvBHXdAjx5FR2VmZladipKFiNiipQMxM2vvdt8drrkmzcEwfDjceScsumjRUZmZmTVdpWMWzMysAnvsAdddBxMnpoTh44+LjsjMzKzpnCyYmTWz3XZLCcMTT8C228JHHxUdkZmZWdM4WTAzawG77go33ABPPZXGM/gpSWZm1h45WTAzayE77wzjx8Mzz6R5GD74oOiIzMzMGqfSpyEhaTFgIDDPAwEj4v7mDsrMrKPYaSe46abUNWmbbdKszz17Fh2VmZlZZSpKFiQdAJwHzAY+y60KYMXmD8vMrOPYYYeUMOy6K2y9dUoYllii6KjMzMwaVmk3pDHA7hGxdEQMyC1OFMzMKvDjH8Mtt8CLL8JWW8F77xUdkZmZWcMqTRa6Afe0ZCBmZh3ddtvBrbfCyy+nhGHWrKIjMjMzq1+lycLpwO8keUC0mVkVtt0WJkyASZNgyy3h3XeLjsjMzKxulX75/zXwO+ATSW/llxaMzcysQxo2DG67DV59NSUMM2cWHZGZmVl5lT4Nad8WjcLMrJPZemv429/S4OcttoD774elly46KjMzs3lVlCxExIMtHYiZWWez5ZZwxx1p8HNtwrDMMkVHZWZm9p2KuiFJmk/SSZJel/RF9u9JkuZv6QDNzDqyoUNTwvDmmylhmDGj6IjMzMy+U+mYhTOArYGfAWtn/25JGvhsZmZV+NGP4M47YepU2HxzeOONoiMyMzNLKk0W9gB2ioh7ImJSRNwD7ALs2XKhmZl1HptvniZrmzULNt00PV7VzMysaJUmC2pkuZmZNdJGG8GDD8LXX8Nmm8GTTxYdkZmZdXaVJgs3ALdJ2lbSapKGA7cA11cbgKRfS3pR0guSrpG0gKQlJN0r6dXs35657Y+VNFnSJEnb5srXk/R8tu4cSWUTmWrrm5m1pLXWgkcegR490gDoB/14CTMzK1ClycLRwN+B84CngHOBB4DfVHNwSX2Aw4D1I2IQ0BUYCRwD3BcRqwD3ZT8jafVs/RrAcOB8SV2z3V0AjAZWyZbhZY5XVX0zs9aw8sopYejbF4YPh9tvLzoiMzPrrCpKFiLiy4j4fUSsHBELZf8eHxFzmiGGbsCCkroBCwHTgRHAZdn6y4Cds9cjgGsjYk5ETAEmA4MlLQssGhGPRkQAl+fq5FVb38ysVfTpAw89BGuuCbvsAldfXXREZmbWGdU5z4KkzSPioez1lnVtFxH3N/XgEfEfSWcCbwGfA/dExD2Slo6IGdk2MyQtlVXpAzyW28W0rOyr7HVpealq65uZtZoll4T77oMRI2DffeGjj+CQQ4qOyszMOpP6JmU7HxiUvR5XxzYBrNjUg2djEUYAA4APgRsk1TdbdLlxBFFPeXPXR9JoUnclevfuTU1NTdlAreXNnj3b7V8Qt33rOuaYLsyZszo///mSPPXU64wY4fYvks//4rjti+X275zqTBayMQS1rwe00PG3BqZExLsAkm4CNgbekbRsdldhWWBmtv00oF+ufl9St6Vp2evS8lLV1icixgJjAQYOHBhDhw5t+F1ai6ipqcHtXwy3fevbYgs48EAYN25FZs/uxjXXLI8fw1AMn//FcdsXy+3fOVU6g/OtdZTfVOXx3wKGSFooe/rQVsDLwARgVLbNKKD2+BOAkZK6SxpAGog8Meuy9ImkIdl+9s/Vyau2vplZIeabDy67DA49FK67bnlGj4a5c4uOyszMOrr6uiHlbVFH+dBqDh4Rj0saDzwNfA08Q7pqvwhwvaSDSAnFHtn2L0q6Hngp2/4XEVH73+UhwKXAgsCd2YKknUhPW/p9U+qbmbUVXbrAOefARx+9wUUX9efDD+GKK2CBBYqOzMzMOqp6kwVJJ2cv58+9rrUi8Ga1AUTECcAJJcVzSHcZym0/BhhTpvxJvhtjkS+fQLqj0KT6ZmZtiQQHHvgG667bn8MPh3ffhVtugcUXLzoyMzPriBrqhtQvW7rkXvcj9emfSnbF38zMWtevfw3XXAP//CdsuilMnVp0RGZm1hHVe2chIn4KIOmfEXFh64RkZmaVGDkSll4adt4ZNtoI7rwzzctgZmbWXCqdwXmOpLXyBZLWlrRfC8RkZmYV2mKLNNszwGabgZ9qaGZmzanSZOEPpG5HeVOBU5o3HDMza6w114RHH4W+fWHbbeG664qOyMzMOopKk4VFgY9Lyj4CFm/WaMzMrEn69YOHH4YhQ1L3pD/9qeiIzMysI6g0WXgJ2K2kbBfSnAhmZtYG9OwJd98Nu+8ORxwBhx8O33xTdFRmZtaeVTrPwm+AOyTtBbwGrEx6tOn2LRWYmZk13gILpG5Ihx8OZ58N06alydwWXLDoyMzMrD2q6M5CRDxCmoPgCWBhYCIwKCL+0YKxmZlZE3TpkhKFM8+E8ePTIOh33ik6KjMza48qvbNARLwl6Qxg6YiY0YIxmZlZlaTUFWnFFWGffdJYhttvhzXWKDoyMzNrTyq6syBpcUlXA18Ak7OynST5aUhmZm3YLrvAQw/BF1/AxhvDvfcWHZGZmbUnlQ5w/ivp6UcrAF9mZY8Ce7VEUGZm1nzWXx8mToT+/WG77eBCT7FpZmYVqjRZ2Ao4LOt+FAAR8S6wVEsFZmZmzadfvzR527BhMHo0HHWUn5RkZmYNqzRZ+AhYMl8gaXnAYxfMzNqJHj1gwgQ49NA0+Hm33WD27KKjMjOztqzSZOEi4EZJWwBdJG0EXEbqnmRmZu1Et27wv/8Lf/lLShw23hhef73oqMzMrK2qNFk4HbgeOA+YD7gYuBX4SwvFZWZmLeiww+Cuu9I8DBtsAPffX3REZmbWFlWaLCwdEX+OiNUjYuGIWC0i/gws3YKxmZlZC9pmmzTweZll0liGc86BiKKjMjOztqTSZOHfdZS/1FyBmJlZ61t5ZXjsMdhhB/jv/4aDDoI5c4qOyszM2opKkwV9r0BaFPCzNMzM2rkePeCmm+D3v4dLLoGhQ2GGH19hZmY0kCxImirpLWBBSW/lF9KTkG5pjSDNzKxldekCJ50EN94Izz+f5mZ4/PGiozIzs6J1a2D9vqS7CncA++XKA3gnIia1VGBmZtb6dt0VVlkFRoyAzTeHsWNh1KiiozIzs6LUmyxExIMAkpaMiM9aJyQzMyvSmmvCE0/AnnvCAQfAs8/C//xPeuyqmZl1LpWOWZgraYyk1yV9BCBpmKRDWzA2MzMrSK9ecPfd8KtfwZ//DMOHw3vvFR2VmZm1tkqThT8Dg4B9SF2QAF4EDqk2AEmLSxov6RVJL0vaSNISku6V9Gr2b8/c9sdKmixpkqRtc+XrSXo+W3eOpO8Nym6O+mZmnUW3bnD22WnQ88MPp3EMTz5ZdFRmZtaaKk0WdgZ+EhGPkj0BKSL+A/Rphhj+AtwVEasCawMvA8cA90XEKsB92c9IWh0YCawBDAfOl9Q1288FwGhglWwZXnqgauubmXVGBxyQkoVvvoFNNoG//tXzMZiZdRaVJgtfUjK+QVJvoKqb0tnjVzcHxgFExJcR8SEwArgs2+wyUrJCVn5tRMyJiCnAZGCwpGWBRSPi0YgI4PJcnbxq65uZdUqDB8PTT8NWW8Ehh8B++8GnnxYdlZmZtbRKk4UbgMskDQDIvlyfC1xb5fFXBN4FLpH0jKSLJC1MmjF6BkD271LZ9n2Aqbn607KyPtnr0vJS1dY3M+u0evWC22+HU06Ba65JCcQrrxQdlZmZtaRKn21xHHAG8DywEPAqcCFwUjMc/4fALyPicUl/IetyVIdy4wiinvLmro+k0aTuSvTu3ZuampqygVrLmz17ttu/IG77YhXd/ptsAmecsTinnLI6667blaOOmsSWW84sLJ7WVnT7d2Zu+2K5/TunipKFiPgS+BXwq6z70aysu061pgHTIqJ26p/xpGThHUnLRsSM7C7GzNz2/XL1+wLTs/K+ZcrLHa+a+kTEWGAswMCBA2Po0KENvEVrKTU1Nbj9i+G2L1ZbaP+hQ2HkSNhrL/jDH1bn/fdX56yzoHv3QsNqFW2h/Tsrt32x3P6dU0XdkCStLulgSccCuwKrNcfBI+JtYKqkgVnRVsBLwASgdhqgUcCt2esJwEhJ3bMuUasAE7OuSp9IGpI9xWj/XJ28auubmVmmTx944AE44gg47zzYbDN4882iozIzs+ZU752F7IvzONIX9mmkq+19gOUkXQEc2Ax3GH4JXCVpfuB14KekJOZ6SQcBbwF7AETEi5KuJyUUXwO/iIi52X4OAS4FFgTuzBYk7QSsHxG/b0p9MzOr23zzwZlnwsYbw09/CuuuC+PGwS67FB2ZmZk1h4a6IY0GhgJDIuKJ2kJJGwDXAAcDf60mgIh4Fli/zKqt6th+DDCmTPmTpLkgSssnkO4oNKm+mZk1bNddYe21U9ekXXdNT0w66yxYcMGiIzMzs2o01A1pP+CwfKIAkP38q2y9mZkZK60E//gHHHUUXHABbLghvPRS0VGZmVk1GkoWVgcerGPdg9l6MzMzAOafH844A+66C955J836PHasJ3EzM2uvGkoWukbEJ+VWZOWVztNgZmadyLbbwnPPwaabwsEHp6cmffhh0VGZmVljNTRmYT5JW1B+HoJK6puZWSe1zDLpDsOZZ8JvfwsTJ8IVV6SnJpmZWfvQ0Jf9mcDFDaw3MzMrq0sXOPpo+NGPYJ990r9HHw0nndQ55mQwM2vv6k0WIqJ/K8VhZmYd2IYbwrPPpjkZTj893XG48koY5GfQmZm1aR5zYGZmrWKRReD//g9uuw1mzID11kuPV/3mm6IjMzOzujhZMDOzVrXDDvDCC7D99nDkkbDllp752cysrXKyYGZmra53b7jpJrjkEnj6aVhrLbj8cj9i1cysrXGyYGZmhZDggAPgX/9Ksz+PGgW77w7vvlt0ZGZmVsvJgpmZFap/f3jggTSZ2+23w+qrw3XX+S6DmVlb4GTBzMwK17UrHHVU6pI0YACMHAm77QZvv110ZGZmnZuTBTMzazPWWAP++c90l+GOO9JdBo9lMDMrjpMFMzNrU7p1S3cZnnsuJQujRqUnKE2bVnRkZmadj5MFMzNrkwYOhAcfhL/8BWpq0l2HsWM9L4OZWWtysmBmZm1W165w2GHw/PNpEreDD4bNN4cXXyw6MjOzzsHJgpmZtXkrrgj33ZfmZXj5ZVhnHfjtb+Hzz4uOzMysY3OyYGZm7ULtvAyvvAL77AOnngprrgn33lt0ZGZmHZeTBTMza1d694ZLL4X770/dlIYNS8nDzJlFR2Zm1vE4WTAzs3Zpiy3SE5NOOAHGj4dVV4ULL/QAaDOz5uRkwczM2q0FFoATT0xJw1prwejRMGQITJxYdGRmZh2DkwUzM2v3Vl0VHngArrwyzcew4YZw4IHwzjtFR2Zm1r61iWRBUldJz0i6Pft5CUn3Sno1+7dnbttjJU2WNEnStrny9SQ9n607R5LqOFZV9c3MrG2S0tiFSZPSpG5XXgk/+AH8+c/w1VdFR2dm1j61iWQB+G/g5dzPxwD3RcQqwH3Zz0haHRgJrAEMB86X1DWrcwEwGlglW4aXHqTa+mZm1vb16AFnnJHmZthoI/j1r2HdddOAaDMza5zCkwVJfYEfAxflikcAl2WvLwN2zpVfGxFzImIKMBkYLGlZYNGIeDQiArg8Vyev2vpmZtZODBwId94Jt9wCn30GW20Fe+wBb75ZdGRmZu1H4ckC8GfgaCD//IqlI2IGQPbvUll5H2BqbrtpWVmf7HVpealq65uZWTsiwYgRacbnk06C229PScRvfgMffVR0dGZmbV+3Ig8uaQdgZkQ8JWloJVXKlEU95c1dH0mjSd2V6N27NzU1NWUDtZY3e/Zst39B3PbFcvs3zeabw6qrdueiiwZwxhnL8H//9yWjRr3JjjtOp1u3sh/5Zbn9i+O2L5bbv3MqNFkANgF2krQ9sACwqKQrgXckLRsRM7IuQrVT7UwD+uXq9wWmZ+V9y5SXqrY+ETEWGAswcODAGDp0aAVv01pCTU0Nbv9iuO2L5favzp57wtNPwxFHzM8556zCXXetwhlnwE47pTsRDXH7F8dtXyy3f+dUaDekiDg2IvpGRH/SwOP7I2JfYAIwKttsFHBr9noCMFJSd0kDSAORJ2ZdlT6RNCR7itH+uTp51dY3M7MO4Ic/TAOeJ0xICcLOO6dJ3p58sujIzMzalrYwZqGc04BtJL0KbJP9TES8CFwPvATcBfwiIuZmdQ4hDZKeDLwG3AkgaSdJJze1vpmZdUwS7LhjemrS+efDSy/BBhukx69Onlx0dGZmbUPR3ZC+FRE1QE32+j1gqzq2GwOMKVP+JDCoTPkE0h2FJtU3M7OObb754JBD4Cc/gdNPT/MyXH89HHQQHH889PHjLsysE2urdxbMzMxa1WKLwamnwmuvwcEHw8UXw8orw5FHwqxZRUdnZlYMJwtmZmY5yy4L556bZoLeay84+2xYcUU48UT4+OOiozMza11OFszMzMoYMAAuvTSNaRg2LM3TsOKKcM01/Zg9u+jozMxah5MFMzOzeqy+OowfD088AeuvD2PHrkT//vDHP/pOg5l1fE4WzMzMKrD++nDXXXDuuU8zeDAcdxz07w9/+AN8+GHR0ZmZtQwnC2ZmZo2wxhofc8cdMHEibLop/P73KWk44QR4//2iozMza15OFszMzJpggw3SpG5PPw1bbgknn5yShuOOg5kzi47OzKx5OFkwMzOrwrrrwk03wXPPwfDhcNppsMIKae4GT+5mZu2dkwUzM7NmsNZaaTK3l1+G/fZL8zT84Aewxx5pcLSZWXvkZMHMzKwZDRwIY8fCG2/AMcfAvffC4MGwxRZw550QUXSEZmaVc7JgZmbWApZdNs0IPXUqnHVW6pK0/faw9tpwySXwxRdFR2hm1jAnC2ZmZi2oRw84/HB47TW47LJ0Z+HAA6FfP/jd7+A//yk6QjOzujlZMDMzawXzzw/77w//+hfcdx9sskm689C/P+y9Nzz2mLsomVnb42TBzMysFUnpUau33JK6Jh12WBrLsNFGsOGGcOWV8OWXRUdpZpY4WTAzMyvIiium8QzTpsF558HHH6cnKfXtC7/5Teq6ZGZWJCcLZmZmBVtkEfj5z+Gll+Cuu1IXpbPOgpVXhm22gfHj4auvio7SzDojJwtmZmZtRJcusO22cPPN8NZbaVbof/87zdXQr1+aHXrKlKKjNLPOxMmCmZlZG7TccnD88fD66/C3v6XxDKefDiutlBKK667z41fNrOU5WTAzM2vDunZN8zPceiu8+SaccEKaJXrkyDSXwyGH+ElKZtZynCyYmZm1E337pmRhypQ0M/QOO6S5GzbaCFZbDf74xzRY2sysuThZMDMza2e6doWtt4YrroC334Zx42CppdKYhuWXh2HD4KqrYPbsoiM1s/bOyYKZmVk7tuiiaUbohx5K8zYcfzy8+irsu29KIPbaK83p4PENZtYUhSYLkvpJekDSy5JelPTfWfkSku6V9Gr2b89cnWMlTZY0SdK2ufL1JD2frTtHkuo4ZlX1zczM2qqVVoKTTkrzMzz0EPz0p/DAA7DLLrD00unnu++Gr78uOlIzay+KvrPwNXBERKwGDAF+IWl14BjgvohYBbgv+5ls3UhgDWA4cL6krtm+LgBGA6tky/DSg1Vb38zMrD3o0gU22yxN9DZ9ekoQdt0VbroJhg9PA6N//nOoqYG5c4uO1szaskKThYiYERFPZ68/AV4G+gAjgMuyzS4Dds5ejwCujYg5ETEFmAwMlrQssGhEPBoRAVyeq5NXbX0zM7N2pVu3NIbhkkvgnXfSHA5bbQWXXgpbbJESh//3/+COO2DOnKKjNbO2plvRAdSS1B9YF3gcWDoiZkBKKCQtlW3WB3gsV21aVvZV9rq0vFS19ZE0mnQHgt69e1NTU9Pwm7MWMXv2bLd/Qdz2xXL7F6u9t//ii8PPfgajRnVl4sQlePjhJbn22l6MG9eNhRb6miFD3mOzzWax4Ybvs+CCbeu2Q3tv+/bO7d85tYlkQdIiwI3AryLi43qGC5RbEfWUN3d9ImIsMBZg4MCBMXTo0LKBWsurqanB7V8Mt32x3P7F6kjtv9126d85c+D+++Hmm7txyy1Lc//9S9O9e7ojMWJEmudh2WWLjRU6Vtu3R27/zqnoMQtImo+UKFwVETdlxe9kXYPI/p2ZlU8D+uWq9wWmZ+V9y5SXqra+mZlZh9O9e0ocxo6FGTPgwQfT3YfnnktdlJZbDtZfP83xMHEifPNN0RGbWWsp+mlIAsYBL0fEn3KrJgCjstejgFtz5SMldZc0gDQQeWLWZekTSUOyfe6fq5NXbX0zM7MOrWtX2Hxz+POf4Y03UsJw6qkpoTjlFNhww3SX4YAD4IYb4KOPCg7YzFpU0d2QNgH2A56X9GxWdhxwGnC9pIOAt4A9ACLiRUnXAy+RnqT0i4io7VB5CHApsCBwZ7YgaSdg/Yj4fVPqm5mZdVYSrLVWWo49Ft57D+66C/72N5gwIc0e3a0bbLIJbLNNWtZbLyUcZtYxFJosRMQjlB8vALBVHXXGAGPKlD8JDCpTPoF0R6FJ9c3MzCzp1Qv22SctX38Njz0Gt9+eHs36u9+lZfHFYcstv0seVlqp6KjNrBpF31kwMzOzdqhbN9h007ScdhrMnAn33Qd//zvce2+a0wFgwICUNGy9NfzoR2lWaTNrP5wsmJmZWdWWWgr23jstEfDvf6ek4e9/h2uvTYOnAVZbLSUNtUtbeMqSmdXNyYKZmZk1KwkGDkzLoYemLktPPJGesvTgg3DVVfDXv6ZtV1553uRh+eWLjd3M5uVkwczMzFpUt26w0UZpOeaYlDw8++x3ycONN8K4cWnbFVaAjTf+bllrrVTfzIrhPz8zMzNrVd26pXkb1l8fjjgC5s6FF15IicPDD6d/r7kmbbvQQjB4cEocFllkCQYPTmVm1jqcLJiZmVmhunaFtddOy2GHpTEPU6fCP/+ZlkcfhdNPh7lz12LffZ0smLUmJwtmZmbWpkhp7MLyy8PIkans00/h4oufoV+/dYsNzqyTKXQGZzMzM7NKLLwwrLmmp4s2a21OFszMzMzMrCwnC2ZmZmZmVpaTBTMzMzMzK8vJgpmZmZmZleVkwczMzMzMynKyYGZmZmZmZTlZMDMzMzOzspwsmJmZmZlZWU4WzMzMzMysLCcLZmZmZmZWliKi6BjaLUmfAJNa+bCLAdXOd9/YfVSyfUPb1LW+XHmlZUsCsxqIq7lV2/4t0fYNbdfYdZWUtce2b8o+Wvvcr6vc7d+0bVri3IfWb/+22vYNbefPnqbvw589dcfQGvtoq589AyOiRwNxNb+I8NLEBXiygGOObe19VLJ9Q9vUtb5ceSPK2l37t0TbN7RdY9dVUtYe276l2r85z323f/O2f0uc+0W0f1tt+yLa3+d+5dv4s6e49u8onz21i7shtT+3FbCPSrZvaJu61pcrr7SsCNXG0RJt39B2jV3XVtu/M5z7dZW7/Zu2jc/9pu/Dnz3faavt78+eltuHP3ty3A2pCpKejIj1i46js3L7F8dtXyy3f7Hc/sVx2xfL7V+sotrfdxaqM7boADo5t39x3PbFcvsXy+1fHLd9sdz+xSqk/X1nwczMzMzMyvKdBTMzMzMzK8vJQh0kDZc0SdJkSceUWS9J52Tr/yXph5XWtfpV2fYXS5op6YXWjbrjqKD9V5X0qKQ5ko5sTF2rX0Pnrz93mle59pa0hKR7Jb2a/duzjrpl27vS+p1VY9tc0rFZG0+StG0d+6yqfkfWXO0taT1Jz2frzpGkOo5XVf2OoKXbXFJ3Sddl5Y9L6l9HHFXVn0cRj2Bq6wvQFXgNWBGYH3gOWL1km+2BOwEBQ4DHK63rpWXaPlu3OfBD4IWi30t7XCps/6WADYAxwJGNqeulwfav9/z1507LtzdwBnBM9voY4PQy9eps70rqd+alMW0OrJ61bXdgQNbmXcvss6r6HXlprvYGJgIbZZ89dwLblTlWVfU7ytLSbQ78HPhr9nokcF0dcVRVP7/4zkJ5g4HJEfF6RHwJXAuMKNlmBHB5JI8Bi0tatsK6Vrdq2p6IeAh4v1Uj7lgabP+ImBkRTwBfNbau1a+C89efO82ojvYeAVyWvb4M2LlM1frau5L6nVYj23wEcG1EzImIKcBkUtuXqrZ+h9Uc7Z19xiwaEY9G+oZ5OeXP62rrdwit0Ob5fY0Htiq9U1Nt/VJOFsrrA0zN/TwtK6tkm0rqWt2qaXurXjVt699Ly/PnTstbOiJmAGT/LlVmm/rau5L6Nq+62qzS87ra+p1NY9urT/a6tLxUtfU7suZs82/rRMTXpFmee5Ucr9r683CyUF65DKv0sVF1bVNJXatbNW1v1aumbf17aXn+3Gkb3N6to9p29u+pcar9fPHnU+M1pc2q/Z7U6N+Hk4XypgH9cj/3BaZXuE0lda1u1bS9Va+atvXvpeX5c6flvVPbrTH7d2aZbepr70rq27zqarNKz+tq63c2jW2vadnr0vJS1dbvyJqzzb+tI6kbsBjf7/ZUbf15OFko7wlgFUkDJM1PGgAyoWSbCcD+SoYAH2W3liqpa3Wrpu2tetWcvz73W54/d1reBGBU9noUcGuZbepr70rq27zqarMJwMjs6S0DgFVIgzabu35n06j2yj5jPpE0JOvbvj/lz+tq63dkzdnm+X3tDtyfjUv4VrX1v6elRoO394X01JF/k0am/zYr+xnws+y1gPOy9c8D69dX10urtf01wAzS4NtpwEFFv5/2tlTQ/stkbfsx8GH2etG66nppVNt/7/z1506rt3cv4D7g1ezfJbJtlwPuaKi966rvpfFtnm3/26yNJ5F7gg5wUe3535T6nWVpxvZeH3ghW3cu303quxNwclPrd8SlFdp8AeAG0mDoicCKuTrPVlO/rsUzOJuZmZmZWVnuhmRmZmZmZmU5WTAzMzMzs7KcLJiZmZmZWVlOFszMzMzMrCwnC2ZmZmZmVpaTBTMzMzMzK8vJgpmZmZmZleVkwcysAJIulXRK7ucXJQ0tLqKGSXpD0tZt+Tj11S1t8zLrQ9KnksY05dhFknS/pC8kPVJ0LGbWsThZMDMro7W+GNeKiDUioqa1jteQ1n7/bcjaEfHbhjaSdKykO0rKXq2jbGTu5+UkTWtMQJK6Sxon6U1Jn0h6RtJ2+W0iYkvSbN9mZs3KyYKZWSNI6lZJWWfWSdrjIWATSV0BJC0DzAf8sKRs5WzbWtsDdzXyWN2AqcCPgMWA44HrJfWv5g2YmVXCyYKZWQOyq+y/kfQv4FNJ3eooO0bSa9nV35ck7ZLbx7qSns7WXQcsUOYYW2ev69xPbtsjJf1L0keSrpO0QLaun6SbJL0r6T1J5+bqLSfpxmzdFEmH1fF+rwCWB26TNFvS0bnV69Rx3HLtUe/xsu3/k73PSZK2quA4q0mqkfRh1nVrp3p+b/W2eUMk/VbSBbmfe0r6KovlCVJysE62enPgAWBSSdlrETE9t9vtgTtybXZU9j4/ze4eLC3pzizmv0vqGRGfRsSJEfFGRHwTEbcDU4D1GvN+zMyawsmCmVll9gZ+DCweEV/XUfYasBnp6u9JwJWSlpU0P3ALcAWwBHADsFs9xyq7n5Jt9gSGAwOAtYADsivatwNvAv2BPsC1AJK6ALcBz2XlWwG/krRt6cEjYj/gLWDHiFgkIs6o77jl2gj4pr7jSRoIHApsEBE9gG2BNxp4f/Nl+7wHWAr4JXBVtq95NKHNy1kTeDb38zrApIj4IiK+BB4nJQRk/z4MPFJS9u1dhSz+zYF7c/vcDdgG+AGwI3AncBywJOn/6O8ldJKWzrZ/sZHvx8ys0ZwsmJlV5pyImBoRn9dVFhE3RMT07OrvdcCrwGBgCOkq9J8j4quIGE+6Ml1WPfspjWd6RLxP+gK9TrbNcsBR2dXoLyKidsDrBkDviDg5Ir6MiNeBC4GRNE6545Zrj4aONxfoDqwuab7sqvlrDRxnCLAIcFq2z/tJydHeZeJsVJvXoVyy8Fzu5wf5LjHYjJQsPFxS9mBu+82B5yLik1zZ/0bEOxHxn6zu4xHxTETMAW4G1s0HlCUcVwGXRcQrjXw/ZmaN5mTBzKwyUxsqk7S/pGezLjIfAoNIV4iXA/4TEZHb/M26DlTPfvLezr3+jPQluh/wZu7OR94KwHK1+8z2exywdF1x1KHccWvl26Pe40XEZOBXwInATEnXSlqugeMsB0yNiG9y694k3bko1ag2L5XdmVgJeD5XvDbzJg8PAZtK6klKjF4F/glsnJUN4vvjFeYZAA28k3v9eZmfv23f7O7QFcCXpLsyZmYtzsmCmVllor4ySSuQrpwfCvSKiMWBFwABM4A+kpSru3y5gzSwn4ZMBZZX+QHGU4EpEbF4bukREdvXsa9y77ch+ToNHi8iro6ITUmJRQCnN7D/6UC/7EtzreWB/5TZtuI2r8PqpGTjM4BsP0OZ987Co6SuYqOBfwBExMdZnKOB6RExJbf99sDfGhHDt7LjjyMlW7tFxFdN2Y+ZWWM5WTAzax4Lk77wvgsg6aekK8uQvlR+DRyWDfzdle93K6pkPw2ZSPqSfJqkhSUtIGmT3LqPs0HFC0rqKmmQpA3q2Nc7wIoVHreuWOo8nqSBkraU1B34gnQVfW4D+3wc+BQ4WtJ8SvNS7Eg2LqNEY9q8nDWBpSStJGlB4A+kpOaN2g2y7lZPAoeTuhDVeiQry49XGAB0r6Lr0AXAaqRxJJ83tLGZWXNxsmBm1gwi4iXgLNKX1HdIXzZrrzZ/CexKGgz8AbAXcFNj91NBDHNJX55XJg1QnpYdK79uHdKTdGYBF5GujJfzR+B3WReiIys5fh2x1HW87sBpWfnbpAHLxzWwzy+BnYDtsnrnA/uX+wLemDavw5rA3aQBx5NJv4vXgdI5GB7MYs9PhvZwVpbvgvRjvt8FqSLZ3aaDSW35ttITqmZL2qcp+zMzawzN253TzMysc5L0BTAHOAdYH7goIm5spn3fAZwbEU1KGCrY/72kQd0TI2KrhrY3M6tUZ5g4x8zMrEER8e08DEqzLL/cjLuvIc3D0CIiYpuW2reZdW6+s2BmZpaTPcnoHWBhDyQ2s87OyYKZmZmZmZXlAc5mZmZmZlaWkwUzMzMzMyvLyYKZmZmZmZXlZMHMzMzMzMpysmBmZmZmZmU5WTAzMzMzs7KcLJiZmZmZWVlOFszMzMzMrKz/D9a6RXfCZnueAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import pyradi.ryutils as ryutils\n", "import pyradi.ryplot as ryplot\n", "import numpy as np\n", "from scipy.interpolate import interp1d\n", "%matplotlib inline\n", "\n", "rangeTab = np.linspace(0, 100000, 1000)\n", "atmogamma = 0.00015\n", "tauTab = np.exp(- atmogamma * rangeTab)\n", "Intensity=200\n", "Irradiancetab=np.linspace(0.05e-6, 100e-6, 400)\n", "r = np.zeros((Irradiancetab.shape))\n", "for i,Irradiance in enumerate(Irradiancetab):\n", " r[i] = ryutils.rangeEquation(Intensity = Intensity, Irradiance = Irradiance, rangeTab = rangeTab,\n", " tauTab = tauTab, rangeGuess = 1, n = 2)\n", "\n", "p=ryplot.Plotter(1, 2, 1, figsize=(12,12))\n", "p.semilogX(1,Irradiancetab*1e6,r,'Range equation for 200 W/sr target through {} km-1 atmosphere'.format(atmogamma*1000.),\n", " 'Irradiance threshold [$\\mu$W/m2]','Detection range [m]');\n", " " ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "## Absolute humidity" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "The [`ryutils.abshumidity`](http://nelisw.github.io/pyradi-docs/_build/html/ryutils.html#pyradi.ryutils.abshumidity) function calculates atmopsheric absolute humidity [g/m3] for temperature in [K] between 248 K and 342 K. This function provides two similar equations, but with different constants.\n", "\n", "Atmospheric humidity is commonly expressed in relative humidity (RH). Relative humidity is the quantity of water vapor in the atmosphere, expressed as a percentage of the maximum absolute humidity. The maximum absolute humidity is determined by the water-vapor partial pressure at the atmospheric temperature and is given by [this equation](http://www.vaisala.com/Vaisala%20Documents/Application%20notes/Humidity_Conversion_Formulas_B210973EN-D.pdf) (Equation 1 in the function):\n", "$$\n", "q=\\frac{1325.252}{T} 10^{7.5892(T - 273.15)/(T -32.44)},\n", "$$\n", "in units of [g/m$^3$] and temperature in [K]. This equation has an error of less than 1% over the range $-$20 $^\\circ$C to 0 $^\\circ$C, and less than 0.1% for the temperature range 0 $^\\circ$C to 50 $^\\circ$C. \n", "\n", "The [alternative equation](#http://www.see.ed.ac.uk/~shs/Climate%20change/Data%20sources/Humidity%20with%20altidude.pdf) (Equation 1 in the function) is: \n", "$$\n", "q=\\frac{1324.37872}{T} \\exp[\n", "17.67*(T - 273.16)/(T - 29.66))\n", "$$\n", "\n", "The function signature is:\n", " `abshumidity(T, equationSelect = 1)` \n", " \n", "- `temperature (np.array[N,] or [N,1])` in [K].\n", "- `equationSelect (int)` select the equation to be used.\n", "\n", "\n", "In this example the values from the equations are compared to [data from](http://rolfb.ch/tools/thtable.php?tmin=-25&tmax=50&tstep=5&hmin=10&hmax=100&hstep=10&acc=2&calculate=calculate):" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "run_control": { "frozen": false, "read_only": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " deg C Testvalue Fn value % Error\n", "[[ 5.00000000e+01 8.27800000e+01 8.28288052e+01 -5.89229459e-02]\n", " [ 4.50000000e+01 6.52500000e+01 6.53105744e+01 -9.27481918e-02]\n", " [ 4.00000000e+01 5.09800000e+01 5.10475988e+01 -1.32422991e-01]\n", " [ 3.50000000e+01 3.94700000e+01 3.95318519e+01 -1.56460806e-01]\n", " [ 3.00000000e+01 3.02600000e+01 3.03161305e+01 -1.85150640e-01]\n", " [ 2.50000000e+01 2.29700000e+01 2.30097997e+01 -1.72968594e-01]\n", " [ 2.00000000e+01 1.72400000e+01 1.72744162e+01 -1.99232035e-01]\n", " [ 1.50000000e+01 1.28000000e+01 1.28192871e+01 -1.50453884e-01]\n", " [ 1.00000000e+01 9.38000000e+00 9.39702200e+00 -1.81142459e-01]\n", " [ 5.00000000e+00 6.79000000e+00 6.79913283e+00 -1.34323439e-01]\n", " [ 0.00000000e+00 4.85000000e+00 4.85173751e+00 -3.58120391e-02]\n", " [-5.00000000e+00 3.41000000e+00 3.41141647e+00 -4.15215865e-02]\n", " [-1.00000000e+01 2.36000000e+00 2.36126844e+00 -5.37186742e-02]\n", " [-1.50000000e+01 1.61000000e+00 1.60720482e+00 1.73915730e-01]\n", " [-2.00000000e+01 1.08000000e+00 1.07451558e+00 5.10408511e-01]\n", " [-2.50000000e+01 7.10000000e-01 7.04731527e-01 7.47585870e-01]]\n", "Highest recorded absolute humidity was 37.51799299496529, dew point 34 deg C\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvUAAAGKCAYAAACM3r4EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABiKklEQVR4nO3dd5gUVdbH8e9hYEhDzjkIjCIKrApmEFYFDBhAMSuGdc367q6u65p1dYNZ12VFBbMirqiomGbNAQQjsiIKQ0aJQ2Y47x9VAz1Nz9DTMz3dPfP7PE89XXW7btWp7iuevnPrlrk7IiIiIiKSuWqkOgARERERESkfJfUiIiIiIhlOSb2IiIiISIZTUi8iIiIikuGU1IuIiIiIZDgl9SIiIiIiGU5JvYhUCWb2qJndXMHHPNPM3q/IY6bJuU4xsymlvJ9nZufEs6+IiKQHJfUiklHChHOFmdVOdSyRKjMpLy93f8LdD0tkXzNzM+uWvOh2zsx+MrNfpzKGimJm15vZ46mOQ0Qyn5J6EckYZtYZOAhw4OjURiNSOjOrWRXOISKZQUm9iGSS04GPgUeBM2K839zM3jCzNWb2XzPrBGCBO81sqZmtMrMvzaxX+F4jMxtvZsvMbK6ZXWNmO/zbaGadw17qmhFleWZ2jpntBjwI7GdmBWa2Mny/tpn93czmmdkSM3vQzOqWdoHh/ivM7EczGxpRXqx3OrKHNyK2s8wsP6x/vpntE17rSjO7L6Jusb8qmNmhZvZd+NncB1isfc3s3bD4i/A6TzSzr83sqIj9a5nZz2bWJ8a1zTSzIyO2a4b7/srM6pjZ42b2SxjvZ2bWKsYxHgM6Ai+FMfwhLN/XzD4M635hZgOjvqebw/cLzOwlM2tmZk+Y2erwXJ0j9nczu8TM5oTx/S2yTZjZ6PBaVpjZ60XtLKLuhWb2PfB9WHZ3+L2sNrNpZnZQWD4EuBo4MYzrizJ812eb2Tzg7Z3FJCLVg5J6EckkpwNPhMvhMZK+U4CbgObAjHA/gMOAg4EeQGPgROCX8L17gUZAV2BAeI6zyhKUu88Ezgc+cvccd28cvnV7eM4+QDegHXBtKYfqD8wK4/8rMNbMrJT9Y9XvTnB9dwF/An4N7A6cYGYDoiuYWXPgeeCa8Lw/AAfEOri7Hxyu9g6v8xlgPHBqxG7DgEXuPiPGIZ4CTorYPhz42d0/J/iR1gjoADQj+DzXx4jhNGAecFQYw1/NrB3wCnAz0BT4HfC8mbWIqDoKOI3gO9gF+Ah4JNx/JnBd1KmOBfYGfgUMB0YDmNkxBIn4cUAL4L3wuiIdQ/Bd9Ay3PyNoA02BJ4HnzKyOu78G3Ao8E15L7xifWUkGALsR/HcQT0wiUsUpqReRjGBmBwKdgGfdfRpB8nly1G6vuPu77r6RIKHdz8w6AJuBBsCugLn7THdfZGZZBAnwH919jbv/BPyDIPkrb7wGnAtc7u7L3X0NQQI3qpRqc9393+5eCIwD2gA79FaX4iZ33+DuU4C1wFPuvtTdFxAken1j1BkGfOvuE9x9M8GPgcVlOOfjwDAzaxhunwY8VsK+TwJHm1m9cPvksAyC76gZ0M3dC919mruvjjOGU4HJ7j7Z3be6+xvA1PDaijzi7j+4+yrgVeAHd3/T3bcAz7HjZ3N7+L3NI/hMin6M/Ab4S9iGthB8p32iesb/EtZdD+Duj7v7L+6+xd3/AdQGcuO8tpJc7+5rw3PEE5OIVHFK6kUkU5wBTHH3n8PtJ9lxCE5+0Yq7FwDLgbbu/jZwH3A/sMTMxoRJaHMgG5gbcYy5BL255dUCqAdMC4eErAReC8tLsi2Zdvd14WpOGc65JGJ9fYztWMdqS/HPzSO3d8bdFwIfAMebWWNgKNv/QhK972yCXvGjwsT+aLYn9Y8BrwNPm9lCM/urmdWKM4xOwMiizzn8rA8k+FFUpKyfTeRnMJfgcyo6190R51lOMFypXQl1MbP/C4fGrArrNCJoe+UReY54YhKRKk432IhI2rNgHPoJQJaZFSW+tYHGZtbb3b8IyzpE1MkhGO6wEMDd7wHuMbOWwLPA74HrCXqIOwHfhlU7AgtihLE2fK0HFPUgt45436P2/5kgWdw97Ckvr7XhuYu0LmnHMlpE8c/NIrfjNA44h+D/KR/t5HqLhuDUIPgLwWyA8K8ENwA3hOPbJxMMRRob4xjRn3U+8Ji7n1vGuEvTAfgmXO9I2I7Cc93i7jF/uETHF46fvxIYDHzj7lvNbAXb71uIvhaI77uOrBdPTCJSxamnXkQywTFAIcEY5T7hshvBkJLTI/YbZmYHmlk2wdj6T9w934IbRvuHPb9rgQ1AYTjM5VngFjNrEA5XuIJgSEkx7r6MINk/1cyyzGw0wdjsIkuA9uG5cfetwL+BO8MfEphZOzM7PMHPYAYwyoIbUfcGRiR4nGivALub2XEW3AR8CaX/YFhCcP9BpP8QjD2/lGCMfWmeJrjH4bds76XHzA4xsz3CIVGrCX5sFcYZw+MEvf+Hh99NHTMbaGbtdxJLaX5vZk3C4VuXAs+E5Q8CfzSz3cO4G5nZyFKO0wDYAiwDaprZtUDDiPeXAJ2t+M3ZMyjbd13WmESkClJSLyKZ4AyCMdHz3H1x0UIwpOYU2z4jzZMENzwuB/YiuHEWgiTq38AKgqEUvwB/D9+7mCDRnwO8Hx7j4RLiOJegh/8XgptPP4x4722Cnt3FZlY0ROhKYDbwsZmtBt4k8bHUfyb4EbGCoEf7ydJ3j084nGkkcBvBdXUnGE5TkuuBceFQjxPCY6wnuNm2CzBxJ+dbRHCT6v5sT5Qh+CExgSChnwn8lxg/rkJ/Aa4JY/idu+cT3Mx6NUHynE/wPZXn/3EvAtMIEuxXCP9i4O4vENwA/XT4nX5NMOSoJK8TjOH/H0Hb20DxoTPPha+/mNnn4XqZvusEYhKRKsiC4ZMiIiKJC3uge7j7qTvdOc2ZmQPdi4YGiYhkAo2pFxGRcjGzpsDZVMCsQSIikpi4knrb/sCRndkQ76PHRUQk85nZuQRTPj7m7vH+v0JERCpYXMNvzGw9wYNASt0NuNvdG1VEYCIiIiIiEp94h9986O7jdraTmUU/CEZERERERJJMN8qKiIiIiGQ43SgboXHjxt6tW7dUhyFpZO3atdSvXz/VYUiaUbuQWNQuJBa1C4ll2rRpP7t7aU8YL7NSk3oza0Uwo8E3wOcEcwCvBf7m7ktKq5uJWrVqxdSpU1MdhqSRvLw8Bg4cmOowJM2oXUgsahcSi9qFxGJmcyv6mDt7MMcTBE9Q7ABMAd4KX0t6MIuIiIiIiFSynSX1Nd19nLvfB6xy9wnuPgWoUwmxiYiIiIhIHHaW1NeKWP9NxLrG4ouIiIhI3FatSnUEVdvOkvpjzcwA3P0LADPLBq5MdmAiIiIikvk2boTf/Ab694c1a1IdTdVValLv7ks9as5Ld9/k7h8nNywRERERyXTz58PBB8OYMXDccVCvXqojqrrKPIzGzBoBlwB9gZzI99z9sASOdzlwDuDAV8BZQD3gGaAz8BNwgruviFF3CHA3kAU85O63heVN46kvIiIiIsnxzjtw4omwYQNMnAjHHpvqiKq2nQ2/ieU5YCDwNkHiHLmUiZm1I/iBsLe79yJIzkcBVwFvuXt3ghl3ropRNwu4HxgK9AROMrOe4ds7rS8iIiIiFc8d/vEPOPRQaNYMPv1UCX1lSOSG132BZu6+uQJjqGtmmwl66BcCfyT44QAwDshjx3H8/YDZ7j4HwMyeBoYD34avO6svIiIiIhWooADOPhuefRaOPx4eeQQaNEh1VNVDIkn9+8BuwJflPbm7LzCzvwPzgPXAFHefYmat3H1RuM8iM2sZo3o7ID9iez7QP1yPpz4AZnYecB5AixYtyMvLK+9lSRVSUFCgNiE7ULuQWNQuJJbq1C7y8+ty7bW9mDevHuedN4dRo/KZNi3VUVUfiST1ZwKTzewToNhTZd39xrIcyMyaEPSqdwFWAs+Z2anxVo9R5jHKSuXuY4AxALm5ua6nvkkkPQlQYlG7kFjULiSW6tIuJk2Ciy6C7GyYMgUGD94F2CXVYVUriYypv4XgCbOtgO4RS7cEjvVr4Ed3XxYO55kI7A8sMbM2AOHr0hh154dxFGlPMHSHOOuLiIiISDkUFsKf/wzDh0OPHjBtGgwenOqoqqdEeupHAT2KhreU0zxgXzOrRzD8ZjAwFVgLnAHcFr6+GKPuZ0B3M+sCLAjjOjl8b1Ic9UVEREQkQcuXw8knw+uvB+Po77sP6tRJdVTVVyJJ/RygQm6SdfdPzGwC8DmwBZhOMBQmB3jWzM4mSPxHAphZW4KpK4e5+xYzuwh4nWDWnIfd/Zvw0LfFqi8iIiIi5TdjRjDv/IIFwRz0556b6ogkkaT+MWCSmd3LjmPq3y7rwdz9OuC6qOKNBL320fsuBIZFbE8GJsfY75dY9UVERESkfMaPD54Q26wZvPtu8KRYSb1EkvoLw9dbo8od6Fq+cEREREQkHW3aBFdcAfffDwMHwjPPQMsS5xeUylbmpN7duyQjEBERERFJTwsXwsiR8OGH8H//B7fdBjUT6RqWpIn76zCzfOBVguEub7j72qRFJSIiIiJp4b334IQTYM2aoHf+hBNSHZHEUpYpLfsBnwCnAT+Z2RtmdrmZ9UhOaCIiIiKSKu5wzz0waFDwVNhPPlFCn87i7qkPp7AcC4w1s5rAwQQ3rf7HzLIJevAnA++4+8ZkBCsiIiIiybduXTCjzZNPwtFHBzfHNmqU6qikNIk8fAp33+Lub7v779y9J8FDpGYBF4eLiIiIiGSgH36A/faDp56Cm2+GF15QQp8JynyLg5ndWMJbG4FngdfKFZGIiIiIpMTkyXDKKWAWrA8ZkuqIJF6J9NT3AK4EDgG6ha9XAn2B3wJzzExNQERERCRDbN0KN9wARx4JnTvD1KlK6DNNIpMR1QBGufsLRQVmNhw42d33NbMzCJ7oqh57ERERkTS3ciWcdhq8/HLw+uCDUK9eqqOSskqkp/5wYFJU2cvA0HD9cWCX8gQlIiIiIsn31Vew997w2mtw330wbpwS+kyVSFL/A8Ewm0jnh+UAzQHNYS8iIiKSxp56CvbdN5jpJi8PLrwwGEsvmSmR4TfnABPN7EpgAdAOKASOC9/PBf5cMeGJiIiISEXavBn+8Ae46y448EB49llo0ybVUUl5lTmpd/fPzaw7sC/QFlgEfOTum8P33wXerdAoRURERKTcFi+GE0+Ed9+FSy6Bv/8datVKdVRSERLpqSdM4N+r4FhEREREJEk++ghGjIAVK+Cxx+DUU1MdkVSkhB4+JSIiIiKZwR3++U8YMADq1AmSeyX0VY+SehEREZEqatUqOPNMuOAC+PWvg/nne/dOdVSSDErqRURERKqgF1+Enj3h8cfh2muDeeibNEl1VJIsZU7qzewOM+uThFhEREREpJwWL4aRI+GYY6B5c/j44+BpsTXUlVulJfL11gJeN7OvzexKM2tf0UGJiIiISNm4w9ixsNtu8NJLcOutwXCbffZJdWRSGcqc1Lv7xQRTWV4F9AFmmtmbZna6meVUcHwiIiIishOzZ8PgwXDOObDnnvDll/DHP2q6yuokoT/EuHuhu7/s7icRzFffAngUWGxmD5lZuwqMUURERERi2LIFbr8d9tgDPv8cxoyBd96BHj1SHZlUtoSSejNraGZnm9k7BA+a+gQ4CNgNKABeLePxGpvZBDP7zsxmmtl+ZtbUzN4ws+/D15i3dpjZEDObZWazzeyqiPK46ouIiIhkos8/h3794KqrYOhQ+PZbOPdcjZ2vrhK5UXYCsAA4DngQaOvu57n7B+6eD1wBdCnjYe8GXnP3XYHewEyC4T1vuXt34K1wOzqWLOB+YCjQEzjJzHqGb++0voiIiEimWbcO/vCHIKFfvBiefx4mToS2bVMdmaRSIr/lPga6u/sR7v6Mu2+MfNPdtwKt4j2YmTUEDgbGhvU3uftKYDgwLtxtHHBMjOr9gNnuPsfdNwFPh/WIs76IiIhIxnjrrWCozd/+BqNHB73zxx2X6qgkHdRMpJK7L44uM7Mr3P2O8P11ZThcV2AZ8IiZ9QamAZcCrdx9UXi8RWbWMkbddkB+xPZ8oH+4Hk99zOw84DyAFi1akJeXV4bQpaorKChQm5AdqF1ILGoXEktFtYvVq2vyz3/uwmuvtaF9+3Xceef/6NNnJTNmlPvQUkUkktRfC/w9Rvk1wB0JxvAr4GJ3/8TM7ib+oTIWo8zLcnJ3HwOMAcjNzfWBAweWpbpUcXl5eahNSDS1C4lF7UJiKW+7cIfnnoOLL4ZffglmtPnzn+tRt26fCotRqoa4k3ozGxSuZpnZIRRPqLsCaxKMYT4w390/CbcnECT1S8ysTdjL3gZYWkLdDhHb7YGF4Xo89UVERETS0vz5cMEFwZzze+8NU6ZA796pjkrSVVl66seGr3WAhyPKHVgCXJxIAO6+2MzyzSzX3WcBg4Fvw+UM4Lbw9cUY1T8DuptZF4Kbd0cBJ4fvTYqjvoiIiEha2boVHnwwmNVmyxb4xz/gkkugZkKDpqW6iLt5uHsXADMb7+6nV3AcFwNPmFk2MAc4i+Am3mfN7GxgHjAyPH9b4CF3H+buW8zsIuB1IAt42N2/CY95W6z6IiIiIulq5sxgWsoPPoBDDw2S+65dUx2VZIK4knozO9jd3w03H40YilOMu7+dSBDuPgPYO8Zbg2PsuxAYFrE9GZgcY79fYtUXERERSTebNgUPkbr5ZsjJgXHj4LTTwGLdPSgSQ7w99Q8AvcL1sSXs4wRj60VEREQkTh9/DOecA998AyedBHfdBS1jztknUrK4knp37xWxXtYHS4mIiIhIlDVr4Jpr4N57oX17ePllOOKIVEclmUq3XIiIiIhUsldfhfPPh/x8uPBCuPVWaNAg1VFJJot3TP2N8ezn7teWLxwRERGRqmvZMrjsMnjySejZM7ghdr/9Uh2VVAXx9tRHzgVfBzieYDrJuUBHoB/wfMWGJiIiIlI1uMPjj8Pll8Pq1XD99cGUlbVrpzoyqSriHVN/VtG6mT0NnOTuz0eUHYemjBQRERHZwU8/wW9+Ezw8ar/94KGHgl56kYpUI4E6Q4H/RJW9SMQ0kyIiIiLVXWEhTJjQnt13hw8/hPvug/ffV0IvyZFIUj8buDCq7ALgh/KHIyIiIpLZtm6F55+HPfeE++/vxqBB8O23wQ2xNRLJvETikEjTOge4wszmm9knZrYA+L+wXERERKRacodXXoG994YRI4Lk/vrrv2bSJOjQYef1RcqjzEm9u08HugMnAXeEr93d/fMKjk1EREQkI7z9NhxwABx5JKxaBePHw9dfw4ABP+upsFIp4krqzezgiPVBwEFALWAZwc22B4XlIiIiItXGRx/B4MHBkp8P//oXfPcdnHYaZGWlOjqpTuKd0vIBoOipsmNL2MeBruWOSERERCTNTZ8Of/5zMNymZUu4665ghps6dVIdmVRX8U5p2StivUvywhERERFJX99+C9ddBxMmQJMm8Je/wEUXQU5OqiOT6i7ennoRERGRauuHH4IHRj3xBNSvD9deGzxIqnHjVEcmEihzUm9mjYBLgL5Asd+l7n5YBcUlIiIiknL5+XDzzfDww1CrFvzud/CHP0Dz5qmOTKS4RHrqnwOygBeA9RUbjoiIiEjqLVkSDK355z+DqSrPPx+uvhratEl1ZCKxJZLU7ws0c/fNFR2MiIiISCotXw5//Svcey9s3AhnnhncENupU6ojEyldIkn9+8BuwJcVHIuIiIhISqxeDXfeCXfcAWvWwEknBWPou3dPdWQi8UkkqT8TmGxmnwBLIt9w9xsrIigRERGRyrBuHdx3H9x+e9BLf+yxcOON0KvXzuuKpJNEkvpbgA7AT0DDiHKviIBEREREkm3jRhgzBm65JRg/P2QI3HQT7L13qiMTSUwiSf0ooIe7L6roYERERESSafNmGDcu6I3Pz4cBA4I55w88MNWRiZRPjQTqzAEq9CZZM8sys+lm9nK43dTM3jCz78PXJiXUG2Jms8xstpldFVEeV30RERGpHgoLgznme/aEc88NZrF54w145x0l9FI1JJLUPwZMMrOTzGxQ5FKOOC4FZkZsXwW85e7dgbfC7WLMLAu4HxgK9AROMrOe8dYXERGRqs8dJk6E3r3h1FOhXj2YNAk+/hh+/WswS3WEIhUjkeE3F4avt0aVO9C1rAczs/bAEQRj9a8Ii4cDA8P1cUAecGVU1X7AbHefEx7n6bDet3HWFxERkSrKHV57Da65Bj7/HHJz4ZlnYMQIqJFIl6ZImitzUu/uXSo4hruAPwANIspaFY3Zd/dFZtYyRr12QH7E9nygfxnqA2Bm5wHnAbRo0YK8vLwEL0OqooKCArUJ2YHahcSidpEeNm823n23BRMntuPbbxvRuvV6rrzyJw49dClZWc6771ZuPGoXUlkS6amvMGZ2JLDU3aeZ2cCyVo9RVuYZeNx9DDAGIDc31wcOLGsYUpXl5eWhNiHR1C4kFrWL1Jo3D/71L3joIVi6FLp2hQcegLPPrkt29m4Ej9ipfGoXUlnKnNSbWYlz0bv7tWU83AHA0WY2DKgDNDSzx4ElZtYm7GVvAyyNUXc+wdSaRdoDC8P1eOqLiIhIBtu6NbjZ9YEH4OWXg7Ijj4QLLoBDD9UwG6leEmnuHaKWfYDfAbuU9UDu/kd3b+/unQmmynzb3U8FJgFnhLudAbwYo/pnQHcz62Jm2WH9SeF78dQXERGRDLR8efDk19zcYH75jz6Cq66COXPgxRfh8MOV0Ev1k8iY+rOiy8xsCHBShUQUuA141szOBuYBI8PztAUecvdh7r7FzC4CXgeygIfd/ZvS6ouIiEjmmjo16JV/6inYsCGYivLGG+G446B27VRHJ5JaFTWmfgrwTHkO4O55BLPU4O6/AINj7LMQGBaxPRmYHGO/mPVFREQks6xfD08/Df/8J3z2GdSvD2ecAb/9bTBNpYgEEhlTHz1tZT3gZIrPRCMiIiKSsNmzg0T+kUdgxYrgoVH33gunnQaNGqU6OpH0k0hP/eyo7XXADLaPYRcREREps8JCeOWVYIjN669DzZrB0JoLLoCDD9aDokRKk8iYet16IiIiIhVmyRIYOzaYknLePGjXLhgrf8450KZNqqMTyQyJDL/JBs4E+gA5ke+5++kVEpWIiIhUae7wwQdBr/yECbB5MwweDHfeCUcfHfTSi0j8EvlPZjywJ/ASsKRiwxEREZGqbM0aeOKJIJn/6qtgfPwFF8D558Ouu6Y6OpHMlUhSfzjQxd1XVnAsIiIiUkV9801w4+v48UFi37cv/PvfcNJJwYw2IlI+iST18wDNBisiIiKl2rQJ/vOfoFf+v/8N5pI/8cSgZ75fP934KlKR4krqzWxQxOZ44EUzu5uo4Tfu/nYFxiYiIiIZaP58GDMm6IlfvBi6dIHbb4fRo6F581RHJ1I1xdtTPzZG2a1R2w5Ez2EvIiIi1cCyZUGv/IQJ8NZbsHUrDBsW9MoffjhkZaU6QpGqLa6k3t27JDsQERERySyLF8MLLwSJfF5ekMjvsgtceWUwHWUXZQ8ilUYTRomIiEjcFiyAiRODRP6994KpKXNz4eqrYcQI2HNPjZUXSQUl9SIiIlKqefPg+eeDRP7DD4OyXr3guuuCRL5nTyXyIqmmpF5ERER2MGfO9kT+00+Dsj594Oab4fjjNae8SLpRUi8iIiIAfP99kMRPmACffx6U7bUX/OUvQY98t26pjU9ESpZQUm9mhwKjgJbufpSZ7Q001JSWIiIimWXmzO2J/JdfBmX9+8Pf/hb0yOtmV5HMUOak3swuBi4FHgJGhMXrgXuA/SsuNBEREalo7sHTXZ97Lkjkv/02KD/gALjzTjjuOOjYMbUxikjZJdJTfxkw2N1/MrMrw7LvgNwKi0pEREQqjDt88cX2HvlZs4IbWw8+GO69F449Ftq1S3WUIlIeiST1DYD8cN3D11rApgqJSERERMrNHaZN257I//AD1KgBhxwCl10GxxwDrVunOkoRqSiJJPXvAlcBt0SUXQK8UyERiYiISEK2bg1mqilK5OfOhZo1YfBguOoqGD4cWrRIdZQikgyJJPUXAy+Z2blAAzObBawGjqrQyERERGSn5s6Ft98OlrfegkWLoFYtOOwwuP56OPpoaNo01VGKSLKVOal390Vmtg+wD9CJYCjOp+6+taKDExERkeIWLYJ33tmeyP/4Y1DesmUwtOaII+Coo6Bx45SGKSKVLJHZb37n7n8HPg2XovIr3P2OBI7XARgPtAa2AmPc/W4zawo8A3QGfgJOcPcVMeoPAe4GsoCH3P22sDyu+iIiIuls+XLIy9uexM+cGZQ3bgwDB8Lll8OgQXqqq0h1l8jwm2uBv8covwYoc1IPbAH+z90/N7MGwDQzewM4E3jL3W8zs6sIxvFfGVnRzLKA+4FDgfnAZ2Y2yd2/Dfcvtb6IiEi6WbMG3ntvexI/Y0Zw02u9esFsNWedFSTxffpAVlaqoxWRdBF3Um9mg8LVLDM7BIjsD+gKrEkkAHdfBCwK19eY2UygHTAcGBjuNg7IY8ekvB8w293nhDE+Hdb7Ns76IiIiKbV+PXz00fYk/tNPobAQsrNh//3hhhuCJH6ffYIyEZFYzN13vhdgZuGoPToC8yLecmAxcJu7TypXMGadCWbX6QXMc/fGEe+tcPcmUfuPAIa4+znh9mlAf3e/yMxW7qx+WH4ecB5AixYt9nr22WfLcwlSxRQUFJCTk5PqMCTNqF1ILPG2iy1bjO++a8D06U2YPr0xX3/diM2ba1CjhrPrrqvp23clffuuoFev1dSurdvVMp3+vZBYDjnkkGnuvndFHjPunnp37wJgZuPd/fSKDCI8bg7wPHCZu6+2+AYGxtopvl8pRTu7jwHGAOTm5vrAgQPLUl2quLy8PNQmJJrahcRSUrsoLAwe/FTUE//uu7B2bfBenz5w8cVBT/xBBxkNGzYCGhHMQyFVgf69kMqSyOw3yUjoaxEk9E+4+8SweImZtQln22kDLI1RdT7QIWK7PbCwDPVFREQqlHtwM2tREp+XByvCaRp23RXOOCNI4gcMgObNUxqqiFQhicx+c2NJ77n7tQkcz4CxwMyo2XMmAWcAt4WvL8ao/hnQ3cy6AAuAUcDJZagvIiJSLlu2wHffwSuvtGHMmCCRX7IkeK9TJzj22CCJP+QQaNs2tbGKSNWVyOw3HaK2WwMDgBcSjOEA4DTgKzObEZZdTZCMP2tmZxOM4R8JYGZtCaauHObuW8zsIuB1giktH3b3b8JjxKwvIiKSqKIEftq07cv06cHNrpBL69bB01sHDQqWLl1SHbGIVBeJDL85K7osnCv+pEQCcPf3iT02HmBwjP0XAsMiticDk2Ps90us+iIiIvEoLAyG0UQm8DNmwLp1wfv16sGvfgXnnQd77QXun3Laaf00V7yIpEQiPfWxTCF40JOIiEjGKSwMeuCnTi09gT/33CCB32svyM0tPk98Xt46JfQikjKJjKnvGlVUj2Ace36FRCQiIpJERQl8UfI+deqOCXzfvqUn8CIi6SaRnvrZBNNGFvVHrANmENyMKiIikjaiE/iiMfDRCfw558DeeyuBF5HMlciY+hrJCERERKQ8lMCLSHUWV1JvZoPi2c/d3y5fOCIiIqVzh/nzYdasIImfNStI3mfM2P5Qp8gEfq+9giReCbyIVGXx9tSPjWMfB6LH24uIiCRk7Vr43/+CpD0ygZ81a3vvO0CDBrDnnnD22dvHwO+6qxJ4Eale4krq3V0z7YqISIWL1ete9JofMf2CGXTuHPS2H3xwkLTn5gavrVujWWdEpNpLaEpLM+tOMC99O4InuT7l7t9XZGAiIlJ1RPa6R/a4x+p1z82FAQO2J+25udCtG9Stm7r4RUTSXSJTWh4FPAG8DMwFcoGpZnaau0+q4PhERCRDbN0KCxYUT9rV6y4iUjkS6am/FRju7u8UFZjZQOA+QEm9iEgVtnFjMFxm3rxg+fFH9bqLiKSDRJL69sB7UWXvh+UiIpKh3GHZsu0Je6xlyZLiddTrLiKSHhJJ6mcA/wfcHlF2RVguIiJpat26YBhMSQl7fn7QEx+pXj3o2DFYeveGDh22b3fsCO3bQ506qbkeERHZLpGk/rfAS2Z2KZAPdADWAkdXZGAiIhK/rVth8eLSe9l/+aV4HTNo2zZIzvfaC449dsekvWlT9baLiGSCRJ4o+52Z7QbsC7QFFgKfuPvmig5ORKS627IlSMaXLIGlS4OlaD1ybPuCBbA56l/hBg2gU6cgOe/fP3iNTNrbtYNatVJzXSIiUrESmtLS3bcQjKMXEZEyWrt2xwQ9VtK+dCn8/HMw1j1azZpBUt6xIxxwwPZEPTJpb9So8q9NRERSI5EpLQ8BfnL3H82sNcHY+kLgandfXNEBioiku61bYfny0hP0yLK1a2Mfp2FDaNkSWrWCHj3goIOC7aKyyNfGjTUsRkREtkukp/4B4PBw/Y7wdQswBo2rF5EM5g5r1sCKFbByZfAauV70OnPmbtxyy/Zk/eefobBwx+NlZUGLFtuT8W7dYifoRYtuOBURkUQlktS3c/d5ZlaTILnvBGwiGFsvIpJSmzfHTsRLS9KL1leuDHrdS2IW9JDXrduQjh2DqRz7998xQS9ab9oUatRI+iWLiIgklNSvNrNWQC/gW3cvMLNsQLdbiUjCNm+GgoKdL6tWlZ6klzS0pUjt2tCkSbA0bhwk4Lm527cj34t+bdgwSNLz8j5h4MCBSf08REREyiKRpP5e4DMgG7gsLDsA+K6CYhKRNOYOmzaVnnivWRNfgh65RM+PXppGjYon4N27x07EY5VpiIuIiFRFiUxpebuZvQBscfc5YfEC4JwKjUxE4rJ1K2zYAOvXb39N1npRAr5lS/zx1asHOTnblwYNggS7ffvi5dFLgwY7ltWvH5RnZSXt4xQREclIicx+kw2cCpxkZkXz1D8D3FzBsWFmQ4C7gSzgIXe/Lep9C98fBqwDznT3z+OpK1IW7kEiu2lTsGzcuH091lLR72/YUHKivWlT4tdlBnXrBkudOjuuF/VsF5WXJQHPyQkSeiXgIiIiyZfI8Jt/ArnAJcBcghtl/0gwK87oigrMzLKA+4FDgfnAZ2Y2yd2/jdhtKNA9XPqHsfWPs66ksZdfhoULg0R68+bgtWgpbbss+8ZTd+PGAyksLF/iXJqaNYMx3tnZpS9FCXZpCXgi67VqaVpEERGRqiCRpP4YYBd3Xxluf2tmnwCzqcCkHugHzC4a4mNmTwPDgcjEfDgw3t0d+NjMGptZG6BzHHUljd1+O7xfyuPNsrKChLhoqVUr/u06deLfd9GiRXTr1qHEZDuehLykfWrV0swoIiIiUjESSeoXA/WAlRFldYFFFRFQhHZAfsT2fILe+J3t0y7OugCY2XnAeQAtWrQgLy+vXEFLxbj88lpceqmRleXblpo1g9caNbzSkuGCggJycn4oU52iMe4bNiQpKEm5goIC/VshO1C7kFjULqSyxJXUm9mgiM3HgNfM7F6CZLkDcCEwvoJjizUoIPph6SXtE0/doNB9DMGDs8jNzXVNUyeR8vLyNHWh7EDtQmJRu5BY1C6kssTbUz82RtnVUdu/AW4vXzjFFP1gKNKeHR9wVdI+2XHUFRERERGpEuJK6t29y872CW9OrUifAd3NrAvBlJmjgJOj9pkEXBSOme8PrHL3RWa2LI66IiIiIiJVQiJj6osxsz2A04FTgLbljijk7lvM7CLgdYJpKR9292/M7Pzw/QeByQTTWc4mmNLyrNLqVlRsIiIiIiLpJKGk3sxaEPR8nwH0Bt4HLq3AuABw98kEiXtk2YMR604wnj+uuiIiIiIiVZEFeXEcO5rVAo4GzgQOJ+gdfwq4DNjN3ZcmJ8TKY2ZrgFmpjkPSSnPg51QHIWlH7UJiUbuQWNQuJJZcd29QkQcsS0/9EmAr8ChwXcSTWy+oyIBSbJa7753qICR9mNlUtQmJpnYhsahdSCxqFxKLmU2t6GOWZbbvL4HGBDek7mNmTSo6GBERERERKbu4k3p3HwjsAkwBfgcsNrOXgPpAraREJyIiIiIiO1Wm53K6+1x3v8nduwODCZ4iuxX4wsz+mowAK9mYVAcgaUdtQmJRu5BY1C4kFrULiaXC20XcN8qWeACzOsCxwOnuPrRCohIRERERkbiVO6kXEREREZHUKtPwm0xkZvuYWaGZjYgoG2Jms8xstpldVUI9M7N7wn2+NLNflaW+pCczOyX8Pr80sw/NrHfEe2oXstPvUW2gejCzDmb2jpnNNLNvzOzSsLypmb1hZt+HrzEnjSipLcRbX9KXmWWZ2XQzezncVpsQzKyxmU0ws+/Cfzf2q/S24e5VdiF4muzbBA+hGhFR9gPQFcgGvgB6xqg7DHgVMGBf4JOy1NeSnguwP9AkXB9a1u9V7aJqL/F8j2oD1WMB2gC/CtcbAP8DegJ/Ba4Ky68Cbi9LO4qnvpb0XoArgCeBl+P9TtUmqv4CjAPOCdezCWaMrNS2UdV76i8GngciH4zVD5jt7nPcfRPwNDA8Rt3hwHgPfAw0NrM2ZagvacjdP3T3FeHmx0D7cF3tQiC+71FtoBpw90UePo/F3dcAM4F2BN/puHC3ccAxMaqX1hbiqS9pyszaA0cAD0UUq01Uc2bWEDgYGAvg7pvcfSWV3DaqbFJvZu0IbuB9MOqtdkB+xPb8sCxaSfvFW1/S39kEPa6gdiGBeL5HtYFqxsw6A32BT4BW7r4IgsQfaBmjSmltIZ76kr7uAv5AMPNfEbUJ6QosAx4Jh2Y9ZGb1qeS2UWWTeoL/8K5098Kocouxb6y7hUvaL976ksbM7BCCpP7KoqIYu6ldVD/xfI9qA9WImeUQ/MX3MndfHW+1GGVqCxnOzI4Elrr7tESqxyhTm6g6agK/Av7p7n2BtQTDZeJRYW2jZiKV0pWZXQicG242Ap42M4DmwDAz20LwC6hDRLX2wMIYhytpv+w460uaiGoXwwjaw0PAUHf/JSxXuxCIrx2oDVQTZlaLIKF/wt0nhsVLzKyNuy8Kh10tjVG1tHYUT31JTwcAR5vZMKAO0NDMHkdtQoLvd767fxJuTyBI6iu1bVSpnnp3v9/d+4RLF3fv7O6dCT7cC9z9P8BnQHcz62Jm2cAoYFKMw00CTg9nutgXWBX+6SPe+pImItsFwQ/ZicBp7v6/iN3ULgTi+x7VBqoBC3qExgIz3f2OiLcmAWeE62cAL8aoXlpbiKe+pCF3/6O7tw/zilHA2+5+KmoT1Z67LwbyzSw3LBoMfEtlt43Kvjs4FQvwKOHsN+H2MIKZDH4A/hRRfj5wfrhuwP3hPl8Be++svpb0Xwh66FcAM8JlqtqFlqg2ssP3qDZQ/RbgQII/gX8Z8e/FMKAZ8BbwffjaNNy/LTB5Z22hpPpaMmsBBrJ99hu1CS0AfYCp4b8Z/wGaVHbb0MOnREREREQyXJUafiMiIiIiUh0pqRcRERERyXBK6kVEREREMpySehERERGRDKekXkREREQkwympFxERERHJcErqRUREREQynJJ6ERHJKGbmZrbWzG6p4OO+bWYbzOz9ijyuiEhlUFIvIlIJzKwgYtlqZusjtk9JdXyJMrOfzOzXKTh1b3f/U1QsQ83sDTO7M1YFMzvZzKaGn/kiM3vVzA4set/dBxE8PVhEJOMoqRcRqQTunlO0APOAoyLKnkh1fLGYWc0MO8dhwAigQYzzXAHcBdwKtAI6Ag8Awyvw/CIiKaOkXkQkDZhZWzN73syWmdmPZnZJxHs/mdnvzezLcNjJWDNrFfY0rzGzN82sSdT+fzSzb81shZk9YmZ1ynCuK83sS2CtmdU0s6vM7IfwXN+a2bHhvo8RJMcvhb3ffwjL3cy6RRzzUTO7uZRzdCwpnjJ6EPg38G7UZ9sIuBG40N0nuvtad9/s7i+5++8TPJeISFpRUi8ikmJmVgN4CfgCaAcMBi4zs8MjdjseOBToARwFvApcDTQn+Lc8OhE+BTgc2CWsc00ZznUScATQ2N23AD8ABwGNgBuAx82sjbufRvG/Ovy1DJdddI6mwAs7iScu7j7L3U9w9/FRb+0H1AnPIyJSJSmpFxFJvX2AFu5+o7tvcvc5BD3OoyL2udfdl7j7AuA94BN3n+7uGwmS1b5Rx7zP3fPdfTlwC0ESHe+57gnrrgdw9+fcfaG7b3X3Z4DvgX7lvOZ73D0f6BVHPOXVDPg5/IEiIlIlJX28pIiI7FQnoK2ZrYwoyyJI3ossiVhfH2M7J+qY+RHrc4G2ZThXZF3M7HTgCqBzWJRD8BeC8ig6RzzxlNcvQHMzq6nEXkSqKiX1IiKplw/86O7dK/CYHSLWOwILy3AuL1oxs04EPeeDgY/cvdDMZgAWvW+EdUC9iO3WwPwSzpGMa4/2EbABOAaYkMTziIikjIbfiIik3qfA6vDm0bpmlmVmvcxsn3Ic80Iza29mTQnG3j+T4LnqEyTgywDM7CyCITNFlgBdo+rMAE4Ojz0EGFBKnMm49mLcfRVwLXC/mR1jZvXMrFY4BWZZ7gMQEUlbSupFRFLM3QsJbn7tA/wI/Aw8RHBjaqKeBKYAc8Ll5kTO5e7fAv8g6O1eAuwBfBCxy1+Aa8xspZn9Liy7NDzHSoIbdv9TUpBJuvZY57mDYAjRNQQ/UPKBi0qLTUQkk5h7rL+ciohIpjKzn4Bz3P3NVMeSDGa2AdhIcLPtnyvwuG8A+wKfuvvgijquiEhl0Jh6ERHJKO5eZ+d7JXTcQ5NxXBGRyqDhNyIiIiIiGU7Db0REREREMpx66kVEREREMpySehGptszsUTO7uZLO9aCZlXhTp5m5mXWLZ18REZFoSupFJGOY2fdm1j2q7C4zW2FmH5lZu4jyU8zs7sqPMjZ3P9/dbyrrvmY20MyiH9xU7ZhZ5/CHT5WY4MHM8szsnFTHISJVh5J6Eckkk4FhRRtm1g/Yi+CJpe8DfwzLGwG/I3jgkMQpVsJc1iS6qiTd5VEZn4E+ZxGJpqReRDJJsaQe6AK87+4bgbfY/mTTW4C/hU8S3ZkmZvaKma0xs0/MbBeI3TMc2btqZmea2Qdmdmf44KU5ZrZ/WJ5vZkvN7IyIusWG+pjZ781skZktNLPRkQEV7Wtm9YFXgbZmVhAubc1snZk1i9h/LzNbZma1oi/OzGqY2VVm9oOZ/WJmz4ZPmY28xrPNbB7wdtR1LQeuN7NGZjY+PMdcM7vGzGrE+ByWA9dHnb+tma0vOmdY1tfMfg6f6trNzP5rZqvCsmeI7d3wdWX4OewXHmu0mc0M/1rzupl1ijiPm9kF4V941pjZTWa2S/hXndXhZ5Ed7jvQzOab2dVhHD+Z2SkRx6ptZn83s3lmtiQcIlU3qu6VZrYYeMTMmpjZy+FntiJcbx/ufwtwEHBfeC33lbG9FX0vJcYkItWPknoRySR5wN5hsgvwDXBQmMgMBr4xs72BXHd/Ms5jngTcADQBZhP8IIhXf+BLoBnBE1yfBvYBugGnEiRtOdGVzGwIwV8SDgW6A7+OdXB3XwsMBRa6e064LCT4HE6I2PVU4Gl33xzjMJcAxwADgLbACuD+qH0GALsBh0dc1xygJcHncS/BE167hvueDpwV9TlE7h95DQsJnkZ7fETxycCEMN6bCJ582wRoH54rloPD18bh5/CRmR0DXA0cB7QA3gOeiqo3hOCvOfsCfwDGEDzltgPQi+D7L9IaaA60A84AxphZbvje7UAPgiffdgv3uTaqblOgE3Aewf9fHwm3OwLrgfvCz+RPYawXhddyUQnXHC36c95ZTCJSjSipF5GMEfbIf0CQwOPuXwPPAx8TJE63A3cDl5jZJWb2rpk9YWaNSznsRHf/1N23AE8QJEjx+tHdH3H3QuAZgkTxRnff6O5TgE0EyVa0E4BH3P3rMHG/vgznBBhHkMhjZlkEieljJez7G+BP7j4//PyuB0ZY8eEb17v7WndfH24vdPd7w89kE3Ai8Ed3X+PuPwH/AE6LqL9t/4hjRHoyjBEzM2BUWAawmSDxbevuG9z9/fg/Bn4D/MXdZ4ax3gr0ieytB25399Xu/g3wNTDF3eeEf8V5Fegbdcw/h9/ff4FXgBPCmM8FLnf35e6+JjzXqIh6W4Hrwrrr3f0Xd3/e3deF+99C8IOoPCK/lw1xxCQi1YiSehHJNMWG4Lj7ne7e291PJEg+3yP4t+08guR/JnBVKcdbHLG+DtihZ70USyLW14fxRJfFOl5bID9ie24ZzgnwItDTzLoS9PavcvdPS9i3E/CCBUOEVhJ8HoVAq4h98qPqRG43B7KjYpxL0CtcUv1oE4D9zKwtQY+7E3xPEPSeG/CpmX1jUUORdqITcHfEtS0PjxUZW/T3Udr3syL8kVVkLsF31QKoB0yLONdrYXmRZe6+oWjDzOqZ2b/C4UqrCYYPNQ5/hCUq8nOOJyYRqUaU1ItIpnmFYEhKMWbWiqDn9kaCYRVfhsM7PgP2TOA8RcldvYiy1gkcJ5ZFBL36RTqWsu8OTwgMk8dnCYaRnEbJvfQQJIJD3b1xxFLH3ReUco7I7Z/Z3pseGW9p9aPjXUkwxOYEgqE3T3n45EN3X+zu57p7W4Lv7wELp/YsJabIa/tN1LXVdfcPS4unFE0ihnZBcJ0LCT6D9cDuEedp5O6RPwii4/s/IBfo7+4N2T58yErYP572Fv297CwmEalGlNSLSEZx93xgtZn1inrrDoLhD+uAH4F9wvHsAwnGIZf1PMsIEtdTzSwr7EHepVzBb/cscKaZ9TSzesB1pey7BGhmwYw+kcYDZwJHA4+XUv9B4JaiISlm1sLMhscbaDi06NnwGA3C41yxk3PG8iTBWPzj2T70BjMbWXQDKcF4fyf4S0K0ZQRDXLpGlD0I/NHMdg+P1cjMRpYxrmg3mFm2mR0EHAk85+5bgX8Dd5pZy/Bc7czs8FKO04Ag6V5pwU3C0d/xkshrKWt7SzAmEanClNSLSCaaDBxRtGFmhxDcQPkCQDgU5RWCntxDgNsSPM+5wO+BX4DdgUR7gItx91eBu4C3CW7OfbuUfb8juPlzTjjMom1Y/gFBkvt5OM69JHcDk4ApZraG4P6D/mUM+WKCnuQ5BFOHPgk8XMZjTCK4KXiJu38RUb4P8ImZFYT7XOruP0ZXDn+s3QJ8EH4O+4bf9+3A0+EQl6+J8VecMlhM8MNiIcH9FeeHnz/AlQTf1cfhud4k6IkvyV1AXYIe9Y8JhsZEupvg3oYVZnZPWFbW9lbWmESkCrPwL6AiIhnDzAYQ3JBa3hsPM5qZvQ086e4PpTqWTGdmA4HH3b39TnYVEUlLeniFiGSiD4DXUx1EKpnZPsCvgLiH0oiISNWl4TciknHCqRNvTXUcqWJm4wiGWlwWTmUoIiLVnIbfiIiIiIhkOPXUi4iIiIhkOCX1IiIiIiIZTjfKRmjcuLF36xbrmSdSXa1du5b69evvfEepVtQuJBa1C4lF7UJimTZt2s/uXqFPgFZSH6FVq1ZMnTo11WFIGsnLy2PgwIGpDkPSjNqFxKJ2IbGoXUgsZja3oo+Z1sNvzGyImc0ys9lmdlWM939vZjPC5WszKwyf3IeZ/WRmX4XvKVMXERERkSorbXvqzSwLuB84FJgPfGZmk9z926J93P1vwN/C/Y8CLnf35RGHOcTdf67EsEVEREREKl0699T3A2a7+xx33wQ8TekPWTmJ4FHqIiIiIiLVSjon9e2A/Ijt+WHZDsysHjAEeD6i2IEpZjbNzM5LWpSSFNfnXc+t792KnqMgIiIisnNpO/wGsBhlJWV4RwEfRA29OcDdF5pZS+ANM/vO3d/d4SRBwn8eQIsWLcjLyytn2FJeW30rH3z3AW8ufZMZ/5vB+V3Pp4al5vdnQUGB2oTsQO1CYlG7kFjULqSypHNSPx/oELHdHlhYwr6jiBp64+4Lw9elZvYCwXCeHZJ6dx8DjAHIzc113aGeHgYOHMhlr13GvZ/eS71m9Xjo6IeoWaPym6tmLZBY1C4kFrULiUXtQipLOif1nwHdzawLsIAgcT85eiczawQMAE6NKKsP1HD3NeH6YcCNlRK1VIgaVoO7h9xN83rNuS7vOlZuWMnTI56mTs06qQ5NREREJO2k7Zh6d98CXAS8DswEnnX3b8zsfDM7P2LXY4Ep7r42oqwV8L6ZfQF8Crzi7q9VVuxSMcyMawdcy71D7+XFWS8y5PEhrN64OtVhiYiIiKSddO6px90nA5Ojyh6M2n4UeDSqbA7QO8nhSSW5qN9FNK3blDP+cwaHjDuEV095lZb1W6Y6LBEREZG0kbY99SKRTt7jZF4c9SIzl83koEcOYu7KCn8Qm4iIiEjGUlIvGWNY92FMOW0KSwqWcOAjBzJz2cxUhyQiIiKSFpTUS0Y5sOOBvHvWu2wu3MxBjxzEpws+TXVIIiIiIimnpF4yzp6t9uSD0R/QsHZDBo0bxFtz3kp1SCIiIiIppaReMtIuTXfh/dHv06VJF4Y9OYyJMyemOiQRERGRlFFSLxmrbYO2vHvmu+zVZi9GPjeShz5/KNUhiYiIiKSEknrJaE3qNuGN097g0K6Hcu5L5/LXD/6a6pBEREREKp2Sesl49bPrM+mkSZy4+4lc+eaVXPnGlbh7qsMSERERqTRp/fApkXhlZ2XzxHFP0KROE/764V9Zvn45Dx75IFk1slIdmoiIiEjSKamXKiOrRhYPHPEAzes15+b3bmb5huU8edyT1K5ZO9WhiYiIiCSVht9IlWJm3DToJu48/E4mzpzIEU8ewZqNa1IdloiIiEhSKamXKumyfS9j3DHjyPspj8HjB/Pzup9THZKIiIhI0iiplyrr9N6nM/HEiXy55EsOfuRg5q+en+qQRERERJJCSb1UaUfnHs3rp77OgjULOODhA5j186xUhyQiIiJS4ZTUS5U3oPMA8s7IY/3m9Rz0yEF8vujzVIckIiIiUqGU1Eu10LdNX94f/T51a9Vl4KMDyfspL9UhiYiIiFQYJfVSbfRo1oMPRn9Ah0YdGPL4EF787sVUhyQiIiJSIZTUS7XSvmF73j3zXXq37s3xzx7PuBnjUh2SiIiISLkpqZdqp1m9Zrx1+lsc0uUQznzxTO746I5UhyQiIiJSLkrqpVrKyc7h5ZNeZkTPEfzflP/jT2/9CXdPdVgiIiIiCVFSL9VW7Zq1efr4pzn3V+dy6/u38ttXfkvh1sJUhyUiIiJSZjVTHYBIKmXVyOJfR/6LZnWbcdsHt7FiwwoeO/YxsrOyUx2aiIiISNyU1Eu1Z2b85dd/oVm9Zvz+jd+zcsNKnj/heXKyc1IdmoiIiEhcNPxGJPS7/X/H2KPH8uacNzn0sUNZvn55qkMSERERiUtaJ/VmNsTMZpnZbDO7Ksb7A81slZnNCJdr460rEsvovqOZMHICny/6nIMfOZhlG5elOiQRERGRnUrbpN7MsoD7gaFAT+AkM+sZY9f33L1PuNxYxroiOzh2t2N59ZRXmbtqLpfMuISFaxamOiQRERGRUqVtUg/0A2a7+xx33wQ8DQyvhLoiDOoyiLdPf5uVm1Zy4oQT2Vy4OdUhiYiIiJQonW+UbQfkR2zPB/rH2G8/M/sCWAj8zt2/KUNdzOw84DyAFi1akJeXV/7Ipcq4sNOF/OPHf3Dao6dx/i7npzocSRMFBQX6t0J2oHYhsahdSGVJ56TeYpRFPx3oc6CTuxeY2TDgP0D3OOsGhe5jgDEAubm5PnDgwETjlaooD9Y2XcuD0x7k5ANP5ujco1MdkaSBvLw89G+FRFO7kFjULqSypPPwm/lAh4jt9gS98du4+2p3LwjXJwO1zKx5PHVF4nXnkDvZq81enPGfM/hxxY+pDkdERERkB+mc1H8GdDezLmaWDYwCJkXuYGatzczC9X4E1/NLPHVF4lWnZh2eG/kcACOfG8mGLRtSHJGIiIhIcWmb1Lv7FuAi4HVgJvCsu39jZuebWdHg5hHA1+GY+nuAUR6IWbfyr0Kqii5NujDumHFMWzSNK16/ItXhiIiIiBSTzmPqi4bUTI4qezBi/T7gvnjripTH0blH8/v9f8/fPvwbB3Y8kJP3ODnVIYmIiIgAadxTL5KObhl0Cwd2PJDzXjqPmctmpjocERERESCJSb2ZZZnZaDOrnaxziFS2Wlm1eGbEM9TPrs+I50awdtPaVIckIiIikryk3t0LgTvcfWOyziGSCm0btOXJ455k5rKZ/Obl3+Aec7ZUERERkUqT7OE3L5nZUUk+h0ilG9x1MDcMvIEnvnqCMdPGpDocERERqeaSfaNsHWCCmX1E8ITXbV2a7n56ks8tklR/OvhPfJD/AZe8dgn7tNuHX7X5VapDEhERkWoq2T31XwO3Au8As4EfIhaRjFbDavD4cY/Tsn5LRjw7gpUbVqY6JBEREammktpT7+43JPP4IqnWvF5znh3xLAc/ejBnvXgWE0+YSPg8NBEREZFKk/QpLc3sEDN72MxeD18HJfucIpVpvw778bdD/8Z/vvsPd3x0R6rDERERkWooqUm9mZ0DPAMsBiYCi4AnzezcZJ5XpLJd2v9Sjt/teK5880ren/d+qsMRERGRaibZPfV/AA5196vd/V/u/ifgsLBcpMowM8YePZYuTbpw4oQTWbp2aapDEhERkWok2Ul9M+DbqLJZQNMkn1ek0jWq04jnRj7H8vXLOWXiKRRuLUx1SCIiIlJNJDupfx+4w8zqAZhZfeBvwIdJPq9ISvRp3Yf7ht7Hm3Pe5KZ3b0p1OCIiIlJNJDupPx/YA1hlZkuAlUBv4DdJPq9IyozuO5ozep/Bjf+9kSk/TEl1OCIiIlINJC2pN7MsYCgwBOgCHAV0cfcB7r4wWecVSTUz44EjHmD3lrtzysRTmL96fqpDEhERkSouaUm9uxcCd7j7Bnef7+6furuyG6kW6tWqx4SRE9iwZQMnTjiRzYWbUx2SiIiIVGHJHn7zkpkdleRziKSl3Oa5PHTUQ3yY/yFXvXlVqsMRERGRKiypT5QF6gATzOwjIB/wojfc/fQkn1sk5U7sdSLvz3ufOz6+gwM6HsBxux2X6pBERESkCkp2Uv91uIhUW38/7O98suATznrxLPZstSfdmnZLdUgiIiJSxSQtqQ9vlO0KnOfuG5N1HpF0V7tmbZ4b+Rx9/9WXkc+N5MPRH1K3Vt1UhyUiIiJVSLJvlD0M2Jqsc4hkik6NO/HYsY8xY/EMLn3t0lSHIyIiIlVMsm+UvRO4wcxqJfk8ImnviB5H8McD/8i/P/83478Yn+pwREREpApJdlJ/MfB7YI2Z5ZvZvKIlyecVSUs3HnIjAzoN4PyXz+frpbrdRERERCpGsm+UPTXJxxfJKDVr1OSp45+i77/6MuLZEXx27mc0qN0g1WGJiIhIhktqT727/7ekJZnnFUlnbRq04ekRT/P98u857+XzcPedVxIREREpRVKTejOrbWa3mNkcM1sVlh1mZhfFWX+Imc0ys9lmtsPTe8zsFDP7Mlw+NLPeEe/9ZGZfmdkMM5tacVclUn4DOw/kpkNu4umvn+afU/+Z6nBEREQkw1XGjbK9gFPY/uCpb4Df7qxiOCXm/cBQoCdwkpn1jNrtR2CAu+8J3ASMiXr/EHfv4+57J34JIslx1YFXMaz7MC5//XI+W/BZqsMRERGRDJbspP5Y4GR3/4hwakt3XwC0i6NuP2C2u89x903A08DwyB3c/UN3XxFufgy0r7DIRZKshtVg/DHjaZ3TmpHPjWT5+uWpDklEREQyVLJvlN0UfQ4zawH8EkfddkB+xPZ8oH8p+58NvBqx7cAUM3PgX+4e3YtfFM95wHkALVq0IC8vL47QpLooKChIepu4qutVXDLjEo586Ehu7nUzNSzZv7WlvCqjXUjmUbuQWNQupLIkO6l/DhhnZpcDmFkb4C6CXvedsRhlMe8oNLNDCJL6AyOKD3D3hWbWEnjDzL5z93d3OGCQ7I8ByM3N9YEDB8YRmlQXeXl5JLtNDGQgW1pt4ZLXLuGzWp9x5YFXJvV8Un6V0S4k86hdSCxqF1JZkt0leDXwE/AV0Bj4HlgI3BBH3flAh4jt9mHdYsxsT+AhYLi7b/sLgLsvDF+XAi8QDOcRSUsX9buIE3Y/gT+9/Sf++5MmhxIREZGySfaUlpvc/TJ3zwFaAQ3c/fJwjPzOfAZ0N7MuZpYNjAImRe5gZh2BicBp7v6/iPL6ZtagaB04DNCTfiRtmRkPHfUQuzTdhVHPj2JxweJUhyQiIiIZpNIG77r7Mi/DhNzuvgW4CHgdmAk86+7fmNn5ZnZ+uNu1QDPggaipK1sB75vZF8CnwCvu/lqFXYxIEjSo3YAJIyewasMqTn7+ZAq3FqY6JBEREckQyR5TXy7uPhmYHFX2YMT6OcA5MerNAXpHl4ukuz1a7cEDRzzAWS+exXV513HzoJtTHZKIiIhkAE2zIZJmzuxzJqP7jOaW927h1e9f3XkFERERqfaU1IukofuG3ceerfbk1BdOZd6qeakOR0RERNJcpST1ZtbQzP5iZi+b2T1m1rYyziuSqerWqsuEkRPYXLiZE547gU2F8dxbLiIiItVVZfXU3w8UAPcAa4EJlXRekYzVvVl3Hhn+CJ8s+ITzXz5fN86KiIhIiZKS1JvZnUVTSoY6Are5+xTgZmDXZJxXpKo5vufxXHvwtTwy4xGOf/Z41m1el+qQREREJA0lq6d+KpBnZieG288D083sceBzYFySzitS5dxwyA3cO/ReJs2axODxg/l53c+pDklERETSTFKSend/AhgEHGhmrxPMNV/08KhT3f3yZJxXpKq6qN9FTDhhAtMXTeeAhw9gzoo5qQ5JRERE0kjSxtS7+yp3vxi4GhgLnApMcvfPknVOkarsuN2O463T32LZ2mXsN3Y/pi2cluqQREREJE0ka0x9m3CWm5eBE4DhwALgYzM7OhnnFKkODuh4AB+M/oC6Nesy4NEBvDZbD0oWERGR5PXUTwA2APcCBtzr7vcDhwMnmNlLSTqvSJW3W4vd+Ojsj+jerDtHPnkkj0x/JNUhiYiISIrVTNJxdwMGuvtmM/sv8DGAuy8BTjWzgUk6r0i10KZBG/575n8Z8ewIRk8azfzV87nm4Gsws1SHJiIiIimQrJ768cCbZnYLMAV4NPJNd89L0nlFqo2GtRvy8skvc3rv07k271p+8/Jv2LJ1S6rDEhERkRRISk+9u19mZvsAXYAn3f2bZJxHpLrLzsrm0eGP0r5Be259/1YWFSzi6eOfpn52/VSHJiIiIpUombPffObuz0Yn9GbW08x+b2bvmNlVyTq/SHVhZtwy+BYeGPYAk7+fzKDxg1i6dmmqwxIREZFKlLSkvoiZ1TWzI83sn2b2E/AU0Ay4Dvhbss8vUl38dp/fMvGEiXy15Cv2H7s/s5fPTnVIIiIiUkmSNaXlLmZ2sZm9BswDzgFmAAe6e293v8rd33X3wmScX6S6Gr7rcN46/S1WbljJ/mP359MFn6Y6JBEREakEyeqpfxXoCtwBtHP3Y9z9X+4+P0nnE5HQfh3248OzPyQnO4dDxh3CK/97JdUhiYiISJIlJal39x7ufrm7T3H3Tck4h4iUrEezHnx09kfs1nw3hj89nIc+fyjVIYmIiEgSJX1MvYikRqucVuSdmcehuxzKuS+dy3XvXIe7pzosERERSQIl9SJVWE52DpNGTeKsPmdx47s3cs6kc9hcuDnVYYmIiEgFS9YTZUUkTdTKqsXYo8fSoWEHbnz3RhYVLOLZkc+Sk52T6tBERESkgqinXqQaMDNuOOQGxhw5htd/eJ2Bjw5kScGSVIclIiIiFURJvUg1cu5e5/LiqBeZ+fNM9hu7H//75X+pDklEREQqQFon9WY2xMxmmdnsWE+ftcA94ftfmtmv4q0rUl0d2eNI3jnjHQo2FbD/2P35eP7HqQ5JREREyiltk3ozywLuB4YCPYGTzKxn1G5Dge7hch7wzzLUFam2+rXrx4dnf0jjOo0ZNG4Qk2ZNSnVIIiIiUg5pm9QD/YDZ7j4nnOv+aWB41D7DgfEe+BhobGZt4qwrUq11a9qND8/+kF4te3HsM8fy4NQHUx2SiIiIJCidk/p2QH7E9vywLJ594qkrUu21rN+Sd854h6HdhvLbV37Ln976k+ayFxERyUDpPKWlxSiLzjZK2ieeusEBzM4jGLpDixYtyMvLK0OIUtUVFBRUizZxeZvLoQBuff9Wpn4/ld/3+D01a6TzPw+pVV3ahZSN2oXEonYhlSWd/689H+gQsd0eWBjnPtlx1AXA3ccAYwByc3N94MCB5Qpaqpa8vDyqS5sYNHAQN797M9fmXYvXd54/4Xka1G6Q6rDSUnVqFxI/tQuJRe1CKks6D7/5DOhuZl3MLBsYBUTfzTcJOD2cBWdfYJW7L4qzrohEMDP+PODPPHz0w7z949sMeHQAi9YsSnVYIiIiEoe07al39y1mdhHwOpAFPOzu35jZ+eH7DwKTgWHAbGAdcFZpdVNwGSIZ56y+Z9E6pzUjnxvJfmP347VTX2PX5rumOqxyKdxayKbCTWzeuplNhZuC9cKI9TKWz104l9w1ubRp0CbVlyYiIgKkcVIP4O6TCRL3yLIHI9YduDDeuiISn6Hdh5J3Zh5HPHkEBzx8AJNGTeKAjgdU2PE3FW5izcY1FGwqKLas2bRjWcGmgmDfzcH6hi0bypyYb/WtFRZ7kXvuvIeh3Ydydt+zOaL7EdTKqlXh5xAREYlXWif1IpI6e7fdm4/O/oghjw/h14/9mrFHj6V3q96lJ95F25tLKA+XzVs3xx1H/Vr1ycnOISc7h/rZ9alXqx61atSibs26NKrdiFpZtcjOyiY7K5taNbavR2+Xul8Z35v09iS+q/0d474Yx8v/e5mW9Vty2p6nMbrvaHq20CMxRESk8impF5ESdW3SlQ/P/pCjnjqKUyaestP9ixLwBrUbbEvEm9VrRqfGnYLy7O3lkUtJ5fVq1SOrRlYlXGnZdKrfiTMGnsFNg27itdmv8fD0h7n7k7v5x0f/YN/2+zK6z2hO7HUiDWs3THWoIiIp4+4sXbuUOSvmMGfFHFZtXMUF+1yQ6rCqLCX1IlKq5vWa89bpb/HSrJeoYTWKJ+MRyXu9WvWoYel8733Fq1mjJkf2OJIjexzJ0rVLefzLxxk7fSznvXwel71+GSN7jmR039Ec1PEgzGLNtCsiktnWbV7Hjyt+ZM6KOfy48sdtCXzR9rrN67btW7dmXX6792/172GSKKkXkZ2qV6seJ/Y6MdVhpLWW9VtyxX5XcPm+l/Ppgk95ePrDPPX1U4z7YhzdmnZjdJ/RnN77dNo11HPwRCRzFG4tZOGahSUm7YsLFhfbv36t+nRt0pVuTbtx2C6H0aVxF7o26UrXJl3p3LizEvokUlIvIlKBzIz+7fvTv31/7jj8Dp6f+TwPT3+Yq9++mmveuYYh3YZwdt+zObLHkWRnZac6XBERVm1YVSxRj0zc566ay6bCTdv2rWE16NCwA12bdOWI7kcUS9q7NulK83rNlbiniJJ6EZEkqZ9dn9N7n87pvU9n9vLZPDrjUR6d8SjHP3s8Leq12HZz7e4td091qCJShW0u3My8VfNK7G1fvn55sf2b1m1Kl8Zd6NO6D8ftdlyxxL1jo46a7StNKakXEakE3Zp24+ZBN3PDwBuY8sMUxk4fy72f3ssdH99Bv3b9GN1nNKN6jaJRnUapDlVEMtQv637hu5+/Y9Yvs4q9/rD8Bwq9cNt+tWrUonPjznRt0pV+7frRtUnXbYl7lyZdaFynceouQhKmpF5EpBJl1chiaPehDO0+lGVrl/HEV08wdvpYzn/lfC5//XJG9BzB6L6jGdBpgP6ELSI72LJ1C3NWzGHWz8UT9+9+/o5f1v+ybb/srGx6NOvBHi33YMRuI+jWtNu23va2Ddqm5cxiUj5K6kVEUqRF/RZctu9lXNr/UqYunMrD0x/mya+f5LEvH2OXJrtwVp+zOKPPGbRv2D7VoYpIJVuxfsX2HvefZ/HdL8Hr7OWziz3ro1X9VuQ2z+X43Y4nt3kuuzbfldxmuXRu3FmJezWjpF5EJMXMjH3a7cM+7fbhH4f/gxdmvsDY6WO55p1ruDbvWg7b5TDO7ns2R/U4ito1a6c6XBGpIIVbC/lp5U/FetuL1peuXbptv1o1atGtaTd2bb4rw3OHB4l781xym+XSpG6TFF6BpBMl9SIiaaRerXqcsucpnLLnKcxZMYdHZzzKIzMeYeRzI2lWt9m2m2v3aLVHqkMVkTit3rg65nCZ75d/X2xmmeb1mpPbLJejehy1rcd91+a70qVJF2rWUMompVMLERFJU12bdOXGQ27kugHX8eacN3l4xsM8MPUB7vrkLvZuuzdHdj+S3q17s2erPenSuIvG4ItUsvWb17Nk7RIWFywucZm9bDa//Hf7WPcsy2KXpruwa/NdGdZ9GLs233VbAt+sXrMUXo1kOiX1IiJpLqtGFod3O5zDux3OL+t+4cmvnuTRLx7lhv/egOMANMhuwJ6t9qR3qyDJ7926N3u03IP62fVTHL1IZincWsiydctKTdSLllUbV+1Q3zBa1G9B65zWtM5pzT5N9mFgr4Hbxrt3bdJVz6iQpFBSLyKSQZrVa8bF/S/m4v4Xs3bTWr5Z9g1fLP6CL5YEy+NfPc7qqauBILnYpeku9G7Vu1iy36lRp4zu1d/qW1m0ZhHzVs1j7qq55K/Kp1m9ZvRt3ZfdW+6uhEl24O6s2rgqrkR92bplbPWtOxyjYe2G2xL13q17c3j9w7dtRy4t6rcoNlQmLy+PgQcOrMSrlepKSb2ISIaqn12ffu360a9dv21l7s7cVXP5YvEXfLnky23J/sSZE7f16jeq3Yg9W+25rWe/d+ve9GrZi3q16qXqUorZsGUD81bNC5L2lXO3Je9zVwXr+avyi83+EalWjVr0bNGTvm360rd1X/q07kOf1n1oWLthJV+FVJZNhZtYsHoB81fPJ391Pvmr8pm/ej7z18xn0ZpF25L1jYUbd6ibnZW9LRnv1LgT/dv1j5mot8pplTb/fYiUREm9iEgVYmZ0btyZzo07M3zX4dvKCzYV8PXSr7f16n+55EvGfzGeNZvWBPUwujfrvr1HP0z2OzTsUKG9+u7Oig0riiXr25L2sGzJ2iXFrwmjXcN2dGzUkf7t+jOy50g6NepEp8ad6NioIx0admDJ2iXMWDyD6YumM33xdCZ/P5lHZzy67Ri7NNmFPq370Ld1X/q2CZL9NjltMvovFtXB5sLNLFyzkPzVQaKevyp/+3qYwEe3Fwh+uLZv2J52DduR2zyX1vV3TNRb57SmcZ3GagNSZSipFxGpBnKyc9i3/b7s237fbWVbfSs/rfypWK/+tEXTeO7b57bt07hO4+1Jfpjw92rZi7q16sY8T+HWQhauWbhDoh7Z016wqaBYnTo169CxUUc6NepE71a9g/XGnejUKEja2zdsv9PH0jeq04gezXpwwu4nbCtbtGYR0xdPD5L9xdOZvmg6z898ftv7Leu33NabX5Tsd2vajRpWo0yfbapsKtzEojWLWLhmIQvWLGDhmoUsWrMIMyMnO6fY0iC7wQ5lOdk51M+un7Lr3bJ1C4vWLCrWu56/Or9YAr+4YPG2vzAVaVi7Ie0btqdDww70adUnWG/UYVtZ+4btaVC7QUquSSSVlNSLiFRTNazGtidMHrvbsdvK12xcw1dLvyqW7D88/WHWbl67rV6PZj3Ys9WedG7UmUUFi7Yl8PNXzy/2OHqAZnWb0bFRR3o068GhXQ/dlsAX9bS3qNciKb2lbRq0oU2DNgzrPmxb2eqNq/li8RfFkv07Prpj23Ce+rXq07t17yDJDxP+Xi17VerzAQq3FrJ07VIWrllYLGGPXBasWcDP637eoW7NGjVx9x2+g9LUr1V/xx8BtSN+BNSKUVbKD4b62fVxdxYVLCqxd33+6vksKli0w9j1nOycbYl5r269dkjWOzTqoKFUIiVQUi8iIsU0qN2A/Tvsz/4d9t9WttW3MmfFnCDJD4fwfLbgMybOnEjbBm3p2KgjB3Y8cFvvelFPe4dGHcjJzknh1RTXsHZDDup0EAd1Omhb2abCTXy77NttQ3dmLJ7B+C/Gc/9n9wNBotyzRc/tPfqt+9K7dW8a12lcpnO7O8vXL98hOY/eXlyweIdkt4bVoFX9Vts+633b70u7Bu1o26BtsaVZvWYYxqbCTRRsKqBgUwFrNq3Zth69rNkY8d7m7eXL1y9n7sq5xY6xZeuWuK+1htXY4Rrq1apHh4Yd6NCoA4fucmiwHpGsd2gYJOwaDiOSGCX1IiKyUzWsBt2adqNb024ct9tx28rdPeOTsOys7G031J7FWcD2HzGR4/Sn/DCF8V+M31avS+Muwfj8Vn3o26Yvy9Yuw3/0UnvYY92s2axus21Jea+WvWIm661yWpXp4UO1a9amds3aFTrvedEPhWI/BKJ/JIQ/ILb6Vto1aFesp13j10WSS0m9iIgkrKomaZE/Ykb0HLGtfHHB4m2J/owlwevEmRO3V5y6fbVBdoNtSfkBHQ+gbU7xRL1dw3a0zmlNnZp1KvHKEpedlU3Tuk1pWrdpqkMRkRiU1IuIiMSpdU5rhnQbwpBuQ7aVrdm4hi+XfMnrH7/OoH6DaNugLW1y2uhmTRGpVErqRUREyqFB7QYc0PEANs/ZzMDOA1MdjohUU5kxb5eIiIiIiJRISb2IiIiISIZTUi8iIiIikuGU1IuIiIiIZDhz953vVU2Y2RpgVqrjkLTSHNjxsY1S3aldSCxqFxKL2oXEkuvuFTpFlma/KW6Wu++d6iAkfZjZVLUJiaZ2IbGoXUgsahcSi5lN3fleZaPhNyIiIiIiGU5JvYiIiIhIhlNSX9yYVAcgaUdtQmJRu5BY1C4kFrULiaXC24VulBURERERyXDqqRcRERERyXBVPqk3s33MrNDMRkSUDTGzWWY228yuKqGemdk94T5fmtmvylJf0pOZnRJ+n1+a2Ydm1jviPbUL2en3qDZQPZhZBzN7x8xmmtk3ZnZpWN7UzN4ws+/D1yYl1I/ZFuKtL+nLzLLMbLqZvRxuq00IZtbYzCaY2Xfhvxv7VXrbcPcquwBZwNvAZGBERNkPQFcgG/gC6Bmj7jDgVcCAfYFPylJfS3ouwP5Ak3B9aFm/V7WLqr3E8z2qDVSPBWgD/CpcbwD8D+gJ/BW4Kiy/Cri9LO0onvpa0nsBrgCeBF6O9ztVm6j6CzAOOCdczwYaV3bbqOo99RcDzwNLI8r6AbPdfY67bwKeBobHqDscGO+Bj4HGZtamDPUlDbn7h+6+Itz8GGgfrqtdCMT3PaoNVAPuvsjdPw/X1wAzgXYE3+m4cLdxwDExqpfWFuKpL2nKzNoDRwAPRRSrTVRzZtYQOBgYC+Dum9x9JZXcNqpsUm9m7YBjgQej3moH5Edszw/LopW0X7z1Jf2dTdDjCmoXEojne1QbqGbMrDPQF/gEaOXuiyBI/IGWMaqU1hbiqS/p6y7gD8DWiDK1CekKLAMeCYdmPWRm9anktlFlk3qC//CudPfCqHKLsW+sKYBK2i/e+pLGzOwQgqT+yqKiGLupXVQ/8XyPagPViJnlEPzF9zJ3Xx1vtRhlagsZzsyOBJa6+7REqscoU5uoOmoCvwL+6e59gbUEw2XiUWFto2YildKVmV0InBtuNgKeNjOA5sAwM9tC8AuoQ0S19sDCGIcrab/sOOtLmohqF8MI2sNDwFB3/yUsV7sQiK8dqA1UE2ZWiyChf8LdJ4bFS8ysjbsvCoddLY1RtbR2FE99SU8HAEeb2TCgDtDQzB5HbUKC73e+u38Sbk8gSOortW1UqZ56d7/f3fuESxd37+zunQk+3Avc/T/AZ0B3M+tiZtnAKGBSjMNNAk4PZ7rYF1gV/ukj3vqSJiLbBcEP2YnAae7+v4jd1C4E4vse1QaqAQt6hMYCM939joi3JgFnhOtnAC/GqF5aW4invqQhd/+ju7cP84pRwNvufipqE9Weuy8G8s0sNywaDHxLZbeNyr47OBUL8Cjh7Dfh9jCCmQx+AP4UUX4+cH64bsD94T5fAXvvrL6W9F8IeuhXADPCZarahZaoNrLD96g2UP0W4ECCP4F/GfHvxTCgGfAW8H342jTcvy0weWdtoaT6WjJrAQayffYbtQktAH2AqeG/Gf8BmlR229ATZUVEREREMlyVGn4jIiIiIlIdKakXEREREclwSupFRERERDKcknoRERERkQynpF5EREREJMMpqRcRERERyXBK6kVEREREMpySehERyShm5ma21sxuqeDjvm1mG8zs/Yo8rohIZVBSLyJSCcysIGLZambrI7ZPSXV8iTKzn8zs1yk4dW93/1NULEPN7A0zuzNWBTM72cymhp/5IjN71cwOLHrf3QcRPD1YRCTjKKkXEakE7p5TtADzgKMiyp5IdXyxmFnNDDvHYcAIoEGM81wB3AXcCrQCOgIPAMMr8PwiIimjpF5EJA2YWVsze97MlpnZj2Z2ScR7P5nZ783sy3DYyVgzaxX2NK8xszfNrEnU/n80s2/NbIWZPWJmdcpwrivN7EtgrZnVNLOrzOyH8Fzfmtmx4b6PESTHL4W9338Iy93MukUc81Ezu7mUc3QsKZ4yehD4N/Bu1GfbCLgRuNDdJ7r7Wnff7O4vufvvEzyXiEhaUVIvIpJiZlYDeAn4AmgHDAYuM7PDI3Y7HjgU6AEcBbwKXA00J/i3PDoRPgU4HNglrHNNGc51EnAE0NjdtwA/AAcBjYAbgMfNrI27n0bxvzr8tQyXXXSOpsALO4knLu4+y91PcPfxUW/tB9QJzyMiUiUpqRcRSb19gBbufqO7b3L3OQQ9zqMi9rnX3Ze4+wLgPeATd5/u7hsJktW+Uce8z93z3X05cAtBEh3vue4J664HcPfn3H2hu29192eA74F+5bzme9w9H+gVRzzl1Qz4OfyBIiJSJSV9vKSIiOxUJ6Ctma2MKMsiSN6LLIlYXx9jOyfqmPkR63OBtmU4V2RdzOx04Aqgc1iUQ/AXgvIoOkc88ZTXL0BzM6upxF5Eqiol9SIiqZcP/Oju3SvwmB0i1jsCC8twLi9aMbNOBD3ng4GP3L3QzGYAFr1vhHVAvYjt1sD8Es6RjGuP9hGwATgGmJDE84iIpIyG34iIpN6nwOrw5tG6ZpZlZr3MbJ9yHPNCM2tvZk0Jxt4/k+C56hMk4MsAzOwsgiEzRZYAXaPqzABODo89BBhQSpzJuPZi3H0VcC1wv5kdY2b1zKxWOAVmWe4DEBFJW0rqRURSzN0LCW5+7QP8CPwMPERwY2qingSmAHPC5eZEzuXu3wL/IOjtXgLsAXwQsctfgGvMbKWZ/S4suzQ8x0qCG3b/U1KQSbr2WOe5g2AI0TUEP1DygYtKi01EJJOYe6y/nIqISKYys5+Ac9z9zVTHkgxmtgHYSHCz7Z8r8LhvAPsCn7r74Io6rohIZdCYehERySjuXmfneyV03EOTcVwRkcqg4TciIiIiIhlOw29ERERERDKceupFRERERDKcknoRERERkQynpF5EREREJMMpqRcRERERyXBK6kVEREREMpySehERERGRDKekXkREREQkwympFxERERHJcP8PJX90eY4V7KoAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import pyradi.ryplot as ryplot\n", "%matplotlib inline\n", "\n", "data=np.asarray([\n", "[ 50 , 82.78 ] ,\n", "[ 45 , 65.25 ] ,\n", "[ 40 , 50.98 ] ,\n", "[ 35 , 39.47 ] ,\n", "[ 30 , 30.26 ] ,\n", "[ 25 , 22.97 ] ,\n", "[ 20 , 17.24 ] ,\n", "[ 15 , 12.8 ] ,\n", "[ 10 , 9.38 ] ,\n", "[ 5 , 6.79 ] ,\n", "[ 0 , 4.85 ] ,\n", "[ -5 , 3.41 ] ,\n", "[ -10 , 2.36 ] ,\n", "[ -15 , 1.61 ] ,\n", "[ -20 , 1.08 ] ,\n", "[ -25 , 0.71 ] ])\n", "temperature = data[:,0]+273.15\n", "absh = ryutils.abshumidity(temperature).reshape(-1,1)\n", "data = np.hstack((data,absh))\n", "# absh = ryutils.abshumidity(temperature, 2).reshape(-1,1)\n", "# data = np.hstack((data,absh))\n", "data = np.hstack((data, 100 * np.reshape((data[:,1]-data[:,2])/data[:,2],(-1,1))))\n", "print(' deg C Testvalue Fn value % Error')\n", "print(data)\n", "\n", "p=ryplot.Plotter(1, 2, 1, figsize=(12,6))\n", "p.plot(1,data[:,0].T, data[:,2].T,'Absolute humidity vs temperature','Temperature [$^\\circ$C]','Absolute humidity g/m$^3$]')\n", "p.plot(2,data[:,0].T, data[:,3].T,'\\% humidity error vs temperature','Temperature [$^\\circ$C]','\\% error')\n", "# p.saveFig('abshum.eps')\n", "\n", "\n", "#highest ever recorded absolute humidity was at dew point of 34C\n", "print('Highest recorded absolute humidity was {0}, dew point {1} deg C'.\\\n", " format(ryutils.abshumidity(np.asarray([34 + 273.15]))[0],34))" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "## Transformation between polar and cartesian coordinates" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "The \n", "[`ryutils.cart2polar`](http://nelisw.github.io/pyradi-docs/_build/html/ryutils.html#pyradi.ryutils.cart2polar)\n", "and\n", "[`ryutils.polar2cart`](http://nelisw.github.io/pyradi-docs/_build/html/ryutils.html#pyradi.ryutils.polar2cart)\n", "provides conversion between polar and cartesian coordinates.\n", "\n", "The [`ryutils.cart2polar`](http://nelisw.github.io/pyradi-docs/_build/html/ryutils.html#pyradi.ryutils.cart2polar)\n", "function signature:\n", "\n", " `cart2polar(x, y)`\n", "\n", "- `x (float np.array)` x values in array format.\n", "- `y (float np.array)` y values in array format.\n", "\n", "returns\n", "\n", "- `r (float np.array)` radial component for given (x,y).\n", "- `theta (float np.array)` angular component for given (x,y).\n", "\n", "\n", "The [`ryutils.polar2cart`](http://nelisw.github.io/pyradi-docs/_build/html/ryutils.html#pyradi.ryutils.polar2cart)\n", "function signature:\n", "\n", " `polar2cart(r, theta)`\n", "\n", "- `r (float np.array)` radial values in array format.\n", "- `theta (float np.array)` angular values in array format.\n", "\n", "returns\n", "\n", "- `x (float np.array)` x component for given (r, theta).\n", "- `y (float np.array)` y component for given (r, theta).\n", "\n", "This transformation code was written by [Joe Kington](https://stackoverflow.com/questions/3798333/image-information-along-a-polar-coordinate-system)." ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "run_control": { "frozen": false, "read_only": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "5.000000000000001 5.000000000000001 -> 7.0710678118654755 0.7853981633974483\n", "[3. 5.] [5. 3.] -> [5.83095189 5.83095189] [1.03037683 0.5404195 ]\n", "[[3. 5.]\n", " [5. 3.]] [[5. 3.]\n", " [3. 5.]] -> [[5.83095189 5.83095189]\n", " [5.83095189 5.83095189]] [[1.03037683 0.5404195 ]\n", " [0.5404195 1.03037683]]\n" ] } ], "source": [ "import pyradi.ryutils as ryutils\n", "import numpy as np\n", "\n", "r, t = ryutils.cart2polar(5, 5)\n", "x, y = ryutils.polar2cart(r, t)\n", "print('{} {} -> {} {}'.format(x, y, r, t))\n", "\n", "r, t = ryutils.cart2polar(np.asarray([3, 5]), np.asarray([5, 3]))\n", "x, y = ryutils.polar2cart(r, t)\n", "print('{} {} -> {} {}'.format(x, y, r, t))\n", "\n", "r, t = ryutils.cart2polar(np.asarray([[3, 5], [5, 3]]), np.asarray([[5, 3], [3, 5]]))\n", "x, y = ryutils.polar2cart(r, t)\n", "print('{} {} -> {} {}'.format(x, y, r, t))" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "## Calculating coordinate index arrays for a Numpy array" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "The \n", "[`ryutils.index_coords`](http://nelisw.github.io/pyradi-docs/_build/html/ryutils.html#pyradi.ryutils.index_coords) function calculates zero-based (x,y) coordinate arrrays for a numpy array indices, relative to some origin.\n", "\n", "This function calculates two meshgrid arrays containing the coordinates of the \n", "input array. The origin of the new coordinate system defaults to the \n", "center of the image, unless the user supplies a new origin. \n", "\n", "The data format can be data.shape = (rows, cols, frames) or \n", "data.shape = (frames, rows, cols), the format of which is indicated by the \n", "framesFirst parameter.\n", "\n", "The function signature is:\n", "\n", " `index_coords(data, origin=None, framesFirst=True)`\n", "\n", "- `data (np.array)` array for which coordinates must be calculated.\n", "- `origin ( (x-orig, y-orig) )` data-coordinates of where origin should be\n", "- `framesFirst (bool)` True if data.shape is (frames, rows, cols), False if data.shape is (rows, cols, frames)\n", "\n", "returns\n", "\n", "- `x (float np.array)` x coordinates in array format.\n", "- `y (float np.array)` y coordinates in array format.\n", "\n", "This transformation code was originally written by [Joe Kington](https://stackoverflow.com/questions/3798333/image-information-along-a-polar-coordinate-system)." ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "The following example demonstrates the variations obtained for frames first or frames last, and the location of the origin." ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "run_control": { "frozen": false, "read_only": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "framesFirst=True, origin=None\n", "x=[[-2 -1 0 1]\n", " [-2 -1 0 1]\n", " [-2 -1 0 1]]\n", "y=[[-1 -1 -1 -1]\n", " [ 0 0 0 0]\n", " [ 1 1 1 1]]\n", "\n", "framesFirst=False, origin=None\n", "x=[[-1 0 1]\n", " [-1 0 1]]\n", "y=[[-1 -1 -1]\n", " [ 0 0 0]]\n", "\n", "framesFirst=True, origin=(0,0)\n", "x=[[0 1 2 3]\n", " [0 1 2 3]\n", " [0 1 2 3]]\n", "y=[[0 0 0 0]\n", " [1 1 1 1]\n", " [2 2 2 2]]\n", "\n", "framesFirst=True, origin=(3,2)\n", "x=[[-3 -2 -1 0]\n", " [-3 -2 -1 0]\n", " [-3 -2 -1 0]]\n", "y=[[-2 -2 -2 -2]\n", " [-1 -1 -1 -1]\n", " [ 0 0 0 0]]\n" ] } ], "source": [ "import pyradi.ryutils as ryutils\n", "import numpy as np\n", "\n", "data = np.ones([2, 3, 4])\n", "\n", "x, y = ryutils.index_coords(data, origin=None, framesFirst=True)\n", "print('framesFirst=True, origin=None')\n", "print('x={}\\ny={}'.format(x,y))\n", " \n", "x, y = ryutils.index_coords(data, origin=None, framesFirst=False)\n", "print('\\nframesFirst=False, origin=None')\n", "print('x={}\\ny={}'.format(x,y))\n", " \n", "x, y = ryutils.index_coords(data, origin=(0,0), framesFirst=True)\n", "print('\\nframesFirst=True, origin=(0,0)')\n", "print('x={}\\ny={}'.format(x,y))\n", "\n", "x, y = ryutils.index_coords(data, origin=(3,2), framesFirst=True)\n", "print('\\nframesFirst=True, origin=(3,2)')\n", "print('x={}\\ny={}'.format(x,y))\n", " " ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "## Latex utilities" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "### LaTeX code for EPS graphic inclusion " ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "The [`ryfiles.epsLaTexFigure`](http://nelisw.github.io/pyradi-docs/_build/html/ryfiles.html#pyradi.ryfiles.epsLaTexFigure)\n", "function writes the LaTeX code to include an eps graphic as a latex figure. The text is added to an existing file.\n", "\n", "The function signature is:\n", " `epsLaTexFigure(filename, epsname, caption, scale, filemode='a')`\n", "\n", "- `fname (string)` text writing output path and filename.\n", "- `epsname (string)` filename/path to eps file (relative to where the LaTeX document is built).\n", "- `caption (string)` figure caption.\n", "- `scale (double)` figure scale to textwidth [0..1].\n", "- `filemode (string)` file open mode (a=append, w=new file) (optional)." ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "run_control": { "frozen": false, "read_only": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\\begin{figure}[tb]\n", "\\centering\n", "\\includegraphics[width=0.75\\textwidth,height=w\\textheight,keepaspectratio]{picture.eps}\n", "\\caption{This is the caption. \\label{fig:picture.eps}}\n", "\\end{figure}\n", "\n", "\n", "\\begin{figure}[tb]\n", "\\centering\n", "\\includegraphics[width=0.75\\textwidth,height=w\\textheight,keepaspectratio]{picture.eps}\n", "\\caption{This is the caption. \\label{fig:picture.eps}}\n", "\\end{figure}\n", "\n", "\n", "\\begin{figure}[tb]\n", "\\centering\n", "\\includegraphics[width=0.75\\textwidth,height=w\\textheight,keepaspectratio]{picture.eps}\n", "\\caption{This is the caption. \\label{fig:picture.eps}}\n", "\\end{figure}\n", "\n", "\n", "\\begin{figure}[tb]\n", "\\centering\n", "\\includegraphics[width=0.75\\textwidth,height=w\\textheight,keepaspectratio]{picture.eps}\n", "\\caption{This is the caption. \\label{fig:picture.eps}}\n", "\\end{figure}\n", "\n", "\n" ] } ], "source": [ "ryfiles.epsLaTexFigure('eps.tex', 'picture.eps', 'This is the caption', 0.75,'w')\n", "\n", "import os.path \n", "if os.path.sep == '/':\n", " !cat eps.tex\n", "else:\n", " !type eps.tex" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "### Numpy arrays in LaTeX code" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "The [`ryfiles.arrayToLaTex`](http://nelisw.github.io/pyradi-docs/_build/html/ryfiles.html#pyradi.ryfiles.arrayToLaTex)\n", "function writes a numpy array to latex table format in output file.\n", "The table can contain only the array data (no top header or\n", "left column side-header), or you can add either or both of the\n", "top row or side column headers. Leave 'header' or 'leftcol' as\n", "None is you don't want these.\n", "\n", "The output format of the array data can be specified, i.e.\n", "scientific notation or fixed decimal point.\n", "\n", "The function signature is:\n", " `arrayToLaTex(filename, arr, header=None, leftCol=None,formatstring='%1.4e', filemode='wt')`\n", " \n", "- `fname (string)` text writing output path and filename.\n", "- `arr (np.array[N,M])` array with table data.\n", "- `header (string)` column header in final latex format (optional).\n", "- `leftCol ([string])` left column each row, in final latex format (optional).\n", "- `formatstring (string)` output format precision for array data (see [numpy.savetxt](http://docs.scipy.org/doc/numpy/reference/generated/numpy.savetxt.html)) (optional).\n", "- `filemode (string)` file open mode (a=append, w=new file) (optional). \n", "\n", "With `\\usepackage{siunitx}` you can use the siunitx `\\num` formatter to format the output." ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "run_control": { "frozen": false, "read_only": false } }, "outputs": [ { "ename": "TypeError", "evalue": "write() argument must be str, not bytes", "output_type": "error", "traceback": [ "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[1;31mTypeError\u001b[0m Traceback (most recent call last)", "\u001b[1;32m~\\AppData\\Local\\Temp/ipykernel_19740/2572545408.py\u001b[0m in \u001b[0;36m\u001b[1;34m\u001b[0m\n\u001b[0;32m 1\u001b[0m \u001b[0marr\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0masarray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1.0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m3\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m4\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m5\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m6\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m7\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m8\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m9\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 2\u001b[0m \u001b[0mryfiles\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0marrayToLaTex\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'array.tex'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0marr\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 3\u001b[1;33m \u001b[0mryfiles\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0marrayToLaTex\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'array.tex'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0marr\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mformatstring\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'%.1f'\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mfilemode\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'a'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 4\u001b[0m \u001b[0mheaderonly\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;34m'Col1 & Col2 & Col3'\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 5\u001b[0m \u001b[0mryfiles\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0marrayToLaTex\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'array.tex'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0marr\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mheaderonly\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mformatstring\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'%.3f'\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mfilemode\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'a'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", "\u001b[1;32mK:\\WorkN\\pyradi\\pyradi\\pyradi\\ryfiles.py\u001b[0m in \u001b[0;36marrayToLaTex\u001b[1;34m(filename, arr, header, leftCol, formatstring, filemode)\u001b[0m\n\u001b[0;32m 892\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 893\u001b[0m \u001b[0mfile\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mopen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfilename\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfilemode\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 894\u001b[1;33m \u001b[0mfile\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mwrite\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'\\\\begin{{tabular}}{{ {0} }}\\n\\hline\\n'\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'|'\u001b[0m\u001b[1;33m+\u001b[0m \u001b[0mnumcols\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;34m'c|'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mencode\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'utf-8'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 895\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 896\u001b[0m \u001b[1;31m#write the header\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", "\u001b[1;31mTypeError\u001b[0m: write() argument must be str, not bytes" ] } ], "source": [ "arr = np.asarray([[1.0,2,3],[4,5,6],[7,8,9]])\n", "ryfiles.arrayToLaTex('array.tex', arr)\n", "ryfiles.arrayToLaTex('array.tex', arr, formatstring='%.1f',filemode='a')\n", "headeronly = 'Col1 & Col2 & Col3'\n", "ryfiles.arrayToLaTex('array.tex', arr, headeronly, formatstring='%.3f',filemode='a')\n", "header = 'Col 1 & Col 2 & Col 3'\n", "leftcol = ['XX','Row 1','Row 2','Row 3']\n", "ryfiles.arrayToLaTex('array.tex', arr, header, leftcol, formatstring=r'\\num{%.6e}',filemode='a')\n", "\n", "import os.path \n", "if os.path.sep == '/':\n", " !cat array.tex\n", "else:\n", " !type array.tex" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "### Upright micro symbol in LaTeX" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "The [`ryutils.upMu`](http://nelisw.github.io/pyradi-docs/_build/html/ryutils.html#pyradi.ryutils.upMu)\n", "function returns an upright $\\mu$ symbol, as required in the SI system for the $10^{-6}$ prefix. \n", "\n", "The upright symbol requires that the siunitx LaTeX package be installed on the computer running the code. This function also changes the Matplotlib rcParams file. This symbol can be used in graphs or text.\n", "\n", "The function signature is: \n", " `upMu(uprightMu=True, textcomp=False)`\n", "\n", "- `uprightMu (bool)` signals upright (True, default) or regular (False) symbol (optional).\n", "- `textcomp (bool)` if True use the textcomp package, else use siunitx package (optional).\n", "\n", "The [`ryutils.upMu`](http://nelisw.github.io/pyradi-docs/_build/html/ryutils.html#pyradi.ryutils.upMu) function returns a string with LaTeX code for the micro symbol.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "run_control": { "frozen": false, "read_only": false } }, "outputs": [], "source": [ "import pyradi.ryutils as ryutils\n", "print('Using siunitx package:')\n", "print('{} renders in LaTeX as an upright symbol'.format(ryutils.upMu()))\n", "print('{} renders in LaTeX as an upright symbol'.format(ryutils.upMu(True)))\n", "print('{} renders in LaTeX as an italic/slanted symbol'.format(ryutils.upMu(False)))\n", "print('\\nUsing textcomp package:')\n", "print('{} renders in LaTeX as an upright symbol'.format(ryutils.upMu(textcomp=True)))\n", "print('{} renders in LaTeX as an upright symbol'.format(ryutils.upMu(True,textcomp=True)))\n", "print('{} renders in LaTeX as an italic/slanted symbol'.format(ryutils.upMu(False,textcomp=True)))\n" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "### Create clean filenames" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "The [`ryfiles.cleanFilename`](http://nelisw.github.io/pyradi-docs/_build/html/ryfiles.html#pyradi.ryfiles.cleanFilename)\n", "function cleans a string by removing selected characters. It creates a legal and 'clean' source string from a string by removing some \n", "clutter and characters not allowed in filenames. A default set is given but the user can override the default string.\n", "\n", "The function signature is:\n", " `cleanFilename(sourcestring, removestring =\" %:/,.\\\\[]\")`\n", " \n", "- `sourcestring (string)` the string to be cleaned.\n", "- `removestring (string)` remove all these characters from the string (optional).\n", " \n", "It is quite useful to use a figure caption as filename, first cleaning it up. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "run_control": { "frozen": false, "read_only": false } }, "outputs": [], "source": [ "print ('\\nTest CleanFilename function:')\n", "inString=\"aa bb%cc:dd/ee,ff.gg\\\\hh[ii]jj\"\n", "print('{0}\\n{1}\\n'.format(inString, ryfiles.cleanFilename(inString) ))\n", "inString=\"aa bb%cc:dd/ee,ff.gg\\\\hh[ii]jj\"\n", "print('{0}\\n{1}'.format(inString, ryfiles.cleanFilename(inString, \"abcd\") ))" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "## Reading two-dimensional lookup tables" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "The [`ryfiles.read2DLookupTable`](http://nelisw.github.io/pyradi-docs/_build/html/ryfiles.html#pyradi.ryfiles.read2DLookupTable)\n", "Reads a 2D lookup table and extract the data. The table has the following format:\n", "\n", " line 1: xlabel ylabel title\n", " line 2: 0 (vector of y (col) abscissa)\n", " lines 3 and following: (element of x (row) abscissa), followed\n", " by table data.\n", "\n", "From line/row 3 onwards the first element is the x abscissa value\n", "followed by the row of data, one point for each y abscissa value.\n", "\n", "The file format can depicted as follows:\n", "\n", " x-name y-name ordinates-name\n", " 0 y1 y2 y3 y4\n", " x1 v11 v12 v13 v14\n", " x2 v21 v22 v23 v24\n", " x3 v31 v32 v33 v34\n", " x4 v41 v42 v43 v44\n", " x5 v51 v52 v53 v54\n", " x6 v61 v62 v63 v64\n", "\n", "The function has the signature:\n", "\n", " `read2DLookupTable(filename)`\n", " \n", "- `fname (string)` input path and filename \n", "\n", "The function returns the following:\n", "\n", "- `xVec ((np.array[N]))` x abscissae\n", "- `yVec ((np.array[M]))` y abscissae\n", "- `data ((np.array[N,M]))` data corresponding the x,y\n", "- `xlabel (string)` x abscissa label\n", "- `ylabel (string)` y abscissa label\n", "- `title (string)` dataset title\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "run_control": { "frozen": false, "read_only": false } }, "outputs": [], "source": [ "filecont = ['x-name y-name ordinates-name\\n',\n", "'0 1 2 3 4\\n',\n", "'10 11 12 13 14\\n',\n", "'20 21 22 23 24\\n',\n", "'30 31 32 33 34\\n',\n", "'40 41 42 43 44\\n',\n", "'50 51 52 53 54\\n',\n", "'60 61 62 63 64\\n']\n", "\n", "with open('2dfile.dat', 'wt') as f:\n", " f.writelines(filecont)\n", "\n", "import os.path \n", "if os.path.sep == '/':\n", " !cat 2dfile.dat\n", "else:\n", " !type 2dfile.dat \n", " \n", "xVec,yVec,data,xlabel, ylabel, title = ryfiles.read2DLookupTable('2dfile.dat')\n", "print('\\n')\n", "print('xVec={}'.format(xVec))\n", "print('yVec={}'.format(yVec))\n", "print('data={}'.format(data))\n", "print('xlabel={}'.format(xlabel))\n", "print('ylabel={}'.format(ylabel))\n", "print('title={}'.format(title))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Probability Density by Kernel Density Estimate\n", "\n", "Jake vdPlas makes a good point why and how histograms are not really good probability density estimators [here](https://jakevdp.github.io/PythonDataScienceHandbook/05.13-kernel-density-estimation.html#Motivating-KDE:-Histograms). Depending on how the data points fall into the bins strange effects may arise. He writes 'The problem with our two binnings stems from the fact that the height of the block stack often reflects not on the actual density of points nearby, but on coincidences of how the bins align with the data points. This mis-alignment between points and their blocks is a potential cause of the poor histogram results.\n", "\n", "The histogram approach is to have (normally same-size) fixed bin edges, in to which the data points are falling and being counted. The Kernel Density method works the other way around make fixed-width block around each data point and add up the blocks' heights. This way the weights add up around the data points' locations. This is somewhat making a bin centred around each data point and then adding up the bin counts at a fine resolution.\n", "\n", "Working with rectangular blocks the new density function has rough edges at the blocks' corners. The next logical step is to use smooth functions (called kernels) instead of rectangular blocks. This yields smooth functions with more appropriate density distributions than does the traditional histogram. In this approach you must define/select the kernel shape and width. If the kernel is too wide a kernel's influence can be too wide, masking the data points. If the kernel is too narrow it will leave 'holes' between the data points. Just image a dirac delta kernel function which simply repeats the data set.\n", "\n", "scikit-learn provides an efficient library for kernel density estimation. The library provides six different kernel shapes and a method whereby an appropriate kernel width can be found. The `pyradi.ryprob.KDEbounded` function packages the scikit learn functionality in an easily usable single function.\n", "\n", "The gaussian kernels typically used in KDE reaches to infinity and the density function domain therefore spans infinity in x.\n", "Suppose the data set is known to be bounded by zero on the lower end: non-zero values are not possible.\n", "The KDE algorithm will happily return a nonzero density for $x<0$, so some means must be found to eliminate density outside the domain bounds.\n", "\n", "\n", "To correct this non-zero density outside of the valid domain:\n", "\n", "- Make additional datasets mirrored around the boundary or boundaries (if both sides are bounded).\n", "- Calculated the KDE of each set separately: the base set, the lower set and the upper set.\n", "- Add the result of the sets together.\n", "- Select the results in the valid domain and set the values outside the valid domain to zero.\n", "\n", "This procedure ensures that the PDF derivative at the boundary is zero and that the density 'outside' the domain is accounted for by mirroring the density back into the domain.\n", "\n", "The next challenge is finding the correct bandwidth. From JvdP book:\n", " \n", "''The choice of bandwidth within KDE is extremely important to finding a suitable density estimate, and is the knob that controls the bias–variance trade-off in the estimate of density: too narrow a bandwidth leads to a high-variance estimate (i.e., over-fitting), where the presence or absence of a single point makes a large difference. Too wide a bandwidth leads to a high-bias estimate (i.e., under-fitting) where the structure in the data is washed out by the wide kernel.\n", "\n", "There is a long history in statistics of methods to quickly estimate the best bandwidth based on rather stringent assumptions about the data: if you look up the KDE implementations in the SciPy and StatsModels packages, for example, you will see implementations based on some of these rules.\n", "\n", "In machine learning contexts, we've seen that such hyperparameter tuning often is done empirically via a cross-validation approach. With this in mind, the KernelDensity estimator in Scikit-Learn is designed such that it can be used directly within the Scikit-Learn's standard grid search tools. Here we will use GridSearchCV to optimize the bandwidth for the preceding dataset. Because we are looking at such a small dataset, we will use leave-one-out cross-validation, which minimizes the reduction in training set size for each cross-validation trial''\n", "\n", "We will not discuss the matter much deeper here, oher than to show how the kernel with optimisation can be done. In the `KDEbounded` function below, if `bandwidth=np.nan` the bandwidth optimisation will be done. If the `bandwidth` is assigned a numerical value, that value will be used.\n", "\n", "The Leave-One-Out algorithm used here is quite slow for large data sets. So use a smaller data set to determine the best bandwidth and then use the value numerically for the large runs.\n", "\n", "The `pyradi.ryprob.KDEbounded` function requires a scikit learn version higher than 0.19. This work was done with version 0.20.1\n", "\n", "A small writeup of this work with more examples is given in \n", "https://github.com/NelisW/PythonNotesToSelf/blob/master/KernelDensityEstimation/KernelDensityEstimation.ipynb\n", "\n", "https://kdepy.readthedocs.io/en/latest/examples.html#boundary-correction-using-mirroring\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "https://jakevdp.github.io/PythonDataScienceHandbook/05.13-kernel-density-estimation.html \n", "https://jakevdp.github.io/blog/2013/12/01/kernel-density-estimation/ \n", "https://mglerner.github.io/posts/histograms-and-kernel-density-estimation-kde-2.html?p=28 \n", "\n", "https://github.com/tommyod/KDEpy \n", "https://kdepy.readthedocs.io/en/latest/examples.html \n", "\n", "https://towardsdatascience.com/histograms-and-density-plots-in-python-f6bda88f5ac0 \n", "https://towardsdatascience.com/how-to-find-probability-from-probability-density-plots-7c392b218bab \n", "https://github.com/admond1994/calculate-probability-from-probability-density-plots/blob/master/cal_probability.ipynb \n", "\n", "https://stats.stackexchange.com/questions/405357/how-to-choose-the-bandwidth-of-a-kde-in-python \n", "https://www.homeworkhelponline.net/blog/math/tutorial-kde \n", "\n", "https://scikit-learn.org/stable/auto_examples/neighbors/plot_kde_1d.html \n", "https://scikit-learn.org/stable/modules/cross_validation.html#cross-validation \n", "https://scikit-learn.org/stable/modules/density.html \n", "https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.html \n", "https://scikit-learn.org/stable/modules/generated/sklearn.neighbors.KernelDensity.html \n", " \n", "https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gaussian_kde.html \n", "https://glowingpython.blogspot.co.za/2012/08/kernel-density-estimation-with-scipy.html \n", "https://pythonhosted.org/PyQt-Fit/KDE_tut.html \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `pyradi.ryprob.KDEbounded` function parameters are as follows:\n", "\n", "- x_d (np.array[N,]): domain over which the probability density values must be returned. Presumably this domain must cover the input data set closely, but it can be wider or otherwise different.\n", "- x (np.array[N,]): the input sample data set.\n", "- bandwidth (float): the bandwidth width to be used, np.nan if to be calculated. If bandwidth is np.nan, calculate the optimal kernel width, aka bandwidth. Be careful, this can take a while, especially for large data sets.\n", "- lowbnd (float): lower mirror fold boundary, np.nan means no lower bound and mirror\n", "- uppbnd (float): upper mirror fold boundary, np.nan means no upper bound and mirror\n", "- kernel (str): kernel to be used ['gaussian'|'tophat'|'epanechnikov'|'exponential'|'linear'|'cosine']\n", "\n", " \n", " " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pyradi.ryprob as ryprob\n", "import pyradi.ryplot as ryplot\n", "\n", "def make_data(N, f=0.3, rseed=1):\n", " rand = np.random.RandomState(rseed)\n", " x = rand.randn(N) #mean-0 var=1\n", " # create bimodal data by adding const to entries with index > f*N\n", " # this means that some of the date keep old values and some shift \n", " x[int(f * N):] += 5\n", " return x\n", "\n", "\n", "# helper function to set up experiments\n", "def showexperiment(bounds=(np.nan,np.nan),bandwidth=np.nan,scale=1,offset=0,kernel='gaussian'):\n", " \n", " x = scale * make_data(20) + offset\n", " \n", " if not np.isnan(bounds[0]) and not np.isnan(bounds[1]):\n", " lowbnd = bounds[0]\n", " uppbnd = bounds[1]\n", " x_d = np.linspace(lowbnd,uppbnd, 1000)\n", " else:\n", " lowbnd = np.nan\n", " uppbnd = np.nan\n", " tenp = (np.max(x) - np.min(x)) / 10\n", " x_d = np.linspace(np.min(x)-tenp,np.max(x)+tenp, 1000)\n", "\n", " x_d,kde,bandwidth,kernel = ryprob.KDEbounded(x_d,x,bandwidth=bandwidth,lowbnd=lowbnd,uppbnd=uppbnd, kernel=kernel)\n", "\n", " print(f'x.min={np.min(x):.4e} x.max={np.max(x):.4e} x.shape={x.shape}')\n", " print(f'x.min={np.min(x_d):.4e} x.max={np.max(x_d):.4e} x.shape={x_d.shape}')\n", " print(f'Bandwidth={bandwidth}')\n", " print(f'kernel={kernel}')\n", " print(f'Integral of PDF={np.trapz(kde,x_d)}')\n", "\n", " p = ryplot.Plotter(1,1,1)\n", " p.plot(1,x,np.full_like(x, 0.011*np.max(kde)),linestyle='',markers='^')\n", " p.plot(1,x_d,kde,'Kernel Density Estimate','$x_d$','$p$')\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "showexperiment(bounds=[-3.301,7.74481], bandwidth=np.nan, scale=1,offset=0) " ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "## ToDo" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "`saveHeaderArrayTextFile', 'loadColumnTextFile', 'loadHeaderTextFile`" ] }, { "cell_type": "markdown", "metadata": { "run_control": { "frozen": false, "read_only": false } }, "source": [ "## Python and module versions, and dates" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "run_control": { "frozen": false, "read_only": false } }, "outputs": [], "source": [ "try:\n", " import pyradi.ryutils as ryutils\n", " print(ryutils.VersionInformation('matplotlib,numpy,pyradi,scipy,pandas'))\n", "except:\n", " print(\"pyradi.ryutils not found\")\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "run_control": { "frozen": false, "read_only": false } }, "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.8.3" } }, "nbformat": 4, "nbformat_minor": 4 }