diff --git a/examples/03_britter_mcquaid_plume_models.ipynb b/examples/03_britter_mcquaid_plume_models.ipynb index cb5525e..df4fa01 100644 --- a/examples/03_britter_mcquaid_plume_models.ipynb +++ b/examples/03_britter_mcquaid_plume_models.ipynb @@ -1,22 +1,251 @@ { "cells": [ { + "attachments": { + "3da40a83-6628-40f3-a58b-03478e3d3cfe.png": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgwAAAIpCAIAAABBhffdAAAgAElEQVR4Aex9P56jOrM2/n437uzNcC8AHN0Id9I90XR0M9hAT3RPBF4A9gKMo5Od3gBsYCbqnsR4AeAFgKN7st6Av9/1c9969fJHxhiwkEUwo8ZSqfRIqKSqUmlyPB419SgEFAIKAYWAQqAKgf9X9VK9UwgoBBQCCgGFwP8ioISEGgcKAYWAQkAhUIuAEhK10KgfFAIKAYWAQkAJCTUGFAIKAYWAQqAWASUkaqFRPygEFAIKAYWAEhJqDCgEFAIKAYVALQJKSNRCo35QCCgEFAIKASUkZBgDjuPkeS5DS1QbFAIKAcEQmKjDdIL1iGJHIaAQUAgIhIDaSQjUGYoVhYBCQCEgGgJKSIjWI4ofhYBCQCEgEAJKSAjUGYoVhYBCQCEgGgKSCInlcjmfz6fTaRRFokGs+FEIKAQUAuNFQAbDdRRFjuO4rqvr+mKxUKb48Q5HxblCQCEgGgIjEBJRFP3jH/94eXnRNO3z8/OPP/7Y7/e+7y+XyzKak8kIWlRmW71RCCgEFAJiIiC6ugm7hD/++EPTtCiKvn379vDwYNv2arUqC4n5fC4myoorhYBCQCEwUgSEXnfHcfz09ARkj8cjrA5hGGqahp+22y0JhjzPF4tFFEVJkpimOdL+UGwrBBQCCgGhEBBaSGiatlwuoyja7/fH47GgSppMJmEY2rbNAuo4jmEY5U0Gm0elFQIKAYWAQqAhAv/RMN+tsi2Xy/3pKTNgGAa9dBwH6TiOC2KD8qiEQkAhoBBQCFyKgOhCgt8eCljk+/5qtdI0LQgCCAkYM/TTs9vtiI5lWdPplP4cIKHrOunErqkuz/Ovry/TNNHqr6+v/X5fSbCrGiuJV74s1JimqWmacRw/nB4W/8riF720bTvPc9C0LKtb4uAkjmN0WRRFwwwYYEUdOvxCp66ZYOzr64sA0TQNn9XX19fDw4Omae26oK7GiwYDmxmjDp/Jw8MDcdV5Rag0TVPqr8I8U/gcWCa7Tadp+vX1BZqdTDKV7ImubtI0zXGcKIoq1U22bcNEUW5bdHrKv8ZxfDgcyvk7eUPzl2EYNIDyPIdYQo+26Etsp8ptqeOZ5tC6DO3eF2ZkfAlxHIPa4XCgubUd/bOlMGUfTg+bufNZQNd1TdMwTnRdj6II0g4TImYEpFvM5phc+AWjKNJ1vTxQgTmNKBaEypdsBn669XfBgl/JQ+VLmBXLDeQzKfiv6FOs4UhEdctzYVTYto0PUNf1r9NzDaS10+lR+AfQ4/TDx8cH8atpWhiG9GchEYah7/uFlyP90/d9TdNGynyHbPN7vMOK+iZl23bfVQxDfxQNSZLEPz3s7NEtPuHp6ZbmwNQ0TavDZ0zqJsuy/vjjj58/f9LJaqz4upXVAlIzTdN1XQEZUyy1QCBN0xalVJHWCJinp3XxOyloWRbOopXbK/o5CShAYaOGvuXp6Wk+nzuO4/t+C9VNGQLx39i2HQSB+HwqDpsgoPyzm6A0ujxjDwjEsdSOQEgEQYDF13Q6TdM0CIL5fL7dbs/6uY692+g7SdPU8zz6UyXGjgA5XKiGjB0B4l9ircYIhAR1AxJYVjfZQ5DpuEBhdH9GUbTZbEbHtmK4DgHySKnLMJb30jTkesCvsRhfX3uvFMYnJJrDwR6kaF5K5VQI9I2ANMsXaRrSd4+Pmr7MQkKN4FEPTcW8QkAhIAICMgsJEfDthAe+T30nVSgiCgGFgEKgEgElJCphEeulNBZ4sWC9HTfSKEKlacjtxsIIah7TOYkRwNkPi6ZpSuw70Q9mQlNFNAuhWWzGnDrzAZwsy2oG2ChzKSExgm6zT88IGO2ZxSRJ5Jhee4rZ0DP8inwtAtPTU/vzyH9Q6qZxdKDSOGmaZpom58jPODryxKUyMo2os+6H1bqTZ0pIjGAMIKLtCBhVLDZDgKIiNsuucikEhkCgLhaAEhJDoH9lHUrzeyWAohWX+OCVaFArfpojUOeGoIREcwxvlrNOwt+MoQYVSxN5okFbL86SZdnFZVQBhUDPCNQdLDsjJNQatud+aUR+jCpsZZvldK0clhVOA9VPY0SgbrY/IyTGuIZF9+R5jmsYxthbEvA8RsEmAeyqCQqBzhE4IyQ6r28wgjIt1uok/GBgqooUAgqBu0VAWiEhU4+apqkW5jJ1qGqLQmBECMgsJD4/P8s9MdJVefMLrstNVm8UAgoBhQAfgTrXJk3TZD5x/fv378lkUobGMIyRiopyW9QbhYBCQCHQHIHo9JTzx3H88+fP8nvJhYRt2+UFeBzHowvtEEVRnufqcrrKEaxeKgQUAs0R0HW9UncdRVGdm5LMO4nKDVSTK+2aIz5MzjRNV6uVEhLDoK1qUQhIjEDdBFg5WwIHmW0SKt6RxGNdNU0hoBDoEAGOfkVmIVG3e+oQWUVKIaAQUAhIgADnzIDMQoKzgZKgU1UTFAIKAYXAAAjILCRWq9UACKoqFAIKAYWAxAjILCQk7jbVNIWAQkAh0CECnIicSkh0iLMipRBQCCgERonAndokRtlXVUybpum6btUv6p1CQCGgEOgAAc5OQuZzEh0gJwYJdce1GP2guFAISIuA2klI27WqYQoBhYBCoFcElE2iV3i7IR5FUd05yW4qUFQUAgoBhUANAkpI1AAj0us0TdVFbyJ1iOJFIXBHCCghcUedrZqqEFAIKAQuRUAJiUsRu01+XddvU7GqVSGgELhvBKQVEpZlydSzh8NBpuaotigEFAJjQUBaIfH19TWWPlB8KgQUAgoBYRGQVkjs93thQVeMKQQUAgqBsSAgrZCQSd1kmqaySYzli1J8KgRGikDdBTzSnrjmHCAcXRfatq2ExOh6TTGsEJADAWl3EnJ0D1qxXC6fnp5kapFqi0JAISAaAmmaVrKkhEQlLOqlQkAhoBC4LwSWy2Vlg5WQqIRFvVQIKAQUAnIiYJpmpfmhbichrU1Czu5Vrboagbrlkm3b6lL0q9FVBEaAwD/+8Q/btsuM1o1/JSTKWKk3o0cAC6U8zxHzKk3Tsy7R+/0+DEO0PDo9mqZZljWdTvHv6EFRDVAInBD4888/K5FI07RSTighUQmXejkmBKIoYldGURQ5jlNoAE30lChnoDdxHEPMsLtyXdfn87lhGPBIVnF5CS6VGBcClZJA07S695ILiclkUug/zBHqGp8CLKP7M4qi+PRgr7Berz3PQyssy6L+xcx+aeuC00P7D4gKkhygBsqWZc1Pz6VVqPwKgeERyPP8cDg8nJ5CYOk8zwuLrX+xd5T30TTN9/0wDNFE//QYhqFpGr0cS+tHx3BPwCZJwm4aaBwnSdJTjSzZJEkwiqheJFTvsCiptJgI+L5fGLeFP33fr+R8cjweC1ml+XMykbl10nRTk4bEcbzZbGhpjyK0XagUG03IXpMnjuPD4RBFUZ7nYRji8KbneQ8PD6Zp1im1rqlRlVUIXIkAqz4tkAqC4Pv375VuHTJPo9IIiTRNV6sVmVULvSvrnzA7Q61EQXBt2xZZwzOdTolVy7K+f//+9vYm0+F/WQebapfjOIZhVAoJyW0ScvQ9OdvI0ZyGrXh6esKEq+u667q2bYtvK87zHNudIAh2p2e1Wum6DuZvsuNpiLbKphCoQ0ByIVHn1FUHh3ovDgIPDw+2bf/3f//3y8vLNVxFUVQ+JbRarUDTtm2YqdgqTNNsPaGbpwfFl8tlFEX7/X5zegzDCMOwzoeEZUClFQLiICCtuinP89fX1/LsIA70zTmBT6fE1iOYHOBQ0BwWNicsBKzRgk2zOS9Kw4WJjB8XlaXMcIv69esXndUYy8aImqAS0iNwj+qmgoOX9H08xgbC6gDNjKZp5RV9uVF5nk+nU4iEsmAo59c0DYt68oW1LOtwOJQ1V6xNj9LwC6Q/iRri8paJVDJAPrKg9v7+jo2F67pBEFQWUS8VAuIgMLS6yfO8w+EwjA22tcZAnO5hOcH8yL4ZbzqKIpINuq6v12vbtisNvNgQsEciCq1GOA2SAXRMupCN/bOyIna0sGlN08gmxCZAkLYakB9xHHMkx3Q69U4Pdk6cnCy3Kq0QuDEClY6xPb0kR92e6LNkwzB8fn5m34w3DdzGyz/LeZIktGMwDIP9iU1/fHzQaGG/EJS1bbvOp5sl0lO68qgEy6Rt2xednMDxi564VWQVAk0Q4HxTWpPyneSh3cN6ve6EIJ9IGIacOYhfVsBfL5p0BOT/eDyy5+AMw6gcBtvt1nXdwg1LlmW5rismAmC4sPOAwLAsy/f9LMvOdgd9GpZlidnMs01QGcaOwO2FBH0Gg30DqHHsPScN/zSNuq5bnjfDMCxcN2vb9nq93m63I0Igy7IwDH3fp8ZCWjRpyHa7pVJwmS2jNCIoFKujQ0DTajcMtT902Mjn52d8LR8fHx2S5ZNSQoKPz5C/JkmCAVAOnlHQKcFJdEje+qur0DQYLfjVsZutJvn51NSvCoHmCNxSSLiuiwlisD0EcJFJSCRJMjB6zcdWk5zoC7YJYRjSwICpWViFUpMGcvJst9v1es0q0LBJ4hfxfZ+KuK47rh0Vp2nqJ2ERuFhIhGGo6zptePFJt2gedtCGYRCpFkTaFZFJSIx9UUlCAhp8LBo0TdN13ff9O5kBkyRZr9esVo0vF7Msc12XjPzKXNFuHlClGiLQRkjQxIQvXNf1hpVRNkgI27bpzZAJJSSGRJtfF/qCZENziy6f7Eh/rZQWnFUUK1l1XWc3ZCNFQLEtIAIXC4nj8QiNKhnTytpkfjthh7ihf5ESEvwOGuxXVkLIZHK4HsAkSchcp2na8/MzRwCwma+vWlFQCBQQaCMkjscjaUUv3Q1AtLAKqwJDA/wpk5AwDOPSLhgAYX4V0JbQ7oE2pvxSd/grgKJvzbKsSudgILPdbllBwqbvEDrV5A4RaCkkaBvB2QsXuMyyDFrXG551AktKSBS6ZrA/cUIF4sEwDPf0aJp26WZ0MIYFqYg15p81bh+Px+12C5Aty7oTu44gPSUlGxyDQq0LLA1BTdNc122CS5IkkBBdLXDCMGy9gpZJSPi+z1ldNumaYfKwbjy2bdPMJVNf9I1klmXr9RqzPw5McL4m9nTFnVt6+u4X6elzZtpaIUHxDzBezy4DyRe+qz3Ex8cHqm7XPWpiaodbi1JnvftVX1yK6sfHB+3jEfqQ/wGymTlC5VI2VP77QeBiIYExhw0EvnAOiePxmGUZlKr8bBchDgnB0ZTxqWEnxM8zll+FPSdBe0d0Vp3eQwmJdiMNGwUyV/B1UGxmtatoB/g9l+JM3RU7CXzSrNkZ557qVihJkmAcd7WHQFclSZJlmRISx+MR9xkINYK32y25/MOFn7/UrRs8QjVKWGZYPR6C5nJYZYNfua7L7xcOHfXTXSFwmZCA/2vhq64TAJAouC6mD0zrhESBvXLVUH+V34/xDcIBCcI564vZvN/rxo8gjRoFG4U4H/xPgM3MzzmKtism+0bgYiHRkCGSEP2ZVSuFBOrlD/1rdiENmz9YNnFcYMmm2tydAShpmqbkRCcDhu0CqJXqyNKZbfIgqMup3isEehESpJXqdT9bFhKwf3CahP6WSUjcfCeRZRlFE0IsjeZe0egOJSS6nYYKcXObrNLwwapIUN12hDTUODNqhU2iSbOxmWXtFoiTfOncgbrCMKxb7JSFhGVZnPYQ80pIEBTXJNgzcWe14XUVoS/UTqIOn9bvWRfYs8I7yzJygrr0WqTWHKqCY0GAM6m2ERKk7qQ9BL3h1FQHFr9sQeuNMAZ1pNj3WDexb8abvtVOgrrmyn0A+kIJif5GIMUBLHwv5RoLNwPSJ1zOqd7cFQK2bddd5XCxkMAHz45FdippeOyORZ89tVeOF6vrOm0y+E5WLE2ky7uQcp5RvAnDsIk+oau2lHcP7TaIxI8SEgRFf4nCroJvsdtutyRX1K6iv04ZEWXbtutWDJcJCXztuq4TOUzc7FTeDhfaC1uWVUmhxUQjjZCoBKSPlzjiTocewjC8UjyAyRZ910fr7oEmK+DhXECfarn5bObWusQy2bG/oStp1+t1J+N/FIBwlEAXCAnsGAp2CBwH7QpKVFFeNWOWocURwpxxLBk436eERMPRCXsSHdqqOxPXkFo5mxISZUx6fbPdbsnRAHc60Xa8sl72HEbl/bKVpaR8SWYbLJU4EY0ka34HOwnSKbELEwBap8lqB6JxetiyOPFAUccp/gd6kaPplkZI9GqToJ6FvO+2N9GP6DJOT7HdrdJdIcCaH85G4WVjqxiGwX7mXfEjPh32goOPj4/yXCR+E1pzaNt23RfaaCdB0pXW8rRU52xS2rGLVSe78NF1nbVVWKcnPD22bXNEvUxCovO2sF6tmEFYzNv1XV0pmJ3qhmBdKfW+EwRYZ1nDMMrbdLYWytz5d83WImYao7RO3S0mzx1y1V5IZFkGA1d5LsZsXqloInUQlNqXqrYR/BJyiC8Gjsejpml1s1vnE2uHXXIRKSz2LyrCycwaHqCJruxEDoUWP13pH9WiRk4RGp/b7XaAtnM4GewnNozKWfPD9vQMxpsgFRV02oJwNRgb7YUEJETl9rMuoBCwhi6I/r1oYUJzImrn7HxRVx2OSkgUkKEdITRLQy7thxcS/j8fcuOh0chP4GthN80FGMf7p+/7hMbz8zPny2LbOOQ4YesdOI2vQ8p+b4JkSyFBc0rlUl3TtPLWFcemsBWAFKEPsgmjyIN9XxiGZ5c8/IN1SkgAT2wd+rNLN+nZXoUEDO+u69q2TWEHaeBdmcC10tLMHbivggbDWZ91rMPOfolNxoDgeRDvpDynEdtsnF1Mbn0ENqXqBk60ERKQEHU7ACz2y8ID2iF2Cw+xcel8DTsEBybowVhXq3LmSystUxDkje/7LVSlrFoJLi7r9brh4rHzhncuJBA+HR7YBTFg27brurBa4d+LmkMFMRHQUgm14CXfs+6i6m6VGYswNMqyLPabLbBEjlI41F34VZo/sTYt69WpgZCX69NDuGF4cGQthigRETZxmZCg6yfrJAQsAZVTcOX2AjhehA6KlIUQEUEGzsgGk5R/1ImG3k2wRReW0vxLCAaD5fqVOLYLEAC0EKZvtdeQRBAbldsU7GWvb91gHVGoiD72s3sFcpM9m7NQxYj+xNoXoe/pNlnStrHKbXZOw7qhspkoMorhcYGQQKugs65s9vF4hEdjpQhhT2JTcbIx0JuzCT64+PWss2alGDtbtYAZzgoJIIwZk/71ff9W+4ZuMWS9M6l12Jrc5PNLkqQScMMwaELpFoG+qZGoOGuooIYbhnH2A+yb7T7oY/YvDDNUhGkHaVZIcD7PETnRNhUSrITgLNKRrXKHVbmTaCEkoKSqrAJRYKF+ybKMs9tge5rSWPqNLhRBoTtIC8/uG+hKZDmcdqA9L3y07unhdHofEweHJtabhZ2NruujO5LGap9s2y6MNxYBCo7Q5JgeW3AsaeyZEJ+fXYVcJCTIL5SD5E0AIW1qIcHRJf7rnAQdUmNxqWwGJ4YSDjEUSpWFBHwQ+fDVHf8pq6EL1dGf2NZkp+d4PAIUmDd71U4QAx0mIBWwZiGBhwSWAOLMm5Wtvsgcsl6vWdlgWdYoAiRQOAfqIMMwYB2pxETAl3ROAl8f5wstiAo59qz8HrlISGD1dnYu5dfYx69kaSskzgsJ2kNwjPvEMT7gygEE4w9pohCZDt8MimdZxi5+OSDCXEaVsgnyc8dZCvYnNo1tTUFgVrLNlrptGsIgDEPIAzgC0KSDBLtjuC23zWtvaLgu3KjTZDQ252HInKTRpi7zfV9wQU74sLsK/gWorKig4rImWCHBmqPLMxWG8bh0j2fUTSQhGraK1ceVBwS7BmRnN+TErgIu7Hw65f1HuS7+G7b269M4st+CDlTVaHKL4lSE04t8HET49ayQQHejsbTIEIHzK3ko+JiNyG7RvEfCMJTSPtGu60cahIYzvWi0avB9v8kqG9YCzmeM0JI4s4M1L758IF7wBIATemVnQHRxthqVpdiX2EmUd1U07YqZIIbX6zW2QfCrYZs2ujRHSJDbDJQzTQbh6JrPutiTO7L4rWDDxDYM+5gkyejsMd12hGVZfNf8bqvritoZIcG3UxWYYPdchZ/q/mSFRMH9iWNAvnInAY1NHUtN3odhmCQJqbZIbdWkLJuHNEjsy4vS8Pu8qIhQmbGwKOxT2QlIYq/KQkewtz5A7z8KlRoZAiHIC41i/yS1BF9PxRaRJi2ssboJwjwh0aQ8m+d6IcFOFiILCbbVt01zluG3Zaxh7diAsxoJNj4EOx4aEpQgG6vMqQx7I1ob2ZiyfIbJBWbs4/bSLkCfjnQ8dykk+IaESljZnQSrOcECs86ad/OdRGVbbvKyztHrJsy0qJRdWLiuSw6jo1hEt2hvwyIFHdQoLNusicW27bqPl90mNtRTNQRN2GwFnx1h+axjrGMhcWmICN/36cQDnXKAsOGTKuim6ppX+f56dVMl2Zu8xBVjN6m6k0ohJMhTi+Ns10l1oyNCVpmxnDxgGeZbICgnR2cwuv4qM4xpjWOpLRcR7U3HQuJKLLDKgJ6db6VUQgIjaexCgnVs5c8pon05Q/LDLtIbepEMyV65rkJMJ863TCYNRLwokxr7G/30sAiwtsxRtI69lLrA8L8O0xV+qPuz4J5Ul62T91cKiefn507YuDmRj4+P8R5WYpXv423FYGOADL9j0TGyLu98X0TKyTdpDAZ1VxXBOZ4d2+yY1zRtFBMRJ4iRzEJipBakrsbuzelst1u6vUDTNPYrujlvgjNAuJ29S06EhrAuW3yG1+s1NY1dd4vQinY8YJdM6nREdtA0ja5LgEW2znjTrtI+So1VSHBOUZyFCTHLzmZTGXpCgNUwtHB26ImrEZFlbb+jcBEuXH7HmRahWJNASGDbV1iMwrBEIy3LslEsjzhCYkJxvxueLJvP59PpFOjkeb5YLKggbSfxJs/z3W5HvzZJ+L5vmiblZOuilw0TURQ5jnNp6xoSHzjbZDJhFdYD135pdUEQYFRYluV5nm3bMvXFpWhcnz8Igs1mczgcYJpyXXc6nV5PticKaZq+v79vNhtEkn57e/M8r6e6bks2TdPZbAaDPHGCoT6ir5U4n0zqZQFJvIYJ1ibB6k+psmsSBZnMMbif5Ra8nc02igxjUU8fj0daKLDeDTL1xa0GDAF7jaFuMObZmcEwDH690OCfjVLOJzLwr3QWpFAv2nJ261CY6ApEbvInZyfx/66Z02kioEirFHGPTFhNAr4SKMvlssDPfr8vvLnPP9k5QkwE0jSdz+dRFNF9n2LyOVKucP4fSjzHcUzTDIJA2LYgiAPcn/b7/XQ6XS6XeZ5XMvz29qbr+u/fv2ezmciNYpn/8eMHpDX7ktJfX1+ULiSWy+VkMlmtVsCk8Kugf9IE3TDB7iQaFmmd7RrfaplWr4IvHtnIvrRuYDt9u91yboVkc6r0WQRYWwU8SgVX7tNRCajLKlfZuFQRpyz5pu+z+AyQASs2WgezNWKHwdqx2V9h5cYZUjbN5rlVmvOFCu3ddE3MIiUkhhlt7BmIys9mGDbusBYWecHPrhcu1uUc1Sb3J2Ft9Wcnd+z26FtwXdcwDJjx4fJEY3XIBTdVWpeoXN4hs9BCwjAMwrqubXXvlZCoQ6bD9xQ+XUAda4fNFJZUkiTUBaO4T5Q9QFC5pQDUbLbWM0BPvXb2RlIyV7C6I7SiEMyKc+9pT8xzyNq2zUZXY3NeZZNgUegjrQwShGocx5QWIZHn+XQ6/f37NywQZWMSy2Se57K6uLDNHD5tmubn5yeuk9vv99++fXMcR7ShwsKyXC4RS1zX9dfXV8dxKg0Vy+WSHGodx/E8rzIbS3mw9Pfv33/+/MmpzjRNYh4HJuhCSXj6UVnR5re///6bePu3BCsxmqQRgYdy1jkg4ZoKytYucY0uXqadxNn4zO3gbV2KFnoNg8zL1BetQeu7IA4G4dseRZhuuqGSE6klSRJy2bjGQtk3+A3p40OAJQnNF0dP2LG6iXWWwqAsY1T3vpyT80YJCQ44t/qJTsw2VzHh26jbzN6qIVLWywb05nz2grSdVhsFPUyBPYwfTCniN6rAfOFPUg+KFq6jS3UTxXlGn6F3Ka0SciPgOA72yK7r8lVMZRxeXl7KL9WbbhEwTfPnz58wlkZRxPc97bbqFtSWyyU5a8ErtNIFlhxqW1QhWpHPz0+cCthut5+fn0Kx15m6qaA6qFsCYJN4jXMe7D+tKRT4LMjzcf1p2/bNt6VJkmB90EKVIVNfjGjk0M3E4kcRZl1gwW3lhz+iU6UjGidglbNFu9hwjZkiiiKSgWyaXpJSgt5cmsCKVeQIBJe2qHX+KIouDXDSuq7Kgp7nzWazw+EQhmEQBGzolMr86qUICNi2nec5zNqbzUbkXQXxBhfYzWbz+PjoeV6apiIgeSc81KF9sZCYz+eaprHkyLJUhvKaKZ6tokxZvRkMgeVyiVA85KQxWNWqousRsG07juMwDA+Hw2q1enx8rFzVXV9RJxQgGLD13Gw2s9nsUq1mJ2zcJ5G6gXGxkAB8rPNWHekWQKdpSg58bBUtSKkinSDgOM5qtcKFJFd+rkrqd9Ij7YiQWl/XdcdxEEClHakBSoFbnEpbrVaCczsAILetoo2QsCyLFQycncRFbXMcZzabOY6DUlEUdUX5IjZUZiAQBMFkMkEv5HmuVExjHxik0nFd93A4OI4znU7ZD1moBk6n0yAIjscjTjg7jjOZTDRNY8+ZC8WwxMy0ERJQImHJr+s6f5zFcYwlZBzHlDPP8yiKaNMAfDEf4T5xnJ0pe1I174lryjavRdacjuMg3Pd6vcbe//qWqq3h9Rh2QiEIgjzP4X4qfqxAz/OgLoOZ8/HxUfADg530kVhEWljhMWvQQXP22ARRwxBkY0ej2eyRH0ScRxHkp3DBbBVE89JEJWOXEhEh/8BOHfgaO3ThrOYAACAASURBVAzzgN5sfq5CBMzvhAfarI+id9ipU7RwHWMfMHXHov/3Sp4WbcNiHy5TGGRlfzVMNMfjESv6wmV+ruvSvRyUwAiwLItES5lsc24xMTXPL3LOwU5c016+4VHqhqApIdEQqJtkQ5wM3Kcm+Myrnx4EfUDEC8EZvkmHtqvUtu26aFpthMTxeDQMA6FlITAQsjH854MT54iXi6kf0z2OPlC/orPxE+iEYYg/IVratRalsiyTZidxDQ4NyyZJgl67JqhiXV1KSNQhI857OlRhWRbn5tHbMmxZFtamWZZRBHJh48XeFqtLa+dEPWkpJKAdguT5+Phgt4FIPz8/g0v29DnpoCAt2Cur8AZFkL5+/6uERMOBgn7pVamladrwYTloOdIQB5WNgkZwjlbdECXDMAqM0dDtUDt6wwbesOoCsCwnLYUE7o6mY8BZlkFr5J4ediXCHsnGihJLgIJCAzsJcAalR93eh+Wek5ZpJ+H7fk/zHXtfEPUmB1WhfsLGFUH5KVpceb1Cb5ATx9dRVqjmCMIMG8GUE3rvJtwW7mMADxTYAxozdvK5CZMjrbRjmwRQsCyLc5kRIcUKCYgWdsVKcx9rQrBOD1FonZBmJ8Ei1hqNQsEsy8ho2Vo8hGHY8Jukji6w0fDP7Xa7Xq+hnyS2afanBIkBNkG/ViYw70AMX2MDa9iQsWSjy38w9gRBBt1aiSECe6CLxb/brrIJt33Zi5AgKwK/bQUhAWGALxNmKIw/EhLYAbSetlhmlJBg0WDTOKaEldcA33/DocJyeDwet9ut67oFkYAdA3TTcM8lO1ahOOdPdgtSoI85sQVNTnXj/Um0XQVHSBDINLaVrYIwaZLoRUhgW8DRZIGzsq6Qohmz9nQSEvhoUTY5PU1aWJlHCYkyLB8fHxRWazAjATq3uZGJFM209udY1cptbPemXGlhfdOOrASl2DHTvBP7aHgTIYF6Wdl/5S62j4YISLMvIQGh3dW42W63CDVKgic7Pa0B1TStoTKkdRXDFOxE3VRQ3Q6wgSBwmggJMmuxggF7hYE7MQzD9XpNRg5cvXelhYygGGmCVU7e0P2puZDAThSdiB4cKfKDsd2XkKBjEDStX9MkTCWapnU1f6mdBLqDjcNMd7Jf01MtynJW5axOA+JwYKnAaU4Yhqz6Aud7OPnl/okVFawaYLBWXyQkwBX5duu63tVydrD2DllRj0IiyzKoL2BmaDG/h2Ho+z4t3DrcG6qdBELf0HnGG06+ZSGBnQ3pvqBQajF+BviQttstO0R1XXdd94ZgDtBkThWsqBh4V9FCSByPR3YP3dBQAV9NmL44m8jtdguWcCaMBQ1+FpXuYTCJiTZ+ehQSwIU+9RZTPKkXOnel1zSN08Fsjwqe/vj4aNEQcVTJMF/ROi5JEtYAcJM1abseT5KEThJg3Pq+36Jr2tUuVCna919zx/ClLWonJFALTl+h18qGUpYTNifys79SGufDDMPAkKDhfTwe2eHNngZjQet8uiPG2iV6FxLH4xGyt8VKEB4stm13LlqlUTdd1Os4jEo7s8J5lItIdZgZOwl2TYcLyEY6w7JGC3KI6hCuUZBKkoTiuAyjxrxGSABSOqeNXqucr9haMK1Xrn3Zk33Ihirgn4ki2+2WPdth2zYiGmAL0qFq/foBM4SQuJ7LbilAD9YtTcGpFaJg9SF3WyOAYI5Yl0E7XPl9tqZ/q4K4943a1Ynr9q3a0q7egvapcj5tR7lcip2+y782fFNYqZQd5wqG7rrIaYUtFH1uGBLEDCs/6CVW1WMREm1CheOTEPzf3W7HusEJzi2fvZeXF4qyXs6ZpulyuTRN03GcKIoMw4AOJAxDXCNYLjLwG9xWhFDhvu/neb5cLq+5tXBg/jnV4d43KNAOh8NisZhMJo7j3M8NS9PpNAxD6DZ3ux1ij3MQu/lPdFMF5ocoipzTQ4wdDgfSn2PDcTgc6FdOgrI1GdsYIU1ycmrs8Ce2yQWy0gqJQjvH+2ccx79//y7zn6ZpEAS4qWm1Wu33eyj3SWaUiwz/Bl/garUiddOV19sN34QmNZqmuVwuaYkaRdFsNguCoElZOfK8vLykaQrt036/x0TcedNg8u2KLOIFVE6OlVdsFW7EuYaNPM/f398LJoprCF5flnfdC7sDkikN3YsELcJxZThfwiOzsEMq75dFaDW0sRi718f0FaFFDXlgNfWwu8ihWGvYfLpLDgrGXrVPzVni5yxYQwtekbja/Xg8whmaSBWy0WdYmHnK6qYsywoaLaJ5wwTHJvEf14sgYSlEUYQrD1kOLcvibPEwnR0OB8uyDqeHyhqGsd/v8X4+n8dxrOv6brejDHUJXdcLap88z6mgruu0S2UppGnKyvbN6UEGmMIMwzBNsyAwWAq3Ssdx/OPHDzCPaHrQQtyKn4HrNU+Pd3p+/fqFjnNd9+3trXJ9OjB7A1Tned58Pvc8D/1uWVYYhpyPrjlLcRwXPqXmZTk5CzTxjVN++gyD00PvNU2r/HLZDJXpxWIxn88F/HIruf3flzeUXb1WLdRWrhb9xj/A0tBfONiu+oKFnTyXyouprqoTnw4LSCdnTsVvMnHIeoKyHqKU4dJE+bTNpRSa5C+sqesqLbw3DAMDHmGHqKLC4H9+fmadYinbzRMcWSDzTgLNLszDcRxXyn9YhrGIeHh4oOVDoTj/z/KmQdM0dt/A2cfouk71stnyPH98fFwul+IvPYIgwM3Yvu+/vb11snjkAy7+r3B6wbI6iiIsscXvyk6AXS6Xb29vi8UiiqLVavXr168gCArL9ksravdhXlSLZVmbzUbTNNM0MS18//69TCEMQ8/zHh4eptNpFEUPDw/YKZqm+fX15XleEARxHHueR93tOM7v378tywJZKAPKlIV7c3MJ1hMDBQHeUy3DkC042w1TafNaCk6QZRW8TH3RHJZCziRJaLJoeO63QGG8f9KxgCsNFfA1GgAHNhYLbYjL9bLX+bHDHnZEzPXs9hHabFYGlGne6g3n3gfJ1U23QrzbeskFu1uynVCj70TTtLpTAlmWFWyDnVQ9RiLsdAll/Rhb0Y5n9oYry7LYWbUhwcGEREN++Nk4AecFvPPKsqy65ighUYeMen8eAYpRYdv2YIHHz7MlfA7WUCGmhro/COm24xYNH5eQ6A/DPihzukPmcxKVHtDsXm8saQFPZuV5Pp/PcYADy6KXl5c6PKGZrfv1Dt8jVgp0Gv2dKhAT2JeXF+gn9/v9ZDLxPC/PczFZvSuuvr6+atvbh1ASgaZMenDRbBK+70O72jAwlEx90fnYpvBH0Ne1UMJ0ztIwBNmGl6OoVvKgdhKVsHTykqNuknknUSsY1Q9tEcAGAieosYFQLkxtsfy/cp7n0VHtxWLx+Pjoed6VNEdR3PM8OqG22Wwmk8ldnVEXrY84H7ISEqJ1ltD8vL6+4hjgz58/yVdHaI7HwFwhmtBms5EyeEllVwRBQPcZLxYL0zQ/Pz8rc6qXt0JACYlbIT+yetM0dRwHEaKOx+OdHB4euJOwObMsa7Vawft+YAZuUp1pmmmawk1uv99/+/YNEQ3KzKh1SRmTDt/U2T6VkOgQ5B5J3da4h4h1URQhQE2P7bx70ggr6/v+4XBARNU4ju8BFdu28zwPw9AwjN1u9/T05HleYdq67VcgfS/UrfykFRLSuDZhaPJ8D3oevMvl0nEcBCFgfTd7rvauyS+XSyhh9vv909MT8L8HRGzbxq5C07TNZjObzVjNG0dvfg/g3KqN0gqJOql4K6CvrHeAaASVHHqet1qtELSS/VwrM6uXHSIAJQwcyaIomk6n94M/vOagXFqtVqR9UjuJDgfYBaQ6cZ8SkAguERSQsRYs+b4//IllXO0Hp8Mr/TKVC2yLTmeLsFEi6o61s/mlSSdJQm3HpMZGuZCmmSI0hAOsOnEtQgcJxwOplRo6sPMboIQEH58mv2ZZRsdTcP9rk1Jy5CEXYU3TBLyJQQ6QlZAYdz9+fHxwoox13jaK8NxVpA0lJLrqI/IWxbJ6yFHRVRNa06F4HogSeFdtbw1a84IcISGtTeICjZvwWb99+4bza31zmue54zgwQmRZxom00Tcnin4lAjBU0Fnl2WzmOM6duD+9vLwgPCKu/3p9fVWH7yoHSeuXCGBeLq6ERBmTO30TRdHj42MURa7r5nmuPEmEHQc4q4yDBVEUPT09zefzezDqzufzMAzJ/WmxWKjoTx2O0oLDMVFWQoKguN9EFEWmaTqOY1lWkiRqgTaKoUAHC3CN7uPj4/3sKtD29XptGMZms3l8fFSD9vpBW+cRqoTE9diOmwKOQeAodRzHdQNl3I2Ul3tMl/AWxa5CYk/Z5emhzsRpO5jQsKuoU5hQEZXgIFC3k5D2+lJ8NvP5vKA2wc2g7P2gHNTk/ilN09VqFUWRrutBEPQX8wCXJskN5m1bF4ZhHMc/fvzY7/er1Wq/3/u+L5+8rzwttFwuv3//jgtisRv+66+/5Gv7AAOsDjRpdxJQ0dql53A4RFFUec31AN0gSBUwUM9msyiK1ut1nuf9SQg0+cqbjQXBTWQ25vN5mqbb7VbXdYRRoTNoIrN9KW+VVnq0NEkSy7J2u91sNpOy7ZdiVc4fxzFivUyn08m/P1EU1e0kpD0ncTweNU2Tw0+uWxfYJEkwep6fn+XAp7mf353kZEW+7/vStNq2bc4FamgmHfFRl09U9rvv+x8fH/7pYT9/wzDqXN6lFRK4i7wSprt9yV6wPPDB3eFPjN9tL6PhWFZjNSDNATToBc72bOGcdhiGZ4uoDJxzEtIKCewk5Oj7JEmuDIyRJAmtLhteJ9chdOowXYdgXkSKXRZYljX26bKhkABEuCQVYtJ13Su/oItgH2Nm27brNp3S2iTg51CpwSyr6gR/M5vNFotFOyZZ84Nt20mShGFYMOa3o6xKiY8ADhZst1so6x3HGfs1Fc2Pg0yn0zAMsTzabDZPT0/L5bJJ8TzPo9PDzxxFkXN6yi5VKF5Q8RNZ/Cr+4PkXh2MUeg15Fu1q6IZsl7O11q6y+tkbriLVTqLcp8O/oWgrrYfT8DwXasRuuPCyyZ+0jT47J5DFDlMkq7VnK2LBRBR9+pX9iV2bFy4vuOH3SKyyCc5OQnJ1k2g9wfZK8/SlXzUComFQ6ro+sPmh3C4lJMqY3OQNq6w3DOPmA+NSEFoLiePxyGreXNetm/1t24YJJwxD27Yty6pkUtd1UuJDKiAbjsHjhkG8JzWXpmnr9Ro/hWFI7yvpD//yHoUEQoXfm5DAyKZ9YicxXK8fr0pIXI9hhxSyLLMsC4NkXDbta4QEAGTbbllW2Z9C0zR6ybluoOA5iRsbj8ejYRjsR8fagTRN6BU5R0hIa5O4N7V7EATT6dRxHFwyikWiilVA8lIlCIHpdBrHcRiGlmXhktTpdMrXv1PZsSeo7XRJqmmaBcslnempm0NghKg8erbf76m4pmmojkDDyYT5fD6uD1NaIUEdI31iuVxOJpPFYnE4HHD0IQxDz/Okb7hq4DUI4DJtbPIOh8Pr62vB0HoN8Z7KFtT6rWthL0nFBbFl4zMR//z8pHSLRPnc7m63WywWI4qeooREi34fuojv+6QfoLqjKPI8bzqdIoo49s6fn5+VCxwqpRIKARYBKNZd193v96+vr47jiLylqAzLwTbnojSixUDwFDYTLB0IiTiO6zYWbGZOmrRYx+PRtu3393dOZqF+kjZ2k1AoX8kMu+iA/xwtfGCX/v79+x3KhjzPd7sdsCVA2Jcc2NngXdB0sw4wnIJS/hQEwXw+9zwPo8t13XHpQ1p3CmKZ8APj4+tDfDNUxB8qdVKWVUOFYTiZTFqzPXTBDg1iQpGSwFgahiFOz/u+Xx6Xz8/PYzHLd9IXgMIwDNbFsNuvBba7Or8XoYZ3H8wkSfL8/AxIOUEa+qi6IU2c8mmYuV22gkW6ztpccHulUgVHxDprMBxt23HYU6k6Vo/Ho+Q7iYK41k/PdDrVdb2gKzQM4+vrS9M0vMfho07mIFR6OBweHh7q9ssN17/wz5vP52WZ0QmrQhGJogiwFMBhL+mjDQEAMQyjyY6qvO1I03S/3+M9LuYDyLquswtAofDpnBnTND8/PxEDbr/ff/v2zT09V6pZOuTTPD0dEiyTsm17tVphWeM4jq7r5Tyaprmuu1qtHh4eptPpcrnUdR0Dz/f99/d3bE3SNI2iCPGBNE0DtSAIoigKgqCsQK6sSISXE4SvEIGVbnnAeUiaRIh4HMcF8UA/DZ9oMgdZlrVYLNbr9Xht0eiLsyMtTdNfv37tdrtCHwEliAHbtqMoKndrJ30HIUHCiWiGYdhQ/FCRsSeCIKBD/r7vv729iSMqesU2z/OnpyeaIgo+5VR1nufv7+9YrxiGsVwuad0GD0PkhGss0jAigrKu69vtVihIHcdBQ6iN/0r0tHm5OdmGKg72SEuWZVDghGG43W5xpp8aAi0E5aH3ZxMowpqtzhYpZDh7TLSQX7Q/OX2RJMl6vaYPDOPSsix8XdeAdj0IUG4Qbzg/NRYVXyfNd10XPSLIcQqcRLu+aWcpoCJ2cjhbhM1Qnj3oV8wt9Kc4CY66SejzHdcgyJmYriF7k7JSCgl00L9WK5r2/Pzs+76AJgHf90lZD4YRBesmg2HgStkwFXTGeGAeqDrDMG7OAzEjWYJjhVJCYgR9PXYhcTweaQFe3r/jVGrrVdtg/bfdbgubHhyvFZ/z6yFar9fQzt82mAe2mNc3R1EoI1Bnov9fLXE5txxv1E5CqH6EbCAzIKbX22qTrsGnIOosy/J9X3ppwWqfbtJ3SkhcM2j5ZXVdr8ughEQdMgK9H+9OIgxDmlk0TbMsS7LI/mzYODRQwNht3Q5l6lDDMGiD2G0VddQuuk+ijoh6X4mA2klUwjKal3XXCorcAHgEkclBelVy4fTGwLPnwCPhVoYKJST662iOkFBhOWgeEzfx8vIiLnP/zlme557nmabpOA4OhbiuG4bher3+94yy/bVcLo/HIxktcL1Pw1tuRoeFaZpZlvm+r+t6FEVjv8hodPgPzXB/oum2lGWySVx/fWnffQE3X1JEQPFCaut7WwBmWUaWXhzBldVcgZtLMGfB97/XkXZvA6lXMAvE1U5iaNHbbX2z2UzkcGCO4zw+PjqOs9lsNE2DL10cx/dzVrnQ3dPp1PO8PM+hg4qi6PHxkQ3AVcg/3j+n02kQBNA+HQ4HXBE63uYozisRUOqmSliEe9kk2sTATOPOXlbVgEi0aZqW9WN1Uc8G5nng6pbLJa21V6sVplTxI3JfipJpmlC16bpOzbyUiMp/WwQsy6obmUpI3LZrmtbOCWXclER3+dI09TwPu4fD4QCrQ5IknN0DxTnojotxUIJggAZf07TFYjGbzebzuVAd2gmU7OZpsViwq4dO6CsivSLA66+CZkqaP2WySYjjArvdbsnwoOt6w8MBiGwhzdC6piGsyyx7t+U1NAUsyw6SDh29lE2iv76mG1jLVahzEmVMhHtTiEt8K/7YQBoXubQqIVHosiRJ6JI1wzAEjERSYLjFn2zg8YtGC6cuJSQ44Fz5k2EYdRSkVTfR4d5e92iDEb+tTSIIAni1IkgyIt81b7tkfdG84XU5TdNM03S9XhuGsd/vZ7OZ53l1GuE6IoK/R+DxMAzhJjuZTKCPEpztu2WPM8NIKyTuVgne7SjHrY2LxWK/3+OwNATGRbWovqiEC4Jhu93qur7ZbKS0Vdi2nec5TslsNhu4P92nF0PlGBjFS2mFxCjQb8hkkiTDX1ECYfD09IQbUo/HYxAEQkXAb4ie4Nnm83me59vt1jCM3W739PQ03otD6qD2PA/uT4fDYbVatW4j7u+qq0W9vwYBnuSu00Nd+t73fdGiR3COh1zaunvLTxrzTrTJyibRcPyQ1ef5+VlKQ8XxeKSg677vN4RFZRsAAc6XrgzXA+B/bRW+73foIsLhJssyumOnQ98baBs49aqfCIEsy2jX6LqulKIChgqcu6Rj+YSAStwEAXXp0E1g76zSwu3qndFlCCVJQm6LlmWt12vmx2uTyinlUgQRFgnaA9u2pYzqQcsRWRt4aaffNj9HSCibxDV6PBnK4gTcbDbbbDa4pTKOY/nU4uPqquVySUdSENVDPvcnXA9sWRYa6DjO2QOGeZ6fzTOujh4Ft0pIjKKb+mJyuVw+PT3tdjtsVvI8p/Vdt1VGUdQtQemp4aj28XhEj8D9SbIAUKZpxnEMS0wURU9PT/xxslgsfvz4IX3Xi9ZAJSRE65Eh+EFA78lkslqt2h19uJRL0rNfWlDlD8Nwu90CwNVqJd9pA9u26SC6c3o42wWOO78aKn0hcFtFWH+1Z1kmjXdThzYJ1vYA8TCAvjs7Pf319Z1QDsMQoqJ5QJRxIUMatrr46sq41V+H3qNN4uvrqy+5Ok66aZqapgnbAwWDGubow/T0jBM2gbi2bRvKGZw2eHx8lOyQ9nw+D4Lg4+PDMAyJ46sLNKSasSKtuglCgrNvbYaPELkKy/+LeMrzHMfiZrPZfr83DAMHp3uyPVzEm8rcAgFWOTObzaDWb0FH2CIvLy8/f/6Erx0Cj/MNFcI2RB7G+tu/3JbydrvVNO2evbDDMCSvVuzfb4UG1Ai3HQ/y1U6GCtwDOIDacGAMKb46GkjatoHZuJPqOOomaQ/TwSYhx5djGEbDgwsQDGxAvQ7PxLX+WiCrWhdXBTkIsLtMWQ8csGsdaQyNnD69yU8cISG5ugnOnWPf9+33+7MNiaJoMpngDlEE1Ht+fg7DMI7jm2uWHh4ext4FwvJvmmYQBORF+vj4KJ9yBjekUjyPl5cXyYwxwo4uMCatkLgTV7k4jh3HgXhAj+KeuO12+/n5eXPxAJaUkOh7FsAeAldqw4tUsmkUgcfh3PX79+/X19cgCPpGVdEHAtIKCYk7OIqiIAggG3D+COZomByCILBtez6fi4OAih07QF9Mp1NcqW3bdhRFs9nMcRxeaM8BeOq6Cgwk27YPh8NisZjP5/Jtm7rGrAN6/9EBDUWiOwTYQZ/nOWmZEDyD/kSFlmXBc1zNwt31wLgpTadT6Bh//PgRnR7f99/e3mQaIWEY5nkeRdFisXAcR9d113VVIJn+Bu7keDz2R/22lCeTCS7+ZNmwbbuFJgqHDC4NihBFEbxOy2qfz8/Pv//+W9O0/X7Pstckbdu2YRiXMtOEck95oihyHEfikdYTbleSBewgEoZheRBeSX/44o7jRFFEAylN0z/++OP3799w3qNA68MzJkGNjuPUzSryCwmR+0/X9bJeCIpXdunnOI5lWVgrjfFTV0LiVoMwz/PFYoHtKdxXWqyQbsV8ud6CkECGKIo8z4Ozhnt62G+nTES9qUSAIyQkVzet12vEuiFFTRzHLW7TtG07TdOHh4fm4w+3aE2nU1ZrhCVPZSfxX06n0zGKB36j1K99IwDtU57nmF6jKLIsKwiC8tKkb076ow+NaxAEi8Vis9lEUfR2epp/qv3xJgnlm/jkDlBpkiSapslxZ8vHx8eoGwI9wACdrqrgIMBqY4a5w4rDTLufsE7ilGUXUiNtI6d1vf50j+ckoOsf9eaaliEvLy9yNIRapBLDI4D7KXEwzXGcMcbzwKaBA114enD5rpSuwJy29/eTcoHtD9vOKEdRJJkvY2fQKEIXIhAEAS6p3e/3T09P4/IitW2b3Q9VNh2a4e12q+s6uQJX5lQvGyKghERDoG6WDcflyKZyMz5UxbIgAEMFLpre7XaO48jnPzqfz/M8x+XqURRNp1M5Yn32NwbZWD6FWiQ3XBdaO8Y/JbAx9n36N01TOl8Ct+NyR0NbDUUEfh2RD3G5Ode/geoGBu3NZvPw8CAfIJ7n/ed//ucff/yBbdPHx8fLy8v10ElJgRMWQVoXWGncLuM4fnp6GrWfu+d5m82G3NvbfWNwS4N/DvZV7RzVCrWT8DBN07KsO3SJwQH+3W5nWdZff/0lsvXL87z5fM5apwu9WflnnudPT0/waUS8yzvs5Upk2Jfww2TfUFoJCYJC0ESe54+Pj6MWEu0Ednx6drsdRxjQtM5OHJZlHU5PZY+y1Ng0ZcY59vnpoZfSJyDIKSi3mNPoZDJpYpao7KwgCN7f3+HP4rouJ/RTEARYheAsNwcKnGnXNE3XdZYgLIicgpUc3val4zi19p5e3aqGJJ5lGW4DRqXSuF0i5vnY/fmaXGURhuF6vbZtu6Ae1XUdU0MYhujlbsfVdrtFiHUcZWe/VcuyGgZp75alm1DbbrekjhOz1ThmdA046/WaRpfrumVSMGNgocCvDuEUobXDEAW1LMt0Xbcsq0xc5DdwfqvkUJ77JOjb9n3/eDxKIyTQkLELicrBRy9936fpie1HdCVlGyaRJEmBH44L+TAsDVkLu5wUbdTxZ+2GKCVJQoPNMIyPjw+2INvX/DkEIX9Qlj2VRcRZst2m+/guylAQz5IIiSRJMKCxEFBCgjpY2ASuq2PVRIZhYMcgzlVR7MJT0zSEYReHvZ46N8syxIZBk8U5yNmJkABodM9d4V4/XddZ0WhZVuWMjIsvCznxp23blmXput5T7/RE9i52EsBuvV4bhiGTkEBbRj0x2baNrf12u12v1zQBYdNg27aYyg36Gl3XpeUh8cxOEJRTpkSSJCTCLcsSYQR2KCTQU2EYomdpTtc0je1ZdmORZRnJy/L+HusbIsuZc8UcJByGR3NOgnWjxFUK5Quq0jRdLBZ//vknqSzkSNi2PS4jWBn2X79+vby8PD09LRYLWAVt28YnF4ah4H76QRCkacpOmjDFTyYTcr0tN3nsb0zTxAFmTdN2u93j4yP7DY69deAfJ+9s2357ezvbosViMZvN2GyFpQP70xjTdT7QIxASuJhztVoBd5yLsW3777//ns1mdBQ5z/PX19f1ei2fK/Sob4+JoiiO4/1+z4Z0xvZIpBVPHgAAIABJREFUZG/L8keOSROc02bIcZzpdAopUi4iwRvIcgTzwEVGNzyVRrB3C2wYhjQ/klm7XEU5KAgnzr9Uqwcx9z7EFXR/tM2ERwHtfC3LIhcFdCHWPjKpm8obWwJH5ERBrVRQ9YrMeUPe4GfFTlvwwW9YfHTZEMwDUyduSx2+CTQP9Fc1qzU6Ho+suqlQaUExxY5wnGYv5Bf8T05LRd9JzOfzwhEBVvdCceThsAwlgOM4N1zslBcgV77hLG2upNxH8TzPPc+bTCZQK+GTw7UZpODuo97haSJ4exzHWZbBaIEQF5PJBNHph2ep1xoRzCNJEsuyoih6fHxcLpe0j2er7u/rG+ZboE0A7r+r2+8ahkFnI6IoOhwO7IphGFZZ2PtLiy4ksHY4q/srWF0kCGVBXd7i9gsqO3DCcZzHx8fNZkMrPgh4mbqjDCmpm8gKutlsHh8fHccpZx77G8SOxe52tVo9Pj7SlEpN66+7+6NMzPu+D/32ZDJ5fHyExx39yib+/PPP3W43OT2O47CLV1hx2MyjTo8mdlN5LDbBvfyh0hldtrhki1y2aX2n4zjGTS+oCKeK3t7e6tZfffNzQ/pQeMIGA0xwQth13QFmtyEbDnUT7rxzHMd1XbbHcR/7kPx0WJdpmtvtFt0H/WEd8ZeXlyzLFosFTlzDbIPMtm2Tnryu+Ijej0ZIXLp9w9avPPvHcVyOqFrONqIuvBWrURQtl0vY7qBQkm82bIEtREVwet7f36EIxe2zlQuUFlWIUITuvNv886GwSGOXiM0jsgCEyu4QB4SGy+s4jusUNqMREtQTrMaTTVMGJMjPsvB+dPLgUulYaG/nf+Z5/n56KGKaZNdhdoWYd3ryPMcs6jiOfNssqNpc18Wu4vHxEadeevLYzvO8J8pddbpodCi61FWMCW5zB3twajoejx8fH+xpeMMwKs9DKu+mnroVymgMOM45/kLtWFwXXt7bnxjDgK5gQpMDCnZssOfROmwdxwOnw1rukxQH23GE5YD9E52HrQC0SbqukztsoWsxZAsvx/gnXA+bBMjrtXVsGAPDMC6aBZSQQNfAD4rWdIKfM28xnLIsI1lIhyVb0FFFhkdg9EKiANl2u4VvJWfqRPDUQkH1ZwsE2PhFlmVxMK8jroREARn2EMmtzhwUWOrwTxFOVHTYnDshJZuQaNJt0uwkmjS2pzwfHx+sLeui3QPLkhISLBqUHkA/Q3UNn5C7dcPj2XeNHCExgnMStD2/KCGavfci5tnMeZ4PH5YjCILpdPrt27f9fq/rOhQjo7P5szAKmIZlAq6TjuOYptnQEUXAtpRZYuN5oHUcN5Ny8co33ZhhK0mrl/UISCskRnQGrb53/u+XKIrKbrtnS7XIkOf5crmcz+eLxeJwOBiG4bouDlG3oKaKNEEgCALYKvb7veM48/n8+sm0Sb0D5DFNMwiC7XZrWRZumXYc55oogUpIDNBr5SqkFRLsEflys0f0ZhifvzRNPc97fHxcrVa47ni73aZpSoEHrkEM6qZrKMhdFo6kx+NxvV7vdrunp6e6iBdjxAFiDzF0oyiazWbz+bwynscYW3cPPEsrJIaZWyUYImmaOo4zm80QTgNhy7o9NKuERMNx4nne8Xg0DKMu4kVDOgJmu4fA4wLC3glL0gqJTtCRmwjiIc5mM6jCKYipsj3ctt/TNEUXOKdHGu0TRfSCGebHjx8y2WBuO2Z6rV1aIYH9rNrVlkdPnuewSzuOE0WRruuu6yZJEsdxT+IB1ZU5UW/qEAjDEIaKKIqenp5kMlRomgYzzOFwwG0cMknBug4d9XtphQQM10rpRKMTsgGBWmGXxlkTvO81Ht9ut5PJj4Ag7TUBQ0WSJK7rwlCBCKy9VjoY8el0ut1ufd8/HA6QgtcYtAdj+z4rklZIiBNg6+YDK03T5XL5+PiIADvY9SdJUrio4+Z8KgbKCMBBKEkSwzDgIER3qJUzj+vNdDpdLpdo2m63w7V342rCnXArrZCQqf/YkBjN2wXfpOl0OpvNcPkrrA64Uq3XrUNzJlXOJgiYppmmKSJerFar+XwujTbfNM2fP3/CShFFkUxNa9Kzo8gjuZCQ5ltqOJiiKPI8D7GOZ7MZ1EqWZfm+v91uYXW4iQpOmrONDTuij2zL5RKGClyBh8ve+6hoYJoFD2BlqBgYf1THOTMguZC4CdzdVhrHseM4nMN0cRzjBNxkMnEcZ7PZ7E4PgjZnWUYZumXsImqthUSapvcm6TnAYj5FHD1o86+5rDc+PZzqBv7J8zyEZVOGioGRP1Nd3yFBbkhf08YR45YPUZIkhSvXkf/j44MibrJ97Pt+kiR8mgP/atu2YRgtKqX4P63DRrWodCxFoM1H10sWexxjHk1j7wJQQcD6G5xhGLJQsxWN79IhdkI8m4buBfemnc0sVAZaPsOLNwgCehPHMesshNbZti3yxWe4wO4ihLGFgpm9J9/ci/gRLTO0+e/v76vVCtp8z/PkAMo0zSzLcFnTarV6f38Pw1C5otxqBE6Ox+Ot6u673slkUq5imJk0TVOaFs+KKDZzmeHyGyI4lpPMUIm0OLMynU7f3t6k8ecpd2Unb+jyO03TLMuS7KJAz/MQCwD+XZZlqXMVnQybApEoiuAGWXivaZrkQgKKDmr2fr+n9Ti9FDNBqiTwjEC+I3VJwqm9FsuR5ekRs4NE4ypN09lsBq4kc25G5BhadbUYSKJ1loD8cISEtOomWm58fX2xypnmq+9OLtS1bTvP80qzc/M9zWQysW17pBKi9feQ5/nb21vr4vdW0DTN4/EYBMFisXAcxzCM5XIpjfbp58+fUKxpmmaa5tvbm+d599bFt2qvtDsJBCZCWE3otTVNi+M4iqIgCEb08WC7PerVUwt1E1RwI+qmW33A5XpxAQnWJZZl/fXXX9IsLwqKtTAMb+LPXcZ8RG9ozarrOrt6juP44eGhWrXLWrFlSktzMx2cAkfdNZjrL2qCbdst7km9qAq5MydJQm7HlmXJ1FiEKsG8rMbJRT273W758qySmgxOopUNk0ZIVLZuXC9bCAk53Jdv3k10+6xhGKI5Rl8JDrlHVzqIX0n8DotzXGDVYTq+ZL39r+QJentWBuSAZrcB65Swqp8/f8IDYr/fz2azsd9lxDrIYQ+BceI4jud5KkTgNSM4z/NataSsMlOancR6vR77shqReWQdaeK3C8E8oH3SdX28JxMR1r4AOK5Hxfzoum6WZYUM6s8mCHDOY6qdxDXSd4iylZ5RQ1TcXR2u655VhnZXm6JURADBPBCaW9M03KQ9xnX34fQUmofLNrCW2mw2kl3+WmjsTf5UQuImsF9QqQSKl+l0qo7LXtDl/WRFaO6fP38ahiFfaG7c/Grb9uFwwOWvY5SC/fT8tVSVkLgWwb7Lf3199V1F3/Tz09N3LYp+EwQQdZxCc0NyNCk4ijxhGG63W6yrcEGFEhXXd5wSEtdjqCicQWBxes5kUj8PiEAQBMfjESqa1WolTdRxTdPm83mapjBURFE0m83UsbsrR5YSElcC2Htx9sBL75WpCu4JAc/z8jx/fn4+HA4/fvyQadENQwU8uzabjeM499SxLdtafZJO05SQaAmoKqYQkAOBz8/P7Xb79fUF/QzFs5Ggdbge1bbtKIokU6wN2TtKSAyJdpu6XNdV0SnaAHdJmSiK8jxH8EeZZsmGGMznc1yRG0XR09MTdhgNywqezTTNMAzX6zUM2tPp9GyIT4yEs9k0TYv++QgOwrXsNXGhHWMeac5JjBH8As+3uisGt3mHpwc84NaNS7+Z9XoNItL74FP4WF3Xfd8Xrb0IwlYYXQ3/xHV+6Hpc9l5ZkD3Iret6HQJZlrFuh+whA8LQtm0Bj6TUHVpyXbfu0iHJA/yNOi7epXOZsPmhEWY/v/5YRcTjKIoosvTZuigqO+Xkh5QneUP5JUsgOCYalSRJ7UHcwZuNcMjXDCQ26nhloFw0drlc0r25laYa0zT3+z3A+fz8/PbtG2QD0iwwQgGoaVqappUdCoNEtVmiUpxK8FKancR6vW5396c4ndjfTgJ7BWjkKJ4dfaKWZdm27bou9gFhGLaIX5RlWZ1UMAyD6IuDdiecoNVAUpwgetfsJFhY4NaF1rEbJlybSruHLMvqgh1omsauuzEMjsdjkiTs7mFE59txYw2LEqWlvU+CZoqxJw6HQ/NF8dgby+ef1MTYLhRgsSxrPp/jlo7md3Xwa9Q0bTqdFpauYIP9F1en4cZAwzBM0xy7GQmtRtRxKN7d0yNHaG7v9ARBsNlsVqdnvV57noc+pTZS4uwgobDb5umh/OWFC/00ooQSEqJ31nw+xxwkOqP98Jem6a9fv3a7HUkIqse2bczLHYoEIs5JQACwYiAIgt1uR6ZvlMX+g+6a5RAU9iecn0jTdLVabTabKIreTk/z2VPYpmmahivB30/PZrOpPE7B9jJuXCBdDSVwOKM8PpfL5eFwYCmIjAZPntGeQrKENOomaMxH3Tvt1E2F9Ts+MN/3xbQHUgeBQ3Y6EJxh4pyfYLuD1ajwS3X7a1fqpjquyt8aWaShiXp+fkbZgroJBgmWLEi1UG+yRIZMc0KFq50E+zmLmGYXLCLy1ylPuDaroE3SdR1iZhQBoGD6o4aw/pSWZWH12ilmAxGzbTvLMtwh6jjOTS5dx96xvwbbtr1arVj65A/98PDgui57ny77YVI2lI2iaLVahWHI5mHJCphmw7AX2RtSWA1ZlzQ7CWxXh4Su87rO7iS2223hOAhswrdasXaLQJZl6/WaVTtwXDC7rboPagWbNpl5+6hreJqaptGViBzDtWVZrDuJZVnr9RrcIuDx6IYuZyehbqYbfhxeVqPcQsL3fdbfHPqEEW3SL+vL45EVFQWVxaWkbpuf1T7J1F+2bRuG4Z8erFQqcYY2CTmfn581TSMQjNNDf1YWF/ClEhICdkpTluoOvzQtL0C+wk4CfqvsdInLZMIwlGxZWoc9jjfTph73bdRlFvk9BidE+wB9l52eXgHJsoxMuJzDdMfj0fd95DQMg/YNEB7UsyO6WpUjJNRhOrZDBU1HUcROqYJyWc9WEASHw+Ht7e3Xr19RFNE1Sli1ff/+vT9jA2wDmqbFcYxQifRG07Sy6xEOWNU3peNfoiiCZxTmWRgtOq6jZ3JxHHuehz51XTcIgv4qHOxUJkZLO685doChW/sDpEPKMARWHqZTQqJDnBWpCgTiOI6iiLx4YYWez+cDiL04jp+enso8Fapmv2rDMNgTtiRaCkXKNK95k+f55vSASN9T7TWs1pVlpV1/Jyocx4njmGdireNPvT+HQBAEX19flUJC2SR63bzeNfHCxPr8/Ey78s5xwWZZ0zTWnJgkCemX2+mIC02AqrodqSZNHruin5QthmH0gVJBb9kEUpWnIQIcdZMSEg0xvFm2cblpwU+JtUVDb2tZVocIUjQOy7IKKyRd17utCG5XbIugp8LpB/KE6apSuEKhUWN0gkqShAwVbMSLTvBRQqITGCuJ3G+Av7JWEW/K7wtzjTh/ItoahZYUhzHiBAqlh4cHHI3Ge8RbxSnWTlTJFAaDVECoiA2DUVj4E4edJEgrFUVRgQfU67puV8aVNE3f398p2ofrupXngTtpVx9E4jh2HAdGoA61Z50MpD7aKwFNjk1C2sN0WMDatl0IIRDH8W63K69ABe9m8rgQh0+E9ClPl1judchnEASLxYIIwhVqPp+Xzc6Up4/E9PSQKRINh4CkkAxdCQnTNIMgcF0XQSMWi8VmsxlRSIz5fI5ThJ7nUTyPan33hV2lDBIXAvZv2fHN/turf/4RxzF7VPCfr0//V249JHg5Li0NB3A0pA8NL6dSzk+k/adhdFbvDy8mDk32pzAModuhiAi0P/B9XxwcCjz7vs9aXDgaXrZgwzQhMCKXSmoaMX+9oUKpmwjVdgkyGuHjpQ8Nf7JxbVn6yibBoiFiGkKic933pU3F9V7sDgx+SuzMWEezybe93W5932fpu64LglmW3bz5dU2re08NwVnc688QbLdbmm1HZ6ug0Ny4y6gOtLPvmwyks0RUhkoE1uv13QkJHI6vhGNcLyEkrp9lLm11nXHYsizf9y+atRFyp5KBMAzJzonljBzh8MqhOGzbXq/XV/ZjlmUEl2VZF/VCJf5DvmRFRZO1RZk3JSTKmHT1hmO4VjuJrkDui86QerMkScpBTEmtZNt2O21P5bdN62KSDdgLt5s++kL/arqAlDDsJBQHhgRoPj8/t+uUq1vWhkCSJAhi0U5vVjmQ2vChypQQ4Fw6JO1hujzPHx8f4XCS5zmZr/vzayLrJTsjVB7rLWTg/wmnGtu24ziez+d1tRCRKIo4Ft00TQt39ZDTDlFgE3BSmp8e9v1Facdx8jyPTw9cg9hD17g/Dh2Eino9tXsR511lzvP8/f2d9f66/sQZa8/HnTldcds3HXCu6/qlpngaSH1zeIf0Pc97eHiodi4oSRR5XojW04ZhkKpaNN6IH8Mw1ut1t2GUsGlgjxrU6V44iik5xiVML+SrZhiG67rX6KDCMAS1K9X9A8ObZRm+BfiqNdSbYckyMKt3Uh3H1ULanYSmaZPJZL1e0x4C82DBZZMmx3KC3X+Uf23yBit6+AIW2GhSHHmw0mc3QGy6kg5tDjhbChSk2arXK7Q8zzscDoZh7Pd7aAwq2dY0DYdCOvSsr6vo5u/ZOBa6rm+329YjRNM03MR5OBww547iUAV2V+/v7zhO4ft+9TKW6Sq1k2DA6DjJOSchrU0CN0lJsAoY0iZxc7jQWDa0xs1Z6pWBj48PfOud1ELq/nEBSFvMs2xjS9oJVopIAQHOTuL/dSyPhCG33+9p8AnDVHtG7uoM0X6/P2t6aQ+lSCVfXl7wrRJTsNzQnxclPj8/oX3a7/eTyWQspp00TeH4tN/vp9Mpv+vFV9he1GXjyFyQJ9L8yfHoGlcbETtoXDy35pb127lGU9+agZsXxKxhGEZDNX0lw6yb7Fi8xdigT9C/VTZNvewJARxirSQu7U4Cis5xCGoul9Dmc7PI8yPrF0vRxeVpXoOWbLdbXdf3+/3T01MTZ7ZKkkEQHI9H13V3u53jONPpVPxdBSKR4G7Uw+Hw9PTkOA4btr2ypeplhwgU/B6JsrRCglooQeLz81OCVlzahPsUEoh6BFsFpnjTNNvNlUEQ4OjJ4XBYLBZ8Nc6lvdNT/ul0Gp4eeDHMZrNRsN0TGgOTLQTtoNqlFRIy6S5///5NHXZXCfHXvz11B2wV0Brt9/vZbNYOiuVymSQJvgXn9MRx3BPPHZK1bRs7Kk3THMehtnuep2RGhzg3JVWphJLgJcdYP67W3dspU3bgUoy/cXVZt9yu12vyVLZtu52pho5TaJp25cmMblvHp0ZsI1AYQvDyi6hf2yHA+dakdYGVxnDdblJoN1BEKEVCQtf1sRhdB8BtvV6Tt57ruu1Ccfi+z568G8vQIju8EhL9jbR7FBKcUCT9Ad0H5SRJPj4++qAsJk1cEdouto+YLeqQK9b7qzVZEjaaprUTNq2rbl2QDfpUF6y0NXFV8Hg83qOQ4Hh0jWtM3Ju6Cb1jGIbaRlQOVASCvXKiZPclnV8yWsl2Jy/ZXeY1LsKdMCMZEY6QkNZwrWmaaZo0qkadYJd+o25Ic+YR0LB5/vvJCX9WimAxOT2e51103NLzPBxh03V9tVo9Pj6OxSAM4wR8ZNGK++n6W7VUZiHxP//zP7eCtdt6V6tVtwTFp/by8iLNSZde0UYojs1m8/r6eqmrtOd5P3/+BAXHcUYhJ+bzeRiGCLqz2WxmsxnJy15xvgfidUjKLCTu1nNUggH98vIyijnr5lB/fn5ut1vLsvb7/bdv3y7dUpimScE8HMcxTXMUPrKmadL9S6vVqvWpw5t33ygYkFlIjKIDmjDJnkNukl+OPLquj2LCujna8/k8jmPcb7rZbJ6enuhgQUPecBsge9L7IuVVw1q6zQa1G4J54NSh4zhqwLQGmXOwTAmJ1qgOV/A+19Tz+VxpnJoPMuhhcAZtsVhcGuAPJ73DMLQsa7fbPT4+LpdL8UUFBfOwLCuKohYCsjnCd5tTCYm77XrRG24Yxn1Kx2s6BrsKxFdYLBaX+m7AXwDFR2TQnk6ncRwjlgkEpBo5l44iznUmSkhcCubQ+dfrNcx0Q1d86/pM0xR/JXtrkKrrRzQO27YRM/xSQ8Vyudxut/CpQzAPcTqCDp+XW/7y8oL4gAjm4ThOOY96U4cAR6wqIVEHmijvp9PppetBUVi/jg/btjmrm+toy1/aNM0wDDHXbzabx8fHi/T18/k8TVOYxKMoenx8vFTS9AQxXwOJ+IAwVERRNJ1OR6E06wmrrsgqIdEVkn3RieO4zjWtryqFocseMBaGqTExQnM9Lg+/lHUor2CogEl8FHMuGSpwCgSXnl7adpWfEKgQEpx9BxVTicEQ2Gw2d3hOYjB476Ei2LSxLVsul5fGHidDxeFwWK1Wr6+v7UKXXw+1cXoa0oGhwvd92OFvxXNDbm+ejXNit0JI3KfD5c07STGgEBgAga+vL8Qev3RPQIYKFPc8bwBuC1WkaXrprhps67o+m82Uj2wBT/bPh4cH9k82PTkej+zf0qSjKHIcR4LW4SisBA1pMbTyPH9/fy9cmIVgVi2otStCG2vLslgbCb2/ZlGVpmmhdQh02o7VhqWiKPI8D8r9MAwv5T+OY8dxUHy9Xt9EWjRsKWXL83xzegDver1mu5KysQl0zVlw8jzf7XZns7GUxUxHUVQrgyULU0XNgTqb/hxvQpqGtOiCSpsEIpE9Pz8XQpIZhoGwd3WxHeEkdlHcwMIenJrAvjcMg94j4ft+4SVC2xay1Qn+ASKzskFVn5+fy4ydfcNOiwMwDH6uDO388fFBHcePkEitMwyDE4OZvcqtAALus+HXchbkITPUDdHj8VihbhJT0N0tV+J4H96qCwrXHkBylGOusD5g5eW5pmmVLzmNMk/PdrvNsgwrbvRFEAT7/X69Xodh6Lrufr8veFv++vVrv9+z3kT701NZF047485O+COxDakscv1LNhrH79+/p9MpbYwaEgfDOKaLAEoDDNQ8zy/tRLY5Ly8vCGuoaRqCebB9RDmDIIjj2Pd913W/vr6+fftW2bTlcrlarTAMLMt6fX0FBc/zJpOJ4zir01MYG1SLaImvr69aloYUVkPWJc0CfL1ea5q0d0PxhwSnE33fL+wkaJlZVwrvG+4ksiyru9PCsqz1ek2cb7db13XpT9SCC+DoZV28d03TWFKUf7BElmV0E5FlWYXlcBM22Ivz+m5LV5cOodWYEw3DKKxCdF2nhmAY0J8sIIUutiwL0ctt267Mz5YVMM25pU3tJGrFp/pBEASCIIj++VSu6cAn5yd+Q/I8B/k4jpHQNG06nRqGgUVlofhut2M12vP5nI2VtFwuMU0U1uak6ChQK2Qr/Nr3nzhJADm32+1g3b2oUs/zttstJA3OeLNoXERqsMxoNU7e7ff7x8dH1k32cDhQ/06nU/hHVfI2n8/pPRtFhn+Yg4qMJiGgTOuEJSigOyF1WyL3djMdizatyulzIl1/eVHZbidRroI2EAUVNq2yORs7drOiaRqps23brlT9U7sowTZ/4DRB0a7eJElIEBY2ee0IlkvVbcjKOZu/oV4m+wENABBh96y+77PDgN2VQv+GK96oN2m4NufnVjk5Non/oPZIluB4dI2rpQNoqAUHxLIscqGhqAw0H13JPEwCBSJQtUOFjb0F9Mvb7RaLxziO2VUkFYeOG2ZP27a/fftG1um///6bshUSZCal1hUyDPOnbdtJkrRW+pum+fPnz/fTE0XRfD73PI+a1kkT8jynNX4nBDVNe3l5+fnz52az+f79O9GsDIkaRdFqtXp4eDBNk7P/Qyfatg3/KMdxSPoSfQETX19ftXPmrQRX3/WiY/quZQD6iKwwQEUCVoFOrLyo0vd9y7JYni3LIu+mysU+qLGrP7b42bTrulggF1aax+ORlNqYICzLoskR6mn8Wa6ivB8q57nJGwQeb6dbpwBKaB2Bc31DAOz1dPgUNE1jhxy7k2AHT2EYwIJdpqzrOluqnEGQN7QRL/MjrU2CvlIBhfZFLMVx3Hpxd1FFwmauXLObprnb7cgOAXd17Lrq1uN17+saXjijS8Uty2K9YhDaSNO0IAgOhwOFnLJt2zCMzWZTR1/w92maLhaLFo5PCKCExQ3FfeqksZ1vIyq5siyLNSrs93vq+sKsQsNP07TdbldJDZFRKn8S6mXl5gkcSqtuQvPGstfjDJdhPgwOAzf/KYqiwseJJaphGK+vr/gpiiLDMJDGF144mrtcLrGbXi6X7OxfyMY2djabEU04TcIv/vv376vV6nA4QOW1Wq1QLw5VsbqFz8/Pb9++UXVsXaZpEueO47DaMzYby8+QaUxtOMjpOI7v+5dyBQo404qDbEmSjEJ3CrM2eme5XEZRxPYp9YLv++/v79+/f0eYkyiKMDwKp/DSNGU7l4qLlphOpzRQi7yVNxfSvCk29fQ3dqy2bbNuiyI3GQNUZA575S0Mwzp9RZZltPwpbOpp6YcxQPbDwueq6zqH+SRJiH5BL+S6Lo0uy7LAYYEHUHZdF8csKD8SGH7saSzKwOo6OOwN8xPNj7Zt13UEn5MwDKk72umviP4w6qYsy4hhTdPKPrLED7JhUJGtHojBxm5ZFn+MEambJ3D6r5INacNyaJo2mUyyLCusxOM4xkpT1/VKPQZ9q4IkpIkv0hOeMCFi3VeoAtqAwgAo5Lnmz8otzjUEhS3red5ms9F1/e3t7dItBRrFxgJpHcwDB9NIbvUHF7SXmqadnSXIgs2OQM/z4jjGzrJJCJD+GlKgHMdxnf4zjuO6zpVWSKRpOpvNxrLDLfQl+6cSEiwaKn0rBDAOsbKu1UucY840TRjYDMNoQWRANxjTAAAgAElEQVQwIXGuHWP9HQGaEH3AMIz9fk9769VqVatUrNxfSPASRyUlaMidq5sk6EFpmgCXJ6yvW3vskPaJ9DPN8dF1vUWp5vTvOec9ejdBycC6H4xU+rOONCNtgmJbDgRwLwVWLY7jtLNC27ad5/l6vY6iaDKZ0AkYOSCSshXSejdBV9ifPnqw0eC67ihsJz0BslwuyQPYMIzv378PgwbpmmGy7rx1URQVYo93XkVPBG3btixrs9mwfqKX1uV5HpyINptNFEVBELA6/TpqtAupy6Det0aAtdUXici6w8J6p/WmWFZYRteu4njVNAqfQG0h5yUKscCWwhhA0Agqwk8UJiw2XnTZH4lqPx6PrE2V5TMMQwrR8fHxQeyxsTo+Pj4o5APYG8voZVvKB7b8K0H9/PxcaH45s3rTHwIc7yZpw4viUoF2Tnv99UQLyq7rjmWyaNG6s0XI9zTLsiRJcFczWwqzNjxKkYHmHZSFRymcVhtqtKH7pvDdcIRFpTS/w4TLRhOChCDfR2RAKVL4Io9lWQgzTq1D9A5d19kRS6XY9oqW3m63wKT1KKU4smz4VdGaKT0/9ygkdF2XY2Eyipmiv0+InUZRC2ZYqhEzVMEbHXKC8hyPR9pNs7Mwm4FNF6rADI5JkBPcuzDH4QoKkNX+uQEqNAdhKMESGsLGGhlL19OpEUR/YpFsmC7E7q6TNxQZpSFZla05Ali+VOaXNiyHruvtDGvsUlGlBURgOp2SPwIsB9vt9nA4sFaEAtv4CQ5vdX7ihSLsnyhOYqbS0JWmKWJyUEE2dCB7KI9VSZmmaVkWG9Fht9uJH2qb2ohEEARJkriuG0XRbDZrYYum2N10idN8PqdepuqusYIQEZW4FAFpDdeVX/Kl6Kj8YiJAs2oURc/Pz/P53Pd9ztE2BJbA3QBkBuc0DXMWOfKvVivUgiLkbK5pGhtdA7dQFPLYtm2aJmc0sj9hA/T+/t5inuU0Z4CfTNMMguDh4WG1Wm02m1+/fhF6zWufTqdBEPzXf/3XH3/8sdvtHh8fPz4+Xl5emlNQOftAQFoh0QdYiqZQCOCyIGhsEE+pkj3sAxAI2rbt2WzWJOL0fr9frVYgaFnW5+cnES9sWSgYOGXQNM3zPCx7cQcqAvsUCrL5Kb1er5+enkYacwwXLv348WO3283n87/++qvFbh6xuxeLRRRF3759Y8Mp5nk+jG8bdcddJehgXaHVMqub2sUPKACk/hQQAWhvMOeuVivHcaBEqrxPGNobz/Nw7bCmaU00ThTaa71eF7QcsFiEYZgkCVRYmqZhsJGGJM9zVnhAVcUqnepQxWoaXqFErS6zgO9N04zj2HXddvfcoUVsHFncAxgEgeM4kD0CtloClvI8Zz0+2BZJKyQKCmK2zSo9XgRwlSmWk4gmRG3RdT2KosLZQwThMQyDtDrw8adSdQlasULzw4of3KVTqUdi9wrgBOIB1KbTqa7rLKk0TXFggmUDxozFYkFaNfbXUaSDIMiyzDCMKIpYG9JFzCOOLO54x8biouIq80UI7Ha72lV1pTlbgpccj65xtY51sR8X551wWx7o8MrHaQPWDQZ6J/xK3k1wkKUzCmCpcF1MmU9M0/QeRMi7qe5YACrFr2CGjkEYhoH3IIU0LNh0zKLg+MRSI05Gl2Chu4Z50oQQXNdQU2XLCHC8w6U9JzEW98Fyb6k3LAK4FwyHD2zbpljTvu8X3F6Px6Pv+5jK1+s1ZpNKp8n1es0Px13IgKvW4KhKsxUrvfATG7ocMz45YbuuS6IF8yaKs2GocXcFtR0Bq6kUvb/bBOFmGAYNg7tFo/OG36OQYD/LzgEdkmBhVTtk1aouDgK4JSIMw4K8SZKEcwcGEazMUzjDUfiTyo4xgSnedV0SnC1awQpmXdcLyLcgqIoQAvcoJPzTQxCMN2HbNmktxtuKdpwnSVLQFLWjo0qJgABtwlgl4UWM0ZlEMrGqYB4XAcjJzFG9SGu4xl2V7NJjvOm///57vMxfw/nDw8M//vGPayiosuIg8PPnT0zujuO0PjAIn9owDBHM4/fv37PZrDU1ccARmRNphUSLC9yF7SdaNwnLYU+MTafTFo72PTGjyF6JABxbcSnsYrGA81Jrmp7n5XkO09RisZDpe2+NSU8FpRUSOA3bE2qKrEJAIdAOARykWK/XOEhx0dlyXBnN1rtcLvM8d133cDg4jqNEBQtOV2mZhQTrtN4VXoqOQkAhcD0CnufBA3iz2TTfLMZxTOdXWB6CIAA1iAoVyYMF5/q0zELienQEodAk3JAgrCo2FAINEUDIWNu29/t9wx2A4zh1p9DhnIOTd79//25IsCGrd55NCQnRB8B6vSYPcdF5VfwpBC5BwDRNmKAPh0OTgIBRFPFPoXueh6jjSvt0ST/8X966LZ0SEi3AHLpIXecNzYeqTyHQAwKe5x2PR8SEiKIIZoZyPU2kCALxLpfLLMvIUFGpoSrTV2/qEJZZSMixAN/tdkrHqj7gO0EgiqLVavX6+lqesC5yBEecRBgqdrudaZplgncC6fXNlFlIUKjn62G6IYU4jn///n1DBlTVCoHBEHBd1zCM/X7/+vpaMD+0WCrBUAGzx+vra90eZbDWjbQimYXESLtEsa0QuFsEcHgCO4DHx0d2Wm+9FSCzx2q1enp6anLyDlqvJjmpp+I4Rjj6QhxiyjDiBOeg9qh/wjgbdRPAPJRmEjRENUEh0ByBJEkQxsMwDIR7wofQOqQH4vJipubEKUKYSJrQOUFn2cA/LHFN066JT9Ucom5z3mNYDurmsSdkii8y9r5Q/A+GgGmaP3/+9H0fqqeLFvWVTJqmeTweEcwjiiLct1rOmef5arVCVM31er3f79n7Pyj/dDpdrVbkm/7jxw/btrfbbRiGtm2/vr5SzhElaiM7dCuOxKEmzU4Crt/iAKs4UQgMiQBCptNU23onQTwj8DsI6rpeCMYehqFlWZQZX18hHC/ysKHdNU2jPLiskCiMJcHZXSmbBA0/QROFuzMF5VKxpRDoBwFcbOe6blfkKYQUfGRXq9V0OiVDQuH+c1xoWDicYdt2HMe07kZkB7r6kBJdMXxzOkpI3LwLzjDw9vaGW9jO5FM/KwTkRSAIgo+PD9/3Kd74lW2Fugn6hsPh8PT0hIMaBXnA1uI4zufnJ/sG6cpAxa3N7GX6N3+jhMTNu+AMA6ZptnD+O0NU/awQGBsCLy8vy+Wy24Ol0LFgm/L19cWBJI7jKIoqz2pUSg4+NU5FAv6khISAnfJvLAVB0O2H8W/U1R93g0Dh2MHdtPt8Q4MgOB6PsI2TEqlQbD6fH4/Hyl+xBSnnL7wR/8/KhmiapoSE6H232+3IiUJ0XhV/AiMgn668J7BZaYp0pWwo1E7+V2zxQp6R/vkfI+X7fti2bVvFPL+f7lYtvS0Ctm07jhNFkW3baZouFovn5+ezLNm2vVgsoGJC2bNFRpRBCQnRO0u+hYnoiCv+7huBMAzpbAQOTFTiYZomGR7W63Ucx4gDpOs6HGcrS43xpRISovfabrfryqND9KYq/hQCAiBg23aWZbvdTtd1TgRZ+/SAX/jpxnF8OBya6KYEaOUFLCghcQFYN8lq27YSEjdBXlV6twhMT8+lzedIlEtJCZVfGa6F6o4KZpRBogIU9UohoBAYCgElJIZCum09uq4rOdEWPFVOIaAQuBYBJSSuRbDv8q7r4lxo3xUp+qNGIM9ztZgYdQ8Ky7wSEsJ2zf8xtlgsKs95is634m9YBL6+vv78808lJ4ZFXZ7aOPZ2JSRG0M2V5/5HwLdicUAETNP8/PzkfOoD8nIXVaVpSmEBJWgwpy3SejfleS6HU5A0DZHgQ1JNUAgQApIFy9F1nY59UBuRkHYnsdvt5FhVqWgKhSGr/lQIKASGREBaITEkiKouhYBCQCEwagTqthGSB/hTcfFGPWoV8woBhcBgCNypkJDDJjHYKFEVKQQUAneLAOcGTGkN19J0tm3blVdfSdNA1RCFgEJABATUfRIi9EJLHtQ5iZbAqWIKAYXA1Qgow/XVEPZMIIoime7L7RktRV4hoBDoGAElJDoGVJFTCCgEFAIyIaCEhEy9qdqiEFAIKAQ6RkAJif/f3tmjOcpje1x+nhs76zfCvQCo6I2gEruj6g2IDdRE0xF4AcYLMI7ebLwB2MBU1FWJYQGGBRiiqaw24PtcnzsaDQb8BTbIf4JuWUhHRz9RHPR11DBQiAMBEAABlQjASKjUmqgLCIAACDRMAEaiYaAQBwIgAAIqEYCRUKk1URcQAAEQuIRAjZN5lY1ETbUvoXinPI7jqOGp8E78UCwIgMBVBFQ2Emr48o2iCE6ornrGkRkEQOAYAdM0q5KobCTU6EnEcaxGRaoeQcSDAAjcnUDNkQQqG4m7c29KARiJpkhCDgiAwLkElDUSmqYp4wVWmYqc+3QiPQiAwN0JKGskGGNqzEkwxjAncfe/EygAAg9LQFkjkee5Mh/g4/H4YR9QVBwEQOC+BJQ9T8I0zel0OhwOCxMypmkWYu7bAEdLV2nc7GhlkQAEQKBrBJQ1EgQ63l8ydNd18zw3TTOKIjm+s2HO+XA47Kx6UAwEQEBtAoPdbqdkDV3XHQ6HVWct9ajKtm0zxoIg6JHOUBUEQKCDBML9VapYGIZVtkDZnkSe5//6179KcfQuUpnJld6Rh8IgoBgBsZ5e13Vy5ZCmqYgsrayyRoIx9vHxUVpnRIIACIDAAxLgnFd1F2jEopSJsqubOOf4AC9tckSCAAiAwCGBqv6EskaiFy6PoiiqahjDMFzXPWxIxIAACIBAGwSSJCkVq6yRsCyrtMKdiszzXDYSnueJdkrTNM/zTmkLZUAABBQmULXMR1kjEUVR7zxsz+dzbK5W+I8QVQOBLhMQX6gFJZWduO7djuuqfRuO42CfROGpxU8QAIHGCVT5MVK2J9G7bkTV+NhyuZSHpBp/MiAQBEAABGoIKGskaurczVvUk8iyrJvqQSsQAIHHJKDscFNfmjOKIlqhTEZiur/6ojz0BAEQUJ6Aykai+5PAYkwsyzJN0+QRpyRJXl9flX/+UEEQAIGOE1DZSPTCLYewE4UHpSq+kAw/QQAEQKARAt7+OhSl8pzEZDI5rHCnYmisKQzD+qmI7neJOkUVyoAACDRIQOWeRIOYWhKl7S/P89I0NU3z5eXFMIzDPgRWN7XEH2JBAAQEgcM3D91S1lW4bdv9cq/t+34cx1EU5XnO/32J9kMABEAABNojYNu2ruulm65VHm5qD2gbkl3XDYIgyzLOeRiGtm1TB8IwjPf39zZKhEwQAAEQOEoAw01HEd0oQRiGSZK8vb3FcUxDT6ZpMsbSNH1/f+/+/MqNMKEYEACBdghU7biGkWiH98lS6awoeZQpCAL5FG5N004WhoQgAAIgcCGBJElKpyVgJC4E2kg2GlZijM1ms9LRQMaYvHmikUIhBARAAAQKBJIk+fvf/16IpJ8qG4kkSao6UKUsbh9Zc1CUUCbLMpyeJGggAAIg0AaBNE2rxrRVNhJhGHbcSERRtFwuC00uFryu12vLsuI4fnl5KaTBTxAAARBokEDpQBPJV9lINEiwJVE0lFQ4XEi0Ft3FnERL8IXYLMtWq1XVcJ9IhgAIPCYBlY1EVe+pUy19dMoBh0m03V5xHM/nc8YY7ETbqCG/swTEAMahhiobiZ6uHCV3sMJ4VB0XddiWiLmMAOc8CAL02C6jh1zKE8Bmujs3MfkJJyWSJBmNRs/7S8ymuK4rjHySJPVenu5cmd4WzzkXVrm3lYDiINAKARiJVrCeKLTg2s+2bcuydrtdEARpmpJtyPNcdCaenp6m0+mJwpEMBEAABK4nACNxPcPLJYRhKFYuJUmSpinNWnPOx+MxGQnRjbi8GOQEARAAgUsJqDwncSmTm+YTw0pvb2+MMbG0aTKZkIdwEUNqiZPsDrXsl0PDQ/0RAwIg0EECKvckxPu3g9xJJdM0fd9njGVZtlwuZ7MZxWdZNp/PC+aBbhXWy8pVw3SFTAPhzhIojLJ2Vk8oRgRU7kmUvmQ71fAvLy/T6XQwGJBWYgnmarVijJVutKalOJ2qBZQBgbMIdP8P86zqKJ9YZSPR/cYzDGO9XtOma8dxhMKvr6/G/hIxCIAACIBAewTI53SpfJWNRPd9N5H/vsLiyzAMNU0TX1ucc9GlwKxD6UOMSBAAgSsJxHFcJUHlk+mqDlqqYnH7eJqKoGkJxliSJD9//qRZB13XxcpXoViSJKvVqmpaAiZEgEIABEDgLAKDQaUtULkncRajuySO45gWNVHpv379Gg6Hm80mDMP5fB6GIedcPlbw8/Pz0CEg5RW9jbtUBIWCAAioSkDl1U3db7Moil5fX0nPLMs+Pj5oNoJmsMk/B/1LaSaTya7iOux2dL/60BAEQKD7BGAk7tlGeZ6LQ+hoTFBMRei6TsNKhRmLe6qLskEABB6PgMpGQrxwO9uspmmKDdW+73POyWZkWZam6eF6AzrJzt5fhmEMpGs0Gsl9js5WGYqBAAh0kEDN21LlOYnuHzpE+yQmk8nn52eapuItT/skRCdDPFLfvn0TRmU8HovNd5QAfQ4BCgEQAIGzCNSMV6tsJM5idJfEhmEEQRCG4R9//CF20jHGDMOYzWaHtn0ymWy32+X++vj4+PPPP8XKqLvoj0JBAATUIFBzbo2yw019OR6AdlAHQVAwCXLHQn4KR6OR7/vb7TYIgiiKBoOB67pyAoRBAARA4FwCVQvrGWPK9iQ0TTscrjkX3I3T0zaIMAzzPOeck50r3f0w2l+MMdu2wzCkNVE31hbFgQAIPAIBZXsSX19f375960sT0ubwp6en5XKZ5/lsNguCgIycmISQ6xKG4WQysW2bc55lWfddGcrKIwwCINA1AjUzmsoaiTRN39/fu9YSh/r4vm/b9tPTE81PbLfbwmHLYRjK0xVRFFmWZdv25+fner0u7WccloIYEAABELiMgLLDTaZpfn19XQbllrnopDnHcU6ZgqbBJVIvTdPn52ehqqZpcBUuaCAAAiBwFgGxtPIwl7JG4uvrazKZHFa4azGbzYZGmd7e3v7666/6ITJ9fxWmuLtWI+gDAiDQOwI1q5tUNhLv7+/dtxOGYex2uyiKlsvljx8/6NlyXbe0Y+Htr949f1AYBECg4wRq5jWVNRIdb5KCetb+YoxFUeS67nK5fHt78zyPcx5FkXDe53kenWlayE7exdfrde8WdBUqgp8gAAJdI6DsxHVf9klkWUZuNqjrYFlWFEXr9VqMKckVEZGFx4hzXrM4oZAYP+9LIEmS0m7ifbVC6Q9OoGZGU9meRC9mrRlj379/p6czDMM///yTxsdEx4Ix9vLyIgbNaAXUgz/Nfa9+GIal3cG+1wv695pAzSCEsj2JXvwdhmFomuZms9ntdpxzMSchP22e5wkjIccXwjWLEwop8RMEQAAETiegrJE4HcEdU4Zh6LouTRktFotSTU7xuuF53inJSuUjEgRAAARqhptUNhI18/UdfCbiOBYT1LJ6VUfRyWkYYzW9xUJK/AQBEACBAoEa303KGgk6+7MAouM/ezFE1nGGUA8EQOACAjVGQtmJ6wsw3SWLWOhCPppc16XW0nVdnGx6F8VQKAiAwOMQkFdRFmqtrJGoqXMBwR1/cs7D/UU6cM7JQojxQTrTlO5mWbZarUp7G+Qf8I4VuUHRtm0n++sGZbVaxHw+h8etVglD+AUEal6YKhiJJEmenp5ms5nsCK8X+wb4/qppUXmm4evraz6fVyWu2kJRlb6P8aUGso8VUUPnXng0UAP1DWohv2oKxalgJOiwTyXfIPLCVnLgUWi/h/pZOK71oeretcrSceu73a5rit1enzAMsyxTeHmhCkbC9/2aWZfbPzQNlmhZ1nq9blAgRIEACDRLIEkSJb9QBaWeGQk6u40xVnBEoet6oZ3EsL6oak8DvRg3uw3bQhPfplCUAgIPTqDrS2CTJBEzDZ7nkVft5XL5/Pxcf6ZQX9xy1D9/YRj2a7dHfXWuvFt6SN+VMpEdBEAgSZIaCJ02ElmW/fz5k2ZraW2Pruvb/WWa5q9fv+SKycP3jDE1jARjDJ/PcisrE4bBU6YpFahIzWESjLFOG4npdCoWZsVxnOd5kiSj/eX7vvz2NAzDcRy5tZSZpTBNU64XwmoQKHzTqFEp1KKnBOrflp2ek6Dl5IPBgNDLqzxppJ4OfKYpip42z1G1a5amHc2rWAKVVjeJTZSKtRGq01MCNV8tne5JHMVdbwCPZu9Fgvrhwl5UAUqCAAh0mUC+v6o07HRPoqD0WQuWTNMUXRBZjmmaN/421zTtmhVKnuedO359ZYkyrk6F5QHGTikGZUCg1wTqx7T7ZCRkHxXUOaqpW1XvKYqiG/c/siw79y1/5QOXZZnM6kppJ2aXBwNFFqo4TSwNh8Pr3/KDwaC0IFFiVUBmUi+BvkUIoKZpeZ7Xp68qsTSehPu+X/9IiBKTJBkOh3Ec08dNlmWj0ag+b2m5RyPP/Xiiv69rNDm3xKNVuFcCwzCuf7Dvpfwp5Q66v2dyMPh/JQeDwe/fv+kEHs/z5vN595U/pQ0USJMkyenvC8755+enWMFsGEaSJJxzwzDCMNR13TCM9/f3b9++lS7/FUuiW+JGhcqvadu26S1A23GqpkYadDxc43/lrFoX9g+Nx+OPjw/S//39/Y8//hDe6RU+L0/XddGaZ9GTE9Pj/e3bt8lkIrMqEKYsTZVIxTHG5BJlrUT43BJJIOVK05RqV/U67ZORsCzr6+vLMAxN05bLJef8ekdpdDbcjQegRNMiAAIdIZBlWZ7n14yLUkVO9zV5sz693IlshDaZHOrklb462iiR9nVXWbsrS9Q0rbA9+b9A7Tp/LRYL0nG73XLOtf3lOM71ipONCYLgelGQAAJXEpjNZjQuhwfySpLI3iyBHqxuEp6zRqNREATZ/qpaQRhFkVgOlCSJZVmD/eV53uG8d5IkYh/Gf1lO/ACB2xKIomg+n/u+v16vxQN/WxVQGgj8h4A8k9oDI/EfxY+FPM97fn4Wg+NPT09xHM9mM875fD4nZ7GyDM/zru9fywIRBoFTCIRh6Hme+JphjNFiCuryW5YlJmxOkYY0INA4gZ8/f9q2TWLVMRJZlsnTfWQqFouF53lBEHDOV6uV7/vUsShdHds4aAgEAcaY67pi5Jqccdm2PZ/Pn56exCR8lmXyfDiMBJ6cmxGgToP4tqZy5TUj6hiJPM8dxxGLYmmJnui5B0GQ57nrurvdbrPZbLfbm7UBCnpkAq7rLpdLseqaHsjtdrvb7RzHmc/n9KCKBI/MCnW/DQHP80QvwfO879+/2/trMBgIUyGvCVLHSFiW5fu++GSr+aszDEMku02roJRHJiA+XGhYKQgCevxoXo0eVMuyxN9nlmXyd9wjo0PdGyfguu58PqeHLQzD+XwuVgZxzoXxkMvt02Y6We+mwkcPEG2qIMh5TAK0b47+9ugvs9QA0HJ7ShbHcdVKx8dkiFo3SGC5XBakycMtpePw6vQkCjU/cdkS/hoL3PCzPQJyl0KUIiwHdfCzLJN7+iIZAiDQCAFaHVsjSvRoRRplexKWZS2XS3JjwBgTM4Si5giAwG0I1H+viLvo1N6mOVCKIHDou0j+jlksFpRS2Z4E1XY6nTLGoihq0GWCQIwACJxFgKYi5JWvjDEswj6LIRI3RSDLMvGBUipTDEMp25MYjUa/f//+8eOHGGVDL770UUBk2wQKayjEkbQFa9G2GpAPAjIBOr1NjmGMxXEs3HmJW6oZicViIQ4unUwm6/WavDyR2zhRbQRA4C4EaNnr19cX+R/TdV3u4N9FJRT6mATE3IMYZaGvlsOFFT1w8PeYTYhaK0Mgy7LlcikcyYT7izGmaZqIVKayqEgvCAwGA9pibBjG19eX4zh5ni+XS8dxDp9JGIletCmUBAEQAIHGCBiGEQQBdRps26ZexXg8Lt3qDyPRGHcIAgEQAIE+EiDnp1VbjGEk+tim0BkEQAAEbkRA2SWwN+KHYkAABEBAaQIwEko3LyoHAiAAAtcRgJG4jh9ygwAIgIDSBGAklG5eVA4EQAAEriMAI3EdP+QGARAAAaUJwEgo3byoHAiAAAhcRwBG4jp+yA0CIAACShOAkVC6eftZuSzL6DzFfqoPrUFAKQLYTKdUc6pRGcMw0jRljO12OzVqhFqAQH8JoCfR37ZTVnOyEFXVE65mqhIgHgRAoEECMBINwuyiKOFztIvKVei02+02m03pTd/3oyiCe+1SOIgEgTYIwEi0QbUTMj3PG41GNLg/GAwOjyo8qmVLBmY0Gg32l+u6Zx28kyQJeTOu8kR2tEZIAAIgcC4BzEmcS6wf6Q/HZDabzeFxIlWVCcPQtm26q+u68Cpclf6seM/z5vO5yFKl2GBQfDgnk8nn5+dZdkWUggAIgMBlBNCTuIxbp3O5rksO4oMg2P37OmohoigS5424rus4ThAEdNLfarVqsMKe5wX7i47Y/fnz5ynCfd//+Ph4fX09JTHSgAAINEbg3+8Q/K8IATrK2zTNzWZzVpU0TVssFodZ6ASrw/hGYjjnjDHHcQ6lMcbW6zXFb7dbxthsNiskk61g4RZ+ggAINEKg2KNvzPhA0J0IjEajPM+32+1ZA/c0vlSai8adyPa0USfLsuI4Xq/XlmXJ8uXhJrpF0yqHI2lkxuS8CIMACDRGoBFTAyEdIfD79+/SL+6j6o3HY13XD5ORbTi3U3IopyaGFjJxzgtpRC+BEggdGGOlqhay4ycIgEAjBP5vvxIuZQjQ6M12uz2rRjSYczjWtNlsqsagzpJ/NDFNTlQlM01TNiGcc2E/qrIgHgRAoCkCmLhurE92d0FZloVhaJrmWQNNjLHlcskYIwMjahFF0dPTk2VZruuKyJYCjuMwxjzPO5Rv23Ycx4vFQr7lunlz5+8AAAl0SURBVC4t7T13Ea0sBGEQAIGTCDRlbSDn7gRoaOhwdveoYrquj8fjQjJd18kxRiG+jZ9i61xBeGmNCsaMMSZGogrZ8RMEQOB6AuhJnGRKe5GIlr2+vLycq22apr9+/ZJzua6bpinnnPbTXbART5Z2NGwYBr36CwW5rqvremHZK63NpUefehjNrtA9qi0SgMBDEcDqJnWaezAYHHWKR4bk6JAULZESaDRNy7JM/Gwj4Pv+dDp1HEfeqxGG4T//+c/6HR6e561Wq7bVa6PKkAkCvSDwP73QEkoeJUDf4LPZrCqlvIlapKl6+5/4zvU8L01TMjxiVoPvL1HEiQF5kStlWS6Xp+z0DsOwsHb2xBKRDARA4CQC149YQUIXCNDw/e/fv0uVIeOh6/psf8lPRmn6UyKrDJK8EukUOSINaSV+VgVms5mu61RTGqS6YBqmSjjiQQAECgSwBLYApK8/aYFQqfa0wlXTNHlpLBmVi6emt9strVsVC2c3mw3NEFxsJGiqXFaytDpBEFDRZFRgIUopIRIEmiKA4Sb5q7rH4TzPqxxo0wpXx3HkpbGHa4TOqnwcx3mec87FAlljf8VxfJYcOTHnfD6fx3Es6yknoDANZ504uXKYHTEgAAJnEcDqprNwdTdxlmVV79Y8zxljVXcvq9LRKZALxJbuk6iSQ6ai2UpVlYV4EHhkAjASirT+NZ/wFyAgw3O47uh6P95iGvwCrZAFBECgcQIwEo0jvZvAqhEkGoYqbEG4Uksq6/CFPhwOr5SM7CAAAp0iACPRqea4UJn6gRrakrZcLsU7nRx4XFjYPptpmpqmFRyyhmFIPYyLJVdNq1wsEBlBAASuJICJ6ysB9iP7X3/99ePHD3HY3PVKj0aj19fX+Xxu23ZVr+KCUkaj0fUDVheUiywgAAJVBNCTqCLTv/iaz/DJZLLdbsV4lKZph4mzLCPHeaLDUY/A87zZbKZpWs1R2EmSWJY1GAxGo9GJYtM0rS8Xd0EABG5JAEbilrTbLat+qc9oNBIetrMsK0xRJEny/PwchmGSJLZtn7iH2fO8LMt2u504Qk62PeRHlsqlsSnhcqNdEJAOAiDQHAEMNzXHstuSfN+XrQh914t3OhkG2mEXRRFNNoiex9GalRoV2p9BMjnn5PS7XpS8S64+Je6CAAjchgCMxG0437+U6XR6qMQ//vEPikzTVMx+W/vrMPEpMfJWuCRJxuOxyKXr+nw+Fz9LA5ZlkWkpvYtIEACB2xOAkbg98/uUGASBPCugadrr6yttdDjRnd8penPORW9gOBx+fn7KucQtORJhEACBLhOAkehy6zSpW41zVhqGkk1FFEWnjzXJWgqXUOQUdjqdRlFkWVaWZfP5nBxMyekRBgEQ6DgBGImON9CN1NtsNk9PT3Ec67r+9vZmWdZlRkJW13Xd1Wr1/PxMkbquw0jIfBAGgV4QgJHoRTO1rqRhGEEQ0EYKXdflDsE1ZSdJQiumqo6wvkY48oIACNyAAIzEDSD3owhy8d24rjXDXI2XBYEgAAKNE4CRaBxpJwRmWda2y78oisgJh3wYqqZpYjkszUYQjizLvr6+DMMIw1DTNNl7h5y9E+ygBAiAgEQARkKCoVBwub9uUyF50dRlJa7XazIt8uT5ZaKQCwRAoFkCMBLN8ryntCRJhO9uf3/dTBuaeygUZxjGubPfbfd+ChriJwiAwFECMBJHEfUgAdmGr6+vNnSNokjTNPH6Ln3v07F0jZSOvRSNYIQQEGiKAIxEUyTvKYc8p4rJgAZVcV23sAU6CIJSO9FIoVEUYZlsIyQhBASaIgAj0RRJNeXQQdZNrYg9yijPc9nB1NH0SAACINA2AXiBbZvw7eS3NOsrnAC2XROcJNE2YcgHgQsIwEhcAK1zWcg3n5g2aFa/5XJp7y/f91uyQ6QwrZJqbyyrWSyQBgIPQgBG4kEauqSanueVrl7NskycYUebIehYoel0+vz83J6dwHFDJY2EKBC4NwEYiXu3QEPlc85L3/g14ufzeWmW1WolDpsbDoe6rm82m91uFwRBnuer1apG5jW3SpW5RiDyggAIXE8AE9fXM+yEhMI25mt0EpstGGPyPAE52Gj1e/9m8x/X8EFeEHgoAuhJKNLclmXd5ku8pX0MdJxqG6t4FWlgVAME7kQARuJO4Jsulr7BxelyJ4o/OsFg2/ZoNKJkNAbV0nucLFxLwk+kgWQgAAKHBDDcdMiklzGj0UjX9ff397O0P7opgaY6vn//LsS2tPro7e2NMabruigIARAAgS4QQE+iC63QjA4vLy8fHx9HOwdyYUdHqDjnYos153y73crZGwynaarrujwd0qBwiAIBELiYAIzExeg6l5G+8QteNOq1PGWCgewErW462vOoL67qLtmqc8fKqqQhHgRAoEECMBINwryzKBrQp3GbE1WRz3U4MUsbychItDSQ1YbCkAkCj0MARkKptg6CIE1TWinUo4rBr1+PGguqPhoBGAmlWpxzrmnaWSNOVXMYp4xEFdiFYTgYDMRu7cLdqp++7+d5DuevVXwQDwL3JQAjcV/+zZceBIHYL10vnXy7TqdT0fPIsiwMQ9d15VNI64XId23b1jTtrFGjJEmWy6XjOC3NdsjqIQwCIHAJgR0u5QhwznVdP6Vas9ms9KEJgoCyk08OMiez2Wy32202m/F4TLnG4zF57NjtdiSK0pxSNKUhUaenR0oQAIEbE2A3Lg/F3YDAer2ml7jjOEEQrNfrmkIXi4VwhkH9gMViIdLLJmSxWAjJnHPKpWnabrcTEhhjjuOI7DWB7XZLI1onpq8RhVsgAALtEYCRaI/tPSUHQSBvTBM9g3N1Yoxpmrbdbimjub8ovFgs6O5utxNDTKZpHi1LJGaMHU18rsJIDwIg0CwBGIlmeXZLGu2DcxxHDAqdq1+hZyBe6zTPrOu6sB+z2cw0zVPkU7fDcRyR95RcSAMCIHAXAnDLIQ+oqBYmv63X1ErX9UN/SpPJ5OPjgzbZycJPnHwW8+RyXoRBAAS6SQBGopvt0hWthsOhrIppmrTC9dBCMMaOOvmQRSEMAiDQCwJYAtuLZrqbkoUjUX3fJ1VM0/R9/4JdEXerCQoGARC4iACMxEXYHiaTPPvNGLMsi5a6TvcXY0yehX4YKqgoCDwQgcFut3ug6qKqZxLIsuxwpiFJEjqfTrYQSZL87W9/w3zDmYCRHAS6TgBGoustBP1AAARA4I4E/hfL3zhDv328xAAAAABJRU5ErkJggg==" + } + }, "cell_type": "markdown", "id": "294a660d-c458-46e8-b4bf-e663503d6dfe", "metadata": {}, "source": [ "# Britter-McQuaid Dense Gas Plume Model\n", "\n", - "The Britter-McQuaid model is based on the *Workbook on the Dispersion of Dense Gases*, by R.E. Ritter and J. McQuaid, and is a series of correlations relating maximum centerline concentrations to downwind distances that are based upon actual releases.\n", - "\n" + "The Britter-McQuaid model is based on the *Workbook on the Dispersion of Dense Gases*, by R.E. Ritter and J. McQuaid, which uses a series of correlations relating maximum center-line concentrations to downwind distances based upon actual releases. The model of the plume is a series of rectangular slices with constant, average, concentration throughout -- giving a \"top-hat\" model. Points outside the defined plume are assumed to have zero concentration.\n", + "\n", + "![image.png](attachment:3da40a83-6628-40f3-a58b-03478e3d3cfe.png)\n", + "\n", + "## Example\n", + "\n", + "This scenario is adapted from CCPS *Guidelines for Consequence Analysis of Chemical Releases* CCPS, 1999, pg 122.\n", + "\n", + "This example is based on the Burro LNG dispersion results, in which LNG was released at ground-level with the following parameters:\n", + "- release temperature: -162°C\n", + "- release rate: 0.23 m³/s (liquid)\n", + "- release duration: 174 s\n", + "- windspeed at 10m: 10.9 m/s\n", + "- LNG liquid density (at release conditions): 425.6 kg/m³\n", + "- LNG gas density (at release conditions): 1.76 kg/m³\n", + "\n", + "and we assume the atmosphere was otherwise at ambient conditions of 298K and 1atm." ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "id": "b4bdee03-c42f-4425-a66a-1148a57a1edd", "metadata": {}, "outputs": [], + "source": [ + "using GasDispersion" + ] + }, + { + "cell_type": "markdown", + "id": "7eb4678a-27cb-442a-885e-5c58260d5736", + "metadata": {}, + "source": [ + "The first step is to define the *LNG* substance, using the given parameters and otherwise assuming the substance is Methane." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "9bea6d39-ae28-449d-8bcc-4d4d3c887ef5", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Substance: LNG \n" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "s = Substance(name = :LNG,\n", + " gas_density = 1.76,\n", + " liquid_density = 425.6,\n", + " reference_temp=(273.15-162),\n", + " reference_pressure=101325.0,\n", + " boiling_temp = 111.6, # Methane, NIST Webbook\n", + " latent_heat = 8.17/0.0160425, # Methane, NIST Webbook\n", + " gas_heat_capacity = 1.6151, # Methane, NIST Webbook\n", + " liquid_heat_capacity = 2.0564) # Methane, NIST Webbook" + ] + }, + { + "cell_type": "markdown", + "id": "d24b114b-eb3c-469d-ae87-dc60cdcbab48", + "metadata": {}, + "source": [ + "The release is defined for us for most parameters - for completeness I am assuming the vent diameter was 1m and estimating a jet velocity from that - and implicitly we are assuming all of the liquid evaporates. The release is denser than air because the vapour is at the boiling point of LNG (-162°C)." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "57674931-d585-4a88-86ba-56f09a10b5d6", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Release conditions:\n", + " ṁ: 97.888 kg/s \n", + " Δt: 174 s \n", + " d: 1.0 m \n", + " u: 70.81526849717933 m/s \n", + " h: 0.0 m \n", + " P: 101325.0 Pa \n", + " T: 111.14999999999998 K \n", + " f_l: 0.0 \n" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ṁ = 0.23*425.6 # liquid spill rate times liquid density\n", + "Q = ṁ/1.76 # gas volumetric flowrate: mass flowrate divided by gas density\n", + "A = (π/4)*(1)^2 # release area, assuming a diameter of 1m.\n", + "u = Q/A\n", + "\n", + "r = Release( mass_rate = ṁ, # liquid spill rate times liquid density\n", + " duration = 174,\n", + " diameter = 1.0,\n", + " velocity = u,\n", + " height = 0.0,\n", + " pressure = 101325.0,\n", + " temperature = (273.15-162),\n", + " fraction_liquid = 0.0)" + ] + }, + { + "cell_type": "markdown", + "id": "942c716a-b605-4134-b230-ff23a79e6c91", + "metadata": {}, + "source": [ + "We assume dry air (0% relative humidity) with a windspeed of 10.9 m/s at 10 m (the default height for windspeed measurements), for completeness we further assume the atmospheric stability class was \"F\"" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "340db347-8953-43a2-b8f8-f988c2b25a96", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Atmospheric conditions:\n", + " P: 101325 Pa \n", + " T: 298 K \n", + " Rs: 287.0500676 J/kg/K \n", + " u: 10.9 m/s \n", + " h: 10 m \n", + " stability: ClassF \n" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "a = DryAir(windspeed=10.9, temperature=298, stability=ClassF)" + ] + }, + { + "cell_type": "markdown", + "id": "470a15c8-eca6-43ca-b216-9901966eef0b", + "metadata": {}, + "source": [ + "We can now put together the scenario and generate the solution." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "21bef454-8019-45d7-b177-e0a6dc9b1d31", + "metadata": {}, + "outputs": [], + "source": [ + "scn = Scenario(s,r,a)\n", + "pl = plume(scn, BritterMcQuaidPlume);" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "e0765efd-04b6-459a-b1b6-0577c66d62a0", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.08925269675562625" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# centerline concentration at 367m downwind of release\n", + "c₁ = pl(367,0,0)" + ] + }, + { + "cell_type": "markdown", + "id": "7afe7140-02f6-4a94-a366-0858f40e72f7", + "metadata": {}, + "source": [ + "From the reference we expect the concentration to be ~0.05(v/v)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "50c1a607-5133-48dd-b1c0-b3830af439dd", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.050711759520242185" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "c₁/1.76" + ] + }, + { + "cell_type": "markdown", + "id": "e892944d-e5dd-4ac4-a2c4-2d368001ff76", + "metadata": {}, + "source": [ + "The difference here is possibly due to the example calculations not carrying all decimal places, as well as the imprecision that comes with using correlations (especially by pencil and paper)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f8204273-0a71-4d42-b6c1-52947aec4ffd", + "metadata": {}, + "outputs": [], "source": [] } ], diff --git a/src/GasDispersion.jl b/src/GasDispersion.jl index 06cf62e..d021e78 100644 --- a/src/GasDispersion.jl +++ b/src/GasDispersion.jl @@ -4,7 +4,7 @@ module GasDispersion # imports using Markdown -using Interpolations: Extrapolation, Line, LinearInterpolation +using Interpolations using SpecialFunctions: erf using RecipesBase diff --git a/src/models/britter_mcquaid_plume.jl b/src/models/britter_mcquaid_plume.jl index 68071e3..1d5621b 100644 --- a/src/models/britter_mcquaid_plume.jl +++ b/src/models/britter_mcquaid_plume.jl @@ -1,44 +1,34 @@ +include("../utils/britter_mcquaid_correls.jl") + struct BritterMcQuaidPlume <: PlumeModel end struct BritterMcQuaidPlumeSolution <: Plume scenario::Scenario model::Symbol - jet_density::Number - temperature_correction::Number - critical_distance::Number - interpolation::Extrapolation + c₀::Number # initial concentration, kg/m³ + T′::Number # temperature correction + D::Number # critical length + lb::Number # plume dimension parameter, m + itp::Interpolations.GriddedInterpolation + xnf::Number # near field distance + xff::Number # far field distance + A::Number # far field constant end -Britter_McQuaid_correlations = Dict( - 0.010 => ( αs=[-1.0, -0.7, -0.29, -0.2, 1.0], - βs=[2.25, 2.25, 2.45, 2.45, 1.83]), - 0.005 => ( αs=[-1.0, -0.67, -0.28, -0.15, 1.0], - βs=[2.4, 2.4, 2.63, 2.63, 2.07]), - 0.020 => ( αs=[-1.0, -0.69, -0.31, -0.16, 1.0], - βs=[2.08, 2.08, 2.25, 2.25, 1.62]), - 0.002 => ( αs=[-1.0, -0.69, -0.25, -0.13, 1.0], - βs=[2.6, 2.6, 2.77, 2.77, 2.21]), - 0.100 => ( αs=[-1.0, -0.55, -0.14, 1.0], - βs=[1.75, 1.75, 1.85, 1.28]), - 0.050 => ( αs=[-1.0, -0.68, -0.29, -0.18, 1.0], - βs=[1.92, 1.92, 2.06, 2.06, 1.4]), -) - - """ - plume(scenario::Scenario, BritterMcQuaidPlume) + plume(::Scenario, BritterMcQuaidPlume) Returns the solution to a Britter-McQuaid dispersion model for the given scenario. -Currently only implements the max concentration at a downwind distance x, the -other coordinates are ignored. +# References ++ Britter, R.E. and J. McQuaid, *Workbook on the Dispersion of Dense Gases* HSE Contract Research Report No. 17/1988, 1988 """ function plume(scenario::Scenario, ::Type{BritterMcQuaidPlume}) - Q = _mass_rate(scenario) - h = _release_height(scenario) + Q = _release_flowrate(scenario) + ṁ = _mass_rate(scenario) ρⱼ = _release_density(scenario) Tᵣ = _release_temperature(scenario) @@ -46,66 +36,116 @@ function plume(scenario::Scenario, ::Type{BritterMcQuaidPlume}) ρₐ = _atmosphere_density(scenario) Tₐ = _atmosphere_temperature(scenario) - # Setting up the Britter-McQuaid curves - britter_interps = [ ] - concs = sort(collect(keys(Britter_McQuaid_correlations)), rev=true) - - for conc in concs - αs, βs = Britter_McQuaid_correlations[conc] - f = LinearInterpolation(αs, βs, extrapolation_bc=Line()) - push!(britter_interps, (c=conc, it=f)) - end + # initial concentration + c₀ = ṁ/Q # relative density g = 9.80616 # gravity, m/s^2 gₒ = g * ((ρⱼ - ρₐ)/ ρₐ) # critical distance - Vr = Q/ρⱼ # volumetric release rate - D = √(Vr/u₁₀) + D = √(Q/u₁₀) # temperature correction - T′ = Tₐ/Tᵣ + T′ = Tᵣ/Tₐ # correlation parameter - α = 0.2*log10( gₒ^2 * Vr * u₁₀^-5 ) + α = 0.2*log10( gₒ^2 * Q * u₁₀^-5 ) # setting up correlation for constant α # calculates the points for the linear interpolation - concs = [ elem.c for elem in britter_interps ] - βs = [ elem.it(α) for elem in britter_interps ] + concs, βs = _bm_pl_c(α) + + # near-field location + xnf = 30 + xmin = 10^(minimum(βs)) + if xnf < xmin + # linear interpolation between the near-field correlation + # and the main correlation + βnf = log10(xnf) + cnf = 306*xnf^-2 / (1+ 306*xnf^-2) + βs = [ βnf; βs] + concs = [ cnf; concs ] + else + # for α > 0.605 the near-field overlaps with the correlation + # stop the near-field correlation early + # (this can lead to discontinuities) + xnf = xmin + end + + # linear interpolation + itp = interpolate((βs,), concs, Gridded(Linear())) - # linear interpolation between the short distance correlation - # and the main correlation - βs = [ log10(30); βs] - concs = [ 306*30^-2 / (1+ 306*30^-2); concs ] + # far field correlation + # starts at last interpolation point and decays like x′^-2 + xff = 10^(maximum(βs)) + A = minimum(concs)*xff^2 - # linear interpolation, extrapolates past the end with a straight line - interpolation = LinearInterpolation(βs, concs, extrapolation_bc=Line()) + # plume dimension parameters + lb = Q*gₒ/(u₁₀^3) return BritterMcQuaidPlumeSolution( scenario, #scenario::Scenario :brittermcquaid, #model::Symbol - ρⱼ, #jet_density::Number - T′, #temperature_correction::Number - D, #D::Number - interpolation #interpolation::Extrapolation + c₀, # initial concentration, kg/m³ + T′, # temperature_correction + D, # critical length, m + lb, # length parameter, m + itp, # interpolation::Extrapolation + xnf, # near-field distance + xff, # far-field distance + A # far-field constant ) end -function (b::BritterMcQuaidPlumeSolution)(x, y=0, z=0) - ρⱼ = b.jet_density - T′ = b.temperature_correction - D = b.critical_distance - x′ = x/D - if x′ < 30 - # use short distance correlation - c′ = x′ > 0 ? 306*x′^-2 / (1+ 306*x′^-2) : 1.0 +function (b::BritterMcQuaidPlumeSolution)(x, y, z) + + # plume dimension parameters + Lu = 0.5*b.D + 2*b.lb + Lh0 = b.D + 8*b.lb + + # domain check + if (x < -Lu) || (z < 0) + return 0.0 + end + + # plume dimensions + if x < 0 + # Britter-McQuaid model does not give a profile for x<0 + # assuming crosswind length is Lh0 + Lh = Lh0 else + Lh = Lh0 + 2.5*cbrt(b.lb*x^2) + end + + # within the crosswind extent of the plume + if (y < -Lh) || (y > Lh) + return 0.0 + end + + x′ = x/b.D + if x′ ≤ b.xnf + # use near-field correlation + c′ = x′ > 0 ? 306*x′^-2 / (1+ 306*x′^-2) : 1.0 + elseif b.xnf < x′ < b.xff # use linear interpolation β = log10(x′) - c′ = b.interpolation(β) + c′ = b.itp(β) + else + # use far-field correlation + # where A is a function of α + c′ = b.A/(x′^2) + end + + # non-isothermal correction + c′ = c′ / (c′ + (1 - c′)*b.T′) + c = b.c₀*c′ + + # vertical extent, from continuity + Lv = (b.D^2)/(2*Lh*c′) + if z ≤ Lv + return c + else + return 0.0 end - c = ( ρⱼ*c′*T′) / (1 - c′ + c′*T′) - return c end diff --git a/src/models/gaussian_plume.jl b/src/models/gaussian_plume.jl index 292042c..08d0072 100644 --- a/src/models/gaussian_plume.jl +++ b/src/models/gaussian_plume.jl @@ -18,9 +18,6 @@ end Returns the solution to a Gaussian plume dispersion model for the given scenario. -The Gaussian plume model is per *Guidelines for Consequence Analysis of Chemical -Release*, CCPS, New York (1999) - ```math c\left(x,y,z\right) = { \dot{m} \over { 2 \pi \sigma_{y} \sigma_{z} u } } \exp \left[ -\frac{1}{2} \left( y \over \sigma_{y} \right)^2 \right] \\ @@ -30,9 +27,13 @@ c\left(x,y,z\right) = { \dot{m} \over { 2 \pi \sigma_{y} \sigma_{z} u } } where the σs are dispersion parameters correlated with the distance x +# References ++ CCPS, *Guidelines for Consequence Analysis of Chemical Releases*, American Institute of Chemical Engineers, New York (1999) + # Arguments - `downwash::Bool=false`: when true, includes stack-downwash effects - `plumerise::Bool=false`: when true, includes plume-rise effects using Briggs' model + """ function plume(scenario::Scenario, ::Type{GaussianPlume}; downwash::Bool=false, plumerise::Bool=false) # parameters of the jet diff --git a/src/models/gaussian_puff.jl b/src/models/gaussian_puff.jl index de9028c..8be9ca0 100644 --- a/src/models/gaussian_puff.jl +++ b/src/models/gaussian_puff.jl @@ -12,13 +12,10 @@ struct GaussianPuffSolution{S<:StabilityClass} <: Puff end @doc doc""" - puff(scenario::Scenario, GaussianPuff) + puff(::Scenario, GaussianPuff) Returns the solution to a Gaussian puff dispersion model for the given scenario. -The Gaussian puff model is per *Guidelines for Consequence Analysis of Chemical -Release*, CCPS, New York (1999) - ```math c\left(x,y,z,t\right) = \dot{m} \Delta t { { \exp \left( -\frac{1}{2} \left( {x - u t } \over \sigma_x \right)^2 \right) } \over { \sqrt{2\pi} \sigma_x } } @@ -27,6 +24,9 @@ c\left(x,y,z,t\right) = \dot{m} \Delta t + \exp \left( -\frac{1}{2} \left( {z + h} \over \sigma_z \right)^2 \right) } \over { \sqrt{2\pi} \sigma_z } } ``` +# References ++ CCPS, *Guidelines for Consequence Analysis of Chemical Releases*, American Institute of Chemical Engineers, New York (1999) + """ function puff(scenario::Scenario, ::Type{GaussianPuff}) diff --git a/src/models/intpuff.jl b/src/models/intpuff.jl index 5489be5..345f3aa 100644 --- a/src/models/intpuff.jl +++ b/src/models/intpuff.jl @@ -13,7 +13,7 @@ struct IntPuffSolution{T<:Number,S<:StabilityClass} <: Puff end @doc doc""" - puff(scenario::Scenario, IntPuff; kwargs...) + puff(::Scenario, IntPuff; kwargs...) Returns the solution to an integrated Gaussian dispersion model, where the release is modeled as a sequence of Gaussian puffs, for the given scenario. diff --git a/src/models/simple_jet.jl b/src/models/simple_jet.jl index 230ce55..910e235 100644 --- a/src/models/simple_jet.jl +++ b/src/models/simple_jet.jl @@ -14,14 +14,11 @@ struct SimpleJetSolution <: Plume end @doc doc""" - plume(scenario::Scenario, SimpleJet; kwargs...) + plume(::Scenario, SimpleJet; kwargs...) Returns the solution to a simple turbulent jet dispersion model for the given scenario. -The turbulent jet model is per Long, V.D., "Estimation of the Extent of Hazard -Areas Around a Vent", *Chem. Process Hazard*, II, 6, 1963 - ```math c\left(x,y,z\right) = k_2 c_0 \left( d \over z \right) \sqrt{ \rho_j \over \rho_a } \exp \left( - \left( k_3 { r \over z } \right)^2 \right) @@ -31,10 +28,14 @@ where *r* is the radial distance from the jet centerline. Assumes a circular jet with diameter equal to the jet diameter. Ground-reflection is included by method of images. +# References ++ Long, V.D., "Estimation of the Extent of Hazard Areas Around a Vent" *Chem. Process Hazard*, II:6 (1963) + # Arguments - `release_angle::Number=0`: the angle at which the jet is released, in radians - `k2::Number=6` parameter of the model, default value is recommended by Long - `k3::Number=5` parameter of the model, default value is recommended by Long + """ function plume(scenario::Scenario, ::Type{SimpleJet}; release_angle::Number=0.0, k2::Number=6.0, k3::Number=5.0) # Density correction diff --git a/src/source_models/jet_source.jl b/src/source_models/jet_source.jl index d4ddb1b..f00afb6 100644 --- a/src/source_models/jet_source.jl +++ b/src/source_models/jet_source.jl @@ -7,8 +7,8 @@ Returns returns a scenario for a simple jet from a circular hole. The jet can either be a liquid or a gas (in which case it is assumed to be an ideal gas and the jet is isentropic). -Liquid and gas discharge models are per *Guidelines for Consequence Analysis of -Chemical Release*, CCPS, New York (1999) +# References ++ CCPS, *Guidelines for Consequence Analysis of Chemical Releases*, American Institute of Chemical Engineers, New York (1999) # Arguments - `phase=:liquid`: the phase, either :liquid or :gas diff --git a/src/utils/britter_mcquaid_correls.jl b/src/utils/britter_mcquaid_correls.jl new file mode 100644 index 0000000..5fd4aa7 --- /dev/null +++ b/src/utils/britter_mcquaid_correls.jl @@ -0,0 +1,116 @@ + +function _bm_pl_c_1(α) + # corresponds to (Cm/C0) = 0.10 + if α ≤ -0.55 + return 1.75 + elseif -0.55 < α ≤ -0.14 + return 0.24*α + 1.88 + elseif -0.14 < α ≤ 1.0 + return -0.5*α + 1.78 + else + @warn "α= $α is out of range" + return -0.5*α + 1.78 + end +end + +function _bm_pl_c_05(α) + # corresponds to (Cm/C0) = 0.05 + if α ≤ -0.68 + return 1.92 + elseif -0.68 < α ≤ -0.29 + return 0.36*α + 2.16 + elseif -0.29 < α ≤ -0.18 + return 2.06 + elseif -0.18 < α ≤ 1.0 + return -0.56*α + 1.96 + else + @warn "α= $α is out of range" + return -0.56*α + 1.96 + end +end + +function _bm_pl_c_02(α) + # corresponds to (Cm/C0) = 0.02 + if α ≤ -0.69 + return 2.08 + elseif -0.69 < α ≤ -0.31 + return 0.45*α + 2.39 + elseif -0.31 < α ≤ -0.16 + return 2.25 + elseif -0.16 < α ≤ 1.0 + return -0.54*α + 2.16 + else + @warn "α= $α is out of range" + return -0.54*α + 2.16 + end +end + +function _bm_pl_c_01(α) + # corresponds to (Cm/C0) = 0.01 + if α ≤ -0.70 + return 2.25 + elseif -0.70 < α ≤ -0.29 + return 0.49*α + 2.59 + elseif -0.29 < α ≤ -0.20 + return 2.45 + elseif -0.20 < α ≤ 1.0 + return -0.52*α + 2.35 + else + @warn "α= $α is out of range" + return -0.52*α + 2.35 + end +end + +function _bm_pl_c_005(α) + # corresponds to (Cm/C0) = 0.005 + if α ≤ -0.67 + return 2.40 + elseif -0.67 < α ≤ -0.28 + return 0.59*α + 2.80 + elseif -0.28 < α ≤ -0.15 + return 2.63 + elseif -0.15 < α ≤ 1.0 + return -0.48*α + 2.56 + else + @warn "α= $α is out of range" + return -0.48*α + 2.56 + end +end + +function _bm_pl_c_002(α) + # corresponds to (Cm/C0) = 0.002 + if α ≤ -0.69 + return 2.60 + elseif -0.69 < α ≤ -0.25 + return 0.39*α + 2.87 + elseif -0.25 < α ≤ -0.13 + return 2.77 + elseif -0.13 < α ≤ 1.0 + return -0.50*α + 2.71 + else + @warn "α= $α is out of range" + return -0.50*α + 2.71 + end +end + +""" + _bm_pl_c(α) + +Britter-McQuaid plume correlations as digtized in TSCREEN + +# References ++ EPA, *User's Guide to TSCREEN* U.S. Environmental Protection Agency EPA-454/B-94-023 (1994) + +""" +function _bm_pl_c(α) + concs = [0.10, 0.05, 0.02, 0.01, 0.005, 0.002] + βs = [ + _bm_pl_c_1(α), # (C_m/C_0) = 0.1 + _bm_pl_c_05(α), # (C_m/C_0) = 0.05 + _bm_pl_c_02(α), # (C_m/C_0) = 0.02 + _bm_pl_c_01(α), # (C_m/C_0) = 0.01 + _bm_pl_c_005(α), # (C_m/C_0) = 0.005 + _bm_pl_c_002(α) # (C_m/C_0) = 0.002 + ] + return concs, βs +end diff --git a/src/utils/monin_obukhov.jl b/src/utils/monin_obukhov.jl index 8021220..eab581a 100644 --- a/src/utils/monin_obukhov.jl +++ b/src/utils/monin_obukhov.jl @@ -2,8 +2,10 @@ _monin_obukhov(roughness, StabilityClass) returns the Monin-Obukhov length for a given Pasquill-Gifford stability class and surface roughness (in meters) -Curve fit from - Pasquill, F., *Atmospheric Diffusion, 2nd Ed.*, Halstead Press, New York, 1974. + +# References ++ Pasquill, F., *Atmospheric Diffusion, 2nd Ed.*, Halstead Press, New York (1974) + """ _monin_obukhov(zR, ::Type{ClassA}) = -11.4*zR^0.10 _monin_obukhov(zR, ::Type{ClassB}) = -26.0*zR^0.17 diff --git a/src/utils/pasquill_gifford.jl b/src/utils/pasquill_gifford.jl index 633f929..4f917ba 100644 --- a/src/utils/pasquill_gifford.jl +++ b/src/utils/pasquill_gifford.jl @@ -1,8 +1,12 @@ -# Plume crosswind dispersion correlations -# Reference: -# Spicer, T. O. and J. A. Havens, "Development of Vapor Dispersion Models -# for Non-Neutrally Buoyant Gas Mixtures--Analysis of TFI/NH3 Test Data" -# USAF Engineering and Services Laboratory, Final Report, 1988 +""" + crosswind_dispersion(x, Plume, StabilityClass; avg_time=600.0) + +Plume crosswind dispersion correlations + +# References ++ Spicer, T. O. and J. A. Havens, "Development of Vapor Dispersion Models for Non-Neutrally Buoyant Gas Mixtures--Analysis of TFI/NH3 Test Data" USAF Engineering and Services Laboratory, Final Report (1988) + +""" function crosswind_dispersion(x, ::Type{Plume}, ::Type{ClassA}; avg_time=600.0) δ, β, tₐ = 0.423, 0.9, 18.4 δ = δ*(max(avg_time, tₐ)/600)^0.2 @@ -39,10 +43,16 @@ function crosswind_dispersion(x, ::Type{Plume}, ::Type{ClassF}; avg_time=600.0) return δ*x^β end -# Plume vertical dispersion correlations -# Reference: -# Seinfeld, J. H. *Atmospheric Chemistry and Physics of Air Pollution*, John -# Wiley and Sons, New York, 1986 + +""" + vertical_dispersion(x, Plume, StabilityClass) + +Plume vertical dispersion correlations + +References: ++ Seinfeld, J. H. *Atmospheric Chemistry and Physics of Air Pollution*, John Wiley and Sons, New York (1986) + +""" function vertical_dispersion(x, ::Type{Plume}, ::Type{ClassA}) δ = 107.7 β = -1.7172 @@ -85,10 +95,15 @@ function vertical_dispersion(x, ::Type{Plume}, ::Type{ClassF}) return δ*(x^β)*exp(γ*log(x)^2) end -# Puff crosswind dispersion correlations -# Reference: -# CCPS, *Guidelines for Consequence Analysis of Chemical Releases*, American -# Institute of Chemical Engineers, New York, 1999 +""" + crosswind_dispersion(x, Puff, StabilityClass) + +Puff crosswind dispersion correlations + +References: ++ CCPS, *Guidelines for Consequence Analysis of Chemical Releases*, American Institute of Chemical Engineers, New York (1999) + +""" crosswind_dispersion(x, ::Type{Puff}, ::Type{ClassA}) = 0.18*x^0.92 crosswind_dispersion(x, ::Type{Puff}, ::Type{ClassB}) = 0.14*x^0.92 crosswind_dispersion(x, ::Type{Puff}, ::Type{ClassC}) = 0.10*x^0.92 @@ -96,15 +111,30 @@ crosswind_dispersion(x, ::Type{Puff}, ::Type{ClassD}) = 0.06*x^0.92 crosswind_dispersion(x, ::Type{Puff}, ::Type{ClassE}) = 0.04*x^0.92 crosswind_dispersion(x, ::Type{Puff}, ::Type{ClassF}) = 0.02*x^0.89 -# Puff downwind dispersion correlations + +""" + downwind_dispersion(x, Puff, StabilityClass) + +Puff downwind dispersion correlations + +References: ++ CCPS, *Guidelines for Consequence Analysis of Chemical Releases*, American Institute of Chemical Engineers, New York (1999) + +""" function downwind_dispersion(x, ::Type{Puff}, stab::Union{Type{ClassA},Type{ClassB},Type{ClassC},Type{ClassD},Type{ClassE},Type{ClassF}}) return crosswind_dispersion(x, Puff, stab) end -# Puff vertical dispersion functions -# Reference: -# CCPS, *Guidelines for Consequence Analysis of Chemical Releases*, American -# Institute of Chemical Engineers, New York, 1999 + +""" + vertical_dispersion(x, Puff, StabilityClass) + +Puff vertical dispersion correlations + +References: ++ CCPS, *Guidelines for Consequence Analysis of Chemical Releases*, American Institute of Chemical Engineers, New York (1999) + +""" vertical_dispersion(x, ::Type{Puff}, ::Type{ClassA}) = 0.60*x^0.75 vertical_dispersion(x, ::Type{Puff}, ::Type{ClassB}) = 0.53*x^0.73 vertical_dispersion(x, ::Type{Puff}, ::Type{ClassC}) = 0.34*x^0.71 diff --git a/src/utils/plume_rise.jl b/src/utils/plume_rise.jl index b9043a8..f9fd391 100644 --- a/src/utils/plume_rise.jl +++ b/src/utils/plume_rise.jl @@ -29,9 +29,14 @@ Base.isapprox(a::MomentumPlume, b::MomentumPlume) = all([ if typeof(getproperty(a,k))<:Number ]) """ - plume_rise(Dⱼ, uⱼ, Tᵣ, u, Tₐ, ::Type{Union{ClassA, ClassB, ClassC, ClassD}}) + plume_rise(Dⱼ, uⱼ, Tᵣ, u, Tₐ, StabilityClass) Implements the Briggs plume rise equations for buoyancy and momentum driven -plume rise as described in the ISC3 model guide EPA-454/B-95-003b +plume rise. + +# References ++ Briggs, G.A. *Plume Rise* U.S. Atomic Energy Commission, Oak Ridge (1969) ++ EPA, *User's Guide for the Industrial Source Complex (ISC3) Dispersion Models, vol 2*, U.S. Environmental Protection Agency EPA-454/B-95-003b (1995) + """ function plume_rise(Dⱼ,uⱼ,Tᵣ,u,Tₐ, stab::Union{Type{ClassA},Type{ClassB},Type{ClassC},Type{ClassD}}) # physics parameters diff --git a/src/utils/wind_profile.jl b/src/utils/wind_profile.jl index 0350f14..c5b7dfd 100644 --- a/src/utils/wind_profile.jl +++ b/src/utils/wind_profile.jl @@ -29,16 +29,16 @@ end returns the windspeed function u(z) for a given Pasquill-Gifford stability class, `z` is assumed to be in meters and `u` is in m/s -# Arguments -`u` friction velocity -`zR` surface roughness -`λ` Monin-Obukhov length -`stability_class` Pasquill stability class (A, B, C, D, E, F) -`k` von Karman's constant, 0.35 +# References ++ Businger et al. 1971 +# Arguments +- `u` friction velocity +- `zR` surface roughness +- `λ` Monin-Obukhov length +- `stability_class` Pasquill stability class (A, B, C, D, E, F) +- `k` von Karman's constant, 0.35 -# References - Businger et al. 1971 """ function _windspeed(z, u, zR, λ, ::Union{Type{ClassA},Type{ClassB},Type{ClassC}}; k=0.35) a = (1-15*(z/λ))^0.25 diff --git a/test/models/britter_mcquaid_plume_tests.jl b/test/models/britter_mcquaid_plume_tests.jl index dda0a67..db8732b 100644 --- a/test/models/britter_mcquaid_plume_tests.jl +++ b/test/models/britter_mcquaid_plume_tests.jl @@ -1,3 +1,22 @@ +@testset "Britter-McQuaid plume correlations" begin + include("../../src/utils/britter_mcquaid_correls.jl") + + @test _bm_pl_c(-0.75)[2] ≈ [1.75, 1.92, 2.08, 2.25, 2.4, 2.6] + @test _bm_pl_c(-0.40)[2] ≈ [1.7839999999999998, 2.016, 2.21, 2.3939999999999997, 2.564, 2.714] + @test _bm_pl_c(-0.20)[2] ≈ [1.8319999999999999, 2.06, 2.25, 2.45, 2.63, 2.77] + @test _bm_pl_c(0.0)[2] ≈ [1.78, 1.96, 2.16, 2.35, 2.56, 2.71] + + # test domain warning + @test_logs (:warn, "α= 1.1 is out of range") _bm_pl_c_1(1.1) + @test_logs (:warn, "α= 1.1 is out of range") _bm_pl_c_05(1.1) + @test_logs (:warn, "α= 1.1 is out of range") _bm_pl_c_02(1.1) + @test_logs (:warn, "α= 1.1 is out of range") _bm_pl_c_01(1.1) + @test_logs (:warn, "α= 1.1 is out of range") _bm_pl_c_005(1.1) + @test_logs (:warn, "α= 1.1 is out of range") _bm_pl_c_002(1.1) + + +end + @testset "Britter-McQuaid plume tests" begin # Britter-McQuaid example, *Guidelines for Consequence Analysis of Chemical # Releases* CCPS, 1999, pg 122 @@ -12,8 +31,8 @@ liquid_heat_capacity = 2.0564) # Methane, NIST Webbook r = Release( mass_rate = (0.23*425.6), duration = 174, - diameter = 0, - velocity = 0, + diameter = 1.0, + velocity = 70.815, height = 10.0, pressure = 101325.0, temperature = (273.15-162), @@ -21,15 +40,44 @@ a = DryAir(windspeed=10.9, temperature=298, stability=ClassF) scn = Scenario(s,r,a) # known answers - # set 1 is covers the short-distance correlation - # set 2 is the answer from the above example, which uses the interpolations - x₁, c₁ = 2.26*20, 1.1827612798486664 - x₂, c₂ = 367.0, 0.08922951841378725 + # initial concentration + c₀ = 1.760006673092362 + # in the near-field + x₁, c₁ = 2.26*20, 1.1827642935170182 + # example, in the interpolation region + x₂, c₂ = 367.0, 0.08925280058298538 + # far field + x₃, c₃ = 1200.0, 0.008660197082795383 - # test type inheritance + # test overall solution pl = plume(scn, BritterMcQuaidPlume) @test isa(pl, GasDispersion.BritterMcQuaidPlumeSolution) @test isa(pl, Plume) @test pl(x₁,0,0) ≈ c₁ @test pl(x₂,0,0) ≈ c₂ + @test pl(x₃,0,0) ≈ c₃ + + # test plume extent + Lu = 0.5*pl.D + 2*pl.lb + @test pl(-Lu - 2*eps(Float64), 0, 0) == 0.0 + @test pl(-Lu + 2*eps(Float64), 0, 0) ≈ c₀ + + Lh0 = pl.D + 8*pl.lb + @test pl(0, Lh0 + 2*eps(Float64), 0) == 0.0 + @test pl(0, Lh0 - 2*eps(Float64), 0) ≈ c₀ + @test pl(0 - 2*eps(Float64), Lh0 + 2*eps(Float64), 0) == 0.0 + @test pl(0 - 2*eps(Float64), Lh0 - 2*eps(Float64), 0) ≈ c₀ + + Lv = pl.D^2/(2*Lh0) + @test pl(0, 0, -1) == 0.0 + @test pl(0, 0, Lv + 2*eps(Float64)) == 0.0 + @test pl(0, 0, Lv - 2*eps(Float64)) ≈ c₀ + + # test large α near-field + a = DryAir(windspeed=0.8322, temperature=298, stability=ClassF) + scn = Scenario(s,r,a) + pl_nf = plume(scn, BritterMcQuaidPlume) + c_nf = 0.9346808405962113 + @test pl_nf(pl_nf.xnf*pl_nf.D,0,0) ≈ c_nf + end