Commit 2d62dcaa authored by Markus Holzer's avatar Markus Holzer
Browse files

Minor changes

parent ad081175
Pipeline #38773 passed with stages
in 45 minutes and 19 seconds
......@@ -18,7 +18,7 @@ class LatticeBoltzmannBoundaryHandling(BoundaryHandling):
"""
def __init__(self, lb_method, data_handling, pdf_field_name, streaming_pattern='pull',
name="boundary_handling", flag_interface=None, target=Target.CPU, openmp=True):
name="boundary_handling", flag_interface=None, target=Target.CPU, openmp=False):
self._lb_method = lb_method
self._streaming_pattern = streaming_pattern
self._inplace = is_inplace(streaming_pattern)
......
%% Cell type:code id: tags:
``` python
import pytest
pytest.importorskip('scipy.ndimage')
```
%%%% Output: execute_result
<module 'scipy.ndimage' from '/opt/local/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/scipy/ndimage/__init__.py'>
%% Cell type:code id: tags:
``` python
from lbmpy.session import *
from lbmpy.phasefield.analytical import *
from lbmpy.phasefield.eos import *
from lbmpy.phasefield.high_density_ratio_model import *
from pystencils.fd.derivative import replace_generic_laplacian
from pystencils.fd.spatial import discretize_spatial, fd_stencils_standard, fd_stencils_isotropic
from lbmpy.phasefield.fd_stencils import fd_stencils_isotropic_high_density_code
from lbmpy.phasefield.cahn_hilliard_lbm import cahn_hilliard_lb_method
from lbmpy.forcemodels import *
from lbmpy.macroscopic_value_kernels import pdf_initialization_assignments
from scipy.ndimage.filters import gaussian_filter
```
%% Cell type:code id: tags:
``` python
A, B, C, D, F = sp.symbols("A, B, C, D, F")
tau = sp.Symbol("tau")
omega = sp.Symbol("omega")
eq = sp.Symbol("eq")
fq = sp.Symbol("fq")
test = fq - omega * (fq - (eq + ((1/omega) - sp.Rational(1, 2)) * A + B))
test.expand()
```
%%%% Output: execute_result
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAToAAAArCAYAAADov0TWAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAJYElEQVR4Ae2di3HVOBSGczMUwLIdkA54VEDoAHYrWOggzFaQyXaQbAUMdACpgEcH0AFsOmD/T9HxyM/7kmXf+GhG15YsS+f8sn4fPay7+vXr11FJt1qtHqq8dyr3cclyvSxHwBFYLgKrCYjui+B+pHJXy4XdNXcEHIGSCByXLEzW3AuV96hkmV6WI+AIOAJFiU5wP5d/D+yxC8upO0fAEXAERkWgGNGJ2C6kCf571Oj+qJp55o6AI+AIRASKEF203n5oXA6S+xHLfui14Ag4Ao5ACQSKEJ0UeSOS+ycqZBbdgy4FRYqXXfEWB2nKf5P3sT4DxY+OgCMwiMDoRCdCOpUEHxIpjOhOkrhwqrRMVqRpm0ksjDVIvu4cAUfAEViLwOhEJwleypoLExBRmpt4hKya7nkjbfP6Uez+Yh1aPq00HuEIOAKOQIrAvTSQ+1wW2pnypKuZdkety1qbjIiW37tUBsWRhsXFzNam7psCH9MIP3cEHAFHoA+B0YgukhQWWJOkWFrC5xhNiw5r7k1D0CcKd1luJ0p71UjrQUfAEXAEOhEYjehU2oXI6HVnqbdLTMyysyRN4iOeCYe3loCjSLIrXZrEzx0BR8ARqCFwXAtlCoiMmCjgU68+h5VW67oS1n3VTGo8Z91dZREqLnRlFZd2hRV05wg4Ao5APwLZLTqRERYX42q/9Rd79JNrpFU6m4X9qqhrxdElJQ/Ikg//v8Q8SfeHPEtV7B4F3TkCjoAjMIxANotOZIRFxtIQJgo4h6BqVpvCL2IaSAz3TmGsNty5/Gf5MIGh42MRGuRH95exOkjuXHGzH5uTTkzADFm0UsWdI+AIlEKg+O4lpRSbspxIcrPYoUWyvBIWL+Xt5ZIu9QEmrGcsbCxlXizuDgAB1SvDNwwB4TAA7Pw25o7+7qp39q7rHcV3Y7VUEbPaoSVawFeS6z8p8VlhSK/mdI1lPVjgZkXXrntgXgionugFMZkH2V3Lf5JvvsAUdbfcPnpn67oeEqQC7FQeS2cMl2WHlpwyKi+sNoYR+r46scmdvlny7Djl1C+7cPPPkGeX1QiQHcM9s1tTOlL97qz3Ui06Gj0+q1Pl8qbFG2HsU0ZOGa3b2tcgIEIc46ulXE79Ssk8eTkQiIQAu4+xu9pXp1PLmrV+99V7kRbdGE+AKgKymOsOLWGJjhpG3xgcxHwjf+e7P2PUfeE8qcvvkeQKFz1pcXvp7USXr+423qElX5Eb54QV0Hrzi5yZHWd8jrfvMzWeatmO4q0721kIxC4/2S4yKhvZL+TPor8krims4pCT2X3SkZ77OE/XbM5aV3SSvK/kqasw5MK5fK/curY4vYd0XmrXtdke9goLYIgkHf8ywmjt0LJXQTvcLNlo0IEAdM7SHXNPdYIVeimCq01QKN22u8j0WYpWVtZjxJtGzoYRoeyo59+Kqz4jTNLxeWGoE8VBFi8UDv9ZovCsdTXgJC/LqphU4vNJZllt2zNLUh2XqPc6nZ3oqsdjrxManI3LkRHdQBxEMrWDhHGvJaMRcIiIjfxfHfkmOV2fCDGk+oT06Q956b7iu8ioTDDlpYI+KcHWXjZKB7lDal0LzNP7ZqurZK+5qDtxLevcEi5R7010vhcTMUXdMvsNvI5j9Sa1a8qn7P8mxoL1sPf+m5hk4q1vDd1E5chsFd2Brsb8VXnWLBzS9jnlERY4x7IsWchfgbWYFpCRsY0biMmEs6Pi3qt8LDu6fSw9+aojeEEQlVNcIA1dD2N91YXbyYveRkc63Zu7DsgPF7qgt6dHv+v4QfKlsjApdF9xKYGTHAs3jEXm1DVitHc7QsABh+y8lFKibiYvrXfu+m3qQ3id3mt1BrTFeQFHd+VsX92VBwTQmY/iIf5vu5ahe3PJiBx8ktdZz7rGmA9pgh46shlDLa3iAvl1xLfSNtP0hZXnTvrpvkF9rDylY90g5FfpojDWIPefEi/Xkp9r8i28utKmeZc4Rwb5wWdK12eht+TYqX67cFyn9yY631Mid7sjsO0OLbuXtMOdsjLCm1C3puOHzZzMSjMroau7TT5v0xuVd1e6NMmY55+GMo/WFS+hpt6QGMRnll+XDnPTNVUV2aye0vhwvkS9N9X5uIWWR2yEgACm0Qx9z3qj6zS2KV1o2BLAGnZNFulAw+HNS3fd0tAlJD64eI4lYYRIdxS96N5aNzKkLfRDF5yuastJLnRJXbO7zpBEShRz1zXVhfMn8oMkH29Yot6DOrtFF5+MbQ5qUFgCu+zQsk0xOdL+qUw6x+ekAyQIWdHwn8mbIzznXWQg3drYaiRe4vFYbDeKg7gri01huujonM5Wzl1XiXvrJD+68IJJiTpevT0sUe9NdXaiqz0qw4HYoCAHGgyWDRYd68+w3oJTHFYFDTGk0ZH1Tqxir5Y9hIQj/qg8ZKRRBMsshtMSaTQ/5f+SXM1FwufxPiZZaFR8/8oMKzpBJFgVLG9oDvIrenxHuZIFSwxZ7EsOJh1q5KdrWG/MKNuSGnvjp93ZWevaQDPUpfQ0y7txuQouUe+1Ok+ye4kePhqaNXwaDo2u2O4ZKh8yYi+89O2uqPm4Q5BxH7RK66fysOhYM9g7S7+PPmPfK/kZJngg+WnUG7up9M5Vv7vo3aXz8caIZUooISA5HjjWQeHZXBPLgd0zzArKVFpvNlhglRXWm2raC4cg4z4IldaPMUaes4Nxag/s34h1jsMgeBvOtvuZSu+d6zeD3i2di1t0seJaizgVz7Q44ytDOxNvV8We2hGICMTn60rPl/UkZo+NZKZrTpcbK+46GgVbyb1Evbt0Lm7RqZaw2vhG8n6jxhh7YOwFi8+dI5AFAT1PfN/KWCrPGxZSmLDIkvn4mSArY4oct+2yLk7vobqeyqJj993ad6ASEhOdsTPf/FEguHMEHIF8CBQnuj7RRXSY6UwQHORgcZ9eHu8IOALTIzBF17WltUiOqXO6rAczftJSwiMcAUdgtgjMwqKL1txWH9PPFlEXzBFwBGaHwOREJ5JjfRBd1uoTo9mh5AI5Ao7AQSMwaddVJMciThZBOskd9GPkwjsC80ZgMqITyTHDeiKSq6bNFcf2z4zVuXMEHAFHIBsCkxCdyIzJh6ciuebkA+TH52DuHAFHwBHIhkDxMbposbEIsuvjZDZErK2vy6apZ+QIOAKLRWAKogvr5XoQZ+aVb1/dOQKOgCOQDYH/AeOPjnn3R/QnAAAAAElFTkSuQmCC)
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAToAAAArCAYAAADov0TWAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAJYklEQVR4Ae2d25HVOBCGzdQEwEIGkAGXCBgygN0IdskAat94oyADIIIpyACIgEsGkAGzkwH7fx63S9bx5fhYln1wd5VGlixL3b+O2t2SrLn269evIic9f/78ltp7p/huzna9LUfAEdguAicLiP5Obd5ZoF1v0hFwBDaKQFZFJyvukXB2JbfRH5uL7QgshUBWRSchHyq8R9jKheXSyRFwBByBWRHIpuik2F5KEsKPSqLrs0rmlTsCjoAjUCGQRdFV1ttPxSi5n1Xbt7wXHAFHwBHIgUAWRSdBnknJvaoEMovuRpuAKve6Ld/yUJoK3xV8rs9A8dgRcAR6EZhd0UkhnYmDDwEXpuhuB3nlpcqyWBGWjYtYGmuQep0cAUfAERhEYHZFJw4eS4GVCxAVN5dVjLKK6WFUNr5f6D6KEuvQ6tkp4xmOgCPgCIQInIaJ1NdSSk9VJ65m6I6ay9pYjFAZLDT22NWkPMqwuZjV2pC+K/ExzPBrR8ARcAS6EJhN0VVKqmhRUuTxOUZs0WHNPYsYvad0m+V2W2XfRGU96Qg4Ao5AKwKzKTq19lLK6Elrq1dbTMyysyKx4iOfBYdzK0CsOtvKhUX82hFwBByBBgInjVSihJQRbujXnuqw0hquK2k9V6+kVtfsu6vdVuWVrqzyQldYSSdHwBFwBLoRuJb6o34pIyyur4r/6GpW91hZRRnigparsIpRav8o4JJSB/cfKKAwmY+j3J8KbFVxt1VAODkCjsB+CCRTdFI+Zm2hoKBvCg+Uj/VWkq7ZPoI7G5b5qHyUV/g8z7Ja+0MB5YcShF4obfvxrnJW+Fc8oqj9hJYV9o2ztE0Ekim6bcLXLrUUHVboHcXX2kvkyxUPvCgeK9jLJdzqAyMo5QsFXja8YJyOAAH1FdM3ZkRgANj1EXB/OIuHyn16eJP+ZBsC6ohVndAifnDz3yj+T/EXxSi9BimPbT1MN9xVcGXXQGd9CfURHg6LeSi7TwqfFeIXmLJ+L5oi98nvBcV+0giwMwUsnTkoyQktKXlUXVhtTA10fXViiztdq+TJcUopX3Lm1l8hv91zBZTdF4XV7SmdqX8PlnurFh2DnpCU1Lm8aQmmMKa0kZJHc1u7BgSKEGIjdi5KKV8unhdvBwUiJsCOue1L4sWZamcgaf9OlXuTFl17v0zLVUegLNZ6Qku5RUc8drmlKGYGzW/v/kzr5VU8TV+ySEd/bYkmyb1Vi26OHwiT+WbJ2cEFuBZrIKyAnTe/+OWt+1aBmBVy47vQ9WsFk0e3m6R7KHZcYVbHuxRo86GEKbUJz/8q2LFfHBJBHzQUQMUnVjbzWDcVXijgAmERlXwrXrWs4rUQj/DMYKcvL5RmXpW4tY+UT/9sSu4+mV3R6dcwlQQwP75w/ssUxs4JLVPbGvu8eGMTNkqh0DXfHhvd1wWDgUHeWKBQeuwpMlkVnfgDb+YVayWrPORE8dWfEQbl+Lyw7BPFKIhHissVccWrllW8liQ+bVGJzyd7t1mprOGzGbmHZD6pcPRoGgJjTmiZ1tL4p/nRQ0/0Y3gVBJQb1g2f6mEthLTaU2TEq1mS8B0q2MbLRvdQ7ig1ytmLx2QMn1utrMasxZXsJHes86DM5uTep69Pq0IsUZdvfQNsIK7fpFZO9eT9v4lVw2q3c6+a7vHWt4FurBKXLqXut5n935TfsHDCB+Nrld37hJb4WdJ6fm4ecXcu1U482Gn7vQKWHVYdW0+QHbxQEDUprxw8iqkrpMFTZPRMavmoD+KTQbNQcUk/KB0qANw2ysRf0WD5lXORupdMVtUFRpPHkeroI3gv1FaoqOPyueVO3b+xPKSH5B6UGUV3qYom/49V1dOpcNo4z5EnntoUWaF83BWOj5r0lYWeL18OimMFUCiv7YSWHbFVblYe1SCDuW+RwVZaKccAwsKp3T+loXsK/E5iGjxFZgb5SnlU71Df8blgqPgKPYM1aBYhsiSTVXUnGUcw1UO8lHZeWFH53HLP/ftFvCG5B2U+jUDy5DgExp7QMq72iaU1+Mo3oaoJ5w/jWk1Jm5WAIoiJes7DzEpphFk5rz/3NSbeeAERYrlRkoXumwI8Bllh2Yh+sH6yvDreotz7ynxSo+QXoxAQwAyasSe0jGojQeFyYKseG9iNKiUDAwfrFpfVyqz9FBksGlzVHZIMyBJSbP0wJREqirXLGsrCNZZ1r5KvHtii3L0yu0VX/TLGRBpQWAJ8tN95QovuX1AnZRXiTuBWDvpLjVy2ta88lCBzcQx8TokxIv1J95nbQk7KMbXBJ2KkkQVXga0cS8jFfEzDXRIfWG/kEwqlkRnFDb8lKc2CC7KELu/aZb1iXn/FP7IgJzy3kspsTu59ZXZF1/qTac8UqPzQUA4MmEJpLLqhE1pQiOUJLTyTg9QePMJr6bpW6bBpBs2Fwt+6F8/fvaiee6qYQcX3r2xQRbmgSLAq2N4QT/Ire36iXQUsMXix+UXSDeWne1hvb5WPHJAp5dCdXbWsV2zXf60vzfKub0QXW5R7UOZFTi/Rj4+BZhPeDBwGXbbTM9Q+Ls7kxQjVMRsdA49ThM8tn9rDomN1eXWLZvvgKL5Z3byhmEG9Ny0lt9pNMsYOkbtN5pO9EUtUUEyg5Mqd6LpmbxduEZYDrlFpKSVqqq+aS90krJmOgccp+OWWj0UXfmdHQxoPbGzGOocwCM7Lq3F/lpL74P5NIPeOzEu4rrgcDTdDgmHN8calU/vmvcZ1cUdptTVk/nc8mS/7GHicgsYC8vESXcTdnoATY4Vpg+vUoTieZtin6kXknti/U+XekTm7RaeegYnv1nlBT6F8mGvB4nNyBJIgoN8TW4CYS0VZYCExiI6F4JU5ReKxLuvm5O7r6+xzdGIGq43TdxvfgVb5+PV++KNAcHIEHIF0CGRXdF2sS9GxgsYCwVFOFnfJ5fmOgCOwPAJLuK47Uku5sXSOy2orsTtlPMMRcAQcgUMRWIWiE/O4s3xgHm7mPFQmf84RcAQcgQYCi7uuUm7sD8JltW8uGwx6whFwBByBqQgsatFJubGlhE2QruSm9qQ/7wg4Ap0ILKbopNxYYeWYn3rZXNdYdszVOTkCjoAjkAyBRRSdlBmLD/cVx4sPKD8+B3NyBBwBRyAZAtnn6CqLjU2QbV8nnOl+Y39dMkm9IkfAEdgsAqcLSI6Swz1lfi6mo/oWMWbe046AI7BOBP4HJQMqc3w9Y1oAAAAASUVORK5CYII=)
$\displaystyle - \frac{A \omega}{2} + A + B \omega + eq \omega - fq \omega + fq$
A⋅ω
- ─── + A + B⋅ω + eq⋅ω - fq⋅ω + fq
2
%% Cell type:markdown id: tags:
# Implementation of high density difference model
According to *"Ternary free-energy entropic lattice Boltzmann model with high density ratio" by Wöhrwag, Semprebon, Mazloomi, Karlin and Kusumaatmaja*
Up front we define all necessary parameters in one place:
%% Cell type:code id: tags:
``` python
a = 0.037
b = 0.2
reduced_temperature = 0.61
gas_constant = 1
κ = (0.01, 1, 1)
λ = (0.6, 1, 1)
χ = 5
φ_relaxation_rate = 1
ρ_relaxation_rate = 1
external_force = (0, 0)
clipping = False
domain_size = (100, 100)
stencil = LBStencil(Stencil.D2Q9, ordering='uk')
fd_discretization = fd_stencils_isotropic_high_density_code
target = ps.Target.CPU
threads = 4
threads = False
```
%% Cell type:markdown id: tags:
## Part 1: Free Energy
The free energy of this model contains a term that is derived from an equation of state. The equation of state and its parametrization determines the density of the liquid and gas phase.
Here we use the Carnahan-Starling EOS, with the parametrization from the paper
%% Cell type:code id: tags:
``` python
ρ, φ, c_l1, c_l2 = sp.symbols("rho, phi, c_l1, c_l2")
critical_temperature = carnahan_starling_critical_temperature(a, b, gas_constant)
temperature = reduced_temperature * critical_temperature
eos = carnahan_starling_eos(ρ, gas_constant, temperature, a, b)
eos
```
%%%% Output: execute_result
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAkQAAAA4CAYAAAAcn+fJAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAYaElEQVR4Ae2dgbXctBKGSc4tIIQKCB2EpAKgAyAVJHRATirIgQ4IFQDpIKSCQDpIXgUkt4O8//PV+MqybMu2fNe7O3OOd21ZMxr9I1vjkSzf+vTp02dOjoAjcLwI3Lp164m0f69r+e/jrcX5aC573VNtv9X2QRv7D7U9l/3e6t/pgAjINndU/LOgAraBHss2l1e7/nusCMi2v8mOP43pfzF20s85Ao7AvhHQRf69NPxaF/qLfWvq2kUI/KJ9HNjGZrLhzzp+re3zKI/vHgaBX+JOk05Uavyr7avDqOOlVkTgF9nzlez73ZDMWzo5dM7THQFHYMcI6OLmafa1ruGvd6ymq5YgILvdJ0l2ayJCOsZB+l7H3ukmWN30oWxBh/idbNFEW3VMlOidNh46PIJ30wapXJ7sSTT9jmz5a0707VyipzkCN4GAGmfTMcRlKe1OuAk1yexre0J6nG9v+9Jvsi4b6PyXZPIE67QhArVtS8eadK5E+XCKnGYiUNs2Kp4hlX9mqnES2TfAcne46LojKvuT6prtTzxCtDuTzVdIxrWb6X/i5imTsO/7UklL+MPF85PK6YzJKp0nqqeh7Af6Z57E06QDaE4r70ft0DDtycsaKU9jl2RSHjoLOv4hulTeZqhBeeG38X/yc/yXzvfm1gQ9TXfyoTe45fJO4ltSF8mvRkH/f63u1QRHglTGZL2j7L3dUv7a+VBEMnFQaTdtWyLdSOertlOTO+dfOjBU9kjbH7Jj9ol1jry1eUvtMFROKX9JPuUpso/yTd5DhvQtSQ+6Hjx6V4LZWH1K+LfGcky/WudUh9HrnnKUh+vuoa65H3rlKpHQrW9HioEMyvg2F2xjQ+3TuRPivWdpY//Kt4g/lIGz0bYdytTGGG2cRqdKGPrbOJ19EXpyQ+OAffISzkz5XymdSEi6wRPXnUlzLW8og06xzRPSwKiTlzza0CPNW4SP+Cbrkuq25ljlgQUOXKe+tY4lu6jeQ+WV8tfMJ1mNA6x/sEF/7NlpT+grqt5Oh3CYSkc/bbTvn6fybnle5e/G3tKl2D7Ku9l1J9l0rsjvtaEtbZHKVvk3ZZvNsEzrVPMY+2izaPngdW9lhvz5e4Nl8v9tOpYtcZVhGQ99l5ahNByLjmOS5uFYtIhffHjYXDypQ0Sj7DliSsPp+ZjqkPKn54OO2U6fcrS1To32qcuTVIbSuFiIpLSOg47RHwBa50f75COtzav9YnyUt4NFXN4W+0HXnpNZo6w59c6VV8pfO1+si2SbjXudmc5Vb6dx2XP3pQ8dL23v/lxeyy9e3lrrtX87P/YPn7bN7yOl5ShfsX3IO1a3peckt3HK9N9rP3NlSsbubUOdRJtgWYLXGoxi+ZIzeN0n+XCcetfLbSU6HS8ChPxsuCmuxRsdfKvQIJ38GM3mDyHJSwllS4kL/12mXIahOnODUsaRY+qSIxwlG/LiPEOFg28PJALArFMHycrVZzY+STmbHAYbIHuruQ5r613KXztfKd5btNOisrk2tH2MbAiftT2GmJfSHTGyLaFSOwzJLuUvzXcw+1BB2QZniGF+JldfchzShuo/lX4Mtpmqw9bn12C0RDfunbTHDrlD1IHj6A64cTBHJyWbP8T5MVrC/0g3iaFXvHF8eJ34cqDQ2TdsyXqZytLNiQjY8yQdx+l7nXulLS6HvAyhtCSZf2v7nH9LFA9DZlCcdwk+V1K2/UWvMZzXlr623qX8tfOV1rt6Oy0tWG2Oa4Nr1q5RWMGB9D85OACV2mFItVL+0nwHs4/uAzhD3AN+0/59Nu0zJzJ3n1Xy5lSK2ZAia/mH5B57OiMcvQeQi2Ov1bnqrws17vSHYLg7dGIJv3gIR8YOQ0e8bvY9jztk4KbCkBWRmQ5JJuF66vKFNm5GowvUKT95GFqwids6bGS/1DmcJxwbnsA5T9SIocOeU6X0lpSXmwaOE5PEbW2Y2fhKzqy6oEAoG9zoEKkbT6ZxZ6mkHj1UCvmrk/SZXe9YiVL+2vliHab2he9NtFPa1JBdSX8mDP4LumLP7OTvcH6zv1I7DClQyl+aj3Lm2keyZ113yj9mG4ZSuAb4b0k6xdHoNn3LnTmY5fRYwl8Zy5xae0mjLyJay5Boey91h2gv5pmvhzk7rTEzIsY6t1n8ajh01pdqPFOddUcN8eEMNR1958TVAfr9aQ0ylMGQW7sOSIYHx4WtR5Lzg3hx2LhBkgdds3mVbm8icXOkQ+ICIYxqNAsfMc2ui3RFNy7I5marY3ThRjy1QB9lbfXEOrfeUqVDpfy183WUmHsg7O+Lp0o7nbKr7E1bY9sDldphSNdS/tJ82XJG7DPruiuwzdS1l9Vvo8RVmEmnufxVsdwIk1pi7f7JNd9ei7drSXc5u0SAqMsaivnb6MlMgUyQfKlOoPdasdKaMXqTp2McGMLl2SiUbmY0XiYSt0Ndxsu/zhMdutRGZIg85MfBsuEwHV6T5LAezK/aeGL/Q9u/Q3mvuTp7LT6SMbcu6MSk7vjJE4eMpxYcozHiRkc9D0VtvRcqUMpfO9+YulXaaWg/S+06pt8hz5XaYUjHUv6xfFn7zLnu3DZZ87SYnxmWdv80p7EBxyNE2TZyFInm4eaUNSNbSD6Xp5hfNxIiLlknJSfY0sQHD3NdcDhKCaeICeG8rZZGo3Ae0rRGbtCRYQdzML5TGk4HN9Lftc+8IbsIGp74R+cYcuP8X/rnKbEYn1hOsj9Wl9+VtzPsp2McOOjO1d/gL+cH9ZP+nH+tbUpOXMAPwuDtmFydq9muBvVPyinNF9dl1r7wqtlO19h1Uu+ga85hbmyj89b+Y1k4/kPX4Fp8S/lL88V6N/sL7DN03bltrtEtuZbJXRXLDdrvdY3m7WXboztE80DcTW46dzUu9Ml1epaWdR5gKuVXGXTSDOsMykJeSuLDiborvuybXzrP2iuc/zrlDcdWh/g0Ds6QHgw/fRlnlmycHKJFTKCjE2nmEintPvl0/i3/ERGhIR9RKHg5ldPD0hpd5tZF+SkDGekkWtKhVK+r1OvfS+3aDe06NexJd84P4drLHyfAW1rvmM/2S/lr57Py5/6rrtXaaQW7Tqov3HIOj0VHeYjoRWLHhJbaYUhGKX9pvrScMfvoXPE9xG2TItve12bfw9ZgqXZQtf32alWeYPfPjmN0Uc5/vjnVAOi8ngUE7Cn+MRf5gVGxYaFUDTN2dmgpylzC/0D5HwoDIi0x4VTwOirpRIHaaIfScFy+Ulr7VKq0BjelmUOD3E5j1DHU6K58HadA/NgAGZ10GMK5zuQ40iHK03kcIcOE5GbCpNJ502zMhiX4IG9WXZQfXHhqT8t+pHSwNIx0eBAqrfeQcqX8tfMN6ZNNl/1rt9O92zWLgxJL7bCWf1Y5BfaZc925bbrWs/shNoHOAcurml790p9AnXvw7as0/51AgDVvePuHjQuLjrzzFsIE/1ancUZoyCkRHch1uGm+SX7Vl6EmhlM6mwTRkOxc7AzhKLEsepsWCqXziR2gF8pD9CYloiR2kcbnrJ6xjOa85KALkQ1zVmM+9mn8sUzyM6+J/5isDMs7iU9gnlsX6kg0qiXpDm5stK8pwmGiTltRab2Hyi/lr51vSJ9eesC7djtda9eenjeUUGqHIXVK+Uvz8ZBTch+Zc925bbrWS/uIc8AyRuBuOOj2J+oUeIr2bQQDAccynu2qwNqn4121smwtzKUHw0HpisusDN2ueqt9Ok/0bVdhtvKVNslveeP/IK+zGrbSwAV5zMlIt85KuCFvu9I0skW81t98myguK5zDoSLT0MrV3PBwUokUte1ZxwyJdD6LENI6q5QqzeSn6ZP4iJd6F9VF+cwW7aqwIY1yOmXH9Yj3la/BKU6rva8ySuptdVncrkrKoW6l+QwH5WcIFcbcyunV26nKMSwW29V0X/Kv8mm/nXY+R454d2Nv6VJkn5Bv8rpTPrfNdB9RdA/bCkvJXdV+ra1LzuB1b3n4F9Ev9L6ecMEZp0kEGPfsPNFPctxcBjz9X/RE9VD/TKLm/xsZ/a3+G9I+0ROiCrk6TPIHMc2f5ODocPFATH7mqe+NymDuAmP6nKOxpdTqwwnlZygLvZEH4bHjrX+JviQkBD/pb5L05lA8f0vWYx0wgRo5RnRQFvFp0nT8QnnQ3comHb1zr/tP4iN5c+qC4wY9Vfk4NhCRMiJwHYyaM/kf8t0Rf3aYMM8yO7Wk3jXa1WQ5QfOifMKE9ggZzkySp+3jvNuColu0UytvjV0bxQ/0M4mv8LspexfZZ8Z157aZ7iNK72G7xLLwuo8vLdp7p1/g5Ml97V7A0LFhtHZ9Gypak1QGXihRmdyQT82iXNaJIaC2gxP2QG2HC3IxSQ6RNJyo3kW9WKgzLkagll2XKqDyecKePal6aXnHxOe2qWetrbC86far8hhNILpoD0kNSBe1oFIBOAkQUQocBYY2iieHlvArD86OzVInDMox5cSdAmPPdDosva6/Hl0q/+LFtyQT+dx8VnVoPa084VwQwFlnkvda4g013uCL2/5amc6/HIFadl2qAdHTXGR1qbxT4nPb1LPmVljeWPtVH47vQD+evuX72UUNnFQA3hafXGhu9KFAFrljCGLSKSrhDzKZ1GwOkb1qyrereFK2ToYhIzqJXLkYM53sq6QyUjmNA6bcB1lmv0xLz7VXBEIbpg1lh/1m6o3T/1rb4vY8szzPPoBAZbsOlDKerPufO8YZiNw2GVAWJm2J5Q233x8FQe6lmvUOkUBivghzGcwhYX4IY80cc9POrkOj9IZm8FPOE+WPv01lN4FnOheX3ytTfHRE6Nbm47iUAj8OWSM7kpdzvCbFih/njBB3J2Q3yegZjhmBB0F5a7eL66J281ZtiHlTDN0uatOLC3fGFIFqdk0F+/FqBNw2qyFsBZwKljxE9nwEanm7reryHV4Tzk0G5SmYiauEp8aolJ8yOmE1dQQcpzT09M3QWhtdSpnGjlUHnCkbhqv1BWRwmcJmTC0/d3wI/COVO5/4WFmFx+LnYcDpsAjUtutha3Napbtt6tnz6LFUX05ghblD2UBGDYeISMeHDOZWIOfHqIhfFWDNGxbSa5+ueToOgnFWGtL5l7Zv/8rH/KbndrzgnyFB9OTfNl6RzjlkC8Q7yzkgQHuJ2+/aOof29zy077XinH8hArXtulANZ8sg4LbJgLIw6dix1H2SwAbTawZXc1/lEKmAkggHr1NnaQ2/eHFQcHRGPzoaQGBNnlwUq9ULedqYiM2r4LyqC3gNiRdHjDfyOpud939H4FAIqE3yAMAHbKcePA6lopfrCDgCjsAeEGCojBGpQboYPFN2wpydsUjJmNM0m183fmaHc/Nn8jRODmG8McJpsjfgsvkkk/PMg2qG1ELnQiRo8dto2YI80RHYAAG1W5+HtgGuLtIRcAROBwHr38dqtNYhGpNt576wnYX/HX5VCieoifbIcWHIjLfZ4rfM2mKUTpSHFaYHvcIgI11TCCeLhe/gbYfoWsG+4wg4Ao6AI+AIOAInhcCFOn0iOLy+OxbJSSuNA4JT8iE9ER1b9Oe/KC3dXcUvHfgiOdEphrhyH+ok4mNzmdKy7fh37aSvLttw2RxMTF7nX3oxvyk3nNHgo/O5id58h2zQiaMA8TXrj3cK8wNHwBFwBBwBR8ARWITAhTpeHIpFiwzCq46ZgnOOg6UNOiRz+FUOQ2W8Nt9Eh9gPRDQHh4ON+RQxEUEaLF8y4UHPdIEmc2DSsmLZRfvSN+fw2BpKi1eWldzsqpNFSnkmR8ARcAQcAUfAEeggsGpSdZDEkJJFVGLhFiGaGnIq5W/e7pITY45WXFZvP+RDr8veyesEojC5r8I/Uvp7OR2DztS1CN9zBBwBR8ARcAQcgWNH4KJCBfiYYm7SMlGnnLORFlnKj2PDq/epgzO0WJSljw3LEQnqOGwhEkU0alHULK2cHzsCjoAjcIwI6F7IAyX3SO6h7PMiC18kWB05lxwnR2B3CNxeq5EuDt5w+aCLh+GphkJ0huWxWTiuIdK0fdJGpKelUn4xMM+HryC3JFmUScSIV+9TR4l0KE1vEtFHO1zkFsliGIs0HDTk+UUvIJwcAUfgbBHgQfcr3Qv5zAFrt7DoLfNNnRyBk0SgRoQIYIimsH4PTxBMoub/m9ip0D7zjRiCyr0mX8L/QvzNWkGSYYRDw8q/nShPOIlDgzPERZwjnnygp5L789Vu81FamzAekvzPEXAEHIGzROB5UusvdEy0yMkROEkEWGjwJCs2VSk5Qbz99UD1P8jQmMonurV4UvVU/fy8I+AIOAI1EdA9653k8QkkX/eqJrAuazcIrB4y201N5ivSmz80X8QqDqJXbE6OgCPgCOwWASLo2pjqwDeg3BnaraVcsbUInGWESBc3c4U+assu6LgWVOd3BByBOgjoWiWK6m971oGzlTIX13DPZH7lK9lj8FtQbQG+4wgcIQLnGiGyN9Byc4+O0IyusiNwegioE2ZuX7P+2OnVbl6NhEUPB5wUbcyjXEL3Ar5FvHKCiGbzYgtzRXu6FAnxTI7AzhGoNal659XsqcfEbiZj+5BVDxpPcAQOj4A6XebYfTEWjQgdM1GLrw95LUsPW3aEF0q+0sY8m+KoViH/a+Ujsm1vv7IPLZoDKf3+Bj9tT7TfGwYLZf1P8uOXY+x+yQOl6YEOTo7ASSBwlkNmJ2E5r4QjcKIIhM74tTrqXmcfzvG5Hd52omMmWpH7bI+Styfpw9wa1uZpVskP+pHGA9ekU1TKr3xMaL6rDUcIuZRHueak6HA+hfJxenpyQpmts6njJyoB5+/LXP75pTuHI7AvBM41QrQvK7g2joAjECNAp8tboD0KHTErzLNu2EGH1IKDcEc6Nc4QOqGf0jlG/+9IG6KZ/JPfNxwqZyIdPcE794khcH4mPYl8QSyn0jpITYr/OAInhIA7RCdkTK+KI3AiCPwoxyLXQe+tejgMuaEj1j7jzSycpV7kJarEWv5I1LJd6cf6bh+1PU111TF1y9VvWWHO5QjsHIFznVS9c7O4eo7AeSKgjpm5Q5NDTTtBh6U7GLpLyfTn/Bit5R+TPecc+vJlASdH4KwR8AjRWZvfK+8I7A4Bhpl2//Yn0Z8C5Jjzk6Ul/OJhDg/lsmI0b5cxh2gwgqP8OFxEoYhSkZ8okDlrOmwJvMG9N7m6zeE7jsAZIOAO0RkY2avoCBwRAkyUzs4f2lkdzNkZGxIbc5rm8iPrTxvWkrODg/NO/9lPFymdeUEM2TVDj8E5YrL359pSYsI2zpOTI3DWCPiQ2Vmb3yvvCOwOATr+3DDU7hQtUIhIzhpq+eXYdJYJCZEeIjs951HOD8OO35szFBT4R/+sW5RzfMAbB8vJEThrBDxCdNbm98o7ArtDgMjJWNRlscJyBnC2+Fo7/6U09LHnMafNoj/2dlaurLX8yGT4iw9ep6t5sywBiyjGZA5Pru7IyaXH/L7vCJw8Au4QnbyJvYKOgCMAAoqYXOqvt7bREnSQJUcE1pwjYWm5+TpNcXP4Vc4rMd0Vz5DuVh5LERAB4vjPpqDrH4sMvb1Oavc2c0LbEnzHETgCBHzI7AiM5Co6AmeEAJGTtoPfeb0ZsrLIS6yqRYimJoeX8jOvKodJU44cpdjJaV7lx+GKFdL+I23vlZ5z0pA9FrFKRPmhI3CaCLhDdJp29Vo5AseKAB12zsnYY334bAjOSkpEclhIMXVK0nyl/C8ki0+CpETUJ3W6SGO+UEuKGt3XARvOUo5wrHKOUi6vpzkCJ4uAO0Qna1qvmCNwlAgQ7XhYqLlNOm4iJYU81bLJSeE19Q9yOJjE3JD2ibawps/jq5RmRW0mM3/SxlteLZXyi+E38XYmT+uYVbqh1slRGmXjTLZ4hDQcr59UXhxJUlJLjQPXHvmOI3CmCFycab292o6AI7BPBP6QWnTgg6RO3s7bvJi/lEaE41VwMgZ5NziBM8EX4HHimETNf/xBVPucB/p1Ijc6hkr430s+ZZhThMPDEFf6TTHD46nymsNEZGloYrhONQTfMawMHtT1P0dgGwT8467b4OpSHQFHYCEC6sxZF2eqE18o/XTZgsP0QE4hTlYRiYeIEo5kbkiuSIZncgROBQEfMjsVS3o9HIHTQWDoY6OnU8NtakKkJ51TNFUSr+eDt5MjcPYIuEN09k3AAXAE9oVAGPa6F6IX+1Jup9oIK5s/9KZUxYAvaxj5JztKQfN8J42AO0QnbV6vnCNwtAgwWdjmzBxtJW5QcXvbbU6ECHx97tANGsmL2jcCPodo3/Zx7RyBs0UgRDD4BMWvZwtCYcVDhIj5Q0UOkfIz6fql8vvr9oUYe7bTR+D/2Yoz0PlgbuQAAAAASUVORK5CYII=)
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAkQAAAA4CAYAAAAcn+fJAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAYP0lEQVR4Ae2d4bXcNBOGN/fcAkJSAaSDcFNBoAMgFSR0ACf/+JdDOgAqSKCDkAoS0gF8FRBuB/d7H1+PI8uyLdvaXe965hyt7ZFGGr2SrPFI9t65ubnZOTkCjsDpIvDTTz89k/b/6Pjn6dZiO5qrnb5Qbb9S+KjA+SOFF+J/0NHpiAioDe6q+Oe1CrQN9FT869tT/z1VBNSGvyh8P6T/5VCkxzkCjsC6EdAA/0Yafqnjr+vW1LULEPhZ5xiwVZvp+IOu3yp8FqTx0+Mg8LPao5k0df6L1PhL4cFx1PFSCyJA275R+Lovz4u+COc7Ao7AuhHQwK6eZnVsbuDr1ti1qxF4oeOrAI37Osdb5HR8BJ5pPOG9M8J4/UK8h8bw42kioDb8R5r/riMPIElygygJizMPgUDqJiPeXQVzVe84V+AmxeS/WpJ+nRsmOis0ddmD8r8rT55gnfaIQOm2VX4fCIHKePmYeJ0mIlC6bVQ8DxfvJ6pxFsn3gOXqcFEd8cp+r2NyPvEls9U12XSF1Lh2M/1X0rh2cQ1iDWfRHHnJYADQsVreCV1jAPxYF3ylI0++P4ofTgB19O6t+HRMi7NO+qUl0JFymPRZ/9WhQ9fiV0sNdV62/k9C8uOJoLO3Rjz0NN1JxzW4pdLm4JtTFxVRhmr9r3Tsdf8uLUl559S7t5hc+dLpUEh50m8wGFlOvIYXknil+2mYfda5dOBJ9YkCffvoS57SYTXtPaF9io67RDtwj2B50+5RWW1bOtGB2qYolqUxyMlPOA2O+zoP5pPfFL6N87yIGX59WgioA7C+/U5HjI6XOscYYZ2UG/4oLZBnsrkXFlCXWW1c0znGEoYNN5K/dB66oU3so06YrOjEGCV/KMQTGBtOMVKYMOKA0fdUwQiDBhwscDNDD57AG9I1ZVkajqRjkIBbnDYX35y6NDoUOKGd9zaJCofceierkitfMp3ywiOHAWw3vOQYUDz80v00icMQU3owXh8rfK3zXjf+UB6l4lT+atp7YvvsbdxJD+5L3A/CB7RSkGfnc6i2kUJ7wzK7sjMSCp+scR9kzX3zG+QCXnXqHqIYkRO6VoPydhGdAUOiIp3jMeGaSWHQezBXXnJ9N2+eMM3rYvpgcKAnBlS8aZSlg46VXgkGP0rTqYd41WSnY1V3HSmDm3pMGExsWG0w0jlpWYbDADK+eYaeW9o6z1x8s+qivEsRdejgUiLzifXuFJkrv4d0GNdVf1Le9FEmtBTtpZ+mChrjSU/GK8YtDw1/KnwYk0nFS44HDpaXJxvJkqEv5fbzTvG58rnpVMCU9tnLuJOu3F/QI35A69R/jKG8Vt82dR32guUYPsQvxCh33FeqqCzGHOPsO4XWeLmoUvjPqSLAzT91A30n/ldq9I4FHFV0srzyZJKhAxJiYuD/nSgXY4MbbmXExEIj19QlRXiDQuOLpcJcAwHMWnVQXqn6TMYnpWhpnnS1if596bzr/JbWO1e+dLpcOPbRT7PKVtsxDv4L2hA563ssMc+luxIkzKHcdujLO1c+N93R2ocKqm24T/Egh+eOyRNDc869y/A6hbYxXY91XILRHJ25d9IfW3TRuvKLU0OAG8fHhNK2f4j4IZoj/0Q3h5ZVHRSA4cN6+3XAC08n37CVl3lwmnzE48ntRcO4PcFwwg2K1ycsh7R4yxpSPE/in3E0ps5tqSxMOwcfy3KfR/Qawnlp2UvrnStfOl1uvYv309yC1c8YG4xZG6OIggP811wcgXLboU+1XPncdMdsHwwf7gEsqT4k6BwPXuo+K/beKRezPkWWyvfle+r8v1WBzgPI5anXaqv6a6CGk34fDPf6IubIS4ZliNBgaGWv+I7FXSeoPBqKxzPTIvEqd72YvHrMzehFKp0JKY403Ki4STWk6z8IYmDY8AROPF6jcFlMl11SWm4alZte55Wxp+NkfCUzqS5oIhnKBjcmROrGk2k4WYrVoUfikL44qezJ9Q6VyJUvnS7UYexcZR+inw61K+U/lx7/1rrSnouXZsbqnYrPbYeULLxc+dx0dZ6T2kd5Txp3Sj/UNiy7MwZay++SCb3RqLl3moJZSpk58pIpiWVKrbXwmIvw1hKae+nlWrRzPSYjYMZO05iJHIYmt0ny6jRM1riPxybrlhpKjzFUTfStiNsL9HttHbIugyU3XNWN9yaSw3AhdEgy3ypgsDGoSYOuybTi75QW3bg5MiExQMIlqEn4SHZyXVQ+ujEgq5utjujCjTjeayVWiyhrX0+sU+vdUkwXufKl08V6TLqu+0KRfjrWroqnrxHWQLnt0KdrrnxuumQ5A+0zadxltM3Y2EvqtyfmIsyk01T5oljuCZNS2dr9kzHfjMWLUrl7PqtEAK/LEgrleVurb6lsqAw2U+O94Y2aFolXrdEbU9cYMBhCSS+U4um8X+mYNJbExzt0rYBniDSkx8Cy5TBdfiLx2UT4UoEn0lcKbGxNpv0k1Tpr8JHc1LpQDkt84ZMnBhkGEobREHGjo57HoqbeMxXIlS+dbkjdIv207j9z23VIv2PG5bZDn4658kPpku0zZdx52ySbp8F8Y1ja/dOMxgoc9xAl+8hJMM3CTSlrjWwu+VSabHkNFDwuSSMllbHxJIcMe136XOCWNDxiFGH0sJEx9kZhPMS8SlZp0ZFlBzMwMFAwOriR/qZz9g3ZIKhkwh/FYbQRz2vbPCVm4xPmE50P1YXvYLSW/XSNAQfdvT30/hLfq5/0J54368byCQvAu/ZhKF/FlexXvfpH5eSmC+sy6Vz1LtlPl7TrqN61rimDuWobxVv/D/MaentoKb658rnpQr2r8xnt0zfuvG0+oZszlkldFMs99N9PNZp2luyPbhBNA3E1qdWxWL5Cn9SkZ7yk8YBQrrzSMUnjtejNi/xiUnoMlHs6Jt/8Ev9NHf9lLFtfWx3CaAycPj1Yfvo8TKwyMHLwFrGBjkmk2qAt3kPS6fiBY0B4aEiHQYYsUSk9jFfponST6qL0lEEe8SZa+FCs1y330++1Tu2G9olbnyl/4vtw7aQPGcgqwLI6htHG62uDXa586XShklPOpUexfqq8lrbrqOoqI2Xw7MRnbPAQ0fHEDmWq9Ktq71hX6dfbPorLHndK623TBrc1lg+Fpcop2n/bVZp0ZffPlmF0OSmLjSZWI9J5ntfVt6f4p+JfHxkSWxaK1bDGTi4tBYlz5K+U/pHqiqclJIwKbsDw8QI13g6dc3N+oGPjGdJ5hZuONpmSb6sz6hqqdFe6llGga9qAPFp8BOo4jLZOe4iHbhhChgki1YZJ8XnTrCNDgppy8CHppLooPbjw1B6X/UR89DWMdHkUyq13n3K58qXT9emT5Avn0v107e2axEHM3HZYKj+pnIz2mTLuvG3arWf3Q9oE2gKWtzW9/WU+gVr3YDeIbkEZ+21980YDFRc7kyreh2MSxkhq0zDegdSEG+uaI8+AsUHTyAuD/+Dr2Bg9ROoaQwkDqjGQ4IuYfMI9SL8m0pCOJ7lOeeIxYKGOEaV8eMolpJbZkKHzh3kyCFJLaFaGpc3Bh/yn1qVTR+kOboQczw4Gkw1onRan3Hr3FZwrXzpdnz4dfo136X66tF07eh6IkdsOferkyuem22W2z5Rx523Tbr14jtgCliECZhC25pOLMIWf9yLwTAOUAWWEEcLkywR2NFL5GBgfdcTYqEjnTJR8gZMvNFcET+FGofKMBPwseUsfHSmnNSkrfzw43PQoj+94NEE8NmVf62hUxdkFR8XzWj/UMrJuWU1ZYR51VHVAhv0/sU643Hn1PvS6YKzhbm9I8WCIbKOneLn4ZNel1g+cbEDuah64UXbHAyZ+TO/EMOMtjlt8LR2y6o3eCrP71YRysvSJKn6/vm5wtniVW7yfgoXyX9qupuJBjxPa4VDtnds+WePO22Z8jlCHOxcse8d9NKjoY9fqG6355DJK5JdpBFj3ZH/JGglLHw/WIx3ZRM3xsa6biZVGV8AgSNVhVF5yDSkfvGN0Joi9Nkzk/JcaexcwMojDCImp0YcIpWdpCL3JD2Liwlr/XLxWJyVShDx8jIEOSQaPD0YgG6jJxwgjyTw+FU/XPA2hu5UNH73ZiN1KK94oPpKZUhczrPnekBmAeBptU7NORwksmJwIKaxGM8hIkFPvEv1qtJxa16x0woP+CBnOtD99H6MYwwraRz+18pa06612x/kdxZe+VmO55D4yWk5u+6CLQs49xNtmfI44aSzVD3LGfTiy6IfxvX535+bmJkx08ucChomNAdB836Z0pRiEypNXa4+9ZFa6ap7fnhFQn8EIu9KRATmbJM+SJUZUZ1DPztQFZyNQql3nKqDy8XBO3lQ9t7xTkvO2Kdda+8Ly0P1X5bFaglfMHpIqkC5LQaWMbS8LXgoMBSz3cJlisKgceaXB2LFd6nd1zjXlhJMCy1hMOlRWhw7xlDP741uSJX9uPosmtI5WztgKAhjrbPJeSryhxht8Yd9fmqfLz0egVLvO1QBP4b68hXN1Wouct025ltgXlgfrv5rDsR2Yx+O3fHeXJXBSAVhbL3SsbvR1gXzkjiWIUaNIaUbl6zxxR5tBtNM5hgmucJ6UbZJhyYhJIlUujRlv9hUrj1RGZYAp9VE+s5+npadaKwLqPwxE+lBy2W+i3hj9fGtodn+eWJ4n70GgcLv2lDLMlg5uGCcg8rZJgDKTtU8sD9x/2WNr351roXHRuppxoYqwX4S9DGaQ7HSOtcd1uEcjmfsEecphczNGkJHdBJ4bg6PSYIixQbUJYuPBYp9Joydpc0lyTGQYZOSNlwn3NLxZJFn2sKT22szKz4VOAoGrWkvrt7OVVt9hHxH9ORwPs/NzwUUIFGvXRVq4cAoBb5sUKvN454IlD5HJB8nFBpEy5u0ebs4x8RTMpM9T8RDlylMGhhahIuXdnBtPx76nb5bWGu9SkH70VHIYPrYMV+ofkMGF4LQdBN6rqpVBXajKbCJvPQwUytezmYZA6XadVrqnHkLA22YInWlxJ4+l5nKcEGynSa0gFVkyYxmqtTGpxtgKHFtzzJJXBXiqbu39Ec+ejhtPlHh/1OU3B/HwDr1oGNNPWNLDeOHYkPKdZWA1GfjJphBQf8GAX+wdMtDIT4Glaoz95BOPpfXj/hCgHZR7sXbdn6bby9nbplybnzqW0h/HBttr2HuZpMskN5OpjHM8HLxOnaQl8pLFkMLQYVksZZBVZSoOEPDqDE4YdX54q7i5IcPyWGXU6dgyxBTn5AisAgH1TdbC+YsUvLE+Ka+iVVwJR8ARWCEC2ADM8b20yCBSrmbsYET00ZDRNFleN312h2MMPVJgGQ033hBhNBF6SXkSzz6oyuOjI/njDXJDqBc1j1gLAuqvvQ8Ea9HR9XAEHAFH4JgI2Pw+pMPFUGShuPsL82nJq1L8JcVLBSy9Vwq8zWZLZ62ixMfTM/jkXMvyTaFw+QsjCwMJw8jJEXAEHAFHwBFwBM4cgUtN+nhweH13yJMTw2Bf1f0YRwTX5v35N+DFp4vkpbu9OsfXaFN/1ImRY3uZ4rLt+jedxMtpGFLQFExuJaJf6cX+ppRhVeGj+NAQM2mMvkHXnuLP64uaVnM/OgKOgCPgCDgCR0AAg+ha5c76yCCyCqidMhyM12uQTJFXWpbKdjqyTBYS3hwMDkK8oRrP0VD5yKBn/IEm+FBc1i13wq/0TRk81APdZn9ZVvJ3JqjhSR0BR8ARcAQcAUdgAIESS2Zs5DSPSliUeYjGNnrmyrOnh+UxM7TCsjrndTr0uu5EfmLghUn9K/wT8flvl15j6lMWfuYIOAKOgCPgCDgCp47AZYEK8KdqqU3LeJ1SxkZcZK48hg0foosNnKs6w9jwMv7QshyeoJac8scTRZjlNat18YMj4Ag4AieNgO6FPFByj+QeyjkvsvCZh8Wec+Xj5AisDoGLpRppcPCGy0cdWQKqSOd4cfg8Nh+Oqwiewo1C/C2fLHllwj4f/qW6IeVFmZTFq/exoQQfivkVU+mJZ5CbJ2tX8zDQyM8HfYWU/zgCjsBGEeBB94HuhezVfKlzPnrLflMnR+AsESjhIQIYvCl8HI4nCDZRc3wcGhU6Z78RS1Cp1+Rz5H+VPG+MNR9hVF4YNHz5t+XlEQ/CoMEYYhCniCcfiO8N/XB7Wv0prW0Yr1l+cAQcAUdgkwi8iGp9X9d4i5wcgbNE4M7NzTZfVqoNqysdj7I0pnLxbs3eVH2WvdEr5Qg4AqtFQPesv6UcD77+3avVtpIrtgSBUh6iJTocSxYPUfxW2iF1wXuVXM47pBJeliPgCDgCQwjIAMKDzosm/AeUG0NDYHncSSOweA/RKdZeg9r2D/Utp+29WtKBDeJ+c9k70l7AKSOgMcKyuFNhBKbgqrTsH3qswPYE215QWCPPzhE4PgKbNIgEu72Bltp7dPxWcQ0cAUdgV0++vPG5eRIWHRzE40WVuQYjy/XZxo3S4s3mxRaWzDq6bL6BHICzQGCrS2Zs7OZpx5eszqIbeyXODQGNTfbY3dcR70SS6omZt0K/POZYVtn22RFeKHmggNGQ/Q2zTPm3Sodn295+5RyatQdSeeGh5k+vnyl0PNXikf//FMKXY+x+yQOl6aFTJ0fgPBDY7Kbq82g+r4UjcH4I1JMxBkBnsq/j+Lsd3nZiYsZbkfrbHrH3T9KHz4jwbZ5qP2KtHzweuEaNolx5pWND8z0FDBXypTzKNSNFl9OpLh+jp5NPXWZjbOr6mUrA+Ps8lX566S7hCKwLga16iNbVCq6NI+AIhAgw6Yaf12ji6omYL8zvdM6Sz9GWb1Q+BgLLVs3LGTrn8yJco//XCr00UX70/w17CxqOQE/wTv3FEDg/l554viA+p9IYSBXHfxyBM0LADaIzakyviiNwJgh8p0k4NUGvrXoYDKmlI17W+EF1wFjqeF6CSiyVD7Kadyr9+L7bfwp8j62lq66pW6p+8wpzKUdg5QhsdVP1ypvF1XMEtomAJmH2Do0uNa0EHT7dwdJdTKY/8UO0VH4o7ylx6PvdFAFP6wicIwLuITrHVvU6OQKniwDLTKt/+1OG290MiNnzk6Q58pKpluiUIV+M5u0y9hD1enAUh8GFFwrPD+nxApmxpsuGwBvcO5urmxR+4ghsAAH3EG2gkb2KjsAJIcBGaTYQr53M2GktM0VKDxlNU+XJ67UMmpcKvP5O+Ks2eqJiq/1V7Avib4j4X0bSsleIzd4pAm8MJidHYNMIuEG06eb3yjsCq0OAiT+1DLU6RTMUwpOzhBp5GTWtz4ToGk8Pnp3O5nPFsez4jY7hPqz34rGnKbWMB95uEAkEp20j4Etm225/r70jsDYE8JwMeV1m6ytjAGPrrQLHXOr7s+cho828P/Z2VqqspfLkiVHEH17zkcVwKYzPEuAVCskMnlTdkU3xQ3k/dwTOHgE3iM6+ib2CjoAjAAIyGq516HzbaA465KWAaMqQMF5opLSKmSKvtG8kfE/HPt2tvJ3S4AHi+nWrQBlO9fWHiM/l3ozQRFnOcgRWi4Avma22aVwxR2CTCOA5aSb4lSPAkpV5XkJVzUM0tjk8V559VSlMqnJkBIVGTvUqv3jXoUI6f6Lwj/gpI428hzxWUVZ+6QicJwJuEJ1nu3qtHIFTRYAJO2VkrLE+/G0IxkpMeHL4kGJslMTpcuX5VhB/CRITXp/Y6ILHfqGGJPtQFwSMpRRhWKUMpVRa5zkCZ4uAG0Rn27ReMUfgJBHA2/EoU3PbdFx5SjJliiWTocFr6h91ZBNzRTrH28I3fZ7ecqqlOjYz3yi03vLSdZa88vlFaVubp3Vtf8zaGDniUTbGZINHzcPw4m2z0JMkVkOVAddc+YkjsFEELjdab6+2I+AIrBOBV1KLCbyXNLFbvO2L+V08PBxvdDz0t3QwJvgzV4w4NlFzDP8Qdac49huhX8tzo2soR56lLsowowiDhyWu+D/FDA++N2QGE56lvo3hiqoIufCNtJrtB0dgWwj4n7tuq729to7A6hHQZM53ccYm8dXX49AK1gbTlY4YWVmktHiUMCRTS3JZeXgiR+BcEPAls3NpSa+HI3A+CPT92ej51HA/NcHTE+8pGiuJ1/PB28kR2DwCbhBtvgs4AI7AuhCQt4JlL76tg/fCKQMBYWX7h95lJK+S1PiC86GXGXNV9HSOwEERcIPooHB7YY6AI5CJAJuFbc9Mpsimk9nbblM8RODre4c23W288iECbhCFaPi5I+AIrAIBeS14ZZ03o2xz8Cr0WrESbNhu/b3HkK41ruDrr9sPAeVxm0Lg/wpVaQzKjlJqAAAAAElFTkSuQmCC)
$\displaystyle - 0.037 \rho^{2} + \frac{0.042578305 \rho \left(- 0.000125 \rho^{3} + 0.0025 \rho^{2} + 0.05 \rho + 1\right)}{\left(1 - 0.05 \rho\right)^{3}}$
⎛ 3 2 ⎞
2 0.042578305⋅ρ⋅⎝- 0.000125⋅ρ + 0.0025⋅ρ + 0.05⋅ρ + 1⎠
- 0.037⋅ρ + ──────────────────────────────────────────────────────
3
(1 - 0.05⋅ρ)
%% Cell type:markdown id: tags:
Next, we use a function that determines the gas and liquid density, using the Maxwell construction rule. This function uses an iterative procedure, that terminates when the equal-area condition is fulfilled with a certain tolerance.
%% Cell type:code id: tags:
``` python
ρ_g, ρ_l = maxwell_construction(eos, tolerance=1e-3)
(ρ_g, ρ_l, ρ_l / ρ_g)
```
%%%% Output: execute_result
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAiQAAAAVCAYAAAB2SvVRAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAP4ElEQVR4Ae2di5XeRBKFsQ8BGBPBQgbGzsDOgEcEQAZwiMCHzYDdCMBkABkAzgAyACYD9n6aqnZLaqlva2bW9oz6HI1a3fW8Vf2Y1j//3Pvnn3/eqcu9e/c+UNsfddtZPxE4ETgROBE4ETgROBG4LgRae437tXARfKXnR3Xb21qXLys/1PYAEI74BJ+uL5BxhP+u8AifE/e7Euw3wM+bGJfXncNvAEz/NxNuO3Y3kW8jwXnd+kdsNWhZU9lzlHIvT0jU8bFan+j569IbFfV9G9U/df9Q17cjpygj/BVtqHzne+l6mQ/cRfNd9fxQ9c9Fc1G1QfO3ntk8JG9uJD5q0LJJSb8fq/4Xz7VeyQOfF7q2yoXo36NTtMj7MgjRyzOY/Rxt5SZa+r8pDZc2v2jRQiP6jEWyrPBZ0G3GbMTOkMlmAwxWGNJPkUwL99C9i3nIs3GfDKh+SAf2fiksMxZV79Wqkj0Ut5Y2ychY9mLUxQn5ktcdF2mHoztp8y6eXTyjfzc/Ulbe92Sqr+uPaIbyQ/QO5lYO44PkdeeOyteu7pDZ9TtlXtddfnTHdti2Syc5NnaO7WHXagyrfWj8ib6LvWi6sRTNaL51ZdY4VHZm82x+H9XvxCwV5V06Nse5+ix/gi7n3c01UHRsSFg7/zPpV4U7DL9RX1606/o424P2d915tbOiX7a5/MgL2qcpQ89Mbr9Xz9iJ7q+qNsBjEMzsCTracRAeEvJB8uVdbej9KZ+5q0BLpbaFtp90MVksL+RPGOmOjd8t5JHEyCs4Zv+SNvTj94xWz118KpndmEmeZWfQYQ8+Ixc/VjhWusFiF3f1W5gHFhbuqb++Sw+2sLnr5ukojeTOYgw/unTN4rYlV3ROjCycJMseF2FnV3fLbulZ4Rm67fxYyt2RiS5nnNv5IXmW32HTbg4HjlZ8XMyldyiOSyxHn0NfN3YuXfjZHf8jdkYsVmNY7fb4E2037qKxYim6kXyzZAZu0GJnveYQm7L+BZ2lX3zkUje2rVigE95ln9osf0RnrS0pX/T4Pa0p00StBxabL5Ig77TpmgFCn8oEStJt3Uf4RYtRZQIKPWwAymZBdQD+e6kPmpoueFeALvmSTryzzUy0MyEVXapzwrFa2ODVVQaH6uz4ZpsPPRMg2mabPj2Dbwt36Je0XXzCbitmkm/bmX5XPHsbki7ukkMcu5iHPxbuaWPew9bmwEqao3fJtuPW0hH83XElOgunoCu5mjrV3hoXVn6kjLxLFvmyi2fQkOeb+ZHyuAf9SqbaR8a5lR+SafuN/trOrXrY2c1jV3fIs+K4ZdPRdunO+WA3dj06fDhqw5IvdLXywx5/kmHFHbt1ObG08g1fXJlB687vtv7EU3ZYsQ07Nse560+lr/xyprbmGhg6idG0ht5XhfKpOi6PTC6f8+cnquQrj2zj/ouupzpuQclesfglhxMETjpmNsimZ1yVAuhaH7jFRseeSlSpPlXt94YvvF6pP3OCz61CguTRFP3YchEXz2xieG4VXn/V/rVoOBZ28YHfwlx0I3Y27bpCo4s5KlzciznCi1yaxaB0Xk/FituOKjdGLk4j48LVXcy/CTw7Mkf8cfNj2O8CwHbFjY+re8TvbatuQU8nP0bGn4u9G0s334iCJVO+jszvI/qHMqGDue2PCIfWFq2PrPufSv+D+wFGa5FPA/6isihJD+B7hX6HnwX9YmfhZlHOzU9L3p9hxOM9Yzb62Hj8saN70qv+H5f8somToud1u+h+1vUe92wPjHnkJKouJNfH6v+p8o9+5Na0XXxgimJhPmhnyr6uu4U5ymSnhfvCsM/EN9vcLvqv+ujGbUuPFSMxd3Gq8sYdF67u2vabwLMpc9Sfgfw44neNQavejU8wdXWP+t0y5pa1NfMjfBwZf13sQ6YVy4F8Q6wlU3T2/D6oP1yzb3uYI8TyRzaOrIFpHHuKT9/VD35DR9GsVANk1r54eLh4Lo+D/Gwk/hAPv9l+posNBrtgjv8m23S/UL+a3mnpfJ8OFV6flCJ6joLYUNBP33PJeal7KXpmB90q2MKCOKNPQslG3iP15wcOs2t2Fx0Dgg0GH8yaLZJ6/lH9LLjskP9WHVn4zWuqeiHu4iOeetPG41Zp4Qfvpp1bgrbaJWsXd/l2CHP0SfYu7urnyLHezG2Zebh9IG4rHbIvN9arvqphipGLk2TC2orrbFyM6E5bbgLPPZnyeXicp63cJXuVHwf93s1hdDnxcXVf1W/seVOKfO5it2er+HfHsLCy5k0Xe2xxYtmyWTpW+ZZ0AzKt+T3l1vc9/TVdr97DHP4Bf2bqJNtZW1jnn7EhAYzWBJ4T3MVM+vxhb3Id4U85j+V0WeDlCIs0f0GTizN3nFuWR9GQcnik/oN4L3iQHBKHVzO8BlptwKDJIhrkQV9syb7qziaDq1lCBrY+0cWm5tcWoWz5RLTgzyBGHjvFpdz0q4fPCOZSM+GCr107J2LvxyHcTcyxYBN3ySBmnLTlCZ5n8QEqM24tycMxqoVs4OSOiyHdN4GnKdP1p4Ym6638GPJbgg7lMAY04jOi+yp+p/+v+34Yu8DPGsPm+BvBfoVbI5YrGjW08q1FN7VtyHTn95bcIf0tAeaYbLG28r3Qha/u2sJnhZ7e1w/AaB33FsE7lfwNbIdkt+t9GZ3B4LRhdoIgzh90/bei+Rxpei6bknD6gnaVshBJFhuPbGd3Rx8bkdbmC966vNDDj+L5d92Ydelk0PBp6M2Njfpewq+L04Dvdf0mPk5CZiXasJOTEeQhm43TRKv7CD5i7ZZZzFw7u1KDQPKO4r6LOeKFRQ/31SmUa/coXS9uo/IW9LMYLfpaONnjYiGr9Vjrvgk8HZmH/DHyo+VvthW/r5DDyGrFJ3Vs3VP3Ib+3hL6O9itih8lOfjAXMD9uzpsDvif2LZbdWB7Mt5lMyTg8vx/U3/LTwrzFqLaZPzXN4NrCHuQDNiTsIgnssuxtUnLnmZ/dWPLyPMpfNhOVMD55TMA4xWFTgZ3/0sWpwle62B2ySP2ii9KScdlz+ZN+vowFnmZRHxsWPlOy9VoBPt759XRBNxXJ4jcfbH8h+fgzFdU5FWEB/1oXOp/pOfXWGzHoW/pm+IhmFHPklrJlZyE4XtnF3cQc7Zu4B5bOZvO4F8E5GLelvsMx2sJJcXPHha37JvB0ZQ74s8R2Kz9sv5cCq+fdHIZuIz627iv4XZn5Rla72AV+zIfdMRx55MybNvZL1DZiuSTbyrcl3fTckenM70u5Q/qXzDwHll3MN3idtXJiNdYW/H/wbksRbQwOGUu1LKA8RMm2FogTictf0TGpbpWygYBeRASiFNnJxoQy2aNn/tzxoWg/mlrXP9L+WU8EBz42BnuF3XnTd8l4BKNkvFwI+FXPnOxwsUGhYDcbrFLEx/tRTkumI6x4pr+Lj2jtmA3aWezbq0jmMO7iYRJyMEd1E3fJID/4c8VmTGC85tKNm/RljGeqR2JUM/ZwQq7od8eFq/sm8ByV6fhT4xP1Zn64fiNDdg7ncPA183hEN3IO+g3ray9HscPwwfywxt8o9gmgbGnGMvurezPfqv5S3ZJZ2did34uwVxVb/yuWV7VBzF8xqrblD0TqG1kDU+50MMKGhF1kc4FWe75CSKa85wnJ5iuLIHT5t+hSX2+hAQA+2ZtBfazn1u54slt0y80CIBLcD9WXJxQ5SJgkin7RgRUL4EqG2iicWsDLX9qkPTTNSshhEV3RoE/9LGg1zmVTNhN0+VDs0+MWlrUsuCw7L8XbP4dwl48W5mgPvLZwp/2JaDg+rAt5wYkY7ZxA7X0mqObbrIcdbty25Lgxmvil08ZpoRD/63FBt6ObOF43ntcRo5Y/+NTLD0gcv6EbymEYjPi4uhHXKpt+t4hfY9swdpWtVn6I/rmukfE3hL0Ry8nkmAe25qPKrSvlR8qp5/dh/SmkcbcwX86bBkZH1hbW1WnN5jeC2ReSyQAWYexnp9j6oh6OamZf3JU89d3lFx0TrljnXzymNvQU/UHHF5aRkGkjjmAsn0HJtq0vj4G3fNFaRc+AX/GojU97F13Qq3DCQaV8GVrKiX50rL4cSG20w1fkRdvqy3hCDnGZ+nS38Ak+K2ah27Iz/RMPeMx8yL7QvcIw2le4S46NecjYxb22I+s7PpYYJO3IPeR247YlU/xWjOBX6eIkGvIDjItf1HUhoIyLkGfrXtovWc3cTjr17+ZH0tX3lky12f6kLPHs5of6Lb9FZ+cwulWc+Li6h/yW7hLvxOHoXbKs2O3RqW8IO8dWyVzlXLRZ40+0FvZuLNNmyd3Nt4rOyQ97fq/kWvrDLyu2lewV5lWf40+TP+LWXD8iTqx502uD1cJUGcCrg+U3rqGwTHSq5+S32qSor8sfoGFMSeiQiZ5aN8AirwxE1eGbfdupntn5zTYMeoYXeYU39EKLTDY/y6v1bZpT8oi22JpYhTwGwNKe5Fm2k1TT50AaMmabRNF18UkZou1ijo1cyRO2N+1MGtFzVEpCbU0GFu7w67IxX9jWxD1trO9h62wDqrbNXK159+qSYcVtT5f6nBhZOEmWNS7SJ0d30tZ38RH7GZ6L/t38qGmz3pKptiF/kKWSubuZH6JxMb/2ucPUbfsteVfO44xB4GfFTno36dTnjn/bdslc5ZzarPGX/onejbs9J0mmk2/W+A387fk96Lv6K/83Y5Y09b2Feei0/BH/kbWFtfdbPiTCjocNCZ9bWJU4msKhC11/6nqi67noZ68sREcwOR5evscm+br8ouHoEzroKQ91tfRAQ0k6bOdYblYkC/DyiB5ZHAfx4VH8KCXshrZV+EuZ2edQQi6biM/V1/ycgGgYMOXVj+rIZ6Js2Qn+3+iajqt0p2z51MUHZum3MHftFF2+CsEvZBN7jhFZnGZ/GSXaLu6iIVdszEWLT9Dv4g4dRbQkN/TYSyFOv8jW6a+mQj/tm/8kkM69IhlW3ELXoXERvBZOorXGBT6J1sqP9F/0PTzt/BiQafsTPnXzw/VbdMi61rljQLftt2QyjihXyWMrdtLl0o1gtxoXl+5YY9gaf8hzsA8srbEWMp18G5rnZIM1vw/ot2KGPIr098a57Y9k2Wtg6Eb2l9N/+41gfKIJe7bJgPAsJwK3EQHlPL9hLD9fcRtdPX26xQiceXyLg3tHXFMOs7njl9sP74fP7MxmJxt3BIvTzbuLwBMNgNlp2d2F4vT8LUbgzOO3OHin6RMCnEZOJ4PThkQTM8fuu9/PcQJ3InBbENCOnFcWvH48y4nAW4vAmcdvbehOwwOBOB3hM4nTq/88IaGbzzzwDuksJwK3HQE+zNv8Ft7b7vjp361C4MzjWxXOO+kMe47ydmb6DEnCELsV/qrlnKwTlPN+InAicCJwInAicCJwrQhov8FflfEvWsr3rPwPu9hiX+CPL58AAAAASUVORK5CYII=)
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAiQAAAAVCAYAAAB2SvVRAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAQFUlEQVR4Ae2d65UdtRKFD7MIwJgIgAywnYGdAY8IDBnA4hf884IMgAh4ZABEADgDuBFgJgPf/WlUfdVqdWur58y1PW6t1UevUj12VbdqdHpm3nj+/PmpLF999dW7uv4ux472gcCBwIHAgcCBwIHAgcC5EGjlGhclcxF8pv775dir2pYtCzs0dgcQ9tjEOl2fwGPP+tdljfA5cH9dnP0S2HkT9+W5Y/glgOn/psJtx+4m4m3EOS9a/oiuBi17KjnHVN6MliY+UPtt1d/EWNQa+zq3/1H9nq6vNWafooysL2hD/A8aexodavW/Lfp31X6ssctijOZvGiN5iLWRSNyr6OBHkvJ5Hr+v+hl9jcdapthokfutxlUtyqXG32JUNfw+zRTIpQ9mv+axqdIY819MA6cT/Z9atNBoPHwRSxb4VHSrPhMvW8/MEwx+0nVPay8ZaxQL9yy7hznsbdxrXSSDtZ+qDl/UJLv74jnkt5Yg8Qhf9nzk4HQSP+e+SKo4smudtWYTzzzfi48Z2y2emnPsGYoP8exiLgWtGMYQ8XOeHclmUzY8HbsTz3N9SCY4dn1n0NnYObpneYt7WOND95/ou34XjePL0XhzeE5QFHrGWP18H5IPE/G0fBsCM/0C88zLskc8oIvnLr6iP9sDRfMrsnTxg/538E8JiTrJuapbm/WfonuiuZ9ZkGn/VP1IVzcpEY21XnQozA1BIpA2btX0uUiCTuqjJ/xIClLipBqw/6OaTbLUh6SCwjzj6I8dl6qnoj5y4fcoBtUmeMPGSCIeaIx2KSOWPFQjbRpai47YEM44qU2y94vqD3UlHGOhapw00TKuPgkJpzkTrdpdfIKnaLuYw1/0XT0z3feiBU+SNfTYKl3csy0O5sixcF9RiNgpk8oVsl3Dlt/WOJs+smIz+8i9L06O7BW9F3hm2SPxUbNe4+naY8fHgN3dGMYI8bP8k2nde9K1u8ZxuO/6zqXLCljYDSi7iI+81r7/pL+DvevLkXhzeZ6kI7TYurr/Zbst+YM+y6ynqol51rH73M6yu3sL0kT7jS722h91XcYJCRtwmZUnzUTwiRr1xshJABsl9NMmnhZUH4PrAeEHrYkEAG5smmUCwIPvrmimUxy1n+r6Q+O1Pox/qPFewfY6IQBMbEendOoBE40t7NUYgcRcJA+sI+MjAYmxsOkLzcUYa6DlZqnLYw38pmuiVdvBJ3g6PrP0lI4kcAlHtZ2v9BzcbcwlG5sc3CGdStZ16p+zId4jfluIzusdH7k42ffFgOyZ3mt4anw0Pia+azxFYNsDM/Hpxodo8JmDOSydGIbO8s+A7CG7UeA6RXpZvnPpsi4udl3VJXd2nB8LMp7WczPTOn63fIkO4tmNt6yrzVP01vPdlS8dLd9mPadK65qYZwLXHmtvmYRe7d2J90Ue/EiKpCOTgogmG1HrJ8zfNf5Qa+5AtFGs9eLDCQInGTMdNM4pTOl86MoEJUSjo6NP0Jc1pxt/NWwhiSCQU8KhNja3Sp2powvBwJWKeEztGMs1Jz+lfdX0VVfrXXxYYGEuuhE9rxQ536eLORJd3CfthBexNPPBNHmehuW3DVGuj1ycRu4LV/ak/k3g2eE5Yo8bH8N2TwCsN1z/uLJH7F7X6hbMdOJj5P5zsXd96cYbXrB4ytaR5/uI/KFI6GBu2yPCob1Fctn3yUHuXGQwWpt8KPCMRlWCHsC3CvPOek4oLqXL2sZ9QtksqMXvnzx3P9cjVfoaZkN2kqv58rQi8dcYWd2TUpjG+F7sLeoYV5uAo9SnUATXB5rnNCXsgw6+JW0XHxblYmGOfrpcPYP3uWoLc4RJRwv3SrGPtW6W3Fbz1+26fluTY/lIi7s4yc6IG/e+cGWXut8Enk2eo/YMxMceu0sMWu2uf/KiruxRu1vK3LKxZnxkG0fuvy72mafly4F4g63FU3T2831QfjbNrrYwh4llj3Tcs7eQU3z0pj74CR1Bs1LcILPxqnO36k/dwfUkEn9rDT/ZfqyLBIMseHq5U3MkLBo6tWS+zYRKnGakjujTMa06zDPHOyRP02T+UJ8MulXQ5VTTB6HG4ccLOfHCYUzNas1zQ6TjKLVnm6T6P3NpnoTlX7Xhhd3l1z3qpnc3NvGBSOtjc6K7Vlr4sXZVzzVGa+PitYm75ndhjjyt3cRd8xw5lsncmpq7xyXD9dtChtbaPhKthZPokNPy6+y+GJEdimvN2fHc4qm54fs8dKXW+kV8aMzGPHhpzWYMZ1ld/7iyRXctu0Pvl6F2sNvSU+s3Y07z1v0nOtvvou36sqWz1i3iLegGeHb3v+BZ11vya9qtvvhsYs7aAXtmorTO2VvIQR5d6AMw/ppxuOrEA+6yMRdDWw4fWR987kv5z3XxogtZIwkJm3UUNm8CoC7v54HgQ5c2L8rAi42ei5dnAGeziAZ+yNlKNkgyuJoFHrpwMnaQBP3RIhQNN0IkKvBDv1nSpH7Y1cNnBHOxTUFm6ZmIvY9duJuYo8Eq7uKBz3iwxwmep/EOKtNvLc7DPiqZrODk3hdDsm8CT5Ona08JTbRb8TFktxjtimEUaPhnRPZ17A77X3S9G7uMn3UPC2fnuTmC/QK3hi8XNBpoxVuLLo2t8HSf7y2+Q/JbDKSThfnK2vc13twrsVVXdw/MfMlB3r3QB2C0jnsz3WYVP4FtEm1M8mvG4QyUj405lvyoxvcFzWMm1J+SCrUB5JJxlWkj0jjvn8T4SW3myMKcn55/Eh1Z+PTyrPpT0TgOeKh6cbIURJrj5S6SIW6cH3SRDJXJVSLNY+jJyQj84M07LYlW9Qg+WtotM5+Jv6Vnl2smEL+9uG9iDnvx7uHOr6rVMeSqPkQnOfhn1W9DzJbEMx9V0y2c7Pui4tXqlrJvAk+H5y57jPho2Rtjk93iszeG4dXyT8hYq0P2LrvXmL6I8Wtih8pOfJwk51z3X2DfgmvTl9Kh9zzq8hSP3c/3nfJbOlmYtxZqbBUj6Teyt5CDpISELJIHa122kpTIPOPdjXot/dH1UzJRMONNahx2nzEZiJ7v6OLXZz/TRXZIUPyui9LicTVz9ck8f4yFNc2iORIWvh5ZO8JjHacePVnQpSJe/OSD7unXea9Gkz0cC/Pw41QImY/UD7llIsaSlrwZPqIZxRy+U1nTcyLY39jEXXIdzJG+irt4gKWTbO63Iq/Msly/1fJ2+2gNJ42794Ut+ybwdHkO2FNjuxYftt01w6K/GcPQSe9WHNuyr2F3oeZL2exil/Gz7mHh5D43bexr1FZ8WZOtxVtNl/odns7zveY7JL9eTD9jueu52bFnJk60zT2wIML+OxfFwKwpBjzkKHeuqtlnjLVATITu+oIu5M0E5c6UQECvi4wufRWjGkMjy036aIx3MNis10roP5vXGgL9rmoSg61Cdt4Mdq3lpIdTm7r8kQem0x31SahmXwtpLfZwWoKOD9UPXKLW8KIkfAraln0xFhiN6LkQ2BqQ/GHctcbFHJFN3MUD+++oXo3Hlr7XGOv6bY23dAw/hj9K0hhb2NHDCb66eveFJVt8zo7nKE/HnhK43F6LD8tueEjucAzndc04xo6sW/g2d1MVY5O/d9pd8nxh7b3YobDWjsScdf+BZQYjcC6xibEJ+5jUuqYvY76om/FWzE/NNZ6FjqHrtKZopOd70Y+mLT8WlLVkj2BeLj2t2QOR5vbsLelg5E2tZ2MN58CvLPEVQjlGO05IVr+yyAvc9Wt0IXcRNDGRaxIA3uwNp95Xv5UwJL1F97RaD4g49z3VcULBWAoE1ZN8tcGK8QWPzDMlQqLjN1hCnzz1vyrzYRNd0GiM0xISkxLntaCE6aSf2mtYlrxYY+kJ4UAZwl02WpgjX7RbuIPNA9FwfFgW4oITMcbBdJb8lYRuO+vh+m2NreujtF4ybZwqgdhf3hdMO7Lx47nxPIePWvZg0yn7Zeu+dOyG1VAMZ9k9/7iyYdcqq3a3iF/g2DB2ha5WfIj+ia6R+28Ie/deM+JtMs3guaZj8Cif72lsRH4wadQW5pI1e24a9uzZW3i+PyMhwVgUaxUe5GSjdbmnAb4fWmymFaG7niOjejOBFXL4yQ+HnTIQ/PGgd0K2agzh1AHaKN9pfAZinoAu8QpCatFyw/MArtfwoKnfSeCmo7QSHsbBpN4EGI91YQt2cbFhLgJO9NgVulr4IETFxdzS84ql/WnjLptHMEeBwG+Bu3iBU2A1KavxfxlXPSWZTKrffKBNCzcaWjvitzVOro/QtYuTaIhT575AH0f2EJ5rRpbj0tHmOWhPiFmNj0zg2A2pHcMQS9euf0RmyR61W/S74xjdb6AMYVfKly0j8eE+NxFhYQ+hdHB8CSmlF2+JyOQ58nxPfF35QdyqRzCP9aY9e/aWu5Lx94U++En/QQgsawlnM36mmgdeKmqzUX6kK72AxSBjup7rmn1Nor67ntMANo4p+VF7IUc0JE71hkTAcVSNHVH487az78XU521fSr05wRMe2JDWRa0x+AJuWdCLUo9fjV59BfNLdKjFA/xYV/NDl9l7JepDz5EhR8cpUVHt4sNaC3PxJ/ly9RRpKvHVGMHTKhbu0nEUc2T1cG/pw5pYl+Ylmz6/Yj2L1dbijTHXb9e9L1yc3PviJLvd+GiZv8CzIurFR0Weui2etj0Fw/DzZTE2NQfstmIYxuJp+WdAtm23eGLvdeN4wkcN13dbdBZ26K5rsV+UyhTtVnxY9x88JMeKd9FZvqz0otuMNyZcnqKzn+/wzQVcKKvyr6bT55bPCrKp2cLctkdc9uwt97Tu6RtffvklWSGbIu8tLIrGUY5EAcP/0fVA1xONlwkAyvJrOyQVvGgzFfWt9SwQLXKgp7DpteRE0hJ06N766ZgAixMPeJHI8PIodkxF/fTrRtPAvMEpEEBNRX34spk9VptAWhSNcxLDTROFNfxF15ae4P+FrjLRWrOpiw8CJcfC3NVTdCRsFOyCN74nWSJp4oafivpd3EUzhDnMM99N3EMJ0ZKMogf6UvDT7xqP/3+EfMrWPwm8olj5FC/Lb6LbfV/ktdjRKrPYFK11X8BItFZ8hFDR9/C042OAp21PtgmcNuPDtVt03RjOMu04HpBt2y2e54hjy3eS5dKNYLfYLwbiw7r/sp+68Z6xtO61zNOJNzs+Mk/r+T4g3/IZ/CjCoHef2/aIl70HZtnw/vSN58+fowgdfnNllmRAeJQDgduIgGKdUyseiJe30b7DptcDgSOOXw8/32YrFcMkd/xw+95FNpTMbHaycZsBOGw7EBACvDN0JCNHKLzqCBxx/Kp78NCfbzLSyWBKSPRg5th98+9zHJgdCNwWBBTvHOHy9eNRDgReWQSOOH5lXXconhFQDHM6Qu6RvvqPExKmeeeB75COciBw2xH4RDdA86/w3nbDD/tuFQJHHN8qd76WxpBzTN/OpHdIAoacrfDfZ4+HdYBy1AcCBwIHAgcCBwIHAmdFQHkGv/nKv2iZ/uzFfwEZmXhdZ0r+4QAAAABJRU5ErkJggg==)
$\displaystyle \left( 0.0695273860315274, \ 8.02904149705209, \ 115.480272671423\right)$
(0.0695273860315274, 8.02904149705209, 115.480272671423)
%% Cell type:markdown id: tags:
With this information we call a function that assembles the full free energy density:
%% Cell type:code id: tags:
``` python
free_energy = free_energy_high_density_ratio(eos, ρ, ρ_g, ρ_l, c_l1, c_l2, λ, κ)
free_energy
```
%%%% Output: execute_result
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAABO4AAABMCAYAAADa4/hFAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAgAElEQVR4Ae2difHltNLFYWoCYIngQQYwE8GDDAaIAMgAaiKg+DIAImDJAF4EwGQAGTBMBnzn579kZF/5uuXlrkdVvrbl7lbrqLW1Jd/X//nnn9ccjIARMAJG4LIReP3113+Whl+qzX5x2Zpau0tAQPbyjmzlz0vQxToYASNwfwi4Dbq/MneOjYARMAK3joD6tjeUxx91fKRx9qtT5vfRKRNzWkbACBgBI9COgDoJOogf7bRrx+4eOWQvXyjf722VdybgOj7TwWDF4UoQUHkd2ABlqOOdK8mC1bxuBGg3aIscjIARMAJGwAjcBALJWfelMvM7Y6pTZup1JX7K9JyWETACRsAINCCgTgGn3Z9qq+kkHK4EAZUbThPK7v3UyYc1F+/Xifgvnd/V8bVkhFbPifeZ6J9O2csSvZJM8jIVXim9N/PDiP6iwXn0eeJh4MM9+fwlxQ1OEb2jMhNdrk9PlNBLHdXVrC20A4Uv4Ea6/y01wDav0s0DzKpNin53m1Ua3xTQvKXrT1Xmr4q41xLmLbYxKzPLl+xct3LU90o/45Pj0KGV7mhdLeQdpUOBKO3WdC1pJ9pZe5GOOO5oH76Fx8EIGAEjYASMwC0goP6N8TbjVsbppwlK7DUfxsA2YBuwDVyeDagX+EzH7y6byyubWpmorHCM4ODCkfC7Dt6MvVGjnYpLfM/yc/h1/KGDra9H++tEe2AvKX6xXuLHicFWbfI1PtCt1Jd8l/cH+us5cd+U+YFHB3iNeUN6iy8qEwfhz6O0yR9pfzCKD9OWfJdyrfxQNjjvyBvX5HNgj9zrCGE8lS/xR8scHb7IcnSN4wf9etvWdagckZFoZ2UmWsoSPfsy1jX5/iPr00KXaGfz3ULXQpvyMq4rYNFj2SIvSiv5zfaSdB3YHen5MAa2AduAbcA2cM02oP6NcRU7ok7Sp50kkVNlxumcxmiMs3G2DexvA+oImGgOJrXGfX/ct8JYZcdKE7wL4QmraHHUDhwJ6KPQOc7mdBMdTrXPjtHp+RK9eKN4MF6QLGy0d8CRto5Z/UWTdRg7HsjsgeMxYZB5qnhGZYoOZ83AuZHkU9f+LvPZQlvyXco1+rfoIvqjGNdkiSda5uA+wBd5CjiEe0eqrrMOs7Yh2pDMlA5Ott5pWEu7kS6a7xBdSjtEq3xsSteSNrT5KMqqWicLOvTt24kc7/O/WBoLY2EbsA3YBq7TBtS/8dLs6Nh7q7J9pIQcjIARMAJG4PIQYEL7lRr70BbJy1PfGi1A4CPxHGzbU9yvOj7QsnxWuxwLH8te9tiSRvq1gEMvb2nkeVR/8vgqHfDhDOB+TYjK/ECJ/FHBki264++/tdCu0f2aeaNlzorKWltGuZW2HS1HMAvJTNtZWN03qBuyuQ85EESI0j1Qh209ig9io7Rb07WknbIfPwljcP+4UufiQkxpBIyAETACRuAyEeDTK1+foo+z4+4yDcBaGQEjcMcIqPFnhcJbmvD83x3DcI9Zx1H0spLx7PDgeTUkp0Omq9IsjZQd/jTmVXqsAvxqFB/SX/J+0fEm58yf9OeWVYPNoUEmafLNyClHYekcbaFt1vlGGGbLXGWbMa3Z9l8Jhyeco+XYIlNicS6/OlLmJE2I0kE7m2+IFKJ0LbRRmVG6lrShXRJomz5ewmgeI2AEjIARMAKXioDGFoyR6eMYF+8aHu8q3cKNgBEwAkZgCQI0/oPVIUuEmOd6ECgcEceU5mP+U4GVQ70jbIpoi3jpyhbZ9zRYyX/wwGql7Jw5lkRVf/HiYMDmP5fMTex+Sqbks1qpFliRheOIFV9daKHNPJxT2qSDcxCs+OOLXZyqkn00SBdeAlA2b+tAF1bx9nnU/eIg2aEyV3qvREs6tfJHLwK6HYSE5YFtNMrEKfinZFHGn+jAWcjHpNlKXNaZEF0031E66RGuP1GZUbqWtKFdEcCZNmqT+r1CD7MaASNgBIyAEdgaAcbDP6vvZSfKbuM9O+62LjbLMwJGwAisQECNPt94YkI8Xs20QqpZrwCB7NTA2TMVjjlKcDosWq02ldiReBwpHGVo1l+2jiMFp91THTiTftOxKiyRmXg6B9tc4nO0eg4ufPOL1VvZicf31fp/3SX+RAF7+QEnF+lJN/LINmG2iJYOKx4vCS1lzhtpynocsAHCwLYTznO2EZWZZT9Rvktn899Kh3+1RQ4hShfNd5SOtKO0W9O1pA3t0sA3gGrlv1Se+YyAETACRsAIXAQCjKk0nmAcyxhw6gXxal0frZZgAUbACBgBI7AlAkz4f1In0E22txRsWVePQF6dVMsIToeXtQdbxmlggvOHf+Zc4vgZ6C8ZL3T8nw4GOd/r+F3yn63Rd6FM/uSAOhfZmj5Jm3TnTxU6p13KB85Ivp13cqeF9MBB17cjuuYtMOV2KgcvEOQy/5SbEgddv6eorN/gDbV0jdjGrEylQb0gsEJ0vNrrB8V/B02UrpMU+8n5nqOO0iEnSrs1XUvaU/mlbaLtcDACRsAIGAEjcIsIsODiWTGe2DyPjzeXaIFXiUAysudJ+Ty44k10HlRfZb6s9OUhYFubLhNhw0SW+rfb25rp1P3kzAgcc7rlVTb5e2A1VaE5RXuNU2rgZEnKrNJffc1Psn/0/1Fnvn+3Oi8RmUoLJxbfvJutcwHa7ySrX9GVcMn96Rvp/uAkuTz7n45JmgMmtRHS+UUlfi6KsuPPIPhn3Vo5zvGXz8NlrrTYLvsfMfMBZ9o5nEv86QkHztpJXcRbtY1GmTX5rIRkKzGrVXGwEubocHxG8x2lI90o7dZ0LWlDuzSAa4t9L03HfEbACBgBI2AEzoFAfqHNuCLyIrhZRzvumiG7WYbBvwNqYM1khkEt34FxMAJbImBbm0YT5zkfUX8xTeInt4hAckKQtdrkNsfVnAqnhqPqZGnRPzluBt+TS5nAecLKNI68fTE9On5aIlM8+U9g+PbW0TBHq+foTDmxiqsMxBMm6zTY6fn7HdVGP9KHf6XmD26m5GabWpxiS5mTSMpnuRqRFXh5y3Vn27rHqRe2jTmZhY5gPBVwYrLNhedH6SAoZNYwzHHdH6AkmTkO9hxyXJfvrWVG5aFMC21WfsH5VC8WFqhmFiNgBIyAETAC6xBIfSnOO8Y5uzjuHq1TcXtuDXLYssAbWQ7evHPkAc72CVpiRuAz4ZwnGMQxmH5Hcd0gOhPdyhmb0mE7O0+B3pWtNUJMHcxvbBpZTX4DCFD2eYVWmR0mvYRjtsFKnF37StpNpYF+r3TUQlR/XgqxLXZLfZtkKm0ckO9qoNWvtFMcfd4B/kFa5LC9c4wNf4aAE+fUTldWkdXw7WxJ+kw6EsXXEqJlPiWTMQb/MpxxayrHCaFjmehYwyKz57KJ0sEXzXeUbg+Ze6SdMWs9gz9tlIMRMAJGwAgYgVtFgM+pVMeSW2T44hx3yhSrcfgHNg4GwnT0DOQc9kUA73DeKrJvSpch3XZ2vnK4N1sLIS3nAJNNJjeslHG4TwTo8HG4jAOrpmpOoZIO58OB06kk2OA66zY1AY/qj5Om9h3HLB+HQ2sIy0x17SnjjFEiz3Q/yFsDLU73QR+aeKnXvXNwlN6et98qf7UV81u/HAiVubDguy/8GQRtXBfSNfqU5dBSjlGZ7CCo1Q3qFSucs71F6dA/lO8Guj1kRnVsSRvaJeEtMWUH6RJ+8xgBI2AEjIARuHQE8q4Lxjabh0t03O26GkcDxdrgbXNgL03gXL41cGWQz4A5B5wrrBJ4kSNu7Gw726lAbWuLgc2NfJ5ELhZkxotA4O2kBRPWQcBhoeMfHYOXUrTDInypeBxIXYBWFx/r6D7E/xBb/aWtflp9Moxs1qtgz06Xsq/oHzfoj6Nm4KBOeUb+56O+KMuf1DsRhGQqHcYAODQog2/KQ3GDtKO0okNv5PZlneJIB5nn6Ee7vCn9Pkgn/rGaMOVInMSY/OhYY7PgM3CK6r6GT6gcyYRCSKbwZ9s1W2Hztly26FJmg3oVpSNh0YbqapRuD5l7pI2eRZi0l4ImX76vi3PUg5y+z0ZgcwTUjtAGXWy4dP0uFjgrZgQWIqB+l/Exx+wnWJYk8boSWMJ3wKPGgX/sGnTKaWDEN1bCb9nEwzdnfkgZZ3BFo8jfyL8/ln+gxEyEZDFoxRnV9O2cGbFX8Vh5xylAGc3uuRYtKwQYUIM5xncxAd3GdqA4BuC2swsppVuxtVPDKdxwZDyRfb956rTPnZ7yTjtPG9W3/efWaWn6ygttJ4H80DbRL9IH/qyy7f/VUnT0a2wR5CVJHxQPDw4G2t6/dOCM+2rc7iluEMTXtduiq62yoi9dpReJSQblhLORPy6q9qNR/UUHPqUDCdmshB44rqN6ixf9ZmWKBtxJqxZY1YiDoQtRWtHhaAVfsM9OV66/kbzBuEhxJwvSi3zm1Ww4FV9yL52wrT6IbpVtiD9ks6LLjjPoCT+Oy5tI0c2WI3SEqMyCNqcNHtV6lWRG6KL5DtElHUO00nFTusa0Q/aCzBykL/UOJ/agfufnPs8jIAxnx7+ioc5Tfy6+L43kZx6V81FI/4ufU0pHbCE09zsfktulfO02tR0S1y3pnOW4RdqSwXyOb+dWx+NrSmfguFNCeVDFZIHEGESHnG7i/Vv0DCTyIJVrwirnT9Lp2drMSw6D6drWmE5JPb9YZ1Wn4AY/yiOdDNtC+snjWKxo6PTZLsI/1r0aP9/iXmnYzlbWiy3KYU8ZN2Jr1IXsVKE94/7AsbAVjsKMNvQ31btd3tLU9FxTF5EX4RfNLI6iyY6PmprE0W7dnUNzCoypeOHI5Hjpv41OiXX8DALCnT4Tp3vv9Jth8WMjcDcIpD6AFxebT2KWgCh9Fo9BSW8Jv3iYY+C4zGOKTvWETXauP1Fkdq7nuVRHx49oZ+dZogn3paJlXPO8T+BhDjflTJ/tx7McyQ3hG8lPlslZ9GD4nQ504aVX+eJHUacL0mW3OWUUv2huJW927heVtYZubb4i/KKZrSPkQXRhe070UZtmLJADL4d4yTmYS6e0o3UePUO0JCrZYV+GaHOeYCV8L10P2p2HR73sgzYsP8/npMNRujk99TxUjqQp2lBZJrpZLFvSznkenyUDbL8Qnt0/Xo2fr7qXUJb8I4O36DjI8j0NevdmOscdOydagEYYfCj9xjGeuWfixwCRtVYOefl9nB5ydfDmkIpG/tF9VVrjNC7tPuWzmkc9w/hZIZBtgHs8xt39FueUvu3sxu0MW0llfZW2Jt1pG/q6kPLDQI2M9fa7RZ1IskkP2YM0t5Jfk6O01rb5s/xKI4Sj6OgveENFWzw+6AM2x7yGybXHCSdWrJ/Mhq4dr630F+bdmGcreZaz3ZjDWJ4fS9ok2qZLKAvpMdtvHdNzKb/4aCNwivXjad0zxsahWcbRFzIW+KCM51oBGUfnWXoe7ktFe9BXKI450aC/1X2oH086hvGV3Nn8jDEo0hhgWaPbKy7hscucUrLD+LXkL8mtjsdb5CylXZuvKL/oZm1KNGF7Jr+RtJNM0sZZk+fQ+DCor/08mmsdoTofpRUd+Qn7MpJc7KxvYxL/H1n32lk05G+23k3RKT6sZ5JxtK1LZRMqS8lrwX3Whmr4lHFKj7F4tR0v6ZZcZ+MigYMCU1zXAUQEi3a2MCNyMo3kdSDrvLqhkYzZQYNoeCMByKvTy3m4xLPyV53YKT7jTUOTD3CbxUM0LMWeHZRBo8N2dgd2hu2n8q4NDC/e1qR7bg/6Aazi6CDI2MGAbW1dl0zqHLL7Tn+tzGP8SmdVXYzykx8d5OsojnrOSsZ+8pKvFY+tHNhQfu5zFbO8RP8AT+N1iNdaTGSfuV3obXytTPNvX07G9DyYpjZ8MFFdWxaSGRpzjtMR30n6vUq69IMHk17FMdnuJ/WZT3FMWP/O9/kMfb6eOosm1JeKDiwOxu2Koz0bjHF0H+3Hm/CV3Nn81PIpPvq4Rbw1ea1xSnuXOaXkNuHXoneSfZax1Np8tfCLdtYuRBOyZ/CNpk26Omp1Flvt279EF6rzLbTZFsST8zY5dxcNTrvBXEP3Az2zvHxOcg/asPw8nxvpALiqp+Jny5E0U3qROUa4rY2mnfNcO0vG4jmdeI/2b49EQGC5cW155K+K/0BL/mjITxbyckaB8aGOV9ynuKU6fCw5k9tDlwq9Rr6Ew8eVMu2878oT53zQqb8K5BP7iNiI7SwA5q2QXLmt0R5i+739B+vC0uLDQUUIfZrggXTV79q6GOWP4khfUwtMQgZbi2pEjhsgQNkwsXA4DQJsbyP4210POPjXCJQI0BZt3YZHx5ylHlxH+60xX75v5tdYmwncYCyRhenMBO2PynictoQ/g8njgoJl9jLal7JtOfpZjmg/3ozPbG4uk2CvOeVu+B0Zj58C4bX5Wss/zmPUnuGLps2OnNr4nbRKP0pLnW+hRdfZoDYFPWmTBj4R2Qf+lmp7MNOG9WlG6XqGbS6iZbk5ljPqZ1tY8nmIo/1bdtyRoZcVJXLCPD9JSB0VHT3/iPZeMgT2JNf0m9VJ/FOVaZb3hgko14/L/KnCvqmDbx4OjpJmg2vb2QYgXpmIq7Q11QO+n0Kd6CfjqS0B/j2cIm+lcs1t7t7FvLYuhvijOIru4I8OhDcrvr/aG4hbky8smSR+Lvx48+qwPwK/KYnuJeP+STkFI3A9CKQ2iG8dnapfmwMn1G8dEbKE/xPlfzBJLuQzvuAP82iza4EJXFOQrGhfioPvmcroZx1lOvS7gzGOZEbHQ0vwacrfnsTCYdZRKpo955R743cwHt8Tz0L22nyt5S9U0RaQuD3DN5t2UX9qfoq/UuJP0rmlzrfQJvGzJ16ivDrS5tQEHGvDSvooXcmz6rqhLPfAclL3At88t5ukbX3wuDC4Y7yhhCWLZb50APxFPA3gV1Ieb+hBEC2dQzZoOi0+yg4tq72QwbkPerb0jR0eZArsasIRbLbKA3iAy9RgYqt0ejm2sx6Ki7qwrcWKQzjRedNmMQmZrDeJjjd0tGm0gfyDY2TSQptHgG/XsLYuruFvwBHs+Ce0/CHZSUxWYD4p89ofJJub/Qfxa8/nJegvrKmzVzXGuATcrMPtI6C6cTFt0Jp+i5Jawi8eXp4MnGBlqQsfxgq1wIoYHAwH8yfJDM+zkCH6al8q2T/pGU4+HFF/65q+ltUhbOs7cP4pvg+iPRgPKS6PYXq6ysXBXFJ8TfmpyOyjkg6M09jSRyA/bLkbtM9Jf+aUjM3Ah22C6P+nnpH/Y6tkdplTJt2lwtFwgB/UKT+RcefVzf2W4CKeJptK+B2M76Npy17YGUhR1MoHfwgBO6NOh+t8Cy2ygwEHInZOG/OJDvww1XqieGzraBsGDSFK90Ad+5XMpnJMehy0TcS3Yrkk7YlcRdrFCdZ69GNFZ0Nj8DkVIglD0//9uDKNkbIEnDfRfaOZ4mkkmdB2nYPi6NgwZmi3/udAjHSy49Sziwlz2GyoKJ0axn3KYDs7JdozadnWQlu76Yzo3KgrT3UwiP5NRzWIlo6f7zV0Lxl0Dx8vICJtWu7cX1aFbxu5ti4287fgmLIKlhxHw0rMj8r2QyNgBIyAEbgZBJr7rVHOm/jVNzEHYmVL5MVdn1TqK+GtvbQKzbN6YQ8Xk32pdPtI6TE/YoIMHbpO9rtJt6nxUBM+SoewJD8PnKPfpNuPimYe2WOueFYU8hKwcyInOuag7Kh4hRjF8U3Bbk6q6y6O+Imw15xyCX7oTnlFx53XOPdrxSVsU8kWpuyZ4m9JG38GssbhvRSBXtWQ9Jiq8wOeFtoB4783WY8nsv++jZFcnPefKq532us+1IZF6f5VIXQVLkekSQdwPlaWB4kmnhruTWkfCP434pUusw39G7vy6lGQP08qJ8lV2DSWKNkFXdNw4rAbO81oMH8qjUP3bNskfo9AAZxiMryF7kexkZHxzYsfdXRfZCwTTM8+03mwUrGkKa7BA2O9tGA7O12JNNtasrGvdebADjmoX8fC1doamVI79ULH/+ngxcL3On5XnnlDPQgpjo/TlyuDcfJRZ2ud+YAfunHEme9n6+KMfgP+KI7IFF60TfzbVf/Cp5bWBpjXxDrOCBgBI2AE7hOBQb+1AIKS/+jq/COycT4xRzpYrai46DyrEz/Xl6Y+lHkbK27ob+l7WXBxMMZRfHg8BO1EKPFBXlN+JmTm6Ixb77RLD3BOMGbNDpTnuh9vTybvnRNjbtwhunPOKQf4pXJqGXde9Xhc2E+FHpcWmxJtaHw/lWiKz2l/yr3KpB/vJ5ujfhHGdvkQ+/CbbfegzpdE6bqFdsAuffI8A0f2eOfQDyL+rqCBN9qGRekG+hy7aSlH5CwsyyqWrWkfyQf1LWN+hKzt0WORI3gqZE8hSymXBAyVjzLyDyoszcSJR8fwVSlMzyIrUkqWlmvykCtOC99R2mTc/xNRS6F8hHHVBEewES/LcXEeDJxuiqNDyo1FRB/KJUI3UDXpmNMpn3V2ouel0yI/p2HE6WE7y4g0nIUp5bSZnZF0Ksej9XDC1gZ/FJDk4Cg+tq3gZmxNmLC15JXyi8Oyf1sLpgrf6ejfXnUx/9bTSF3r6hC4J96jJ6W/xi7W1sVV/DM4km/akWODnIzNIsxXYpfTPjhL7sELlQMiRxgBI2AEjECHgPqCbn/ZFBxqU5kzLB1zjsWu6rckLMwvvVnBNl60MNbn4D7lF6cSY+ZoGMyzRkyTfWnS8X2llcftHyoOhx2TWSbwfNducjyiZ4PxkHjC+Ix0HN8ey8+YtruXrsyBGNP+OiaQni/0nOhPdLwYP19wz1htEpcF8jLLEvxax0Ch8bjwWjO+zPnJ5yX5yryc1/IjY9amxvacbD+cNvTC7T9KKzuJcehhjxzUK3Q4COKhnQjV+Rbag4SGETVdmMvRbj3R8YvSCrVhUbph8ovvZssRyRNlOUh0AZahtAeJPNx0c7tKfJ6LN/dvj5OxIZOKOg45rlbIPa0AYPXOW5L1fh85vMhyWFl3tDOQLGhpjHiLMOjU0zNk4N2dSkuP9w/gplS21GEWm5Srg+8rSJcXekbnVH1LVkFjUcejdHIHPxCZ0sU5O/m2ALxEB1+2hVJGjju7nSUbe56Uyw5Slg9T3icPKd0t7Yw8LLU1VnSW3wthiT5xvL3BBmvhWm3tPTJTyRer6GhoObol5co/19gwb6zKkBvkKWxK2qbrNXYBr3QmvVzvyrRz3GRdbOHHNhAewbFQYnKQk2nWYL4Gu5x+7Sy5g/6qRuM4I2AEjIARiCGgNnXxmHOcAu3+Kfo9pcG4ka2Lk33oWDfuxcdEmXkUY/yDoOfReVbJe6wvZfz2n5JYaeOM40Vs3lKZxziz/XjiRVweQ5Sic1yPycL8lDLL6zxWf1VGjq67PCgOR0n3hxzSOdPzDDzOFtBFmJB+xqrUJceV+C0Zd74loTnPpfzBdcJlk3lHa74GiuimhT9qU6KbtWclzarXpjJJuA3aLKWV7aovu5xHPTta5zMd5xbakq+8LvJzzAbeUVqhNixKV+oQuZbccFsn2lBZlumKZxL3lrRLmUuuVR4DW8kypAPt9qRP5XEizEukM18+U8kJR7cs6fkTHTXPdMcv5XAq0fBwUCCTIRnW9yLIDXFHK34KJ0+EkRMN6NVCH5W7GV0Um5QgGFQLu0Eh8KiVV4OIRaQXb2fK1ZJVZYBx8XaGkittDbvDcdUSrtXWui3nwmu8sq6Wd96Os7L01eghb3h5k3bQYY/ouO3qI+VTkVMhXx21ti5G+VtwzPZJ2/9iJodbYD6ThB8bASNgBIzADSEQ7bemshzhZz70VH05K9fKwByGSTHxjAv6FfqKY6L2ruL6lXaK6+ZAxfhhdp5VJiZ+xl7VvjQ9q441SE/Pcdjl+R9io/14BB/kEZry88Ay+ZvHWOR5KmQaxrBsE2RVIXFgxLh/vHVQ0dWw51i/Bb8lY6BrHY9HcYnaVNSeMYBo2lVjUST1/mDBkmwvUuc7mS20U0oU8VP5ySS5Tsy2YWLAnzNLp7rVt3U5kZlztBwR01KWzDHmcG9JeyYb3eNubhchjNI8SoR0JCg7DnjcaxPSMd23KpjadjmcTBhJ5zVPTOPJbRctMLNTjvvqqjKlwYqu3Ph2fIGfbIQB0vOQKF8Zk3weKJKx0ZlGFy9sh+mAqO2GDrkVx7YU6tQXb2dSmxVkpS3ytoTBFo3vsXDxdobya2xNvNTz0kZx5DEAPeZkuVZbI5/dGzdwK0JuJ8s6iL0MHJrJXrCZfiBeyKhdZlzB6xRhbV2M8rfgSL4zvnOd3RaYnwJnp2EEjIARMAKXgUC035rSdpZf4yEm6XwWZ3BIIH1hftZPZNNY4ano+7iU+DOdy35wdp6V+PJpsi9VWujCaqLOOZgZivMbui7HONF+fBafIo3W/BSsw0vlhzEoOh6sVlQe83ge3Qjcs3OE8uEPKThHnXbw7znWb8FvyRjoWsfjUVyiNhW1Z8o7lLbs7JkO/uCButOFdE05Deq24pkbROo8jqYw7UOqs7+sOK3Ve/w9r1QXaKNyO0Xd6A+e68jPqDv5uqeBfkyn+9YQLUfkhssyiGVL2nP5or6h36ahc9wJaBqtl8oUHUUXdI3xsaXu04eYh9U6iv9HR+fhzPE6f6M4jKEPuv8i3ZQTV9Ip7zsS0eIcKTsoDP3oyryOMfZDg/40QPp2ogHoc4QINuBCftaGziG7Vkgr/5XY2ZJVZUBxLXaGrqttTXWWzoT2Als6Fq7V1uhoB22Q8kx+aRfZqt81xorjnk6wbzdSHJ09dNH6+pfoCcjbPayti1F+ZSSEY5HhnP8O3yK+v9wQ816mL+4PAexIR+uf7dwfUM6xEbgRBKL9VmobDuY6Uf4JuOjbcv/WkSgdxg6MFWiLunlUPiuuH2d0xPF5ViLv05rqS5mLHfzBmNJnG9nPymv5cj/UjzfiE5035vzMnf8rgo+lP0X7XV8AAB4TSURBVGPTMjC/5A/GsiMSPLrvkIm22xao86BcSubK9eqxPunpWGxfSd8l486rHI832FXUpkL2TNk3pE15vBzZy8E8QGUXrvMttEW6b6frfk5SPCM/rKblG3bUiy4kexr4e/Kz0Zl6EqkrEbpjekbLEfVCZdmAZUvaI3iqt2U7WiVojXxcMFChacxwcjGJ5PxfFXI/8dQ1b2hQYrC6RPEsrYY3O+8wGAz4P/Do3AVdf57oMBi+oUCggL/KdHrO/RarypBN+F4HlacalF5+hlOMQEdGHum4cHCcJESwkSIHKxEXKkde1263XZh05+i5WDurlHlkVRlYXIWdoehaW1P9oOOhDvNh475+I7sSrtXWvlU++WOd3KaRNfLNv6DlASBxud34UrT5ZQWrj3kD1bedEM6EORxn2Bc9Xtzmp9Qi/FEccwbADCx+zRGV81aYV0Q76o4QWPpZhDuCyFk1AjeHQKTfqs51EhKz/CViaQzB2IHAmII5x68aH7CDiJeDPMNZNg6D8YPoQ/OsQsjRvpRxjHRhYQZbRl8WfOV3jLto0bb04yF8FuQnrz56LqW61YQJS1b+dDs/dE/ajMXyeApsu5VBXUb0k/LN7WABiniY932j55Pf6oZJYfVYX2msta+lY6CrHI8/wB6aO4bqiPBvsWeSn7Vp7EY29LZon+v8RtKZMUY5XyA6XOdbaJVm2JchnfgjmrHPZuDvSfp3J9EyD5pqw3rSCF1ET+kXKkcSFm20LEO4t6TdZ3z6Ajso29ZpyoYnr0vJBvL9SVWorGp5Lr2oKAchPacy1LbmHtATIR6chK0T6aqsc0amfHQrwpT/3DH1KkWwEQ2VD6dkGL8+gYmLlO7khxQn2M4anXSetDOUEw1v7mgMIw6qm7GzlHfqzIGtJfthINQ5ftM9jefBW4X07KZtTXmkQ3ui/FfbK7CMBMmh3cPWxo7BCPtd0WyF+V2BdkGZVfk1taul6uLlhQGBl4v0YYwFDtoeCOaCZDH46eub7ukbafdo7weT5jlZezxPOLHqpmtryzTW4rCWv9TlnNdbYLQWi7X858TvWtMW5vSXVzXmvFasb01v2Q4TeBx0rDxi3M7kmsUm9Ev0L7POO/GcdU6p9JvHneLZfO4nrByMgBGoIJDalb/1iNW+rAoMB/Ee7d8ehyWdjnCrVWWlxjTGDH4PBsAl0RVcdysRVais7Jl7KzSVHQwoT36maFrjcSIeOBJbhZyY/qidpU4OnEJOu6T7rdgZ2TmwtYQJAwbeaDLIIVCnphqle7A13mB2A0DAWBGy8yG/1Voh6uZZt8L85oG6lAyqvWBy9J0O3j4+0dFs55LBCglW55cTrt8Vj/Mt1x+RhANt12D3QJjzNIQ48g8ciGtxWMt/mqyHU1mF0Vos1vKHc2nCMQLXOOYc58H3J0ZA9ZWxAy+a+zGbrrMtscoHBx5biOfmV+ce6y8ZA+0xHgdOByNgBA4RyGNcXjK3htwmVfkeVWPPG0mD1P1V91ZqqGFmy2v+i+OtxJ5DDktCcdoNBvOKI290JM91dNeJbqAjdDxPeAyerbmRPD5QebJtxWt0LXgn7SzhxKoyJoQsaQ/ZTsIgRFvocamXNVtj4gxunPPxGRiNM3EPtqY8vqF8U6d+Hee/9V4Y5jq9auVea7rXRr8l5teW92vWlzZCB6vecZZ935oXlTtbyfgXxPGEi3teJjQHyeIjxGXbhW5zf7bTnM4SBuU3b7sfsK/FYS3/QJkFN0qfrYK1bYHN0iRnFUZrsVjL35xhM/QIqN5e45iz198XZ0OAFzWM0fOL57Ei9AGz/Ql9R5KTJ+djObvdS/fmcad4dpn77ZZJCzYC149AbhuaXyrP9W+PLxCbg5U+G+nYfYhVslhpdZVBhUmnchAUj2FMrXoq6emQqjJKoju5rtpZ6uDAKbqqbAzX1dsZGarZmuLeHGf2yP092Borhwjjb1g8xLb/Uo+zzHbu++DI+GyF+X2gdv25pF3Nzu0yNzjNv2Ayo/apdMKVNLPX4mci90zH2R3nSRfyUsvPWhzW8s9iOUPApJNjVdgIo7VYrOVfhYGZjYARaEOAPkJtB2083yH7ROe8GuZtXdMufSqaWj+jRweB+s+K31PPKZeMge5hPH5QQI4wAmdEIDvuou1JWNVHYcrTEdZW+vAdgm4lmdR4rqO7Vlz1jWtN1TSo588xwjw1Odcal/LN93Kavb/XmucZvat2Jp7wqrKa/Hu3MzC5I1vj7W23KrNmCwviaOBzY7+A/S5Ytsb8LkC7gUyy0vdlJR+5P+P5oqD2ijrHivWWzyIsSivI9In6kakV7GtxWMsfzMLuZFtgtBaLtfy7g+QEjIARGCLAGF0HO2o4+P5U9w0qnZkfhSfZouXFyjnmlE1joDsajw8L2ndG4LwIPCV5tRN5jLqZNo83k7SRIBrPmqiU+ciqshp7F5dkzH27YJL/mh8o73eZ76kyO2JnLavKquLv2c4A5F5sTflk4Lblyi8+mvxMA61Vq4cog1sNO2B+q1DdTL6oD4HM8G2igyBenCusjKCu4qBjstYPpPQ8x3WrJtL9LoMtpT0blD4vFlkdcRD0bDEOCFvDL96jOB4ou2PEFhitwWItljtCY9FGwAicEIHUl5x0btU6BhL9SfU7IfxOyghcMgKMmbacH/Z5vTjHXa+ZL4yAETAC94UAb1IJNPg/dVf+MQJGIDvlcL5NhQOnlpwzrKLDCd69DEzOJ1ZUdy9ndI/TDifZ0s8iiHW7kPRhNUjvWBxJX4RDIWMR/xyOhfzdLzfEaBEWRQbX8heifGkEjIARMAJGwAjcAgJpnMKYlMUYmwc77jaH1AKNgBEwAu0IaML+Qg0+zgm+vWLHXTuE5rhfBPhGUR9Uj/he3TPVqXf7yId/kH1Dzz5QPG9CceIxuOLcBz2rrvrvCfa7YKvWql0FUm2AwwJVB/xBHBcks5jllBgNsFig8Vr+BUmaxQgYASNgBIyAETgjAiy+IOwyj7Pj7gFc/xoBI2AELgGBH6TEx5egiHUwAheCQO3bdlm1vPLprxyRzt/pPHaCscKO0K3Ok5Ns0WcR5MyC/39ZDgIDgX/Unfx+kmTyT6vVLbKF7CU4FOzVbwTm54txzALKs/JDXvLgtXzUpaPnNefoC2HEtuZq2Bijc2BZzZcjjYARMAJGwAgYgZtBgE+v/KnxzNTuiVUZteNuFXxmNgJGwAhsigD/UvaZJql5VdCmwi3MCFwbAhr88E+AqN053Eb657h+gETdSbQ4wcuQHUmTDrSSeOoaffTs/annrfHSF4ciW3r7PNRktOIwltHKvwZHpVVzzPGdPVZC8o/uTd9d2hqjVizWYjnm970RMAJGwAgYASNwkwgwzhm/ON4so3bcbQalBRkBI2AE1iGgCeUvmqQygWflyS4fNl2nobmNwFkQoC7kFXOlAnmlWFlXqDus3sLBVga2oO/2FrRMqPGafD1VvcdpX4b3dPNOikdvBoItOJSy8nUL/yXhuAdGLVhk/MrzWv5Slq+NgBEwAkbACBiBK0ZA4zV2TxC+fTht/2vH3faYWqIRMAJGYA0CfFT/G3UA/APm2PmwRq55jcC1IoBTi3oxDqx8GzvpWFmHU6UPqks4wTg2WynXC195gbNeIgb6IlI6/028npfbR1twQMw4tPBfDI47YdSCxRhH7tfy12Q6zggYASNgBIyAEbhOBHjB+tOec7dH14mLtTYCRsAI3CYCavB5U4PD7vlt5tC5unME8kf782q5Hg45q/jziH90jP8wgjrxUvFsQegCtLrge5CfPsR0zi7iWJ3Vy050OFn4Y4NV22RzOic6kxeOPqS2IYSD8n0POC7G6FRY9oXnCyNgBIyAETACRuAmEdCYi5edjD932yYLcF5xBwoORsAIGIHLQoCGn1V3X2mC+eqyVLM2RqAdAdkyzjMCgxvCj4pjW/jPyYnyGrae4n7rKIY/rJb7Ws+f6syfUXD+r3hKZ1yW/aXovtBzAv8se/TPITqqC/mR3vyxA4M/wge6B7dflc/8XbhZHG4dxy0w6tB9WIF51KY2wDIl5ZMRMAJGwAgYASNwowjk1XZHv1e8Nu+va1CyVob5jYARMAJGYGMENDll1RFb5XZ9e7Ox2hZnBM6GQHLoPFGdubgtsWcDZUHCe+EouYv+nGJBFsxiBIyAETACRsAIGIHdEdDYhpfGP+t4k5d9eyborbJ7omvZRsAIGIHlCLAFkH+YzatvlksypxG4DwQYPB18L+4+sr5pLvfCkQHtroPaTVGwMCNgBIyAETACRsAIHEeAnRLs7Nh9fPP4uB5+agSMgBEwAudAQB3ACzntvlLadAgfnkMHp2kErgUB1ZU3pCtO7l+vRedL1HNPHNWm2al6iYVunYyAETACRsAIGIFmBDRm4o/T+JO0n5qZFzB4xd0C0MxiBIyAETgFAuoI+K4V3/2q/aPmKVRwGkbgWhB4khS1c2hdiRnHdfiZ2wgYASNgBIyAEbhxBDQ3Y3fCB5qrfXSqrNpxdyqknY4RMAJGYAECqUPgI/X9P2ouEGMWI3DrCPCHFh+qvuy+VeHGgTSON17Azp4RMAJGwAgYASOwHAHNyd4TN4sq/rtcSjun/5yiHTNzGAEjYAROjoA6CT58+qUcE+W/aJ5cDydoBIyAETACRsAIGAEjYASMgBG4NwQ0H+PTLP/T8d9Tvyy24+7erM35NQJGwAgYASNgBIyAETACRsAIGAEjYASMgBG4CgTu+s8pksf0eSqp/M+Nn57ae3oVlrKTki6DnYDdUKzLaEMwLcoIGAEjYASMgBEwAkbACBgBI2AEjEADAnftuBNOX8tJ93nGSw4K/r3xdx3v5jifd0fAZbA7xKsTcBmthtACjIARMAJGwAgYASNgBIyAETACRsAItCNw11tl5aj7R5DxMevuX+h0z6q7P3S8rzh/R6rdnpo5XAbNkJ2cwWV0csidoBEwAkbACBgBI2AEjIARMAJGwAgYgQ6Be/9XWVbb8Q9qDudDwGVwPuyjKbuMokiZzggYASNgBIyAETACRsAIGAEjYASMwIYInHXFnVby8K8cH+hgpdsv517lJn34W99n0uNutsq6DFTiNxb2LtN7rCc3ZiLOjhEwAkbACBgBI2AEjIARMAJGwAhcCQKbrbjTZP69cZ5xIOjAKXcQFI/Djj+GYEvqtzqeKO5Hnc8SlDb6P9Px/lkU2CDRlIeBJMW5DAaIXNfNpZXpLdSTa7CA1nK/hjxZRyNgBIyAETACRsAIGAEjYASMgBFoR2Cw4k6TRVacEf7SwaozPkr/JxFzQbx/i4YVdPnbcFwT+F7cq4fLh1/R8uy54r8cxX+h+6eK/6iM3/ta+uBc5I8pPhrrunfaY/kug+suA8ozWoYROtFcTL2SLjdRT7Yuo6g84de1e9CnwP2PanO6b2zmSM7Rck9lktvRJ2J9qeNLycztMOK6INqm9r2gzyK+n5FbpWvREQGipy3O4S1dzP7Tt3h48fK59Ov/bCgL0LMW3GfTTvnJ6SCbekFfWStHns2WT5I5S5fz5LMRMAJGwAgYASNgBIyAETACd4SAJhuvcSjwb6psE833TEj4o4Z3ctyxc6LFyYAw+JgkvlHjSc86h0T5HHod8L9Xxu95rbQ6Z0ROI92H8px5tjorbZfBgy1SJtdaBqEyjJa16KhLZ69XlIeOb7Ktp3uXUUPbWeJX4Mgq477dLeJnyz2Vwc+Zh7MC7S4XH4ziQ3aZZFDW0PcydI2ef4xkztKJB5qojrnP+SKnI14cctj/UVvTc/DCCdr1X+VZ8b3d5nhodfS46zqUdqIbyEOODjDv5ZGOQijvUbqsu8+HZWxMjIltwDZgG7AN2AZsA7YB28At20B20n2mycNgUkamFZgEDiZdU2CIrjppqtGLlvQGk59Mp3gmaf3ELcfvcVY63cRKZyaH+WDFRdXhOKeD+Nj++9kcXe05fDpcBg/lcJVlEC3DKB12onD2eiUdbqaeRLHfie6gbVA6OIx+H7cJipstd2h0HDi0FEcb+neWqeumtkX0OO0GbbDuf+bIMpNtztKJJ6Rjkgdtr3dOq5Z2fpb4WKldddwpnrzP4i6aUNqiI62Bk073+YXToByTzEj5hDEq8+1rD05tA7YB24BtwDZgG7AN2AZsA/dhA480uSCwNfVga5XiftXxQdpqpMttgozrWx15q9FYKFu9no4jd7rvVpVINud8MMl7tTA9JnAcS4LL4PrLIFqGUbomO9qxXt1SPYlivzUdnx74sKlA54l5UfBHpX1my2b5bctoXtimyuoxXmLw3dE+yLY+5MgRUTrRR3VENGnXPs1A31Tth6QHutJeT7XZUdyjaaPLIL0j/UU071E6Je1gBIyAETACRsAIGAEjYASMwL0hkB13TBxwmI1DnkTxfNOgCddnOr7R8XU6mDgRWN2z1PnVCYj+aML1pg6+8zc4ovwb07kMinLYGNuouLVlEOWP0kX17un2qFc3Vk+i2G9Nx0uQZyqfn3WU7RurmsvvqvVlGbjAQfenymfKaZXTieaFJHmh8uqITGgIUbqQjgUmtX7or4ckX3uSzuXpE+k6cDKWD3U9i3tL2kqLfz+n3+i/Zyf+3HeNyzGUd+kYpRtlzbdGwAgYASNgBIyAETACRsAI3AMCj4tJy7H8vnXsYX4mWWxLYrL4tg4ccF9pgsMKhUEQHduumHD2q+4Ux4oK+AnZYfhwl371nElunsQxWf2tJn/AdAU3yleeYB/T1mVwDJ2Vz9aWQZQ/SldmRzy71ivSUho4d1iZRb2i7vIHB9V6qGdnCUuwKxWN8m9Nhw7C8ifJ/UmXOHn+1jV/RMBqMLafEn8QRHO03MVHedUCq9BI80U0L4UQnGN/ig8Zn+igvUVPtu72zirdh+giOkoWur5SmlzW2jn6EwJ22QfRs2117Czrn3MhubO4L0k7JyIdqDf0S/wxxsCBqPvZ8kFOlA5aByNgBIyAETACRsAIGAEjYATuD4HHynKeKDFhnwoRxxI0PzAJQogmNEyy2MbFFqt+wqd7nHZMVnqnXbr/Rc/YkkfI5+4myYIPZ0I3yVUcEzYmRv32rY74On9cBucvt7VlEOWP0mVEdqtXOQHVJRwPfNexq5O6xxlBHXwz01zIuRW7sdpR/q3pOj2E70ep3cIhB+Y4RjnXQqjcx4yS/57iOsdrehbNSxZFuoQn0jf/yyntOc5G/t01OxmjdA/Sit+KjvkpsrG9cSBPhJxm7l9YGTjrXA7iHk4bRVIe0PWpDl5O/aZjNhzJ+4A3Sjdg8o0RMAJGwAgYASNgBIyAETACN4nAo2Cu8oqHSXJNjnDQvcoEaUKFw65fEaHJCBNWJjv9hDDTp3OemP0wisdp95Nk5kkjjz/W0TkBJZfvOf2oo/uaPw8JKZ4tuQNH4MPTq/u9ujJI+Oet0JQPRy7jqysAKTxbBjOZivL3dDvXK+oIK8D4N8zSkY4TgjpVc6LMZPHsj3vsFmoS5W+mS1jTRrKCjbYxv9ygDAYhUu4Dhn9v+KMD2sr/+zdq9qrLS1E3+VfvweoxSaBN/g6aKN2RVKd0/BQeye/tTtc47XK/UjrpDla4TaUnGeA7h3tL2rx4egHGOnh59L2O31M6U2rk+Km85+f5HKXL9D4bASNgBIyAETACRsAIGAEjcKMIPFa+at8UytnNqzX+yhGNZyZabIHln/W4xmHHKomD7bNJLhNZJp15osYkDscf8V8lmu4kmn41EPSiY/IEXRd0z4QvTwDfSNGbnJJOWXYps8NLz0snSH7ORG9q69RNloEy/rXy3GORcMOJiuNiVbjAMoiWYZTuGD6r61Uh/Dtdjx3puR6tqjd3UEbhshQWvLR4v6gPHyoOhxIOGhxifDutb/cUVwvjch/QJLz5BEHZzoR1LISRzjhQb8nDEx15ddkcXb/SOgub0LF7TP71/D+6wdlP+41D8dd0gFWXnp6hR/9CSNeTIdHO4h5Nu5aQeNmOS9nxYoLv31XLUc/QeVw+ByKjdAeMjjACRsAIGAEjYASMgBEwAkbgJhF4nCYsZK42Sc9xtQlaD4gmGqx8e0uy3u8jhxd5lUbnmBs+eriTjOwIG0/IWFkXmdSyZbafKEqXF7rnG09M+DYNkt07o0rBKS2clC2rXVi9wYQVURnvUmyOu7oyUCZY7Vh+G4ttgcSxoofyWRwurQyiZRilAxjhtGe9yiubsK/xCtdcF11Ghxba18eWspQYbB+nVB/Ej8MHJ/YfOsC8W1EcKfdeSLoQD84s2uDBpwNadCxoX43lF/e0b3zWgKijdAVPdzmlY0mHDroftK/iAzsC396jD2Fb99H2sKN++AnjPpc24pQ+DsVuxR3nIuDMpAz7ciyewVctn5KG6yjdmM/3RsAIGAEjYASMgBEwAkbACNwuAo9T1vK2rXFOuxVkiuwdYmOCdM8qjNrKjo5fEyIcaHkVDysoaoFVP6xK69MSD5Nkjm5LbI2piGPCNJjwFc+u4fIWy4DyyKtz7qEMomUYpdulXhUFwcos6tzYAfOJ4lkZFHWOFCJ3v4xiN6VIlH8zutyOVXDGAYQzCoddbmvRe7bcy8yJn5cT70pWv9JOcV17m8owmhfETtHmJLNNROk6voCOWX7tjLOse3kjOWDzVGdWKpYBmndSPLb7pa67/kPXY/uewr2Ul6/7tFNE99kFyZ5cWZcZ81m0c+XTkUbpslyfjYARMAJGwAgYASNgBIyAEbgPBB6lbDIJYkI0Dqygq03sx3TfanJU2/6IM40JXjdR0ulgAsWzNGGBtp94El9MuKb44IGfCVq3EoT7Kw03VwYqP+yiLDsceUyqV63k2rF815ZBlD9Kt0u9KvCj/gwcq6pLOCo4BnWx4Dn3ZRS7KT2j/JvRpTrAqtr88mKsG+1X106mB7PlngWk8nqqNMbbnXEW5Zcp0bwgNn+aICeRz/QFfOYg6xmlo33GnuZ07PoB0fInGODRhXSNnXb5I30dH40PPaedyc8yLXEh3JXOs7m0JYuAzMHnHLrYf/vPjE8XLZnRvIfoUlo+GQEjYASMgBEwAkbACBgBI3BPCGgChIOMLLNdi4/U53smT3/rYFtjGQfx7zku8TIh/WYU90XiZ1tT5me70B/5PvEyYRmkM3rOBPHnMi7xsQWq001nJqkDnTJ9ejZIMz/b+pzS+mKpXPHfZBmAhwLlTP56e1iK0zE+yccWzlIGKZ+zZRilUz72rFfUbwqGrcy5fhKH/p/luD3Okn8zZdRQljifWK01sH+w1jGwV91Hyx06yos2cnwM2rxEd7R9z2UtWlY4833K0i5oo3v+lO9ZOvG06EifMWgjdE8as/YoGmy51k+EcBdvKG104cjYJBywZ9Ifx4fyLr4QXZmmrx9s0zgYB9uAbcA2YBuwDdgGbAO2gXuwgdfJJCGtbMAZxoqCv3Q81fGVng9WR4mOiRUrGwbbUhXP5COv+nhL16z2+FJ0yOuD6JjkfKIjrwbh2QEdkTmIJ3/jiLQJOBj61Vx6zqSV1SA5fWi6kNJjElpbEZjJNjmntJq/cZcTFz/5usUywDYoI1bKDOwh532r8wWUQbQMo3S71KuEEyuxqBfUSQLXOOAHdb57suHPDZZRtCxxXj/XUbZ95Tcge5SF0Wy5i4b2ELpaYKV0/81R0YZ0zIJETzsED4H2/KAv4MEcXYuOhTwuc9pVfCAgSD7tChjgoCOw7fhX5b3/zqhoQrinvCDjaNqiI61yRSrp08eMV9uFykfyQnQo5mAEjIARMAJGwAgYASNgBIzA/SHQO+6uOetp4tN9T02Tp1dlXvQMp8RVOO5Kva/teqoMFN85IFQunaM33bOSJ38ra9OspvJe7DzdVJkLFiaccHg8UTn0zp1TqesyOhXSTscIGAEjYASMgBEwAkbACBgBI2AErh2BR9eegaR//r4d25jOGXAaDhyH51TmxGkflEFy0uEg+kbX73HomlWR5YqjrdW85zJowZJVQ4MVQi3MK2ldRisBNLsRMAJGwAgYASNgBIyAETACRsAI3AcCt7LiDucQ2436f6VNTiNWeeGgwGHE1qm/tMKo30Kle4eNEBDetTLgu1hvjJNQGbw+jvP96RBQWVEmlA1bl9la6GAEjIARMAJGwAgYASNgBIyAETACRsAIXCACN+G4u0BcrZIRuFgE5LjDmc1H/9+U447Vbw5GwAgYASNgBIyAETACRsAIGAEjYASMwAUicCtbZS8QWqtkBC4Wgd+k2Yd22l1s+VgxI2AEjIARMAJGwAgYASNgBIyAETACHQL/D0X+v14viAG1AAAAAElFTkSuQmCC)
$\displaystyle 0.5 c_{l1}^{2} \left(1 - c_{l1}\right)^{2} + 0.5 c_{l2}^{2} \left(1 - c_{l2}\right)^{2} + 0.3 \rho \left(- 0.037 \rho - \frac{1.0 \left(1.7031322 \rho - 51.0939660000001\right)}{1.0 \rho^{2} - 40.0 \rho + 400.0} + 0.042578305 \log{\left(1.0 \rho \right)} - 0.0530922164415325\right) + 0.5 {\partial c_{l1}}^{2} + 0.5 {\partial c_{l2}}^{2} + 0.005 {\partial \rho}^{2} + 0.00085206629489322$
2 2 2 2 ⎛ 1.0⋅(1.7031322⋅ρ - 51.0939660000001)
0.5⋅cₗ₁ ⋅(1 - cₗ₁) + 0.5⋅cₗ₂ ⋅(1 - cₗ₂) + 0.3⋅ρ⋅⎜-0.037⋅ρ - ──────────────────────────────────── + 0.042578305⋅log(1.0⋅ρ
2
1.0⋅ρ - 40.0⋅ρ + 400.0
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAABPkAAABMCAYAAAAfkWLUAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAgAElEQVR4Ae2d69UdxdGFD1oKQIgIbDIAKQKLDDBEIDkDvPTP/1h2BpgIwGQARCBMBvYXAUIZ6NvPvF3jmTlzqbmd66615nRPd3V19e57Tc+cD96/f38wGQEjYASMwGUj8Le//e1HafhXub9etqbW7hoRULv6o67/XqPu1tkIGIHrRsDjz3XXn7U3AkbACBiBYwQ0tz1R6L90/Vn+d8cc+4U82k+0JRsBI2AEjMAWCGhiYIL4l1wb+LYA1DJaCKhdfaWAT1qBK24kD4PhK10sbkxXgoDq66gNUIe6/nglRbCa14sAYwbjkMkIGAEjYASMwE0goHkNw95fdf2b9dQpC/WBT/KdEm7nZQSMgBGYh4AmBQx8/5XLJGG6EgRUXxhMqLtP5Z/19E78fy/F/E3ux7r+rrBdTtlJ7ueS/1xub/tS+OxyFJmUfYjeiefDiJR/srziwdD0l5KGhRL34PJTCWs5Ch/Ve468whv4PFNGb3UdnarN8rUUvaAb6f+71AHbeJgQC9LeNiz+UYzHiqa0k3VOevF905DzVP6XCmv1J92n20ZGXiO/g/hDzwj+TmGBT4Sl+BqyJvt1ljfLh6JZ3q35AiTJHW0visfIx9jwz0hj1wgYASNgBIzAtSOgeY21NmtW1vQnIRv5TgKzMzECRsAIzEdAk8ErpfqL3E/np3aKUyOgesIo8q0ujEAYg9jUfqjwllFCYYMk3n8r8mu5P8AkF5mEfSb/poa+Ivtnua32VcIXl0PpMYxQ9j59XygcA1mUb7K8RR8WR2HkO8jPgglDIq9AhKwU/uKHb1KeeMgHA9I3cj/jHpKf8mGQoE4qI6PcFB/pL5VUhv9INwxp4EPdgStt8Z3ciuRPYRz8fa5kZOscPrD/B3Lk0qZ+1oXRsWpbclN1Wfgm5ZEPJH7qk/ZFW4065v4T3deL9Bl8k2WW7IokM8Wb5UNolncHvlntpeT/J7l1m3tAxb9GwAgYASNgBK4XAc1rrB05tf7nU5Ti0SkycR5GwAgYASMwDwFNAmwymRBOMhnM087cfQiwMWXy1oUx6rs+nrEwpcOoy+uRldEKXvnZ7HLfPNFE1BZE+zqSS566FpcDxZQeAxgG6vpSMPn9pLAwymXLCx+v/2LYC4oTfK8jQPFZvVPyilx0ro2LhCkfTvVRL83Tilk+RFwq/aqyYZT+QNfHlFNXy9jCva7FbUNps3WOkfmp+CsDH4DJzwm6X3Q122y2LrPyyAqibjm1F+2MsCe6uobrST7JyJaZMqZ4s3woneXdmq/kPbe9ULf0JZMRMAJGwAgYgZtBQHMsa0ceFDLP70428u0OsTMwAkbACCxC4Eel4hRNd1O5SJgTXQUCGHSPXgVU2BtdL9QWMDJsSV9I5h6vxqFvH7VOz4khW14wwdhUG5ykd+3vy2gibI48Th7+pwd7jD9PFI4xHsryPXDf72+2zjHo9o191F2zL2TrMivvoDqFl1ODrb6hcAzXzROdKT7JyZZZrGnec8qckzdlSpPwBXPGpa3HurQOZjQCRsAIGAEjsBMCGPpYC+8+xz3eqQAWawSMgBEwAgsR0ODPU57WKZaFopzsuhDAUNQyLBT1w9hBfH3Kb03R1MaGjB5rxFZpJftIR4VxOufrjvBUeZUWg1r9DT9kKAz9oeaproeQid+Z8sibJ69DRsVYqGX5JrS7+ejJOi/1AxBve9D4rYQ9k8up0Mm2IZ6oo0l5RTYnNzmBNlTnha064ZnhmyyzBEafyfJm+dA1y7s1X+A012W8+0JX31g4V5b5jYARMAJGwAhcBAJaV/ygizdQjt7+2FrBx1sLtDwjYASMgBFYjQCDvzc4q2G8HgENQ8SY0nwvbSviRBIGkt1JZeO0G4ay+POKw5ryKi3GiGqBJP/qfjImT3GcWuojTnodFM9JMtwUH7xdUlrKQ3qMSmDFa7Jh2NXt6Uj5Vq+LKsePdKELp4mrMq7VQnLC2DYmiocbGM7g6Wvv6AWh2xEp3VHbWCAPAyJ/dkQdf6kLwyLf4eMfzpt9ZpJP/KkyS/Yhy5vl20PmnLzJfyGBMePT6r69MH8nMwJGwAgYASOwFwKshX/UfMqJvt3Weo/30t5yjYARMAJGYD4CGvD5oD8bw+6pp/nCnOKaEAiDBoaeIcoYDIbSdsMxUMw+BdcVkrzHIMfVpNnlVd/A6IIR57kuDE98n20xLZVX0lXGuLHMM3ziARde+62++yeX8vHHC62Ti2P5bBhH+/peOlRtUC5l5FXl+g9GVuY1p8452QYWXaINQK2+IB2n2sYceSH7meQ2DdO/6/6lrjh1l+GbU+Ysb5YPnLK8W/OR91LiD2D66n6pPKczAkbACBgBI3ARCGgNwVsIrGFZ/w09IF6t66PVEizACBgBI2AEtkSAzT7HuauN9paCLevqEYhTTFsUBAPF2y0EjclQO8ZQ9IJFzRjfQFyrvCyKdP1DF4ui73T9W/7PB9JOBq+Qx58t0EfrP4UYyGyUr+j+udzKwFdkYLjE6HdyI4fyxJhXjzvy84SZejuVMRgIos5fciMdahzk/0RBoV/r6bfiptpGSp7k0C8gTp52T5J9r/Bv4cnyVZKmf6LM05z/w2eK95wy5+TdVw7GJcYNkxEwAkbACBiBW0SAgxys/2LNsXkZH28u0QKvEoHSyHhHHIrFFU+sY0H9EONfI7ASAbe1YQCFDZtY+t9uT3aGc3fMmREYM7jFKZv4HtmRqqVf/a4IjB2fHjEcByDzFOM7BqyWQaaosra8YQjnFUr+EXZVWZQ+JU98GLx4lXO0jyb5+MfX+qRYwSXm38GFn2QT97OuQZ4iq+nwj7i/NgOSfuoOI+0fdfXVY1JMxZauc+XFK7t/UCpeZ2FcxHDEH7pwYdgd1EX8R3W5QF6ffE5Y8jozp2AxxkJZvgfu9m+3X2fxyfKRW5Z3a752Sefdgemctj1PurmNgBEwAkbACJwXgXjwzZpi6oHxIk1t5FsE200mav3roRbEbGRY0PIdGpMR2BIBt7VhNDG0s7n9dZjFMbeIQDFCULS+zW2E9RkUKjhKeuLD+FCFX8BPr0FmTnnFi5HnILfbLygrJ724ftCVoqXylC7+EIfvhQ1Shk886Ey9cjqsSYRD3bI+hOpXad/JyRhy6zRTHsnk37z5Ht6Q3GiDU6IG49FbF/F9siKsbuPwi7d5yvGgsHjtu+LTfbptJOWFjuQ9RBg8ed2G+CxflK8pM8KiLJF3hPfyUo6S9ygfibO8W/M1FV/gx/j5bkE6JzECRsAIGAEjcPEIlDkXQx9rnPsw8qnQLFp8ouz0zfOVsG9+VJqFNGG8sjK42Ti9mtvk6Ha2DY4LpdxVW5uJERv8eLozM6nZbwAB6j5OcjWLEyd+RtuGxrU5D2XeKoM+I0Ez31X+Ms5SnqE5JFteHjgdJG/1ib1SoNnylDfGyo/l1if45K/qSm7TMJXikyzkcOqya8zgjx44KVjL1P0piNNptIkuVW1P+gzVYZd/6j5b50NyMOphYAvcZtdlR3BXHtFDOkbSqJu1fH39ekhmlzfLN1aepTLn5B2YzXEZl/ra4hwZ5jUCRsAIGAEjcMkI8EmXb7Se2eJNiaNyPj4KOX+AT/mcpw6wJF/aCZA9kXA72xPdcdn31tbG0SixGuTZbLK54USN6T4RYMKPk0pNBDhd1WcQavLM9WOo6DMozpUzxo/hCBrasGfLi0GnadiphOon5GN0mEOz5JW++Vxu99VaDHr1d9uyfEXRI4N+Sc84MHSabk4Z5/L+s6d8yDjSc67gDn+qzqUL2PI68x/kp74Ochkf0aeJT6ouZ8gjK95kQM8ukS+n6KK9ZflSZS6ZZXmzfIjN8m7NV4o028H4GIbU2YmdwAgYASNgBIzAFSDAmxysI1jX1GvJrfR+tJWgDeVwyofCBrHhwcLJwnc1Sc7em5rVOu4hYKrcimeBXy2kS/4YYjhNsNXT+z2KtUam29ka9EbSuq2NgDMeFeNebCDHuR176Qh8VBRkw9oi9ZEnut7rqk4hRaTumeTfysXAUZH8GDa+0FX9ccBD6Ca/jO3PE5Jml6MhE92h5tzyEKLfGeXFuNYyfheMkP8X+fvkD+qtNGl5ks2aAeMHdcYT1/pq5p3lU5qDeNEbuXXbKGHkQ3nOMe9W5VL+NUkP/ukbqk8vPtzWv4MYUx5da9o4+HSNw334ZOsyK+8gvXn1G6NybXCXnzpr9cMZfOl+LZkp3iwfNZXl3ZqPvDs02F46fJ/q/hx9oKOGb43AdgiofzEGXSxdun4XC5wVMwILEVCfY+3KNfoJmIXiD4+XJuymk6JHr3UqjEUR33iZ80Rut1M+0oMFK7rM0adb1Gu9x1DKv7hMvvctHgyqbDJZaF0UoZuu1uJP925nF1VLlVH+6tvaGSBlkOeUyN2NTyozi1+MnN/Lz4R3tST9MURAYbTlMwjU6Y9yqyd1lLGE9Z2eZtzlpDEGuN904f5J961xT2Fr6TsJCF2PZCm/iFtaDmSiM/X5hpsBmiyvdOEh1AtdPPEMos3wb7Ato7juM3qn5SkPjIvkxff4utSskywfMgLTv0rfMKTxqvXSP8fo6jX7XnrwUI92FxhjgHyrqz5JF0LFk8F4VRtXHvyL8kfK87Vc5ngI/br1narLrLyHbCrDGG2ri8dRPxRPik9yJ9t55D2D95wy03kLo8n20ig7XvoHewHTQgSEeWqtLL6rmHuz5VkI1+7JpP817D/T+8TdATtRBtferk4E08Vnc8563CBv9gHYXTanD96/f18LlaLx1JKNBQtOFjipDaf4+Fc/FmKx6I1F2aeKY5G/iIpOGAzmfGvoKC+lx2jV97pNxat4AGYhskrfStiF/qiMTDIsvKuNZp+aimPCZ5HPZmNxvfXJjjDJdTu74XZGPauOb6Gt0Rdio8F4xv3RJlNhm5AwYwz9Re4uT3T6lFRei/si8jLpxTOJo3gYn2Mj2Kcq49aHfREOW46AMP2PUp/NsLRc8+tOKdyZY5/JxVBiMgJGoCCgPsF8wQORVWv+LQGVLrvPk119lSd7Ek71xhqkYtE9+MRnA57JjyGehwWx9wq+1J5M6dJzr3hZB72uMnj44b75Le8qtOgYek+uncQ/ia94UuVp6HZQGjDklX8w42Tu0InkZrLN/cp3t/2nZE9iN6dAkje5dp8jbw3v2rJl0osn1a7EN7mOjbJm8m3wxkM1gp7qeqn07yIet+Sd6fOpsaEjO2X7aJQpkn+nsNaYExHhKr53DIv4cKf4SvyofUY8m9aj5KWxzOYd5e26Sk8f/kruB924tfePQ4CE89rQ13J5TeEgl4H533J5Spkx9DHRQFQq/MhBXquxKixNSossBsdVC2HJoSw8DW7JKeFMAOjOZEml3iypvDwZp057T8sovGrUcisjQ7k/yM3Ufwo3yXI7u/F2RkNQPV91W5P+jBksnGOhSpkYi9h8YBSpxknKugVJHvlxbdbXpvRSnqv6YiZ9KVcGR06rcUKnr/yc6ogFzlSxHD8PARYXtPG6nc9Lbu6FCNCmNx1DFurhZEbg0hBgrA/Dxdl1y8xzY0quSM+mtrWJlizW6LxWXz8IlB+sYq/WPOWa3ZPNmXt5yNmaK3SPke+JrubeMTPnV7ApXXYdki1PJZcfyQY/Dm6Qx1lIee+2/yzlWrNvP8JEMkfX7kcJdgpYW7YZ6SfblWSl9wPZfItM2iX9uXrDTi42j/+TS5ut1sJyU30+y0d1iZfypGwfRS5jEX26Gl/kcs819SDmaAxTmj464pujYxG4WT2WMmfHWbKfzLuv0I0wHrYflO8LXc0xvMGyzFsZ+SSU11DqQRpRCuPkBIM2VuZ6QiFugPgo+GZPSSSLhs3ktcXJOuQ0reVVESijPJXO8vP0gg526wQO4NGdqKuBROF05MABnsnNtfjZsHDMe/CEoOIPinc7u592RpVfc1ujrb5Sm8WoF5vxGHxfKy7CKOcWRP+DqsH+wbvf79q+OCN9GkfJPJpnFFbhIndrvPcD94okC1dedaSNM373GVivqDTXoapwZoFNu35zHRpbSyNwGgQYh5QTY1Frfbomd8lKrU/78lDaVWvWpemVjv1IH/Wt3Vmzoycb5eZp9/SeTOkn596SR5+x7KXy/VlXzNFz5vw5+KbLI126FJvwbvgp7nfZf5b6WLtvHyp/79p9iHnr8LVlm5k+065SbXpmvhjZ+JxZ/Qkt+dHlF4WDf/TJbJ/P8h2UxxzbB+MKp/Zi/6PbA2uY0fWi+IfGMNLXNMQ3U0fkbVaPkpXGshQkk3dh7XWocwjbSxPnKnDsRziNzm+PSmIMXb/2CGIRimWRCj0ZKT8meiYuThFibGTSJ2wpfaH0owaopYKvLV3BATy6dcrkTWPBjQsjB4PBFCGrK68vjdtZHyo3GnblbY3xkLZft/9kX1hamzG+jU6cS4X3pFvbF7PpszgOGTyOTg70lMVB6xCgLllUmk6DAG8NQLMWcw9J/GsEbhoBxqHNDHwFqez6tA/Y7DzXl5aw2em1zmCj11p7NISzRv+PeChTkxhLMPjEOqIZN+XPzr2c2gnDw5jM7JyPjNn4jGV8oXF77T93w07tiP1y3z7xVBCvLdva9N1yZtv0nHx5M6hvvU9eTbtLts9n+bplG7xXO0BHxqOW/UTh2GYGxwLFjY1hdX5ZvjrBek+2HjfHckL1aAdTJyP7xIzOb49LCgrUqsQSHhkTH09pStQ+jiqdSYqJfvaJsj6NJG+oI/Wx30sY9fqFrrrOhVPzCeBeOLid7YXs5cq9yram/sCiudUnylgC0nsYRJ6WKowxt9zu5qzti6n0WRzFdzS/KIynaV/vhoAFVwgIZx6k8e2nr3TVT5UNz24I8NS2eoC5Ww4WbASuDAHGH6nMOHSqOTCDUGqeGxG0JP2XwoD9T5+xk3UJf6iBEbCPusa/Pp5WmGRl516MgcwRP8ptfrObebpeEyl+ztppCT4t/c95o7KOnoBX/J77z72xO1q7nxDrtWVbm75V1BltOpWv5EU/fdvK6OHmtxL2TC59Kdvns3xFfMphDGJ9ODTeDAkZG8OaabJ8zTSL/TPqcQ8sB/UGX13Exz5wkHduxGMJjsY2ljaVsWRVR68l6CNdGOv4VgCW0yNSOBNDNGYaEB+ch5dTZOiEW5Pi+ia8On7Eg7WZCrsaGsFmqzKAB7jURr6tBA/JUZnczobAOWO421oOfOHE5F0d4ZZ/sN8UPp7mMaYxBrJYz2xaon+QbleSPpHXWD6DY/6a9AWfDI5gx2Ym+7mAJZiPlf+u4oQzbdQGvhPUurCmj1/VmuQEsDiLO0dA/eKixh/pc/J5Unli6KwNZt0moXjmuT76hEDFt/Zbuk/vyUKo0vTOvQr/gUt8GK5+l5+5mZMnzU+a6LZN4utdOyl8Nr5KM7s8bW3ad0UH1iPxmRTK0/dHIpSBPSjzJPhg6ER//pGc8g+dwNll/1n0Vvaj1LuGU1rKklkvnXyfSGnWlG1peuU5q10VDFvr2Dl6i3fMqIP9BKKdHcSb6vNZPmTOIAyNtHHGly91YbPp7SMKr0i8o2PYXL7gz7jKe3U9ko/kpDBv6jQ372bahj8zJjbYp72PxRIDwbsR9kzG8NR/6KAC00A5Vs4T63pBW8IZINn8Vk+Q5DKpASq8rdMzCltLNNLBSXOt8C3TT2GzYV5MaAz0pyS3s1OiPZGX21q1QJtAqRrsmdzoK891sYDmFE4vCVMm/SdyqwcScknHw4rMmBYT+9te4dsGru2Ls9MLizSOpahgyTVKKzEfle1II2AEjIARuFsEZs9zHaRmpddcxp6JzX/moWCdVZlbSdt9IJbak9WC/ucZnHuVFyf42E+xmYYPXXvn6aLX2NppFj7KZ2l5lPSYin7/Ugz7zhpz+THa8YCxMjrjFw971g/lf4ckufyTZ7WHjTDCe2iv/edc7CrVpCt1lV2jnmOfiJ6LylYV8OFnbvp0uxJ+tIWhNj03X+wfyOoSeUDo1UtFj74+3+LP8rUStW9Ch2eSVY8v8mPkf6mrsuFEEt2nxrAsX8hNulvVY2920pl6GcI8nXev8IfAd3KiDY2wzYt6lGSPDegguwBgoETJiuRn0MS41zWwMVjGE6GKVz+8Okr4HgT4p9g4b6H7KDbClMGZp0zvu5mVuFdyWycgu3zlHjxorJdGbmenq5HZba20Mb6RxkU7rP5VbULlq21rlEtl5IOq/NsYDyG+08W/2PEku0Ul7HO5zRPHGATps30TeSs9fN2AM99P9sUJ/VrphUEKR2SKl7Hphdz64VBfXoqnHtZg3ifWYUbACBgBI2AEMgi05rlMgg5PMz2vKg++JdBJ17zFUMWeqnUSUvfZPVktS2lG517FM+eyz+M0D/Mz/Bzm6FsTped8yRiiGp8l5RkSWsIDt9rAV8IxZrDGDWPLa91zmqne3+qesldGD4WPrVPOuf+ssZOuB+k5d710qWt3itMqGwEzqU4vXNL9RLxr23Sdr/R9ic6SWe8P5KfNRTvrtkvYg6Lttvp8RDbcLF8jyYNXutB2IQze3XHpe4V/2+CpGPWTHcOyfCF30pUue9fjIJZz8h4pCP0tMB9hmxf1WOwIHqKwKnJEcwnRSNmsVd8tkIvBj0nh66YwhWdOujSTzPFThug0c9KN8kpnKoN/lJpTKTwF+7VPcAYb8fCUD0MDGNakMAaGGCgy+lAvGb46DzzKh/qLfJpxVTtRfNPAEfEMihhI3M4CkRmusNu0nZG1ZE72Q/H0tbXWnyAUORiVh15VILubaWsqLwvpdyoTxs36qS6FFH2rq37SVYX8r59m+lr0IeRPkvJf0y7W9sVV6SdwpOyMI2MLnMBnEeYrsYu8j1zJPXr4csTkACNgBIyAEZiNgMbXD8YSKX7N+rRP9Kp5TgLT6aU7J+PQfxaVMmOAYo2dodaerCfB4NxbdPxUbqzz2VRjOGLzy4b/J1296xeF962d0vj06BlBU+UJvpYrfdgzsY9604rQjeLYsxD8pa5f8awg1nW9mKyQSdIl2M1dL6XW7sJqzVq0D4YlZWvKWZseWZPtSuVutWmlmZWv0rPH+oPShUEZAyDtkYt+hQ5HpDSME5N9Pst3lMFxQJ8e7PsYs57pqozcyi81hmX5jtVYFDK7HqXfUX9diOVk3j0lqvaBPeGHosNs+8tjCqQLmXTULkVYXyXXvErPqaCncj+tA9uekMOJvcGJgCSSAS+DESc0WpN6iUMGVuChvBS9Pyn/d8plSx0msSmlOvrGg3T5VXFMTAwMGVo08Uh+TO6tPEq+GHIHnyoo7iramfSk/b0uBQxj6kv0bxX6RDcl3y3bGZovbWuvpA8GrnhyydF/wnjSQxvso2tta59QmJ5ycTqPgZarOqouHvy0G55uNYlwaAibh9gFv8pz8fhDWl3kis5dirDBMX9OevGmcWwoMrjACR7JXYw5+kvO1n3qILmt+Sp0tWsEjIARMAL7IqDxd/H6tE8z5gldRMWc2GSLsNXzpPJgnflE7qCsZsbhFz+bavZd7AlapLDsnqyVTjdjcy/rPYwSNSkfDB085I1XO7lPzfniS+Mr3qXlqXXteGJtz1pgiKpyKBKjCq/wUkfBTxx4nIXmYIeC4l+yXnqqpFHewXIWTDZbT80tW1exOenFm2pX4pts0+Kh7aNOjA1N1SKs1cfFD76tcUth0a5avAhT3GCfb2aW5Wum6folI/rnWBuo+pF4U2NYlq+ry9S95G5Wj8qr+wryKObZvKfKMBWvfFrtJPgVzpg9aH95XBjZtMegF2lx6eRQbOof7o5/nymoz4pdpZcSGKBo5FxUxiCJj4b1nRha+iiMTsZABSEnS+g1hz8rdzO+LDYlQzDorewZCoFHX33NELGI9eLbmUq15LQaYFx8O0PJlW2NdoeRaw5da1urXnsXXt0Te31l5yk6Y9y7TiRPgnnidjRZd/i4rfoj9dMjp4d9ddDavphNPwfHA+VXyRj7f50o4RaYT2ThaCNgBIyAEbhjBLLz3BBEmfTsn55r7uNEXJPY87B5I5x1RP2mgPxs7D6WW5/gk7/aM8llvTG5JxNPi5RucO4tcb1rE/LTxcY49otz5vwMPug5uzytwh3fxJqMMg9R8LDm5XVFTisSBs7sE7qvMCr4iPbcF2SxQ6kl6yWwqdalR6XaP2BO2fq0yabPtqtsm87m26dzhNHvjw5Dqb1N9fkqfZYvMptwh8oTyaKP0CcmxzDxYP+Z5FMZ6rEuMppwt67HKrskltm8J4pQRW/e38LIxyQS1uOmIljn+zavTR78/xyoFAxSNJKD4ketwop/oaviFfsWp9XIFopB+eHuAn+z2IivmoQbOC0tDZNxdM6lMpaku/h2pkK9Er5zT6uBxcW3M5RU2VL9UHxHbU1h3UUNRj8WeGMGmWtta+9UtqOJVmEM6FCMVfjrcY4bSJgwUXNln3CSHwRe4a8CdvpZ2xez6efgSFED36nJbgvMd4LWYo2AETACRuAGEMjOc0NFzaRnLdFcT1SytIbgzx1Yg9SGPCLK2oKNcncj/LmiY402uSdDVocG517lxbqRq/r0Uicdt090RRnmzPkZfJC/pDyk6yWVI/a17DW7J3dYW0DoBnHf3BNUgcmfPfcFWexQdcl66Vxrd/SdUzb4u5RNn21X2TadzfegNkh//VbXH+RHPmH0I+qqtW9QOHuJqT5P+hSfZGWJU6yUqUvox3hQ9fniRv+veRXeN4Zl+Wo5Cc/W9XiQ7lkss3lPFWOX/vaIXFUYJoa3cml0FclPY+O1vpcPIQ8NUOHvdVVW7QiX+43CaAw16f6rctOcoMineV+xiBcDY3NTRyMfPfFXJcz9YIB4nmD9qPAA9Dkogw24jBlUsnrTQbeQk82v4lM9X0M7W3JajfJdSztD19VtTXXJAMh40ZqMEN6ha21rLKBbY5DKTHkZF/lcQHNS5ilWPW4oDodwUxsAAB8USURBVB4mRviy/ew38UOk3Z2k16q+mE2vgqRwbBQ4yl/h2wivvQXfLTCvZdpjBAIB2pcuTmpwsbnjinYZbHaNgBG4cQTU7081T/YhyZjTGnekD/MeawvGqGrfFa7C6nWJ/Nk9mVhrirze1SFtD3u3o7FQ+fM6G6+zxsGB9JyvNCl8JX9JedraH9/9SUFfSAfWsk1iP8qfrYUxAjyYC/gsDUZOrsCqma7Pv3pfQF66jvbdCkthV3Rdsl46y9odEOeUTbyLsVFW2XaVatNZvUtDoU7eFn84R/sGyUz1+SxfZNRwPyr+eg8TcZKJAZwHDfUhMPlp+y3bUPD3uPBm+soU36COJc+t6zGF+cy8e+A5Coox9ChiacAH798/fCu8VBwV+U4XG04MY18rvLVJ1T3fXqDSW6+M6h5Q6AgQjYXGy1+MI68m3UdjQQ5E5WIJrfjkcs/fM/d+30jhbLQZcMc+9i+WBxIfAzgTUy+/wsNKjQGNvCkvQDNpMYiejJTfFDYYUrGeB84t3RSewkZ8YM+CICaxlpy5NyVfJr7Bb/KFTPGA8cW3s4a+6Mr3IXvbT4PvatoZOqs8i9ua0tLXaYt/lr/qt8jsI8Vfc1tjTGg+lKDcjD11v5GfPscYQvvAD+Fn0mmNnVXMwI94WSiDKR+2TqcbEJcKVj5r+2I2/SSOobB0AmMeIvEdzNYT9gbPJpiHPLtGoImA2h19t17fcK/4F3JH54CmDPuNgBG4DATUb5kvUuvTPo2VPjvPDe2NUukj7zLeMA8yb0LMg28UjuGJPIjrI06nYZipSH74Yq8wuCcr7IfCPzX3ss59rYv9XRD7q3pNRKDu58z5KXyKfunyFD1C38ASPdmXVpvphsx38IvAjPG/Wx5wQVaTkAHv4L5HcaRZvf+UnMVtS2kXrZdKnpvtE5vAZfzKP9suFmODHson1U/El2rTWb1L3rEPo6xQX19K9Xnlm+J7yKYqd9r2IdnoGTo+lf/INhRyccXPmglco9/VYxjxQVN8ip+j42b1qHznYpnKO8rd5ypPjHHYwuq1Zx9fN0z8o/NbbeTrJjzXfVH4tdx6smrqUuLTRj7SKg0VhkHiJJvnpr5b+ks5qpNm8sekVGehMCp7FBvx0BgxYG62YSn5Ll5E1QU4oafoPNjOUEU81QQtL4aXI7y76ornJtpZKTtlOWprKmM1mMmtBqJyf5B79ASixN10W1MZmcyeye0dr7ptZOhe6WMhxj/WtRaYQ2nuNXwrzO8Vv0stt+p11njbLIfSxmKZB5TMbcyDR2NSM82QX+lYbNX9UPeMeYyHJzPAD+lGuPQBJzZfR4tBha3CYW36Mb1PGbcFRmuxWJv+lHjdcl6qh9FN0C2X3WVbj4DaD290YMyrHjrKxdiBoYNxmPF2ytB31n2B9J29RlWazfeJwslkBIxADwJlTOHVZh7kxIOMHs7jIPGPzm+PjpOcPeToe3wbaMRAfLQg3kDuqUVgSMMAwKmfpUQDio3AUhnddBjAJo1g3URnvh9tZ8KZSQ6cUga+UpZbaWcU56itFUxYMLCo4dUFFjm0p+ZTXd3WdA9tjSdVWxjlwiBBuzONI7AV5uO5OHZ3BDSG8CoST68ZV/g+zez2r7ScsngjlxManKpg3OHhwmxZSgdVDzcevBf5y9NtNpktWovD2vQtZc5/swqjtVisTX9++G5Kg2tcn95UBVxrYdSPq5NIcuu3CuTnbar4sxHW/Owlxujc+4Il66U91u5jGDnOCNwzArFW5SH1XBqd3x7NlXYCfgYkFuhPtspLsnjtFqNFALmV6FPL4Sgn3zpsnUikXLqYSF7rqvyFr6UffCV+09eQJZfXtzeV2VJ8n5vBdlZwYsPIaQ4m9FTbEd+ttDMQ72trbKbBDTeuV2Ck+xaBmQLAbdN2IXkX09akC2MU5XzTKvyCG8mKPr3qROCCrK8qyZaYX1XBb1RZ1SfjK6fsMax9N7eYSscDLwyFrU2YwrjHcDibJKv+fEhJjG5TfzA0O58lCaRbfOu4lXwtDmvTt5RZcKP8XxQdFqRuJ5GcVRitxWJt+nZpfLcWAdXHxawZ1pbF6U+OwC/KkXUsD7T7iLlhdJ5R2rPtC5T37DWq0uyydu8Dz2FGwAhUCNDnoDjs8XCX+FV/HZ3fHidknJqFARWlWagNfutggVLVR2OVbuqpywLRp0kiTJhQjkjhNAyevEwRk1GvjKmENxjf286EJZ0NnDDyxcQOZhl8genq2xmFUNmP2onCPiQuSffQ1p4VLLY4yYco+nHILKLtdBAIfLbCvCPet1eGAONtGMibqmN4/0pjFgbAo4cQTcYxv9IyB/A6xNmN70UXytJXnrU4rE0/BmMmjs0o1yraCKO1WKxNvwoDJzYCRmAbBJg7dDH282mfL+XGSZuP5Ge84tvBffOPolrEmMDp4lPvP5esl+5h7d6qHN8YgTMjEEa+zFgyS9VHs7hPw9x3guiggRSjzORptSEVlZaFMd+x6X3KO5TuVsJLuSn/bEvxrWDQKUdvOxNP+rRaR151e+/tDBDuqK3xlLc67VlV/vofBvgY7NdLu00JW2N+myjdT6k4Wfy2p7gxzxG/iDSO0RdZc8z5ZMOivJKJvpROQyej1+KwNn2yCLuzbYHRWizWpt8dJGdgBIxADgGNuRj6eOjPxTezqu9myWU/ldqUi+9c+89Z6yXpyf7Y+8Rc0zCXEdgKgecIUv+LdetWcg+PN5O0kSAV8ugEEaJL4bOnqXq1KTK2PB3Ym88lBqrsd1nuoboQHkPtbM5ptV7x99zOAORe2prKycJtyxNlfOCZf3JedfqIOrhV2gHzW4Xq5stFP0kU8ujbdaRRWgwxnK6gD2PMYwNXL7Dkj7Dq5EW5PzR5lOZkpHzZfHHC4ogUtxgHhK1Jr7SjOB4pu2OAdFmN0Ros1mK5IzQWbQSMwJkR0NjC/HLSfZjynLVGFf9J9TtzlTh7I3ApCLCO2nIvWZfr4ox8tWb2GAEjYATuCwGeukIM+D9UPv8YASMwhEAY8NjIDNGRAUwbGU7nYUivHvTIpb9xgrt6wKN7DHxrPtmg5NtR0YfTJLURsiN9EQ4NGYvSS59RHBvyd/duiNEiLBoFXJu+IcpeI2AEjIARMAJG4FYRKGsX1qkc8ticHm0u0QKNgBEwAkZgNgIa7Hn1A4MF334xGQEjsB4Bvp1Uk/oY39fjtGzzJDfGdYx+GPugVZ9seBCx6S+vTw29ppvNqIVDNlGDr5U+iWMj+e7eU2LUwmJBydamX5ClkxgBI2AEjIARMAIXhkCsO3c52PH4wgprdYyAETAC94zA9yr8F/cMgMtuBJII9H2LL5LGiar4UHqEfytP97MfnNyDqlN/MmAt+mSD0pH+55CDwATxz8KD33VSHP8e3PuabkP2EhwayXu/aRjxi3EMAU1X5aEssahtRlX5KL5pfI34XxXOq9W9tDFG58Cyt1wONAJGwAgYASNgBG4aAT4J81+tY4be1FhVeBv5VsHnxEbACBiBTRHgH9heacB/oWuXbzRsqq2FGYEzIaD+wSus5F4Z5zpqRFi9cKJPFV4M6U0Ko9Ogsa3JPORHH8V9OhQ/N1zyMD5ywrAuQ58M8tVFVJS5yRZhgzLmphf/YhyVts+Id1A4Jyz5c7VZ34QS/6YYSd5JsWxWlP1GwAgYASNgBIzAXSHA2qf74HkzAB5vJsmCjIARMAJGYBUC2mT+pIsNOSdXbORbhaYT3wEC9JE4idcsbpxAa/Yh+hSnwjDGNYnX43d7ktrMaKafcj2Xvhj+m/SJbjCIEY7eLBDn4NCUFf456S8Jxz0wmoNF4Nd016ZvyrLfCBgBI2AEjIARuDEEtHbjTQ1o7edYHqT0/NrI1wOKg4yAETACZ0SAD9p/w+ZdV9cgcUa1nLURuDgEMHTRX7rEibquQY8TaBhgalL/wmDGtdkJvFr4So90Q9eWvohU+O+Ey22+wjoHB8R0aU76i8FxJ4zmYNHFkfu16ftkOswIGAEjYASMgBG4HQR4QPuD1jG77fMe3Q5WLokRMAJG4PoR0IDPUx0G/dfXXxqXwAikEYg/JIhTeHVC9QleW32viz/FqKn0lbdyeeWhIvmfyMN3LV8+hFSGMcI49VXLLnwYZPjThlWv6kY+J3IpC1dN0p8xI4WDeO8BxzUYnQTLuvLsMQJGwAgYASNgBO4GAa3DeFjKmnS3V3UB0yf5QMFkBIyAEbgsBBj4Oc33ta53l6WatTEC2yGg9o2hDWLRA/1LYbyy/qPc6jUG+kAJ459wu8QpvL8r/rlc/mgD90+6bxruQjanY79SPPSxrtE/vqi4LuRHevOnFSwKIb7ZCW5v5MZ37CZxEO9N47gFRhW6Dyc7R9vUBliWrOwYASNgBIyAETACd4RAnOIb/F7yFlh88P79+y3kWIYRMAJGwAhsiIA2kZxa4rW8XZ/0bKiyRRmBi0SgGH+eyb2413IvErABpfbCUXIX/fHGgJoONgJGwAgYASNgBIzAxSGg9Q4PnX/U9aH8ux7ieHRxpbdCRsAIGAEjAAK8bvhKk0Cc3jEqRsAILEOARdXR9+2WibrrVHvhyEJ318XuXdeaC28EjIARMAJGwAhcAgK8lcFbJLuveR5fQmmtgxEwAkbACLQR0ATAHwd8rVAmhM/asb4zAkYgg4D60BPxYSh/k+E3Tz8Ce+Io2TbA9sPuUCNgBIyAETACRuAGENBahz+KY2/3wymK45N8p0DZeRgBI2AEFiCgiYDvbfEdrb5/EF0g0UmMwN0h8KyU2IakdVVvHNfh59RGwAgYASNgBIzAHSKgfRxvQvA95T+fqvg28p0KaedjBIyAEViAQJkQmBjqfxBdIMZJjMC9IsCfdXym/rP7qxE3DrBxvPEKdvGMgBEwAkbACBiBbRHQ+vMTSeSwxp+2lTwuzX+8MY6PY42AETACF4GAJgk+1Mq/gzb/NfQidLMSRsAIGAEjYASMgBEwAkbACBgBI/CAgPZsfDLmZ11/kv+kD5tt5HMrNAJGwAgYASNgBIyAETACRsAIGAEjYASMgBEwAleOwF3/8Uaxrr4udciHuaGXp7a0PmR7n7+ug8uvd9fR5deRNTQCRsAIGAEjYASMgBEwAkbACBgBI3DXRj5V/99lwPhLNAP5+RfLf+v6OMLs7o6A62B3iFdn4DpaDaEFGAEjYASMgBEwAkbACBgBI2AEjIAR2BeBe//jjVcy7PFvJ0F8FPGPCuMDiabTIOA6OA3Oa3JxHa1Bz2mNgBEwAkbACBgBI2AEjIARMAJGwAicAIF7N/Jxio9/jDOdDwHXwfmwz+bsOsoiZT4jYASMgBEwAkbACBgBI2AEjIARMAJnQuCsf7yhE3P84wgn6fge3k+6P+u/Rip/TvJ9LvduXtd1HajGb4z2rtN77Cc31kRcHCNgBIyAETACRsAIGAEjYASMgBG4QQQ2O8mnjf/RK64YG3RhwDsihWPc408vMOz9U9czhf1L7llIeaP/57o+PYsCG2RaytCSpDDXQQuR67q5tDq9hX5yDS1gbr1fQ5msoxEwAkbACBgBI2AEjIARMAJGwAjsi0Drjze0seQkG/SbLk6z8cH9/xKQoJ/Fy8m8OI2HHzoymhW+z+T+9YGl+v0n4br+pevPjfDdvcoPQyRl/1T+d7tnOJKB8ncdXHcdHLJ1mOS7mH4lfW+ln2xdRyl5wo8xMf7Nm1GAe8a7n7jpUKrelZY6iXH0mfxvuVd4jMO1WIXNGlsa/CHjuwm5vXxKk9YRAeLnD5CCnsoz+Y/nSsNDmr/Irf9IKQQoLI17Jm/xUJ7IB9lVv1D4UT0W3sn6yfJFmewaASNgBIyAETACRsAIGAEjYAT6EKhf19Umg3+V/VruDzDKZfNCGMa4SUOfeP4jXjZkpIMfOch7J7dFCmOz+UrXH5rx8pP2d10Y2442qQrfnJRPtQGVW23ayv1B7mSZt1ZGeboOBKpwoE5wr7EOUnWosmX5LqJflTrBeHQL/SSL/dZ83wR+MXbontPLGM+qcbcRPlnvSkM/QeZnjXSMrV/pYtyujU7yp8qCHPEiF72o70qGXO4/kVt/yiDDV3iyOsacA/8/lN9BLsa7n3UxJwyOB4oDr1/lHj0gUtgk7uJJ5V34uv82zQlw8Pmz4ut6lD9VP1k+yTcZASNgBIyAETACRsAIGAEjYARGEXhErDYZGNw4RdfcoGCc4755qkK3g8QG60NdH+j6WBcbRGT0ERuy77vx5Z40zX+87Uu/SZjyqzZhEsYmkA0sG0pOXXAaZjYp/QtdYDmbSjrXwRXXQbYOs3ylEZ29X0nfm+knWex34sPQ1qWXCmie7ov4TL1j0IsTZVU66c34xRiK0amibFmCXy5pMTzWRkLdP9HVNbJl+FI6Sjb0ra6nyrcy8BEgPw97ftE1OA+JB6NmLymO8TiDezZv5PFv0xj2ggKnbj1my57li/zsGgEjYASMgBEwAkbACBgBI2AEehGojHyK4fRD38m5NwrHcMUGbzOSPF7NbW1OG8IxsD1v3O/pZfOHQRE3LjZwbJKXEDgtxcp1cP11kK3DLN+sNrhjv7qlfpLFfms+TsDVJ+5mVewwM2PXf1Tv3TEHo1PzW5zZshwkC+MVDzv4TmpNCudkYK1/lk8CsjqSF3l3DYmEMzf1zkPSA10Zr7n6KIt7Nm90aeUnHYbyzpY9y9dXPocZASNgBIyAETACRsAIGAEjYARqBB4XH5uM1qauhMeGi/j6lF+JW+VoY8SJCL7XFxukNwojD04NRb7y7kfK78P9pM+W7DqYDdnmCdbWQTZ9lm92AffoVzfWT7LYb83HA5OvhOWPcnmtM8Y9TnENnlKbaAAY8ziBHLK67E9KQLYssPPw5d2IzCIyzZfSUfmFrjzk6dJvJeCZ3Dg1FzxfKm39GnkENtxJ3OfkLV7yb80bCsNACHXrMVV2pcvyVZn4xwgYASNgBIyAETACRsAIGAEjMITA48YGZ4iH8KdjkREnWRju2Kx9pAtjHd/k4+RDixTGRve/cuvTfPJzUoP0UK+RT/FsiGPDx8b2F4UdyUfANZHKEBvcMbVdB2PorIxbWwfZ9Fm+ZnGUZtd+RV7KA0MQJ77oV/RdDCe9/VBxZyHps6qfZNNvzQdYkvkDl7wYhH6Xn9dqOWX2YwmXt00KH613xVNffcTptoPieeV3LmYY0hibkfGlLsZb9Oz+QUiKT3ImdZTsg/gwLOLtG+eYTyDaZU3i5zXdrmGtjscjnkncxTM778hEaek31eu28rcelOk+W/YUX+Rp1wgYASNgBIyAETACRsAIGAEjMITAI0XEporN/RBlNorw8J29f+hiA8v1b/nZBNWkewx8B7m1ga/cc5ohwng9sCbx/lEX3/HjtB/y+WYTpwDZXN0CuQ7OX4tr6yCbPssXiOzWryID9Sf6EafL+HdS+i2Gk1YfDN4zu3Ox66qbTb81X6WHsMWYE4YgMGdsHHpIkar3SnDjR3l8otvKSFuCs2UJKTHWP6Mt6GK8ZVzGyBcn1uDN8oXc2pWcro4RFye54z5c+KHI8yAZlBHj3KQhWjwZ3NN5o4hkcoISIyPYUIe/6Jok0ompWT+9abJ8vYkdaASMgBEwAkbACBgBI2AEjMDdIvA4WfI4STHIrk1J/b0mmHTPaRAMdxgMOAlCGCdT2NhioOuj2MR934nEMBgnMiLqC3m+5kZySfetrs/l/4AwqITDh/FiKM+K9wp+rq4OCv6vC7ZsbKGXCn/34L2638k6mChRNn3NJ6z27FcHycdwQ7+p+mjRH4PFE4W90EUfviaqsVuodDb9bL6CNW0frBkXGQv5pl7rX1kVdlDYZL3D10P8GQZjZf3nFT083aCqLEoT4y8GrDBGBi9j8rcKb7aHUT7xDvXzIR1fKo//U7q63cmPUSzkNA16YZAO/QZdyaCNT+E+J++DZGLYqwy0RT4PtI7qsUepobJ3WbN83XS+NwJGwAgYASNgBIyAETACRuCOEcDI93ak/E9LXLwiO8LaG8WmjA0bJ/Hwc0qI0xfV5qgnBYYgNqixqTvIz2aY8MqgF2kUXn8XCX5d3xW+ikX3bA7ZRENPHpxtfiU7NuhdgRVeiud0R5d4dW7otaybrAMB8PcmFgU3Tog1DUpdnFL3RVbUbzPNueogW4dZvmaZuv7V/aohEOM4/bJJ9DdoVb+5gzpK16Ww4AHHp3JjbOCPLDA+YcypjGe6r8c9hfVRt95bPErPuMTDleY4k9axIYx8ukS/pQzPdMWptSm+pkGwkjegY8Qxjv9BN4wbjN8YH9+UC6yq/BSHHpR1kgrvJO7iS+Xdl6HSxpzFaUf+Yb63HhXeVz9HIrN8RwkdYASMgBEwAkbACBgBI2AEjMDdI/CYDYkugOjb0EdY32auBk/pOWn3VO6ndWDb80RxyKqMeO2ohzvFh8Gmu3njJN5Piu/dODVkcfKl3lSK/1fdY1hjc7gpSWZs1FtyS14YNOecojmI/ybrQOC8Utma3/LiFUXCOAFE/Swmpb+oOpA+qTrM8gGMePfsV8inz9Evuydnoy+6jqiINtVjovBL1XlJTtvHgFWT0mMcwuDNpwjAnFdGDwqbrHf4mqQ0GL4Yg7snANM6Km3wjo21jG+Mx2Q/ygdDk5SmV8cODzJbfVvpwA7CgMkc8kTu6JxUcT/8zMF9NG/EKV+Mj7jdvoHhkzqs6xG+IPFPlh3eLF/ItWsEjIARMAJGwAgYASNgBIyAEWgi8KjcYByL0zvN+OpUlAJq41kzsuHndMeTxn1441QVG6KQxcmMPuI0EUa5Oi/5kcnFpneK2Fxl+KbknCv+FuuAzXqc+jkXrnPyXVsH2fRZvl36VQMQTnzR5zBuNOlL3WBQyRpSmmn39mexG9Ijm34zvhjHenA+FIwx7sX4iN6Zeq/LJxk8yPhYbn2CT36McTGmZ8uCTHj7xnLioGgTWb4qUULHim/g5xOFx4MeyvRc8nh4UF8Kg+dFCauMgvJX84fcbvs+KIxydHFX0BE18yaSE428mjuGUUuIeKfqp+LP8rWE+8YIGAEjYASMgBEwAkbACBgBI9BA4FHx88oYG8sucTKvzwjQ5funNih9r2BieGMzeFA8m6qjzVaJYxMEb71JLeHBHy7BNUkmaZDNhqs6YVJHXp/n5upA9UK7aNYdRj+MR91TMJdSW2vrIJs+y7dLv2qATf9pGWFVNxg1uFp9sZHm3N4sdkN6ZtNvxlf6AKfkwujW1Y3xqxonS8RkvYeAUl8YvbqvXDOmvi182bLAzknqPj2ZCyhD6JnlOyR1hI9vQ/LPw+BRUfHTTqvy6R5jH9++a12KZ5yJuOAlLIW75E3mLVkQMuP13Cqg/MT8GfhUwZJLX5qqn0OWr5mh/UbACBgBI2AEjIARMAJGwAgYgS4ClZFPGww+sv5WLhvDiuRno8WrsnyQvCLCdL3XxWmGJn2jsNZrtrr/qjA0jQVsvrqvYrEJ+lYX302KUyIlaeWgW1NGFSheTmvEJpZN4KUajip9p35UnpuuA5WPeqZ9Db3SPQXR7vFr62BG+lRdq8C79SvpSv/GmPM0gC1hGIT4U4OL7E8zMO4dq2akT9VRVp4wZQzj5Bm416R7XuP8UW5z7EvVu9JQf9QXZa3ShKsw6hCD1EFuqiyFl9NtGMviFVnS980FWb6UjuQtgjfG9CpAP9n2iI4tbIuALO7ZvJnDWifGhQ/jGnnXmJO3wlNlz/Ih02QEjIARMAJGwAgYASNgBIyAERhD4IP3799X8dposElhY/dO12+6nuv6WuGtzb7u+X4Um8CusY4NTXWCQi6GAzZrfxUf8mrSPRuiL3U1N3NHfHUCeZQmNpzkDaFrfUpM8RgYObER+cNTkcLIjw+59500LFzbOCWv2d/ki9yV/lbrgLZBHXH6ptUeouxbuRdQB9k6zPLt0q8KThhQ6Bf0EQg/xqJWn69iNvy5wTrK1iWG7te6mmNf85uVNcrCaLLexcN4CF8fcQK7NqjLn9IxBImfMZc0EOP50VxAxBSf4tM6NuThjbx78YEBknzGFTCoTnXLxfj4RuH1d1HlT+EuvphnRvMWH3k1HzyRP3NM9xRfquxKl+JTHiYjYASMgBEwAkbACBgBI2AEjMAoArWRb5TrwiPLJqn6/pv875rq6v5qjHxNva/NP1QHCq+MFXIro3C5P8htnlzarLiSS30vNrRupsiFCxJOGEeeya0NQadS2XV0KqSdjxEwAkbACBgBI2AEjIARMAJGwAjcEwKPbqSw8T0+Xn07J2FgbBkZz6nMifM+qgMZczDwYUzidBj/qMuJGk5bNk8y6XZTuuc6mAMkp5FaJ4/mJF7J6zpaCaCTGwEjYASMgBEwAkbACBgBI2AEjIAR6CJwKyf5MCTxyhOvqFWGi2Jg4vQYxgyMS7y+9ZvC69e4dG/aCAHh2lcHv0v8k24W4v2gG+b70yEg/KkT6obXp3m90WQEjIARMAJGwAgYASNgBIyAETACRsAIXDkCN2Hku/I6sPpG4KQIyLCH4Zs/D/hQfk7VmYyAETACRsAIGAEjYASMgBEwAkbACBiBK0fgVl7XvfJqsPpG4KQI/KLcPrOB76SYOzMjYASMgBEwAkbACBgBI2AEjIARMAK7IvD/EuqkdQdqycYAAAAASUVORK5CYII=)
$\displaystyle 0.5 c_{l1}^{2} \left(1 - c_{l1}\right)^{2} + 0.5 c_{l2}^{2} \left(1 - c_{l2}\right)^{2} + 0.3 \rho \left(- 0.037 \rho - \frac{1.0 \cdot \left(1.7031322 \rho - 51.0939660000001\right)}{1.0 \rho^{2} - 40.0 \rho + 400.0} + 0.042578305 \log{\left(1.0 \rho \right)} - 0.0530922164415325\right) + 0.5 {\partial c_{l1}}^{2} + 0.5 {\partial c_{l2}}^{2} + 0.005 {\partial \rho}^{2} + 0.00085206629489322$
2 2 2 2 ⎛ 1.7031322⋅ρ - 51.0939660000001
0.5⋅cₗ₁ ⋅(1 - cₗ₁) + 0.5⋅cₗ₂ ⋅(1 - cₗ₂) + 0.3⋅ρ⋅⎜-0.037⋅ρ - ────────────────────────────── + 0.042578305⋅log(1.0⋅ρ) - 0.
⎜ 2
⎝ 1.0⋅ρ - 40.0⋅ρ + 400.0
⎞ 2 2 2
) - 0.0530922164415325⎟ + 0.5⋅D(c_l1) + 0.5⋅D(c_l2) + 0.005⋅D(rho) + 0.00085206629489322
⎞ 2 2 2
0530922164415325⎟ + 0.5⋅D(c_l1) + 0.5⋅D(c_l2) + 0.005⋅D(rho) + 0.00085206629489322
%% Cell type:markdown id: tags:
This is the free energy expressed in the order parameters $\rho, c_{l1}, c_{l2}$. Next we have to transform it into coordinates $\rho, \phi$.
%% Cell type:code id: tags:
``` python
transformation_eqs = [ c_l1 - (1 + φ/χ - (ρ - ρ_l)/(ρ_g - ρ_l)) / 2,
c_l2 - (1 - φ/χ - (ρ - ρ_l)/(ρ_g - ρ_l)) / 2]
transform_forward_substitutions = sp.solve(transformation_eqs, [c_l1, c_l2])
transform_backward_substitutions = sp.solve(transformation_eqs, [ρ, φ])
```
%% Cell type:markdown id: tags:
To do the transformation, we use the substitutions dict.
After the substitutions the differentials have to be expanded again.
%% Cell type:code id: tags:
``` python
free_energy_transformed = free_energy.subs(transform_forward_substitutions)
free_energy_transformed = expand_diff_full(free_energy_transformed, functions=(ρ, φ))
free_energy_transformed.atoms(sp.Symbol)
```
%%%% Output: execute_result
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAADMAAAAVCAYAAADrVNYBAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAC8klEQVRYCc2XjVFbMQzHG44BUkYIG+TKBmSD5DoBZAMYoRc2gI5AN0g7QRs2gA2AbBD+v3eWK+t95AXSS3Xn2Na3ZFl+GWw2m08RBoPBULgv4EX/GemH2suvkWyPNR7k11P04ygiJHAl3K+E/xPph9ynAAhiIT+XKel/XeJkbAg71eCohobrM4v/sg/fPnlkc6lx73XGkzkTw5MY1pp3gdNdmPfE+yA9lFyGGEwm7Lh43pH/n7B/OBjV7bk8I0sHhw8HowjGKsv/ouMdvyedOg1q9asG5TXRHjXcm+t33DdkM6STngnBvR1poLPWhrOAW8RgeF86L7+M3SIvA3NmAtH6RjOGV5onfY0j70GyC+3ppKabEl5pfPZ8aU0iTwq8BKvWLCTOvGq0tlnR7jVunQzGKDPTQaBFuzTatllyPAuPnk97kstTce7xrAX4y2JqNEPiJBmoCTnGS9EJNr9BWi+MngyQ2cIhT+9aJ91FIoWjnAuHvQ7RCNZ8n1pGESKrPETZ2SD4CE/AxWCQX3qePmvJcMI4XdjWnq8R8KMmPcKTvHwIVTDGKELtVXU0lObMaU1W4r41i6anaZYeErmKNHAajSctfO1r5UhID7VX1RGhEYAB2fQt+bv2d3LohzHoQnt+QzfN6Cq+AyVLtTDobE3A18pa9tZGjMEYvmm+ELLqMol4JkVVy5Rh6vZF+0xPgbxqJrutkPi4zCfGlHDonEtn14P8YjLMx37TtUapjMw0KAnuz1jrqpVq/S0a1X4tOsEONXMXcgaF88CpANfi444AvFmzqLOidPz0DgYdSflcRskkR3zXoRv+U/FS210wEdH+n9x0MW6j7VJmXhe17O+Lp8U15dh2KvByMn11Rd3FPgbzLGqu3YKz3OT7UqLLHeUlDDobIdE55d+NDN1IdJdJUtZyexYRxVvbq3iuvFzbehuf6JwK9or3pU2f4eHXqH2t5EAcIwFVnyWa86eC0fc5J6davzqiLfHjG02HTlfz7Q2bYtPrSxjjxAAAAABJRU5ErkJggg==)
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAADMAAAAVCAYAAADrVNYBAAAACXBIWXMAAA7EAAAOxAGVKw4bAAADA0lEQVRYCc2W61FTQRSAgUkB0Q6EDjBUoHQQtAJjBzL+Cv8c6ACoQKADtAKFDqQDkQ7i9+3sLpe9ezM3IU48M5tz9rz2vHZvNmez2UYJR0dHQ3gj+dDfSvm69sSyzdm7rFvouzKOrZKB0id43yP/Zylf5z4mYBLH0Ncsi/4Idiat6XQ6Zs1Yw8Trg9Gf9NFbpQ5nXrMumz7LzuyR5h0ZPzym24va6aW1WqVb3DlyGcpksmBB4veC+v9E/dnJ0MW3RGaV1g7PToYMdknov3jxBsuUk+Cd1fcsx2ufPWjDe3MIveh90zYD9nb6gKWfbZY+W88w/BaUnfGpmxsMjk/R+Qj2kBNon0jxMesG2gCWAmz1cQAO/qE966bDmYV82ZTlZGIQ7xDqoAroXCrwsIitYhgxeFZP2oAWBuzHGI2T7+jA79wQnueUcAVDmXYBBmyspAGIrUp1/uFPkOv0FSuBI3aYNmC7Wju4odJJniNp+lIxddmJeQKc6yfkBcxz8Gfwly2ZEiyr4Oi0DOELHnSBfN4Yepd6zbcOE+DTAnjuReJFnArT9VqahAkb99WWRhD+1wmjw9YK1UCjPL/oe/ivpBj3Ht45pkm3gr3wxlAWykfGDrQKBM/x8q/XG+gwTSEZGAlaX9UkACszgQQG3hxJi3CGY2c5AHRTP7FrWF9ORgZs7bLLRGvgv5UH9HIBymRqRon3ASJ1T94ejkLFwD4M9+AshzaRP+DcTY1KiHp2Pb9MkadPX7WuEdPVvT8JNv2jlgBDHwJflOp/LfhWyoAdr31W6tZp7VB4aQxfQ+cKYpcBvuNi4J6ZXibpqk/4AbBrxTpIwj4YBwZvtaykLT6bZ4d8h5UC7FINRUHPLvu9WhoWGbPmIXaoeV+aspJ2HKtdiYrl3Svte+/LZFpf1Q5P+b50yAObJLw3+qxClNvlH1WF+Ux9PylSmUzrq9rhrzPAQn9CwPNGZxT1+3Y5qMcitP6tPHkA1ETRSvmB9HX5yj4/texXCjGoEbhXMjE2HyBjbMX2FyjLoHFGrT5SAAAAAElFTkSuQmCC)
$\displaystyle \left\{\phi, \rho\right\}$
{φ, ρ}
%% Cell type:markdown id: tags:
Now the free energy depends only on ρ and φ. This transformed form is later used to derive expressions for the chemical potential, pressure tensor and force computations.
%% Cell type:markdown id: tags:
## Part 2: Data setup
%% Cell type:code id: tags:
``` python
dh = create_data_handling(domain_size, periodicity=True, default_target=target)
# Fields for order parameters
ρ_field = dh.add_array("rho")
φ_field = dh.add_array("phi")
c_field = dh.add_array("c", values_per_cell=2)
# Chemical potential, pressure tensor, forces and velocities
μ_phi_field = dh.add_array("mu_phi", latex_name=r"\mu_{\phi}")
pbs_field = dh.add_array("pbs")
pressure_tensor_field = dh.add_array("p", len(symmetric_tensor_linearization(dh.dim)))
force_field = dh.add_array("force", values_per_cell=dh.dim, latex_name="F")
vel_field = dh.add_array("velocity", values_per_cell=dh.dim)
# PDF fields for lattice Boltzmann schemes
pdf_src_rho = dh.add_array("pdf_src_rho", values_per_cell=len(stencil))
pdf_dst_rho = dh.add_array_like("pdf_dst_rho", "pdf_src_rho")
pdf_src_phi = dh.add_array("pdf_src_phi", values_per_cell=len(stencil))
pdf_dst_phi = dh.add_array_like("pdf_dst_phi", "pdf_src_phi")
```
%% Cell type:markdown id: tags:
## Part 3a: Compute kernels and time loop
We define one function that takes an expression with derivative objects in it, substitutes the spatial derivatives with finite differences using the strategy defined in the `fd_discretization` function and compiles a kernel from it.
%% Cell type:code id: tags:
``` python
def make_kernel(assignments):
# assignments may be using the symbols ρ and φ
# these is substituted with the access to the corresponding fields here
field_substitutions = {
ρ: ρ_field.center,
φ: φ_field.center
}
processed_assignments = []
for a in assignments:
new_rhs = a.rhs.subs(field_substitutions)
# ∂∂f representing the laplacian of f is replaced by the explicit carteisan form
# ∂_0 ∂_0 f + ∂_1 ∂_1 f (example for 2D)
# otherwise the discretization would not do the correct thing
new_rhs = replace_generic_laplacian(new_rhs)
# Next the "∂" objects are replaced using finite differences
new_rhs = discretize_spatial(new_rhs, dx=1, stencil=fd_discretization)
processed_assignments.append(Assignment(a.lhs, new_rhs))
config = ps.CreateKernelConfig(target=target, cpu_openmp=threads)
return create_kernel(processed_assignments, config=config).compile()
```
%% Cell type:markdown id: tags:
#### Chemical Potential
In the next cell the kernel to compute the chemical potential is created. First an analytic expression for μ is obtained using the free energy, which is then passed to the discretization function above to create a kernel from it. We only have to store the chemical potential of the φ coordinate explicitly, which enters the Cahn-Hilliard lattice Boltzmann for φ.
%% Cell type:code id: tags:
``` python
μ_ρ, μ_φ = chemical_potentials_from_free_energy(free_energy_transformed,
order_parameters=(ρ, φ))
μ_phi_assignment = Assignment(μ_phi_field.center, μ_φ)
μ_kernel = make_kernel([μ_phi_assignment])
μ_phi_assignment
```
%%%% Output: execute_result
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAArQAAAAyCAYAAACtSEg2AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAgAElEQVR4Ae2di5EdNRaGsWsCYL0RLGRg7AgWMsBFBAsZLEUELsgAiIBHBkAEBmcAGQCTwe7/aXRktVrdfdTT994Zz1FVj9TSeenX67Rat+fR//73v3ciBAKXRODRo0fvSf+Huv7SRfq5rpfqm68VRwgEAoEHhoDmhHdV5S9ytZkTCP/RnHB9k4y/gUAgEAhMEbia3sZdIHARBL6U1j+0WH2Ddi1m/1X0i65/cB8hEAgEHhwCX2o++MxqrTnha6V/0/W+5UUcCAQCgUCNwOP6JtKBwIUQeCm931W6/6k0u7URAoFA4GEi8KmcWN7aWOCh9z3lPbWMiAOBQCAQqBGIHdoajUhfBIHO0YKPZQgLWIRAIBB4mAiwO/vrw6x61DoQCAT2IBA7tHtQC56TIMBRA128Vvzajh+cRFEI5VgHu13sgnFW8ZAgWbPdM+TrsjOQh+jpCUGHrkPr09MTeedBgPGvqz4vi4PLsaQ4Vy8gHnJ/V90vNs+cp/eHlr0IPNIEsZc3+G6BgAal7UD+KTGcC+PM2B9ekV7+o+la+/Lk8plsL+fdWpr6XvSfijadla3zLa1yHKwfdP0kuq8s/xKxF7sl27z8R9PV9kg2k/+sfZTPLjg4L4Vr4Z/OMIsWh9Tal/bhnv76c80sur91T7k5HaQJH4j2+iZ589crE2rRbo4V0bjr45U5SAcmn8Oj8EwXR2Y+V70NC/KTI0J+ulmhy7Sb9c504Gw/oCKL+x867eOiQ4CFCnvL+q6tkxXUsfhWx3lNu5WWLPowfXXWj7Z4jy6v8LjovC07hvo7OFS2GyyztqxoVusnOteckPVy/tnCEyW6P+6TTC+da57JNnrH2vDYoELS0Z1frbIWb9GpfHOs5/pszsPZLtd8VNln+i1r1jesoI5l02HjvJa7O63J6Z24zouBGotdyI8Nd6UZTL/res/y1mLRufiPpuvZlO1m8XT1I9HjCK3SioZJgietp1u0pyqXbhfGS/q9/EfTtfZIPv1q1j7KYwL7SReLSHvBk/qnYvomO+alzSjTRfuUPky5AnwsNtyQRse7NW+mG5Hpaoesa7M+Wb9XppeOxYMHsBoj6g4OH1q+0i66ERsz7aR9ch4OYNs+LrrMj63Uv7Yfmb9TvnWJbnOcb8mo7KBdZ/3Iw38kjWxw9YclnV5+D51oXOO3wnCzLT16szzX+JU8W9f+a5goj7mdOaKsdV66SsbmPIN8XZtjspLpHhvGQ6yALbP5tabZohP/Zr8SjQvzrMtdd8m9E+O8xWvv/ebEtFdw8PUnfXWgT3XNFgXlpQlqCzcv/9F0Pbukg68RuAa08cNjaWIFBioTXHFelWaQUcjT364+Kl5+ULKLHz5d976NVIfF9lFZ1+FQPtiXyT3LoC2Kc6Q0bUbeb3X76H5zYs9tjl2bMkXjbgfReuvjkjmoG0evLNCGifLo139X9146l40ZS2hn/Vx5tFFpH2g8dJWtLLLtWMWxnDgJRt/GLW9b7rmXjLYvcj/D2SPrtjTS626Tni4v/wCdq79ji8JmW3r1Znne8Ut/L/3fcFHepB/p3kVX8W/OM1nmrK9gD5fJyvUZGhvGKzmL86vRZPmLdJLh6ldZ1+acmfW55plMu9k36rq0aexq8055L32r6/pjEUToIMBWui6eJo8OLyRw8hoyK3il+EPpZCFaC17+o+kmNmVsrpXJ5QrioTNO6q7ODz+vZ+vjFtCR/72uveFdMXLtCV7slmR7+Y+mK/Y42of+1gsslPZai3Lai7bgSiG3md3uib0yvfhgg7c+XpleOnTTX3/vjF2OZNRniL10I7o5rvSRrq3gpeM1KjvwzH2To0Fq94+4thT1xvkWT1suGcmZVf7XSj/lUprXx8wVlwgjbdKzz8vvpXP1d+HmbUuvXurmHb/orud1wwX+eq3z0hm/J/aONWS5x4Ypzv1xMi9aWR076Ly4ezFHvavuA32jrlJJix892HXOwJrO1Q3h0HZhSZks6r3BuMzhK6ET/NUhNV2UrwUv/9F0rU2faHGbLHgtQeeeXdjJuctMw6D+QgOEH4XxNMuiecnzcl7ssvmzyMt/NF1tyGr7qB1+rIlJC3veErys82kvXf+o2010LECE+szbTY7j74BMLz7s4rvqI/O8Mr101Jg+zQ+WitNPZhVsAvbSjejGsflYbfKTLtODatqybh8vHbzMfZyhXqoPNGthaZyv8bRl7ByBA7Fd7ETvtamVP3o/0iY92V5+F91Af/e2pUsvFZPuzTmh6ou9te7PDNAzL10P0I0871hDzMjYMLWr86sRKd6ic+HuwbzS6a27t29UoifJI8b5ROBtb65qAepcPBXzmmryQXvlAzrGX/RHOrWtp0zn+q4tULvUV4N3jf/JUqGX/2i61h7Jx+msF8uWpNyLlp2VT3QxiX2ke8p4Iv5c/SktTop5yjv3kx52zIIXuxljzvDyH01X2yPZ7vYxPvEw9hnj9iMKK5rEomMuwFnih2azBxqV8woNx4pvCSNz8z++9WQqDxlbYW2szOrjlemlM+OEAw9kvUDfxwFIfdtDt0P3j+LBmech42+laT/GF0cDipNP2kMnPsIzXX+Ivh67yORVL4vlLDS03XHeMomHvgR2zAO0F3NCerBXPFmDVHaxIDtv2xdd/LfRI95Zf8+AbbblbfSiQ/yzOUHtd618intjlLmBwJGAnz10N+Rv/opndZ6RXNeYRKJoR8YG9XXNr1t0Knf1ize1fpMS7wxzKx2o+2bfMJkWS289J7jGObzZ3u5YN9lHxFeNEDz23hMVC9ivDe1beZs7GfX99wkqaIObCXwprHVyL//RdMVW4cPEye6N7SiXsjYh2uT0ipZ+Radm8vgqy/hNMa8wN+W0ck9878VuyQwv/9F0yZ6Mrat9mgrQ57m6QXKfqoBJ9LkuHLTefEDf/V5tmvp3toVX8bTzzBHakOnFRyq7oVcfr0wvXVcxmbluyUlbJOrTDesWti+kj7HGIk+9GVOztvTSidfmoGfiKQ840oHDzC/Ui6MsWuo6PM7Fg3380MvmBvoWO7F3xpGVLRaG28QYc+zl99I14tNtr79T4GnLXXrVhk8lf21OoJ9Q3gb4CGabl+6G64bPPc8YU7a3Oya9Y0MyXOufk24YdwfmVt1JvFB3w/9k4xwjpPtsY/3xpNY3nW+y8MgYKk0H5CD3ahAtE+q9DbnRf1EFmLSvL1QRe3rdq97Lv5euuzPXGissfyCvWbBS31IeCy7p2aILzz0IXuyWquLl30Pnap/aMLUVkzS/Zp+M/ZpGZa91faWLp+zvdPFAwq5gCSrDcS3jRmlr5+5uvso3ZRbh/UQXH099+uJSbldmh36Ljv7/o+q49VbLS1ebMNGd2wHc2UWlDWlPHiQm7eOhEw3zPYHd+nYHnjPt31Y0LFbD4zzbwY8MkzObtN08IHHeuOcAZZI7HU3aZIelXv4ZnTDrjl/lD7Xlhs0zvY7x+x9k1m2qNL6EzRG2meGlSyZK79A8k5hu/iyONdnFWNkcQ6Lxzq9eusq8bnKCuwPzrhBlTuo+2jdEPzzOMSTjeraxftXUns72ssmzCWZxsavomVDvZRDwOFe8SqBT/6D7vfXgtdlkB6MS1Nv9tmJ7WrPzRZZfx17+o+mSDcKEB5auc1IbmenoN/+q8pmEym6P8sHZ+lZFNpaULuzpyUl4qrxeNE04jtTSKykvdiarjb38R9MxebjapzVY92Bki0uneJpF/5YuGyecryW9FJDLD0B4vbioo5UpHi8+Pb1L9fHK9NL1dNMO9EmOLC31scS3QDekWzJoc86bWz/nNSCLMwsQzifnHXn966JLht386bUVO6jIeaaLV8Wk94zzb8VXzwW6TU448bv8uU2QXchgY2JE1gvh9HpB71CbdGR4+b10rYql/m50W235qxF2Ys+6ZK/tmQdYO9OckPsda8CXysO3wDl7lS/6aLLLSyf6tYCsxXlG+hfHpMpcYyPTede/TTrZu7e9Ew7CzTUPr9Vdgrb6xm3GOXbuGuvZ5uF1/Sohoz8SYMyt48oPdLznSdecMVN1J2N1js+FAZ+gYnL4N4PsaEPzwEVsb6K1vF4HS6Z4+Y+mQ7mwYReAV4SL9iUjb/6wWJVXQlV+nWSC88iqeWZp2WML+aRM9jJh4kRt7ZBN+LzYTZiqGy//0XSD7VNZnJJlcWkLJJd2YsF63ZSxCDJncDGx8gbnieg+UNwL1r/pS5syJQeZyCl8lVDLW+o/3focjXllT0nKZhZHcGDeXAxLdF4bK8E8iNcPjuZcsLnAfJbaR7GLrtJ/LZ6lwFxAGB7nqjf20H7s9taBfELbz25yB/5SB5Ev9cMBSTekFSbW72oZlrfUF2kPO09qtF1+L13NnNNb/R08lkJ9jnXVPhOgNtwcv6JNmzq5LSZztPjpi4SCmZdOvO555kZFWru2xuTm2JBe+uXm+ic69zo50t6S68bc6k0svm7dK92rfSPLGh7nWffusS77Jn0m20F9Vtf1KyNUvOS4YlTqnBJGg84AUD6NiKL3SYumdFTlTUI2iB2yRZoJw+CN5PNdtCEnxlSI7xvxs1DzdH/YhGjyc2yvBJvsd55U5W1Zfe/lP5qONn4ufNKrh8ogBhr/dpR8Hnzo/NCym5OCypgoWVxTyPf0q9UdrEx+iciL3ZJtXv4j6Z7JGG/7FLtzW9BeS45EakfRbe3Eov+vIvhNIvVr9YtavlemF5+izVEfr0wvXa07zYGqa+nXsgdscWrKfKe8LTqX7lzX7pyMPpUzbz/x0pWKvDm2UGVNklaXPeMcbJj/23XkE+Uzf5jsicI7cONqkxU7vfxeuqQqt+3a+F2SZ6Ya3kt0vXXJO35NRxuzZqQ3B21Bc9+jG5lnzAF6X/2qOyYHxgYYb86vosHh3qSTPayTBC/uw5irbnvnmRvL3jxw7BnnyDj/WBewTLYoB7DJB4t1T0UoTP8tRvHkY866p8PhyAAcTznQ89p+Qlfp4Gmh/OeZrBc+rkU+47c40y/yqLx8GN54RmLxY2e3DiNyerRZ9uTDztAp8IqifAi9x5vpsG2TXzSH0q3Yw4eq235DXyofXFaa/sGDjvU1+syt2shkLcVZZ7Fhia6XL14Xdj1e8rz8R9P17JGOWfvUdCrnwQKju+2xxJ/z4cOhUtQfL5lu8jH+nDfpM1kGttYyh9tB/Fv1ccmUHBedYSl65sLZnKE85rWEUa7jJp3o3bpFC2ZlbJk9WReLayrz0mU+xquS039qojzmqDL3KD08zsXDg+2kr+keTGh3zu2mOeKuxbLN3SY92738XjrTIfqt/u5tS3f9pLM7p+T8evyiG9q6/7O5MWlr3bvoqLPCbIzlfPS088zmWKt4N8eQYV7H0tnFoqapdPTmPBfuS3pyfsHc9Cp/s+6i8faN4XGe63z4WM82L67r5mBYJ5v8dyQxM4EBFo4qV/mPNEoDWHHAlC5KyKfcwM2VQ8ekQXVPY9Z8DM5Jp6xlZDmbPJKBnG7Hb+Ut3Yt/gsUS3Z58ZOvioHSNPwOjYKa0tUnBuKLf5IdW4VA601/H0oGiSZvpnr5R8FO6tIXSJ3dmc90ZrKVv1TZ70tiv6963keowa5+6/tQx05Q2asoZb2XcU1bxlHzlMT+0jgrO3GRBy/wumZl2sx0ae1frMyJTtrt057pDy3zZXvU4AKNNukEbmeuYb4vTkPl78+QmnWEpeTjDpU8gXxdtWY+JoXGeZdCByjqQ88Ck9CWz4a7FsnGzP+T6UMeTz9u5nT39fbMts6zN+mU61/gVBox/ZJa+qTS2TNpa9y66rNs1z0jmyFhzjSH0t5f0rM6vRr9Gp7JN3EXjwrzCCJntXMR9mY8y7WbfEM/QOM9yzX85dKzLltV1/RHKtTWdHEkR86vWJ7pe6aLx+MQSnY3D3H9yTz5B+QDBAmbHEXAeUrnKmDj5cRTb6SkoDzCpXJ3HBMl51fI6UnR0EF4R2KuQxG9/VO7iyTbsPgsrfnaAcdSKvWbDbWPJprGRf63rT13Pdb2scdA9GNMpeTUzOU8ywO/V46LDJguygfZk0qDvEOgHr2Sr9QEGAXZTB46z0Mbooc+U9tb9SYLso+Pz1F367IiiAYzvZBtttY9hITrakPHKlz3SWLYyi0VDG5dXdkrDg7MzGRtZlr1KYx75S1f53rDSJQzIHOqbzvq4ZEqWl44+ACa9wOv1dHxJ8lx0CPHqzrSMtS90gbeFyVxLpmS66EyA6JmjwIBAe/bmKGS6xrnkMSZ5oOV8L2kC6bPMCUnbLf5420R0Z5kTqIp0bY7fTOdpS1d/z/K8cwJ6CdaPZv2SQtXDRZdpqfPqPJPbYHNMIo8g+tGxsbr+3UhNcjfppNuFu+i8mLvnmVx3T99wj/Ms8yRjXRisruvm0AL6M5t4rTHWYgnGsfyXeDjwTmX5IQS/iEsDTOnJ9wSVz5NB+QpCpgP4ydk85SO3u7iO8IiWRsLB6i7SKotwBgRym3HMpP0E0Em1S+9qxz+p8hAeCDwwBDzjXDTD68wDgzGqGwjcaQQ845wKnGqsS+7quv44o4fnP9ltyflrEfTwEdI3LKWMJw0mrX+TaSGDcG33ObYntiY77TI8aTN38LDLvPpL4wUdkX0sAjzsjPatIyygv7V97gi5ISMQCATmCHjG+Z51Zq4pcgKBQOBSCHjGObadaqyvrutX2QllV5UjBO6gHbcX4uX7crwuf6rYePmeX+tIIL9+HWa0vRhn9t1ewUpejwcb0Bvhsgg8v8QuuXRewom+LNKhPRC4HAKr47xaZ9JxtsuZGZoDgUDgFgisjnPknnKsb63rV9L/LFfu1xy7IwlP51hUgXJ+doEZB7V1ctt7Y4W2e362I2ONx+tAm4yIT4MAZ4QjBAKBwNuNwNY4t3UmHjTf7n4QtXu7Edga59T+YmOdIwc4spN/JYdF3iBnli3o1xv0OKg4qiXIGSYPp7a3i9qVN8iD3CXHuNgRidMioDbb9aOs01oV0gOBQOBIBBzj/FbrzJG2hqxAIBDYh4BjnCP4YmP9sQy81nWbp2Z+TLbKr3Ic1J7j+lL5dg6XrWrS/P/z5Iiyda3Lfv2oohRWeYxIMfr40VmEQCAQCAQCgQsioDn9tuvMBa0P1YFAIOBF4JJjPX3lwGvobejkmHJG9zNzVk1W5bCylc1XEMpnflTG7u8vuiZfPVjjEW0KS/qsPOJAIBAIBAKBQCAQCAQCgbcDgXM6tOy+crTBvh/nRlDOKR/z/tHLIHqON/C9u/jKgRe0oAsEAoFAIBAIBAKBQOCeIvD4XHbLueRYAkcIekcPFs0QPY7w6pGGDjPHFPjYd4RAIBAIBAKBQCAQCAQCgbccgbM5tOAopxYn87O8g+qFlv/2xI/HXEGy+Rdx7M7GD8JciAVRIBAIBAKBQCAQCAQC9xuBsx05qGHCoR1xUmverfQpZW/pjvJAIBAIBAKBQCAQCAQCgfMjcBGH9vzVDI2BQCDgRSC/QeGoD8eDftbDZ/czerW8PTw1f6QDgUAgEAgEAoHbIHDWIwe3MXSJVwspX0KYBBZXXUNndScC4iYQeKAIaNzgyH6hCyf2G13PlPeD4sWwh2dRWBQEAoFAIBAIBAI7ELiqebQw2Tdf7RNaX46cRfXwiwZH0750wH+U+Iv73i6QR554fxEdXzWwXSTShA9uojd/R3TDVek3Id/VdmZ59uMz9FI3MOv+iE30X5sgxfy7Xj5H1j0fvKUbOaJBJ86HBe45P9zVb0Ti4yGAT6iZ7VZU13mzD3jqIxowMT0ejLANB+qDJWzM2FwPF63xEIvvU8nGWbtokB13bbzRPu2XSL6RnTwg0q9etIBRNsrTkXEn2qO1K+4DgUAgEAgE7hECWqT4sRYW/6aLz2PZPQsV/5iAH2WlvLXYw48sXT/VcnTPoo4BHzb5LnvEh41/ZxmkkccZ3YnNyhvRDS36i01K4zj9bnKVBp+v7Z5Y4WNdJAqOOd+w5F8EG744btg9wZd7Xau6KxkT/VkXdk70G73FKgcnHJQWIy/mrvpIhwujTIfdOPzYAIazNsz1Q6aLtq2f3Yufh45J3c99LxtcWC/Z5eEXjbvPo0eBsUOfnGDPvS4Inrb2KG+YpyejzYv7y/bPwD/wjz4QfeC+9QFzrvgyQHHWrBJ5sZo4oFZWx6Jz8YsOR2TiwCFHgUX0b5OptEte5p05ZianjiXTpTvLxNkozmfO4x9DFCwo14XxxXlU2hb+3zq6S/2sTPQTmV7dmQ6M2NmaOGbZhon+mkbl2D1zaJFFfk2b9eCwlHrnPLDcrI9o3BiZ3opn4lRZeR2P0LZ89f2507LbjXXPNi+/6Nx9Hj1Z7uwhKZcxRidjYi9PW6ee3JYm7qfjPPAIPKIPRB+IPjDtA49ZxRR4lWiv7FNG/vNK8Yf5tWKd36a9/JzP+70jr/1GrVdea8favUu3bGOXld3TyStpdRxexdb/qAG8rvOlSN7twvEBFSGz9xkxZBR8B3SLLf1Xtdoe8laD5FOvic0VwwjmrvpI9ghGlSmnSwoD+kGvr59O6VzyCNZzbv94dfV5U6D++40uOx5i2Rb/pcRzu7F4D4/xEt+R9qhNinQgEAgEAoHAPUTAHFoWPhasNpgTRvla8PLjuP6x4vixw0nwyruh9v316mZBv16xMWlTOb/+/gexqc8OKbflrKzyrE49fP/MvM9y7NKdaXnY+Fjyf6p0UMSOatFPRhU+kb0TR70qc2Fe6dqsjxejyoZzJHltXtrsHAo7OlxYd/gsy8vv7fMmN50vVht/revLfPHwQuD4gvXllGF/RPfpKI/xKr4L7VGZE8lAIBAIBAKB+4jAlRai7iLVVOZJc19uR/jlSLAz1QvsHLLD+XpEngkSD69wqcc/dbHwvkSW4hI8ujMxzuUfkolNn+jC6Xxf1+qPrUSPk4EzyY+titOo9LXKlJ1+BEZcB+wlYDPBrVtyf5Rc/h0wDsffSvNDO+zkeMDs3wSrnNf/XUdXZe4+MFgfqXwTpKeL0RuK06Skt27Lj3SPIrDix4jXp9HalzqCdU/CCL/qtjneah2SzREYHjjLLq3yeIPA+CLYA+7Nnf7u5Lkz7VEqEolAIBAIBAKBe43Alax/kmuwtrCvOTy34teCyOKGQ2dfPhiVh23fm2MiecjiWANHBFZ34jq6xVp2oZ6J32xi4cZp5KsEE2cxy8BR43UsTvSvutoADzRtoO4Ew9dil27Z8kL6cVJxOHCmcTiIJ0E0YMKu88whyYSjmHvrk8Q7MZrYfNRNxoeHpeSk6Z70VxmT3xTTT5ZwOcqMWs4o1jUv6Vvxq770uXq8JfnKx5ktOKXMm/ufVcaZcoLF6WYnT3qoukPtkeoSfwKBQCAQCATuNwKPnebbTqKTfEa2xs8PV37EyZhxLWcUeeLDIbk2UqVxTnBku7uRRpfjiW4t0OZQ8hq07LJm2u8Vf1vRpGzRvdb1lS52w77ThZNkr2kTjf78h4Tyi1OrNI6F2c2O8LDurAcZ7DZSZ3PmW/2TXWPR7QkFczFv1qdW4MSoZjkkLXxo3+Kk6R7800OObLJ+MnsAOET57YTUWO+RtMY/6fMIFy48EIFNeYAjvwrWNxkDKezkua/tYdWOOBAIBAKBQOCOIoBD2zsLaebabpCd9bT8Ot7Nr0URp5NXnPWr0d3yKqNwVt6TfBy8bljQbbS9HTt2p1jYORbQDaoHO5c4mD9IvjkBOFTk/UsXO6r/1YUThW2vdBFqfXX6pvRmZ2yiWzJwQnDmeW0OhvxAzHAsjnem23LuhzAfrI/VIcXi7WI0ITrgJte7ddLA63UlnnbhwWIxSA7fYOUhhV1/77UmcwjrjmG7+WV/b7yhAkeWHfwam1o1fZWHTvCyMMQj3T2nebg9THnEgUAgEAgEAoFAjcAVi5QWG/JwmNpgeT0nK9Hu5c8L3BPxT36pPyJPMnhNiowPWsPzvdk/KXborhfuCa9ukpMsGclpke7WCfhVNDhSXDhvKVAvJcrZRDIlw3YH0w/ldE/2pm6IFODFSS5BOjhXy24tn+Xi7CO28fmrxfaDGduy7h5eljeRsVUf5ErmEEbwHBhwuMpRlAW52DepV0uX67nUv1ryzfs9WNdC9/KrLXAoZ+NN+bRvclhrPZZWOf2YUB6K9vCI/5D2SJbEn0AgEAgEAoFAoEHgKt/b6+qm+J0nVXlbVt8P8WtB5JX4+1qcbUcR5yc5isrDwfDKY7e0t2OV7Jas1tlEz17dVl9zgNJ5QsnjSwdrTqjx9WIcKr6WYPxL9TbepFs6cUJwVI3PynFOOb6AIw0GYPpc9+lVbyG62ZVkB5t8nGmcjSXd3j6A+LY+R2CE3D2Buif9MGfMcPRTyPc4a6UP5qJzRLfFeohfdV3s81Vl7W1BlZWS9A2O1aDTgvWJEZ673B5Wr4gDgUAgEAgE7ikCj7PdODY4h21gZ4rFbOY4NYRufi2uOD3PsxNVi2HRNefUK4/vZrIj2QYclXoBTuVO3exEsfi2ASx4LWtywaR9DQuP4Wh0yYmWbn5UhiOaQk5jJw6DBZfu3B7sqvbsRBZ6cJS5XrSXyrDdyky/F/OR+rgxkj1Hh7Q7XQlt+8S3KqP/lF30ivbUSTfWC4a4+dVHVseb6s9DEu00C+JlTM6c/j08knOX22NW98gIBAKBQCAQuGcIaHFiVw+r2b1q/+sV/x2IH0gZDY4SxLP/RKU8Dz8OGHQ4bu01+S9VA/Im/9lIfHyeqvfvO0d0c5Sh/HtUpak3Mmt8eIU7+U9dlOsCnzbf/jsXu6qGJTomdLkdNnVnOhyNdK7XZOZ87Jr9R6eGBhsn//0r8262YaZz1Yf6cTW6uxgZjeg5SoF9s/8oZzQWr9GqDEeu9CloKz4cwkm/sbJzxdimq+5P1sfOPt5kB+1UsAIDBfCbjP8aG5UN8WR5RYfu71R71HWL9PS/7wQegUf0gdS/lp8AAAFdSURBVOgD96EPPMJIQt4xxJlgt+ZPXc91zb7nKjoWYnb32vOgLMir/Jl3aVeRneByVnHAHuTZLuMTpdnlnX1fdES3+MGDulAnAnJ7WLS7V9jCQl12Z2EmZHkkTebid209urNMnI4vdNnONtlrcm0HGLsJ7E6+kr3pCxPSu9mGiUt/vPURnQsj0eFkEqDHDnb02D3E8Z58ccJLKzrwoZ/SZzmrbbuEOLOkLxZkmwtr0Z1rvPGg8Ymuui/NxlENmGwb4rnL7VHXK9KBQCAQCAQC9w+B4tDeP9PD4kDAh4AcKR40PmwdYx93UB2NQLTH0YiGvEAgEAgEAoHHAUEg8AAQYKd2tmv+AOp9V6sY7XFXWybsCgQCgUDgniIQDu09bbgwewgBfoRoX6cYYgzikyAQ7XESWENoIBAIBAIPF4FwaB9u2z+kmnMmPMLdQSDa4+60RVgSCAQCgcBbgcD/Aep4w6Lk0PdcAAAAAElFTkSuQmCC)
$\displaystyle {{\mu_{\phi}}_{(0,0)}} \leftarrow 0.0004 \phi^{3} + 0.000473530700220886 \phi \rho^{2} - 0.00760399528440326 \phi \rho + 0.0205263968409311 \phi - 0.02 {\partial {\partial \phi}}$
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAArQAAAAyCAYAAACtSEg2AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAgAElEQVR4Ae2d7ZUcNRaG23MmADNEAGQAJoI1GSzrCGAzgON//PNZMgBHYCADdiMwdgaQgc1kMPs+al2NqlpVdVXd0z3tuTpHLZXqfunV1y2VuvvRzc3NJkIgcEoEfvjhh0+l/6nie0XyXyq+UPlbpRECgUDggSGgsf9YVX6eq82cQPhG5dfbbHwGAoFAIDBE4HJ4GVeBwEkQ+I+0/qXF6me0K/1Oyf8UP+I6QiAQCDw4BP6jeeDfVmvlf1L+jeJnVhZpIBAIBAI1Ahf1ReQDgRMh8EJ6X1W6P1ae3doIgUAg8DAR+FZOLG9tLPDQ+6nKPreCSAOBQCAQqBGIHdoajcifBAEtUuOjBf+UISxgEQKBQOBhIsDu7B8Ps+pR60AgEFiDwKM4Q7sGtuC5CwTk2HLU4JniK+V/vAsdIXOLgPC1c8u/KH+Qc4mS87ni4OFE15yFvFL6111iL/kHr89d2huy+xBQ+/KA+0+lceRAQDzk/q66n2ye6eu1QX1sBMKhPTbiWZ8Gpe1AvlMRkzRnxtyLvpf/0HRjuCSfV4D/VlrOu41p6mvR8SoxnZWtyy2vezhAvyr+rvxJnVrpP/s2Uh2a7aNydsHBeSpciyadYVaKs2jtS/twTX/9b82s6791zX1zaskTvtC96212++mVCbVoF9tBNO76eGV20oHJ9/AoPFHkyMz3ssuwoJy6uOgy7WK9Mx042xeoKOL6V+kat4+LDgEWJMNssCIeNgd1sht1KprZcV7TLuUliz5MX93pR0u8h75f4XHSeVt2dPV3cKhsN1h22rKima2f6FxzQtbL+WcLV8o0v9wnmV461zyTbVwckxgm2u6xkfma8yv36iD5s3S6b+NsEnfR9GDunmews9JvZu/0DbtRp+I72Div5a7Nh0O7Frk9+NQJ3oidb/H/hhilDCbKvlJ+0akVjYv/0HTYOg7S8afK3ir9enyvdS06HCGbZFokG91n8FNHFrDFxbMpZM9C6XVhPKXGy39ourE9kt9sH5UzgYJzq79xdhFn7DdF+iZtZg7tRnlbTL9WPvVh0VCOritFeJDLPfr5tdISdN0j09UOkumqD0aI1ivTS8fi8ZPkfmWVzPbwxoExnRxLpS46ZIjWpTvToru0Ty7DAWRRqtvHRZf5sRUZ9AOzn2t2xxZ3SUWzOM7RsxQkJ2EmOvraoB8t8R76vvS726Sl28vvoRNNT393taVHL/USnWv8Zjowo9+lzQmlzDl84Ze5Pc09XjrxpCD6xXlGNO6xhlDRu8fG1ortZ7Zlcf2bo9O9xX4lGhfmuS7uukuuq2/Uda7z4j/IOK9l7pO/2Ic5ePsRUAf4VlyPldYLDRM11/UTalO4l//QdC1jpIMFuze8qxkkAyz+VmSis2ALFztdq4LkPVUE6+6Q+c6+jVSP2fbRfZyttLtuqcBiofyvrq1/giFP4TixFmznr94V5B4T+0eKjxQ/U8QZsrY0XlKXTPF2jRXRL9bHK9NLlysFZmOHkoc26o4TaMFF16M707IgjsM3Kijt46WrhJhDbG3NrceKrQegiq1kB+O8lHZkZDOLLX2IduWNAV8Ko+zoQXq7+uLYQC+/lw75ol3s79mOxbbs0SuZrvErupeKHDcqb9qUZ4PiD8V6rfPSiS0FzzzjGmtIy3VfHENb1bef4pudX41yji7r9qw1XsxR6667aBf7htVjIt17nE/IbRYLr9l1/aLJFYUbOppi7WQdChV2Mlu7jq9VTmOxaMwFL/+h6QY2ZWxYsImuIB52/gZ1Vxn87xXrhRI6yn9RXBsei5G4Jnixm5Lt5T80XbHH0T70t1YY7MaKgPaiLYgp5DazyzWpV6YXH2zw1scr00uHbvrrn8Jl3N9wBlmszAnz0vXoZre07AwrPxW8dBvZy8MLc9/gaJDKcaAWdYlmZ5xPGTVVLhlghtPDzhm7wtjDQwJzxSlCT5u07PPye+lc/V24edvSq5e6eccvuut53XCBv17rvHTG70m9Yw1Z7rFhinN/HMyLdq9OHXRe3L2Yo95V946+UVep5MW/9zgvwvwZ5tjxPFu4L0ouMmME2HFpDcYxXe81naA1KZsu7s8FL/+h6cY2PVOHHix4Y4LGNQtTveNjJAzq57r3HVF5Fs1TnpfzYmf2j1Mv/6Hpajtm20c42w5s4VEZT/YvSoEyKmO3ll3X0m7KswAR6l2WbYnjs0OmF5+NZLrqI/O8Mr101Bhs+B1lFrhWsAnYS9ejG8eGL0tx5tz0YANtWbePlw5e5j52RKfqA81cmBrnczzje+yYgQOpxW/3sGksv/e6p01asr38Ljrh4O3v3rZ06aVi0r04J4jG+mJrrbNdvSdeuhagC2XesYaYnrFhamfnVyNSukTnwt2DeaXTW3dv36hED7KHGOcDgfteXNYCBBpPxW+UDn7QXteAjvHl1UHN96Hlc33nFqhVVZZcG+Rz/FdTN738h6Yb2yP5OJ31YjkmKdeiZWflmSKTGDs8StIT8ffKpwVTKU+fxJMH2XL2baQ6uNvHABcPY58xvnS+mbkgvdIS7c4DjcrSq1nR8FvCyOQM7Wzb6v6OTJXt2w479fHK9NKpbimIngeyVqDvb3Q/1d9DJ5queouec844NzxkcHSH9mPHCQe3OD1eOvERnigy/9VjF5k7XzSDmDCibY7zLeXtp3hod7BjHqC9mBPSg73SwRqkeycLsqWrTcaGevm9dGP5XIt3p79nusW23Edv1t0avzwQcbu1njE3EDhCgnNMfpYOgjqIZ3ae0X3XmESmaF1jyPSL3jW/LtHp/up+Jd4dzCv7vHVf7Bsm01LprecE1ziHN9vbHOsm+xDpxUgIHvv7URmXLGBMaB98EPB0MurLGbRDBxu0TOBTYa6Te/kPTVdsFT5MnExWtqNc7o0zosHp5YwmCxUPQ/bLBeDLgxOy7lvwYjdlt5f/0HTJnoypq31GFaBNiM0guTi7TOTMEThofzQI6bv8DNiPijhWRNqZyXcnqHxOphefHbm5oFUfr0wv3ZTuDXXTTfr30gPCmK5bt3SxUNjDBfUG752HCC+deG0OYgctjV2ltDsOLY7zIKise5yLBzv5oleaH5RHBjux9zF0t8moEl5+L91IfLps9XdueNpylV613dz4RTcPVK05/nNuKphtXrot15bPPc8YE/Yq3xyTuucdQ/Avzq+S56Hrxt2BuVV3kE7U3fC/s3GOEdJ9tLF+Maj1diJku7oEGUOl6Qi/l8KJjGh5ajrbIPupJ9/A/Eb56xNVxJ5e16r38q+lYwGyxXPSRtH8yk2lLISkLLKpbymPM0yejn6OwYvdVN28/GvoXO1TG6b2YPJ9qnQw9kc0fBEDR5WJ/5UijurAudE1T+xl3Chv7dzczdf9RZm1DY18Ex/JXaxPQ5YVNWXazSpdoqP/s/Oz9FbLS1epTrvf5Vo6aAdwZ9OBNqT+nOkdt88inXiY7wk4K+Nxzpn2lxXNRvnucS4e7OCYRJobUKbAA9JjlTUffhLF/f5Y6g9L1nv5d+iEWbO/q7yrLRcMbOldGr9pU6huU+VZY22OsE0RL10yUTK65pmqXpNjTTIXx0aW451fvXSVec3sAHfZuYR5U4gKB3WXnK6+IfrucY4h4jvqWL8c1Z7O9mJUZhPM5GJX0Z/tLq6Ax7liB4rBxk6EklWBHQ2eOFuhtfttdFc5Y+eLrLxOvfyHpks2qF48sDSdk9rITEe/+aQqZxKqd6vA2fpWRdaXlUzsaclJeOp+vWiacCYFHLNW8GLX4qXMy39ouo3q5GqfhuFgZItL4/awSHpw1GyccL6W/FRALs4yrxcndejeQKZ4vPi09E7VxyvTS9fSvVFd6JO8sp/qY4lvgq5Lt2TQ5pw3t37OOGMRYQHC+eSVLjtKLrpk2Paj1VbsoCLniSJyya8Z53yrvZ4LdJmccNLHfOwTZBcy2JjokcVu8dsJvV1t0pDh5ffSjVVM9XejW2pLHiamgmdd2gi7wfjVNX2OyBrAF03xLXDOXudIH012eelEPxeQNTnPSMfkmNQ919jIdN71b5FO9q5t74SD7NnBvAXQXN1Fv9Q39hnnmLNqrGebu9f1SwNAAox57LjyBR3vedI5Z8xU3ctU9ccR5fftmBz+ofzcIr2qDshUhLc10VpZq4MlfV7+Q9OhXDLZBXisdNK+ZOT2g8WKV0JzGDLBeWRVYnez0mEL+eCmypkwcaKWdsjGfGfZRqpnT/sM6qyLsriMb0gu7bRR+nZ0j0WQOYPIxMobnCulXyhtBevfyFqUKRpkIqfwVUKtbKr/NOsjea629dJV9pSseFkcwYF5czJM0a3QzYN4/eC4kQywY3OB+Sy1j1IXXaV/buzS1wjd41zysYf2Y7e3DpQTxv1sW9rxSR1EPtUPOyRtSStMrN/VMqxsqi9uvPxeulp5zi/1d/CYCvU5VqtLTWtlpX6yc3H8SkDa1KFOyg/maJXRFwm1TC+de57Zqkhz19KYXBwbspl+ubj+ic49D4ONImYaxuQtWFnCSHRuzE0AqfiadVe56Qb3qbB6nCNQOlaPdfEO+owZqPLZdf3SCJVOOa4YZZ2TBt0BQGVUHEX89uTSTgx07JCVzqzrgwXJ5ZvyXU6MKRffz4os1DzdH2xCNPk5tVeCo+JyKH78QDGm8/Ifmo42/lL4pFcPlVEMNNqcch58WOSgLefhVMbgZHFNIV/Tr2Z3sDL5KRIvdlO2efkPSfdExnjbp9id24L2mnIkUjuKbmknFv3vi+DbzBVZ8dfyvTK9+BRtjvp4ZXrpat02B5Z+LXvAlvqX+U75JTqXbslhXE3NyYxF5m2caxcdduYwpd/uW13WjHOwYf4fryPPVI7NJtt03Zd0CpPUv2Uk9+eCl99Ll3Tltp0bv1PyzFbDe4quVT/v+DUd45Q1I705GN8YXbfoeuaZjfCZHWsZv8UxJLvAeHF+FQ0O9yKd9LJOEry4d2O+VPcZ3ckwfVjfWDPOkXH0sV47tDgY9aKzESBUhEgjEZ4rWkNwnw5H2StFXifwpMPZqI8VC53KUlAZTwuDSUtl9rT2TveafFvu2885Ht3jnN/Ov37ccs/nxMtkC/9d/QMGjp/VuTYGB7o10dc05L38B6UTHgw84iConL8hZHIqi7iu6UcsohboWzXvS13z8JAelIzoHqVe7KZM9vIfkq6nfWq7WSAILWeUchyP1uJjfNautOfOmBfvuO17ZHrxQaYFs2uqPl6ZXrqkV3VnLmQhG2PAglrOojrpXLoli10W4tQmAmMwtZ2HLlVk+8HrUmwYB+Yo9FmbrxnnO/1B8sCOeFebCON6rLl2tcmMYC+/l85ULfV3b1v26HXNCWpX+j5z/Sf0PwxWSp+kD5S29tLBr+CeZyR3cUyKxjuGcO6s3ydD+BB/a/3z0iHCi7sLcwQSZNdi3UXm7RtrxjlmHH2sX6BVlaeT2aRCkQWbnHFCcWzrXTboOaP1tSKOyTul0LE7+lQp90vQNTo431UaW3kcXPjqb8Gb81x464yThzOwLaexFjWZFy8LEIPx4CHLfq+0yFcebP6lmA7Ho5QyxRvF9GRGGUHX2ObhPyhdUt7+wHZiHahH/cqAhT497SllAGN/fb/mPXm+A+NzaKNW+9QYW9td14VVnjlgMCaFD30XPr74YHw8BDJBlqBrzqQT6ocdrl0yxe/qwwiswmx9vDK9dOgVLXMj/Zr+kHCwVGUFow66nnqDLfOd1VuXySbmVn5VxHZZvHQb8TCf4wiXOVR55A/mKF13jfMsA6yuFFPIZWAHTiyc9zLINlebUB/FO5+3K5Cs3a+rspKVLa629NYvC3aNX9HS1uMHy1Zbe+lQ75pnVB/XmESggntsbMkHn+BvbTC4Mbpo0nXg7sV846276Fx9Q/XoGufUW7Kp79HH+qObmxuU40mzaDFomWxeK1LObieL0seKOJ7lVb7y0NO5AGWjtLzqVx4nDCe1dl7T08CojKcbzquWiUz5G5VxdMEmYl3eBpW7eESHDavPwoqfyZwFodTh1or9cpJJYyP/WpGd6S8VX6i84KDrja55gGBhGTh/uvbyH5QOmyzIBtqTDkvfIdAPXqs89RGlPNBgN3X4SpG6YQ99ZlBPlR08SAdOFztXpc/2KBGfF7t72Uayf7Z9DAvR0YaMlW+UT2PZ7lmqctqYSd8CPLzBGIyNLIuJl8A88l6ReeCagjp0yHS1g8nONizVxyVTsrx09AEwaQXeuqTdKKUuOoR4dWdaxtpzRfC2gJM7bh8XnQkQP3MUGBBoz9YchUzXOJc8xiQODed77YGe/FHmBOnZK3jbRHRHmROojHQtjt9M52lLV3/P8rxzAnoJ1o92+iU3VQ8XXaalzrPzjOS5x1qW2Ts2vPPrIp1sdeEuOi/mvXX39A33OM943slYFwaz67o5tIDOb5GV1wAYNRdEi2OZXiUoT2U5q4XzlQaY0sEPY+uan5Epv4KQ6QB+cDZP5chtLq49PKKlkXCwmou07kU4AgK5zZ4q5WHpaEH6Zjv+0QwJRYHAA0DAM85F073OPADoooqBwNkg4BnnVOauxrrkzq7rFxlJPP/B03wun0ugh4+Aw4Izy5MGk9Y/KLSgcpzca7vOqT2xjYrTLsPVuHAFD7vM7AxGOC0CPOz09q1DWEx/G/e5Q8gNGYFAILCLgGecr1lndjVFSSAQCJwKAc84x7a7Guuz6/pldkJxOAfn5JbQEh9nZ3ntyOvyz5UaC+VjRwL59eswo22lOLOPWzdmylo82IDeCKdFgPOzR98ll85TONGnRTq0BwKnQ2B2nGs8MqczH6fjbKczMzQHAoHAHgjMjnPk3uVYX1rXL6X/Sa7cHzl1JxKezrEoLednJ5iZzHAw6zC+tnvQNs/PqryHx+tAm95I7wYBzghHCAQCgQ8bgaVxbutMPGh+2P0gavdhI7A0zqn9ycb6hZTjyPLrA1POIgZOBvGxBf12kmB7AwcVR7UE8VGGztYualNeJw9ypxzjYkdk7hYBtdmqL2XdrVUhPRAIBA6JgGOc77XOHNLWkBUIBALrEHCMcwSfbKxz5ACncp+nZr5MNvuFH93nm74tx/WFdHPWIjmwoiHPv9wkR1QpTvBzpfaNRl2mv+ad5IEgB/TxpbMIgUAgEAgEAidEQHP4vuvMCa0P1YFAIOBF4JRj/cJr5BSdjJ91Zis++y3bUiRedu/4MwXO4vLzYPw0EL95ZgGn9Fvd45ttKTh4jJQvhP1mF5EGAoFAIBAIBAKBQCAQCHyYCKSf7TpG1eSIsqvK0YZ6t9WlGodW0e2cipadXX7vLn7lwIVwEAUCgUAgEAgEAoFAIHC+COy9Q+utupxLjjXwLyrsurqD6HGEe49E8Bu0gz8jcCsMwkAgEAgEAoFAIBAIBAKBs0LgaA4tqMg5xcnkLw7ZQfUG/u2J81euIFr+8pHd2fhCmAuxIAoEAoFAIBAIBAKBQOC8ETiqQwtUcjS7jhyI3ntG11riF/H07ugab6SBQCAQCAQCgUAgEAgEAmeGwNHO0J4ZLmFuIPBgEdADIW9QOOrD8SD+AbD5M3o1QGt4av7IBwKBQCAQCAQC+yBw9B3afYxt8Woh5XdwB4HFVZHFOEIgEAh0IKBxgyP7XBEnlrcj/Czfr0onwxqeSWFxIxAIBAKBQCAQWIHAYIdWCxNfpiLwbxCfKfJzWu6zqB5+0eBo2rED/lHiPdcq39kFcsr7W/zsKBm/nc/9QvyDs7c9uiVvU+nnkvBKZaaH+9TFvnyGXq7BrHnkQeU/6b4F/q73G5UNbLSbKre2sKKBbgpFg06cDwtcc364qd+IdJ+HAM4ym+12a6My07vYB0S7WB/R9GKEbThQO+1XjMwZyXbT1rzi+1ax9yhLLeIgedngxrql0MMvmp7xlvqTeGx8JrW65if1+MtDflZvEFTWzTMQoAvJuBftMbYrrgOBQCAQCATOB4Hi0GpReSOzXyhNP4+llIWKMn5qa9Gp9fCLhsX1J6Xl57SUZ1FnwURPccSUd9kjOv48AecQe7ET+6nHtdISdN2jG1qcKhztZJNSrj9XiqO/UYo+nNfiFCrP7+VC97Xy5WfGMi31oe7pn7OU4oz9TxHHreCr/KJu8aQgWuQV/RTqGv04v0V/Iq4+dA/M+LOLgYOiay/m1jdm6yN5Lowy3UvZxMMNDzlg85HKB22osk0PLfStIBm028Bpa9HdZZn0u7CessHDLxp3n0eP6BmLfKnyE+UL9srTjjw40lfLA52uN7ru5oGvDshQPGl71PZEPhAIBAKBQOD8ELjAZC0mLGK8pi9OkPIsaFzXu3C63A0d/Cx+YweMhQxdOGIpdMiDHscM5+eR4meKOKFlMd5KTJ8u3ZnenMLiYKucRb04nsqDGTtL5U8fdG309a6pijc4a1eiLX8DqzyOwR+KY3w9ujfiRz9O0TjwxxRj/YVGfDw87IQsz9sHvPVxYSTd14o8BNA3Xu0YVxX00FZs4+y7ccExr1WHezXeqrrzoMOXKgfjJ19TxnGEcVjDM5Zx0vYYGxPXgUAgEAgEAueHQHJoZTY7dYOdl1yV10qfakHDmZsLXn4WxD8b8nAE63OvXnlzNo3vuXTLNhxUdggHr6RVzg5y2VnWffBikS+Lv+6XvMrrgMzaGbZ7yCj4duiGn53i2h7KZoPkU6+BzRVDD+au+kh2D0aVKXeXFQb0g1Zfvzulu5J7sN7l9o9XV583BcLmZ8XBA6fdU/pe8cvqOmXX8NQy7kl71CZFPhAIBAKBQOAMEbjMNrPwDRy4XG5OGPfL7m2+VydefhxXXttPOX6Ps1CvvNqGpbxXNws6O4ZTNiY9uo+8j2qlKrPd2rLrqjKrEw7BONjO1BPdQJ5LdxbCw8Z3kv+7UnY3zV52oov+TGvJM9Gxg91yWlyYi9ddH9G6MDLjjpTS/8pO+ZF0jtW4sB4zVddefm+fL6KFDbvHXyhaf3qtMsb+p4o2Hyh7G9bw3HKn+eDU7VGZE9lAIBAIBAKBc0TgUouROShz9l9N3ezhFy07U63AzuFG9zk+0G2PeNIrXIn4WJGFlzO0g104XS/qFh8B5/Iv0WPTM0WcTnZDZ79sJXqcjHSsQfnycKA8zrFupXO+pHXAXgI2E9y6JfM3onhwov9WnqMb2Pl7Llf2NqiMowZNR1f33JiLtqc+twYoJ94mRgOiO7iQ3rot2WlHC1hNHU/h/p0E6XZj3TKgh1+03j6fVImehyP6fnngUZ43CIwvwo5Du5Ln3rTHtlrxGQgEAoFAIHDuCFyoAle5EtczlZlbhPfi14LI4oZDZ18K6ZWHbZz7+1ERGcQ3yuM8zYaGbuitrvxcEQ4PclngcWhxHgdBZez44SxCgxPNudhxsB2ucTl1J5hOS726cVjMecaZps4DR17XG9kHvjiiOw4J9xV6MffWJwmXXg9GifbQH9KNE88vOqS2VB6Hnx1B8KKfgM0xQy/WY9v24ld96XP1eEvyVY4zu1FanNl8bW8OuHzDh4WVPPetPaw6kQYCgUAgEAicMQIXTtttJ9FJvkM2x8+XoNht7HntWOSJjx23a9OoPE4bi3BzN9LocjrQLV5zKHHAzFE0ll+UeVnRpHJds6uM04tz+UoRJ2ns+PJFrY3Ki5OtPI6F2c2u2Brd6EEGu43UGUeFM8pj/Th04/qItCsUzMW1WJ9asnR7MKpZDpKXXtp3ozQ5aUrBH5wos36CY3vfQo31Gtvm+Ad9HuHCgh1YsLGHSorrYH2TMZDCSp5zbQ+rdqSBQCAQCAQC9xSBS9nVOttp5tpukJ31tPI6Xc2vRRGnE2eufjW6Wl5lFM7KU8n9VLG5Kzmh20S0eNidYuF/opicIiO2VDJxzHEw2c0tPztFmeInKufniXBkcThe54jzWeur87qVwo5uycEWfkbJdtRw7JGF04DjzT88oRe6Jee+C/Ms11sfqb8N4m1idEtxmFyuN04adloAo9ppo62gmQyix5nj59XMqZukrW5wpvltdV1nu7CuGXN+Nb9sao03xIIJfWXKZh6UrN2gJ3TxSLY5zXu1x1Z1fAYCgUAgEAgEAkMELrXQsJBR2lqwrazlZCVJa/nzAneldPBN/R55ouU1KTK+SMbsfpj9gzuiZ3Gd0309YBhesLhvJAPHlHTsBHDkACeJyKv5FESHTHM+rcx2B3HqrR0WdW8lplfmtXOALTgd7NbyU0o49Nj2WOlk+yFL9013Cy8rG8iAR6yT9clyuzCC54ABh4ujKHN4Yt+gXmP9mX+qf43JF6/XYF0LXcsvvqk+T/smh7XWY3nxmcNfHopU1s0jeQdpD7Mr0kAgEAgEAoFAoEbgMl/Y6+r6Hvmr6v74Xn3dxa8FkZ1EfjO27Mwqb44iDoZXHrulrR2rZLdkjp3NjcrW6rb6mgOUzhNKXtmJNYKOFIcq7aRmnql6m8ikWzpxKHBUd5w1leEc40iDAZh+qev0qld5C+hl95py6HE2pnQnLPN9JbNhXJ9DYDSrcOYmdU/6oVEdwQxHP4V8jbNW+mC+dYxkX6y7+FXXyT5fVZY3Bq1A3+DICDotWJ/o4bnP7WH1ijQQCAQCgUDgTBEwhxbHxnYL66qwM8VituM41UTKu/klC6cHJ4uFsg4sunbO0yuP380cy0Emjkq9AFO2cepmJwr94wAW7GSaXDCpnVGjx8kmGB16qdtLxfIPTCrDwcLOevfPpVu82EGcOlKBbGyzhwNd3gaV869P3K+dOS/mG/F56+PG6Na6g+XS7nQlbdwnaA/6T9lFr2jvOuvGesIQN7/qNzvedD/1pZae3M7gxq5/CSpPbxRKQZWZ4hHJfW6PqgaRDQQCgUAgEDhHBC4wWosQjuT7vBileiiPU/QvxfQFIAopU7xRLDtfudzLzy4NizFy+NvUElXGF5dwgDZKXfJEmvjhsSBe+yes2llDplc3Dg7OXnHwld/BQjQ40umb4ZVuHM65TNIAAAHkSURBVD1oS13yPXSPd5LBAbqyi6y8VzdiqR9nddFXgq55tcw3+W0nudyrMvCM+byYI8ZVH9H1YIRcwsfbpLwdyJfNZI6WflsfieAhyna4wZ7+Xt9vKriLQul1YS26o4w31ZF2GmAh3TjCOP2Dv2bWtYVennvbHlahSAOBQCAQCATOF4FHNzc3yXoWT2Vw4q4V3yl+qfhC5cXh0vVG17y2xeEbL4CL/JkXZ6gV2Akuu5XKL8pDiOiQx+JK4FUojuPO74uKDrtdukWHXLDABgJyW1iwe1U7zsjni19ld1bXKWR55E3m5O/aenQjSHQ4Hc8Va2d5Ti47wNiI3QQc6NeSk35hQqkLcxhFaw7/bH1E58JIdDiZBOiRSb/DAcU5xwEswUsrOvChn9L2nNVGJrJ5ECJ/siD9LqxFd6zxxsPYM8W6L+2Moxow2dbFI/p72x51vSIfCAQCgUAgcH4IFIf2/EwPiwMBHwJypJITr3TgGPu4g+rQCER7HBrRkBcIBAKBQCBwERAEAg8AAXYGd3bNH0C972sVoz3ua8uEXYFAIBAInCkC4dCeacOF2V0IlPOzXVxBfFcIRHvcFbIhNxAIBAKBB4pAOLQPtOEfWLU5Ex7h/iAQ7XF/2iIsCQQCgUDgg0Dg/5cbX+BvwiqcAAAAAElFTkSuQmCC)
$\displaystyle {\mu_{\phi}}_{(0,0)} \leftarrow 0.0004 \phi^{3} + 0.000473530700220886 \phi \rho^{2} - 0.00760399528440326 \phi \rho + 0.0205263968409311 \phi - 0.02 {\partial {\partial \phi}}$
3 2
μ_φ_C := 0.0004⋅φ + 0.000473530700220886⋅φ⋅ρ - 0.00760399528440326⋅φ⋅ρ + 0.0205263968409311⋅φ - 0.02⋅D(D(phi))
%% Cell type:markdown id: tags:
#### Pressure tensor and force computation
For the pressure tensor a trick for enhancing numerical stability is used: the bulk component is not stored directly in the pressure tensor field, but the related quantity called `pbs` is stored in a separate field.
$ pbs = \sqrt{|ρ c_s^2 - p_{bulk} |} $
The force is then calculated as $ \nabla \cdot P_{if} + 2 (\nabla pbs) pbs$
In the following kernel the pressure tensor field is filled with $P_{if}$ and the pbs field with above expression.
%% Cell type:code id: tags:
``` python
# Bulk part
pressure_assignments = [
Assignment(pbs_field.center,
pressure_tensor_bulk_sqrt_term(free_energy_transformed, (ρ, φ), ρ)),
]
# Interface part
P_if = pressure_tensor_from_free_energy(free_energy_transformed, (ρ, φ),
dim=dh.dim, include_bulk=False)
index_map = symmetric_tensor_linearization(dh.dim)
pressure_assignments += [
Assignment(pressure_tensor_field(index_1d), P_if[index_2d])
for index_2d, index_1d in index_map.items()
]
pressure_kernel = make_kernel(pressure_assignments)
# Force kernel
pressure_tensor_sym = sp.Matrix(dh.dim, dh.dim,
lambda i, j: pressure_tensor_field(index_map[i, j]
if i < j else index_map[j, i]))
force_term = force_from_pressure_tensor(pressure_tensor_sym,
functions=[ρ, φ],
pbs=pbs_field.center)
force_assignments = [
Assignment(force_field(i),
force_term[i] + external_force[i] * ρ_field.center / ρ_l)
for i in range(dh.dim)
]
force_kernel = make_kernel(force_assignments)
```
%% Cell type:markdown id: tags:
#### Lattice Boltzmann schemes for time evolution of ρ and φ
- ρ is handled by a normal LB method (compressible, entropic equilibrium)
- stream and collide are splitted into separate kernels
- macroscopic values are computed after the stream, but inside the stream kernel
- velocity field stores the velocity which was not corrected for the forces yet
- the φ collision kernel corrects the velocity itself, because then the updated forces are used for the correction. When u is computed, the updated forces are not computed yet
- when ρ and φ are updated, they are clipped to a valid region, this clipping should be only necessary during equilibration of the system
- exact difference method is used to couple the force into the ρ-LBM
%% Cell type:markdown id: tags:
The following cell handles the clipping of the order parameters:
%% Cell type:code id: tags:
``` python
if clipping:
def clip(ac, symbol, min_value, max_value):
"""Function to clip the value of a symbol which is on one of lhs of the assignments
in an assignment collection"""
assert symbol in ac.bound_symbols
for i in range(len(ac.subexpressions)):
a = ac.subexpressions[i]
if a.lhs == symbol:
new_assignment = Assignment(symbol, sp.Piecewise((max_value, a.rhs > max_value),
(min_value, a.rhs < min_value),
(a.rhs, True)))
ac.subexpressions[i] = new_assignment
break
# TODO: how can this 'densgin' be derived automatically?
tred = reduced_temperature
densgin = -67.098 \
+ 549.69 * tred \
- 1850.6 * tred * tred \
+ 3281 * tred * tred * tred \
- 3237.3 * tred * tred * tred * tred \
+ 1687.6 * tred * tred * tred * tred * tred \
- 361.51 * tred * tred * tred * tred * tred * tred
ρ_clip_min, ρ_clip_max = densgin * 0.5, 1.2 * ρ_l
φ_clip_min, φ_clip_max = -χ * 1.5, χ * 1.5
```
%% Cell type:markdown id: tags:
Next, the collide and stream kernels for the ρ lattice Boltzmann are created
%% Cell type:code id: tags:
``` python
force_model = EDM(force_field.center_vector)
lbm_config = LBMConfig(stencil=stencil, method=Method.TRT_KBC_N2, compressible=True,
relaxation_rate=ρ_relaxation_rate, kernel_type='collide_only',
density_input=ρ_field, force_model=force_model, velocity_input=vel_field)
lbm_opt = LBMOptimisation(symbolic_field=pdf_src_rho,
symbolic_temporary_field=pdf_dst_rho)
config = ps.CreateKernelConfig(target=target, cpu_openmp=threads)
# Standard collision step, that does not compute ρ and u from pdfs, but reads
# them from fields - this is necessary because ρ may have been clipped before
# the velocity field is not force corrected, which is the correct for the EDM model
# but might be wrong for other force models
ρ_collide = create_lb_function(lbm_config=lbm_config,
lbm_optimisation=lbm_opt,
config=config)
lbm_config = LBMConfig(stencil=stencil, method=Method.TRT_KBC_N2, compressible=True,
relaxation_rate=ρ_relaxation_rate, kernel_type='stream_pull_only',
output={'density': ρ_field, 'velocity': vel_field})
# First the assignments are created, then the density is clipped
# then a kernel is created from the clipped assignments
ρ_stream_ur = create_lb_update_rule(lbm_config=lbm_config,
lbm_optimisation=lbm_opt,
config=config)
if clipping:
clip(ρ_stream_ur, sp.Symbol("rho"), ρ_clip_min, ρ_clip_max)
lbm_config = LBMConfig(update_rule=ρ_stream_ur)
ρ_stream = create_lb_function(lbm_config=lbm_config,
lbm_optimisation=lbm_opt,
config=config)
```
%% Cell type:markdown id: tags:
The φ lattice Boltzmann solve the Cahn-Hilliard equation.
%% Cell type:code id: tags:
``` python
φ_lb_method = cahn_hilliard_lb_method(stencil=stencil,
mu=μ_phi_field.center,
relaxation_rate=φ_relaxation_rate,
gamma=1)
corrected_vel = vel_field.center_vector + sp.Matrix(force_model.macroscopic_velocity_shift(ρ_field.center))
lbm_opt = LBMOptimisation(symbolic_field=pdf_src_phi,
symbolic_temporary_field=pdf_dst_phi)
config = ps.CreateKernelConfig(target=target, cpu_openmp=threads)
lbm_config = LBMConfig(lb_method=φ_lb_method, compressible=True,
kernel_type='collide_only',
density_input=φ_field, velocity_input=corrected_vel)
φ_collide = create_lb_function(lbm_config=lbm_config,
lbm_optimisation=lbm_opt,
config=config)
lbm_config = LBMConfig(lb_method=φ_lb_method, compressible=True,
kernel_type='stream_pull_only', output={'density': φ_field})
φ_stream_ur = create_lb_update_rule(lbm_config=lbm_config,
lbm_optimisation=lbm_opt,
config=config)
if clipping:
clip(φ_stream_ur, sp.Symbol("rho"), φ_clip_min, φ_clip_max)
lbm_config = LBMConfig(update_rule=φ_stream_ur)
φ_stream = create_lb_function(lbm_config=lbm_config,
lbm_optimisation=lbm_opt,
config=config)
```
%% Cell type:markdown id: tags:
#### Time loop
Now we can put all kernels together into a time loop function
%% Cell type:code id: tags:
``` python
op_sync = dh.synchronization_function([ρ_field.name, φ_field.name])
p_sync = dh.synchronization_function([pbs_field.name, pressure_tensor_field.name])
pdf_sync = dh.synchronization_function([pdf_src_phi.name, pdf_src_rho.name])
def time_loop(steps):
for t in range(steps):
op_sync()
dh.run_kernel(μ_kernel)
dh.run_kernel(pressure_kernel)
p_sync()
dh.run_kernel(force_kernel)
dh.run_kernel(ρ_collide)
dh.run_kernel(φ_collide)
pdf_sync()
dh.run_kernel(ρ_stream)
dh.run_kernel(φ_stream)
dh.swap(pdf_dst_phi.name, pdf_src_phi.name)
dh.swap(pdf_dst_rho.name, pdf_src_rho.name)
return dh.cpu_arrays[φ_field.name][1:-1, 1:-1]
```
%% Cell type:markdown id: tags:
## Part 3b: Compiling getter & setter kernels
%% Cell type:markdown id: tags:
The setter kernel computes ρ, φ from C and sets the pdfs to equilibrium using the values in the order parameter and velocity fields.
%% Cell type:code id: tags:
``` python
init_assignments = [
Assignment( φ_field.center, transform_backward_substitutions[φ].subs({ c_l1: c_field(0), c_l2: c_field(1)} )),
Assignment( ρ_field.center, transform_backward_substitutions[ρ].subs({ c_l1: c_field(0), c_l2: c_field(1)} )),
]
init_rho = pdf_initialization_assignments(ρ_collide.method,
density=ρ_field.center,
velocity=vel_field.center_vector,
pdfs=pdf_src_rho.center_vector)
init_rho = init_rho.new_without_subexpressions()
init_assignments += init_rho.all_assignments
init_phi = pdf_initialization_assignments(φ_collide.method,