diff --git a/experiment/wood-bilayer_PLOS/PolarPlotLocalEnergy.py b/experiment/wood-bilayer_PLOS/PolarPlotLocalEnergy.py
index 6379e849fdaf9814bb1bfe56dd0bfd2d2c8f34b5..cfe226377313c755ef192b8a24237b4104a454c9 100644
--- a/experiment/wood-bilayer_PLOS/PolarPlotLocalEnergy.py
+++ b/experiment/wood-bilayer_PLOS/PolarPlotLocalEnergy.py
@@ -54,6 +54,7 @@ show_plot = False
#--- Choose wether to perforate upper (passive) or lower (active) layer
dataset_numbers = [0, 1, 2, 3, 4, 5]
+#dataset_numbers = [5]
number=7
for dataset_number in dataset_numbers:
kappa=np.zeros(number)
@@ -105,23 +106,23 @@ for dataset_number in dataset_numbers:
kappamin_neg=r[imin+N_mid,jmin]
alphamin_neg=theta[imin+N_mid,jmin]
Emin_neg=E[imin+N_mid,jmin]
- kappa_neg[n]=kappamin_neg
+ kappa_neg[n]=-kappamin_neg
#
- fig, ax = plt.subplots(figsize=(6,6),subplot_kw=dict(projection='polar'))
- levs=np.geomspace(E.min(),E.max(),400)
- pcm=ax.contourf(theta, r, E, levs, norm=colors.PowerNorm(gamma=0.2), cmap='brg')
- ax.set_xticks([0,np.pi/2])
- ax.set_yticks([kappamin_pos,kappamin_neg])
- colorbarticks=np.linspace(E.min(),E.max(),6)
- cbar = plt.colorbar(pcm, extend='max', ticks=colorbarticks, pad=0.1)
- cbar.ax.tick_params(labelsize=6)
- if (show_plot):
+ if show_plot:
+ fig, ax = plt.subplots(figsize=(6,6),subplot_kw=dict(projection='polar'))
+ levs=np.geomspace(E.min(),E.max(),400)
+ pcm=ax.contourf(theta, r, E, levs, norm=colors.PowerNorm(gamma=0.2), cmap='brg')
+ ax.set_xticks([0,np.pi/2])
+ ax.set_yticks([kappamin_pos,kappamin_neg])
+ colorbarticks=np.linspace(E.min(),E.max(),6)
+ cbar = plt.colorbar(pcm, extend='max', ticks=colorbarticks, pad=0.1)
+ cbar.ax.tick_params(labelsize=6)
plt.show()
- # Save Figure as .pdf
- width = 5.79
- height = width / 1.618 # The golden ratio.
- fig.set_size_inches(width, height)
- fig.savefig('./experiment/wood-bilayer_PLOS/wood-bilayer_PLOS_dataset_' +str(dataset_number) + '_exp' +str(n) + '.pdf')
+ # Save Figure as .pdf
+ width = 5.79
+ height = width / 1.618 # The golden ratio.
+ fig.set_size_inches(width, height)
+ fig.savefig('./experiment/wood-bilayer_PLOS/wood-bilayer_PLOS_dataset_' +str(dataset_number) + '_exp' +str(n) + '.pdf')
f = open("./experiment/wood-bilayer_PLOS/results_" + str(dataset_number) + "/kappa_simulation.txt", "w")
diff --git a/experiment/wood-bilayer_PLOS/PolarPlotLocalEnergy_5_6.py b/experiment/wood-bilayer_PLOS/PolarPlotLocalEnergy_5_6.py
new file mode 100644
index 0000000000000000000000000000000000000000..05aa178229e225ccae67dbd267edd265db3e6384
--- /dev/null
+++ b/experiment/wood-bilayer_PLOS/PolarPlotLocalEnergy_5_6.py
@@ -0,0 +1,148 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Jul 6 13:17:28 2022
+
+@author: stefan
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.colors as colors
+from matplotlib.ticker import LogLocator
+import codecs
+import re
+import json
+
+def energy(kappa,alpha,Q,B) :
+ G=kappa*np.array([[np.cos(alpha)**2],[np.sin(alpha)**2],[np.sqrt(2)*np.cos(alpha)*np.sin(alpha)]])-B
+ return np.matmul(np.transpose(G),np.matmul(Q,G))[0,0]
+
+def xytokappaalpha(x,y):
+
+ if y>0:
+ return [np.sqrt(x**2+y**2), np.abs(np.arctan2(y,x))]
+ else:
+ return [-np.sqrt(x**2+y**2), np.abs(np.arctan2(y,x))]
+
+# Read effective quantites
+def ReadEffectiveQuantities(QFilePath, BFilePath):
+ # Read Output Matrices (effective quantities)
+ # From Cell-Problem output Files : ../outputs/Qmatrix.txt , ../outputs/Bmatrix.txt
+ # -- Read Matrix Qhom
+ X = []
+ # with codecs.open(path + '/outputs/QMatrix.txt', encoding='utf-8-sig') as f:
+ with codecs.open(QFilePath, encoding='utf-8-sig') as f:
+ for line in f:
+ s = line.split()
+ X.append([float(s[i]) for i in range(len(s))])
+ Q = np.array([[X[0][2], X[1][2], X[2][2]],
+ [X[3][2], X[4][2], X[5][2]],
+ [X[6][2], X[7][2], X[8][2]] ])
+
+ # -- Read Beff (as Vector)
+ X = []
+ # with codecs.open(path + '/outputs/BMatrix.txt', encoding='utf-8-sig') as f:
+ with codecs.open(BFilePath, encoding='utf-8-sig') as f:
+ for line in f:
+ s = line.split()
+ X.append([float(s[i]) for i in range(len(s))])
+ B = np.array([X[0][2], X[1][2], X[2][2]])
+ return Q, B
+
+# Number of experiments / folders
+show_plot = False
+
+#--- Select specific experiment [x, y] with date from results_x/y
+data=[5,6]
+#DataPath = './experiment/wood-bilayer_PLOS/results_' + str(data[0]) + '/' +str(data[1])
+DataPath = './results_' + str(data[0]) + '/' +str(data[1])
+QFilePath = DataPath + '/QMatrix.txt'
+BFilePath = DataPath + '/BMatrix.txt'
+ParameterPath = DataPath + '/parameter.txt'
+#
+# Read Thickness from parameter file (needed for energy scaling)
+with open(ParameterPath , 'r') as file:
+ parameterFile = file.read()
+thickness = float(re.findall(r'(?m)h = (\d?\d?\d?\.?\d+[Ee]?[+\-]?\d?\d?)',parameterFile)[0])
+energyscalingfactor = thickness**2
+# Read Q and B
+Q, B = ReadEffectiveQuantities(QFilePath,BFilePath)
+Q=0.5*(np.transpose(Q)+Q) # symmetrize
+B=np.transpose([B])
+#
+# Compute lokal and global minimizer
+kappa=0
+kappa_pos=0
+kappa_neg=0
+#
+N=500
+length=4
+r, theta = np.meshgrid(np.linspace(0,length,N),np.radians(np.linspace(0, 360, N)))
+E=np.zeros(np.shape(r))
+for i in range(0,N):
+ for j in range(0,N):
+ if theta[i,j]<np.pi:
+ E[i,j]=energy(r[i,j],theta[i,j],Q,B) * energyscalingfactor
+ else:
+ E[i,j]=energy(-r[i,j],theta[i,j],Q,B) * energyscalingfactor
+#
+# Compute Minimizer
+[imin,jmin]=np.unravel_index(E.argmin(),(N,N))
+kappamin=r[imin,jmin]
+alphamin=theta[imin,jmin]
+# Positiv curvature region
+N_mid=int(N/2)
+[imin,jmin]=np.unravel_index(E[:N_mid,:].argmin(),(N_mid,N))
+kappamin_pos=r[imin,jmin]
+alphamin_pos=theta[imin,jmin]
+Emin_pos=E[imin,jmin]
+# Negative curvature region
+[imin,jmin]=np.unravel_index(E[N_mid:,:].argmin(),(N_mid,N))
+kappamin_neg=r[imin+N_mid,jmin]
+alphamin_neg=theta[imin+N_mid,jmin]
+Emin_neg=E[imin+N_mid,jmin]
+#
+E=E/E.min()
+print(Emin_pos/Emin_neg)
+fig, ax = plt.subplots(figsize=(6,6),subplot_kw=dict(projection='polar'))
+levs=np.geomspace(1,E.max(),1000)
+pcm=ax.contourf(theta, r, E, levs, norm=colors.PowerNorm(gamma=0.4), cmap='brg')
+ax.set_xticks(np.array([.0,1/4,2/4,3/4,1,5/4,6/4,7/4])*np.pi)
+anglelabel=["0°","45°", "90°", "135°","180°","135°","90°","45°"]
+ax.set_xticklabels(anglelabel)
+ax.set_yticks([1,2,3,4])
+#ax.set_yticklabels(["1$m^{-1}$","2$m^{-1}$","3$m^{-1}$","4$m^{-1}$"])
+#
+ax.plot([alphamin_pos,alphamin_pos+np.pi], [kappamin_pos,kappamin_pos],
+ markerfacecolor='red',
+ markeredgecolor='black', # marker edgecolor
+ marker='s', # each marker will be rendered as a circle
+ markersize=5, # marker size
+ markeredgewidth=0.5, # marker edge width
+ linewidth=0,
+ zorder=3,
+ alpha=1, # Change opacity
+ label = r"$\kappa_{2,sim}(m^{-1})$")
+ax.plot(alphamin_neg, kappamin_neg,
+ markerfacecolor='blue',
+ markeredgecolor='black', # marker edgecolor
+ marker='D', # each marker will be rendered as a circle
+ markersize=5, # marker size
+ markeredgewidth=0.5, # marker edge width
+ linewidth=0, # line width
+ zorder=3,
+ alpha=1, # Change opacity
+ label = r"$\kappa_{1,sim}(m^{-1})$")
+colorbarticks=np.linspace(1,15,10)
+cbar = plt.colorbar(pcm, extend='max', ticks=colorbarticks, pad=0.1)
+#bounds = ['0','1/80','1/20','1/5','4/5']
+cbar.ax.tick_params(labelsize=8)
+#cbar.set_ticklabels(bounds)
+fig.legend(loc="upper left")
+if (show_plot):
+ plt.show()
+# Save Figure as .pdf
+width = 5.79
+height = width / 1.618 # The golden ratio.
+fig.set_size_inches(width, height)
+#fig.savefig('./experiment/wood-bilayer_PLOS/wood-bilayer_PLOS_dataset_' +str(data[0]) + '_exp' + str(data[1]) + '.pdf', dpi=300)
\ No newline at end of file
diff --git a/experiment/wood-bilayer_PLOS/auswertung.ipynb b/experiment/wood-bilayer_PLOS/auswertung.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..bbfc92325bb2f5efa10b8a6211d44690207541e9
--- /dev/null
+++ b/experiment/wood-bilayer_PLOS/auswertung.ipynb
@@ -0,0 +1,530 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 83,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "%%capture\n",
+ "#!/usr/bin/env python3\n",
+ "# -*- coding: utf-8 -*-\n",
+ "\"\"\"\n",
+ "Created on Wed Jul 6 13:17:28 2022\n",
+ "\n",
+ "@author: stefan\n",
+ "\"\"\"\n",
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt\n",
+ "import matplotlib.colors as colors\n",
+ "from matplotlib.ticker import LogLocator\n",
+ "import codecs\n",
+ "import re\n",
+ "import json\n",
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt\n",
+ "import math\n",
+ "import os\n",
+ "import subprocess\n",
+ "import fileinput\n",
+ "import re\n",
+ "import sys\n",
+ "import matplotlib as mpl\n",
+ "from mpl_toolkits.mplot3d import Axes3D\n",
+ "import matplotlib.cm as cm\n",
+ "import matplotlib.ticker as ticker\n",
+ "from matplotlib.ticker import MultipleLocator,FormatStrFormatter,MaxNLocator\n",
+ "import seaborn as sns\n",
+ "import matplotlib.colors as mcolors\n",
+ "\n",
+ "def energy(kappa,alpha,Q,B) :\n",
+ " G=kappa*np.array([[np.cos(alpha)**2],[np.sin(alpha)**2],[np.sqrt(2)*np.cos(alpha)*np.sin(alpha)]])-B\n",
+ " return np.matmul(np.transpose(G),np.matmul(Q,G))[0,0]\n",
+ "\n",
+ "def xytokappaalpha(x,y):\n",
+ " \n",
+ " if y>0:\n",
+ " return [np.sqrt(x**2+y**2), np.abs(np.arctan2(y,x))]\n",
+ " else:\n",
+ " return [-np.sqrt(x**2+y**2), np.abs(np.arctan2(y,x))]\n",
+ "\n",
+ "# Read effective quantites\n",
+ "def ReadEffectiveQuantities(QFilePath, BFilePath):\n",
+ " # Read Output Matrices (effective quantities)\n",
+ " # From Cell-Problem output Files : ../outputs/Qmatrix.txt , ../outputs/Bmatrix.txt\n",
+ " # -- Read Matrix Qhom\n",
+ " X = []\n",
+ " # with codecs.open(path + '/outputs/QMatrix.txt', encoding='utf-8-sig') as f:\n",
+ " with codecs.open(QFilePath, encoding='utf-8-sig') as f:\n",
+ " for line in f:\n",
+ " s = line.split()\n",
+ " X.append([float(s[i]) for i in range(len(s))])\n",
+ " Q = np.array([[X[0][2], X[1][2], X[2][2]],\n",
+ " [X[3][2], X[4][2], X[5][2]],\n",
+ " [X[6][2], X[7][2], X[8][2]] ])\n",
+ "\n",
+ " # -- Read Beff (as Vector)\n",
+ " X = []\n",
+ " # with codecs.open(path + '/outputs/BMatrix.txt', encoding='utf-8-sig') as f:\n",
+ " with codecs.open(BFilePath, encoding='utf-8-sig') as f:\n",
+ " for line in f:\n",
+ " s = line.split()\n",
+ " X.append([float(s[i]) for i in range(len(s))])\n",
+ " B = np.array([X[0][2], X[1][2], X[2][2]])\n",
+ " return Q, B"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Parameters from Simulation"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 84,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "materialFunctionParameter=[\n",
+ "[ # Dataset Ratio r = 0.12\n",
+ "[0.12, 0.0047, 17.32986047, 14.70179844, 0.0, 0.0 ],\n",
+ "[0.12, 0.0047, 17.32986047, 13.6246, 0.0, 0],\n",
+ "[0.12, 0.0047, 17.32986047, 12.42994508, 0.0, 0 ],\n",
+ "[0.12, 0.0047, 17.32986047, 11.69773413, 0.0, 0],\n",
+ "[0.12, 0.0047, 17.32986047, 11.14159987, 0.0, 0],\n",
+ "[0.12, 0.0047, 17.32986047, 9.500670278, 0.0, 0],\n",
+ "[0.12, 0.0047, 17.32986047, 9.005046347, 0.0, 0]\n",
+ "],\n",
+ "[ # Dataset Ratio r = 0.17\n",
+ "[0.17, 0.0049, 17.28772791 , 14.75453569, 0.0, 0 ],\n",
+ "[0.17, 0.0049, 17.28772791 , 13.71227639, 0.0, 0],\n",
+ "[0.17, 0.0049, 17.28772791 , 12.54975012, 0.0, 0 ],\n",
+ "[0.17, 0.0049, 17.28772791 , 11.83455959, 0.0, 0],\n",
+ "[0.17, 0.0049, 17.28772791 , 11.29089521, 0.0, 0 ],\n",
+ "[0.17, 0.0049, 17.28772791 , 9.620608917, 0.0, 0],\n",
+ "[0.17, 0.0049, 17.28772791 , 9.101671742, 0.0, 0 ]\n",
+ "],\n",
+ "[ # Dataset Ratio r = 0.22\n",
+ "[0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0 ],\n",
+ "[0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0 ],\n",
+ "[0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0 ],\n",
+ "[0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0 ],\n",
+ "[0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0 ],\n",
+ "[0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0 ],\n",
+ "[0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0 ]\n",
+ "],\n",
+ "[ # Dataset Ratio r = 0.34\n",
+ "[0.34, 0.0063, 17.14061081 , 14.98380876, 0.0, 0 ],\n",
+ "[0.34, 0.0063, 17.14061081 , 13.97154915, 0.0, 0],\n",
+ "[0.34, 0.0063, 17.14061081 , 12.77309253, 0.0, 0 ],\n",
+ "[0.34, 0.0063, 17.14061081 , 12.00959929, 0.0, 0],\n",
+ "[0.34, 0.0063, 17.14061081 , 11.42001731, 0.0, 0 ],\n",
+ "[0.34, 0.0063, 17.14061081 , 9.561447179, 0.0, 0],\n",
+ "[0.34, 0.0063, 17.14061081 , 8.964704969, 0.0, 0 ]\n",
+ "],\n",
+ "[ # Dataset Ratio r = 0.43\n",
+ "[0.43, 0.0073, 17.07559686 , 15.11316339, 0.0, 0 ],\n",
+ "[0.43, 0.0073, 17.07559686 , 14.17997082, 0.0, 0],\n",
+ "[0.43, 0.0073, 17.07559686 , 13.05739844, 0.0, 0 ],\n",
+ "[0.43, 0.0073, 17.07559686 , 12.32309209, 0.0, 0],\n",
+ "[0.43, 0.0073, 17.07559686 , 11.74608518, 0.0, 0 ],\n",
+ "[0.43, 0.0073, 17.07559686 , 9.812372466, 0.0, 0],\n",
+ "[0.43, 0.0073, 17.07559686 , 9.10519385 , 0.0, 0 ]\n",
+ "],\n",
+ "[ # Dataset Ratio r = 0.49\n",
+ "[0.49, 0.008, 17.01520754, 15.30614414, 0.0, 0 ],\n",
+ "[0.49, 0.008, 17.01520754, 14.49463867, 0.0, 0],\n",
+ "[0.49, 0.008, 17.01520754, 13.46629742, 0.0, 0 ],\n",
+ "[0.49, 0.008, 17.01520754, 12.78388234, 0.0, 0],\n",
+ "[0.49, 0.008, 17.01520754, 12.23057715, 0.0, 0 ],\n",
+ "[0.49, 0.008, 17.01520754, 10.21852839, 0.0, 0],\n",
+ "[0.49, 0.008, 17.01520754, 9.341730605, 0.0, 0 ]\n",
+ "]\n",
+ "]"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Load data specific simulation:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 85,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "#--- Select specific experiment [x, y] with date from results_x/y\n",
+ "def get_Q_B(index):\n",
+ " # results_index[0]/index[1]/...\n",
+ " #DataPath = './experiment/wood-bilayer_PLOS/results_' + str(data[0]) + '/' +str(data[1])\n",
+ " DataPath = './results_' + str(index[0]) + '/' +str(index[1])\n",
+ " QFilePath = DataPath + '/QMatrix.txt'\n",
+ " BFilePath = DataPath + '/BMatrix.txt'\n",
+ " # Read Q and B\n",
+ " Q, B = ReadEffectiveQuantities(QFilePath,BFilePath)\n",
+ " Q=0.5*(np.transpose(Q)+Q) # symmetrize\n",
+ " B=np.transpose([B])\n",
+ " return (Q,B)\n",
+ "\n",
+ "Q, B=get_Q_B([0,0])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The following function computes the local minimizers based on the assumption that they are on the axes with alpha=0 or alpha=np.pi and |kappa|<=4"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 86,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from scipy.optimize import minimize_scalar \n",
+ "def get_local_minimizer_on_axes(Q,B):\n",
+ " invoke_function=lambda kappa: energy(kappa,0,Q,B)\n",
+ " result_0 = minimize_scalar(invoke_function, method=\"golden\")\n",
+ " invoke_function=lambda kappa: energy(kappa,np.pi/2,Q,B)\n",
+ " result_90 = minimize_scalar(invoke_function, method=\"golden\")\n",
+ " return np.array([[result_0.x,0],[result_90.x,np.pi/2]])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 87,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "n=len(materialFunctionParameter)\n",
+ "m=len(materialFunctionParameter[0])\n",
+ "kappa_0=np.zeros([n,m])\n",
+ "energy_0=np.zeros([n,m])\n",
+ "kappa_90=np.zeros([n,m])\n",
+ "energy_90=np.zeros([n,m])\n",
+ "kappa_exp=np.zeros([n,m])\n",
+ "omega=np.zeros([n,m])\n",
+ "for i in range(0,n):\n",
+ " for j in range(0,m):\n",
+ " Q, B=get_Q_B([i,j])\n",
+ " minimizers=get_local_minimizer_on_axes(Q,B)\n",
+ " kappa_0[i,j]=minimizers[0,0]\n",
+ " energy_0[i,j]=energy(kappa_0[i,j],0,Q,B)\n",
+ " kappa_90[i,j]=minimizers[1,0]\n",
+ " energy_90[i,j]=energy(kappa_90[i,j],0,Q,B)\n",
+ " omega[i,j]=materialFunctionParameter[i][j][3]\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 90,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 99,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAANcAAAB4CAYAAABhN2eOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAG60lEQVR4nO3dsVLb6BrG8Wd3TrdrqT8zSrHFFrEotzjiBsjMFiduVGWGxpMbgC5QZCZVrgAo08SN6UDFmdkK3YBEuwWa2d5iU+sUGRuMsYMEr/2S/f+q8DE2zwQ//myZV/qhaZpGAJ7cj5sOAHyvKBdghHIBRigXYIRyAUYoF2DkX6u+Wde18jzXZDJRmqYL65LU7/cVRZFtSuAZWrlzBUGgfr+vuq7n1kejkZIk0c7Ojk5OTkwDAs9Vp5eFRVEoCAJJUlVVTxoI+F6sfFnY1ng81unpqX777T/6/ff/PuVdd/bly9+SpJ9++nnDSW54y+Qtj+Qv05cvf2tr69dWt+lUrq2tLdV1rSAI5t5vDQYDDQYD/fnnX+r1el3u2oy3PJK/TN7ySD4zPdQ3y5XnuYqiUFVVCsNQZVkqTVOdn58rDMO5Ax0Abvxg8Ye7nnau6+trSb6eAb1l8pZH8pfp+vpav/zy71a34XMuwAjlAoxQLsAI5QKMUC7ACOUCjFAuwAjlAoxQLsAI5QKMUC7AyIMmkaXFiePpehiG6vf7hhGB52nlzrVs4jjLMoVhqCRJVJaleUjgOVq5cxVFoeFwKGl+4jhJEu3u7iqOY+3t7S3cbjro5oGnLFPeMnnLI/nL1CVPp/dcVVVpb29PvV5PR0dHs/XxeKw3b97ojz/+1+Vuge/Kyp1r2cRxnucaDodKkkQfP36crTOJ3I63TN7ySD4zPdTKct2dOK7rWmVZKkkSZVmmKIq0vb29rqzAs8Ik8gZ4y+Qtj+QvE5PIgCOUCzBCuQAjlAswQrkAI5QLMEK5ACOUCzBCuQAjlAswQrkAI5QLMNJ5zH80GimOY1VVpZ2dHduUwDPUecw/iiL1+30lSWIeEniOOo3553mufr+vLMsUBMFCwTyNaHvKMuUtk7c8kr9Maxvzl6Q4jhd2NMb8gRudxvxv//s2xvzb8ZbJWx7JZ6aHWjmJXNf1bMw/iiJFUaSyLBXH8Wz9vpeFnsrlbaJV8pfJWx7JX6Yuk8iM+W+At0ze8kj+MjHmDzhCuQAjlAswQrkAI5QLMEK5ACOUCzBCuQAjlAswQrkAI5QLMEK5ACMry1XXtbIsU5Zlc8OSUycnJ6rr2iwc8JytnOcajUZK01RBEOjw8FDv37+ffa+qqnsLJ/maIvWUZcpbJm95JH+ZnnwSuSgKBUEgSQtFqqpqYWiSSWTgxsqda5k8z5UkiS4vL+fWmURux1smb3kkn5keauXONR3zl+ZH+8MwVJ7nKopiduo1APNW7lxpms7G+dM0VV3XKstSSZLMSjeZTNYSFHhuGPPfAG+ZvOWR/GVizB9whHIBRigXYIRyAUYoF2CEcgFGKBdghHIBRigXYIRyAUYoF2Ck0wXH67pWVVWzy7dyXWRgUacLjp+fnyuKIg2Hw7l1ADc6XXA8TVNJ0uXlpV6+fLlwO08j2p6yTHnL5C2P5C/TWi84LklnZ2fa39+ffc2YP3Cj0wXHJSnLMr19+3buXBqM+bfjLZO3PJLPTA/VaRJZ+vp+7OzsTFEUze1eAL5iEnkDvGXylkfyl4lJZMARygUYoVyAEcoFGKFcgBHKBRihXIARygUYoVyAEcoFGKFcgBHKBRjpPOZ/3zqAG53G/JetA7jRacx/2fqUpxFtT1mmvGXylkfyl6lLnk4XHF9mPB7r9PRUL1680IcPH57yrh9lPB5rMBhsOsYcb5m85ZH8ZWqbp9MFx5etDwYDffr0SVdXV61CWzs9Pd10hAXeMnnLI/nL1DZPpzH/u+t3vX79ul1qY97ySP4yecsj+cvUNs+jx/y9HVH0diLTb/0/nJycKE1TBUGwljzfyjQajRTHsaqq0s7OzsbzTNfDMFS/319LntuZJpPJ3AbS6nHdPNLx8XEzmUyapmmag4ODb65bW/ZzP3/+PFvf3d3deJ6maZqrq6vm4OBg9v1NZzo/P28uLi6apmnWmmlVnrIsm6b5+vtbt6urq+b4+Hhurc3j+tEfIhdFMXvWvXtE8b51a8t+7nR3WHYi03XnmX69ic8Il2XK81xVVSnLstlZvjaZJ0kSHRwc6PDwUK9evVpbnlXaPK7/cX+hcfdEppuS57nLc+zHcezm88uqqrS3t6der6ejo6NNx2nt0eVqe0TR2qqfe/tEppvOE4ah8jxXURSz1/CbzrSpv7RZlmf6BOThyXCqzeP6SQ5oTI8cRlGkKIpUlqXiOJ5bX9eb0WV5pK8HD3q93lpPZLosT5Ikquta79690/b29r1HXded6fbvLAiCtR70uS9PGIazl86TyWTtO/1oNNLFxYX29/cVhmHrx7XJSUEB/APfcwHrQrkAI5QLMEK5ACOUCzBCuQAjlAsw8n9hJ/xnXRcf2AAAAABJRU5ErkJggg==",
+ "text/plain": [
+ "<Figure size 225.682x139.482 with 1 Axes>"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "plt.style.use(\"seaborn\")\n",
+ "mpl.rcParams['text.usetex'] = True\n",
+ "mpl.rcParams[\"font.family\"] = \"serif\"\n",
+ "mpl.rcParams[\"font.size\"] = \"8\"\n",
+ "mpl.rcParams['xtick.bottom'] = True\n",
+ "mpl.rcParams['xtick.major.size'] = 2\n",
+ "mpl.rcParams['xtick.minor.size'] = 1.5\n",
+ "mpl.rcParams['xtick.major.width'] = 0.75\n",
+ "mpl.rcParams['xtick.labelsize'] = 8\n",
+ "mpl.rcParams['xtick.major.pad'] = 1\n",
+ "\n",
+ "mpl.rcParams['ytick.left'] = True\n",
+ "mpl.rcParams['ytick.major.size'] = 2\n",
+ "mpl.rcParams['ytick.minor.size'] = 1.5\n",
+ "mpl.rcParams['ytick.major.width'] = 0.75\n",
+ "mpl.rcParams['ytick.labelsize'] = 8\n",
+ "mpl.rcParams['ytick.major.pad'] = 1\n",
+ "\n",
+ "mpl.rcParams['axes.titlesize'] = 8\n",
+ "mpl.rcParams['axes.titlepad'] = 1\n",
+ "mpl.rcParams['axes.labelsize'] = 8\n",
+ "\n",
+ "#Adjust Legend:\n",
+ "mpl.rcParams['legend.frameon'] = True # Use frame for legend\n",
+ "# mpl.rcParams['legend.framealpha'] = 0.5 \n",
+ "mpl.rcParams['legend.fontsize'] = 8 # fontsize of legend\n",
+ "\n",
+ "\n",
+ "#Adjust grid:\n",
+ "mpl.rcParams.update({\"axes.grid\" : True}) # Add grid\n",
+ "mpl.rcParams['axes.labelpad'] = 3\n",
+ "mpl.rcParams['grid.linewidth'] = 0.25\n",
+ "mpl.rcParams['grid.alpha'] = 0.9 # 0.75\n",
+ "mpl.rcParams['grid.linestyle'] = '-'\n",
+ "mpl.rcParams['grid.color'] = 'gray'#'black'\n",
+ "mpl.rcParams['text.latex.preamble'] = r'\\usepackage{amsfonts}' # Makes Use of \\mathbb possible.\n",
+ "# ----------------------------------------------------------------------------------------\n",
+ "# width = 5.79\n",
+ "# height = width / 1.618 # The golden ratio.\n",
+ "textwidth = 6.26894 #textwidth in inch\n",
+ "width = textwidth * 0.5\n",
+ "height = width/1.618 # The golden ratio.\n",
+ "\n",
+ "fig, ax = plt.subplots(figsize=(width,height))\n",
+ "fig.subplots_adjust(left=.15, bottom=.16, right=.95, top=.92)\n",
+ "\n",
+ "# ax.tick_params(axis='x',which='major', direction='out',pad=3)\n",
+ "\n",
+ "for i in range(0,n):\n",
+ " ax.xaxis.set_major_locator(MultipleLocator(1.0))\n",
+ " ax.xaxis.set_minor_locator(MultipleLocator(0.5)) \n",
+ " ax.yaxis.set_major_locator(MultipleLocator(0.5)) data=np.zeros([3,m])\n",
+ " data[0]=omega[i,::-1]\n",
+ " data[1]=kappa_0[i,::-1]\n",
+ " data[2]=kappa_90[i,::-1]\n",
+ "\n",
+ " # relative_error = (np.array(data[dataset_number][1]) - np.array(dataset[dataset_number][2])) / np.array(dataset[dataset_number][2])\n",
+ " #print('relative_error:', relative_error)\n",
+ "\n",
+ " #--------------- Plot Lines + Scatter -----------------------\n",
+ " line_1 = ax.plot(np.array(data[0]), np.array(data[1]), # data\n",
+ " # color='forestgreen', # linecolor\n",
+ " marker='D', # each marker will be rendered as a circle\n",
+ " markersize=3.5, # marker size\n",
+ " # markerfacecolor='darkorange', # marker facecolor\n",
+ " markeredgecolor='black', # marker edgecolor\n",
+ " markeredgewidth=0.5, # marker edge width\n",
+ " # linestyle='dashdot', # line style will be dash line\n",
+ " linewidth=1, # line width\n",
+ " zorder=3,\n",
+ " label = r\"$\\kappa_{1,sim}$\")\n",
+ "\n",
+ " line_2 = ax.plot(np.array(data[0]), np.array(data[2]), # data\n",
+ " color='red', # linecolor\n",
+ " marker='s', # each marker will be rendered as a circle\n",
+ " markersize=3.5, # marker size\n",
+ " # markerfacecolor='cornflowerblue', # marker facecolor\n",
+ " markeredgecolor='black', # marker edgecolor\n",
+ " markeredgewidth=0.5, # marker edge width\n",
+ " # linestyle='--', # line style will be dash line\n",
+ " linewidth=1, # line width\n",
+ " zorder=3,\n",
+ " alpha=0.8, # Change opacity\n",
+ " label = r\"$\\kappa_{2,sim}$\")\n",
+ "\n",
+ " #line_3 = ax.plot(np.array(data[0]), np.array(data[3]), # data\n",
+ " # # color='orangered', # linecolor\n",
+ " # marker='o', # each marker will be rendered as a circle\n",
+ " # markersize=3.5, # marker size\n",
+ " # # markerfacecolor='cornflowerblue', # marker facecolor\n",
+ " # markeredgecolor='black', # marker edgecolor\n",
+ " # markeredgewidth=0.5, # marker edge width\n",
+ " # # linestyle='--', # line style will be dash line\n",
+ " # linewidth=1, # line width\n",
+ " # zorder=3,\n",
+ " # alpha=0.8, # Change opacity\n",
+ " # label = r\"$\\kappa_{exp}$\")\n",
+ "\n",
+ " # --- Plot order line\n",
+ " # x = np.linspace(0.01,1/2,100)\n",
+ " # y = CC_L2[0]*x**2\n",
+ " # OrderLine = ax.plot(x,y,linestyle='--', label=r\"$\\mathcal{O}(h)$\")\n",
+ "\n",
+ "\n",
+ "\n",
+ " # Fix_value = 7.674124\n",
+ " # l3 = plt.axhline(y = Fix_value, color = 'black', linewidth=0.75, linestyle = 'dashed')\n",
+ " # --------------- Set Axes -----------------------\n",
+ " # ax.set_title(r\"ratio $r = 0.22$\") # Plot - Title\n",
+ "\n",
+ " # Plot - Titel\n",
+ " ax.set_title(r\"ratio $r = 0.49$\") \n",
+ " ax.set_xlabel(r\"Wood moisture content $\\omega (\\%)$\", labelpad=4)\n",
+ " ax.set_ylabel(r\"Curvature $\\kappa$($m^{-1}$)\", labelpad=4)\n",
+ " plt.tight_layout()\n",
+ "\n",
+ " # # --- Set Line labels\n",
+ " # line_labels = [r\"$CC_{L_2}$\",r\"$CC_{H_1}$\", r\"$\\mathcal{O}(h)$\"]\n",
+ "\n",
+ " # --- Set Legend\n",
+ " legend = ax.legend()\n",
+ " # legend = fig.legend([line_1 , line_2, OrderLine],\n",
+ " # labels = line_labels,\n",
+ " # bbox_to_anchor=[0.97, 0.50],\n",
+ " # # bbox_to_anchor=[0.97, 0.53],\n",
+ " # # loc='center',\n",
+ " # ncol=1, # Number of columns used for legend\n",
+ " # # borderaxespad=0.15, # Small spacing around legend box\n",
+ " # frameon=True,\n",
+ " # prop={'size': 10})\n",
+ "\n",
+ "\n",
+ " frame = legend.get_frame()\n",
+ " frame.set_edgecolor('black')\n",
+ " frame.set_linewidth(0.5)\n",
+ "\n",
+ "\n",
+ " # --- Adjust left/right spacing:\n",
+ " # plt.subplots_adjust(right=0.81)\n",
+ " # plt.subplots_adjust(left=0.11)\n",
+ "\n",
+ " # ---------- Output Figure as pdf:\n",
+ " fig.set_size_inches(width, height)\n",
+ " fig.savefig('WoodBilayer_expComparison_local_'+str(i)+'.pdf')\n",
+ " plt.cla()\n",
+ " \n",
+ "\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "6"
+ ]
+ },
+ "execution_count": 36,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Sampling of local energy"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "N=500\n",
+ "length=4\n",
+ "r, theta = np.meshgrid(np.linspace(0,length,N),np.radians(np.linspace(0, 360, N)))\n",
+ "E=np.zeros(np.shape(r))\n",
+ "for i in range(0,N): \n",
+ " for j in range(0,N): \n",
+ " if theta[i,j]<np.pi:\n",
+ " E[i,j]=energy(r[i,j],theta[i,j],Q,B) * energyscalingfactor\n",
+ " else:\n",
+ " E[i,j]=energy(-r[i,j],theta[i,j],Q,B) * energyscalingfactor\n",
+ "# \n",
+ "# Compute Minimizer\n",
+ "[imin,jmin]=np.unravel_index(E.argmin(),(N,N))\n",
+ "kappamin=r[imin,jmin]\n",
+ "alphamin=theta[imin,jmin]\n",
+ "# Positiv curvature region\n",
+ "N_mid=int(N/2)\n",
+ "[imin,jmin]=np.unravel_index(E[:N_mid,:].argmin(),(N_mid,N))\n",
+ "kappamin_pos=r[imin,jmin]\n",
+ "alphamin_pos=theta[imin,jmin]\n",
+ "Emin_pos=E[imin,jmin]\n",
+ "# Negative curvature region\n",
+ "[imin,jmin]=np.unravel_index(E[N_mid:,:].argmin(),(N_mid,N))\n",
+ "kappamin_neg=r[imin+N_mid,jmin]\n",
+ "alphamin_neg=theta[imin+N_mid,jmin]\n",
+ "Emin_neg=E[imin+N_mid,jmin]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Local minimizer in area of pos. curvature: kappa = 1.819639278557114 , alpha = 0.0 , Q = 0.1299383783013474\n",
+ "Local minimizer in area of neg. curvature: kappa = 2.797595190380761 , alpha = 4.709241091954239 , Q = 0.09277799602960847\n",
+ "2293.056976484917\n",
+ "3317.472225915116\n"
+ ]
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGdCAYAAAAMm0nCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAB9bUlEQVR4nO3deXxU5d3//9c5kwUQEgkJBGSbGFYXUASMehtiCGBTKi3VW+u+9Y6EIORW+Fq1bjfFG9oIBGLanxW0oha0aDEKTcMw3JZUEAyKC1smQYVAAMkgYpY55/fHcE5mkgESSDLb5/l48NAzc0iuIcOcN9f5XJ9L0XVdRwghhBAiiKj+HoAQQgghRGtJgBFCCCFE0JEAI4QQQoigIwFGCCGEEEFHAowQQgghgo4EGCGEEEIEHQkwQgghhAg6EmCEEEIIEXQi/D2A9qJpGvv376dbt24oiuLv4QghhBCiBXRd5/jx4/Tp0wdVPf08S8gGmP3799OvXz9/D0MIIYQQ5+Drr7+mb9++p30+ZANMt27dAPcfQExMjJ9HI4QQQoiWcDqd9OvXz7yOn07IBhjjtlFMTIwEGCGEECLInK38Q4p4hRBCCBF0JMAIIYQQIuhIgBFCCCFE0JEAI4QQQoigIwFGCCGEEEFHAowQQgghgo4EGCGEEEIEHQkwQgghhAg6EmCEEEIIEXQkwAghhBDirMrLF+Fw5Pt8zuHIp7x8UYeORwKMEEIIIc5KUVTKHQubhRiHI59yx0IUpWMjRcjuhSSEEIGk9M+PgGoh5d7/bf7csjmguUi5//d+GJkQLWO15gBQ7lhoHhvhJck603y+o8gMTAsscBwgr6LK53N5FVUscBzo4BEJIYKOaiGlstAdVjyULptDSmUhqBY/DUyIlrNac0iyzqTcsZD1tmF+Cy8gAaZFLIrCfEdVsxCTV1HFfEcVlrPsmCmEECn3/i+lA7K8QowRXkoHZPmcmREiEFmtOShKFLpeh6JE+SW8gNxCapHcgYkAzHdUmcdGeJltTTSfF0KIM0m5938pXQYplYXUPfUSKUqDhBcRdByOfDO86HodDke+X0KMouu63uHftQM4nU5iY2OpqakhJiamTb6mEVqiFIU6XZfwIoQ4J3VP9SBKaaBOjyDqmSP+Ho4QLda05qU9amBaev2WW0itkDsw0QwvUYoi4UUI0Wqly+Y0hheloVlNjBCByldY8ayJOd0S6/YiAaYV8iqqzPBSp+unLewVQghfPGteop450qwmRohApuuaz5kWI8Toutah45EamBZqWvNiHAMyEyOEOCtfBbueNTGly5BaGBHQkpIePu1z/qiBkQDTAr4Kdn0V9gohxGlpLp8Fu0aIQXP5Z1xCtEBNcSWKqhCT3r/Zc86SfeiaTmzGgA4dk9xCagHXaQp2cwcmMtuaiCs066CFEG0o5f7f88kVAyjcXtj8uXv/l+1XJVNQVuCHkQlxdoqq4CyuxFmyz+txZ8k+nKfCTUeTGZgWeNTa+7TPycyLEKKlVEVladlSALJGZJmPF24vZGnZUrJHZvtraEKckTHz4iyuNI+N8BKTMcDnzEx7kwAjhBAdxAgtniHGM7x4hhohAo1niHGu3wcu3W/hBSTACCFEh/IMMX/69E/Ua/USXkTQiEnvb4YXLL5rYjqK1MAIIUQ7qc5fQnVB87qWrBFZ3PwvmGL/kUg1UsKLCBrOksbwgktvVhPTkVoVYF588UUuv/xyYmJiiImJISUlhQ8++ACAo0ePkpOTw5AhQ+jcuTP9+/dnxowZ1NTUeH0NRVGa/XrzzTe9ztmwYQNXXnkl0dHRJCcns3z58vN7lW3IZrNht9t9Pme327HZbB08IiFEwLKoHF6c3yzEvP/be7l5YwOoFuq1ep+FvUIEGs+al75zryMmY4DPwt6O0qpbSH379uX5559n0KBB6LrOK6+8wk033cQnn3yCruvs37+f3//+9wwfPpzKykqysrLYv38/b731ltfXWbZsGZMmTTKPL7zwQvP/HQ4HmZmZZGVlsWLFCkpKSnjggQfo3bs3EydOPL9X2wZUVTVDSmpqqvm4EV7S0tL8NTQhRIBJmDYNgMOL883j9397L9aV/8Zxy9X89tll9DxVAwPITIwIWL4Kdn0V9nak894LKS4ujgULFnD//fc3e27VqlXccccdnDhxgogId1ZSFIXVq1czZcoUn19vzpw5FBUVsWPHDvOxW2+9lWPHjrF27doWj6s99kIyeIaV1NTUZsdCCOGpuqCAw4vz0SIsqA0uHLdczU+eXWY+L4W8ItB1ZB+Yll6/z7mI1+VysWrVKk6cOEFKSorPc4xvboQXQ3Z2Ng888ABJSUlkZWVx7733oijuNeSlpaWMHz/e6/yJEycyc+bMM46ntraW2tpa89jpdJ7Dq2oZI6TYbDY2btyIy+WS8CKEOK2EadM48mIhan09WoTFK7xA48yL1sGt2IU4m81rylFUhdGZ1mbPbSlyoGs6YyYn+WFk51DE+9lnn9G1a1eio6PJyspi9erVDB8+vNl5hw8f5rnnnuPXv/611+PPPvssK1eupLi4mKlTpzJt2jTy8xs3gKqqqqJXr15ev6dXr144nU5Onjx52nHNmzeP2NhY81e/fv1a+9JaJTU1FYvFgsvlwmKxSHgRQpxWdUEBen09SmQkaoML229XsqXI4XVO1ogspo2cxpYiB5vXlPtppEJ4U1SFzWsczd6v7vepwy8N7AytnoEZMmQIZWVl1NTU8NZbb3H33Xdjt9u9QozT6SQzM5Phw4fz9NNPe/3+J5980vz/K664ghMnTrBgwQJmzJhx7q8CeOyxx8jNzfUaQ3uGGLvdboYXl8uF3W6XECOEaMa4fRQ/I4eEadOoLijA8W4pXxyKB/D6l61xURgzufm/doXwB+P9uXmNwzz2fJ/6mpnpKK0OMFFRUSQnJwMwatQotmzZwqJFi/jjH/8IwPHjx5k0aRLdunVj9erVREZGnvHrjR07lueee47a2lqio6NJTEzk4MGDXuccPHiQmJgYOnfufNqvEx0dTXR0dGtfzjk5XQ0MICFGCGFqGl7AfTtpNAXw7ho2r5kMBNZFQYimPEPMxx9UoDXoAfE+Pe9GdpqmmbUnTqeTiRMnEh0dzd///nc6dep01t9fVlZG9+7dzfCRkpLC+++/73VOcXHxaetsOpqvgl3PmhjPYyFEmHNpXuHFYISYLgcOsXkNAXVREMKX0ZlW832qRviuielorQowjz32GDfeeCP9+/fn+PHjvP7662zYsIF169bhdDqZMGECP/zwA6+99hpOp9MspE1ISMBisbBmzRoOHjzI1VdfTadOnSguLuZ3v/sdjzzyiPk9srKyWLJkCbNnz+a+++5j/fr1rFy5kqKiorZ95edI0zSfBbvGsaZJEZ4Qwi0hZ/rpn5s2jXHAl9NtAXVREMKXLUUO832qNehsKXL4/f3aqgBz6NAh7rrrLg4cOEBsbCyXX34569atIyMjgw0bNvDRRx8BmLeYDA6Hg4EDBxIZGcnSpUuZNWsWuq6TnJxMXl4eDz74oHmu1WqlqKiIWbNmsWjRIvr27ctLL70UED1ggDP2eZGZFyFEawTiRUGIppre3jSOAb++X8+7D0ygas8+MEIIcb5Od1GQ20gikJzufdme79d27wMjhBDi3Pj68Pe12kMIf3P3eWkeUoxjXfPfHIgEGCGE6GCBfFEQAqCgrABVUcma3LwzdOH2QjRdY1rmNB+/s+NIgDkPm1atQFFVUqbe1uy50rffQNc0rrn5dj+MTAjhL635XPDVnt0IMW3dnl2I1lAV1eceXZ7bXvhbqzvxikaKqrJp5QpK337D6/HSt99g00r3h5gQIry05nNBURWfu/kaG+f5s8upCG9ZI7LIHpnN0rKl5m7pgbZnl8zAnAfjX1ibVq4wj40PqWtuud3nv8CEEKGtNZ8Lvnbz9bXrrxD+YISUpWVL+dOnf6Jeqw+Y8AKyCqlNGB9OlogIXA0NEl6EEK36XDBCCxYFXLqEFxFQrvzLldRr9USqkWy7c1u7f7+WXr/lHkcbSJl6m/khZYmIkPAihGjV50JMen8zvGBRJLyIgFG4vdAML/VavXk7KRBIgGkDpW+/YX5IuRoamt37FkKEn9Z8LjhL9pnhBZferCZGCH/wrHnZdue2ZjUx/iY1MOep6b1t4xiQmRghwlRrPhea1ryYt5NAZmKE3/gq2PWsifE89hcJMOfBV2GerwI+IUT4aM3ngq+CXV+FvUJ0NE3XfBbsGsea7v99/yTAnAdd03wW5hnHumzsKETYac3ngq75Ltg1jqWhnehoLxTvwqIqzEhv3qRuccluXJrOrAxZhdSuZC8kIUSwKC9fhKKoWK05zZ5zOPLRdY2kpIf9MDIRbhaX7CaveBe5GYOZkT7orI+3B9kLSQghgoSiqJQ7FgJ4hRiHI59yx0KSrDP9MzARdoxwkle8yzzuyPDSGhJghBDCz4zQ4hliPMOLr5kZIdqLZ4hZsn4PdS4t4MILSIARQojzZrPZUFWV1NTUZs/Z7XY0TSMtLe2MX8MzxDgqCtD1Ogkvwm9mpA8yw0uURQ248ALSB6bNVOcvobqgwPdzBQVU5y/p4BEJITqKqqrYbDbsdrvX43a73Qw3LWG15qAoUeh6HYoSJeFF+M3ikt1meKlzaSwu2e3vITUjAaatWFQOL85vFmKqCwo4vDgfLPJHLUSoSk1NJS0tzSvEGOElLS3N58yML+6C3TozxDgc+e05bCF88qx52TX3RnIzBpNXvCvgQozcQmojCdPcS84OL843j43wEj8jx3xeCBGajJBis9nYuHEjLper1eHFs+bFOAZkJkZ0GF8Fu74KewOBBJg25BlijrxYiF5fL+FFiDCSmppqhheLxXLO4QV8F/YK0V6MpfwubWKzgl2HI5+fWjXI+AmuAOpNJPc12ljCtGkokZHo9fUokZESXoQII3a73QwvLperWU3M6ei65rNg12rNIck6Ez0Aup6K0GYs5Z+SvK5ZeCl3LERR3IW8szIG+3GU3mQGpo1VFxSY4UWvr6e6oEBCjBBhoGnNi3EMnHUmxmhS19gFtfECYoSaxi6ogXMBEaEjGJfyS4BpQ01rXswCXpAQI0QI81Ww61kT43l8JhZV8Vln4FmXIER7Cbal/BJg2oivgl1fhb1CiNBj9HlpGlKMY62F+6IFUxdUEZqs1hwzvAT6Un4JMG3Fpfks2DWPXXIPW4hQdaYmdS0t5DUESxdUEZp8LeUP1BAjmzkKIUQAGvz4B2YjsV1zb/T3cEQYON1S/o6+jdTS67esQhJCiAATDF1QRWg53VL+JOtMd01MADZVlFtIQggRQJrWvBjHEDgNxERoqSmu5Ef1O5KSms+09Ci/iR/1moBcyi8BRgghztECxwEsikLuwMRmz+VVVOHSdR619m7x1wumLqgidCiqQrfiNGIyBoC18XFnyT6cxZVclHEXMUn9/TfA05BbSG3NNg/s830/Z5/vfl4IERIsisJ8RxV5FVVej+dVVDHfUYVFUVr19VyaTuqYbUTFl3g9PiN9ELkZg9las5KCMt+bxgpxrmLS+xOTMQBncSXOkn1AY3iJyRhATHrghReQGZi2p1rANtf9/6mzGx+3z3c/nva4f8YlhGhzxszLfEeVeWyEl9nWRJ8zM2cyK2MwhdsTWFq2FICsEVnmc1HxJWz75k1SlOw2Gr0QjYyQ4iyuxLl+H7j0gA4vIAGm7RmhxTPEeIYXz1AjhAh6niFmYcVB6nT9nMKLwQgtniGmcHshS8uWkj0y2yvUCNGWYtL7m+EFixLQ4QVkGXX7MUKLJQpcdRJehAhx/Tdsp07XiVIU9o0bcd5fzwgtkWok9Vq9hBfR7ozbRlgUv87AyDJqf0ud3RheLFESXoQIYXkVVWZ4qdP1ZjUx5yJrRJYZXiLVSAkvol151rz0nXtds5qYQCQBpr3Y5zeGF1fd6Qt7hRBBzbPmZd+4Ecy2Jvos7G2twu2FZnip1+op3F7YRiMWwpuvgl1fhb2BplUB5sUXX+Tyyy8nJiaGmJgYUlJS+OCDD8znf/zxR7Kzs+nRowddu3Zl6tSpHDx40Otr7Nu3j8zMTLp06ULPnj159NFHaWho8Dpnw4YNXHnllURHR5OcnMzy5cvP/RX6g2fNy5PV7v/a5kqIESLE+CrYzR2YeN4hxrPmZdud28gemc3SsqUSYkS70DXft4uMEKNrgVlp0qoi3r59+/L8888zaNAgdF3nlVde4aabbuKTTz7hkksuYdasWRQVFbFq1SpiY2OZPn06v/jFL/jXv/4FgMvlIjMzk8TERDZt2sSBAwe46667iIyM5He/+x0ADoeDzMxMsrKyWLFiBSUlJTzwwAP07t2biRMntv2fQFvzVbDrq7BXCBH0XKcp2DWOXedQYuirYNdXYa8Q52vzmnIUVWF0prXZc1uKHOiazpjJSX4YWcucdxFvXFwcCxYs4Je//CUJCQm8/vrr/PKXvwTgq6++YtiwYZSWlnL11VfzwQcf8NOf/pT9+/fTq1cvAAoLC5kzZw7V1dVERUUxZ84cioqK2LFjh/k9br31Vo4dO8batWtbPC6/FfHa5rmXUvsKKfb5oLkg7bGOG48QIqgUlBWgKipX7pvQ7OJSuL0QTdcY/e2NAX9xEYFvS5GDzWscjJls9Xqfne7xjtLuRbwul4s333yTEydOkJKSwtatW6mvr2f8+PHmOUOHDqV///6UlpYCUFpaymWXXWaGF4CJEyfidDr5/PPPzXM8v4ZxjvE1Tqe2than0+n1yy/SHjv9DEvqbAkvQogzmjZyGlkjslBUhc1rHGwpcpjPZY3IYvS3N7J5jQNFbV2TPCGaGp1pZcxkq9f7zN/hpTVa3Qfms88+IyUlhR9//JGuXbuyevVqhg8fTllZGVFRUVx44YVe5/fq1YuqKvd94KqqKq/wYjxvPHemc5xOJydPnqRz584+xzVv3jyeeeaZ1r4cIYQISMbFY/Mah3kcTBcXERw832cff1CB1qAHzfur1QFmyJAhlJWVUVNTw1tvvcXdd9+N3W5vj7G1ymOPPUZubq557HQ66devnx9HJIQISR14mziYLy4ieIzOtJrvLzXCd01MIGr1LaSoqCiSk5MZNWoU8+bNY8SIESxatIjExETq6uo4duyY1/kHDx4kMdFd0JaYmNhsVZJxfLZzYmJiTjv7AhAdHW2ujjJ++dsCx4HTrkLIq6higeNAB49ICHHejO1Cmq4qNAr4VUubfrvRmVbUCCXoLi4ieGwpcpjvL61B97ptGcjOuw+MpmnU1tYyatQoIiMjKSlp3IRs586d7Nu3j5SUFABSUlL47LPPOHTokHlOcXExMTExDB8+3DzH82sY5xhfI5i09UZvQogAkDq7eWuEdtwuJFgvLiI4eN6WfGhJWrOamEDWqltIjz32GDfeeCP9+/fn+PHjvP7662zYsIF169YRGxvL/fffT25uLnFxccTExJCTk0NKSgpXX301ABMmTGD48OHceeedzJ8/n6qqKp544gmys7OJjo4GICsriyVLljB79mzuu+8+1q9fz8qVKykqKmr7V9/O2nqjNyFEgPBsjbBxQbttF9K05sU4BmQmRpw3XzVVvmqvAlWrAsyhQ4e46667OHDgALGxsVx++eWsW7eOjIwMAF544QVUVWXq1KnU1tYyceJECgoat363WCy89957PPTQQ6SkpHDBBRdw99138+yzz5rnWK1WioqKmDVrFosWLaJv37689NJLwdEDxoe23uhNCBEgUmc3hpd22C4k2C8uInCZS/W1Cc1qqgq3F6JdpDFm8o0B28DOIJs5dpC23uhNCOFn7bxhq9FkbHCUiqJ67wxsNBkb2iUCXdOJzRjQZt9XhL7T7W4eKLuey2aOAaQ9NnoTQvhRB2wXMmZyEqMzrSiq0mw/mtGZVoZ2icBZXCn9YESrZY3IarY9RaCEl9Zo9TJq0TpNa16MY0BuIwkRjDp4uxBj5sVZXGke+9p8T4jW8Nye4k+f/ol6rT6owgtIgGlXp9voDZAQI0Sw0ly+bxcZx5qrzb+lZ4hxrt8HLt+b7wnRGlkjsszwEqlGBlV4AQkw7ao9NnoTQvjZmZrUteNGrTHp/c3wgkWR8CLOW+H2QjO81Gv1FG4vDKoQIwGmHT1q7X3a52TmRQjRGs6SxvCCS8dZsk9CjDhnTWtejGMInt3OJcAIIUSAa1rzYhwDEmJEi71QvAuLqhAVX9KsYLfucDpXdjsSVCFGAowQQpxFdf4SsKgkTJvW/LmCAnBpJORMb5fv7atg11dhrxBnY1EV8op3kTqm2iu8LC7ZTV7xLnIzbiHl4h5ouubnkbaMBJgOYrPZUFWV1NTUZs/Z7XY0TSMtLc0PIxNCnJVF5fDifACvEFNdUMDhxfnEz8hpt2+ta+6C3cPW1RxxqFit7u9lhBZd03E48tF1jaSkh9ttHCL4zUgfBEBeMYyKHQx4hpfBp54f5McRto4EmA6iqio2mw3AK8TY7XZsNpuEFyECmBFaPEOMZ3jxNTPTVowmdUccKuWOhQBeIcbhyKfcsZAk68x2G4MIHY0hZhdL1u+hzqV5hJfgIgGmgxihxTPEeIYXXzMzQojA4RlijrxYiF5f3+7hxZMRWjxDjGd4MZ4X4mxmpA8yw0uURQ3K8AKylUCHM0KLxWLB5XJJeBEiyHx12eXo9fUokZEM/ezTDv/+RmhRlCh0vU7Ci2g147ZRlEUNyBkY2UogQKWmpprhxWKxSHgRIohUFxSY4UWvr3cX8HYwqzXHDC+KEiXhRbSKZ83Lrrk3kpsxmLziXSwu2e3vobWaBJgOZrfbzfDicrmw2+3+HpIQogU8a16GfvYp8TNyOLw4v8NDjLtgt84MMQ5Hfod+fxG8mhfsum8nBWuIkRqYDtS05sU4BmQmRogA5qtg11dhb3trWvNiHAMyEyNOq7x8EYqi4tImNrtd5HDk81OrBhk/waUFV0WJBJgO4qtg11dhrxAiALk0nwW75rGr/ftmGGFl45Fn+Fwfywxr88LeNeWTcGk6szIGt/t4RPBQFPcKtinJ3kHXMxAHUg1MS0mA6SBGn5emIcU41rTgaBwkRDg6U5O6jlqFpOsaSdaZfK6PJa94F+Ce/jcuSK9ujWXZx+7bA0J4CtUVbLIKSQghgkzTWgZftQ1CNBUsK9haev2WACOEEEEo0JfCisC03jbMLAK/Ie1Lfw/HJ1lGLYQQ52HTqhWUvv2Gz+dK336DTatWdPCIvM1IH2SGl2BuRiY6TqitYJMA4weB/sEohABFVdm0svnf1dK332DTyhUoqn8/PheX7DbDS51LC7olsKJjeda83JD2JUnWmZQ7FgZ1iJEiXj8wPhgBUqbeZj5ufDBec8vt/hqaEOIU4++m599Vz7+jnn93O9rpamAAmYkRXmqKKzmg/oVvlWVeNS9Waw615TVBvQxfAowfBPIHoxCikeff1Y/+9ldcDQ1+/zt6umZkgIQY0YyiKvxYfpSLku71CinOkn10LU7logz3CrdgJAHGTwLxg1EI0VzK1NvMv6OWiAi//x11aTqpY7YRFf810BhUjNCytWYlBWXdmTayY5Z3i8AWk96fJGbiLK7EyT5i0vvjLNmHs7iSmIwB9E1/wt9DPGdSA+NHKVNvwxIRETAfjEKI5krffsP8O+pqaDht/VpHmZUxmJSLE1hatpTC7YVez0XFl7Dt+Juoiny0i0Yx6f2JyRiAs7iSbx7/0AwvMen9/T208yLvcj8KtA9GIYQ3z1u7M1e8wzW33O6zsLejZY3IIntktleIKdxeyNKypWSPzCZrRJZfxycCT0x6f7Ao4NLBogR9eAG5heQ3TWtejGNAZmKECAC+6tJ81a/5ixFSlpYt5U+f/ol6rV7CizgtZ8k+M7zg0nGW7Av6ECMBxg8C/YNRCAG6pvmsSzOO9QDY/iNrRJYZXiLVSAkvwifPmhfPGhggqEOMBBg/CIYPRiHC3TU3n76dQaD8A6Nwe6EZXuq1egq3F0qIEV6ahhdoDC3BHmIkwPhBMHwwCiECW9OaF+MYkBAj2LymHEVVGBylNivY3VLkQNd0hmYMQNeCdzchCTBCCOHBZrOhqmqzneMB7Ha7ubO8Pxlh5aH63zLqm6tghHdNDMCobyaiazpjJif5c6jCTxRVYfMaB0y2MjrTaj6+pcjB5jUOxky2Bu3Mi0ECjB9V5y8Bi0rCtOb9GqoLCsClkZAz3Q8jEyJ8qaqKzWYD8Aoxdrsdm83m9/ACoOka2SOzGfXNVe6LFDA602qGGNfH3dm8xX2REuHJCC2e7w/P8OIZaoKVBBh/sqgcXuzeh8IzxFQXFHB4cT7xM4KvtbMQwc4ILZ4hxjO8+JqZ6Whmk7oR7v94XqRGfTPRDC+hcJES584zxHz8QQVagx5S7wsJMH5khBbPEOMZXnzNzAgh2p9niNm4cSMulytgwktToX6REudndKbVfF+oEUpIvS9a1chu3rx5jB49mm7dutGzZ0+mTJnCzp07zecrKipQFMXnr1WrVpnn+Xr+zTff9PpeGzZs4MorryQ6Oprk5GSWL19+fq80QCVMm0b8jBwOL87nq8sul/AiRIBITU3FYrHgcrmwWCwBGV4MozOtqBFKSF6kxPnZUuQw3xdag86WIoe/h9RmWhVg7HY72dnZ/Pvf/6a4uJj6+nomTJjAiRMnAOjXrx8HDhzw+vXMM8/QtWtXbrzxRq+vtWzZMq/zpkyZYj7ncDjIzMwkLS2NsrIyZs6cyQMPPMC6devO/xUHoIRp01AiI9Hr61EiIyW8CBEA7Ha7GV5cLhd2u93fQzqtUL5IiXPnWfPy0JI0xky2snmNI2TeH626hbR27Vqv4+XLl9OzZ0+2bt3K9ddfj8ViITEx0euc1atXc8stt9C1a1evxy+88MJm5xoKCwuxWq384Q9/AGDYsGF8+OGHvPDCC0ycOLE1Qw4K1QUFZnjR6+upLiiQECOEHzWteTGOgYCbiWlamGkcAzITE4YKygpQFdVdC9WkYHdr33Uwujub17jPDfb3x3nthVRTUwNAXFycz+e3bt1KWVkZ999/f7PnsrOziY+PZ8yYMbz88svoeuNa9NLSUsaPH+91/sSJEyktLT3tWGpra3E6nV6/goFnzcvQzz41bydVFxT4e2hChCVfBbupqamkpaVhs9kCaibGCCs3XB7HkE4WwH1R8vyXtrNkHzWnGpaJ0KcqKkvLlrLlwMde4cVYem+56jvGTLYGdf8XwzkX8WqaxsyZM7n22mu59NJLfZ7z5z//mWHDhnHNNdd4Pf7ss89yww030KVLF/7xj38wbdo0vv/+e2bMmAFAVVUVvXr18vo9vXr1wul0cvLkSTp37tzse82bN49nnnnmXF+OX/gq2PVV2CuE6DhGn5emMy3GsRZAnbLdfV6sXNTJ4tVV1bhoXVBZg3PfcWIyBvhzmKIDNfYDeha1bzajyWq+0ecIPw+yjZxzgMnOzmbHjh18+OGHPp8/efIkr7/+Ok8++WSz5zwfu+KKKzhx4gQLFiwwA8y5eOyxx8jNzTWPnU4n/fr1O+ev1yFcms+CXfPYFTgflEKEizP1eQm020dNm9R5hpghnSxmeAn2hmWidcJlo89zCjDTp0/nvffeY+PGjfTt29fnOW+99RY//PADd91111m/3tixY3nuueeora0lOjqaxMREDh486HXOwYMHiYmJ8Tn7AhAdHU10dHTrX4wfnalJncy8CCFaw3N/G+d6987DEl7CVzhs9NmqGhhd15k+fTqrV69m/fr1WK2nLwD685//zM9+9jMSEhLO+nXLysro3r27GUBSUlIoKSnxOqe4uJiUlJTWDFcIIcJKTHp/sCjg0sGiSHgJY742+gw1rZqByc7O5vXXX+fdd9+lW7duVFVVARAbG+s1M7Jnzx42btzI+++/3+xrrFmzhoMHD3L11VfTqVMniouL+d3vfscjjzxinpOVlcWSJUuYPXs29913H+vXr2flypUUFRWd6+sMfLZ5oFogdXbz5+zzQXNB2mMdPy4hwsACxwEsikLuwOYrI/MqqnDpOo9ae/thZK3jLNlnhhdcOs6SfRJiwsQLxbuwqAoz0gc1q3m5Z/XckNzos1UB5sUXXwRg3LhxXo8vW7aMe+65xzx++eWX6du3LxMmTGj2NSIjI1m6dCmzZs1C13WSk5PJy8vjwQcfNM+xWq0UFRUxa9YsFi1aRN++fXnppZdCcgm1SbWAba77/z1DjH2++/G0x/0zLiHCgEVRmO9w/4PMM8TkVVQx31HFbKvvlg+BxFmyD2dxpXnbyDgGJMSEAYuqkFe8i23OlWx1vmmGl8Ulu9nw0UjGjSXkQoyie65fDiFOp5PY2FhqamqIiYnx93BaxjOspM5ufiyEaDeeYSV3YGKz40BmhJXjGTY6JXXHas3xejwmYwBHkt5F1zWSkh7282hFe1lcspslnyzlmosTWP7zx1lcspu84l3kZgw2Z2Y0XWvcSytAtfT6LQEm0BihxRIFrjoJL0J0ICO0RCkKdboeFOEFoKa4EkVVOJL0LuWOhSRZZ3qFmG/1V9mvLvN6XIQmI7REWVTqXJoZXoJJS6/f59XITrSD1NmN4cUSJeFFiA6UOzDRDC9Rp6mJCUSxp24bWa05JFlnUu5YiMPh7iV1JOldCS9hZEb6IDO8RFnUoAsvrSEBJtDY5zeGF1ed+1gI0SHyKqrM8FKn6+RVVPl7SK3mGWLW24Y1m5ERoW1xyW4zvNS5NBaX7Pb3kNqNBJhA4lnz8mS1+7+2uRJihOgAnjUv+8aNYLY1kfmOqqANMYoSha7XoShREl7ChGfNy665N5KbMZi84l0hG2LOuROvaGO+CnaN//panSSEaDO+CnaN//panRToHI58M7zoeh0OR76EmBBVXr4IRVFZUz7Jq2AXYHLSWr67Kpa84l0AIXc7SQJMoNBcvgt2jWPN1fFjEiJMuE5TsGscu4JorYPDkc8LxTvpEfc0T/3yThyOfModCwH3zMzikt24NJ1ZGYP9Ok7RNhRFpdyxkCNHu5ObMdYMKcbP/a5RM+nePQVXCGze2JQEmEBxpiZ1MvMiRLs6U5O6YJt5KXcspEfc0yz7OI7u3XczI90981LuWMirWy9k2cdx5Ep4CRmNM2tPkWSdCQwy3wdG7dOM0zfND2oSYIQQ4SkEu1/rukaSdSbpN9xJ9+67PW4d5HiFl1C7lRDujBBT7liIo6IAXa8Li8JtCTABKFTamgsR0EKw+7VnkzojpOQV72LJ+j3UuSS8hDKrNccML+FSuC2rkAKQ0da86eoHo9DQoih+GpkQISR1dvOVfiHW/TqceoKEO1+F26FOZmACkK/VD8HU1lyIoOG50m/jgpDrfu2rJ4iEmNBxug7MRg3Mj3trGDb+CX8Ps91IgAlQniFmYcXBoGprLkRQSZ3dGF5CqPt1031wjGMIveW04UpRFcr3LuaIstqr5qVH+U3U7K1kf/IyOjliQ/Z2kgSYAJY7MNEML8HU1lyIoOKr+3WQhxhjU79xYxOYkZ4JeNfEbHOuZIy1e8Bv6ifOLCa9P9HE0GPPz+mh3wTWxg08kzJmEGsdgK5r/h5mu5EAE8B8tTWXECNEGzrdDvAQ1CHGpelcc3ECW51vUri9B1kjsgB3iNnmXMlW55tcrWT7eZSiLQxL/y1O3KHFuX4fuHRiTu2NFUNozrwYJMAEqKY1L8YxBFdfCiECVgh3v3Y3qXucwu09WFq2FICsEVkUbi9kq/NNskdmm6FGBL+Y9P5meMGiEJPe399D6hASYAJQqLU1FyIghUH3ayOkLC1byp8+/RP1Wr2ElxDkLGkML7h0nCX7wiLEKLoeRD2yW8HpdBIbG0tNTQ0xMTH+Hk6rSB8YIURbuvIvV1Kv1ROpRrLtzm3+Ho5oQ0bNi3HbqOlxMGrp9VtmYAJQqLQ1F0L4X+H2QjO81Gv1FG4vlBmYILd5TTmKqjCkk6VZWNn5o4sL+neD4kqAoA0xLSGN7IQQYaM6fwnVBQW+nysooDp/SQePqH0Vbi/kozV7eajhSbbduY3skdksLVtK4fZCALYUOdi8ptzPoxStpagKm9c4+HbnUa/w4v55OjgxIJaYjAHoIbiBoyeZgQlwNpsNVVVJTU1t9pzdbkfTNNLS0vwwMiGCkEXl8GJ3h9KEaY1LiKsLCji8OJ/4GaGzaqNweyFLy5byUO8n0bfEsyXRQVZmY02M9nEc+pZ4xkwO0Z3+QtjoTPfPbP0aB2MGxDKaxvAyZrLVfD7USYAJcKqqYrPZALxCjN1ux2azSXgRohWM0OIZYjzDi2eoCXaarp0q2L2FLYnuixtAVmaWV3gJl4tdqDF+bpvXOPj4gwq0Bj3sfp5SxBsEPMNKampqs2MhROsYoUWJjESvrw+58OKL8S90NUIJy4tdqHpxug2tQUeNUHhoSWj8g7al128JMEHCCC0WiwWXyyXhRYjz9NVll6PX16NERjL0s0/9PZwOEYoXu3AWqqG0pddvKeINEqmpqWZ4sVgsEl6EOA/VBQVmeNHr609b2BtKthQ5zPCiNehsKXL4e0iilQrKCpoUYLtrXh5akoYy+jCb1zjC6ucqASZI2O12M7y4XC7sdru/hyREUPKseRn62afEz8jh8OL8kA4xW4oc1BRXcsOIHjy0JI0xk61eFztnyT5qTi27FYFLVVSWli2l4M8rvQp2C7cX8mLEc2EXYqSINwicrgYGkJkYIVrBV8Gur8LeUGL8S/2GET3oVunEWbLPqwC06z4n3SqdxGQM8PNIxdkY/Xs+WrOF0aNHMzrzBnO1mWfBdqgvnzZIgAlwvgp2jf9KiBGilVyaz4Jd89gVejv36pq7NmJYptXs0gruVSye4SWUG56FEneIKWRp2XO89Jfnm20PEQo1MC0lRbwBTvrACCHakhliTu2bI+ElOIXy9hCylUCIOFM4kZkXIc5u06oVKKpKytTbmj1X+vYb6JrGNTff7oeR+Ue47lwcSmR7CDcp4g0im1atoPTtN3w+V/r2G2xataKDRyRE4FNUlU0rm//dKX37DTatdIebcOJr52IR2F4o3sXikt0AXjUv2+7cxqiYW722hwgnMgMTRIwPYsDrX5PGB/E1t4TPvyKFaCnj74rn3x3PvzO+ZmZClbNkH+XlC+mUEcfQ9Ce8amJi0vvjcOSj6xpJSQ/7eaTCk0VVyCvexTbnSrY63zRrXhaX7GbDRyMZN9a9PQQQVjMxEmCCiHwQC3FuPP/ufPS3v+JqaAi7vzNGWOmUEce3yjKiHbFY0917PzmLK9nPq3yrLCPJOtO/AxXNzEgfBMCST4q55uJbzfCSV7yL3IzBzEjPpHB7DzQ99IrQz0SKeIOQEVosERFh+UEsxLlaePsUXA0NWCIimLniHX8Pp0PVFFeiqIo501LuWEiSdSZWaw5flfyPGV6s1tDZ0DLUGKElyqJS59JOhZdB/h5Wm5NOvCEsZeptZnixRERIeBGiBUrffsP8O+NqaDhtPVmoivVYbWS15pBknUm5YyHrbcMkvASJGemDzPASZVFDMry0RqsCzLx58xg9ejTdunWjZ8+eTJkyhZ07d3qdM27cOBRF8fqVleV9T27fvn1kZmbSpUsXevbsyaOPPkpDQ4PXORs2bODKK68kOjqa5ORkli9ffm6vMASF+wexEK3leat15op3uOaW230W9oYTqzUHRYlC1+tQlCgJLwGovHwRDke+eby4ZLcZXupcGs+89aofR+d/raqBsdvtZGdnM3r0aBoaGvjNb37DhAkT+OKLL7jgggvM8x588EGeffZZ87hLly7m/7tcLjIzM0lMTGTTpk0cOHCAu+66i8jISH73u98B4HA4yMzMJCsrixUrVlBSUsIDDzxA7969mThx4vm+5qDWtObFOAZkJkYIH3zVifmqJwsnLxTvwlnzEdfF1ZkhxuHIx2rNYXHJblyazqyMwf4eZthTFJVyx0IA1pRPMmteJietZdE/v2TZx5l07747bGdiWhVg1q5d63W8fPlyevbsydatW7n++uvNx7t06UJiYqLPr/GPf/yDL774gn/+85/06tWLkSNH8txzzzFnzhyefvppoqKiKCwsxGq18oc//AGAYcOG8eGHH/LCCy+EdYCRD2IhWk/XNJ91YsaxroVX4SOAs+Yjln0cB1c9zVO/vNOsiXl164Us+ziOXAkvAcGYFVv0zy95Z2+SGV7KHQt5ePxMkpIGk1e8CyAsQ8x5rUKqqakBIC4uzuvxFStW8Nprr5GYmMjkyZN58sknzVmY0tJSLrvsMnr16mWeP3HiRB566CE+//xzrrjiCkpLSxk/frzX15w4cSIzZ8487Vhqa2upra01j51O5/m8tIAkH8RCtN6ZmtSFY+B3OPK5Lm4hXPU0yz6OO/Uv+BwzvNx71dGwvBgGKqs1h9gLX2PKxUVcqj5CuaPOrFeacWrXAFeY7H3U1DkHGE3TmDlzJtdeey2XXnqp+fivfvUrBgwYQJ8+ffj000+ZM2cOO3fu5G9/+xsAVVVVXuEFMI+rqqrOeI7T6eTkyZN07ty52XjmzZvHM888c64vJyjIB7EQ4nzpukaSdSbpN9xJ9+7uVS1L1u+hzuUOL3deedTfQxRNPHPzHay3zfVZrxTOYfOcA0x2djY7duzgww8/9Hr817/+tfn/l112Gb179yY9PZ29e/dy8cUXn/tIz+Kxxx4jNzfXPHY6nfTr16/dvp8/VecvAYvqc9fc6oICcGkk5Ez3w8iECByyj5hvnk3qZqQPOhVe3IWhT/3yTj+OTJyOu8Fg83qlcHdOy6inT5/Oe++9h81mo2/fvmc8d+zYsQDs2bMHgMTERA4ePOh1jnFs1M2c7pyYmBifsy8A0dHRxMTEeP0KWRaVw4vz3WHFQ3VBAYcX54NFVscLoaoqNpsNu93u9bixw7saZlsI+NJ0VYvRrl74X01xJc6SfV49e25I+9Jc/v7lP//H30P0u1bNwOi6Tk5ODqtXr2bDhg1YrWfftrusrAyA3r17A5CSksLcuXM5dOgQPXv2BKC4uJiYmBiGDx9unvP+++97fZ3i4mJSUlJaM9yQZcy8HF6cbx4b4SV+Ro7PmRkhwo0x82Kz2cxjI7ykpaWF/Waoi0t2s+STpYwbm8Dynz9uNkkD98xM4fZCNF1j2kj5PPEHRVUo37uYI8pqrx49PcpvomZvJfuTl9HJERvWMzGtCjDZ2dm8/vrrvPvuu3Tr1s2sWYmNjaVz587s3buX119/nZ/85Cf06NGDTz/9lFmzZnH99ddz+eWXAzBhwgSGDx/OnXfeyfz586mqquKJJ54gOzub6OhoALKysliyZAmzZ8/mvvvuY/369axcuZKioqI2fvnByzPEHHmxEL2+XsKLEE14hpiNGzficrkkvNDY0XXc2AS2Ot+kcHsPZqS7+3U13XNH+EdMen+iiaHHnp/TQ78JrI3bQSRlzCDWOgA9zLYOaKpVWwkoiuLz8WXLlnHPPffw9ddfc8cdd7Bjxw5OnDhBv379+PnPf84TTzzhdUunsrKShx56iA0bNnDBBRdw99138/zzzxMR0ZinNmzYwKxZs/jiiy/o27cvTz75JPfcc0+LX1gobyXg6avLLkevr0eJjGToZ5/6ezhCBKTnnnsOl8uFxWLhySef9Pdw/O6F4l1YVMWcaTF2N84akcU9q+d6bRgo/MvccPPU7uExHh2VQ1VLr9+yF1IQM24bKZGRMgMjxGkYt40sFovMwJyGEWIi1UjqtXoJLwHmm8c/BJcOFoW+c6/z93DaneyFFOI8a16GfvYp8TNyfBb2ChHOPGtennzySdLS0nwW9oa7rBFZZniJVCMlvPjR5jXlbClymMfOkn1meMGl80VBmf8GF2DOq5Gd8A9fBbu+CnuFCGe+CnZ9FfYKKHhpFZcdTOOz/jbqtXoKtxeaIWZLkQNd0xkzOcnPowwPiqqweY07wAzpZMFZXElMxgB2/uji2D8qGbbvOM6SfSF/G6klJMAEI5fm83aReewK78IuIQCzz0vTkGIca9K5GnDfPvro4G7GfJ3J9JHT2dp3HUvLlgIw6puJbF7jYMzks684FW1jdKb7z/rYPypxdraY4cX4OcScCjVA2IcYqYERQoSMBY4DWBSF3IHN92LLq6jCpes8au3th5EFJs8CXs+wsrXvOkrfc4eaMZOt5kVVdJwvC8r4Zucx9rh0tAbd6+fgLNmHrunEZgzw8yjbR0uv3zIDE+xs80C1QOrs5s/Z54PmgrTHOn5cQviBRVGY73C3d/AMMXkVVcx3VDHb6nuT2XCl6Vpjwe4I92Ob1zhQI4YwpmEwjK5mdOYN/h1kmBo2bSQbptvQGnTUCMUrRIb7zItBiniDnWoB21x3WPFkn+9+XLX4Z1xC+EHuwERmWxOZ76gir8IdZDzDi6+ZmXA2beQ0r4Ld0ZlW1AjFvGhm3/+ffhxdeCkoK6Bwe6F5vKXIYf4ctAadpX/+qx9HF5hkBibYGTMvtrmNx0Z4SXvc98yMECHMCCnzHVUsrDhIna5LeGmhLwvKSLYo7AG0Bp0tRY6wuW3hb6qi+qw92tp3HZvf282YLZlsSXTI7TwPEmBCgWeI2bgAXHUSXkRYyx2YaIaXqNPUxAhvW4ocHNt5jGGdLYz1KBwF79Uwon0YM2Gl7+3G9fVgM7wsLVtK9k+zGfWN1fx5SIhxkwATKlJnN4YXS5SEFxHW8iqqzPBSp+vkVVRJiDmDLUWOZqtchmQMgMlWr9UwUnvRvrJGZKFtWcVHFPHnozbqD9c3q1HStZBcd3NOJMCECvv8xvDiqnMfS4gRYahpzYtxDEiIOQ13nxfv1UbO4kp6WxR6d7bg7N+NvhJeOsS0B27mpb/M89lUUGZevEkRbyjwrHl5str9X1+FvUKEOF8Fu74Ke4W3MZOTmq9yOdX5FYvC8Gkj/Te4EPdC8S4Wl+w2jwu3F5rhpV6r5+6/zfXj6AKbzMAEO18Fu74Ke4UIA67TFOwax67QbHvV5r4seZbaAU7iK6eAS/fq/Opw5KPrGklJD/t3kCHCoirkFe8CICq+xOzLU3c4nSWfvMg23LuFy/YOzUmACXaay3fBrnGsuTp+TEL4yZma1Mnto5Zxluyjdq+TI8mric0YQI/ym8zOr0eS3qXcsZAk60z/DjKEzEgfBMCST14kOqHYDC95xbvIzXiIqPjB5uokCTHeJMAEuzM1qZOZFxEOpJljm3GW7MNZXElSxgxirQModyyEJOiRcRPlexdzRFlNknUmVmuOv4caUmakD2JrTXdKyzPIWzmAOtcucjMGnwo37oCj6bL1RVMSYEKItFEXYclo5gjeIcbz9qpoEV3TzdVGMbhDSrljIQ61AD25jj7avRJe2skrv3iSwY9/QJ1LI8qimjMzIDMvpyNFvCHEaKPetFDRKGy0KIqfRiZEO0qd3bxwXZo5npPYJkulrdYcFCUKXa9DUaIYNv4JP44utJSXL8LhyDePF5fsNsNLnUvjmbde9ePogoPMwIQQzw6kxrG0URdhQZo5trkXinfhrPmI6+LqzBDjcORjteawuGQ3Lk1nVsZgfw8zaCmK6r5FB6wpn3Sq5mUwk5PWsuifX7Ls40y6d9/tNRMjvEmACTHSRl2ELWnm2KacNR+x7OM4uOppnvrlnTgc+ZQ7FvLq1gtZ9nEcuRJezotxK27RP7/knb1JZngpdyzk4fEzSUoabK5OkhDjmwSYECRt1EVYkmaObcbhyOe6uIVw1dMs+zju1ExAjhle7r3qqFxU24DVmkPsha8x5eIiLlUfodxRZxZJzzjVlsclnXdPSwJMCJI26iLsNK15MY5BQsw50HWNJOtM0m+4k+7dd5NXvIsl6/dQ53KHlzuvPOrvIYaMZ26+g/W2uWadkWeRtITEM5MAE2KkjboIO9LMsc15NqmL6FFM54S9nKy+gSiLylO/vNN8rnB7IZquMW3kNH8MM2jVFFeiqAox6f1PNQZsrDP6quR/6K3dKbt+t4CsQgoh0kZdhKUzNXNMe1yaOZ6nzY7viIj/B50T1lPn0sy294XbC1lathRVkctIaymqgrO4kq9K/sdsDHhD2pdcpN/Lt8oyDqh/8fcQg4LMwIQQaaMuwkV1/hKwqCRMm9asSV11QQG4NBJypsvMy3laXLKbDR+NZNxY2MqbjEuKI68YtjlXstX5ZuNOyaJVYtL7s59X+VZZxkW6u7eOs2QfXYtTuSgDvlWWEe2IlZ47ZyEBJoR4Nqmz2WyoqkpqairgffvIbrejaRppaWkdPkYh2oRF5fBidw+NhGmNty+qCwo4vDif+BnywX++FpfsNpf2zkjPpHB7D5aWLSV22FtsdTYwKuZWCS/nIcrajYsc99K1OJVv1n8ILncTwb7pTxDtiEWXzrtnJXN/IUpVVWw2G3a73etxu91uhhshglXCtGnEz8jh8OJ894wL3uHFM9SIc+PSdI929u5usJFqJBoNqEQwouvNfh5h8Nm8ppwtRQ7AXWc0NP0Jc9dvXYGvfmgA3KuTZLPMs5MZmBBlzLzYbDbz2AgvaWlp5vNCBCsjpBxenM+RFwvR6+slvLShpk3qCl5axWUH0/isv416rZ7OPdcD7nO2FDnQNZ0xk5P8MNLgoagKm9e4A8zoTCvOkn1meFF06PbNcT+PMLhIgAlhniFm48aNuFwuCS8ipCRMm2aGFyUyUsJLOyncXshHB3cz5utMpo+czta+68wdkkd9M5HNaxyMmWz18ygD3+hM95/R5jUOuu5z0q3SyfEBMazffoQbRvSgW6UTZ8k+r+0cxOnJfYQQl5qaisViweVyYbFYJLyIkFJdUGCGF72+3rydJNqOsdoo5aeDGDPZyuY1DkZ9M5HskdmUvrfbDC/GxVmc2ehMqxlWvqp1sX77EcZMtjLsoRHEZAzAWVzpnpkRZyUzMCHObreb4cXlcmG32yXEiJDQtObFOAZkJqYNabrWuNpohPuxzWscqBFDGNMwGEZXMzrzBv8OMsAVlBWgKqpZ9Nwn+UK27PqOnSc1UDW2XPQBo5lmzrzo0n23RSTAhLCmNS/GMSAhRgQ1XwW7njUxnsfi/DRtUjc4SqWms8rOkxpqhMJD9/+n+ZyzZB+6pksTtiZURTVvuWWNyGJXnWaGFzQVtsbDSPe5cvuo5STAhChfBbu+CnuFCEouzWfBrnnskiWo7WX/nmMMjbYAsPOkxpYih1mQ6iyuJEbCSzPGzMvSsqVoH8ehb4lHGX2YFyOe46GGJ9G3xLMl0SG34VpJAkyIMvq8NA0pxrGmyQe8CD6bVq1AUVVScqY3e6707TfQNY1rZOal3WwpcrD5VMHp0EonFw3uznqPgtSYjAEyg3AaWSOyzPCytf9atkR8cOrW3C1sSXR4rU4SLSMBJkQZTerMD/ypt5nPGSHG/MC/+Xa/jFGI1lJUlU0rVwB4vadL336DTStXcM0t8l5uL1uKHGbB7jBjCXBxJT/rHolyajVNXwkvXl4o3oVFVcxeOlf1vIrC/oVsuegDVCI4echdO2SEFql9aZ1WrUKaN28eo0ePplu3bvTs2ZMpU6awc+dO8/mjR4+Sk5PDkCFD6Ny5M/3792fGjBnU1NR4fR1FUZr9evPNN73O2bBhA1deeSXR0dEkJyezfPnyc3+VYcz4wC99+w2vx40PfEUa2okgkjL1Nq655Xav97RnePEMNaJtufu8NK42OmxdzeGL/46ig67A8b7dzHMdjnzKyxf5a6gBw6Iq5BXvMveP2tb/H2Z40Whg+/erzHNHZ1qlj04rtWoGxm63k52dzejRo2loaOA3v/kNEyZM4IsvvuCCCy5g//797N+/n9///vcMHz6cyspKsrKy2L9/P2+99ZbX11q2bBmTJk0yjy+88ELz/x0OB5mZmWRlZbFixQpKSkp44IEH6N27NxMnTjy/VxxmjA90z3+1yge+CGae7+mP/vZXXA0N8l7uAE0vrnWO4xy5+G8AxO/9GUO7uC8nDke+uUFhuDNmXvKKd5n7R42KufXU/lJlbHW+SeH2HrIlwzlSdP3cd/irrq6mZ8+e2O12rr/+ep/nrFq1ijvuuIMTJ04QEeF+gyuKwurVq5kyZYrP3zNnzhyKiorYsWOH+ditt97KsWPHWLt2bYvG5nQ6iY2NpaamhpiYmNa9sBBkhBZLRIR84IuQsPD2KbgaGrBERDBzxTv+Hk5YMQp2v8+wmxsSdi1ONY+TrDNlI0IP96yey1bnmzQcnsDJ6hvMLRqMHjuyKaa3ll6/z+v+gXFrKC4u7oznxMTEmOHFkJ2dTXx8PGPGjOHll1/GM0eVlpYyfvx4r/MnTpxIaWnpab9PbW0tTqfT65dolDL1NjO8WCIiJLyIoFb69hvme9nV0NDsFqloP56rjYamP0GSdSbfKsvYlfGg1+7K4ay8fBEOR755PMba3QwvkarOT63vA+7C3uyR2WiyceM5OecAo2kaM2fO5Nprr+XSSy/1ec7hw4d57rnn+PWvf+31+LPPPsvKlSspLi5m6tSpTJs2jfz8xh92VVUVvXr18vo9vXr1wul0cvLkSZ/fa968ecTGxpq/+vXrd64vLSTJB74IFZ63QGeueKdZTYxoX7qmm6uNXijexZrySShKFLpSj6JH0lu7E3DvZv1C8S4/j9Y/FEWl3LHQDDENRzLM8FKvKfxlW+M/+rNGZDXrtSNa5pxXIWVnZ7Njxw4+/PBDn887nU4yMzMZPnw4Tz/9tNdzTz75pPn/V1xxBSdOnGDBggXMmDHjXIfDY489Rm5urtf3lxDj1rTmxTgGZCZGBBVf9Vu+6rxE+/FsUmcUqU65OI2fJdvQqeNo8t95pWQSecW7yG2yIWS4MGagyh0LeXXrhSz7OI57rzrKdXFP8+HRp1n2cRzdu+82a2TEuTmnADN9+nTee+89Nm7cSN++fZs9f/z4cSZNmkS3bt1YvXo1kZGRZ/x6Y8eO5bnnnqO2tpbo6GgSExM5ePCg1zkHDx4kJiaGzp07+/wa0dHRREdHn8vLCWnygS9Cia5pPuu3jGNd+ht1qMlJaym/+Eve2ZtJUtIsJietZdE/v+SdvUlmnUe4slpzzPAyJXkt18W9T5J1Juk33En37rvJOzU7Fc5/RuerVQFG13VycnJYvXo1GzZswGpt3nDH6XQyceJEoqOj+fvf/06nTp3O+nXLysro3r27GUBSUlJ4//33vc4pLi4mJSWlNcMVyAe+CH42mw1VVUlNTW3Ws8hut5tNGyWIdyxjtdHD42fyXY9tLPmkhCXrM6hzJTHl4iImJ5UD7kJVTdfC5jZJTXEliqoQk96fmNixTEnOZ3LS+yhKFD3Kb6JmTyUzMtyhxSV9X85LqwJMdnY2r7/+Ou+++y7dunWjqqoKgNjYWDp37ozT6WTChAn88MMPvPbaa17FtAkJCVgsFtasWcPBgwe5+uqr6dSpE8XFxfzud7/jkUceMb9PVlYWS5YsYfbs2dx3332sX7+elStXUlRU1IYvPTx4fuBX5y8Bi2q2W/f8wK8uKACXRoKPDqdC+JOqqj63v/DcLkN0PF3XzNVGKc5Cth1fSS0Q9V0GD48fhq5rXqtswoWiKjiLKwGYkryOcos7vOh6HeV7F5N0sbtUQmZezl+rAsyLL74IwLhx47weX7ZsGffccw/btm3jo48+AiA5OdnrHIfDwcCBA4mMjGTp0qXMmjULXddJTk4mLy+PBx980DzXarVSVFTErFmzWLRoEX379uWll16SHjDny6L63OjOc2M8IQKNrz28fO31JTpWUtLD5v/XHU6ntnoX0QnF1AJryh8iKr4kLJcIG1splO9dzBFlNUnWmfQov8l9nLyaWOsAYpDP2rZwXn1gApn0gfGt6S6+vnb1FSIQGaHFYrHgcrkkvASIxSW7zYJdI7TomgVFdYVVeNm8phxFVRidaTVvr/XY83PiK6eAS8fZvxudb/w/s8lfuC81P5OWXr9lL6QwY4SUw4vzOfJiIXp9vYQXERRSU1PZuHEjLpcLi8Ui4SUAeIaXGemD2LzGwuhvy9ly0QfomoW6w+nmuVuKHKe2IwjNdvmKqpgbMvYY5r69FllyBbh0XLrOiQGxDD8VWnTp+9ImZCOcMJQwbRpKZCR6fT1KZKSEFxEU7Ha7GV5cLhd2u93fQwp7Lk33Wm308aGPGbVvEqO/vRFFdbG1ZiXQuBGkoir+HG67cu9lZGXzGgdHvvwZPcpvMsOLRVEY0skCuFcned5+E+dOZmDCUHVBgRle9Pp6qgsKJMSIgNa05sU4BmQmxo9mefR5KdxeyIsRS3lo9JOM2jKJq3pdxYvHn6Pgz4PQt8R7bQQZSgrKClAVlawRWebrO/aPSpydLXx50sWuyw9y58AUs7A3RnbsbjMSYMLM6WpgAAkxIiD5Ktj1Vdgr/Md7T59b+LK2jG8+VcmqfwFdU1FGH2Z05g2AeysCXdO9GuIFM1VRWVq2FHB31R3SyWKGly/r6oi66jtiRrhDi4SYtiUBJoz4Ktj1rInxPBbCnxY4DmBRFHIHJpp9XoyQkldRhUvXefTUsSa9jPxO0zWvgt2LhsTRbd9xAHbVN8Cow4D3PkqhwnjNRohJ3Xk135wKLxY9glHfTIQRjaFFl94vbUYCTDhxaY3hxTYPVAukzm4MLa5TFwL7fNBckPaY/8YqwppFUZjvcPeZyvXo85JXUcV8RxWzrYmAzLwEiqZN6nb+6OLYSRfDOrvrPi789kacRxrDSyjMQLxQvAuLqjAjfZAZYkrf243r68Fs7lcEw+KYHv0rs7B3dKY1JF53IJEAE0a8mtSpFrDNdf+/Z4ixz3c/nvZ4xw9QiFNyB7oDihliBiZ6hRfjeRF4jILdMZOtHO/0FxL2Oum+8Wc4FcUrvDgc+e5meEFa0GrsAwXupnSjvplohpetfUr4dewqRqe7a2I8Q4xoOxJgwlXqbPd/PUKMV3gxnhfCTzxDzMKKg9TpuoSXAOcZXtz9ULqzX3kFgO57f8bOH12MpnEbgiTrTL+O93wYK6+MEBNxoIzN/f7N1j4lKKqLqPgSYJAZWuTWUduTABPOPEPMxgXgqpPwIgJK7sBEM7xEnaqJEYHL3eelcbVRj/KbqNlbyZHk1QBcUHmHV3gJxmZu5eWLUBQVqzXHDDFLPnnR3YU4OoPbu+XSa9APXoW9MvPSPiTAhLvU2Y3hxRIl4UUElLyKKjO81Ok6eRVVEmICmGeTOqNgd23Sf1J7dATXJT+Noq1BdzSY4WVxyW5cmu61HDvQKYpKuWMh4O7pEhVf4g4v1RloR1K5a+o+M5h5hhjR9qSRXbizz28ML64697EQAcCz5mXfuBHMtiYy31FFXkWVv4cmzsJztVGXiy9k2cdxrCn/CbragKJF0KP8JrOLryXImttZrTkkWWdS7liIw5FP6d5qM7w06BGsKZ8EuENL9shsNOm6225kBiacNa15MY5BZmKEX/kq2PVV2CsCk67pZsHuDOC77/7NigoLjorhzBj4BQv2rOMv5X3MLr6F2wvRdK3ZaqZAU1NciaIqWNPdMyyL/vkl9r2ZTLm4iKwrPuPveyd6FfbKzEv7kgATrnwV7Poq7BXCD1ynKdg1jl2huQdtyPBsUudw5HNd3EK+qJvAF/UVZH/yM36s7sO9Vx01w4vRBC/QKapiNqNbwyTe2ZvElIuL+FnSeroW/39kZfSi08UXeoUY0X4kwIQrzWWGF8+mYWZo0VyAR9Mwa28/DlaEBY/eRM3ebx69iWTmJXh4Fuz+9YYchv3hv4mM/wcA18X9nXkb1vN65ccBv2u1507T4O6oezBpJ1Mu3sZkawm60sD3GXb6pj/BjFO/xyWrjtqdBJhw5dGkzqtpmEeIado0TIh21aQ3kUl6EwUtXde8CnZPVt9A6pFBcMFucr/uhoZ3eAnUHas9d5oenWllP69yg7KMHnt+Tvw/X2L/lWv5Nn4Z0Y5Yr9VJon1JgBHSNEwEBulNFHKMJnVGwW5uxmBSfhzO5jXuVUdlfdd7hRejh0ygMWZeNq9xUBvxGrWRy93hpfwmdAXG3DIXhyPRa3WSaH8SYAQgTcNEgJDeRCHHM7wYNS+b++3mzurJpNUO4Z7Vc8mO+pVXA7xA2fDR107T5bv+Tg/LFOL33YRL0bHoCs6SfWZhry6rjjqMLKMWptyBiWa/DWkaJvwmdXbjsn7pTRT0XJruFV6Wli0l5aeDuGhId8bXDeNX9pRm4cV5arWPvxk7TRduLwTcMzE9yqcwdN8Uvvixng9u/YKYjAE4iyvdIcaaE7RbIwQjmYERJmkaJgKCr95EEmKCltGkznO1UdaILBgB/5y1gZ6XvIdSD1v7DmZIyU+8Nnz0135JjRs1NtlpelMKQ6MtfPFjPbt/hDHfTCQm0723k7E6STZs7DgSYATQvO+GcQzSb0N0IOlNFLI0XWtWsLvzpAb1CurQdxj70c9xlnuHF3/tl+S9UWPjTtPJ1UP5JmYXG8ZWm7e9ANnvyE8kwAhpGiYCg/QmCmmeTeq8N338A3Oet1Lb4wvuSXqXvum/9wova8on4dqzq8O2GygvX8TkJBUyJpkhJuVH907Tf+m3hrK+JWz/eZl5vuw07T8SYITZNOwKx5fYK3eSmpoKeDcNs9vtaJpGWlqaP4cqQlB1/hKwqCQMczUr2K0uKABXFxLSHjd7E4ng1nTHamfJProdG8pfj1nh4iIm/3MIqA1meDEKgNub0WVXSXLvdTQ5CTPEXK98jdJvN59cVIyOTuH2Qq/CXpl58Q8JMMJsGmav3InNZgPwCjF2ux2bzSbhRbQPi8rhxfkwI4eEaY3/Sq8uKODw4nziZ+TIzEsI8dyx2ijYnZExmBhqWfJJMY7KQTw8cI9XeOmI7QaMLrs9Mm6CJMwQUxDv4JOEYuqr0/kg9UHWOSNlp+kAIQFGmIzQ4hliPMOL8bwQbckILYcX55vHnuHFM9SI4Gc0qfPc8DEmvT+THfmUVOyhPKKS5f98mPqT5eROHtxsu4G2bnbnq8uuEWJmrl9HZEIl9dXp/Hg4gzXlg80mdbLTtP9JgBFePEPMxo0bcblcEl5Eu/MMMUdeLESvr5fwEuI8N3w0al4W3jCTn/ytAqXrHq77bgjXlO2ipOpTlja4w8uobyaat5/Ot1eM0eNllDqxWR2Ls7iSvzCUL+P3Mqw7PDxqHTu07GaFvbLTtH9JgBHNpKammuHFYrFIeBEdImHaNDO8KJGREl5CnBE8mhbsnqzexeaE9dCviMH7JjOsOoaFlz1NzDdXmuFlSCeLOXvTWsYS6ah49dSsDoyZ7A4x/y4/wscXreNiRvMy9UyJreNnSXvR9TomJ631Kuw1QozwH2lkJ5qx2+1meHG5XNjtdn8PSYSB6oICM7zo9fXuAl4R8oz9kjxrXr787z9Q1nc9r/Vcw5cnXfT44V+8ur2IfSO6eoUXY/bmmbde5YVTweJMyssX4az5iLziXdQdTid7ZDZLy5ayte86TiQ5UT93cvvGFE5cuIspFxdxj96Fq7R1JFlnnqqJWUtuxmDZqDFAyAyM8NK05sU4BmQmRrSbpjUvxjEgMzEhLinpYZ/bDWg0UNZ3PaO+mQh1CpaEL/jrXiuVJ1aSMjaBh9P/A4cjn0X//JJ39maSm6GcttDXc4XRdXFPw1VPn/p+6TwQU8PSsqXoCRby9ufRc/j7TEheTZJ1Jj3Kb/JZ2Ct7HQUGCTDC5Ktg11dhrxDna9OqFSiqSsrU25qFl9K330BPiGXQjBwJMWHC13YDZs2L7uDLHRO5e+Q/4OIi3q/pxOfO13DlR/H9dxHE/fAL/nRxD451eZZlG7sTH3c13+6ai4rKhKpJAPTvsxvLx73ZmXQtzsu2cF3c0zit/4+8YohgOJ2GWHjEMY3hnSL5yuXi8I6biDn5U6weXXaNECN7HQUOCTDCZPR5SU1N9brAGKFF09x/cUvffgNd07jm5tv9OVwRxBRVZdPKFQAkuzSv8LJp5QquueX2xtDikgtGqPO13YBXwe6xbSR8dhN3W9+F2B94vzqD7bXfMGZ/JoP6dWb43u/5Z+djDOnUjXc+iyU59jsO9llHXFUPoo4mYzt2EBfdebAcXnUmolvH85NBz/NuRR5qDzuj9qfT1zmYf8Z8yY/XDDG/N3h32ZWZl8AiAUaYPPu8eF5gPEOM5wVGiHOVMvU2APd77JbbSZl6m9d7y3heZl7Ci7HdgGd4GdLJgvOznuy87BBf/Tian/b/Le+tn09pDzv0K4KvM9nfcw/jd8xm+JVr4eIi3tmbSUqkystD8hm/7xaSv72Wb5L38SpfcNfhybiqdV64OBK1h52UWoUxX2fy5SVdibxea1bYC9JlN1BJgBE+eV1gTh37usAIca4832Mf/e2vuBoa5L0V5ozalc37yputNkpP/w+SHPks+mc6DXoEEUdS2Tr4KQDGfJ3JV51dDN02iXuS3jVDTCfq+Gf/lXynOhm950ZQ++LqplMYt533nVFEJxQz+FAGJ5KcFHx7ktzD6WSPxCvESJfdwKXouh6SPx2n00lsbCw1NTXExMT4ezhBywgtlogIucCIdrHw9im4GhqwREQwc8U7/h6OCCBG8a2x2qixYHcwTudv+atzB7pm4deb52PRI/hpbAQWRWHX+Af4dcl8GvQIug55HEV18eC//8CwqCjW9/mMN75L4sp+K0ntM5yfjexCuWMhHx59mmUfx5GbMZio+JJ27forzqyl1+9WLaOeN28eo0ePplu3bvTs2ZMpU6awc+dOr3N+/PFHsrOz6dGjB127dmXq1KkcPHjQ65x9+/aRmZlJly5d6NmzJ48++igNDQ1e52zYsIErr7yS6OhokpOTWb58eWuGKtpIytTbzPBiiYiQ8CLaVOnbb5jvLVdDA6Vvv+HvIYkAEuuxVNozvETFl/BX5w5+NeAqbvn6fix6BIM6gUVROJz0Ln93uGdpOsUXo6guRn0ziWFRUSRc8h4nuu9kqnU/y4c+x9SPrqFH+U0kWWdyXdzT3HvVUVyaTtaILAkvQaBVAcZut5Odnc2///1viouLqa+vZ8KECZw4ccI8Z9asWaxZs4ZVq1Zht9vZv38/v/jFL8znXS4XmZmZ1NXVsWnTJl555RWWL1/Ob3/7W/Mch8NBZmYmaWlplJWVMXPmTB544AHWrVvXBi9ZtIZcYER78bwlOXPFO1xzy+1sWrlC3mOiGV3XiL3wGjO8GIW+40/8lh77L8Hau47hnSLZPOgllis/uGtg4u1EJpQwft8t3H4ok2GdLWy/4Cuu7VFP0icXs/NHFzEZA9wrjE6FmDuvPNphu16L83det5Cqq6vp2bMndrud66+/npqaGhISEnj99df55S9/CcBXX33FsGHDKC0t5eqrr+aDDz7gpz/9Kfv376dXr14AFBYWMmfOHKqrq4mKimLOnDkUFRWxY8cO83vdeuutHDt2jLVr17ZobHIL6fw1rXmRGhjRVk73XpL3mDgbcwuAU4W+N4zoQbdKJ99n2Cks3887ezNJjrVxsM867tuZQ9TRZPrHH+HjhgupHLuZrc43eajhSfQt8WadzflsSSDaXkuv3+dVxFtTUwNAXFwcAFu3bqW+vp7x48eb5wwdOpT+/fubAaa0tJTLLrvMDC8AEydO5KGHHuLzzz/niiuuoLS01OtrGOfMnDnzfIYrWsHzQpK8/wjVBQWknFoR4lnYW11QAC6NhJzp/hyuCAI2mw1VVUlNTXUvw/cIKXa73b2M/9SxrsnSaeFb00LfPlEqyuDu1Fm7EXvUPUtz8EgZKr/iqsTRkAh9J3/Otq0nGdH1Zq5O6oGmH2ZM4mj3fkzp/f38isS5OucAo2kaM2fO5Nprr+XSSy8FoKqqiqioKC688EKvc3v16kVVVZV5jmd4MZ43njvTOU6nk5MnT9K5c+dm46mtraW2ttY8djqd5/rSBHhdYDy7ohohRtc0rwZkQpyNqqqNDRE9egh5NlAEZOZFtEjT3ahjeJhnzIceb3L2FTxlroQ+dYtoZLsNTXSQcw4w2dnZ7Nixgw8//LAtx3PO5s2bxzPPPOPvYYQMzyZ1njsFgzvENO2eKsTZ+Orq7Kv7sxBCtMQ5BZjp06fz3nvvsXHjRvr27Ws+npiYSF1dHceOHfOahTl48CCJiYnmOZs3b/b6esYqJc9zmq5cOnjwIDExMT5nXwAee+wxcnNzzWOn00m/fv3O5eUJHzxDjLFjsIQX0VqeIcbY8VzCixDiXLRqFZKu60yfPp3Vq1ezfv16rFbv7oSjRo0iMjKSkpIS87GdO3eyb98+UlJSAEhJSeGzzz7j0KFD5jnFxcXExMQwfPhw8xzPr2GcY3wNX6Kjo4mJifH6JdpWwrRp5k7BSmSkhBdxTlJTU82dzi0Wi4QXIcQ5aVWAyc7O5rXXXuP111+nW7duVFVVUVVVxcmTJwGIjY3l/vvvJzc3F5vNxtatW7n33ntJSUnh6quvBmDChAkMHz6cO++8k+3bt7Nu3TqeeOIJsrOziY6OBiArK4vy8nJmz57NV199RUFBAStXrmTWrFlt/PJFa1QXFJjhRa+vdxfwCtFKdrvdDC8ulwu73e7vIQkhgpHeCoDPX8uWLTPPOXnypD5t2jS9e/fuepcuXfSf//zn+oEDB7y+TkVFhX7jjTfqnTt31uPj4/X//u//1uvr673Osdls+siRI/WoqCg9KSnJ63u0RE1NjQ7oNTU1rfp9wrdDS5fqXwwZqh9autTnsRAtsWHDBv2pp57SN2zY4PNYCCFaev2WrQTEWXkV7A6rAdUCqbObF/La54PmgrTH/D1kESAWOA5gURRyByY2K9jNq6jCpeuM2bdLCnmFEKYO6QMjwoRL8w4ptrkAJEybbT5vPp7WdPmiCGcWRWG+w90e4QpN8wov8x1VzLYmmqFFk94vQohWkAAjzsqrSV3qqdDiGWI8w4vxvBBA7kD3ykJ3WBlG7sBEr/BiPC8zL0KI1pIAI1rPM8RsXACuOgkv4rQ8Q8zCioPU6bpXeBFCiHPRqlVIQphSZ4Mlyh1eLFESXsQZ5Q5MJEpRqNN1ok7VxAghxPmQACPOjX1+Y3hx1bmPhTiNvIoqM7zU6Tp5FVX+HpIQIshJgBGt51nz8mS1+7+2uRJiRCPbPPP94Fnzsm/cCN52vo1mmychRghxXqQGRrTOqfDy4cgZbB5wF7nQrLA3b8BduHSdR629/TZM4WeqBWxz+dex48yPmdpY82Kfz7WfLIYrZjD11OokuZ0khDgXEmBE62guSHuczQPuMpfH5g5MNENM6dEa5mvuf22LMHbq/XCtbS5vXwHXDnzOa+bu2tTZzD7VB0YIIc6FNLIT56zpclhfy2NFmDNCi1ErJavVhBBn0dLrtwQYcV6M0GIUZ0p4Ec08l9BY8P1ktb9HI4QIcC29fksRrzgvsjxWnJGsVhNCtBMJMOK8yPJYcVqyWk0I0Y6kiFecM+P20c2c5JdKHVsHDvEu7AXsdjvaqT1wRGirzl8CFtV7z6xTNS/VBQXg6kKCEWJAamGEEOdFAow4J54Fu6Mqd57aTRhmWxtDTOPjEl7CgkXl8OJ8ABKGubzCi7FruRlaNJcfByqECAUSYMQ5cXkW7J6abfEMMeWOCpx2m7n7sAh9CdOmAbhDzKndyz3Di/G8zLwIIdqCrEISbcZut2Oz2bBYLLhcLgkvYcoILUpkJHp9vXd4EUKIs5BVSKLDpaammuHFYrFIeAlTCdOmmeFFiYyU8CKEaBcSYESbsdvtZnhxuVzY7XZ/D0n4QXVBgRle9Pp6dwGvEEK0MamBEW3CuH1k3DYyjgGZiQlhm1atQFFVUqbeBuBV87KnV3e+/+gjMAp7ZSZGCNGGJMCI8+YZXiIPfUPp22+QeuqC5hliSt9+A13TuObm2/05XNGGFFVl08oVACQf/M4rvGxauYJrbrmd+IGDG1cnSYgRQrQRCTDivBl9XoyQYlzQjBCjaZr5+DW3SHgJJcbMy6aVK/ihfzKXNwkvxvMAuDQ/jVIIEYokwIjz5tnnxfOCBu4Q4xlevC5oIiR4/sw/21+Bq6Gh2c9aZl6EEG1NAoxoc54XtI/+9lefFzQRWlKm3mb+rC0REfKzFkK0O1mFJNpFytTbsEREyAUtTJS+/Yb5s3Y1NFD69hv+HpIQIsTJDIxoF74uaBJiQoPNZkNVVXN1mectwrr4Puz7/FPzFqL8zIUQ7UUCjGhzxgVtZP9kLrcONos6Aa/ltrg0EnKm+3Oo4hyoqmquLos6vN8rvBir0ZJ7xUuIEUK0Kwkwok15/mvcWFabPCMHbrnd53JbEXyMmRebzcbFCXHNwotn3x9dk5VHQoj2IQFGtCld05oV7HqGmO8/+ojDf18n++MEOc8QU3G0BpdrZ7PwIjMvQoj2JAFGtKmmTeo8dyjuHhnJhbK5X8hITU1l48aNsveVEMIvZBWSaHeyuV9okr2vhBD+JAFGtLvqR287/eZ+9vlgm+e/wYkWW+A4QF5FFeC9fcSTTz7J0XGTmF9+QEKMEKLDyC0k0a6qCwo4vKaM+EudJOQ8TPXnXRv3xbnke7DNhbTH/TxK0RIWRWG+o4qKigriNjQW7OZVVLFS78Qt1oHYbGsB2cBTCNH+JMCIduO5M7ERVhLSHocZOe4Qc6mThJzHIXW2v4cqWiB3YCIA8x1V3DJuEqmpV5NXUcV8RxWzrYnkDkzErtSiycojIUQHkAAj2o9La16wa5tLgiUKLo2CgddLeAkyniHmnQ3bqdN1M7yAzLwIITpOq2tgNm7cyOTJk+nTpw+KovDOO+94Pa8ois9fCxYsMM8ZOHBgs+eff/55r6/z6aef8h//8R906tSJfv36MX/+/HN7hcJvEnKme4eX1NlgiQJXHQkj6kj4vbSbD0a5AxOJUhTqdJ0oRTHDixBCdKRWB5gTJ04wYsQIli5d6vP5AwcOeP16+eWXURSFqVOnep337LPPep2Xk9PY1MzpdDJhwgQGDBjA1q1bWbBgAU8//TR/+tOfWjtcEUjs88FVZ4YY7BJKA55tXrOfU15FFXW6ziOVrzLD8bJZ2CuEEB2p1beQbrzxRm688cbTPp+Y6P2vsXfffZe0tDSSkpK8Hu/WrVuzcw0rVqygrq6Ol19+maioKC655BLKysrIy8vj17/+dWuHLAKBfT7Y5vLhyBlsHplNbuWr7gJeMG8j5VVU4dJ1HrX29uNAhRfV4vVzMmpe3na+zbUVf+ZfV8xgqsMdYGQmRgjRkdp1GfXBgwcpKiri/vvvb/bc888/T48ePbjiiitYsGABDQ0N5nOlpaVcf/31REVFmY9NnDiRnTt38t133/n8XrW1tTidTq9fIkCcCi+kPc7mkdnMd1SRN+Au9+oj21ywzzcvjBZF8fdohafU2ebP6V/vPtkYXj5ZDGmPc+1NzzHbmuj+mcpMjBCiA7VrEe8rr7xCt27d+MUvfuH1+IwZM7jyyiuJi4tj06ZNPPbYYxw4cIC8vDwAqqqqsFqtXr+nV69e5nPdu3dv9r3mzZvHM888006vRJwXzeW+CKbOJvfUQ/MdVWC9i9w0KD1aw3ytyqsYVASQUzNk19rm8o36IhFavfnzhMaZF5eu+22IQojwo+j6uX/qKIrC6tWrmTJlis/nhw4dSkZGBvn5+Wf8Oi+//DL/9V//xffff090dDQTJkzAarXyxz/+0Tzniy++4JJLLuGLL75g2LBhzb5GbW0ttbW15rHT6aRfv37U1NQQExNzbi9QtBtjxsUoBpXwEgSeS2isYXqy2t+jEUKEKKfTSWxs7Fmv3+12C+n//u//2LlzJw888MBZzx07diwNDQ1UVFQA7jqagwcPep1jHJ+ubiY6OpqYmBivXyJwyUqWICMF2EKIANNuAebPf/4zo0aNYsSIEWc9t6ysDFVV6dmzJwApKSls3LiR+vp685zi4mKGDBni8/aRCD7GSpYIdOp0vVn9hNGqXnS86vwlPrZ7cNcwVfd4hurvf2LWLgkhhL+0ugbm+++/Z8+ePeaxw+GgrKyMuLg4+vfvD7inf1atWsUf/vCHZr+/tLSUjz76iLS0NLp160ZpaSmzZs3ijjvuMMPJr371K5555hnuv/9+5syZw44dO1i0aBEvvPDCub5OEUCM20e3KD8St2Gtex8dj5UsnvvsCD+wqD63ezC2gYifkQOXXNFsFZkQQnSkVgeYjz/+2OvCkpvrLsu8++67Wb58OQBvvvkmuq5z2223Nfv90dHRvPnmmzz99NPU1tZitVqZNWuW+XUAYmNj+cc//kF2djajRo0iPj6e3/72t7KEOgT4aj1vs63lllMhpuk+O6LjGc0HDy/Oh8kjSfipd3jxak6oufw0SiFEuDuvIt5A1tIiINGxFjgOYGlS82LMuGwbOBSXDrOTekt4CQDGXlbGLuLNwosQQrSDll6/JcCIgPDcc8/hcrmwWCw8+eST/h6OOOWryy5Hr69HiYxk6Gef+ns4Qogw4PdVSEK0lN1ux+VyoSgKLpcLu93u9Xzp22+wadUKP40uPGxatYLSt733pqouKECvr2dP7x7sjOvqXdgrhBB+JgFG+JVx+yi5Zw+6frGF5J49sNlsZogpffsNNq1cgaLKW7U9KarKppWNIca4fbT/pkns6nkhXa8ey+HF+RJihBABo1078QpxJp6rjVJTUyntFc+mlStITp2IzWbj688/5aB9HdfccjspU5sXhIu2Y/z5blq5gh82b6bPu2vZf9Mkyip2m3/+1QMGNa5OkloYIYSfSYARfqNpmtdqI8+LaKeeF7Gvej9pEl46jOef/6cjB6F5hBfwCC0uzV9DFEIIkxTxioCz8PYpuBoasEREMHPFO/4eTkiz2Wyoquq16sv486/veRFjf3Gr9OMRQnQoKeIVQank0ZlmeHE1NHgVllYXFFCdv8SPows9qqo2qzkywsuPPXrz7Zc7/DxCIYTwTW4hiYBR+vYblO3bw6ADR0m59Q729OrOppXu1UfJB79r7AIr2owx8+JZc9QrdSJ7Dh0huWcPDtrXUdorXm7jCSECjgQYERCM1UbX3HK7GVaSZ+TALbezaeUKqg8cJUUaqbWL1NRUvv78U/YcOoIyfDTHDx1pVlgNSIgRQgQUCTAiIOia1my10eHF+XSPjGRQ9650vnqMhJd2lJQQx97qo+i6jsViaVZYrWtSuCuECCxSxCsC1leXXoLeoPnuAmuf796HJ+0x/wwuiMl2DkKIQCZFvCKoVRcUuMOLqqPX13s3ULPPd++ErFr8N8AgZlEU5juqyKtw7wBuhJej4yaxecBQLrYO9CrsFUKIQCS3kETAMbrAxs/IIeGS76nOX9TYQO2S793hJe1xSJ3t55EGJ2PmxXP376PjJrFS79Rkl3AbgMzECCECkgQYEVC8wsupmpcEACPEXOokIUfCy/nyDDERqTfRoCtmeIHG0KJJ7YsQIkBJDYwIKNX5S8CiNivYbXg2nu8+jQYlgoRVlV7P5VVU4dJ1HrX27sihBhfbPPcttybBr/+G7WRXLCdK15h17yI/DU4IIRpJDYwISgk505uvNrLPJ0KrJ+ayWhIuOeaugTklr6KK+Y4qLIrSsQMNNqrFfeutyZ9ddsVy5lS8TJ2imjUxQggRDOQWkghsRsFu2uMsHXAX9Rv+lzm2uQDkDbiL+Y4qr1sf4jSMmRePP7v6Df/LnIqXIe1xIk/9WQLyZymECAoSYETg8ggvpM4mF8gbN4f/3QBzbHOpH3iA2ePmyAW3pTxCTLYyn2i93uvPFpAQI4QIGhJgRODSXM1WG+UOTKTvgLsA3HUbTS60drvd3OU63PmsJ0qdTYN9Ps5Po0G5gISnvf9sAVyhWRYnhAgxEmBE4PLRpC6voooGFBb1uwNNtaBXVJkXXqOfiYSXUyxq4/JzI8TY5/Pdp9Ec3hFD/KVO9yxXavMQI4QQgU4CjAgaRsHubGsioyp38vu9+zFKUkdV7jTDi/QtcTNCi2cPner8Re7wcqrHjlETI8vShRDBRgKMCAqe4SV3YCIMTATs/H7vl8wHRjv280iT8FL69hvuPZZuvt1v4+5om1atQFFVcw8jzxDz78QL0fQBXDPjDu/bShJihBBBSJZRi6Dg0vVmq41SU1O56utdjHZ8ia4qzcLLppXui3k4UVSVTStXUPr2G+ZjCdOmsad3HLt69UBVm9fEkPa4u95ICCGCiMzAiKDgq0md3W5H13Wu+noXuq7zWkE+d0zLMcNL092tw4HxejetXGEel8yZxa6e3Rl86BjJB45QXVDQPMQIIUSQkQAjgpJnwW5qaiqvFeSz59ARns/5LyIPfRt24cVms6GqKqmpqV4h5t+r3uBkj14kJl/K5L8+b27VADRvGCiEEEFEAowIOk3DC8BEl8LX1fv5MaEPgFd4qS4oAJdGQs50v4y3I6iq6rX5YsrU28zwUpdwERefWpnVrLBXQowQIkhJgBFBx+jz4lnz8uk3e4k6vN/9PO4amJSpt3ltDhnKjD8LI8REHd5vhpeo6m/NPxvwCC0u2ahRCBG8JMCIoNO0z0vp229QVrGbkQMH0efdtey/aRKbVq7gh82b6fPuWq+drUPJAscBLIrSbAdpm83G1n6D0a5KZ3ZSb6IO7/eqiQGZeRFCBD8JMCKoNS3YrR4wCBbn80PvOMrYDTdNYpjnxdo+/1SH3+ZN8oKNRVGatf6POryfrf0GsyVpOGMqv/KapWoaYoQQIpiF1xpTEXJ0TfMq2E2YNg0lQiX5wFEGH/qOzldd1XiysbeSavHTaNtW7sBEZlsTme+oMneSfrGhsxlerqz4CrvdDrhDyzW33I6uyW0jIURokBkYEdSaNqmrLihAb9DMEPPDxr+RNzqd3MpXvTaGBHdzPJeu+1yiHbBs89wB7NRrMGZe5juqcG34Xy7XNS5UbmPxPbeaxc6A1+okIYQIBTIDI0KGZ8Hu0B2fE/+zK+hS/AW/enSMz/Ay31GFRVH8POpWUi3u12Kfbz6UOzCR/658hUcrXkZXVBaPuxpwh5a0tDRsNps5EyOEEKFCZmBESPAML0aBasL816GiP4c/vQCA1++5i1x8bEsQTIymcx7t///17pM8WvEyCwbexx8G3I3FY4NLowZGk1tHQogQ0+oZmI0bNzJ58mT69OmDoii88847Xs/fc889KIri9WvSpEle5xw9epTbb7+dmJgYLrzwQu6//36+//57r3M+/fRT/uM//oNOnTrRr18/5s+fjxCn5dKarzayzydheA3xl5/AomnUb/hf+to+8RlePG+3BJrq/CXuXjYGo/2/bS4Hb+7P4FeX8a8rZvDoPS80q4mBxpkYIYQIJa0OMCdOnGDEiBEsXbr0tOdMmjSJAwcOmL/eeOMNr+dvv/12Pv/8c4qLi3nvvffYuHEjv/71r83nnU4nEyZMYMCAAWzdupUFCxbw9NNP86c//am1wxVhIiFnerPwYtw2Sli5j4QZM5hT8TI5la8Sge4zvKiBum+SReXw4vxmIebg57Ec/ewCNFXl2pueA3wX9gohRChq9S2kG2+8kRtvvPGM50RHR5OY6Htq/ssvv2Tt2rVs2bKFq06tEMnPz+cnP/kJv//97+nTpw8rVqygrq6Ol19+maioKC655BLKysrIy8vzCjpC+OQRXsyalwF3UT/wAHMqXgbg7pWVvHLLlGZdfQNpB2tzZ2kf3XNLfj2V7w/15erLvqHXJTXu19yksNel6/4ZuBBCdIB2+Sfnhg0b6NmzJ0OGDOGhhx7iyJEj5nOlpaVceOGFZngBGD9+PKqq8tFHH5nnXH/99URFRZnnTJw4kZ07d/Ldd9+1x5BFKNFcPgt2I8fNgbTHGfz9UdYlDGTysjeahZdA2sHac2fphGnTiJ+Rw+HF+axJu5aymlq6DulNr1X7zNtJTQt7g2p1lRBCtFKbF/FOmjSJX/ziF1itVvbu3ctvfvMbbrzxRkpLS7FYLFRVVdGzZ0/vQUREEBcXR1WVe8q7qqoKq9XqdU6vXr3M57p3797s+9bW1lJbW2seO53Otn5pIlh4NKlrVrA7cDY/T4WXl73BloHDQNfJ/vxzSk91q/XsKeOPPZROtykjQN0ll/DV1Vfx/cnvGHzoO9L/+i/3b/JR2CuEEKGuzQPMrbfeav7/ZZddxuWXX87FF1/Mhg0bSE9Pb+tvZ5o3bx7PPPNMu319EZxcuu6zYHdU5U4ANEXlpR3biLJXcc3o/qTEfw3QfA+lduzg67klQNNNGUtHpbGvczy1RSvd+xpd0InLKr8j+cBRqgsKGut+jNCiudp8fEIIEYjafRl1UlIS8fHx7Nmzh/T0dBITEzl06JDXOQ0NDRw9etSsm0lMTOTgwYNe5xjHp6uteeyxx8jNzTWPnU4n/fr1a8uXIoJQ09sonjUvT6em8nzOf/Fjz74AJB86Bra5VBeVcfjvnzSuavKsqWkHXlsCeOxn9LYezUq9E7d06WZuytj5yEEm2/5lBiygeYgRQogw0O4B5ptvvuHIkSP07u2+kKSkpHDs2DG2bt3KqFGjAFi/fj2apjF27FjznMcff5z6+noiIyMBKC4uZsiQIT5vH4G7cDg6Orq9X44IYr4KdiMPfQvAjz378rfyPaRu64ve8Ak/ZAznLz/5Bbm+CoLbooOvR0ddz266ALl8xA8DqrlX78SYyq+Iq/jKDC+arrl32vZR2CuEEOGk1QHm+++/Z8+ePeaxw+GgrKyMuLg44uLieOaZZ5g6dSqJiYns3buX2bNnk5yczMSJEwEYNmwYkyZN4sEHH6SwsJD6+nqmT5/OrbfeSp8+fQD41a9+xTPPPMP999/PnDlz2LFjB4sWLeKFF15oo5ctwpGmac0Kdo2al9cK8tlXvZ9dcbEMPlrDqB7/5NJXLgW9Hof1V1ibFATPtiZit9vNr9lS1flLwKKScInFq2Yld2AiF694Ff21jdDv//h04H3mfkZoGmO7diV96f9njhswQwwuaVInhAg/rQ4wH3/8sdcHtnHb5u677+bFF1/k008/5ZVXXuHYsWP06dOHCRMm8Nxzz3nNjqxYsYLp06eTnp6OqqpMnTqVxYsXm8/Hxsbyj3/8g+zsbEaNGkV8fDy//e1vZQm1OC/G+7ZpeAGYiIXSz3awu3ccAD2+OE7C8BpqlUhmuK7nEbud1cdrefOCnsy2JjL0r69iO/mD+2va5lG67Rv+1e1y1Ak3cef7f6Ni715++NnPSOUjvt7xCXkD76PP8Eu5x6JiLymh8/djyUx7nNKVr6JvrmDQBZcy+C/LiL/Uyf/+x33kD7iLB+zvgqaBqhLxkykAzQp7U2TmRQgRplodYMaNG4d+hv4S69atO+vXiIuL4/XXXz/jOZdffjn/93//19rhCXFWTXewNupJUmbkkNCrO9+vW8Xhje7l+gnDa7hO3cxtDXegXRDLtZv/Sc8vo9l88gcu/ewzhvfsSenJb9i0uZJ+l9Xx8dJvObxmFfxsMi7bPKCUzy76T97oksBtX+zgi0suYcehQ1z60UeUHOlP2eGBjDz0bw5/ZiP+Uicv//QW8gfcRQMK5cNG8lDESeri+/jclFF2lhZChDPZC0mEHc8mdU33UEqwz4c+JVT/LJPDf/8EBl5nNr9b1O8OLv/xBHtOfk9yzx6kpqdT+uZr7O4dxzVjBpBc+g6X74jh5ck3M+VanRs+KeW12Mk8kjyN236oJuH9ldgSLiItLY2GkzWUVexm8KFj9DlwkvhLncRcVsu3A/+TBza8y9Fxk1jZcyAjPVZQyc7SQgjRSAKMCG+eeyh5bj+QOhsGFlB6xMmX/XqZIebDfmO46utdHLSvY0VEBK7ecQw6cJTuy/dxuD6G+EudzOtWQMQn9fx+4P38fsBdjHZ8Sdy3e8yVRJ/86QVcDQ0MPnSM5ANHUFSdmMtqidbrmVmxnK/TbiM19WoGnqq3gcbVSbIpoxBCuEmAEWHNq0ldkw6+f/nJL8yCXYf9R/odOciWK4YBcO3h/XDoWywREQz+7nv0+nqUyEgSRtSBq55aJZLfD7iLKEVhzLd7cLlcWCwWuhyrxtXQgKqoZnjRNYW1EbO56TpIss0lyZoEpDbbEsDYWVoIIUQ7bSUgRFBKe6zZ9gOzrYmMqtzJK45eXDTqv5htTWSLdRj/umIc9T0vwtXQwK7uXVEiI9Hr66neHkWDGkm0Xs8jla9Sp+tsvigZi8WCy+XihwsTUBUVTdfYf1lnhr54H/Ezchj8l2VUf9612bYAsiWAEEL4JjMwQvhgdPAdVbnTq3dM1NtvcO2eb+liHcSPPXrTveYEu3tDwk1XkFz6Dod3xDDfejNTrtN55JPFJH53iEdGPkJyr3gS3nd30x1cdZTuPXTKavrQ5XA/Uqa561kOL86HGTkkpD0uHXWFEOIsJMAI4YMx62FzfNmsd8wjt9xO8sFvsX/2GZ3HjiXOUsmmzZXsG/kTPrb25r41q/hemcz6LincUbOGuD1duC95Grf95BYydnzC5ssuIy0tjWtO7b8ETXq6SEddIYQ4KwkwQpyBZ88jz+XX1flLSE1Pdxf/2uYB8K9ul9P7jpuItyby/d69WNIeAz7ish2fcNsP1fQZfik/yczgglMN8Jouh5ZuukII0XKKfqamLkHM6XQSGxtLTU0NMTEx/h6OEEIIIVqgpddvKeIVQgghRNCRACOEEEKIoCMBRgghhBBBRwKMEEIIIYKOBBghhBBCBB0JMEIIIYQIOhJghBBCCBF0JMAIIYQQIuhIgBFCCCFE0JEAI4QQQoigE7J7IRk7JDidTj+PRAghhBAtZVy3z7bTUcgGmOPHjwPQr18/P49ECCGEEK11/PhxYmNjT/t8yG7mqGka+/fvp1u3bhw/fpx+/frx9ddfy8aOfuJ0OuVnEADk5xAY5Ofgf/IzCAy+fg66rnP8+HH69OmDqp6+0iVkZ2BUVaVv374AKIoCQExMjLxR/Ux+BoFBfg6BQX4O/ic/g8DQ9OdwppkXgxTxCiGEECLoSIARQgghRNAJiwATHR3NU089RXR0tL+HErbkZxAY5OcQGOTn4H/yMwgM5/NzCNkiXiGEEEKErrCYgRFCCCFEaJEAI4QQQoigIwFGCCGEEEFHAowQQgghgk5YBpiioiLGjh1L586d6d69O1OmTPH3kMJWbW0tI0eORFEUysrK/D2csFFRUcH999+P1Wqlc+fOXHzxxTz11FPU1dX5e2ghb+nSpQwcOJBOnToxduxYNm/e7O8hhZV58+YxevRounXrRs+ePZkyZQo7d+7097DC2vPPP4+iKMycObNVvy/sAszbb7/NnXfeyb333sv27dv517/+xa9+9St/DytszZ49mz59+vh7GGHnq6++QtM0/vjHP/L555/zwgsvUFhYyG9+8xt/Dy2k/fWvfyU3N5ennnqKbdu2MWLECCZOnMihQ4f8PbSwYbfbyc7O5t///jfFxcXU19czYcIETpw44e+hhaUtW7bwxz/+kcsvv7z1v1kPI/X19fpFF12kv/TSS/4eitB1/f3339eHDh2qf/755zqgf/LJJ/4eUlibP3++brVa/T2MkDZmzBg9OzvbPHa5XHqfPn30efPm+XFU4e3QoUM6oNvtdn8PJewcP35cHzRokF5cXKynpqbqDz/8cKt+f1jNwGzbto1vv/0WVVW54oor6N27NzfeeCM7duzw99DCzsGDB3nwwQf5y1/+QpcuXfw9HAHU1NQQFxfn72GErLq6OrZu3cr48ePNx1RVZfz48ZSWlvpxZOGtpqYGQN77fpCdnU1mZqbX34nWCKsAU15eDsDTTz/NE088wXvvvUf37t0ZN24cR48e9fPowoeu69xzzz1kZWVx1VVX+Xs4AtizZw/5+fn813/9l7+HErIOHz6My+WiV69eXo/36tWLqqoqP40qvGmaxsyZM7n22mu59NJL/T2csPLmm2+ybds25s2bd85fIyQCzP/7f/8PRVHO+Mu45w/w+OOPM3XqVEaNGsWyZctQFIVVq1b5+VUEv5b+HPLz8zl+/DiPPfaYv4ccclr6M/D07bffMmnSJG6++WYefPBBP41ciI6XnZ3Njh07ePPNN/09lLDy9ddf8/DDD7NixQo6dep0zl8nJLYSqK6u5siRI2c8JykpiX/961/ccMMN/N///R/XXXed+dzYsWMZP348c+fObe+hhrSW/hxuueUW1qxZg6Io5uMulwuLxcLtt9/OK6+80t5DDVkt/RlERUUBsH//fsaNG8fVV1/N8uXLUdWQ+DdNQKqrq6NLly689dZbXisf7777bo4dO8a7777rv8GFoenTp/Puu++yceNGrFarv4cTVt555x1+/vOfY7FYzMdcLheKoqCqKrW1tV7PnU5Eew6yoyQkJJCQkHDW80aNGkV0dDQ7d+40A0x9fT0VFRUMGDCgvYcZ8lr6c1i8eDH/8z//Yx7v37+fiRMn8te//pWxY8e25xBDXkt/BuCeeUlLSzNnIiW8tK+oqChGjRpFSUmJGWA0TaOkpITp06f7d3BhRNd1cnJyWL16NRs2bJDw4gfp6el89tlnXo/de++9DB06lDlz5rQovECIBJiWiomJISsri6eeeop+/foxYMAAFixYAMDNN9/s59GFj/79+3sdd+3aFYCLL76Yvn37+mNIYefbb79l3LhxDBgwgN///vdUV1ebzyUmJvpxZKEtNzeXu+++m6uuuooxY8awcOFCTpw4wb333uvvoYWN7OxsXn/9dd599126detm1h/FxsbSuXNnP48uPHTr1q1ZzdEFF1xAjx49WlWLFFYBBmDBggVERERw5513cvLkScaOHcv69evp3r27v4cmRIcpLi5mz5497Nmzp1loDIG7ygHrP//zP6murua3v/0tVVVVjBw5krVr1zYr7BXt58UXXwRg3LhxXo8vW7aMe+65p+MHJM5ZSNTACCGEECK8yE1vIYQQQgQdCTBCCCGECDoSYIQQQggRdCTACCGEECLoSIARQgghRNCRACOEEEKIoCMBRgghhBBBRwKMEEIIIYKOBBghhBBCBB0JMEIIIYQIOhJghBBCCBF0JMAIIYQQIuj8/7WOozOaAn7MAAAAAElFTkSuQmCC",
+ "text/plain": [
+ "<Figure size 640x480 with 1 Axes>"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "print(\"Local minimizer in area of pos. curvature: kappa =\", kappamin_pos, \", alpha =\", alphamin_pos, \", Q =\", Emin_pos)\n",
+ "print(\"Local minimizer in area of neg. curvature: kappa =\", kappamin_neg, \", alpha =\", alphamin_neg, \", Q =\", Emin_neg)\n",
+ "\n",
+ "n=100\n",
+ "bending_path_pos=np.outer(np.array(np.linspace(0,2,n)),np.array([kappamin_pos,0]))\n",
+ "bending_path_neg=np.outer(np.array(np.linspace(0,2,n)),np.array([-kappamin_neg,np.pi/2]))\n",
+ "\n",
+ "for i in range(0,n):\n",
+ " plt.plot(bending_path_pos[i,0],energy(bending_path_pos[i,0],0,Q,B), 'x') \n",
+ " plt.plot(bending_path_neg[i,0],energy(bending_path_neg[i,0],np.pi/2,Q,B), 'x') \n",
+ "\n",
+ "print(energy(1/1,0,Q,B))\n",
+ "print(energy(-1/n,np.pi/2,Q,B))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "base",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.10.12"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/experiment/wood-bilayer_PLOS/explorer.ipynb b/experiment/wood-bilayer_PLOS/explorer.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..d3433b5b3b4231339d88b5a4dc850f2a2258e805
--- /dev/null
+++ b/experiment/wood-bilayer_PLOS/explorer.ipynb
@@ -0,0 +1,223 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "%%capture\n",
+ "#!/usr/bin/env python3\n",
+ "# -*- coding: utf-8 -*-\n",
+ "\"\"\"\n",
+ "Created on Wed Jul 6 13:17:28 2022\n",
+ "\n",
+ "@author: stefan\n",
+ "\"\"\"\n",
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt\n",
+ "import matplotlib.colors as colors\n",
+ "from matplotlib.ticker import LogLocator\n",
+ "import codecs\n",
+ "import re\n",
+ "import json\n",
+ "\n",
+ "def energy(kappa,alpha,Q,B) :\n",
+ " G=kappa*np.array([[np.cos(alpha)**2],[np.sin(alpha)**2],[np.sqrt(2)*np.cos(alpha)*np.sin(alpha)]])-B\n",
+ " return np.matmul(np.transpose(G),np.matmul(Q,G))[0,0]\n",
+ "\n",
+ "def xytokappaalpha(x,y):\n",
+ " \n",
+ " if y>0:\n",
+ " return [np.sqrt(x**2+y**2), np.abs(np.arctan2(y,x))]\n",
+ " else:\n",
+ " return [-np.sqrt(x**2+y**2), np.abs(np.arctan2(y,x))]\n",
+ "\n",
+ "# Read effective quantites\n",
+ "def ReadEffectiveQuantities(QFilePath, BFilePath):\n",
+ " # Read Output Matrices (effective quantities)\n",
+ " # From Cell-Problem output Files : ../outputs/Qmatrix.txt , ../outputs/Bmatrix.txt\n",
+ " # -- Read Matrix Qhom\n",
+ " X = []\n",
+ " # with codecs.open(path + '/outputs/QMatrix.txt', encoding='utf-8-sig') as f:\n",
+ " with codecs.open(QFilePath, encoding='utf-8-sig') as f:\n",
+ " for line in f:\n",
+ " s = line.split()\n",
+ " X.append([float(s[i]) for i in range(len(s))])\n",
+ " Q = np.array([[X[0][2], X[1][2], X[2][2]],\n",
+ " [X[3][2], X[4][2], X[5][2]],\n",
+ " [X[6][2], X[7][2], X[8][2]] ])\n",
+ "\n",
+ " # -- Read Beff (as Vector)\n",
+ " X = []\n",
+ " # with codecs.open(path + '/outputs/BMatrix.txt', encoding='utf-8-sig') as f:\n",
+ " with codecs.open(BFilePath, encoding='utf-8-sig') as f:\n",
+ " for line in f:\n",
+ " s = line.split()\n",
+ " X.append([float(s[i]) for i in range(len(s))])\n",
+ " B = np.array([X[0][2], X[1][2], X[2][2]])\n",
+ " return Q, B"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Load data from simulation:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "#--- Select specific experiment [x, y] with date from results_x/y\n",
+ "data=[5,6]\n",
+ "#DataPath = './experiment/wood-bilayer_PLOS/results_' + str(data[0]) + '/' +str(data[1])\n",
+ "DataPath = './results_' + str(data[0]) + '/' +str(data[1])\n",
+ "QFilePath = DataPath + '/QMatrix.txt'\n",
+ "BFilePath = DataPath + '/BMatrix.txt'\n",
+ "ParameterPath = DataPath + '/parameter.txt'\n",
+ "#\n",
+ "# Read Thickness from parameter file (needed for energy scaling)\n",
+ "with open(ParameterPath , 'r') as file:\n",
+ " parameterFile = file.read()\n",
+ "thickness = float(re.findall(r'(?m)h = (\\d?\\d?\\d?\\.?\\d+[Ee]?[+\\-]?\\d?\\d?)',parameterFile)[0])\n",
+ "energyscalingfactor = thickness**2\n",
+ "# Read Q and B\n",
+ "Q, B = ReadEffectiveQuantities(QFilePath,BFilePath)\n",
+ "Q=0.5*(np.transpose(Q)+Q) # symmetrize\n",
+ "B=np.transpose([B])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Sampling of local energy"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "N=500\n",
+ "length=4\n",
+ "r, theta = np.meshgrid(np.linspace(0,length,N),np.radians(np.linspace(0, 360, N)))\n",
+ "E=np.zeros(np.shape(r))\n",
+ "for i in range(0,N): \n",
+ " for j in range(0,N): \n",
+ " if theta[i,j]<np.pi:\n",
+ " E[i,j]=energy(r[i,j],theta[i,j],Q,B) * energyscalingfactor\n",
+ " else:\n",
+ " E[i,j]=energy(-r[i,j],theta[i,j],Q,B) * energyscalingfactor\n",
+ "# \n",
+ "# Compute Minimizer\n",
+ "[imin,jmin]=np.unravel_index(E.argmin(),(N,N))\n",
+ "kappamin=r[imin,jmin]\n",
+ "alphamin=theta[imin,jmin]\n",
+ "# Positiv curvature region\n",
+ "N_mid=int(N/2)\n",
+ "[imin,jmin]=np.unravel_index(E[:N_mid,:].argmin(),(N_mid,N))\n",
+ "kappamin_pos=r[imin,jmin]\n",
+ "alphamin_pos=theta[imin,jmin]\n",
+ "Emin_pos=E[imin,jmin]\n",
+ "# Negative curvature region\n",
+ "[imin,jmin]=np.unravel_index(E[N_mid:,:].argmin(),(N_mid,N))\n",
+ "kappamin_neg=r[imin+N_mid,jmin]\n",
+ "alphamin_neg=theta[imin+N_mid,jmin]\n",
+ "Emin_neg=E[imin+N_mid,jmin]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 28,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Local minimizer in area of pos. curvature: kappa = 1.819639278557114 , alpha = 0.0 , Q = 0.1299383783013474\n",
+ "Local minimizer in area of neg. curvature: kappa = 2.797595190380761 , alpha = 4.709241091954239 , Q = 0.09277799602960847\n"
+ ]
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGdCAYAAAAMm0nCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAy2ElEQVR4nO3de3RU5b3/8c8kYYaLmWAKSeAQMJY2GrlYLOJg5VJDgo1dckqP0ipiBS04wQO0EDnHu+vXULRSVC61tcQeRcTaeCEFmiYkVIlAKRSChaJGoQ2TYCkzATGEzPP7g2aXkQlmkpBkJ+/XWrNw7+c7O8+XnTAf9zyz4zDGGAEAANhIVHtPAAAAIFIEGAAAYDsEGAAAYDsEGAAAYDsEGAAAYDsEGAAAYDsEGAAAYDsEGAAAYDsx7T2BCyUYDKqyslKxsbFyOBztPR0AANAExhjV1NSof//+iopq/DpLpw0wlZWVSk5Obu9pAACAZjh06JAGDBjQ6HinDTCxsbGSzvwFuN3udp4NAABoikAgoOTkZOt1vDGdNsA0vG3kdrsJMAAA2MznLf9gES8AALAdAgwAALAdAgwAALAdAgwAALAdAgwAALAdAgwAALAdAgwAALAdAgwAALAdAgwAALAdAkwTLCn8q54qOhB27KmiA1pS+Nc2nhEAAF0bAaYJoqMcejJMiHmq6ICeLPyroqP4bdcAALSlTvu7kFrTvdd/SZL05L+utNx7/Zes8DJvwpetcQAA0DYIME10doh5pvg9naoPEl4AAGgnvIUUgXuv/5Kc0VE6VR+UMzqK8AIAQDshwETgqaIDVng5VR9sdGEvAAC4sHgLqYk+u+alYVsSV2IAAGhjBJgmCLdgN9zCXgAA0DYIME1QHzRhF+w2bNcHTXtMCwCALsthjOmUr76BQEBxcXHy+/1yu93tPR0AANAETX39ZhEvAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwnYgCzIoVKzRs2DC53W653W55PB6tX79eknT06FHNnj1bqamp6tGjhwYOHKh7771Xfr8/5BgHDx5UVlaWevbsqYSEBM2fP1+nT58OqSkpKdGIESPkcrk0ePBg5eXltaxLAADQqcREUjxgwAAtWrRIX/rSl2SM0fPPP6+bbrpJO3fulDFGlZWVeuKJJ5SWlqaPPvpIM2fOVGVlpX79619Lkurr65WVlaWkpCRt2bJFhw8f1u23365u3brpRz/6kSSpoqJCWVlZmjlzpl588UUVFRVpxowZ6tevnzIzM1v/bwAAANiOwxhjWnKA+Ph4Pf7445o+ffo5Y6+88opuu+02nThxQjExMVq/fr1uvPFGVVZWKjExUZK0cuVK5eTk6MiRI3I6ncrJyVFBQYHKy8ut40yZMkXHjh3Thg0bmjyvQCCguLg4+f1+ud3ulrQIAADaSFNfv5u9Bqa+vl5r1qzRiRMn5PF4wtY0fPGYmDMXesrKyjR06FArvEhSZmamAoGA9u7da9Wkp6eHHCczM1NlZWXnnU9tba0CgUDIAwAAdE4RB5g9e/booosuksvl0syZM5Wfn6+0tLRz6j7++GM99thjuvvuu619Pp8vJLxIsrZ9Pt95awKBgE6ePNnovHJzcxUXF2c9kpOTI20NAADYRMQBJjU1Vbt27dLWrVs1a9YsTZs2Te+++25ITSAQUFZWltLS0vTwww+31lzPa+HChfL7/dbj0KFDbfJ1AQBA24toEa8kOZ1ODR48WJJ01VVXafv27Vq6dKl+9rOfSZJqamo0ceJExcbGKj8/X926dbOem5SUpG3btoUcr6qqyhpr+LNh39k1brdbPXr0aHReLpdLLpcr0nYAAIANtfg+MMFgULW1tZLOXHnJyMiQ0+nUG2+8oe7du4fUejwe7dmzR9XV1da+wsJCud1u620oj8ejoqKikOcVFhY2us4GAAB0PRFdgVm4cKFuuOEGDRw4UDU1NVq9erVKSkq0ceNGK7x88skneuGFF0IW0vbt21fR0dHKyMhQWlqapk6dqsWLF8vn8+n++++X1+u1rp7MnDlTzzzzjBYsWKA777xTxcXFWrt2rQoKClq/ewAAYEsRBZjq6mrdfvvtOnz4sOLi4jRs2DBt3LhREyZMUElJibZu3SpJ1ltMDSoqKnTJJZcoOjpa69at06xZs+TxeNSrVy9NmzZNjz76qFWbkpKigoICzZ07V0uXLtWAAQP0i1/8gnvAAAAAS4vvA9NRcR8YAADs54LfBwYAAKC9EGAAAIDtEGAAAIDtEGAAAIDtEGAAAIDtEGAAAIDtEGAAAIDtEGAAAIDtEGAAAIDtEGAAAIDtEGAAAIDtEGAAAIDtEGAAAIDtEGAAAIDtEGAAAIDtEGAAAIDtEGAAAIDtEGAAAIDtEGAAAIDtEGAAAIDtEGAAAIDtEGAAAIDtEGAAAIDtEGAAAIDtEGAAAIDtEGAAAIDtEGC6kk25Uuni8GOli8+MAwBgAwSYriQqWtr0/84NMaWLz+yPim6feQEAEKGY9p4A2tDYBWf+3PT//r3dEF7G/++/xwEA6OAIMF3N2SFm8+NS/SnCCwDAdngLqSsau0CKdp4JL9FOwgsAwHYIMF1R6eJ/h5f6U40v7AUAoIMiwHQ1Z695eeDImT/DLewFAKADYw1MVxJuwW64hb0AAHRwBJiuJFgffsFuw3awvu3nBABAMziMMaa9J3EhBAIBxcXFye/3y+12t/d0AABAEzT19Zs1MAAAwHYIMAAAwHYIMAAAwHYIMAAAwHYIMAAAwHYiCjArVqzQsGHD5Ha75Xa75fF4tH79emv8008/ldfr1Re+8AVddNFFmjx5sqqqqkKOcfDgQWVlZalnz55KSEjQ/Pnzdfr06ZCakpISjRgxQi6XS4MHD1ZeXl7zOwQAAJ1ORAFmwIABWrRokXbs2KE//vGP+vrXv66bbrpJe/fulSTNnTtXb775pl555RWVlpaqsrJS3/rWt6zn19fXKysrS6dOndKWLVv0/PPPKy8vTw8++KBVU1FRoaysLI0fP167du3SnDlzNGPGDG3cuLGVWgYAAHbX4vvAxMfH6/HHH9e3v/1t9e3bV6tXr9a3v/1tSdK+fft0+eWXq6ysTNdcc43Wr1+vG2+8UZWVlUpMTJQkrVy5Ujk5OTpy5IicTqdycnJUUFCg8vJy62tMmTJFx44d04YNG5o8L+4DAwCA/Vzw+8DU19drzZo1OnHihDwej3bs2KG6ujqlp6dbNZdddpkGDhyosrIySVJZWZmGDh1qhRdJyszMVCAQsK7ilJWVhRyjoabhGI2pra1VIBAIeQAAgM4p4gCzZ88eXXTRRXK5XJo5c6by8/OVlpYmn88np9Op3r17h9QnJibK5/NJknw+X0h4aRhvGDtfTSAQ0MmTJxudV25uruLi4qxHcnJypK0BAACbiDjApKamateuXdq6datmzZqladOm6d13370Qc4vIwoUL5ff7rcehQ4fae0oAAOACifiXOTqdTg0ePFiSdNVVV2n79u1aunSpbrnlFp06dUrHjh0LuQpTVVWlpKQkSVJSUpK2bdsWcryGTymdXfPZTy5VVVXJ7XarR48ejc7L5XLJ5XJF2g4AALChFt8HJhgMqra2VldddZW6deumoqIia2z//v06ePCgPB6PJMnj8WjPnj2qrq62agoLC+V2u5WWlmbVnH2MhpqGYwAAAER0BWbhwoW64YYbNHDgQNXU1Gj16tUqKSnRxo0bFRcXp+nTp2vevHmKj4+X2+3W7Nmz5fF4dM0110iSMjIylJaWpqlTp2rx4sXy+Xy6//775fV6rasnM2fO1DPPPKMFCxbozjvvVHFxsdauXauCgoLW7x4AANhSRAGmurpat99+uw4fPqy4uDgNGzZMGzdu1IQJEyRJS5YsUVRUlCZPnqza2lplZmZq+fLl1vOjo6O1bt06zZo1Sx6PR7169dK0adP06KOPWjUpKSkqKCjQ3LlztXTpUg0YMEC/+MUvlJmZ2UotAwAAu2vxfWA6Ku4DAwCA/Vzw+8AAAAC0FwIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMAACwHQIMbGn5ruVa+eeVYcdW/nmllu9a3sYzAgC0JQIMbCnKEaVlu5adE2JW/nmllu1apigH39oA0JnFtPcEgOaYOXymJGnZrmXWdkN48V7ptcYBAJ0TAQa2dXaIeXb3s6oL1hFeAKCLiOg6e25urkaOHKnY2FglJCRo0qRJ2r9/f0iNz+fT1KlTlZSUpF69emnEiBF69dVXQ2qOHj2qW2+9VW63W71799b06dN1/PjxkJrdu3fruuuuU/fu3ZWcnKzFixc3s0V0ZjOHz1S3qG6qC9apW1Q3wgsAdBERBZjS0lJ5vV698847KiwsVF1dnTIyMnTixAmr5vbbb9f+/fv1xhtvaM+ePfrWt76lm2++WTt37rRqbr31Vu3du1eFhYVat26dNm/erLvvvtsaDwQCysjI0KBBg7Rjxw49/vjjevjhh/Xss8+2QsvoTFb+eaUVXuqCdY0u7AUAdDKmBaqrq40kU1paau3r1auX+dWvfhVSFx8fb37+858bY4x59913jSSzfft2a3z9+vXG4XCYv//978YYY5YvX24uvvhiU1tba9Xk5OSY1NTUJs/N7/cbScbv9zerN3R8K3atMEPyhpgVu1aE3QYA2E9TX79b9FENv98vSYqPj7f2jR49Wi+//LKOHj2qYDCoNWvW6NNPP9W4ceMkSWVlZerdu7e++tWvWs9JT09XVFSUtm7datWMGTNGTqfTqsnMzNT+/fv1z3/+syVTRicRbsHuzOEz5b3SG/bTSQCAzqXZi3iDwaDmzJmja6+9VkOGDLH2r127Vrfccou+8IUvKCYmRj179lR+fr4GDx4s6cwamYSEhNBJxMQoPj5ePp/PqklJSQmpSUxMtMYuvvjic+ZTW1ur2tpaazsQCDS3NdhA0ATDLtht2A6aYHtMCwDQRpodYLxer8rLy/XWW2+F7H/ggQd07Ngx/f73v1efPn302muv6eabb9Yf/vAHDR06tMUTbkxubq4eeeSRC3Z8dCz3XHlPo2Ms5AWAzq9ZASY7O9tafDtgwABr//vvv69nnnlG5eXluuKKKyRJw4cP1x/+8ActW7ZMK1euVFJSkqqrq0OOd/r0aR09elRJSUmSpKSkJFVVVYXUNGw31HzWwoULNW/ePGs7EAgoOTm5Oe0BAIAOLqI1MMYYZWdnKz8/X8XFxee8zfPJJ5+cOWhU6GGjo6MVDJ65pO/xeHTs2DHt2LHDGi8uLlYwGNSoUaOsms2bN6uurs6qKSwsVGpqati3jyTJ5XLJ7XaHPAAAQOcUUYDxer164YUXtHr1asXGxsrn88nn8+nkyZOSpMsuu0yDBw/W97//fW3btk3vv/++fvKTn6iwsFCTJk2SJF1++eWaOHGi7rrrLm3btk1vv/22srOzNWXKFPXv31+S9N3vfldOp1PTp0/X3r179fLLL2vp0qUhV1gAAEDX5TDGmCYXOxxh969atUp33HGHJOnAgQO677779NZbb+n48eMaPHiwfvjDH2rq1KlW/dGjR5Wdna0333xTUVFRmjx5sp566ilddNFFVs3u3bvl9Xq1fft29enTR7Nnz1ZOTk6TGwsEAoqLi5Pf7+dqDAAANtHU1++IAoydEGAAALCfpr5+8yt7AQCA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBggHZ05OlndGT58vBjy5fryNPPtPGMAMAeCDBAe4qO0sdPPX1OiDmyfLk+fuppKZofUQAIJ6J/HXNzczVy5EjFxsYqISFBkyZN0v79+8+pKysr09e//nX16tVLbrdbY8aM0cmTJ63xo0eP6tZbb5Xb7Vbv3r01ffp0HT9+POQYu3fv1nXXXafu3bsrOTlZixcvbmaLQMfV95571Ofe2SEhpiG89Ll3tvrec087zxAAOqaIAkxpaam8Xq/eeecdFRYWqq6uThkZGTpx4oRVU1ZWpokTJyojI0Pbtm3T9u3blZ2draiof3+pW2+9VXv37lVhYaHWrVunzZs36+6777bGA4GAMjIyNGjQIO3YsUOPP/64Hn74YT377LOt0DLQsZwdYvYNHUZ4AYAmcBhjTHOffOTIESUkJKi0tFRjxoyRJF1zzTWaMGGCHnvssbDP+ctf/qK0tDRt375dX/3qVyVJGzZs0De+8Q397W9/U//+/bVixQr97//+r3w+n5xOpyTpvvvu02uvvaZ9+/Y1aW6BQEBxcXHy+/1yu93NbRFoM/uGDpOpq5OjWzddtmd3e08HANpFU1+/W/QGu9/vlyTFx8dLkqqrq7V161YlJCRo9OjRSkxM1NixY/XWW29ZzykrK1Pv3r2t8CJJ6enpioqK0tatW62aMWPGWOFFkjIzM7V//37985//DDuX2tpaBQKBkAdgF0eWL7fCi6mra3RhLwDgjGYHmGAwqDlz5ujaa6/VkCFDJEkffPCBJOnhhx/WXXfdpQ0bNmjEiBG6/vrrdeDAAUmSz+dTQkJCyLFiYmIUHx8vn89n1SQmJobUNGw31HxWbm6u4uLirEdycnJzWwPa1NlrXi7bs/ucNTEAgHM1O8B4vV6Vl5drzZo11r5gMChJ+v73v6/vfe97+spXvqIlS5YoNTVVv/zlL1s+2/NYuHCh/H6/9Th06NAF/XpAawi3YDfcwl4AQKiY5jwpOzvbWnw7YMAAa3+/fv0kSWlpaSH1l19+uQ4ePChJSkpKUnV1dcj46dOndfToUSUlJVk1VVVVITUN2w01n+VyueRyuZrTDtB+6oNhF+xa2/XBdpgUAHR8EV2BMcYoOztb+fn5Ki4uVkpKSsj4JZdcov79+5/z0eq//vWvGjRokCTJ4/Ho2LFj2rFjhzVeXFysYDCoUaNGWTWbN29WXV2dVVNYWKjU1FRdfPHFkXUIdGB9Z2c3+mmjvvfco76zs9t4RgBgDxEFGK/XqxdeeEGrV69WbGysfD6ffD6fdY8Xh8Oh+fPn66mnntKvf/1rvffee3rggQe0b98+TZ8+XdKZqzETJ07UXXfdpW3btuntt99Wdna2pkyZov79+0uSvvvd78rpdGr69Onau3evXn75ZS1dulTz5s1r5fYBAIAdRfQxaofDEXb/qlWrdMcdd1jbixYt0rJly3T06FENHz5cixcv1te+9jVr/OjRo8rOztabb76pqKgoTZ48WU899ZQuuugiq2b37t3yer3avn27+vTpo9mzZysnJ6fJjfExagAA7Kepr98tug9MR0aAAQDAftrkPjAAAADtgQADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADoMW2vfmBthdUhB3bXlChbW9+0MYzAtDZEWAAtJgjyqFtb1acE2LOhJcKOaIc7TQzAJ1VTHtPAID9jcxKkSRte7PC2m4IL1d/M8UaB4DWQoAB0CrODjF/XP+hgqcN4QXABcNbSABazcisFEXFOBQ8bRQV4yC8ALhgCDAAWs32ggorvARPm0YX9gJAS/EWEoBW8dk1Lw3bkrgSA6DVEWAAtFi4BbvhFvYCQGuJ6C2k3NxcjRw5UrGxsUpISNCkSZO0f//+sLXGGN1www1yOBx67bXXQsYOHjyorKws9ezZUwkJCZo/f75Onz4dUlNSUqIRI0bI5XJp8ODBysvLi6gxAG3HBMMv2B2ZlaKrv5kiEzTtNDMAnVVEV2BKS0vl9Xo1cuRInT59Wv/zP/+jjIwMvfvuu+rVq1dI7U9/+lM5HOfe+6G+vl5ZWVlKSkrSli1bdPjwYd1+++3q1q2bfvSjH0mSKioqlJWVpZkzZ+rFF19UUVGRZsyYoX79+ikzM7MF7QK4EK7+5qWNjnHlBcCF4DDGNPt/jY4cOaKEhASVlpZqzJgx1v5du3bpxhtv1B//+Ef169dP+fn5mjRpkiRp/fr1uvHGG1VZWanExERJ0sqVK5WTk6MjR47I6XQqJydHBQUFKi8vt445ZcoUHTt2TBs2bGjS3AKBgOLi4uT3++V2u5vbIgAAaENNff1u0aeQ/H6/JCk+Pt7a98knn+i73/2uli1bpqSkpHOeU1ZWpqFDh1rhRZIyMzMVCAS0d+9eqyY9PT3keZmZmSorK2t0LrW1tQoEAiEPAADQOTU7wASDQc2ZM0fXXnuthgwZYu2fO3euRo8erZtuuins83w+X0h4kWRt+3y+89YEAgGdPHky7HFzc3MVFxdnPZKTk5vbGgAA6OCa/Skkr9er8vJyvfXWW9a+N954Q8XFxdq5c2erTC4SCxcu1Lx586ztQCBAiAEAoJNq1hWY7OxsrVu3Tps2bdKAAQOs/cXFxXr//ffVu3dvxcTEKCbmTD6aPHmyxo0bJ0lKSkpSVVVVyPEathvecmqsxu12q0ePHmHn5HK55Ha7Qx4AAKBziijAGGOUnZ2t/Px8FRcXKyUl9NMF9913n3bv3q1du3ZZD0lasmSJVq1aJUnyeDzas2ePqqurrecVFhbK7XYrLS3NqikqKgo5dmFhoTweT8QNAgCAzieit5C8Xq9Wr16t119/XbGxsdaalbi4OPXo0UNJSUlhF+4OHDjQCjsZGRlKS0vT1KlTtXjxYvl8Pt1///3yer1yuVySpJkzZ+qZZ57RggULdOedd6q4uFhr165VQUFBS/sFAACdQERXYFasWCG/369x48apX79+1uPll19u8jGio6O1bt06RUdHy+Px6LbbbtPtt9+uRx991KpJSUlRQUGBCgsLNXz4cP3kJz/RL37xC+4BAwAAJLXwPjAdGfeBAQDAftrkPjAAAADtgQADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAP+y5ZUXVfbqS2HHyl59SVteebGNZwSgMQQYAPgXR1SUtqw9N8SUvfqStqx9UY4o/skEOoqY9p4AAHQUnsnfkSRtWfuitd0QXkbffKs1DqD9EWAA4Cxnh5itv3lZ9adPE16ADojroQDwGZ7J31F0TIzqT59WdEwM4QXogAgwAPAZZa++ZIWX+tOnG13YC6D98BYSAJzls2teGrYlcSUG6EAIMADwL+EW7IZb2Aug/RFgAOBfTDAYdsFuw7YJBttjWgDCcBhjTHtP4kIIBAKKi4uT3++X2+1u7+kAAIAmaOrrN4t4AQCA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7RBgAACA7UQUYHJzczVy5EjFxsYqISFBkyZN0v79+63xo0ePavbs2UpNTVWPHj00cOBA3XvvvfL7/SHHOXjwoLKystSzZ08lJCRo/vz5On36dEhNSUmJRowYIZfLpcGDBysvL6/5XQIAgE4logBTWloqr9erd955R4WFhaqrq1NGRoZOnDghSaqsrFRlZaWeeOIJlZeXKy8vTxs2bND06dOtY9TX1ysrK0unTp3Sli1b9PzzzysvL08PPvigVVNRUaGsrCyNHz9eu3bt0pw5czRjxgxt3LixldoGAAB25jDGmOY++ciRI0pISFBpaanGjBkTtuaVV17RbbfdphMnTigmJkbr16/XjTfeqMrKSiUmJkqSVq5cqZycHB05ckROp1M5OTkqKChQeXm5dZwpU6bo2LFj2rBhQ5Pm1tRfxw0AADqOpr5+t2gNTMNbQ/Hx8eetcbvdiomJkSSVlZVp6NChVniRpMzMTAUCAe3du9eqSU9PDzlOZmamysrKWjJdAADQScQ094nBYFBz5szRtddeqyFDhoSt+fjjj/XYY4/p7rvvtvb5fL6Q8CLJ2vb5fOetCQQCOnnypHr06HHO16qtrVVtba21HQgEmtcYAADo8Jp9Bcbr9aq8vFxr1qwJOx4IBJSVlaW0tDQ9/PDDzf0yTZabm6u4uDjrkZycfMG/JgAAaB/NCjDZ2dlat26dNm3apAEDBpwzXlNTo4kTJyo2Nlb5+fnq1q2bNZaUlKSqqqqQ+obtpKSk89a43e6wV18kaeHChfL7/dbj0KFDzWkNAADYQEQBxhij7Oxs5efnq7i4WCkpKefUBAIBZWRkyOl06o033lD37t1Dxj0ej/bs2aPq6mprX2Fhodxut9LS0qyaoqKikOcVFhbK4/E0OjeXyyW32x3yAICuyF/4kQJFB8OOBYoOyl/4URvPCGh9EQUYr9erF154QatXr1ZsbKx8Pp98Pp9Onjwp6d/h5cSJE3ruuecUCASsmvr6eklSRkaG0tLSNHXqVP35z3/Wxo0bdf/998vr9crlckmSZs6cqQ8++EALFizQvn37tHz5cq1du1Zz585t5fYBoPNxRDkUCBNiAkUHFSj8SI4oRzvNDGg9EX2M2uEI/02/atUq3XHHHSopKdH48ePD1lRUVOiSSy6RJH300UeaNWuWSkpK1KtXL02bNk2LFi2yPqkknbmR3dy5c/Xuu+9qwIABeuCBB3THHXc0uTE+Rg2gK2sIK+4Jg+S+fuA520BH1dTX7xbdB6YjI8AA6OoaQouiHVK9IbzAFtrkPjAAgI7Lff1AK7wo2kF4QadCgAGATipQdNAKL6o3jS7sBeyo2TeyAwB0XI2tgZHElRh0CgQYAOhkwi3YbfiTEIPOggADAJ2MCYZfsNuwbYKd8rMb6GIIMADQycRNGNToGFde0FmwiBcAANgOAQYAANgOAQYAANgOAQYAANgOAQYAANgOAQYAANgOAQYAANgOAQYAANgOAQYAANgOAQYAANgOAQYAANgOAQYA0CFt2rRJpaWlYcdKS0u1adOmNp4ROhICDACgQ4qKigobYhrCS1QUL2FdGb+NGgDQIY0dO1aSrCstY8eOtcLL+PHjrXF0TQQYAECHdXaI2bx5s+rr6wkvkMRbSACADm7s2LGKjo5WfX29oqOjCS+QRIABAHRwpaWlVnipr69vdGEvuhbeQgIAdFifXfNy9qePuBLTtRFgAAAdUrgFu+EW9qJrIsAAADqkYDAYdsFuw3YwGGyPaaGDcBhjTHtP4kIIBAKKi4uT3++X2+1u7+kAAIAmaOrrN4t4AQCA7RBgAACA7RBgAAC4gD74YKkqKp4OO1ZR8bQ++GBpG8+ocyDAAABwATkcUfqg4qfnhJiKiqf1QcVP5XDwUtwcfAoJAIALKCVltiTpg4qfWtsN4eXSlDnWOCJDgAEA4AI7O8RUfLhcxpwivLQQ160AAGgDKSmz5XA4ZcwpORxO24WXxysO68kPfWHHnvzQp8crDrfpfAgwAAC0gYqKp63wYsypRhf2dlTRDocWV/jOCTFPfujT4gqfoh2ONp0PbyEBAHCBfXbNS8O2JNtciZl3SZIkaXGFz9puCC8LUpKs8bZCgAEA4AIKt2A33MJeOzg7xPz0wyqdMqZdwotEgAEA4IIyJhh2wW7DtjH2+p1O8y5JssKL0+Fol/AiEWAAALigLr30vxsds8uVl7M9+aHPCi+njNGTH/raJcSwiBcAADTJ2WteDo4brgUpSWEX9rYFrsAAAIDPFW7BbriFvW0loiswubm5GjlypGJjY5WQkKBJkyZp//79ITWffvqpvF6vvvCFL+iiiy7S5MmTVVVVFVJz8OBBZWVlqWfPnkpISND8+fN1+vTpkJqSkhKNGDFCLpdLgwcPVl5eXvM6BAAALVbfyILdeZckaUFKkuqNadP5RBRgSktL5fV69c4776iwsFB1dXXKyMjQiRMnrJq5c+fqzTff1CuvvKLS0lJVVlbqW9/6ljVeX1+vrKwsnTp1Slu2bNHzzz+vvLw8Pfjgg1ZNRUWFsrKyNH78eO3atUtz5szRjBkztHHjxlZoGQAARGp+Sr9Gr7DMuyRJ81P6tel8HMY0PzIdOXJECQkJKi0t1ZgxY+T3+9W3b1+tXr1a3/72tyVJ+/bt0+WXX66ysjJdc801Wr9+vW688UZVVlYqMTFRkrRy5Url5OToyJEjcjqdysnJUUFBgcrLy62vNWXKFB07dkwbNmxo0twCgYDi4uLk9/vldrub2yIAAGhDTX39btEiXr/fL0mKj4+XJO3YsUN1dXVKT0+3ai677DINHDhQZWVlkqSysjINHTrUCi+SlJmZqUAgoL1791o1Zx+joabhGOHU1tYqEAiEPAAAQOfU7AATDAY1Z84cXXvttRoyZIgkyefzyel0qnfv3iG1iYmJ8vl8Vs3Z4aVhvGHsfDWBQEAnT54MO5/c3FzFxcVZj+Tk5Oa2BgAAOrhmBxiv16vy8nKtWbOmNefTbAsXLpTf77cehw4dau8pAQCAC6RZH6POzs7WunXrtHnzZg0YMMDan5SUpFOnTunYsWMhV2GqqqqUlJRk1Wzbti3keA2fUjq75rOfXKqqqpLb7VaPHj3CzsnlcsnlcjWnHQAAYDMRXYExxig7O1v5+fkqLi5WSkpKyPhVV12lbt26qaioyNq3f/9+HTx4UB6PR5Lk8Xi0Z88eVVdXWzWFhYVyu91KS0uzas4+RkNNwzEAAEDXFtGnkO655x6tXr1ar7/+ulJTU639cXFx1pWRWbNm6be//a3y8vLkdrs1e/aZ2yRv2bJF0pmPUV955ZXq37+/Fi9eLJ/Pp6lTp2rGjBn60Y9+JOnMx6iHDBkir9erO++8U8XFxbr33ntVUFCgzMzMJs2VTyEBAGA/TX79NhGQFPaxatUqq+bkyZPmnnvuMRdffLHp2bOn+c///E9z+PDhkON8+OGH5oYbbjA9evQwffr0MT/4wQ9MXV1dSM2mTZvMlVdeaZxOp7n00ktDvkZT+P1+I8n4/f6IngcAANpPU1+/W3QfmI6MKzAAANhPm9wHBgAAoD102l/m2HBhiRvaAQBgHw2v25/3BlGnDTA1NTWSxA3tAACwoZqaGsXFxTU63mnXwASDQVVWVio2NlYOh6PVjhsIBJScnKxDhw51qbU19E3fXQF903dX0NH7NsaopqZG/fv3V1RU4ytdOu0VmKioqJCb7LU2t9vdIU/8hUbfXQt9dy303bV05L7Pd+WlAYt4AQCA7RBgAACA7RBgIuRyufTQQw91ud+7RN/03RXQN313BZ2l7067iBcAAHReXIEBAAC2Q4ABAAC2Q4ABAAC2Q4ABAAC2Q4CRtGzZMl1yySXq3r27Ro0apW3btp23/pVXXtFll12m7t27a+jQofrtb38bMm6M0YMPPqh+/fqpR48eSk9P14EDBy5kC80SSd8///nPdd111+niiy/WxRdfrPT09HPq77jjDjkcjpDHxIkTL3QbEYuk77y8vHN66t69e0hNZzzf48aNO6dvh8OhrKwsq6ajn+/Nmzfrm9/8pvr37y+Hw6HXXnvtc59TUlKiESNGyOVyafDgwcrLyzunJtJ/L9papH3/5je/0YQJE9S3b1+53W55PB5t3LgxpObhhx8+51xfdtllF7CLyEXad0lJSdjvcZ/PF1LX2c53uJ9bh8OhK664wqqxw/mWCDB6+eWXNW/ePD300EP605/+pOHDhyszM1PV1dVh67ds2aLvfOc7mj59unbu3KlJkyZp0qRJKi8vt2oWL16sp556SitXrtTWrVvVq1cvZWZm6tNPP22rtj5XpH2XlJToO9/5jjZt2qSysjIlJycrIyNDf//730PqJk6cqMOHD1uPl156qS3aabJI+5bO3K3y7J4++uijkPHOeL5/85vfhPRcXl6u6Oho/dd//VdIXUc+3ydOnNDw4cO1bNmyJtVXVFQoKytL48eP165duzRnzhzNmDEj5MW8Od8/bS3Svjdv3qwJEybot7/9rXbs2KHx48frm9/8pnbu3BlSd8UVV4Sc67feeutCTL/ZIu27wf79+0P6SkhIsMY64/leunRpSL+HDh1SfHz8OT/bHf18S5JMF3f11Vcbr9drbdfX15v+/fub3NzcsPU333yzycrKCtk3atQo8/3vf98YY0wwGDRJSUnm8ccft8aPHTtmXC6Xeemlly5AB80Tad+fdfr0aRMbG2uef/55a9+0adPMTTfd1NpTbVWR9r1q1SoTFxfX6PG6yvlesmSJiY2NNcePH7f22eF8N5Bk8vPzz1uzYMECc8UVV4Tsu+WWW0xmZqa13dK/x7bWlL7DSUtLM4888oi1/dBDD5nhw4e33sQusKb0vWnTJiPJ/POf/2y0piuc7/z8fONwOMyHH35o7bPL+e7SV2BOnTqlHTt2KD093doXFRWl9PR0lZWVhX1OWVlZSL0kZWZmWvUVFRXy+XwhNXFxcRo1alSjx2xrzen7sz755BPV1dUpPj4+ZH9JSYkSEhKUmpqqWbNm6R//+Eerzr0lmtv38ePHNWjQICUnJ+umm27S3r17rbGucr6fe+45TZkyRb169QrZ35HPd6Q+72e7Nf4e7SAYDKqmpuacn+0DBw6of//+uvTSS3Xrrbfq4MGD7TTD1nXllVeqX79+mjBhgt5++21rf1c5388995zS09M1aNCgkP12ON9dOsB8/PHHqq+vV2JiYsj+xMTEc94HbeDz+c5b3/BnJMdsa83p+7NycnLUv3//kB/uiRMn6le/+pWKior04x//WKWlpbrhhhtUX1/fqvNvrub0nZqaql/+8pd6/fXX9cILLygYDGr06NH629/+JqlrnO9t27apvLxcM2bMCNnf0c93pBr72Q4EAjp58mSr/NzYwRNPPKHjx4/r5ptvtvaNGjVKeXl52rBhg1asWKGKigpdd911qqmpaceZtky/fv20cuVKvfrqq3r11VeVnJyscePG6U9/+pOk1vl3sqOrrKzU+vXrz/nZtsv57rS/jRoXzqJFi7RmzRqVlJSELGidMmWK9d9Dhw7VsGHD9MUvflElJSW6/vrr22OqLebxeOTxeKzt0aNH6/LLL9fPfvYzPfbYY+04s7bz3HPPaejQobr66qtD9nfG893VrV69Wo888ohef/31kLUgN9xwg/Xfw4YN06hRozRo0CCtXbtW06dPb4+ptlhqaqpSU1Ot7dGjR+v999/XkiVL9H//93/tOLO28/zzz6t3796aNGlSyH67nO8ufQWmT58+io6OVlVVVcj+qqoqJSUlhX1OUlLSeesb/ozkmG2tOX03eOKJJ7Ro0SL97ne/07Bhw85be+mll6pPnz567733Wjzn1tCSvht069ZNX/nKV6yeOvv5PnHihNasWdOkf7Q62vmOVGM/2263Wz169GiV75+ObM2aNZoxY4bWrl17zltpn9W7d299+ctftu25bszVV19t9dTZz7cxRr/85S81depUOZ3O89Z21PPdpQOM0+nUVVddpaKiImtfMBhUUVFRyP91n83j8YTUS1JhYaFVn5KSoqSkpJCaQCCgrVu3NnrMttacvqUzn7Z57LHHtGHDBn31q1/93K/zt7/9Tf/4xz/Ur1+/Vpl3SzW377PV19drz549Vk+d+XxLZ24ZUFtbq9tuu+1zv05HO9+R+ryf7db4/umoXnrpJX3ve9/TSy+9FPJR+cYcP35c77//vm3PdWN27dpl9dSZz7cklZaW6r333mvS/5x02PPd3quI29uaNWuMy+UyeXl55t133zV333236d27t/H5fMYYY6ZOnWruu+8+q/7tt982MTEx5oknnjB/+ctfzEMPPWS6detm9uzZY9UsWrTI9O7d27z++utm9+7d5qabbjIpKSnm5MmTbd5fYyLte9GiRcbpdJpf//rX5vDhw9ajpqbGGGNMTU2N+eEPf2jKyspMRUWF+f3vf29GjBhhvvSlL5lPP/20XXoMJ9K+H3nkEbNx40bz/vvvmx07dpgpU6aY7t27m71791o1nfF8N/ja175mbrnllnP22+F819TUmJ07d5qdO3caSebJJ580O3fuNB999JExxpj77rvPTJ061ar/4IMPTM+ePc38+fPNX/7yF7Ns2TITHR1tNmzYYNV83t9jRxBp3y+++KKJiYkxy5YtC/nZPnbsmFXzgx/8wJSUlJiKigrz9ttvm/T0dNOnTx9TXV3d5v01JtK+lyxZYl577TVz4MABs2fPHvPf//3fJioqyvz+97+3ajrj+W5w2223mVGjRoU9ph3OtzHGdPkAY4wxTz/9tBk4cKBxOp3m6quvNu+88441NnbsWDNt2rSQ+rVr15ovf/nLxul0miuuuMIUFBSEjAeDQfPAAw+YxMRE43K5zPXXX2/279/fFq1EJJK+Bw0aZCSd83jooYeMMcZ88sknJiMjw/Tt29d069bNDBo0yNx1110d6ge9QSR9z5kzx6pNTEw03/jGN8yf/vSnkON1xvNtjDH79u0zkszvfve7c45lh/Pd8DHZzz4a+pw2bZoZO3bsOc+58sorjdPpNJdeeqlZtWrVOcc9399jRxBp32PHjj1vvTFnPk7er18/43Q6zX/8x3+YW265xbz33ntt29jniLTvH//4x+aLX/yi6d69u4mPjzfjxo0zxcXF5xy3s51vY87c6qFHjx7m2WefDXtMO5xvY4xxGGPMBb7IAwAA0Kq69BoYAABgTwQYAABgOwQYAABgOwQYAABgOwQYAABgOwQYAABgOwQYAABgOwQYAABgOwQYAABgOwQYAABgOwQYAABgOwQYAABgO/8fT8g/Pg6Om+8AAAAASUVORK5CYII=",
+ "text/plain": [
+ "<Figure size 640x480 with 1 Axes>"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "print(\"Local minimizer in area of pos. curvature: kappa =\", kappamin_pos, \", alpha =\", alphamin_pos, \", Q =\", Emin_pos)\n",
+ "print(\"Local minimizer in area of neg. curvature: kappa =\", kappamin_neg, \", alpha =\", alphamin_neg, \", Q =\", Emin_neg)\n",
+ "\n",
+ "n=10\n",
+ "bending_path_pos=np.outer(np.array(np.linspace(0,1,n)),np.array([kappamin_pos,0]))\n",
+ "bending_path_neg=np.outer(np.array(np.linspace(0,1,n)),np.array([-kappamin_neg,np.pi/2]))\n",
+ "\n",
+ "for i in range(0,n):\n",
+ " plt.plot(bending_path_pos[i,0],energy(bending_path_pos[i,0],bending_path_pos[i,1],Q,B), '-'\n",
+ " ) "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 20,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "array([[0. , 0. ],\n",
+ " [0.31084391, 0.52324901],\n",
+ " [0.62168782, 1.04649802],\n",
+ " [0.93253173, 1.56974703],\n",
+ " [1.24337564, 2.09299604],\n",
+ " [1.55421955, 2.61624505],\n",
+ " [1.86506346, 3.13949406],\n",
+ " [2.17590737, 3.66274307],\n",
+ " [2.48675128, 4.18599208],\n",
+ " [2.79759519, 4.70924109]])"
+ ]
+ },
+ "execution_count": 20,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "bending_path1"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "base",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.10.10"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/experiment/wood-bilayer_PLOS/perforated_wood_lower.py b/experiment/wood-bilayer_PLOS/perforated_wood_lower.py
index 3627cf878b980eec37461757cfa50b5dd6cf7c6c..69e8ffd690a0171f570b6155da081d1c5a286127 100644
--- a/experiment/wood-bilayer_PLOS/perforated_wood_lower.py
+++ b/experiment/wood-bilayer_PLOS/perforated_wood_lower.py
@@ -230,7 +230,7 @@ parameterSet.gamma=1.0
## numLevels : Number of Levels on which solution is computed. starting with a 2x2x2 cube mesh.
## {start,finish} computes on all grid from 2^(start) to 2^finish refinement
#----------------------------------------------------
-parameterSet.numLevels= '3 3' # computes all levels from first to second entry
+parameterSet.numLevels= '4 4' # computes all levels from first to second entry
# parameterSet.numLevels= '4 4'
#############################################
diff --git a/experiment/wood-bilayer_PLOS/perforated_wood_upper.py b/experiment/wood-bilayer_PLOS/perforated_wood_upper.py
deleted file mode 100644
index 06be5119e77f13d37f0950bce810822c59fded37..0000000000000000000000000000000000000000
--- a/experiment/wood-bilayer_PLOS/perforated_wood_upper.py
+++ /dev/null
@@ -1,280 +0,0 @@
-import math
-#from python_matrix_operations import *
-import ctypes
-import os
-import sys
-import numpy as np
-# import elasticity_toolbox as elast
-
-class ParameterSet(dict):
- def __init__(self, *args, **kwargs):
- super(ParameterSet, self).__init__(*args, **kwargs)
- self.__dict__ = self
-
-parameterSet = ParameterSet()
-#---------------------------------------------------------------
-#############################################
-# Paths
-#############################################
-# Path for results and logfile
-parameterSet.outputPath='/home/klaus/Desktop/Dune_release/dune-microstructure/experiment/perforated-bilayer/results'
-parameterSet.baseName= 'perforated_wood_upper' #(needed for Output-Filename)
-
-# Path for material description
-# parameterSet.geometryFunctionPath =experiment/wood-bilayer/
-
-#---------------------------------------------------------------
-# Wooden bilayer, https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6191116/#pone.0205607.ref015
-#--- define indicator function
-# x[0] : y1-component -1/2 to 1/2
-# x[1] : y2-component -1/2 to 1/2
-# x[2] : x3-component range -1/2 to 1/2
-#--- define indicator function
-def indicatorFunction(x):
- pRadius = math.sqrt((param_beta*param_r)/(np.pi*perfDepth)) # perforation radius
- if (x[2]>=(0.5-param_r)):
- if(((x[0]**2 + x[1]**2) < pRadius**2) and (x[2] >= (0.5-perfDepth))): #inside perforation
- return 3 #Phase3
- else:
- return 1 #Phase1
- else :
- return 2 #Phase2
-
-# def indicatorFunction(x):
-# factor=1
-# pRadius = 0.25
-# if (x[2]>=(0.5-param_r) and np.sqrt(x[0]**2 + x[1]**2) < pRadius):
-# return 3
-# elif((x[2]>=(0.5-param_r))):
-# return 1 #Phase1
-# else :
-# return 2 #Phase2
-
-# # --- Number of material phases
-# parameterSet.Phases=3
-
-# def indicatorFunction(x):
-# factor=1
-# pRadius = 1
-# # if (np.sqrt(x[0]*x[0] + x[1]*x[1]) < pRadius):
-# if ((x[0] < 0 and math.sqrt(pow(x[1],2) + pow(x[2],2) ) < pRadius/4.0) or ( 0 < x[0] and math.sqrt(pow(x[1],2) + pow(x[2],2) ) > pRadius/4.0)):
-# return 1
-# else :
-# return 2 #Phase2
-
-# --- Number of material phases
-parameterSet.Phases=3
-
-
-# Parameters of the model
-# -- (thickness upper layer) / (thickness)
-# param_r = 0.22
-param_r = 0.49
-# -- thickness [meter]
-param_h = 0.008
-# -- moisture content in the flat state [%]
-param_omega_flat = 17.17547062
-# -- moisture content in the target state [%]
-param_omega_target = 8.959564147
-# -- Drehwinkel
-param_theta = 0.0
-
-# Design Parameter ratio between perforaton (cylindrical) volume and volume of upper layer
-param_beta = 0.3
-# Depth of perforation
-# perfDepth = 0.12
-perfDepth = param_r
-# perfDepth = param_r * (2.0/3.0)
-#
-#
-#
-# -- increment of the moisture content
-delta_omega=param_omega_target-param_omega_flat
-# moisture content for material law
-omega=param_omega_target
-
-# --- Material properties from https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6191116/#pone.0205607.ref015
-# --- for European beech, moisture content omega = 15%
-# --- L=direction orthogonal to layering and fibres = orthogonal to wood stem cross-section
-# --- T=tangential zu layering
-# --- R=orthogonal zu layering
-# --- in MPa
-# --- Properties are defined by affine function in dependence of moisture content omega via property = b_0+b_1 \omega
-# --- coefficients of affine function are contained in the following array
-# --- data taken from http://dx.doi.org/10.1016/j.cma.2014.10.031
-
-properties_coefficients=np.array([
- # [b_0, b_1]
- [2565.6,-59.7], # E_R [MPa]
- [885.4, -23.4], # E_T [MPa]
- [17136.7,-282.4], # E_L [MPa]
- [667.8, -15.19], # G_RT [MPa]
- [1482, -15.26], # G_RL [MPa]
- [1100, -17.72], # G_TL [MPa]
- [0.2933, -0.001012], # nu_TR [1]
- [0.383, -0.008722], # nu_LR [1]
- [0.3368, -0.009071] # nu_LT [1]
- ])
-# Compute actual material properties
-E_R = properties_coefficients[0,0]+properties_coefficients[0,1]*omega
-E_T = properties_coefficients[1,0]+properties_coefficients[1,1]*omega
-E_L = properties_coefficients[2,0]+properties_coefficients[2,1]*omega
-G_RT = properties_coefficients[3,0]+properties_coefficients[3,1]*omega
-G_LR = properties_coefficients[4,0]+properties_coefficients[4,1]*omega
-G_LT = properties_coefficients[5,0]+properties_coefficients[5,1]*omega
-nu_TR = properties_coefficients[6,0]+properties_coefficients[6,1]*omega
-nu_LR = properties_coefficients[7,0]+properties_coefficients[7,1]*omega
-nu_LT = properties_coefficients[8,0]+properties_coefficients[8,1]*omega
-# Compute the remaining Poisson ratios
-nu_TL=nu_LT*E_T/E_L
-nu_RT=nu_TR*E_R/E_T
-nu_RL=nu_LR*E_R/E_L
-#
-# --- differential swelling strain
-# --- relation to swelling strain eps: eps=alpha* delta_omega with delta_omega = change of water content in %
-alpha_L=0.00011 # PLOS paper
-alpha_R=0.00191 # PLOS paper
-alpha_T=0.00462 # PLOS paper
-# Umrechnen
-#alpha_L=(1-1/(1+delta_omega*alpha_L))/delta_omega
-#alpha_R=(1-1/(1+delta_omega*alpha_R))/delta_omega
-#alpha_T=(1-1/(1+delta_omega*alpha_T))/delta_omega
-# --- define geometry
-
-
-
-# --- PHASE 1
-# y_1-direction: L
-# y_2-direction: T
-# x_3-direction: R
-# phase1_type="orthotropic"
-# materialParameters_phase1 = [E_L,E_T,E_R,G_TL,G_RT,G_RL,nu_LT,nu_LR,nu_TR]
-parameterSet.phase1_type="general_anisotropic"
-[E_1,E_2,E_3]=[E_L,E_T,E_R]
-[nu_12,nu_13,nu_23]=[nu_LT,nu_LR,nu_TR]
-[nu_21,nu_31,nu_32]=[nu_TL,nu_RL,nu_RT]
-[G_12,G_31,G_23]=[G_LT,G_LR,G_RT]
-compliance_S=np.array([[1/E_1, -nu_21/E_2, -nu_31/E_3, 0.0, 0.0, 0.0],
- [-nu_12/E_1, 1/E_2, -nu_32/E_3, 0.0, 0.0, 0.0],
- [-nu_13/E_1, -nu_23/E_2, 1/E_3, 0.0, 0.0, 0.0],
- [0.0, 0.0, 0.0, 1/G_23, 0.0, 0.0],
- [0.0, 0.0, 0.0, 0.0, 1/G_31, 0.0],
- [0.0, 0.0, 0.0, 0.0, 0.0, 1/G_12]]);
-materialParameters_phase1 = compliance_S
-
-def prestrain_phase1(x):
- # hB=delta_omega * alpha with delta_omega increment of moisture content and alpha swelling factor.
- return [[1/param_h*delta_omega*alpha_L, 0, 0], [0,1/param_h*delta_omega*alpha_T,0], [0,0,1/param_h*delta_omega*alpha_R]]
-
-# --- PHASE 2
-# y_1-direction: R
-# y_2-direction: L
-# x_3-direction: T
-parameterSet.phase2_type="general_anisotropic"
-[E_1,E_2,E_3]=[E_R,E_L,E_T]
-[nu_12,nu_13,nu_23]=[nu_RL,nu_RT,nu_LT]
-[nu_21,nu_31,nu_32]=[nu_LR,nu_TR,nu_TL]
-[G_12,G_31,G_23]=[G_LR,G_RT,G_LT]
-compliance_S=np.array([[1/E_1, -nu_21/E_2, -nu_31/E_3, 0.0, 0.0, 0.0],
- [-nu_12/E_1, 1/E_2, -nu_32/E_3, 0.0, 0.0, 0.0],
- [-nu_13/E_1, -nu_23/E_2, 1/E_3, 0.0, 0.0, 0.0],
- [0.0, 0.0, 0.0, 1/G_23, 0.0, 0.0],
- [0.0, 0.0, 0.0, 0.0, 1/G_31, 0.0],
- [0.0, 0.0, 0.0, 0.0, 0.0, 1/G_12]]);
-materialParameters_phase2 = compliance_S
-def prestrain_phase2(x):
- return [[1/param_h*delta_omega*alpha_R, 0, 0], [0,1/param_h*delta_omega*alpha_L,0], [0,0,1/param_h*delta_omega*alpha_T]]
-
-#Rotation um 2. Achse (= L)
-parameterSet.phase2_axis = 1
-# phase2_angle = param_theta
-# -- Drehwinkel
-parameterSet.phase2_angle = param_theta
-
-
-# --- PHASE 3
-parameterSet.phase3_type="isotropic"
-epsilon = 1e-8
-materialParameters_phase3 = [epsilon, epsilon]
-
-def prestrain_phase3(x):
- return [[0, 0, 0], [0,0,0], [0,0,0]]
-
-# # --- PHASE 3 = Phase 1 gedreht
-# # y_1-direction: L
-# # y_2-direction: R
-# # x_3-direction: T
-# parameterSet.phase3_type="general_anisotropic"
-# # Drehung um theta um Achse 2 = x_3-Achse
-# N=elast.rotation_matrix_compliance(2,param_theta)
-# materialParameters_phase3 = np.dot(np.dot(N,materialParameters_phase1),N.T)
-# materialParameters_phase3 = 0.5*(materialParameters_phase3.T+materialParameters_phase3)
-# # rotation of strain
-# def prestrain_phase3(x):
-# return elast.voigt_to_strain(np.dot(elast.rotation_matrix_compliance(2,param_theta),np.dot(elast.strain_to_voigt(np.array(prestrain_phase1(x))),N.T))).tolist()
-
-
-
-# --- Choose scale ratio gamma:
-parameterSet.gamma=1.0
-
-
-
-
-#############################################
-# Grid parameters
-#############################################
-## numLevels : Number of Levels on which solution is computed. starting with a 2x2x2 cube mesh.
-## {start,finish} computes on all grid from 2^(start) to 2^finish refinement
-#----------------------------------------------------
-parameterSet.numLevels= '3 3' # computes all levels from first to second entry
-# parameterSet.numLevels= '4 4'
-
-#############################################
-# Assembly options
-#############################################
-parameterSet.set_IntegralZero = 1 #(default = false)
-parameterSet.set_oneBasisFunction_Zero = 1 #(default = false)
-#parameterSet.arbitraryLocalIndex = 7 #(default = 0)
-#parameterSet.arbitraryElementNumber = 3 #(default = 0)
-
-#############################################
-# Solver Options, Type: #1: CG - SOLVER , #2: GMRES - SOLVER, #3: QR - SOLVER (default), #4: UMFPACK - SOLVER
-#############################################
-parameterSet.Solvertype = 3 # recommended to use iterative solver (e.g GMRES) for finer grid-levels
-parameterSet.Solver_verbosity = 0 #(default = 2) degree of information for solver output
-
-
-#############################################
-# Write/Output options #(default=false)
-#############################################
-# --- (Optional output) write Material / prestrain / Corrector functions to .vtk-Files:
-parameterSet.write_materialFunctions = 1 # VTK indicator function for material/prestrain definition
-#parameterSet.write_prestrainFunctions = 1 # VTK norm of B (currently not implemented)
-
-# --- (Additional debug output)
-parameterSet.print_debug = 0 #(default=false)
-
-# --- Write Correctos to VTK-File:
-parameterSet.write_VTK = 1
-
-# The grid can be refined several times for a higher resolution in the VTK-file.
-parameterSet.subsamplingRefinement = 2
-
-# --- (Optional output) L2Error, integral mean:
-#parameterSet.write_L2Error = 1
-#parameterSet.write_IntegralMean = 1
-
-# --- check orthogonality (75) from paper:
-parameterSet.write_checkOrthogonality = 0
-
-# --- Write corrector-coefficients to log-File:
-#parameterSet.write_corrector_phi1 = 1
-#parameterSet.write_corrector_phi2 = 1
-#parameterSet.write_corrector_phi3 = 1
-
-# --- Print Condition number of matrix (can be expensive):
-#parameterSet.print_conditionNumber= 1 #(default=false)
-
-# --- write effective quantities to Matlab-folder for symbolic minimization:
-parameterSet.write_toMATLAB = 1 # writes effective quantities to .txt-files QMatrix.txt and BMatrix.txt
diff --git a/experiment/wood-bilayer_PLOS/wood_bilayer_test.py b/experiment/wood-bilayer_PLOS/wood_bilayer_test.py
index e319d8a4784a6afcf5c3a77e052684516a8ee78a..3226fe34cc57be012a3e392b0bfddf293bcfa184 100644
--- a/experiment/wood-bilayer_PLOS/wood_bilayer_test.py
+++ b/experiment/wood-bilayer_PLOS/wood_bilayer_test.py
@@ -159,12 +159,12 @@ for dataset_number in dataset_numbers:
[0.17, 0.0049, 17.28772791 , 9.101671742, 0.0, 0 ]
],
[ # Dataset Ratio r = 0.22
- [0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0 ],
- [0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0 ],
- [0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0 ],
- [0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0 ],
- [0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0 ],
- [0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0 ],
+ [0.22, 0.0053, 17.17547062, 14.72680026, 0.0, 0 ],
+ [0.22, 0.0053, 17.17547062, 13.64338887, 0.0, 0 ],
+ [0.22, 0.0053, 17.17547062, 12.41305478, 0.0, 0 ],
+ [0.22, 0.0053, 17.17547062, 11.66482931, 0.0, 0 ],
+ [0.22, 0.0053, 17.17547062, 11.09781471, 0.0, 0 ],
+ [0.22, 0.0053, 17.17547062, 9.435795985, 0.0, 0 ],
[0.22, 0.0053, 17.17547062, 8.959564147, 0.0, 0 ]
],
[ # Dataset Ratio r = 0.34