From de58091144464436bec495939d9d8467e2ce115d Mon Sep 17 00:00:00 2001 From: Kuan-Fu Feng Date: Thu, 22 Feb 2024 16:27:19 -0800 Subject: [PATCH 1/3] update attenuation notebook to use StackStore --- src/noisepy/monitoring/esyn_plotting.py | 127 +++++++++ src/noisepy/monitoring/esyn_utils.py | 25 +- .../attenuation_singlestation.ipynb | 242 ++++++------------ 3 files changed, 207 insertions(+), 187 deletions(-) create mode 100644 src/noisepy/monitoring/esyn_plotting.py diff --git a/src/noisepy/monitoring/esyn_plotting.py b/src/noisepy/monitoring/esyn_plotting.py new file mode 100644 index 00000000..5dbfab3f --- /dev/null +++ b/src/noisepy/monitoring/esyn_plotting.py @@ -0,0 +1,127 @@ +import matplotlib.pyplot as plt +import numpy as np + + +def plot_waveforms(ncmp, wav, fname, comp_arr): + fig, ax = plt.subplots(1, ncmp, figsize=(16, 3), sharex=False) + + for n in range(ncmp): + absy = max(wav[n][1], key=abs) + ax[n].set_ylim(absy * -1, absy) + ax[n].plot(wav[n][0], wav[n][1]) + ax[n].set_xlabel("time [s]") + ax[n].set_title(fname + " " + comp_arr[n]) + fig.tight_layout() + # print("save figure as Waveform_readin_%s.png"%(fname)) + plt.savefig("Waveform_readin_%s.png" % (fname), format="png", dpi=100) + plt.close(fig) + + +def plot_filtered_waveforms(freq, tt, wav, fname, ccomp): + nfreq = len(freq) - 1 + fig, ax = plt.subplots(1, nfreq, figsize=(16, 3), sharex=False) + + for fb in range(nfreq): + fmin = freq[fb] + fmax = freq[fb + 1] + absy = max(wav[fb], key=abs) + ax[fb].set_ylim(absy * -1, absy) + ax[fb].plot(tt, wav[fb], "k-", linewidth=0.2) + ax[fb].set_xlabel("Time [s]") + ax[fb].set_ylabel("Amplitude") + ax[fb].set_title("%s %s @%4.2f-%4.2f Hz" % (fname, ccomp, fmin, fmax)) + fig.tight_layout() + plt.savefig("Waveform_filtered_%s_%s_F%s-%s.png" % (fname, ccomp, fmin, fmax), format="png", dpi=100) + plt.close(fig) + + +def plot_envelope(comp_arr, freq, msv, msv_mean, fname, vdist): + nfreq = len(freq) - 1 + ncmp = len(comp_arr) + + fig, ax = plt.subplots(ncmp + 1, nfreq, figsize=(16, 10), sharex=False) + for n in range(len(comp_arr)): + for fb in range(nfreq): + fmin = freq[fb] + fmax = freq[fb + 1] + ax[n, fb].plot(msv[n][0][:], msv[n][fb + 1], "k-", linewidth=0.5) + ax[n, fb].set_title("%s %.2fkm %s @%4.2f-%4.2f Hz" % (fname, vdist, comp_arr[n], fmin, fmax)) + ax[n, fb].set_xlabel("Time [s]") + ax[n, fb].set_ylabel("Amplitude") + + for fb in range(nfreq): + fmin = freq[fb] + fmax = freq[fb + 1] + ax[-1, fb].plot(msv_mean[0], msv_mean[fb + 1], "b-", linewidth=1) + ax[-1, fb].set_title(" Mean Squared Value %.2fkm @%4.2f-%4.2f Hz" % (vdist, fmin, fmax)) + ax[-1, fb].set_xlabel("Time [s]") + ax[-1, fb].set_ylabel("Amplitude") + plt.tight_layout() + plt.savefig("Waveform_envelope_%s_F%s-%s.png" % (fname, fmin, fmax), format="png", dpi=100) + plt.close(fig) + + +def plot_fmsv_waveforms(freq, wav, fname, noise_level, twin): + nfreq = len(freq) - 1 + fig, ax = plt.subplots(1, nfreq, figsize=(16, 3), sharex=False) + + for fb in range(nfreq): + fmin = freq[fb] + fmax = freq[fb + 1] + absy = 1 # max(wav[fb], key=abs) + ax[fb].plot( + [wav[0][0], wav[0][-1]], [noise_level[fb], noise_level[fb]], c="blue", marker=".", ls="--", linewidth=2 + ) + + ax[fb].plot([twin[fb][0], twin[fb][0]], [-0.1, absy], c="orange", marker=".", ls="--", linewidth=2) + ax[fb].plot([twin[fb][1], twin[fb][1]], [-0.1, absy], c="orange", marker=".", ls="--", linewidth=2) + ax[fb].set_yscale("log", base=10) + ax[fb].plot(wav[0], wav[fb + 1], "k-", linewidth=0.5) + ax[fb].set_xlabel("Time [s]") + ax[fb].set_ylabel("Amplitude in log-scale") + ax[fb].set_title("%s @%4.2f-%4.2f Hz" % (fname, fmin, fmax)) + fig.tight_layout() + plt.savefig("Waveform_fmsv_%s.png" % (fname), format="png", dpi=100) + plt.close(fig) + + +def plot_fitting_curves(mean_free, intrinsic_b, tt, Eobs, Esyn, fname, dist, twin, fmin, fmax): + numb = len(intrinsic_b) + plt.figure(figsize=(8, 2)) + for nb in range(numb): + plt.yscale("log", base=10) + # plt.xlim(0,120) + pymin = np.min(Eobs[nb][:-2] / 2) + pymax = np.max(Eobs[nb][:-2] * 2) + plt.ylim(pymin, pymax) + plt.plot(tt, Eobs[nb], "k-", linewidth=0.5) + plt.plot(tt, Esyn[nb], "b-", linewidth=1) + plt.plot([twin[0], twin[0], twin[-1], twin[-1], twin[0]], [pymin, pymax, pymax, pymin, pymin], "r", linewidth=2) + + plt.title( + "%s %.2fkm @%4.2f-%4.2f Hz, mean_free: %.2f b: %.2f~%.2f" + % (fname, dist, fmin, fmax, mean_free, intrinsic_b[0], intrinsic_b[-1]) + ) + plt.xlabel("Time [s]") + plt.ylabel("Energy density Amplitude") + plt.tight_layout() + plt.savefig("Fitting_fmsv_%s_F%s-%s_MFP%.2f.png" % (fname, fmin, fmax, mean_free), format="png", dpi=100) + plt.close() + + +def plot_fitting_result(mean_free, intrinsic_b, tt, Eobs, Esyn, fname, dist, twin, fmin, fmax): + plt.figure(figsize=(6, 2)) + plt.yscale("log", base=10) + + pymax = np.max(Eobs[:-2] * 5) + pymin = 10 ** (-6) + plt.ylim(pymin, pymax) + plt.plot(tt, Eobs, "k-", linewidth=1) + plt.plot(tt, Esyn, "b--", linewidth=1) + plt.plot([twin[0], twin[0], twin[-1], twin[-1], twin[0]], [pymin, pymax, pymax, pymin, pymin], "r", linewidth=2) + + plt.title("%s %.2fkm @%4.2f-%4.2f Hz, intrinsic b: %.2f" % (fname, dist, fmin, fmax, intrinsic_b)) + plt.xlabel("Time [s]") + plt.ylabel("Energy density Amp") + plt.tight_layout() + plt.savefig("Final_fmsv_%s_F%s-%s.png" % (fname, fmin, fmax), format="png", dpi=100) diff --git a/src/noisepy/monitoring/esyn_utils.py b/src/noisepy/monitoring/esyn_utils.py index bbe1cbfd..87107de8 100644 --- a/src/noisepy/monitoring/esyn_utils.py +++ b/src/noisepy/monitoring/esyn_utils.py @@ -3,7 +3,6 @@ from typing import Tuple import numpy as np -import pyasdf ### ----- # These scripts are aim to perform the 2-D radiative transfer equation @@ -15,28 +14,6 @@ logger = logging.getLogger(__name__) -### ----- -def read_pyasdf(sfile: str, ccomp: str) -> Tuple[float, float, np.ndarray, np.ndarray]: - # useful parameters from each asdf file - with pyasdf.ASDFDataSet(sfile, mode="r") as ds: - alist = ds.auxiliary_data.list() - try: - dt = ds.auxiliary_data[alist[0]][ccomp].parameters["dt"] - dist = ds.auxiliary_data[alist[0]][ccomp].parameters["dist"] - logger.info(f"working on {sfile} (comp: {ccomp}) that is {dist} km apart. dt: {dt}") - # read stacked data - sdata = ds.auxiliary_data[alist[0]][ccomp].data[:] - - # time domain variables - npts = sdata.size - tvec = np.arange(-npts // 2 + 1, npts // 2 + 1) * dt - return dist, dt, tvec, sdata - - except Exception: - logger.warning(f"continue! no {ccomp} component exist") - return None - - ### ----- # Function that Calculate Mean Square def msValue(arr: np.ndarray) -> np.ndarray: @@ -252,7 +229,7 @@ def get_SSR(fnum: int, para) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: return SSR_final, mfpx, intby -def get_optimal(fnum: int, para) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: +def get_optimal_Esyn(fnum: int, para) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: """ # Getting the optimal value from the grid searching results (the SSR output from the get_SSR) # Return with the optimal value of mean free path, intrinsic absorption parameter diff --git a/tutorials/monitoring/attenuation_singlestation.ipynb b/tutorials/monitoring/attenuation_singlestation.ipynb index c2008bb1..b691a5d9 100644 --- a/tutorials/monitoring/attenuation_singlestation.ipynb +++ b/tutorials/monitoring/attenuation_singlestation.ipynb @@ -5,7 +5,7 @@ "id": "58f7c398", "metadata": {}, "source": [ - "# Welcome to the Tutorial measuring intrinsic absorption parameter!\n", + "# Welcome to the tutorial measuring intrinsic absorption parameter!\n", "\n", "This Notebook calculate the scattering & intrinsic absorption parameters of the Rayleigh waves following the instruction proposed by Hirose et al. (2019).\n", "\n", @@ -30,20 +30,16 @@ "outputs": [], "source": [ "import os\n", - "import sys\n", "import glob\n", - "import obspy\n", "import numpy as np\n", "\n", - "import matplotlib.cm as cm\n", "import matplotlib.pyplot as plt\n", - "import pyasdf\n", - "import scipy\n", - "import math\n", "\n", "from obspy.signal.filter import bandpass\n", "from noisepy.monitoring.esyn_utils import *\n", - "from noisepy.seis.noise_module import mad" + "from noisepy.monitoring.esyn_plotting import plot_waveforms\n", + "from noisepy.seis.noise_module import mad\n", + "from noisepy.seis.asdfstore import ASDFCCStore, ASDFStackStore" ] }, { @@ -57,30 +53,36 @@ { "cell_type": "code", "execution_count": null, - "id": "ea154d07", + "id": "6f435225", "metadata": {}, "outputs": [], "source": [ - "# detremine the h5 file for processing\n", - "path=\"../example_data/\"\n", - "sfiles = sorted(glob.glob(os.path.join(path, 'CI.ADO_CI.ADO.h5')))\n", - "print(\"The h5 file: \",sfiles)\n", - "sta_pair=sfiles[0].split(path)[1].split(\".h5\")[0]\n", - "print(\"Processed station pair: \",sta_pair)" + "path = os.path.join(\"../\", \"get_started_data\")\n", + "\n", + "os.makedirs(path,exist_ok=True)\n", + "data_path = os.path.join(path, \"STACK/\")\n", + "print(os.listdir(data_path))\n", + "wave_store = ASDFStackStore(data_path)\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "42a5f808", + "id": "d243ac22", "metadata": {}, "outputs": [], "source": [ - "# detremine the component for processing\n", - "comp_arr = [\"EN\",\"EZ\",\"NZ\"] \n", - "num_cmp=len(comp_arr)\n", - "fnum=len(sfiles)\n", - "print(fnum, num_cmp )" + "pairs = wave_store.get_station_pairs()\n", + "print(f\"Found {len(pairs)} station pairs\")\n", + "print('pairs: ',pairs)\n", + "\n", + "stations = set(pair[0] for pair in pairs)\n", + "print('Stations: ', stations)\n", + "\n", + "sta_stacks = wave_store.read_bulk(None, pairs) # no timestamp used in ASDFStackStore\n", + "stacks = sta_stacks[0][1]\n", + "print(\"Target pair: \",sta_stacks[0])\n", + "target_pair=sta_stacks[0][0][0].network+\".\"+sta_stacks[0][0][0].name" ] }, { @@ -90,32 +92,31 @@ "metadata": {}, "outputs": [], "source": [ - "# information of the input data: lagtime (lag), sampling rate (samp)\n", - "lag=200 \n", - "samp=20\n", - "leng=int(lag*samp*2+1)\n", - "npts=leng\n", - "print(\"Lag-time: \",lag,\", sampling rate: \",samp,\", total data length in points: \",leng)" + "print(\"Processed station pair: \",sta_stacks[0][0])\n", + "\n", + "params = stacks[0].parameters\n", + "dt, lag , dist = (params[p] for p in [\"dt\", \"maxlag\", \"dist\"])\n", + "\n", + "npts = int(lag*(1/dt)*2+1)\n", + "tvec = np.arange(-npts // 2 + 1, npts // 2 + 1) * dt\n", + "print(\"Lag-time: \",lag,\", sampling rate: \",(1/dt) ,\", total data length in points: \",npts)" ] }, { "cell_type": "code", "execution_count": null, - "id": "c045328e", + "id": "42a5f808", "metadata": {}, "outputs": [], "source": [ - "def plot_waveforms(ncmp,wav,fname,comp_arr):\n", - " fig, ax = plt.subplots(1,ncmp, figsize=(16,3), sharex=False)\n", - " \n", - " for n in range(ncmp):\n", - " absy=max(wav[n][1], key=abs)\n", - " ax[n].set_ylim(absy*-1,absy)\n", - " ax[n].plot(wav[n][0],wav[n][1])\n", - " ax[n].set_xlabel(\"time [s]\")\n", - " ax[n].set_title(fname+\" \"+comp_arr[n])\n", - " fig.tight_layout()\n", - " plt.show()" + "# detremine the component for processing\n", + "comp_arr = [\"EN\",\"EZ\",\"NZ\"] \n", + "num_cmp=len(comp_arr)\n", + "fnum=len([target_pair])\n", + "\n", + "# define the array being used\n", + "stackf=np.ndarray((fnum, num_cmp, 2, npts)) \n", + "vdist=np.zeros((fnum,1)) # S-R distance array" ] }, { @@ -125,25 +126,21 @@ "metadata": {}, "outputs": [], "source": [ - "# Ready to read in h5 file\n", - "stackf=np.ndarray((fnum,num_cmp,2,leng)) \n", - "print(stackf.shape)\n", - "vdist=np.zeros((fnum,1)) # S-R distance array\n", - "fname=[] # file name array\n", - "\n", + "# file name array\n", + "fname=[] \n", + "stack_name=\"Allstack_linear\"\n", "aa=0\n", "# loop through each station-pair\n", - "for sfile in sfiles:\n", + "for sfile in [target_pair]:\n", " ncmp=0\n", - " sta_pair=sfile.split(path)[1].split(\".h5\")[0]\n", - " fname.append(sta_pair)\n", + " fname.append(sfile)\n", + " # read stacked waveforms accordingly to the cross-component\n", " for ccomp in comp_arr: \n", " print(aa, sfile, fname,ccomp)\n", - " # read stacked waveforms\n", - " if ( read_pyasdf(sfile,ccomp) == None):\n", - " continue\n", - " dist,dt, tvec,sdata = read_pyasdf(sfile,ccomp) # read waveform from pyasdf \n", - " stackf[aa][ncmp]=[tvec,sdata]\n", + " \n", + " stacks_data = list(filter(lambda x: x.component == ccomp, stacks))\n", + " print(stacks_data)\n", + " stackf[aa][ncmp]=[tvec, stacks_data[0].data]\n", " vdist[aa]=dist\n", " ncmp=ncmp+1\n", " plot_waveforms(num_cmp,stackf[aa],fname[aa],comp_arr)\n", @@ -152,43 +149,39 @@ "fnum=len(fname)\n" ] }, + { + "cell_type": "markdown", + "id": "12acf76a", + "metadata": {}, + "source": [ + "### Step 1 --- Calculation of mean-squared (MS) envelopes --> observed energy densities (***Eobs***)\n", + "--> normalized MS envelope is referred to as the observed energy density Eobs " + ] + }, { "cell_type": "code", "execution_count": null, - "id": "225e782a", + "id": "710ff292", "metadata": {}, "outputs": [], "source": [ - "def plot_filtered_waveforms(freq,tt,wav,fname):\n", + "def plot_filtered_waveforms(freq, tt, wav, fname, ccomp):\n", " nfreq = len(freq) - 1\n", " fig, ax = plt.subplots(1,nfreq, figsize=(16,3), sharex=False)\n", - " \n", + "\n", " for fb in range(nfreq):\n", " fmin=freq[fb]\n", " fmax=freq[fb+1]\n", " absy=max(wav[fb], key=abs)\n", - " #absx=max(tt, key=abs)\n", - " #ax[fb].set_xlim(absx*-1,absx)\n", " ax[fb].set_ylim(absy*-1,absy)\n", " ax[fb].plot(tt,wav[fb], \"k-\", linewidth=0.2)\n", - " #ax[fb].plot(wav_fold[0],wav_fold[fb+1], \"b-\", linewidth=1)\n", " ax[fb].set_xlabel(\"Time [s]\")\n", " ax[fb].set_ylabel(\"Amplitude\")\n", " ax[fb].set_title( \"%s %s @%4.2f-%4.2f Hz\" % ( fname,ccomp,fmin,fmax ) )\n", - " \n", - " fig.tight_layout()\n", + " plt.tight_layout()\n", " plt.show()" ] }, - { - "cell_type": "markdown", - "id": "12acf76a", - "metadata": {}, - "source": [ - "### Step 1 --- Calculation of mean-squared (MS) envelopes --> observed energy densities (***Eobs***)\n", - "--> normalized MS envelope is referred to as the observed energy density Eobs " - ] - }, { "cell_type": "code", "execution_count": null, @@ -196,7 +189,7 @@ "metadata": {}, "outputs": [], "source": [ - "freq = [0.5, 1, 2,4] # targeted frequency band for waveform monitoring\n", + "freq = [0.4, 0.8, 2,4] # targeted frequency band for waveform monitoring\n", "nfreq = len(freq) - 1\n", "\n", "MSE=np.ndarray((fnum,num_cmp,nfreq+1,npts)) # filtered two-side averaged stack CF\n", @@ -216,7 +209,7 @@ " dafbp[fb] = bandpass(data, fmin, fmax, int(1 / dt), corners=4, zerophase=True)\n", " \n", " MSE[aa][ncmp]=[stackf[aa][ncmp][0],dafbp[0],dafbp[1],dafbp[2]] \n", - " plot_filtered_waveforms(freq,stackf[aa][ncmp][0],dafbp,fname[aa])\n", + " plot_filtered_waveforms(freq, stackf[aa][ncmp][0], dafbp, fname[aa], ccomp)\n", "\n" ] }, @@ -227,14 +220,14 @@ "metadata": {}, "outputs": [], "source": [ - "def plot_envelope(comp_arr,freq,msv,msv_mean,fname,vdist):\n", - " \n", + "def plot_envelope(comp_arr, freq, msv, msv_mean, fname, vdist):\n", + "\n", " nfreq = len(freq) - 1\n", " ncmp = len(comp_arr)\n", - " \n", - " fig, ax = plt.subplots(ncmp+1,nfreq, figsize=(16,10), sharex=False) \n", - " for n in range(len(comp_arr)):\n", - " \n", + "\n", + " fig, ax = plt.subplots(ncmp+1, nfreq, figsize=(16,10), sharex=False) \n", + " for n in range(len(comp_arr)):\n", + "\n", " for fb in range(nfreq):\n", " fmin=freq[fb]\n", " fmax=freq[fb+1] \n", @@ -242,18 +235,17 @@ " ax[n,fb].set_title(\"%s %.2fkm %s @%4.2f-%4.2f Hz\" % (fname,vdist,comp_arr[n],fmin,fmax))\n", " ax[n,fb].set_xlabel(\"Time [s]\")\n", " ax[n,fb].set_ylabel(\"Amplitude\")\n", - " \n", + "\n", " for fb in range(nfreq):\n", " fmin=freq[fb]\n", " fmax=freq[fb+1]\n", " ax[-1,fb].plot(msv_mean[0], msv_mean[fb+1], \"b-\", linewidth=1)\n", " ax[-1,fb].set_title(\" Mean Squared Value %.2fkm @%4.2f-%4.2f Hz\" % (vdist,fmin,fmax))\n", " ax[-1,fb].set_xlabel(\"Time [s]\")\n", - " ax[-1,fb].set_ylabel(\"Amplitude\") \n", + " ax[-1,fb].set_ylabel(\"Amplitude\")\n", "\n", - " plt.tight_layout() \n", - " plt.show() \n", - " " + " plt.tight_layout()\n", + " plt.show()" ] }, { @@ -548,96 +540,20 @@ "\n", " # parameters for getting optimal value from the sum of squared residuals (SSR) between Eobs and Esyn \n", " para={ 'fb':fb, 'fmin':fmin, 'fmax':fmax, 'vdist':vdist, 'npts':npts, 'dt':dt, 'cvel':c, 'filenum':aa, \\\n", - " 'mfp':mfpx, 'intb':intby,'twin':twinbe, 'fmsv':fmsv_mean, 'SSR':SSR , 'sta':sta_pair}\n", + " 'mfp':mfpx, 'intb':intby,'twin':twinbe, 'fmsv':fmsv_mean, 'SSR':SSR , 'sta':target_pair}\n", " # call function get_optimal\n", - " result_intb[fb], result_mfp[fb], Eobs[fb], Esyn[fb] = get_optimal(fnum,para)\n", + " result_intb[fb], result_mfp[fb], Eobs[fb], Esyn[fb] = get_optimal_Esyn(fnum,para)\n", " \n", " plot_fitting_result(result_mfp[fb],result_intb[fb],fmsv_mean[aa][0][:],\n", - " Eobs[fb],Esyn[fb],sta_pair,vdist[aa],twinbe[aa][fb],fmin,fmax)\n" - ] - }, - { - "cell_type": "markdown", - "id": "d0f8ab81", - "metadata": {}, - "source": [ - "## Plotting SSR (for inter-station measurement)" + " Eobs[fb],Esyn[fb],target_pair,vdist[aa],twinbe[aa][fb],fmin,fmax)\n" ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "374cf19a", - "metadata": {}, - "outputs": [], - "source": [ - "def plot_grid_searching(sta_pair,freq,SSR):\n", - " nfreq=len(freq)-1\n", - " \n", - " fig, ax = plt.subplots(1,nfreq, figsize=(16,4), sharex=False)\n", - " \n", - " for fb in range(nfreq): \n", - " fmin=freq[fb]\n", - " fmax=freq[fb+1]\n", - "\n", - " loc=np.where(SSR[fb].T == np.amin(SSR[fb].T))\n", - " #locx=list(zip(loc[0], loc[1]))\n", - " print(\"%4.2f-%4.2f Hz \" % (fmin,fmax),\"loc \",loc)\n", - " ymin=y[loc[0]]\n", - " xmin=x[loc[1]]\n", - " print(\" intrinsic_b %.2f \" % ymin,\"mean_free: %.2f \" % xmin)\n", - " result_intb[fb]=ymin\n", - " result_mfp[fb]=xmin\n", - "\n", - " grid = SSR[fb].T\n", - " im=ax[fb].imshow(grid,extent=(x.min(), x.max(), y.max(), y.min()), aspect='auto',cmap = 'viridis_r',interpolation='spline16' )\n", - " #im=ax[fb].imshow(grid,aspect='auto',cmap = 'viridis_r' )\n", - " im.set_clim(1,1.3) \n", - " cb=plt.colorbar(im,extend='max')\n", - " cb.set_label('SSR/SSR_min', rotation=90, labelpad=14)\n", - " ax[fb].set_title(\"%s SSR @%4.2f-%4.2f Hz\" % (sta_pair,fmin,fmax) )\n", - " ax[fb].set_xlabel(\"mean free path\")\n", - " ax[fb].set_ylabel(\"intrinsic_b\")\n", - " ax[fb].invert_yaxis()\n", - " ax[fb].plot(xmin,ymin,\"+\", markersize=20, color='red')\n", - "\n", - " plt.tight_layout() \n", - " plt.show() \n", - " \n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "73654905", - "metadata": {}, - "outputs": [], - "source": [ - "#plot_grid_searching(sta_pair,freq,SSR)\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "a6b364c4", - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "6681804f", - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { "kernelspec": { "display_name": "noisepy", "language": "python", - "name": "python3" + "name": "noisepy" }, "language_info": { "codemirror_mode": { From 5704fe51a133be9a22d834e00fd190edea35f902 Mon Sep 17 00:00:00 2001 From: Kuan-Fu Feng Date: Mon, 11 Mar 2024 17:22:31 -0700 Subject: [PATCH 2/3] convert to seis-io and update some settings for monitoring script and notebook --- src/noisepy/monitoring/esyn_utils.py | 2 +- .../attenuation_singlestation.ipynb | 118 +++++++++--------- 2 files changed, 59 insertions(+), 61 deletions(-) diff --git a/src/noisepy/monitoring/esyn_utils.py b/src/noisepy/monitoring/esyn_utils.py index 87107de8..9ea93747 100644 --- a/src/noisepy/monitoring/esyn_utils.py +++ b/src/noisepy/monitoring/esyn_utils.py @@ -229,7 +229,7 @@ def get_SSR(fnum: int, para) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: return SSR_final, mfpx, intby -def get_optimal_Esyn(fnum: int, para) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: +def get_optimal_Esyn(fnum: int, para) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: """ # Getting the optimal value from the grid searching results (the SSR output from the get_SSR) # Return with the optimal value of mean free path, intrinsic absorption parameter diff --git a/tutorials/monitoring/attenuation_singlestation.ipynb b/tutorials/monitoring/attenuation_singlestation.ipynb index b691a5d9..1306b4ca 100644 --- a/tutorials/monitoring/attenuation_singlestation.ipynb +++ b/tutorials/monitoring/attenuation_singlestation.ipynb @@ -19,7 +19,7 @@ "Step:
\n", "0) Data preparing and filtering
1) Calculation of mean-squared (MS) envelopes --> observed energy densities (Eobs)
\n", "2) Calculation of synthesized energy densities (Esyn) via a grid search
\n", - "3) Determination of best-fit parameters: intrinsic absorption parameter *b* (for single station)
" + "3) Determination of best-fit parameters: intrinsic absorption parameter *b* and intrinsic *Q-value* (for single station)
" ] }, { @@ -30,16 +30,15 @@ "outputs": [], "source": [ "import os\n", - "import glob\n", "import numpy as np\n", "\n", "import matplotlib.pyplot as plt\n", "\n", "from obspy.signal.filter import bandpass\n", "from noisepy.monitoring.esyn_utils import *\n", - "from noisepy.monitoring.esyn_plotting import plot_waveforms\n", + "#from noisepy.monitoring.esyn_plotting import plot_waveforms\n", "from noisepy.seis.noise_module import mad\n", - "from noisepy.seis.asdfstore import ASDFCCStore, ASDFStackStore" + "from noisepy.seis.io.asdfstore import ASDFStackStore" ] }, { @@ -119,6 +118,28 @@ "vdist=np.zeros((fnum,1)) # S-R distance array" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "03f5bbd7", + "metadata": {}, + "outputs": [], + "source": [ + "def plot_waveforms(ncmp, wav, fname, comp_arr):\n", + " fig, ax = plt.subplots(1, ncmp, figsize=(16, 3), sharex=False)\n", + "\n", + " for n in range(ncmp):\n", + " absy = max(wav[n][1], key=abs)\n", + " ax[n].set_ylim(absy * -1, absy)\n", + " ax[n].plot(wav[n][0], wav[n][1])\n", + " ax[n].set_xlabel(\"time [s]\")\n", + " ax[n].set_title(fname + \" \" + comp_arr[n])\n", + " fig.tight_layout()\n", + " # print(\"save figure as Waveform_readin_%s.png\"%(fname))\n", + " plt.savefig(\"Waveform_readin_%s.png\" % (fname), format=\"png\", dpi=100)\n", + " plt.close(fig)" + ] + }, { "cell_type": "code", "execution_count": null, @@ -189,7 +210,7 @@ "metadata": {}, "outputs": [], "source": [ - "freq = [0.4, 0.8, 2,4] # targeted frequency band for waveform monitoring\n", + "freq = [0.6, 0.8, 2,4] # targeted frequency band for waveform monitoring\n", "nfreq = len(freq) - 1\n", "\n", "MSE=np.ndarray((fnum,num_cmp,nfreq+1,npts)) # filtered two-side averaged stack CF\n", @@ -360,56 +381,23 @@ "metadata": {}, "source": [ "### Step 2 --- Calculation of synthesized energy densities (***Esyn***) via a grid search \n", - "### #The 2-D radiative transfer equation for scalar waves ***(Shang and Gao 1988; Sato 1993)***\n", - "Assuming isotropic scattering and source radiation in infinite medium to calculate ***synthesized energy densities Esyn*** :" - ] - }, - { - "attachments": { - "image-2.png": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4kAAAB9CAYAAAAPxSbHAAABXWlDQ1BJQ0MgUHJvZmlsZQAAKJFtkL9LQnEUxY9lSSXYEEG/4A1NYiHPoGgzoQgaRAuqJZ5P0+CpX7++iLb+gaKhxraG/oAwiIbGpgiCwKZoqT2QqOR1vlqp1YXL+XC493C5QJvXEMJyA8jmbBmbm9GWV1Y1zxM6MYgueOAzzKIIR6MLHMG3tlblDi6lt2Mqy2+Uh97u/Y/XpeyetR96/jvfUt3JVNGkfrB1U0gbcAXJ0S1bKN4h90keRT5QnK7zieJEnS9qM4uxCPmG3GtmjCT5gRxINPnpJs5am+bXDep6byq3FKf2s4cRh4kUBGxIqoZZ6Jjkj/7fmajtRJDnxjY3NpBGhrsawnQErFrGPHJMHUeArCPIDqlf//5hw8tPA1NXQPthw0vwn2dlYKCn4Y2eAr41oLQrDGn8fNZVcRfXQ3qdvRLoeHWclxHAcw5UpeO8HzlO9Zj5zLssfAKCWGV7JZf6/QAAAFZlWElmTU0AKgAAAAgAAYdpAAQAAAABAAAAGgAAAAAAA5KGAAcAAAASAAAARKACAAQAAAABAAADiaADAAQAAAABAAAAfQAAAABBU0NJSQAAAFNjcmVlbnNob3SXMOzLAAAB1mlUWHRYTUw6Y29tLmFkb2JlLnhtcAAAAAAAPHg6eG1wbWV0YSB4bWxuczp4PSJhZG9iZTpuczptZXRhLyIgeDp4bXB0az0iWE1QIENvcmUgNi4wLjAiPgogICA8cmRmOlJERiB4bWxuczpyZGY9Imh0dHA6Ly93d3cudzMub3JnLzE5OTkvMDIvMjItcmRmLXN5bnRheC1ucyMiPgogICAgICA8cmRmOkRlc2NyaXB0aW9uIHJkZjphYm91dD0iIgogICAgICAgICAgICB4bWxuczpleGlmPSJodHRwOi8vbnMuYWRvYmUuY29tL2V4aWYvMS4wLyI+CiAgICAgICAgIDxleGlmOlBpeGVsWURpbWVuc2lvbj4xMjU8L2V4aWY6UGl4ZWxZRGltZW5zaW9uPgogICAgICAgICA8ZXhpZjpQaXhlbFhEaW1lbnNpb24+OTA1PC9leGlmOlBpeGVsWERpbWVuc2lvbj4KICAgICAgICAgPGV4aWY6VXNlckNvbW1lbnQ+U2NyZWVuc2hvdDwvZXhpZjpVc2VyQ29tbWVudD4KICAgICAgPC9yZGY6RGVzY3JpcHRpb24+CiAgIDwvcmRmOlJERj4KPC94OnhtcG1ldGE+Cry305sAAEAASURBVHgB7V0L/BZT+j9IouLPYlHEZq1CSKiEcsld282lhETJPZbFLrKEyFpLkSQtKRRyL6UsFX+s22LpZim5bLLK7v632vk/37POdGbeed93Zt6ZeW/f5/P5/WbmzLnNd+adM885z/N91nNEFIUIEAEiQASIABEgAkSACBABIkAEiIAgsD5RIAJEgAgQASJABIgAESACRIAIEAEiYBCgkmiQ4JYIEAEiQASIABEgAkSACBABIkAEuJLIZ4AIEAEiQASIABEgAkSACBABIkAE1iHAlcR1WHCPCBABIkAEiAARIAJEgAgQASJQ9whQSaz7R4AAEAEiQASIABEgAkSACBABIkAE1iFAJXEdFtwjAkSACBABIkAEiAARIAJEgAjUPQJUEuv+ESAARIAIEAEiQASIABEgAkSACBCBdQhQSVyHBfeIABEgAkSACBABIkAEiAARIAJ1jwCVxLp/BAgAESACRIAIEAEiQASIABEgAkRgHQJUEtdhwT0iQASIABEgAkSACBABIkAEiEDdI0Alse4fAQJABIgAESACRIAIEAEiQASIABFYhwCVxHVYcI8IEAEiQASIABEgAkSACBABIlD3CFBJrPtHgAAQASJABIgAESACRIAIEAEiQATWIdBg3S73iAARIAJEgAgQASJABIDA6NGjVb9+/dQmm2ySKCBjx45VQ4cOVStWrEi03lqv7Pnnn1edOnUq22UuW7ZMvf766+r4448vWx/YMBHIEoH1HJEsG2RbRIAIEAEiQASIABGoZATGjBmjbrnlFvX++++rDTfcMNGu7rLLLlrxPProoxOtt9YrO/vss9UOO+xQtsv829/+pvbee281depU1bZt27L1gw0TgawQ4EpiVkizHSJABIgAESACRKDiERg/frz6xS9+oWbMmJG4gjhr9mw1f/589eijj6pevXpVPBbs4DoEttxyS3X11Verbt26qXnz5qnmzZuvO8k9IlCDCHAlsQZvKi+JCBABIkAEiEA9IrB48WJ1zTXXqIYNG6p77703MgQffPCBateunbruuuvUJZdcErl8sQInnXSSevnll9VfP/1UNdhgg2LZeb4CETjmmGPUkiVL1Jy5c1STxk0qsIfsEhFIBgES1ySDI2shAkSACBABIkAEyoTAp6J0nXHGGWq33XZTDzzwgPruu+9i9eSss87SdaShIKJP8Kvr3bs3FcRYd6cyCt13333qs88+U5ddellldIi9IAIpIUAlMSVgWS0RIAJEgAgQASKQDQJTn3xSde/eXe23336xG7z77rvV3Llz1ZAhQ2LXUajguHHjtPIK37ooAuqIm2++WXXo0EFts8022h/u3HPPVX//+9+jVOPm/ctf/qJ69uyp4Bu5/fbbq2OPPVZNmTLFPZ/lTiX1Jex1//jHP1Z9+vRR999/v1q0aFHYYsxHBKoOAZqbVt0tY4eJABEgAkSACBCBIAQuvfRSNWLECL1a98gjjwRlCUzDKl/Lli1V06ZN1YIFC9T66yc/hw6yk0aNGmlFNLATeRKhxD3zzDNqvfXWUzbXYJs2bbTf5FZbbZWnZG7yK6+8omAuieu162vQoIG67bbb1HnnnZdbKKWUSupL1Ev861//qpXs444/Tk1+dHLU4sxPBKoCgeTfglVx2ewkESACRIAIEAEiUGsIQPGJIzAhBHslVvnSUBDfeusthb/+/ftH6h5Wq+bMmaMVOIRg+Oabb9Tll1+u63j33XfVXXfdFbo+KJgDBgxQ7du3VzNfnKn+7//+T/3xj39Uu+66q1qzZo0mZVmzdm3o+krJWEl9iXMdLVq00KEwHn/scfXGG2/EqYJliEDFI0AlseJvETtIBIgAESACRIAIpIkAlMSNN95YDR48OJVmYMq6+eabq1NOOSVS/X/4wx/U8OHD1UUXXaRg5og6brzxRs2wiYpmzZoVuj4om1DOnhTT3EO6HKKZWw888EDNtArlGnEb54m5bRZSSX2Je71gwP3Pf/7jKu1x62E5IlCpCFBJrNQ7w34RASJABIgAESACqSMAP8T33ntP7b777trcNOkGsWL32GOPacUOimgUgQ/imWeemVOkS5cuOg1KSlhZunSp+s1vfqM22mgjTxFc97bbbqvT1ma0klhJffGAEeFg//3314o7VmOxCk0hArWGAOMk1tod5fUQASJABIgAEahwBJYvX65XtRB7rtyCVT7IAQcckEpXJk2apJWIQYMGRa7/oYceCizzxRdf6PTOnTsHng9KPPHEE4OStdkpSHA222yz1DDwN1xJffH3LcoxFEWszE6dOlWb8kYpy7xEoNIRoJJY6XeI/SMCNY5AqXHNahweXh4RqCkEvvrqK70yBtISmDeCkOXoo49WUKQ+/PDD0Ne65557qh49eoTOXygj+gKJonAh/9dff63JX4opujBlxWodfAGTEphrQoBdqTJXAsN///33qmvXrtoEtdT6SilfSX0Jcx0HHXSQVhKfeuopKolhAGOeqkKASmJV3S52lgjUDgKIazZ06FD9cfjPf/5TsxHWztUFX8ns2bNVu33bMQBzMDyppX777bcKygko/ynlQwAEH1DsoOB9+tmn6oLzL9DEKVB07rzzTk3QErZ3J5xwQiJKIshgMFG14YYbKmPCGaYPL7zwgurXr586+OCD1cMPP5y3CEIkQAm96aab8uaJegLKKUhrfv7znyusZJUqTzz+uGrYsGGifYzbp0rqS5hrgGINefnll9Xq1avLrmSH6XOaeTjGpYlu/rrTGuPok5gfc54hAkQgRQSSiGuWYvcSr3r06NF6pn7535YnXjcrLIwAPujxMY+YbJTyIPDnP/9ZHXXUUap58+baP69J4yYKacbEc+TIkZqEBUQsYf5A3pKEGOKXVq1aqU033TRUlVD4EEbiyy+/VNOmTVMrV67MWw6mrPBDPOOMM/LmiXriggsuUFi9HDNmTMGiAwcOLHgeJ7F6e++996phw4apvffeu2j+NDNUUl/CYAcs9thjD+2XCNbZF198MU14Kr5ujnHlu0WpjXHCdEUhAkSACJQNAWGIc+TV6vTu3btsfUi7YSGLcIQ90GnSpIkjqwBpN8f6fQg8/fTT+hmTwOGOrOz4zvIwbQQkppwD7PEnipVu7qqrrnIk8L2zes2aRJuXOImR3ifnnHOOzi+xCEP145JLLnF22GEHp127droc3l133HFHYFkhgdHX3L1798DzcRLRlrCcOqJQFSwuiqwjymnBPLL64Ej4C0fCfhTMl8XJSupLGOxsTGSiQz8LV199tZ1cV/sc48p7u9Ma47iSWD7Fny0TASIgCMSNa1Yt4F1//fU6/hj6O2rUKL0CUC19r5V+YtXn8MMPV5999pk69NBD1ZIlS2rl0ir+OmCCBxNTsFnCN08mSnTA9okTJ6rJkyerBhtskOg1IN5fFHnnnXd09i222CJUMYSiWCgmpPAJbNmypS6DMBVBAj81PHNhVqVglgqTxUIsmSBIue666xTqRWzDfPLggw/qdw6YT8Gsij/5hPVkR9pxxx2n2rRpo99LnpMJH2AF/9///rcOFwEmWZDk2JJlX0y7+fAOg52pw2ybNWumdw2ZkEmvly3HuPLf6dTGuPLqvmydCBCBekcg6sx/NeElHxyOBObWs8y//vWvq6nrNddXrOBiJUuGc0dIRJx//OMfNXeNlXJBYnrnfPLJJ7o7559/vsZcFBJHPiYdCULuYB8rR2nIPvvso9sTv8dQ1WMlDc/EkCFDQuW3M91yyy26LMp/8MEH9im9Lx9uzo477uiIspZzziQ8/vjjjiibekWvbdu2ur4jjzwyx+Jg5oszna233trB1i8zZs5wkyQcgyOmqG6/sI8/CfHh5sHqrfiB6j//Sq74ijsfffSRmzfujvgAOyNGjHD23Xdf3ZcnnnjCEWVK7yPNSBZ9MW1hWwjvMNjZdZl9Yw3TrVs3k1Q3W45xlXOr0xjjMLtEIQJEgAgUReDjjz8u+LFTtII8GWpVScQHR+PGjfVHkaxi5bl6JmeJwEsvveQ0aNBA35OTTz45y6brqi1hLNWKDxQ1mFlDiYLpI8wu582blwoW48ePd/baay/dFtrDH8wAb7vttoLtSXxAnRcKbFSRVT9nk0020eUvvPBCT3FZVXIkHqFz5ZVXetLtAygsmES6+eab3WTxC9T1Qbk28uqrr2oTUwmzYJL0FgqW+HI6nTp18qTDvNdg7jkhB1BYhfTHEVZO51//+pfn9Oeff+6Ij50j/k2e9DgH+K1JmAvdD/QFOEu8R+dHP/qRA2UYklVfTP/D4F0IO1OPfzt8+HB9nR06dPCfquljjnGVd3uTHuOoJFbePWaPiEBFITB27Fg9Gy7sf07Hjh0T71stKomY0dtuu+30h8Pmm2/uYHaeUhkIwP/KKBH4wKakgwBWauE7B6xvv/12B/55lShGyROCmVjdw+oRrhG/d/saoXRiQgJKR5AAHyhMEhbDc1pMcnV9vXr10ulYAcQKItoAnvafEO3odJSxpZCiA0UNdf3P//yPpy7UC4V+5513tqsqaV9Yq11LCqwqQr777jsH6ZAs+xIW70LY6U4H/Pu9+IkC05/+9KcBZ2sziWNc5d7XJMc4hsCQXzaFCBCBXATkFaiE1EHTuz/wwANKiBqUmIppvxKZ/S5rXLPc3lZWinz8KJmV152Cv4aYOVZWB+u4N7Jqo2PzwVfssssuU2L6pv/qGJJULh3so2D/lFUuBTbOSpRV369SojzoronCFquLAwYM0IHU8XtHQHVZLdX1wE8RjLpCchNYL+JDLl++XA0ePNhzvn///uqwww7TLLA4IaanOnwL9oGnX8By2qdPH39y4PEdEmYETKYQUOYHyamnnpqTHDeGZaNGjXRICPgcitKr623atKneZt2XsHjnXHyIhMabbKJz4X7Wi3CMq9w7negYV7m6MHtGBIhAGgicddZZoaoVynY9Ew5fEgjYs8DoZ8QwusmrUs+iFtvCxClIam0lESuvBgv/KkHQ9TMtewSmTJni3iOJnejIR2z2najhFiUwu179wkran/70p4q9UpkocJ+D6dOnx+qnYTDFb94wpMLkC8cTJkzIW+e5556r89x6661588Q9EWc1rFBbpbzrYXILLPKtqBZqN+hc3L6ExTsOdjIhoK8RK7GF/E+Drqca0zjGVf5dS2qM40qivL0oRKBeEBDfCQX2tnvuuafgJYsfjWYivPbaa5WYU6mlny/Vq4piNuaWQ1yzFStWuMfFdvLNqBcrV03nMdMvSq/b5euHXe/uc6dyEADb5oEHHqjZJMXXVmG1VyjcK6eDVd4TvBuEuET17dvXE3vv7bffVuI3WDFXt9VWW7l98TNuuieK7MCqQibAlCh7aubMmUrM8BRiI4r5ppKwPnlLg+0TApbNSpdKetfH7UuaeJtxEHE2a52tm2Ncpf9a/9u/xMa4yteH2UMiQASiILBw4UIHDuXwGbBFTEadhg0bajIFEBbgL2jWU2jptW+KfETr4mDTg6/F0KFD7eoS26+llcTTTz9dzyjLa9oRU7PEMGJFySOA3wjuE/4Qv3LBggXJN1KnNWIFHauIEvrARQCx58IyjrqFMtiR0Bf6GZDwNLFbAxuoIegBi/Fmm23mYOWqkEgoC91u8+bNcwhkUO5Xv/pVYHqhOs05sxqGVbxyS9IriXGvJyzecbBDfES8R8CeXOvCMa567nASYxzjJMovm0IEagEBMQtVQjqgEC8HsbwwS37UUUfpuFuIvyUU7zpWFfxD5MNE/73//vv60oUcQfsaygCpxOlZwW8Es+M9e/VUfU7uo2644QZ1zTXXpAJT1LhmqXQigUoRb03My9yasDpFqVwEsJKI2ImQVatW5fiGVW7PK7tnMvGkJNC79kX82c9+pjs7ZswYNXr0aPXcc89VXOeNL2Kh+ITFOi0my0oUY51t2LBhOg6g39fQX4eYpuokxOw87bTT9LvZ5BGmVCVmukoULJMUaStMsjo/3vVBfoyRKksos3xaJ1RTvGrC4h0HO/PsCBlQvM5VSSmOcVVyo37oZiJjXPXoxOwpESAC+RAoldobDH077bSTg1lteb/oP8QPA7X3ylUr8zWbSHrUuGaJNJpCJYcccoiLnZ+SPoXmWGUCCPzv//6ve8/w3MOPg1IaAnhfgAlZyGp0RVhBxLslidh7pfUsuDTCFuDe+0NYBOfOnypKsPss2XEA85dwnB49e7hlEENQCGs0SyqYRsUMtVDRoucQnxHXJcqqDrGBWH5Zy+LFi93rmzhxYtbN57QXFu+o2JlQH8b6JqfhGkngGFd9N7LUMY4hMKrvnrPHRMCDQFLU3qBsx0cFlEOZffa0kcZB3LhmafSl1DoRxwzYmb/777+/1CpZPiMEYBpp7hti1FFKR0AsEByEfoGZeicxW//kk09KrzSlGkwIi1NOOaWkFqAcI6wEnqWwpqsIBSGriFqpNs8gFMx33323pL6gsPiTuzEcobykPdnn7zDe7+KH7v62oPgKc6o/W6bHYfGOil3Xrl31dYofWKbXk2VjHOOyRDvZtkoZ49ZDV+TlRCECRKDMCMSlGRcfQk2QIL4wSvwuPFcBmn+YlsKZXmL16RAWMKcxtO8ms8SuUi1bttSmqRJIWVOwm3PcFkcAtPEIEwKBmS9wj2sqVrw15kgSgXHjxilh8nWrBFX+0Ucf7R5zJzoC+Kx48803lfh6Kpl0il5BhiUGDRqkibxgmv/ss8+W1LIomv8Nr7LkM9WkcZPQdYE0R/w3lawmuqEvQhcukBGEY/9Z+x+G4PFhFAbvKNjtt99+6vXXX9cEWOJL6mutNg45xlXvfSxpjEtWX2VtRIAIxEWgnNTewuyoZ0KPP/74uN2v23LywaGJKmQI0RgikC2lehAAgZMo9vre4R7id0ipHwTkA0rfe/GfLPmiYSIaN5RGyY2zgrIhsM022+hnCEQhtSgc46r7rpYyxjEERvVODrDnNYZAuai9QW4AynZhPlW//e1vXVSx8ti4cWNliB3cE9zxIHDfffdpogqTeM4555hdbqsAAaz4gqQJvz/InDlzlHzsqYMOOqgKes8ulorAz3/+c4XVRIRCAckLQlfEFfG9VPij1A8CIGn64osv9Fi5//771+SFc4yr7ttayhhHdtPqvvfsfRUgIHNQ6s9//rNmvgRjHdj/gkTo4VXnzp1D//3kJz/R1Zj4g1OnTlVQ+PwCM1Q73d/+7Nmz1eeff67wsQSTUwhinB166KFq1qxZ/up47EMAphxGwC67xx57mENuqwSBAQMGeHp6yy23eI558F8Evv32W22SCSZlw+hY7diAkVL8AMHPoGQVsNovh/3PGIEZEhsTsttuu+mJ1oybz6Q5jnGZwJxqI3HHOCqJqd4WVl7PCCB478CBA/VKHBQH+KsIk6eSuFw6RAU+SpKQUqm9582bp7tx3nnn6S0CA4sjvu5vr169kuhizdYxd+5cJSQT7vVBsaYkg4CQnShh7VUSn1P16t1LXXzxxbpi/K4Q+B703kn5/whhjfbXNT1/8cUXFRQiyn8RQAid1ru11u8uhNjp3r272nrrrXUaaPGrXYwPKifFqv1OZt//P770km5UmD+zbzyDFjnGpQdyVYxx1W1py94TgcpEQFbmHDAlggpeYhY6YsLmPP30085ee+2lfRfktePIyl1gMPs4V1QKtTeo6tEfML8tXLjQadOmjVMOuvQ4113uMn379nXvJzDEPaaUjoCs5Dnbbruts/7667v4Xn755Y4obg7YH4G1+VuwYEHpDUoNMtPq1om677jjjkTqrfZKbr/9dkfMlTTr8aOPPurIpJKD4OFIA04IRi/KVVVfppgM6mtJwi+xqoFg5yMjAAbXBg0aVDSDb+SLsgpwjLPASHC3WsY4hsBI8KazKiJgEDBxk8TU0yTp7bJly5xNN93U/RgVRkzP+bgHpVB7g9oadPUIfYEPc3wUUoojAIVls802c+8l9sWst3hB5giNwC9/+UsXX2HddUA1j9+PsGbqdGHtdYRJNnR9hTI++eSTbltQfkhg4zhQnsRX2UH4AjEv9cAnvrcuXq1atUpswsvTSIYHmMDD85TUpEOGXWdTZUIAsT/xrsB7qRaFY1z6d7XSxzgqiek/A2yhzhB466239MwiPmShvPnFxFTC4ILgtEkKXupQ+gp9OC9ZusQRUhpPs+KD6IhJWeaxtDydqLIDIftxP5JxL4844ogqu4LK7y4CnANbKIaXXnqpnsz44IMPnGnTpjkdO3Z0brzxxsQuAgq+iXOHNrGK+de//jWx+quxIrAdAwus4vrFH6T5lVde8WepquPnnntOX+tZZ51VVf1mZ8uHgJkomTJlSvk6kWLLHONSBPeHqit9jCO7qYyAFCKQJAJjx45Va9asUatWrVJdunTJqVpmH920+fPnu/tJ7MhqlirGsNZsu2Y5TSG2H/4o4REAA6YtorTYh9xPAAH4wkHAsHvnnXeq2267Tcmqlf6D32ySAnZf+AzP/IGIAgRPiH2ZlN9jkn3Noq5V369STz31lG5qwoQJCgRXtuAdZwt8E2X11U6qqv0jjzxSv68feeQRJZMPZHWuqruXfWcRW3jixIkKZGUgfatF4RiX/l2t9DGOSmJKz4CsJmla7auuukodd9xxKbWSfbWyWqID8+JjrWnTptl3oApaBJU6BB+2+KD1i50GEhtKtgiAqAdss1DQxZdEB6/u0KGDAnlJFHnjjTc82RFQmZIcAsuXL3dJgRYvXqzZdhGqIE2xlUS0A8WoXpXEDz/4UDN+AgcxRQ8M8m4z+bZu3RpZq1pGjBihmU7F5F6TI1X1xbDzqSKASSuQvD300ENKrA5SbatclXOMSxf5ahjjIimJ+LAS8xO9SpIkdPhoFvOhJKssW12Yfb7uuuvU8OHDlfg31Bwl8j/+8Q+FlbK58+aqac9P0wpj2cCu0IbxQQvBhxXiC1HKj8CatWvV2HvvVfj4Q1yrIIGSh5lhE1okKI9Jw8vdXgXGb719+/bmNLcJIPDss896xhpMTKUtfkX/zTff1CFravUjsBCeYtbrnr7sssvUYYcd5h7X6k7btm0VGJ3HjBmjLrnkEgXLDAoR8CMAhuXRo0er3r17K6xAJyX33HOPuv7660NXB6shIZPy5AdLufhXe9LyHYCR9f777w88zTEuEJZEE6thjIukJIpfiBJ7ffXee+8pIRFILE5SMfO4RO9KipXBPKfb8d0U6NMRe2nSpEkKK2+1JDDF6tO3j5oyeYo2LRI/Dh0fqJausdRr+f7773UV4s9UalUsnwACUBB7CGW/MZ1DsGy8y3B/MNgbER8rNXjw4FATVi8J7bm4FJiiascdd9S/eTeBOyUj8Pzzz7t1dOvWLZP4k36TYawUvP7660VNuN2O1tCO+De7V/OXv/ylLpREXLD4Yen7ffbZZ+tJIxcE7hCBHxAYMmSIwoQ5VhOTlN13312dccYZ6u2331YzZsxQ5lvC3wYmoI866ih18MEH+0+pnjLJgfENCy+IyxwksHJC6KyePXsGndZpHOPyQpPYiaoY4+K6ZRpWJ0FLO3ubbb9+/Rx5wB2EAPjyyy/1dtGiRY48rJoe/uabb3Y6HXigp8xNN90UtxsVU+7rr792xFRJXxdIFsSWu2L6lkZHJOafvlbxY9OEJ2m0Ua11tmvXzn2+QQZTSOTjywERByU9BEx4EBmAHTBkGkG4D4n15t4rvMNkwsqcLriVjwRPORB8UJJFoHnz5i7GCLuQldjt4pkYNmxYVk1XVDs222sYgi3kX7t2bUVdQ9zOiMuAZnRFyA8KEbARAEkNmKxnzJxhJye+jxA85rva3voZ0ws1jN+tXRb7LVq0cPC9Wkw4xhVDqPTz9lhTqWNcbHZTmWHNefhAlQ12xTAifh5u+RerPMbS3//+d6dV61b6ehAXb/LkyWEgqOo8q9esccT8SF8zlOKZL86s6utJsvNGKcEL+fDDD89btTi+Oy1btkyUoTFvY3V6AgM57gNwxjvLL2Iy5L6HkE8CtvuzBB6DWRP5zd8111wTmI+J8RB49913XWzxbs1SxHzMbRv3t15ZazGBZZ5vbCWodt7bgLA5CJNRSyFgZCXFQQy8WrqmvDeQJ0IhgDAw+LAHE27aMn78eM/vD7/BsJOYpm9iDptTx6hRo8zpgluOcQXhKflktYxxsb1t/Q6t8gBrZriwNvzwcYCfB4gj9q9iwgd5UpTExFNw8ofAHrzQEr7OVCX/YE4B0+IgabDBBkqUYQXTPbB89T+9v94G5a23tMMPO9y9ZJhlS9BU99jswDQZ5h4rV65U551/nknmNmEEnn7qaV0jWGZhAu6XM888U98fmPiA0fCGG27wZ8k5hvkq/LNtESXUPkx0v9DvMNGGKqgyYxqMLp3Q+4RMewa2QltkMLcPM9uXlW5tIg0sbP/ArDrw05/+VO20005uc/itSOgc99jsYBy44oortKk2GGJrRQ466CB16623ap/UWrkmXkdpCEhIKyVKVqJ+iPl6BH9ov2CciiL+73T8Pk866aSiVWQ9xqFD9TbOVc0YF1cdFtbOnBmKqLPpwuyoA3jH7UMllLvyyitdHJo1a1Yzceb06qgESMasWSERNjj3+vv27Vsoa92cEz83R3zUXFzk/ecIRbYjhD8OVs2FgMMRchR9HvuU9BAwq7rCwphYsG+YheCe2n9pmZeH/R2mh2B5au7cubOLL1wbshS4RNj3FvtY9c9a8Myafoh/XNbN6/ZGjhzp9gF92W677ZyhQ4c64kvjPPHEE87JJ5/swHpGlMmyYFQWUNgoEcgAAeO+ZN4BQo7mLFu2LHTLMJk2Zc027EpklmMcLqgex7lqGeNim5sKZXzOAyizFqEfYGSUVUTnhBNOiFSmkjLLaoIDE1vzAxw3blwldS92X4Sh1TEmV4ceemjBeqAQycy7i8Fjjz1WMH+9nBTWKk9gbvOM2FsE4qWki8BFF13kPpvYT0KClAhZYUmiak8dUX6HnoJVfiCz9U6TJk30fROChsyvRuLkuc+M+b0icHzW8rOf/cztR7mURLzfYW5rcAjawswUpqkUIkAEkkHgm2++0ZMv9u8NPvVRxJ7AN/VgUSOMZDXGoS/1OM5V0xgXy9wUJnJ+1iSZYdTmpvIwhhZQ94oTbuj8lZZR/JdcdkTQ359++umV1sXI/QHbY48ePZRhXSpGey6zyB7KZpjbghGw3gXMY2C53XfffRUwsgUmXKCdlll6O5n7KSAA028TvuB3v/udEkWx5Fb8JncbbbRRYAy5UhqK+jsspa1KKwuW2VWrVulu4XeUtdgmlqZtsHvWo+DdBQZrhIOAa4EtcBWRlUTN/ioKrX2K+0SACJSAAEIjrF692lODrDx5josdTJ8+PSfL0UcfnZMWlJDFGId263Wcq6oxLsysgj9P0Eyr2Dn7s3mOQRCB5e9aEbAkyjPu/kl8m6q/NKyGHHDAAe414frEJzHUdYExy+BRr2yA+YACYcrUqVMdrC5K2IV82ZieEgLyges+m3hGYb0g/g+xW4OJnXnWsQW5RZJSyu8wyX7Ua10gp7DvL/ajMAomhVslrCT6rwXWQmD8lLAgjtDz+0/zmAgQgQQQQJQA/zvomWeeCV0zyJbAwGrXseWWW4ZmH057jMOFcJwLfTsTzxhljIu1koiZRb8UmqFAYM9BgwaFij/mr7dSj/1B0uPGQ0TgdTHRVOLPqQOzIiaXEcwkPfzww+rCCy/UhDj+2R2Tr9QtVv5A2LH33nurOXPmuNXJS0UHE/7ss88U/sQswD3n39mn3T5u0oQJE9x97ihNmCIhEnRcI1EoCEnGCIjZjWcFUSa5FIKmz58/P1ZPxBTIU27bbbf1HMc9SOJ3GLdtlluHAGKI+QnYFixYsC5DHe+Jn5QONi9hftQmm2xSx0jw0olAegiIj7uncjG/Vwh8H1ZmzZqlxM/Pkx0xYI1VjedEwEFaYxya4jgXAHjGSZHGuDgqKpzU5ZrcPziuw4Y6nyBuHJxuMXNQC4LwD5iVMRiAhCSsrFy10sEskZghOptvvrlbh6nLxMyDgzKcjE06tqDxT1pAQCTmcp527DbtfcS+zCd++/dXX301X1amE4GyICDspY4Mku6zjhiJcQhn/L9LmSAr+XqS+h2W3BFWoBHwE08htm/WUokriVljwPaIQL0h8OGHH7pjlPn+AslJFDn//PNz6kDcxbCSxhiHtjnOhb0D6ecLO8Y1kIcwkgjTnMLqly177bWXEoXHTnL3L730UgUaXnG6Vdtvv72bXs07L4ittyzXupcgJD7ufrGd1f9erfEDHb/fVw1+na1atVJLlixRmPXB6p0toEQH9kE+M3a+KPvwLzX+JJgtBw2xEVDBm9li+F0VWjHx28uPGTNGyYvGVBV5i37Bpy9tQegCPJuU2kYAM7N4d9mr4V999ZWCBYSwzioxQQ0NgMSC9eTdeJONPcdxDpL6HcZpm2VyEWjUqJEn8fsffCQ9iTwgAkSACCSMgB0awVRdjBvC5DPbmTNnml13e9xxx7n7xXbSGOPQJse5Yshndz7sGBdZSQx6gKG0vPbaa/rqYCKJ5eSlS5eqSZMmKQlIq9O7du2a3dWn3BLM1WyBQhdWoEy//PLLOrvYmOtYeaZs+w7tXUdeKIhQ0GylDccwAU1SJASDW90uu+zimuDBLEFokEO3B1NVKL7m5SI+eEr8NEObN7id+GEHCrGEjfAnJ34MMoYvvvgi8XpZYWUgANN48ZH1mFHbPQNByimnnKK23GpLdUiXcOY8fjOejRuVriQm9Tu0r4378RHwD6D2ezh+rSxJBIgAESiMwIwZM3IyYKHl8ssvz0kPSlizZo2S1UjPKcQ8Fd4IT1qhgzTGOLTHca4Q6tmeCzvGRVYSERzcL1Ca/IqTP0+3bt38SVV77A+kLfERY11L06ZNPeW6dO6iBg4cqNnievbqqR74wwMKyvUrr7yi8yH4ur+Mp4ISDuDvaPtoYXUtikIKpRL5jZKIlVbUKUvaJfQq/aJ+BrH0W2QLWSDwr3/9S+GdYzO8YYIGLI1YNRRzaIXg4GBXwzNw2qmnqbfeeivUM49JMFvMarudFne/1N9h3HZZzotA2AHUW4pHRIAIEIH4CGAsAvOlXyQmqT8p0nGXLl0i5U9zjENHOM5Fuh2pZA47xkVSEvFBZVYMTa9hMmmoscWKVq98QVHAvpGtttpKderUyRxW/fbrr7/2XAPMROOIUf5MWTgLjx8/Xu25555qwoMTFEw8QYU8ZcoU1UKUrS6dO5usiW9NyAtT8cEHH2x2Q2/9H8sLFy2MrSQ2btw4dLvMSARsBKAgInTC7Nmz3WTxG1aY4MIWgplVWDnA1BQCE2+E47n77rv1cb5/ElRdCXOc5/TGG5e+kmgqTOJ3aOriNj4C/gFU4lrFr4wliQARIAIhEMCYZCbaTXZ8P8NKK6zYk/2mjMS9NrtFt2mPcegAx7mityH1DGHHuEhKYhBjkgRb13GU7CvCbAhMKREzD2anBx18UGyzQ7veStm3/RHRJz8TXth+Ak8j22yzjf5gxQcn2EGhIEKwcphF/EX/CnEctla/Yrdo4SJ16CGHmkuMtMVH/JKlS9R/1uZnVI1UYZ7M62+wfp4z8ZLxgoXvGyUdBKDkgZkrn0BBxLNrs8M1b95cK4TY2oKJK6MkIt2Ygdt5/Pv+CSKc90+O+MtEOU7idxilvUrL+84777ixZ8P2Dc+DUf7DlimWz6/416O5KSZOhECtGFQ8TwSIQAEEsJAC3o4wgkUBv9x6661KyA79yYHH8LnHoo39jYp3WZTvubTHOHS8nsc5mAKbOMCBNzEgEXrArrvuGnAmflLYMS6Skhj0AAfNUOBHAX8yhG+AX+Kxxxwb/0oqrOTy5cuVxIfy9MoodJ7EIgdYlbXDXRi/uMsuu0zttttuRUonf3ru3LlupfjojUK3bAr6lcRPPvnEnIq1bbZdPDPeWI0lVAgzgQh3QUkHARAiwUw0nyCUi60gNmzYUE+6+BVElAeNvy0ImI7fdyElNGgA9b9s7Tqj7ifxO/S3CZ9M2xfEfz6J48mTJys/eVWcekEkJCzKkYr26NlDTZk8JVKZYpn99xRK4pq1a1WDDTYoVjTvebgpRMHI9gvCCrfEJ8xbt30C/uH+jzD7fNj93/3udwofqBQiQATiI4AFgLCTLX6yPrjxwComrGD8sBVElGvbtm2kicy0xzj0KelxrprGuFNPPVWTeQKHsCLREALNkMOWD8oXdoyLpCT6H2A0XIhcBNoyFMYorEpBF1NJaUaZs/u0QYwPh6BVWfg2Dh061K46k318vNgfZnip+B+gMB3xK4kSOD5MMeYhAokggMFx5MiRnrr69u2rDjroIE+aOfA/45iFBeNpISXRNqM39QSlmXNRtkn9Dv1tYnUVym+aArKEWhZ8rK2/3nolXSIwKuU+hC0rYZRK6icLEwEikD0CUM7ee+89T8N77LFHKD95UyhoIScqM2rQeBaUZtqMuk1jnOMYF/Uu5ObPN8aFVhIxE/L+++97asbyZyHGpI8//ljPYhT66PJUWAUHfkUIXfavLIa5jKAfM3yi/B+uYeoqNQ8YIG2J6uRsymLVxpa1MvNebyKx91Rc/OoNqzjXG7QiaOq59tprlR0EGCvihSZdMFj5pRhZUxBxVFL+akn9Dv3XVE3HHTp2UN8s/yZSl/dss2ek/GEy+81LQXqEQbSeBKGR+C6rpzvOa00DgbD+hPgmtEM0oS9Rf39BzKjHHHNMpMtKc4xDR+p9nMOqYBDGhW4SQuMlLWHHuNBKIkJf+GcT4I9YSGAahpXEWpIddthB+wva5BV+R+Mw1+tflUUYETCblkP8fQkyIQ7TL79SCDOLehPEhvTjWW8YlON68W6CL68tWEXE7zWf+Ce9oAQUm9AKern7X7b52iuW7n9u4v4O/e3AHyWOSby/nkLHCN+ThEx+dHIS1ZRcB2ambdl0003tw1j7sBS57rrrQpe988471Zdffqnzt2/fXoX92At6RkM3amU866yzFP4oRIAIpI+An8wFLUYxNYV5Ony6bYF/ot+twj4ftB/0/khqjEN7aYxz1TTGjRo1Kgj2zNPCjnGhlcRp06blXESxQSufmVdORVWUgA9JBJW3/e1AVhJFYLLq/0AdPHhwWWaqMXNlUy7jIxkfJHEEfpa2lKIkGrZJu7409rHaVCk/2jSur17qXLBggY7Paq53PTENvOqqq8xh4HbOnDmedFhGFFstClIWklhJTPJ36LkoOQBBTy2xS/uvL41j/z0NuxpQqC9gKYS1SFh58MEHXSURxBdRyoZtg/mIABGoDAT8xGkgROwcgdEeSqb/G+yAAw5QGAujSFpjHPqQ1jjHMS7KHf5v3rBjXCglEbP0/nANeIDjkJugexdccIGCqZd/RiH6ZZanBGaEbSXRJhgI0yP/qiw+QAafMzhM0cTzzJs3z2Oit99++xX9UM7XCT8OcUODoH74lyEcSNqyxRZbUElMG+QM6v/Tn/7kaWXHHXdU22+/vSfNPoC/rJ/cAwNqMWkasKKUxCxrkr/DYtfA88UR8M+yJqEkFm+VOYgAEahHBOCLiEgAtuyzzz7K78Jjn/fv+804cT4Kq6mpL60xDvVznDMol38bdowL5WSBwNMgdLAFD3AcEyawnd5xxx1qdRUTHfg/Pv3KkY1T0P7z0573JHft2lU1adzEk5bVgd/EoZgJcaF+2f5gyNe8wEd6oXp4jghERcDv8I9Yo4UEzI0I1WPLiSedaB8G7oPd0h/y4p//Kj2GXpK/w8COMzESAv4BdPMtNo9UnpmJABEgAmERePrpp3OyRv0Ws1m9URlWEI85Npo/IsqlNcahbo5zQKEyJOwYF2olETEP/RLVoRblX3vtNXX22WfrqvqffrreVuM/vxNpFDY5rMrOecVr5nbaaaeVDQY74Dg6AQp6CFZKQX8+fPjwnI9inSHgn02djJXmjh06BOQKl4RVSFC5+x25w5UOnyvI/j58aeasFAT8JjWFwsjAzOKhhx7ydB0zrmFjeoK8yl49/Oc/SlcSk/wdei6MB7EQ8A+gP9oif2zOWA2wEBEogADeLzA3hrWV/3ujQDGeShCBSy+9NPS3T6nN+q1aUF8xdy67TZBELl682E5Su7baVcUNI5bGGIfOcZzz3KKyHoQe40RpKSpigujI1Xj+JM5J0XJ2hokTJzpisqPrEBM/Rz7U7NNVtS+mak6DBg1cPMSHL3T/ZVXWLQdMxW/PEcKX0OVNRiHLcS6//HJHgs47srLp9OjRwxETTXM61FYeEkdeBm5/UA8E90ZWYhwh/XCEuTVUXRIg1K0H1yVhT0KVYyYikAQC9913n+f5O++88/JWK/6/nrxi0uOIj3De/P4TQkvuKS/EAP4skY6T/B1GapiZ8yIgfjmee4x3bdYi7KJuH2RyNevm2V4ZERg9erS+92Kt5T4D/m8wHnu/SZPGQ0ijMnkChNPC8x2G65BJ8khty2R+znNyzjnnRKrDzpz0GIe6Oc7ZCJd/P+wYV3QlcennSwMDP4ZlTJo1e7Yac889CmamAos8/0oHG2/UqJHeD/q3aNEiNWLECCWKh0JsKawKdOzYUT3xxBMKBC+gqYcP0rvvvqsQABtEMvfff7+nqvnz52vfx8ZNGqu777rbjTWDWFMIYv/mm2+qt95+S+36s13V9ddfr8vCXvree+9VWBmEOe1Nw28KNAMFYyIcig3dMEzdwHYaxvzWvyqLFYxiZBmeC/vh4Pe//7266aab3FOfffaZjrEDMg4QJIQROErb4Tvgk4V7dOKJJ+q6pr8wPfRM2gszZniaLBQ/05ORB0QgAQRAkiUTN/p9gepmvjgzsFYEfZcPMPccViDxW2rdurWbVmwHTMS2easdY7RY2aDzSf4Og+pnWjQE4A/tJyMrxnobrQXmJgKFEZBJL53hwAMPLOhbXbgWni0FgazCkWFMsr/D0GdwQ0SRRx55JCd7VHNVu4KkxzjUzXHORri8+5HGuGL6rFBg58xQyOU5wrTmLFm6xFMcq06YkX/88cediy66SK9EIa//r9Aq5FtvveUYDVdMH53f/va3TocOHdw6UDdWuWTQdtNQP8rZcuWVV7rnzz//fH1Klu/dNNMn8U3S56655hoHKwomHdvf/OY3dpWefVF6PXnlo9RzPt+BfS1o48knn8yXtWD67rvv7mnf9Puuu+4qWM4+KXHlPHWIYud069ZNpw0ZMsTOWnS/Z6+ebl3ycnVWrFhRtAwzEIEkETj22GPdZxC/B5nw8VT/wAMPOOJP6MmD90tUufjiiz11wKpAWOWiVuPmT/J36FbKndgICNuz5/7iWcJKddbClcSsEa+M9mSS2pHJK0dcNhysMlFqFwExK3Z23nnnnPdN9+7dQ1/0LbfcklMe7yxRPkPX4c+Y9BiH+jnO+VEu33GUMQ4rR3lF2EsDHz6jkGArcRD1UjleanZ6vv1WrVvlbQ8nzIeerIY5Qm6j82ILs0q0YcwfkWYUGrTlNy+DWSfMWnHu1FNP1fXAjPLtt992JDCyTkd9UHSvvvpqB2YdQpnvQDkyfe/Xr58uF/QPH4UwFTV5xX49KJsnDS9823wE5reyAunJE/ZAZhjdtk0fsIWyFlaAi13W7ANX8QUMW43OB1NVU/7MM8+MVJaZiUASCMiKnmfAxbtJYo86N9xwg+N/lwlDsTNu3LhYzWIixjzrZishOGLVhUJJ/g5jd4IFXQRkVj7n/s6aNcs9n9UOlcSskK6sdjCpjfdKnz59Kqtj7E1iCIgligMT0SAFEfce36aYtB85cqQj1nI57eJ9dPPNNzud8nwHog6xuNPjHyZL4Q4URZIe49A2x7kodyDdvFHGuIJKovkASnKLD7ZCAmUQ7UHBM0oi8otpZY6dNhQ8s+po/OnsuiXMhq4LZY1AKcMMHdoQYhQ9Q4zVBaEP1lnsFUIJfGyKBW4vu+wyXQ/qwsqmELcE5jOJWEG1sYRCHFcwS2SUYLtOYUoNXeVtt93m6Y+YvTr9+/ePrLjC39T0AYqvxIEM3QdmzI8AFPW7775bK/5iDql9T7EajudSTLLzF6zjM59++qnTq1evHKsA83xi4LziiivcyaY4UL0oA7Spz2ynT58epypdJqnfYewOsKAHAXx8mfuKLSwj4E+TtVBJzBrx8rdnTz6/9NJL5e8Qe5AKAueee67nHWO/b/z7Qf7Q++67b+jyqE/IESNdR9JjHBrnOBfpFqSaOcoYV9AnUXopz1e2ggDs8CdEOIWbbrzRDSCMmIz+WGhgbjr44IMV4g7CJ09MTjUjpumxiTFj+8fJi1eZkBVgDZOVQwU6/COPPFIXg6+iEfkhmt3ArXysa0pf+EbC1xHBu2UGJjAvEsXU1PXLzJsp5ImePXtqv0hgAl9IXAf2o/hWiUmwgn/lO++8o3784x8r+D+Iw3LIHvw3m3w8KXmJuWVkFlTX5SZwJxYCS5YsUSf36aNeEb9RW/Ccw69VZhj1tm3btvbput+XySL16KOPKsRBROge+DeDIVfMs/W7Ac97qdJa3ht+seOm+s8VO07id1isDZ4Pj4BMNHgywyc+jL+5pxAPiEAMBMSdRskkq5IJAgU/a0ptInDnnXcq/MUVMReMWzRUuaTHODTKcS4U9JlkijTGpaquxqhcwi64MyQwzYTZaCERwho3P8xFbcGKAWZcbIGvpNwF/YdVhZNOOsk+7QhBjj6H1cYws8fwHzArk5hxFipiT31ZHGD20ZgtTJ06NYsm3TaAscGzU6dOsZha3cq4oxHACqKQCGlcmzRpon1w4YNqVs0N3lhJlokBolYGBHbccUf3ucf9iOrDW4Yus8mQCBx11FGee1sKS2DIJgOzwS8J1i74u/HGGwPzMLG2EDj88MP1sydkerV1YbyaqkOAY1zV3bLQHY4yxhU0Nw3dYoIZofD85Cc/cQdpmJ/CuTefIOyDCUchjKNuNpQBjTCUSFtAemM+slu0aOFxDF+5aqU2LcJ5YYayixXch0mgqRNkO1kL/CHRPkKVZCkLFy50oMSg7ebNm9PMNCHwf3/HHRpTfDDAz84Ink+/mYohZTJ5uM0GAWEAdn/zeP7FoiGbhtlK6ggg9I95n2MLc3oKEUgbAZjKw48ak+N02UgbbdZfDAGOccUQqt7zUca49WUQrCiRl6S6d+xYhS0EpqcwYcwnoCZv06aNPg2zSZiWQRBCQ3wNlZDP6GP8++qrrzzU9bfffruyg6m/MP0FhUDbEPlA19sw/wYNGqR69Oyhsz777LNKmFLDFEskjzCwKmG3UkLGo0aNGpVInWEqgXktgr2uWrVKhxcBjTNMVimlIzBewrnAFBohXxDexUiTxk20iYr43JkkBfNpSvYIIASOLcLqbB9yv0oREL9y5TfFgUsDhQikjcA9Eips9erV6rDDDuNYmjbYrL8oAhzjikJUlRkij3GVqgvbZqFgeioU2NTOi1AYWF2EKam8dD2XBzZDuav6b//99/ecw4Eoe+55kN5EEaxcnnDCCbo8+uun349SV5i8WHEF+xWuB+auYUNwhKm7WB7x6dRmvGgbq7Xik1msCM+HRGDZsmUOCIQQ6iWfwMTUPMfIG8YsOl9dTI+HANhMzT0w23KYmsfrPUvlQwAhicz9xFbiheXLynQikCgCeNbwzMnkYCL1wm0BDJrt27d3ZAJXmy3DdPrbb78tqf606i2pUyELg+WzR48emi0f1k8ggislVETIZqsyG8e4qrxtRTsddYyrCHNTIZlwYCNrCz58d9llF3fARoyzfIKXqhnYoRgOGDDAEQKWnDAO8D80+YLqM2xyQUyp+dr2p4M1CPEW8fEeJwabv758xwjnAUUYZm5ZfpzC/BHYAkfcn1Ko//NdWz2ng+4akxWFBAy9eL5wDzBBUEwQNqbY39q1a4tVw/M+BOywL7gXftN2X3YeVgEC8Gs3YwS2UUIKVcHlsYsREBCSMOePf/xjSTFQwzb3/PPP6+cO75Sk3sUmLjQmre1nWiyvHLGqCtu1nHxp1ZvTUMIJEszd9eu3MYG70h3i4kHJRYBjXC4m1Z4SdYyrCHNT+TDWLKEwYTQCNjmxiTaH2qzRPfDtdOnSRckPXadOmDBB/eEPf1BCgKPkReDJKS8JfQzTTNtkD4lgk/zoo4/0ecMq9otf/EKBwTSKiH+geu7559TWW2+tJCCpGj9+fJTiofM2atRIyWqnmj17thK/zdDlSs0ohCq6XVHElcScVC1btiy1Spa3EACrnfi4Wim5uzCplheVPiEv8ZwMsqqthg4dqk1WJSSJaty4cdG/F154IacekyBkSOrKK69Uxx13nMJvDXX7zSthpj1v3jzdd5mt1vnwu4YgvXfv3kr8fHPKmTaqcYvfgi1pM87ZbdXbPkxA4R6A944QZKlTTjlF4V0vEyaJQgGGaFsO6Oi9x/Y57sdHQPzZlcQ21uzkeIeBpRluGw8//HD8ShMqiWftiCOO0Mzf+BYQaxnt9pJQ9YHV3DPmHp2Obx6ZAAzMEyVRJqzUnDlzlIQdUGKdotniDQs52NgLsbAXaietegu1mcQ5jJd4d8iqqhKrK/3ekAkAteuuu6o1a9YoiZWt1qxdm0RTNVUHx7jsbmfFjnGVoBWbwJ3iR+jpzrBhw9wZsNdff91zzn8gP3Y3b1D8QVFo3PNgjfOLKHPu+fvuu0/PLIky6YC9NI7IDdfx2qIGMY3TVpZlBg8enJg5TJb9rqW2YFItry79J/6vnkvDiqFNzmTyFdtiJtsvmEWXjzfdjoSI0Svjhn0YhEWi/LlFwISG34tpB+dhDYDfElh/Tbooi26Zat95+umn3evC9YkfabVfUsX1H88z4tHieRL/XEcm3hyZBHTN3bt165ZYn2FGB7M886xii3GDkiwCWLUxjOA21mb/jDPOyLECSrYH+WubMXOGfgbwvF177bUOjmGxUyxmcv4ai5+RcF+OTORp65CkLHNkMs8ZPXp0TuP4vQBnm+QvJ1OBhLTqLdBkIqewiggSRL9rBoLam1VFjHcULwIc47x4pHFU6WNcRZibmqD0EiPRMb6ACExvzBolRmLRe2P8AfFCtxkhTUFb4YRvol9uuukm9+MA5hj4sJVZTX82HhOBsiNgfi9QyuDDaAuC5pqPLbDkYR9bsFnhz6ThI82kYesfIBH4FoMnPpYefPBBtwkhZnLr94cGwMvOhOlAOBQEg0Z7sqruljn11FPduqp9B2ZhUFwM3tiC8ZeSDAIYA4wLAD5ObRM5MDkDb5hdl+pjZXorKy+ee4nJFkqyCGAyynyU472DEFV+kzbc1759+ybbcIjaMFkNVxG8V+fOnatLwPcfx0IQF6KGeFmEeE4/dxKnOF4FAaVOPvnkQLNViQmt25IV0oBSxZPSqrd4y6XlmDRpUl6WYvAq4JmbNWtWaY3UYGmOcene1GoY4ypCScTqgrCU6hk7fGTutdde7uoDnIwLkdaYW4hQAPih5yOMweoizmMQWLFihSnmbrFSiXPIIyZ6egXEPckdIlAhCMAX1XxUYVXPliVLl+jfzcCBA7WyAqUNzzNWA42YFfezzjrLJOVsoRSasDL+lUr48qJO/MH/1haE6ACFO86hfgy+/fv31z49UA7hyxJ3Zd5up5L2zz77bBePIEwqqa/V1Bc8u+3atdPYYvIQqy22mOcYSjo+5JOQX/7yl557KeatSVTLOn5AAKs4sDjAxBNWg20B8Zq4Lrj44z0CgrSsREyX9TsP7z2bNOy5557TfXrjjTdS64qZDA+avE66UTE51dcj5pWJVp1WvYl2MqAyPJNYxcWkaVLvkYBmqjqJY1w6t69axriKUBIRAN7MwGMLtqlnnnnGgeN4GMFsB4KNd+3aNUz2vHlAGjJt2jRcfLktAAAR2UlEQVTNjpo3E08QgTIiIH6yepDH6jrM42zBhxZWxI3g9wPFxWbyNeZ0fuXPlMFqioSO0eVat26dM3CCNAdm4TBlwr4tNoEUPuJh3gOltpbl1Vdf1VgBZ/wdzHiJidxuxJs1mPoVCjSAj/YrrrjCgblYUoKg9aZNTBRi0oOSHAIY54HvmDFjAisFS7Ztmv7oo48G5ks6EWO+eefBxNQIZvl33nlnvdpp0gpt4QawePHiQllyzpn3Byyg/KaQOZkTSMBqJe4B2g0jYNLGpCQmIAtJ1HoL1ZXluRdl9RB4lPrtmGWfs27LPKPm3cgxLpk7UC1jXEUoiaVCjmDyoDP2m96VWi/LE4FKQgA+gPiIatW6VahZdiFS0QMgTD8hmEwxK4QwJw0SsAxjMIAZ3yuvvBKUJW+af8YRfiD1IMYkErjh/mS5AlKL+L755pv6GQSeEjMuk0tE8HLDGIx2YbZNSRYB+BpCES8k9ofTkCFDCmVN5BwmGYwvtc1lgMmtDh06OFtttZXzySefFG0LjIF4fqIqG7C0wPNWyLKjaOMhM8BcGytmCJ0VRjBpjlV8YOCfkLTLR63XLlvu/QsuuEBbkEEZpuRHgGNcfmzinKmmMa50Gi15w5VTEMQdQWixlRdaObvCtolAaghImBjVs2dPJWbZ6tlnnlVill20rUWLFuk8Yt6lt2Kyp5nccADGU7/ILLgyLKdgrfUzm/nz+4/Fp8NNEnIozULpJtTwzoUXXuheHVhehazHPeZOdAQkdJBbCCymWcjEiRMVWIONgHmTkiwCEgdQ3XjjjQUrBbuzETAzFxKwOBf7s++pvy78VoXLQIn7iQ5eP3bsWJ1FVhaVWCYpMLBOnDRRtWjRwl/UPf7uu++UmNFrFnS0hXcg2ETDCNoHczREJtjCFCkpjyhESlYslazkFq1HlGfN8CqTJ6pjx445TPF2BVHqtcuVe19IBZW4Jynhq1AyeVHu7lR0+xzjkr09VTXGxdGCy1UGs3swMYKpHExBYCYCs6AsbPnLdc1slwiAmGO33XbTrHsSeiI0IPidyKvNQXxQCEyhcIy/oJh+xpQV54MYgAs1jFinpm5so/SzUL3VcA6+LMZPFNcOHyNKPATko1v7rJlnKSsfVpsR2PbhjXcVLBUXAazsmHvvJ46DDw/eaTBlL8SOaspjG8TabPrWr18/ty2ZJNBEW+BHALEO/GHxTiskeMfZsZxNu7bJaqHyMLtFmTgESegbSGiwOiuTeZroB77k+UxWwSgLV4NibOuIuQxCMvjpmeuR0EcOSH3gn+mXsPX6yyV1jDHtscce088FYpqa1We4QuA+wIpGwjflNIcxFc8RrF8oxRHgGFcco7A5qm2MqypzU3zYmhcXfAjgAJ+PqCbsDWM+IlDJCODDCD6FMImC+WgUMaajCFsCAXOw+f0E+QUhzIU5H5XNDx8Lpuzhhx8epZs1kRckPub6sYVvJyU6AvA1NDhiArCQmVv02oNLgLTMtIntqFGjgjMyNXUETKB2ENlhUtgI3oO2Im/fr0L7tpKIia9O4pOHdx9CWtjlzCQPfBChfPn9rU0/zBbvUryTYYaHOm1TZSiOYQTmrOhDPtP/oDqAQ9SQMPAFBcN0mPED5EH4rrKxwTH+YBprS5R67XJJ7YMVFsRVNvYg0IECCNZc+xrs0CJQpOFDCUb8LN4vSV1vuevhGJfMHai2Ma6qlET80M0Ml5iuOIjhQiECtYoAPlRA/Q/GXzjY5xMQPPg/bDFbZWbbf/3rX+uitgP60KFD3epQHis2UEbNwIrZ2SiC2WZTFkyB9SZi9uaJsdenT596gyCR6zVskniWwHKdhWA1xjy7+PDHhzglewTsdxZWFG2Bj6i5RyaMT9TQPngvBSmaUBDBoo6Vy7BKA37vmBiDnzfE74/tDylkXwv24e+HFUu828OGcIlDlw8yMyiI2PoFMSDziQntAyU4SOLWG1RXqWk2K7G4S2i/UOBqlF3gbEgQV69Z48DvFX/YtwWxrT/66CM7ifsWAhzjLDBK2K22Ma6qlETcFziR4wUc9mVewr1kUSJQNgTw8QEiBXwIYcY2n4BRD8QC48eP92RByAnzUQVTKghif5k0zLRC8EEMWnoMriA0MOfzmXDjw8If+gJ9BTsfym6xxRZ5TZ50gzX87/rrr3fxwwdKMUbAGoYi9qXBhM88g2DHDRKsBMA02l4dCMoXJg3PP55Z02aQaVqYepindAQQ0gf3QXwAPcyySYX2MT00AeXRFtxXkhAoelj5Ns8RTFkLyUUXXaTzwrw1jOA5jRoSBpOCMDH1jx9QjkaOHKlNMYPafvvtt93rOPPMM3OyxK03p6KEEsRfTvcXiiFIDDHRg0lPuCOJP6UjfrC6JXwzYvUQMSL9ZrmIrQ03AYyDlPwIcIzLj03YM9U2xlWdkhj2RjAfEahWBDCYwY8QHxwwq4YJVtAfQlTA1AZKojGPwgvImDGZDxYhsNFQYGDEQGrSW7VqpT/IsDqPOHS/+tWv3HNHHHFEDnww0wJ7pzhde85h0sbUCb+QehX4bZj4fcAjC8bCWsMaz7FZAQeGfhM5MClixRsKJFZWShX4LZlnFwzZmC2nZI8A2GUx0YRJMf+HehKhfcwVwfzU3G+YlvqVBZMvzrZv375u3VBU8j1LUNJMAPdCq3l2H2zW1zAhYcDaihVEXCsURfvPjAFCsGU34e4jjJLB6JFHHnHTsRO3Xii5YH0u5c+/8mc6Bh9i9BfWZRif7r77bnPKs4XCi3xQ5m08sI/VRjwPlMIIcIwrjE+Ys9U2xlFJDHNXmYcIZIjAwIED3UHaDNaFtsbnEF285JJLPGXbtGnj6TliHPnrgrkdBB/d9qrKgAEDtJkrSAtgTgqFNCgIs23u88ADD3jaq7cDfNDigwMYQ8EPQ59fbxgVu16sBphnFKaAIOSAiQ6e7R122EF/zIGwolTBBysmWExb9G8vFdH45WE1gfuQ7wPfrjlOaB9Tvn379u79fvLJJ01yItuXXnrJrRvXEuT3jYZg3o/zYZWSOHT5zZo18/TFPONmWyguI8J4IN+GG26YYwobt174bZq2426x+uqXv/3tb25YJ9R76KGH+rPo499bPvP52o/qhx/YUB0kcowr/SZX0xhHJbH0+80aiEBiCDz77LORB1OQKBjBhwpmUzEQQlnx++2CUMVWBDGo2gQROA/yAv9AilVLv9mSadMonlCK4FdU73LyySe7+MHslxINAcy0nnvuuZ6PPzyPWAE5/fTTE3vG7JXzffbZhy4M0W5TYrnNBxMmm8LIpEmT9O8LRFsQTG6Z95WEw8pbhYSncPPlUybyFg55Aivcpi8mPq2/qFHCwNQaRuwVyiBW6jB1hM2D1TqsguIaYN6alCShJAaZvmJV1eCN7bvvvptUl1lPAQQ4xhUAJ8Spahrj1sP1yI+LQgSIQI0gIM736qmnnlKHHHKIElOcnKtauXKlmjFjhhJzUyXmkTnnkSBmq0r8HXU8RWHrUxKCIzAfE3MRkI9WJUq1klluJbPxSnxBlXxw5WZkSkEEZMJByaSFWr58udqu2XZqn7b7KDFFLVgm7EncIzz/qFv8R9Vrr72m71nY8syXDAKjR4/WMQLFekJhP4wI6ZYSM2ElJvkK8S1ltV4hBiNElCglJDeB1YgJvZo+fbpq0KCBkuDpSnzQAvOVkoh+oX8QmaRT8+fPV4g5a2Tp50vVTjvupM+h30IQY04FbvGuFtNUtWrVKn1efO30cxuYOYFEcR1QMumnaxK/XyUMognUqtSs2bPVlxJzsRTBOOQfz0SBVg899JCuVnxN1RNPPFFKEywbEgGOcSGBKpKtKsa4EEovsxABIkAEiEAEBKZMmeKanYKl07AgRqiCWVNEAL6zMn7rvxEjRqTYEqvOhwBWBGHSeOKJJ0ZaxY0a2gftwzzZmIGDoMsvjz/+uLPNNts4YLgsRUCyg2syzxbCVdgCc32cC/L5tvOZ/azp8kHcZPru9w01faqkLfyITX/nzZtXSV2r+b5wjKvsW5zUGLe+/MAoRIAIEAEikCACQqmvLr74Yl2jsAUqUUQSrJ1VlYKAhHdRUyZP0VVIiBklvo6lVMeyMRB45plnlMTdUxJTVT04YYJeWfNXs3r1aoWVIqymGcE+VuYhEktRb80qGw6WLl2q0/Bv8uTJSoLH62OskMknnd4XRUhvzT/UJ2bMSlgvlfjAmuRY22bbNVOdO3d2y4oi7LaLRAnDoc+J6aSbp9AOVmyMCAt1IE7mfBLb2bLiB8GKPfCoZBESHbVkyRLdxVatWynxN63k7tZc3zjGVe4tTXKMo5JYufeZPSMCRKCKEZBQIUqCNusrEBp29fHHH1fx1dRG11esWKEkBp++GJgBCntjbVxYFV0FFBEofxKGR+FjpsEGG+T0HgrdoEGDlKykqaZNm7rnzzvvPCWEQ/rYmMALG7R7HsonRJhFlfg6ukojzIkhEidRt6sP5B9M6rt3765ktV8JQZJJLmkLhdOIrEwqYVTVhzDxX7hwoTYf/bm0GUaEqMnNJgyh7r69I/5N+lpRdykivukaD9QhPrqqYcOGCgr1ddddV0q1qZWFS4WRE3qfYHa5zRABjnEZgh2yqcTHuMpeMGXviAARIALViwCo/Y1JFJhmk6Tcr15UytNzkHIgRpqMtZrVFAHNKdki8Prrr7vEWeIzHRjaR1Z3HVkl1Pdp2LBhuoOlhPZBBbJiqesDc7QR9AXstogZC6bbpASkFCb8BJ41WXHRVffq1Uv3IYilM1/bqCurkDB2kG9RDHXsXZjgLlu2LF/3ypouK7YaT2Asfvhl7Us9N84xrnLufhpjHNlNK+f+sidEgAjUIAKILYZYXPiYQVgRSnkQkBUefQ+EqMYB0yUlWwTAwmx+B/gtFPsDSzPiYkJKCe2D8vANRHsItfHll186w4cPd/AcgIF05aqVyJKoyCqoe31oB7FqsYVfZNTJCcP+iv6nGRIGbKvmnoj5pg7hk4/ROlGwYlQGRm7gif4iNi2lvAhwjCsv/qb1NMY4KokGXW6JABEgAikhAAIKszoCcgtKtgiYIOoI/4I4X5TsEdhvv/1cJcQoI4W2wlbpdrLU0D5Q0uzwFBtttJFWHLFSl4YIe6rnWnfffXd93KFDh8jNoY9ZhIQRf0m3zwg4P27cuMh9zaqAHZNyyJAhWTXLdgogwDGuADgZnEprjGMIDBmlKESACBCBtBFAOAdZyVAgsmnRokXazbF+CwFZEVGDBw9WEjc0h0bfysbdCkag1NA+8FMU5UL9+9//VrKC6BLfpHXJ8HN85513PNWPGTNGhSWt8RSUgzTp8k1bYnKq4JsIfMRk1iRzSwRCIcAxLhRMqWRKa4yjkpjK7WKlRIAIEIFcBISKX2XBUpjbMlPEd0aJjxWBIAKZIHDrrbcqxBo0IqvYmo1TzGhNErdEoOYQ4BhXvluaxhhHJbF895MtEwEiQASIABEgAjWIAFgGEVLDsJKC9VRMOGvwSnlJRIAI1CoCDIFRq3eW10UEiAARIAJEgAiUBYHNN99cx4E0jZ999tlml1siQASIQFUgQCWxKm4TO0kEiAARIAJEgAhUEwLG/1CIa9T+++9fTV1nX4kAESACiuamfAiIABEgAkSACBABIkAEiAARIAJEwEWAK4kuFNwhAkSACBABIkAEiAARIAJEgAgQASqJfAaIABEgAkSACBABIkAEiAARIAJEwEWASqILBXeIABEgAkSACBABIkAEiAARIAJEgEoinwEiQASIABEgAkSACBABIkAEiAARcBGgkuhCwR0iQASIABEgAkSACBABIkAEiAARoJLIZ4AIEAEiQASIABEgAkSACBABIkAEXASoJLpQcIcIEAEiQASIABEgAkSACBABIkAEqCTyGSACRIAIEAEiQASIABEgAkSACBABFwEqiS4U3CECRIAIEAEiQASIABEgAkSACBABKol8BogAESACRIAIEAEiQASIABEgAkTARYBKogsFd4gAESACRIAIEAEiQASIABEgAkSASiKfASJABIgAESACRIAIEAEiQASIABFwEaCS6ELBHSJABIgAESACRIAIEAEiQASIABH4f9gx2r9HbtM/AAAAAElFTkSuQmCC" - } - }, - "cell_type": "markdown", - "id": "7edc1730", - "metadata": {}, - "source": [ - " ![image-2.png](attachment:image-2.png)" - ] - }, - { - "cell_type": "markdown", - "id": "7fd29983", - "metadata": {}, - "source": [ - "* l scattering mean free paths \n", - "* c is the Rayleigh wave velocity\n", - "* r is the distance between the source and receiver\n", - "* H is the Heaviside function" - ] - }, - { - "cell_type": "markdown", - "id": "8f06cd85", - "metadata": {}, - "source": [ - "The sum of squared residuals (SSR) between *Eobs and Esyn* was calculated, and the values of the parameters l and b that minimize SSR through the grid search were determined.
\n", + "The 2-D radiative transfer equation for scalar waves ***(Shang and Gao 1988; Sato 1993)*** assuming isotropic scattering and source radiation in infinite medium to calculate ***synthesized energy densities Esyn*** :\n", + "\\begin{gather*}\n", + "E_{syn}(r,t)=\\frac{e^{-ctl^{-1}}}{2{\\pi}cr}\\,\\delta(t-\\frac{r}{c}) + \\frac{e^{l^{-1(\\sqrt{(c^2t^2-r^2)}-ct)}}}{2\\pi l \\sqrt{(c^2t^2-r^2)}} H(t-\\frac{r}{c})\n", + "\\end{gather*}\n", "\n", + "* $l$ scattering mean free paths \n", + "* $c$ is the Rayleigh wave velocity\n", + "* $r$ is the distance between the source and receiver\n", + "* $H$ is the Heaviside function\n", "\n", - "An energy density with the spatially homogeneous intrinsic absorption is obtained by multiplying e^-bt with the right-hand-side of Equation (1), where ***b is intrinsic absorption parameter***\n" - ] - }, - { - "attachments": { - "image-3.png": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqYAAACGCAYAAAAPSEiAAAAMQGlDQ1BJQ0MgUHJvZmlsZQAASImVVwdYU8kWnluSkEBCCSAgJfQmiNQAUkJoAaQXwUZIAoQSYkJQsaOLCq5dLGBDV0UUOyB2xM6i2LAvFlSUdbFgV96kgK77yvfm++bOf/85858z587cewcAzRNcsTgP1QIgX1QoiQ8LYoxOTWOQngIcGAEaoANTLk8qZsXGRgFYBtq/l3c3ACJvrzrJtf7Z/1+LNl8g5QGAxEKcwZfy8iE+AABexRNLCgEgynnLSYViOYYV6EpggBDPl+MsJa6S4wwl3qOwSYxnQ9wCgJoGlyvJAoB2GfKMIl4W1KD1Quwi4gtFAGgyIPbPzy/gQ5wOsR20EUMs12dm/KCT9TfNjEFNLjdrECvnoihqwUKpOI875f9Mx/8u+XmyAR82sGpkS8Lj5XOGebuZWxApxxoQ94gyomMg1oH4g5CvsIcYpWTLwpOU9qgxT8qGOQP6ELvwucGREBtDHCrKi45S8RmZwlAOxHCFoJOFhZxEiA0gni+QhiSobDZKCuJVvtD6TAmbpeLPcSUKv3Jf92W5SSyV/utsAUelj9GKsxNTIKZAbFUkTI6GmAaxszQ3IVJlM7I4mx09YCORxcvjt4I4XiAKC1LqY0WZktB4lX1ZvnRgvtjGbCEnWoX3FWYnhivzg7XwuIr44VywywIRK2lARyAdHTUwF74gOEQ5d+yZQJSUoNL5IC4MileOxSnivFiVPW4hyAuT8xYQu0uLElRj8eRCuCCV+nimuDA2URknXpzDjYhVxoMvAVGADYIBA8hgzQAFIAcI23oaeuCdsicUcIEEZAEBcFIxAyNSFD0ieE0AxeBPiARAOjguSNErAEWQ/zrIKq9OIFPRW6QYkQueQJwPIkEevJcpRokGvSWDx5AR/sM7F1YejDcPVnn/v+cH2O8MCzJRKkY24JGhOWBJDCEGE8OJoUR73Aj3x33xKHgNhNUVZ+LeA/P4bk94QmgnPCRcJ3QSbk0Qlkh+inIU6IT6oapcZPyYC9wGanrgQbgfVIfKuD5uBJxwd+iHhQdAzx6QZavilmeF8ZP232bww9NQ2ZFdyCh5CDmQbPfzSJoDzWNQRZ7rH/OjjDVjMN/swZ6f/bN/yD4ftpE/W2Lzsf3YWewkdh47gjUABnYca8RasaNyPLi6HitW14C3eEU8uVBH+A9/A09WnkmpS61Lt8sXZV+hYLL8HQ3YBeIpEmFWdiGDBb8IAgZHxHMexnB1cXUDQP59Ub6+3sQpvhuIfut3bs4fAPgd7+/vP/ydizgOwF4vuP0PfefsmPDToQ7AuUM8maRIyeHyCwG+JTThTjMEpsAS2MH5uAJP4AsCQQiIADEgEaSC8TD6bLjOJWASmAZmg1JQDpaAlWAt2AA2g+1gF9gHGsARcBKcARfBZXAd3IGrpwu8AL3gHfiMIAgJoSJ0xBAxQ6wRR8QVYSL+SAgShcQjqUg6koWIEBkyDZmDlCPLkLXIJqQG2YscQk4i55F25BbyAOlGXiOfUAzVQHVRE9QGHY4yURYaiSai49AsdCJajM5FF6Gr0Wp0J1qPnkQvotfRTvQF2ocBTB3Tx8wxJ4yJsbEYLA3LxCTYDKwMq8CqsTqsCT7nq1gn1oN9xIk4HWfgTnAFh+NJOA+fiM/AF+Jr8e14Pd6CX8Uf4L34NwKVYExwJPgQOITRhCzCJEIpoYKwlXCQcBrupS7COyKRqE+0JXrBvZhKzCFOJS4kriPuJp4gthMfEftIJJIhyZHkR4ohcUmFpFLSGtJO0nHSFVIX6YOaupqZmqtaqFqamkitRK1CbYfaMbUrak/VPpO1yNZkH3IMmU+eQl5M3kJuIl8id5E/U7QpthQ/SiIlhzKbsppSRzlNuUt5o66ubqHurR6nLlSfpb5afY/6OfUH6h81dDQcNNgaYzVkGos0tmmc0Lil8YZKpdpQA6lp1ELqImoN9RT1PvUDjU5zpnFofNpMWiWtnnaF9lKTrGmtydIcr1msWaG5X/OSZo8WWctGi63F1ZqhVal1SKtDq0+brj1CO0Y7X3uh9g7t89rPdEg6NjohOnyduTqbdU7pPKJjdEs6m86jz6FvoZ+md+kSdW11Obo5uuW6u3TbdHv1dPTc9ZL1JutV6h3V69TH9G30Ofp5+ov19+nf0P80xGQIa4hgyIIhdUOuDHlvMNQg0EBgUGaw2+C6wSdDhmGIYa7hUsMGw3tGuJGDUZzRJKP1RqeNeobqDvUdyhtaNnTf0NvGqLGDcbzxVOPNxq3GfSamJmEmYpM1JqdMekz1TQNNc0xXmB4z7Tajm/mbCc1WmB03e87QY7AYeYzVjBZGr7mxebi5zHyTeZv5ZwtbiySLEovdFvcsKZZMy0zLFZbNlr1WZlajrKZZ1VrdtiZbM62zrVdZn7V+b2Nrk2Izz6bB5pmtgS3Htti21vauHdUuwG6iXbXdNXuiPdM+136d/WUH1MHDIduh0uGSI+ro6Sh0XOfYPowwzHuYaFj1sA4nDSeWU5FTrdMDZ33nKOcS5wbnl8OthqcNXzr87PBvLh4ueS5bXO6M0BkRMaJkRNOI164OrjzXStdrblS3ULeZbo1ur9wd3QXu691vetA9RnnM82j2+Orp5SnxrPPs9rLySveq8upg6jJjmQuZ57wJ3kHeM72PeH/08fQp9Nnn85evk2+u7w7fZyNtRwpGbhn5yM/Cj+u3ya/Tn+Gf7r/RvzPAPIAbUB3wMNAykB+4NfApy56Vw9rJehnkEiQJOhj0nu3Dns4+EYwFhwWXBbeF6IQkhawNuR9qEZoVWhvaG+YRNjXsRDghPDJ8aXgHx4TD49RweiO8IqZHtERqRCZEro18GOUQJYlqGoWOihi1fNTdaOtoUXRDDIjhxCyPuRdrGzsx9nAcMS42rjLuSfyI+GnxZxPoCRMSdiS8SwxKXJx4J8kuSZbUnKyZPDa5Jvl9SnDKspTO0cNHTx99MdUoVZjamEZKS07bmtY3JmTMyjFdYz3Glo69Mc523ORx58cbjc8bf3SC5gTuhP3phPSU9B3pX7gx3GpuXwYnoyqjl8fmreK94AfyV/C7BX6CZYKnmX6ZyzKfZfllLc/qzg7IrsjuEbKFa4WvcsJzNuS8z43J3Zbbn5eStztfLT89/5BIR5QraikwLZhc0C52FJeKOyf6TFw5sVcSKdkqRaTjpI2FuvBHvlVmJ/tF9qDIv6iy6MOk5En7J2tPFk1uneIwZcGUp8Whxb9NxafypjZPM582e9qD6azpm2YgMzJmNM+0nDl3ZtessFnbZ1Nm587+vcSlZFnJ2zkpc5rmmsydNffRL2G/1JbSSiWlHfN8522Yj88Xzm9b4LZgzYJvZfyyC+Uu5RXlXxbyFl74dcSvq3/tX5S5qG2x5+L1S4hLREtuLA1Yun2Z9rLiZY+Wj1pev4KxomzF25UTVp6vcK/YsIqySraqc3XU6sY1VmuWrPmyNnvt9cqgyt1VxlULqt6v46+7sj5wfd0Gkw3lGz5tFG68uSlsU321TXXFZuLmos1PtiRvOfsb87earUZby7d+3Sba1rk9fntLjVdNzQ7jHYtr0VpZbffOsTsv7wre1VjnVLdpt/7u8j1gj2zP873pe2/si9zXvJ+5v+6A9YGqg/SDZfVI/ZT63obshs7G1Mb2QxGHmpt8mw4edj687Yj5kcqjekcXH6Mcm3us/3jx8b4T4hM9J7NOPmqe0Hzn1OhT11riWtpOR54+dyb0zKmzrLPHz/mdO3Le5/yhC8wLDRc9L9a3erQe/N3j94Ntnm31l7wuNV72vtzUPrL92JWAKyevBl89c41z7eL16OvtN5Ju3OwY29F5k3/z2a28W69uF93+fGfWXcLdsnta9yruG9+v/sP+j92dnp1HHwQ/aH2Y8PDOI96jF4+lj790zX1CfVLx1OxpzTPXZ0e6Q7svPx/zvOuF+MXnntI/tf+semn38sBfgX+19o7u7XoledX/euEbwzfb3rq/be6L7bv/Lv/d5/dlHww/bP/I/Hj2U8qnp58nfSF9Wf3V/mvTt8hvd/vz+/vFXAlX8SuAwYpmZgLwehsA1FQA6PB8RhmjPP8pCqI8syoQ+E9YeUZUFE8A6uD/e1wP/LvpAGDPFnj8gvqaYwGIpQKQ6A1QN7fBOnBWU5wr5YUIzwEbo79m5GeAf1OUZ84f4v65BXJVd/Bz+y8giXxlIATszwAAAIplWElmTU0AKgAAAAgABAEaAAUAAAABAAAAPgEbAAUAAAABAAAARgEoAAMAAAABAAIAAIdpAAQAAAABAAAATgAAAAAAAACQAAAAAQAAAJAAAAABAAOShgAHAAAAEgAAAHigAgAEAAAAAQAAAqagAwAEAAAAAQAAAIYAAAAAQVNDSUkAAABTY3JlZW5zaG90i8LtvQAAAAlwSFlzAAAWJQAAFiUBSVIk8AAAAdZpVFh0WE1MOmNvbS5hZG9iZS54bXAAAAAAADx4OnhtcG1ldGEgeG1sbnM6eD0iYWRvYmU6bnM6bWV0YS8iIHg6eG1wdGs9IlhNUCBDb3JlIDYuMC4wIj4KICAgPHJkZjpSREYgeG1sbnM6cmRmPSJodHRwOi8vd3d3LnczLm9yZy8xOTk5LzAyLzIyLXJkZi1zeW50YXgtbnMjIj4KICAgICAgPHJkZjpEZXNjcmlwdGlvbiByZGY6YWJvdXQ9IiIKICAgICAgICAgICAgeG1sbnM6ZXhpZj0iaHR0cDovL25zLmFkb2JlLmNvbS9leGlmLzEuMC8iPgogICAgICAgICA8ZXhpZjpQaXhlbFlEaW1lbnNpb24+MTM0PC9leGlmOlBpeGVsWURpbWVuc2lvbj4KICAgICAgICAgPGV4aWY6UGl4ZWxYRGltZW5zaW9uPjY3ODwvZXhpZjpQaXhlbFhEaW1lbnNpb24+CiAgICAgICAgIDxleGlmOlVzZXJDb21tZW50PlNjcmVlbnNob3Q8L2V4aWY6VXNlckNvbW1lbnQ+CiAgICAgIDwvcmRmOkRlc2NyaXB0aW9uPgogICA8L3JkZjpSREY+CjwveDp4bXBtZXRhPgpaMOcfAAAAHGlET1QAAAACAAAAAAAAAEMAAAAoAAAAQwAAAEMAABZl7N20PgAAFjFJREFUeAHsnQnQTtUfx3+2LIlC9i1lp0SKyJZCEmWkUSghskQzlckyoSEla9YSIYaRJBmmZEkZS5YiWUIKEZE1y+v+z/fM3Pu/z33v87z3We6zvL5n5n2fe8859yyfs9zfPcvvZDGUERoSIAESIAESIAESIAESSDCBLBRME1wCjJ4ESIAESIAESIAESEAToGDKikACJEACJEACJEACJJAUBCiYJkUxMBEkQAIkQAIkQAIkQAIUTFkHSIAESIAESIAESIAEkoIABdOkKAYmggRIgARIgARIgARIgIIp6wAJkAAJkAAJkAAJkEBSEKBgmhTFwESQAAmQAAmQAAmQAAlQMGUdIAESIAESIAESIAESSAoCFEyTohiYCBIgARIgARIgARIgAQqmrAMkQAIkQAIkQAIkQAJJQYCCaVIUAxNBAiRAAiRAAiRAAiRAwZR1gARIgARIgARIgARIICkIUDBNimJgIkiABEiABEiABEiABCiYsg6QAAmQAAmQAAmQAAkkBQEKpklRDCI7d+6UVatWye7du2X//v3y2WefSf78+ZMkdUwGCZAACZAACZAACfhPgIKp/4yDxnD9+nWZM2eODBs2TA4cOCD169eXdu3aSZs2baR06dJBn0sWhz/++EMaNmwohw8ftpI0duxY6dOnj3VvXrzzzjsCt1OnTmmrbNmySY8ePWT8+PGmF/6SAAmQAAmQAAnc4AQomCaoApw+fVoLoOvWrZN8+fLJhAkTpHPnzglKTXTRIu2vvPKKDqRSpUp61NctRAilRYsWlbS0NPn222+lUaNGbt5oRwIkQAIkQAIkcIMSoGCagIL/559/pF69evLrr79KiRIlZP369VK2bNkEpCQ2USIflStXtgLDkoQmTZpY9+bFlStXJG/evNK4cWNZuXKlac1fEiABEiABEiABEtAEKJgmoCJgZHT27NmC6ew1a9boKfwEJCNmUY4bN0769+9vhde2bVtZtGiRdW9erF27Vo+STps2Tbp3725a85cESIAESIAESIAENAEKpnGuCCtWrJAWLVroWAcNGiTDhw+PcwpiH13z5s3l2rVr+g/CZ/bs2eXQoUN6NNge28CBA2XEiBFy5MgRKV68uN2J1yRAAiRAAiRAAiQgFEzjWAkuX74sFSpU0JuFsmTJIgcPHpQcOXIINhFhirtKlSoC+1Qyly5dkgIFCsjQoUPljjvukKefflonf8iQIdrOnpfatWtr4XXbtm12a16TAAmQAAmQAAmQgCZAwTSOFWHTpk3ywAMPWDFCUNu8ebN1j534r7/+uvTq1cuyS/YLcwR469atUr16dSlTpowcPXpUihUrpgVwjJ7CYONT4cKF5Y033tCjpsmeL6aPBEiABEiABEgg/gQomMaR+YcffmitrcT60gEDBmhBFdPeU6ZMsXazY3of0/ypYPr16yfz5s2T48eP69FejJy+9dZbOukLFy7U6q9wg+v27dvLd999l/JranXm+I8ESIAESIAESCDmBCiYxhxp8AB79+4tkyZN0h4+/fRT6dChg+X5zJkzUrduXb1TP3fu3Hqav0iRIpZ7sl5UrFhRatasKfPnz9dJPHbsmB41vXr1qt7otHr1am3ftWtXWbx4sfz9999601ey5ofpIgESIAESIAESSBwBCqZxZA81SdiFj9HS8+fPS65cuQJix6jpyy+/rO2mT58u3bp1C3BPthuskS1XrpzMmDFDunTpYiUPI6MYIYXZtWuXXjuLZQoQvBcsWGD54wUJkAAJkAAJkAAJ2AlQMLXT8Pm6VatWsmzZMr3R6dy5c+li27dvn94cBQdM848cOTKdn2SyMAVpnPxUqlQpK2k4NAAnQsFgvSxGiqHndNasWSl7iICVOV6QAAmQAAmQAAn4RoCCqW9o0wfct29fmThxonbAbnbniCmEVZwCBYORU3PaX1sk4b/WrVvLnj179PIDZ/KwEWrnzp06P9jwhDWzmOZPheUJzrzwngRIgARIgARIID4EKJjGh7OOZfLkydaO+71790r58uUDYr9w4YIeTYXl+++/L6+++mqAezLdQPVVoUKF5Pnnn7eEbXv6pk6dKj179rSsatWqJVu2bLHueUECJEACJEACJEACTgIUTJ1EfLw/ceKElCxZUrAxaMyYMQGnJSFa+9GeUMPUrFkzH1MTXdDLly+Xli1bytKlSwVLFJwGa2hx3OrZs2e10+DBg2XYsGFOb7wnARIgARIgARIgAYsABVMLRXwu2rVrp4/rrFGjhtZhaur5ROyjRo3Sa0urVasmO3bskKxZs8YnURHE0rFjR5k7d66cPHlSChYs6BqCfenChg0bpE6dOq7+aEkCJEACJEACJEACIEDBNM714PTp09KkSRPZvn27QIUSNhBBOMVoaYMGDeTixYt6g1SjRo3inDLv0eGggHr16klaWpo+XhTK9N2MOQIMwRWjxcksaLuln3YkQAIkQAIkQALxJUDBNL68dWwYZezUqZNWHZUzZ069Yx3HdOK4Uiirr1q1agJS5S1KnHc/evRouXLlin4gT548ei0p7NxM06ZN9YYn6G2lIQESIAESIAESIIFQBCiYhqLjsxs2EG3cuFGvw8TxpJlxxzr0mMIks7DtczEzeBIgARIgARIgAY8EKJh6BEVvJEACJEACJEACJEAC/hKgYOovX4ZOAiRAAiRAAiRAAiTgkQAFU4+g6I0ESIAESIAESIAESMBfAhRM/eXL0EmABEiABEiABEiABDwSoGDqERS9kQAJkAAJkAAJkAAJ+EuAgqm/fBk6CZAACZAACZAACZCARwIUTD2CojcSIAESIAESIAESIAF/CVAw9ZcvQycBEiABEiABEiABEvBIgIKpR1D0RgIkQAIkQAIkQAIk4C8BCqYR8MXxmiNGjIjgydg+Yp6qFKtQkyFfsc5TrNgwHBIgARIgARIgAf8JUDCNgPHp06elRo0acvjw4Qiejt0jhmHELjAVUjLkK9Z5iikgBkYCJEACJEACJOArAQqmEeL9/vvvpWHDhpKWlhZhCNE/5ocQl+h8+ZGn6EkzBBIgARIgARIggXgQoGAaBeXhw4fLkCFDMgyhTJkyGfpxerh27ZpcuXJFj2Li2s34JcQlMl9+5cmNH+1IgARIgARIgAT8IbB161bZtGmT7N+/XypUqCB169aVatWqSZYsWUJGSME0JJ7QjtevX5eHH35Y1qxZE9Rj6dKlZfv27XLbbbcF9RPKAYLa8ePHZe3atTJy5EjZsWOH5d0vIS6R+fIrTxY0XpAACZAACZAACfhG4Ny5c9KnTx/55JNP0sXRrl07mTt3rtx0003p3EwLCqYmiQh/jxw5Ivfcc4+cOnUqaAht27aVRYsWBXX36gChrV+/fjJhwgT9iJ9CXKLy5WeevHKmPxIgARIgARIggcgIdOrUSebMmSNVq1aVSpUq6f04W7ZsEfP93rp1a1myZEnQwCmYBkXj3WHp0qUC0KHM1KlT5aWXXgrlxbNb586dZfbs2VYhe34wTI+JyJfbF1aYyaZ3EiABEiABEiCBBBBYt26d3n8zZswY6d+/v5WC9evXyxNPPKGXJ2Iq/6+//pLChQtb7vYLCqZ2GlFcY9j6gw8+CBpC7ty5ZfPmzfoLIqgnjw4nT56U8uXL6wL2+EjE3uKdr0KFCkWcVj5IAiRAAiRAAiSQOAK9e/eWEydOyMKFC9MlYtKkSQJ3mPnz58szzzyTzg8sKJi6Ygnf8r///pM6deoErAF1hoJFvxBOc+XK5XQK+37YsGGeNl6FHbDjgXjnyxE9b0mABEiABEiABFKEQKNGjQQCKKbxnWbv3r1SsWJFbY2BvF69ejm96HsKpq5YIrPcvXu33HfffXLx4sWgAfTo0UOmTJkS1N2rwy+//CJVqlTx6j0qf/HMV1QJ5cMkQAIkQAIkQAJJSeDo0aNSokQJnbbly5dLixYtXNNJwdQVS+SWM2bMkK5du4YMYPHixfLkk0+G9JNsjpk1X8nGmekhARIgARIggcxIwFx/euutt8qxY8eCzh5TMPWh9Nu3b++6vsKMCqqjoPapVKlSplVK/KZivqC24pZbbnHle+DAAdmzZ0/QrzbXh2hJAhEQCFUPIwgu7EfOnj0r+fLlC/u5VHuAbTrVSiz10pvothwtsXD7gnD9h0pf3759ZeLEiXoZ4tChQ4N6jVow/eabb2TDhg2C6V6oTCpWrJhAdydORcJaAyiH79Kli1YZMHjw4KAJgUOkYZ05c0agiiAjkzNnTilSpIgULVrU107633//1UeWHjp0KGiSHnroIVm9erVky5YtqJ9kc0i1fK1cuVI6dOggn3/+uTRo0EAuXbokONnq66+/lmXLlgmWQ9x1112yb9++ZEOdqdKDr2QcFhGtwRpttN1UM856GI/0X7hwQSu1hio76EEuW7asvPDCCzJu3LiQ+gPjkbZYxsE2HUua3sO6Udt0Itqy91Jx9xlNXxDLvgMjpFAdVbx4ca3bHfJYUKP0SkVkfv75Z6NJkyY4rN36U0KWgT/TTn2hGzfffLO+79mzZ9B4og1LCRmGOl3JiteMP9Rv9erVjbfffttQQm3QdEXj8MMPPxjZs2cPmSZ1alQ0USTk2VTIlzogwBgxYoSRNWtWo1WrVoYSigz1lWYohb7pykMJpgnheKNEqha7G5UrVzZy5MiRjn2o9unm9tVXX6UUNrd6GI8MrFixwkDfC4Y1a9Y01AiPbg9KRYuhDgQx1AhIPJLhexxs074jdo3gRmzTiWrLrgUQhmUs+gK8S2PRdyjF+oaawjfUIGaGOYAuzLCNWsBqKP1T1otGjUYZatRUCwCXL182tm3bZrRp08ZyRwf50UcfucYTy7DUCUwBcSLe/PnzG++9956hVBPozlkdixXgR43CGIcPH3ZNW7SWKFC3F6xpByFenegUbTRxfz6Z84WX7lNPPaW5P/bYYwbqI4ya4jM2btxojBo1ylBfala5UDCNT/VRo+2G0mFncTfbgFIXYkyfPt1Q+muNjz/+2Jg8ebJupx07dtQCrekPv0rvXXwSG4NYgtXDGASdYRCPPvpoAGelj1g/M2vWLP2Cuffee1OKZbAMs00HIxMf+xulTSeyLUdbkrHqC6LtO0aPHq0HhlatWuUpSxEJpvjqNl8Yd999t6FUCqWLLC0tzVBK5y1/EFbdTCzDQvhqKYEVJ9I4c+bMgGghqDz77LMBftTQsh5VCPAYgxswcI4qm9zM35IlSxpqCUQMYotfEMmar9OnTxtKU4Eu2+bNm7vWS1B65JFHrPKnYBq/eoMZCrPe4xcfhahLoYz5EYR2kirGaz30Kz+o+3bOSrG1FRU+0uFWrlw5Q53uZtmn+gXbdGJKMLO36US35WhLNZZ9QaR9h9rsrWfM8OvVhC2YQoiyd3rTpk0LGhekY/jFCBWmU50mlmEhbKXUNSBtmMpVyuid0erp+4IFCwb4/fLLL9P5i4UFOv/bb789IC47P1xjdDnVTDLmC6Ns4KnUURjnz58PitT+YULBNCimmDtgdNRe9998880M48BHL5YB4CM3VYzXeujMD6YLMbppjvI73b3e//jjjwYGDNSmP0MdkJHuMaWiRZcDfjOLYZtOTElm9jad6LYcbanGui8It+9Q+zkMpbfdwIhrOCZswVQt/g14uSDiYAajIVhjWrt2bVcvsQwLESxYsCAgbffff79rvLB0fkmo0wiC+o3WAWtg7S9kt2ulbDbaaOL+fDLl64svvrAYz5s3LyQLs7NBOVAwDYkqpo6YmbDXfbWBwlP4ELCGDx/uyW+iPYVTD51p7d69u+bz+++/O51ieo81/fhoR1lgGUVmMGzTiSnFzNymU6EtR1vq4fYF4fjHfhTIf+PHj3dNJpbWPffcc65uYQumSqN/wMtFKYt3Ddi0xKakYBufYhkW4nvxxRcD0qa0AJjJSPfbsmXLAL/48vPT9OvXLyA++wsa1/iq+Omnn/xMgi9hJ0O+MPKudmtrvkrbQYb55EssQ0Qx94C6ba/zWPt99epVT/Eo7R6GVyHWU4A+eQq3HtqT8dprr2k+SqOJ3dq3a7OvhNB/8OBB3+KJV8Bs0/Ei/f94MnObTqW2/P8Siewq3L7Ai//t27frzd9NmzbVfTf6b/Nv9erVxpgxY/Q+JewtcDNhC6Zz584NeMGo46UMpY7ALWxth7U/wYZxYxkWIlN6QQPSZl9b5Uygcy2q3zvkMSWJTQf2l7PzGusjQ7F05iEZ7pMhX0ollMVVqRzLEEs0LzGMZm3dutVQamoyjCeYByymV+rNDKU+I5iXTGdvrk8y63zbtm1d86jUy6Xs2sdw6yEAYPoefY/JRR3O4col1pZ//vmn3giFeJV+4lgHH/fw2KbjjlxvKjbrLX4zU5tOpbYcbcmH2xdk5F/pBg/YHG+vI/Zr7NBXp2S6Jj9swRSqIuyB4xqFGImJZVhKJ2VAukKNyKijsAL8YloLOzz9NigwU32Wk6F5361bN7+TEfPwE5kvCIomO4w6u23Ec2Y4nJcYRvUgVEHzBBqSGRfWPeJDQx0x62nECQIIZgigFQKqN8xwlF5dA+t24LZp0ya9NhbCGf7wTGYx+HI284zfYFo68AWNckT+U8lEUg8xldW4ceN0XPDhEumHD+or+sKFCxcaSoF1yI8fc6Mg9gBghCiVDdt0/Esvs7bpZGnL0Zaon31BqL6jRo0aAX2avd+3X/fq1StoFsMWTBGSUwUBIsOUbiQmVmGNHTs2AAZUBrmZ3377zShQoECA33gKgzNnzgyI215Q5jVeKqlmEpUve7njJe/FeH2JYT1NrVq1dHnhQwebdbB+FUtElKJgqxzz5MljhFojDD2SEGxRvnnz5jXeffddrRrJuc7ZLH/zN9hMg5c8JpMfzALYVXQhf/jqdhps+oHgnoqbcsKphxglwLp7s5zdfgsVKuTEE/QedbJZs2aanV13Mj6A3DZ/mgFhdNaMGzpBU9mwTce39DJzm05kW462FOPVF/jdd0QkmEK4MxU4mx0bfr3ssnWCj1VYzpe8U1sARgQGDhyod6ra0/z44497XuvmTHuk9/ZpAntazGuMzKlToyINPmHPJSJf9evXt16uw4YN85R3Ly8xjDqZwhRUFe3atSsgbAibdhU1KLtggqS5qQV+lixZYoWDzYGmzlWz7CFM4GsUf3a/1kMpeOGcocDhFk6DEUKzXEIJ+c7nkuU+nHqIujNgwAAD60rBwix7qHCCHf6gJsurOajWiE6YMMFw7pBWpz6FDAJ6Y8248QGWysasO8hPsA2NbNOxK+HM3KYT2ZajLaF49QV+9x0RCaaAh9349q9zs4PDtGe4JtqwMH2LUSszDfjFiSfY8YXphjvvvDPdyTM4IAAbtxIxZQjFxHgJ2dPrvH7wwQfjLjCHW25O//HOF07tMncXg5/XkeaMXmKYQjdHOBEuRoPdDEa+7KrAcLKU81QLrCU1yxZ+nQYHTNhPS8MHVmYz6nxkiwFYYKoesxT46oaaH2xEtC+TQOeaSibSeog82uvZoEGDoso26q19qVBGs1g7duwIKBe3UeyoEhTHh9mm4whbRZVZ23SytOVoS9PvvsDvvuN/AAAA///yx5FAAAAg8klEQVTtnQvcZVP5x9cYIZQ7DTHUkNQgoWLGtSmXmYmkcS1RhJlqKhIZ91suRYnJZSIJTR+aT417imQaZISIJsRg3NXQSJz//m7/Z1tnv/ucs8++ncv7ez6f99377L3Ws9b6rct+1rOe9SxXy0EXXXRRbciQITXnXN3fhRde2DbXPLyuv/76uvTj+bHfSy+9dO3LX/5y7corr6y98sorbeexyAh/+tOfam9729sa5nuJJZaoPfbYY0UmWQmvKst1zz331OF3zTXXpCrj3nvvHcUbMWLEgDg/+clPovcrr7xy7b///e+AMPbg+OOPj8LSzmhfPp1xxhnR+w022MB/Fd3vuOOOUZhll102et4vN+973/ui8llfbHRdb731eq7YWdvhggULaosttliEze9+97tcZX/66acjXuD7q1/9qim/Rx99tC78zTff3DR8N79Un662dvq1T3dLX85bm2WPBWWPHS4vAN///vfrBjcGxEUXXbSGgNIuZeV1yCGH1OXh8MMPr1Ex06dPr3v+7ne/u/b666+3m63Swp9yyil1+fM/1j/+8Y9LS7dsxlWV69prr63D749//GOqorX6iO20004R36222qopzxdeeKHGhMfqbqmllqr95z//ieLstdde0bvVVlsteu7fnHnmmVEY+Lz66qv+656+f/jhh+vKtvjii9fuv//+2ty5c2sPPfRQ7fbbb6/9/Oc/r62//vphOPpyr1HWdnj11VdH2Cy55JK56/2qq66K+C2yyCI12mYzeu6556LwtLuLL764WfCufqc+XV319HOf7pa+nLc2yx4Lyh47cgumADhlypS6AY5Bbu211669/PLLbeObhZd91Ew4uOWWW6J00YjZc65ptWoRgxJv0MS94x3vqMsfefzc5z5XYqrls66qXL5mE9zuvffeVIVr9RGj7VqbSVMX2267bRSeeH/961+jfOy+++7Ru6FDhyYKH2i2LD0mdW+88UYUv9dvpk6dGpWNMn784x9PLNJ5550XhsurNUxkXvLDrO3wG9/4RoTNJz/5ydy5/OY3vxnx+/CHP9ySH/3U2h3X4447rmUcP8BnPvOZ2nve857C/jbffHOffVv36tNtwZUrcD/36W7py7kqKIhc9liQd+xoVb5CBFMS+dKXvlQ3yDHQZdX6tcPriSeeqEsXQe+1116Lyh1fat1tt92id52+OeKII+ryDmYjR47MJNB3uix++lWV67TTTqvD7x//+IefjYb3zT5imHigbbIP9oQJExrysRdf/OIXo/DE+81vfmOvanEsLrvssuid3WD6Yuk1Wu63sL12/fSnPx2VjTJi2pBE3/3ud2vLLLNMXd9NCteNz7K2ww033DDC5tRTT81dtE022STih9Cbhnxzov333z9NlCjM8OHDo/Ss/ea5rr766hHvdm/Up9tFLHv4fu7T3dKXs9fOmzGrGAvyjB2typdaMPWFvSSmSNAf+9jH6gaqpGVQwhXFi3xgm+oPhuPHj6/L3j//+c86QQPbzVZLXHUMSvrBMl7cPheh+m9/+1tJKVbDtspyZZ3dNvuIsbzst6fNNtusJXAnnnhiXRyWUYweeOCBunr+yEc+MqD9+wP96aefblFbXu+6667aWWed1TIc5iuzZ8+unXDCCbWTTjqp9uc//7kSrSz9HGHTx5Nl/CQiX0xIy6I5c+aEJgPHHntsjcnGz372s8KSytIOn3nmmbp2QV3moX/96181NPKG9a9//etU7Pw4X/nKV1LFsUDYGTKeFvWXZ1LWD30aEyD6c6M+YrhzVZ/20Sjuvhv6ct7SVDUW5Bk7WpUxtWCKBN6KGFxtYOSaNAPGjrQoXuRnjz32qEvzBz/4wYBsfuITn6gL86Mf/WhAmCofsKlpxRVXrMsTeF1xxRVVZqPwtKouF8u+fnsrwsb03//+dx3PpDYcB+573/teXZxZs2bVBYm/Zwn0+eefD8P8/ve/r2FfSDnYBNVq0kakv/zlLzWEWSY2rfrSHXfcUWMDFwLiAQccUPvUpz4V2oBvvPHGpU/QMKnx6wcNW6cIgZRJgeVn5syZhWUlSzu8/PLLo7ystNJKLScKTz75ZO29731vbZ111klcUaE8VjY+GC+99FLL8mFqZXG4MjnoVWommHZ7n164cGGN79aqq64a1oe/4pJUH+rTSagU86wb+nLeklQxFpQ9dqQWTN///venwouPoA12yy+//IA4Bx54YK0oXtjiMahbelzRUMWJzRV+mFYf83j8In8jeKCF8/PD/Ve/+tUik6mcVyfK9cgjj9ThmNZ+uNlHDOBWWGGFiC82n3w4mhGaSL8+mbHGCU8QJoASFnOBYcOGhfHYMHXMMcfUbZqKx+f3vHnzap/97GfrNls1a8vY3FIWBNhbb701YnnyySeH6RI3jQATRWzz5jvf+U4dLu0uFaPZhUdRZJshqVOElaIoSzsEC2szacxFMAGx8LSlOH3rW9+K3m+66abx14m/n3rqqSgOvKdNm5YYrhce9mqfZsUPO11/WbSZYKo+XW5r7Ia+nLeEVYwFZY8dqQVT3Jqk0eZsueWW0WA3evToARizUaQoXnfeeWeUFgPrGmusMSA9HrBEstxyy9WFTbtRJpFhjof+hgf70GACgYlDL1MnyvW///0v1P4Zjmk1zq0+YnENOx+PZnTwwQdHbatRG0RoRjBlYxWCIYMH9qe/+MUvQg8SzfjbO3bDoo1lGe8DH/hAmGYzwZQNMGCzzz77GIvwCm5MDnmHJrEsIm9WN1wRDNPSs88+W1tzzTVbaoTT8iPcoYceGuYnaVxqh088bJZ2iDBi2LDxqxWxGoHGGXda8+fPHxD8ox/9aMSPtpWG4mYr7EjuVerVPs24wOqJr0FvJpiqT5fbQruhL+ctYRVjQdljR2rBlEGUZcdWxOBpA26SFhCXTUXxitv27bvvvg2zd9BBB0X5Iv20mwMaMszwwnfhYBixpN+L/kr94neyXL5HhqOPPtrPVsP7Vh+xX/7yl3VtpdUOZxPyqNPDDjtsQLp8aJiMjR07tjD/ueb7tJFgevfdd0dloDxxsokELqzSTDjj8Vv9xobS30TGPS5G0hATyTFjxoT5P+ecc9JESRXGNhu1u/s8DfN22mHc3Q6/IZbSs5gY4A8VLbCNKWlXDv7whz9EcXDjlaTpDzPWA/96vU/jNs3qr5Fgqj5dTUPsZF+mhIx/WWWCqsaCsseOtgRTbLSaubLBxs86F1c2wvgEaLbhJy8v+LJk5af305/+1E+u7t7v+MRhibNKJ/vsGMd5up9fPta9rKUA4E6Xy3dgv80229TVeaMfLIdbPaCVixOCGrZ8Foar74LMD+/bUTLpoo379OCDD4abQ+CBfXVR1Eow/drXvhblH5++cfJdVPmbteLhsv6mL/r4ffCDH0zFiqW0jTbaKIz79re/vfbiiy+mitcqEFpGG3uy+Fhuxb+dduh7YUBzCl1yySVhO8myknPddddFWLMkHG+DjfLu+xveYYcdGgXriee93qf971MjwVR9upqm2Mm+jLLNNowyNrQ7HlQ1FpQ5doQKx7RVbR8ZNC1JGhYGfv80CGawccJmzPhwzcPL78jGs1GHtnzwcbSwXHENUQXhMJ3NJn7a3KfV8FWRxyxpdEO50MzZ6TkslacxiTDNGXWAfWdSe2ZA8B3ns3SOazKfWIIzIQpeLMvH6cgjj6yrd2yi0cDizxNH/nvuuWe4KYnlV+InLdPGefK7lWA6atSoMN1VVlklKXrN19ph31o0mcbT2jwmDo0EJiaI2FDiUcPqknhg04r+/ve/hxtHWKbH1peZfBKxCx+e2L1jCoFnDjZBMgYxDjTSFrJkhSurr3/962E6N954Y406jVM77fALX/hC1CY4gIEJPOVmk1ycyCvaE7SgfAxIP06+S7I0XiQsPsKo1U9W137Gq9PXXu/T/ves0Xes3/t0s77G5j+UOHgtwQVknDCHYcXBFGf0G5QCmEhw4I75VGfco78z3jIJTqJO9WXGJCbj1ie5vutd74rKlJTX+LOqxoKyxo4Ig3jBGv32waKDoBG57777Qjc0zDAY8C0MH/EkTUd8ExLh2+UFX5Ym+dBZenbdeuutw5NkGmlCadQWlivLVzTotEuMjbBp9dy3QbT0cahN5+ll6pZy+dqSRjvzGbAYlFhqtzqwK+VgkIrXB7NPfzPfWmutFQozaNzYKIJ7G3ggbM6YMSOxKuOaQ0uz2ZWBtBW1EkzNhnHddddNZIVQbXnYb7/9EsO0+xD8ONJv4sSJEW9LgyseDvAKwKEDeBXAlpfJLP3QD2f3N9xwQ8MsYIO6yy67hEvYeDlAG019oBVFcLQPlDHAzha+tBWOMI7bnLPjPf6hQghEYGSs4WOIA3h4oBlPojTtkHi+7S0aTlZO8C4SzzN4Yi/G5Mkwia9Cwc8EFsKkaTvEwZbune98Z8iXlRzzEsG7XqF+6tNpBNN+7tPN+hrus9i7YppE2rk/QfPHMmSMSy+9tE6pwLgNMTn1lQ0oFRpRJ/py3G7T+jwKvbRUxVhQ5tgRYZC2wGgWOE2AAdwAi1/5wKAFbLSLGc0Mrkzy8PKPi4ynb7+TbFspJ9oof/ejhUewKIt8o3ZLjw80s7Jepm4q10033RS1SYSVJEIjZPg3usbPuYcP9cQOavso+HGZzSJkMZtvRHRiNhj58dLco/1rRq0EU3xLkg4mM0mE0GNL2wiIRZD/cU1TxmZhMLGIC2qWR9IxW/ULLrjAHof+H40nGyN9wpaWd9ihI4idffbZYb35rrzQjPpkWl/f+T0TyiT/zMRL0w4J5/utpQ6YGDTT9JvWnfHVND/wgZiE+1rmRhOzN0O/9Z+NgoYVmthepH7q037faaQx7ec+naavMRk1jynxE/nskBOWgSEUWOaxB3kAZRbjNWZLeAai7fO7EXWiL5OX7bbbLuqX1j8xu0pDVY0FZY8dIQZpChwPg6YULSmaJpai+PBS4fGlzng8Gkd8h3NWXnHe3fgbZ/nxI0cRjG+77bZuzG7qPHVjuWywoTP7rpFSFypFQFxkMEjgvLxVWzd2LE0gzLA5BWEebR1nkiMQoYVDw8XAycDK8ogNvCxNNqNmgilCjg1qzexubcMMWrleIXDng0L5kkxhTCPi++RkjDE86H/nn39+VFyEX8N83Lhx0XNuzPOB79UA+9BmBwGkaYe4qqItYFdK3lqRTcbxaBInNKhWNoTvRsK8H48PmK04MeFis8VgpW7p060E037u07S9tH2NSSLtncll3AQLhY9t/MMNno1vLN1j0sRqDjRp0qSQR7OxkXBV92XShDAjwye8bapFi5iGqhgLqho7Um9+SgOMwryFABXo7+6zj8eZZ575VqAevOvWcvFxNRti3G91AyGUspRO3fuCUrO8gS8DJpq0RisPxG8mmPLeNto18mlZhsaUdMumLbbYIsQTH7BxgQq87DSSQw45JMqKrxVlAhAnq6P4yUdoV6k76oIJBcSyYbMPRdHtEI27LWEmaTYnT54c5pF88sFNQ1OmTAnjoIFrZ5kwDe9+D1NWn24lmIJrv/Zpypa2r3GSo31LcYZvhB02E1YTVvH1SzhsNjEX8n3/2qbpo446yqInXqvuy34mmDAxiUZgTzPZJG4VY0FVY4cEU781FHjP8px1ILtit9Lr1M3lYqZpxuP+8munMMf2kbpHAGgmZMbzRztB89eMWgmmthER7wJJ5Ntl+RrBpLDd8sw/lSVpo5DvTsf3aWuG+thqIlj4xIYsWyL14xCGlQHTunAqD0uJaajIduh7Okk6thS/pja+pHHnx2oN5UXYjq9epSnbYA9TVp9OI5j2Y5+29tROXzOFj3/4BitP3/72t41dnbbTP9iD/m+T1zRmL1X2Zcs8ecSmln6d9mhh4pY9FlQ5dkgwtdZQ4JUB3z4WdmVQabTzt8CkS2XVC+ViidQEDewIO0XYM1vdM/NNu7kE20g0ZJjJNKNWgim7/km/kf303Llzo/yl3TDTLD9VvGOXPmWiftEsx4klesPcvBvgOcI2D7E0Fyc7GxtBLclWGK2K8Wxkvxznye+i2iHmClaPaETZ0Y/WCPKPgG5k9xoG/P9/mLiw4QnzIl+D5IfRfWMEyuzTaQTTfuzTPtpp+xqrIfQJhDcId3gs4/uTTtsTwNjrn26HOSFx0T6zGpGGiu7LmBU00oLSt82sgYNYkghzpqrHgirHjlArnlRwPcuOADMs/+hJOgG/zb4lO+fOxuylcnGWNDZ0CBu+PWGVCOKqhLq3P3aMk69GhD9YbE4RuthgmCQk+XFbCab+8nWSs2ZcU1ne0mja/LQ7dW+akg996EMDssBAbwO6fbAIZJsYaAvs7vUJbalt5sQdlJG/EYmlQd/VW9KRxxYvfi2iHZonAN+vrpkpsFmPOmRXf5I21c8P/nYRSEeMGJHKrtWPq/s3ESizT6cRTPuxT2fpa2ZLyYQTkyS8gPimUriQs7Et7jv6gAMOCN+1M8mk9ovsy83c4CGM4gkGzwKNyPeGU8VYUPXYER6Z3ajwet4+AnzozGDZOgZXNrv0MvViudhNj60mH+0yHKq3qk8EJXMn5bcFZvJsZmEDDUvoLDObQEU4hKAkQTKeXivBFB4IY/Bkk02c7CQ0BPhGs/d4nE7/Nhs7NlzGyTQhCF8I+UYs74EBKxZxMsEO/H1TC+rA39zGR8nqcNq0aXE2TX/naYfswEfjQ9p2FCX3TGDw12rv/CXMpMxgr8ZHnE14aTX3SXwG+7My+3QawbQf+3SWvsbKoy3HM8HGQ4fvrQL/xPQTvFjE3VaaJnXq1KmhUNvORLOovszKI9pa+mWcsGtttbK68847R+NR2WNBJ8YOMNBSfrxl5PhtS432EePq27fkYN3RqL1aLjo/GlMaeieIJQlm5n57aHSPENvOCUzmC7OZLz52mZMewq9P7Ao3d0tlHM/pp1XkvS1lxm1i8UNs9ndxu0nTdmIb6JNpXTgSeM6cOf6rGs/ip8iZd412bL6MadZ2iCbb2gs+GJlk8RtTD1xE8XE+66yzLJmGVwQq+xg2DKQXqRAoq0/jIsrqupFPZDLYb306a1/zJ/1xky3zYoHPZJ8effTRCGMmr7i/RIPaDhXRl5n0ohhIOqgjTV4Yu8z2veyxoFNjhwTTNC0hRRhmYDaw2BWhwdfEpGDTdUH6tVxVAo19Di6h8M/GjB07P9xB4QMV376cj96O1pIladOGIpzwsUwitAWmwceWEmLpy5az8KfJ714hW8rERsv6FWVkly0mEPFlO+xQTZjjI2R07rnnhuFxRv3444/b4/CK0E7/RYvKPWQn1rHr17djC1+W+M/sX9H8MMGlLumPnFrDSVnxvJeYFbGOIVBkn6bvc8CCfTeaKTP6qU/n6WusOIEXE0+ERSNMb+zgiPjqhgn/jJ145xg5cmSdptV4lHGN92XMkmwMy5IeYz4+nPt1LJBgmqVVxOJg32UbbmxwYdnRX1KMRcn1k8HJ/9DmYtYkcifK1SQ7g/4VyyrYHCa1NbSJ7GyNE23QtKtoGcz2lp3/ndIkx/OY9jf2aDjQxmYbTcv2228fnuDExp9G/kBZ5mYZm5OeMH9AU8wOe5bAzLWMnz624Ggj4Y89JloX0sPGM80uXp9X3vvZs2dHwgof0Xnz5uVlqfhdhgCHOnAykX037EpfbaRR65c+naevsdJDP43v3TCzGyakbIryiZPd7DAKFAO+uY4frox79eX2UB1C8KAziDIiENiDuMD+ywUG13UcgpNrXLDsUvcsz49gduWCndQu8KHoAk0JE4rwLw/PZnE7Va5Ae9AsW3qXEYFg81rYdoKB2QX+Xt3w4cMzcup8tOCD4wKtsQuW110gsLnA1rJppgIBPOw7wa5qFwimLhA4XaBpToxDWHgHqx1hPws2vLhAkHVrr712wziJjAp6SH8PTBVcoBlywfJdQVzFph8Q6PU+nbWv8e0bPXq0C1agXOAyqq2qDDSNLlDsuGBFpPL+rL6cvqokmKbHKjHkrrvu6qZPn574ruyHZc4pOlWuMstUdn2IvxAQAkJACJSLQLDa4YKd4i4wgapcuCy3ZOJuCEgwNSQyXIONBy6wHcwQs5goZQlxnSxXWWUqBnFxEQJCQAgIgSoRCOy6XXBwigtMmFzgf9cF5jlu1qxZLjCXqzIbSqtCBCSYZgQ7sBkJlxMCu7eMHPJHK0OI63S5yihTfqTFQQgIASEgBKpGgKX3wL7bBSfhubFjx7rAztsFm5hC05qq86L0qkNAgmkGrANfgKENWuB+IkPs4qIULcR1Q7mKLlNxaIuTEBACQkAIVIlAsEEpFEwDP6Uu2PjpAr+lbtiwYVVmQWl1AAEJpm2CjuA0fvx4F/gzbDNm8cGLFOK6pVxFlql4xMVRCAgBISAEqkQg2E3vAtdS4UbHKtNVWp1DQIJpm9gH7j1c4DuszVjlBC9SiOuWchVZpnJQF1chIASEgBAQAkKgLAQkmLaJbOA/zbGs0A0UnClbWDa6pVxFlqkwcMRICAgBISAEhIAQqAQBCaaVwKxEhIAQEAJCQAgIASEgBFohIMG0FUJ6LwSEgBAQAkJACAgBIVAJAhJMK4FZiQgBISAEhIAQEAJCQAi0QkCCaSuE9F4ICAEhIASEgBAQAkKgEgQkmFYCsxIRAkJACAgBISAEhIAQaIWABNNWCHXh+zlz5oRnBU+aNKkLc5ctSwsXLnTnnXeeGzNmjFt33XWzMVEsISAEhIAQEAJCoKcRkGBaUPU99dRTbtSoUW7o0KHurrvucksuuWRBnN9ic88997ijjz46PC944403dhwfWjbtvffe7sYbb3QXXHCB23777QtP7tVXXw0F0pNOOslxygfHze2www6FpyOGQkAICAEhIASEQPcjIMG0oDq6/PLL3W677RZyu/LKK91OO+1UEGcXCmyTJ092M2fOdAsWLAj5brLJJqULpo899phbY401wvTGjRvnZsyYUViZYHTxxRe7Y445xpHOa6+9FvKWYFooxGImBISAEBACQqCnEJBgWlB1Pf7446HGdKmllnI33XSTW3nllQvi7BxHss2fP98hjK6//vruvvvuC++r0JjuueeeoUB89tlnuz322KOwMsHo2muvdZtuuqm7/vrr3YQJE0LeEkwLhVjMhIAQEAJCQAj0FAISTHuqupwbO3ZsuNxdhca0KmjuuOOOUNAmPQmmVaGudISAEBACQkAIdB8CEky7r06a5kiCaVN49FIICAEhIASEgBDoYQQkmOasvDfeeCO0AWV5/e6773ZsStpmm21ycm0cvSrB9LnnnnMPPvhgaDYwb948d9RRRzXOVM430pjmBFDRhYAQEAJCQAj0CQISTHNUJELp5ptv7tgt//LLL4ecrr76arfddtu5F154wbFTPwstvfTSbvXVV0+MWoVgeuihh7pLLrnEPfnkk2EeRo8e7W6++ebwHmH19ddfT8xbq4drrbWWW2KJJQYEk2A6ABI9EAJCQAgIASEwKBGQYFpAtU+ZMsUdd9xxbvHFF3fPP/986CqKzUITJ07MxH3FFVd0zzzzTGLcKgRTEmaXPBu5uFI+ds9DK620knv22WfD+3b//fCHP3QHH3zwgGgSTAdAogdCQAgIASEgBAYlAhJMC6j2nXfe2V111VVu2223dTfccEPI8bLLLnMnnnhiJu7LLLNM6EA/KXJVgqnvKgovA1tttVWYnS233DLUBiflrdWzww47LHFnvwTTVsjpvRAQAkJACAiBwYGABNOc9cyy9gorrOBeeukld8oppziWwcukqgRTBOvdd989XHp/8cUXQ21wWeWSYFoWsuIrBISAEBACQqC3EJBgmrO+brvtNrfZZpuFXDjxacMNN8zJsXn0qgTTAw880J177rlu6623dr/97W+bZyrnWwmmOQFUdCEgBISAEBACfYKABNOcFYntJceErrLKKuFmoSFDhoQcsc3kuM0stOiiiyZuEoJXVYLpeuut5+6//3537LHHuiOPPDIqBpu8arVa9LudG45pXWSRRQZEkWA6ABI9EAJCQAgIASEwKBGQYJqz2keNGuVuvfVWxwlJ7GQ3yrP5adiwYaELKuPlX6sQTJ9++ulQ0CbdW265JTzRyvKQZ/PT1KlT3f7772+soqsE0wgK3QgBISAEhIAQGNQISDDNUf2vvPKKW3bZZcOd6xdddFEonLJjHe3p+eefn9nedLnllnNz585NzFkVgun06dPdrrvuGnoXwO3VYostFuVlxIgRoeeB6EEbNyeffLIE0zbwUlAhIASEgBAQAoMNAQmmOWoc357sUoeeeOIJd/zxx4cboVj+Los4W/722293G220kbvzzjtLSWbSpEkO105jxoxx1113XSlp+Exnzpzpdtxxx/DRjBkz3Lhx4/zXuhcCQkAICAEhIAQGCQISTHNUNFrSffbZJ9yx/vnPf97NmjXLzZ49u7Qd7A888IDD9hMbz6FDh7qHH364oSP+HMVyI0eOdPfee6874YQT3OGHH56HVcu4lGWvvfZyl156aRiWpX6W/EVCQAgIASEgBITA4ENAgmmOOkdziQYTQpi75ppr3KqrrpqDY3LU+fPnu1122SXUkC5cuDAKhBkBR6Biz7rOOutEz/PccBQpdqQIjJQP/mXRqaee6s4555xQwPbT2GCDDdz48ePDjVf+c90LASEgBISAEBAC/Y2ABNOc9fvQQw85hDkEOHbT9zpdccUVbsKECW748OHukUce6fXiKP9CQAgIASEgBIRADyEgwbSHKquKrO63337uwgsvdJMnT3ZnnHFGFUkqDSEgBISAEBACQkAIhAhIMB3kDYEjVPEggK9SdtyvttpqbsGCBQ571jXXXHOQo6PiCwEhIASEgBAQAlUiIMG0SrS7MK0tttgi9FWKR4Hll1/eHXTQQe6II44IPQx0YXaVJSEgBISAEBACQqCPEZBg2seVm6Zo06ZNc/vuu2+4259d/hMnTnT4G006oSkNP4URAkJACAgBISAEhEBWBCSYZkVO8YSAEBACQkAICAEhIAQKRUCCaaFwipkQEAJCQAgIASEgBIRAVgQkmGZFTvGEgBAQAkJACAgBISAECkVAgmmhcIqZEBACQkAICAEhIASEQFYEJJhmRU7xhIAQEAJCQAgIASEgBApFQIJpoXCKmRAQAkJACAgBISAEhEBWBCSYZkVO8YSAEBACQkAICAEhIAQKRUCCaaFwipkQEAJCQAgIASEgBIRAVgQkmGZFTvGEgBAQAkJACAgBISAECkXg/wDnMd/hIIvOXwAAAABJRU5ErkJggg==" - } - }, - "cell_type": "markdown", - "id": "d771c52c", - "metadata": {}, - "source": [ - "*Here i is the number of time windows and j is the number of station pairs.* from Hirose et al. (2022) --> ![image-3.png](attachment:image-3.png)" + "An energy density with the spatially homogeneous intrinsic absorption is obtained by multiplying $e^{-bt}$ with the right-hand-side of Equation (1), where ***b is intrinsic absorption parameter***\n", + "\n", + "\n", + "The sum of squared residuals (SSR) between *Eobs and Esyn* was calculated, and the optimal values of the parameters $l$ and $b$ that minimize SSR through the grid search were determined. *Here i is the number of time windows and j is the number of station pairs.* from Hirose et al. (2022)
\n", + "\\begin{gather*}\n", + "SSR=\\sum\\limits_{i=1}^{6} \\sum\\limits_{j=1}^{N} (log_{10}E_{obs}(t_i,j) - log_{10}E_{syn}(t_i,j))\n", + "\\end{gather*}" ] }, { @@ -500,18 +488,20 @@ "outputs": [], "source": [ "def plot_fitting_result(mean_free,intrinsic_b,tt,Eobs,Esyn,fname,dist,twin,fmin,fmax):\n", - " plt.figure(figsize=(6,2))\n", + " plt.figure(figsize=(4,2))\n", " plt.yscale('log', base=10)\n", - "\n", + " \n", + " intrinsic_Q=(2.0*np.pi*((fmin+fmax)/2))/intrinsic_b\n", + " \n", " pymax=np.max(Eobs[:-2]*5)\n", - " pymin=10**(-6)\n", + " pymin=10**(-4)\n", " plt.ylim( pymin , pymax )\n", " plt.plot( tt, Eobs, \"k-\", linewidth=1)\n", " plt.plot( tt, Esyn, \"b--\", linewidth=1)\n", " plt.plot([twin[0],twin[0],twin[-1],twin[-1],twin[0]],[pymin, pymax,pymax,pymin,pymin],\"r\", linewidth=2)\n", "\n", - " plt.title(\"%s %.2fkm @%4.2f-%4.2f Hz, intrinsic b: %.2f\"\n", - " % ( fname,dist,fmin,fmax,intrinsic_b))\n", + " plt.title(\"%s %.1fkm @%4.2f-%4.2f Hz, \\nintrinsic b: %.2f, intrinsic Q: %.2f\"\n", + " % ( fname,dist,fmin,fmax,intrinsic_b, intrinsic_Q))\n", " plt.xlabel(\"Time [s]\")\n", " plt.ylabel(\"Energy density Amp\")\n", " plt.tight_layout() \n", @@ -547,13 +537,21 @@ " plot_fitting_result(result_mfp[fb],result_intb[fb],fmsv_mean[aa][0][:],\n", " Eobs[fb],Esyn[fb],target_pair,vdist[aa],twinbe[aa][fb],fmin,fmax)\n" ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "23bd56f2", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { "kernelspec": { - "display_name": "noisepy", + "display_name": "noisepy_test", "language": "python", - "name": "noisepy" + "name": "python3" }, "language_info": { "codemirror_mode": { From c1c6ad5e94550ffa36569a45520ee318ae43519a Mon Sep 17 00:00:00 2001 From: Kuan-Fu Feng Date: Wed, 13 Mar 2024 10:38:22 -0700 Subject: [PATCH 3/3] Correcting the equation description and missing sign in the notebook markdown --- .../monitoring/attenuation_singlestation.ipynb | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/tutorials/monitoring/attenuation_singlestation.ipynb b/tutorials/monitoring/attenuation_singlestation.ipynb index 1306b4ca..338dccf7 100644 --- a/tutorials/monitoring/attenuation_singlestation.ipynb +++ b/tutorials/monitoring/attenuation_singlestation.ipynb @@ -381,23 +381,25 @@ "metadata": {}, "source": [ "### Step 2 --- Calculation of synthesized energy densities (***Esyn***) via a grid search \n", - "The 2-D radiative transfer equation for scalar waves ***(Shang and Gao 1988; Sato 1993)*** assuming isotropic scattering and source radiation in infinite medium to calculate ***synthesized energy densities Esyn*** :\n", + "The 2-D radiative transfer equation for scalar waves ***(Shang and Gao 1988; Sato 1993)*** is assuming isotropic scattering and source radiation in infinite medium to calculate ***synthesized energy densities Esyn*** :\n", "\\begin{gather*}\n", - "E_{syn}(r,t)=\\frac{e^{-ctl^{-1}}}{2{\\pi}cr}\\,\\delta(t-\\frac{r}{c}) + \\frac{e^{l^{-1(\\sqrt{(c^2t^2-r^2)}-ct)}}}{2\\pi l \\sqrt{(c^2t^2-r^2)}} H(t-\\frac{r}{c})\n", + "E_{syn}(r,t)=\\frac{e^{-ctl^{-1}}}{2{\\pi}cr}\\,\\delta(t-\\frac{r}{c}) + \\frac{e^{l^{-1(\\sqrt{c^2t^2-r^2}-ct)}}}{2\\pi l \\sqrt{(c^2t^2-r^2)}} H(t-\\frac{r}{c})\n", "\\end{gather*}\n", "\n", - "* $l$ scattering mean free paths \n", + "* $b$ is the intrinsic absorption parameter \n", + "* $l$ is the scattering mean free path\n", "* $c$ is the Rayleigh wave velocity\n", "* $r$ is the distance between the source and receiver\n", "* $H$ is the Heaviside function\n", "\n", - "An energy density with the spatially homogeneous intrinsic absorption is obtained by multiplying $e^{-bt}$ with the right-hand-side of Equation (1), where ***b is intrinsic absorption parameter***\n", + "An energy density with the spatially homogeneous intrinsic absorption is obtained by multiplying $e^{-bt}$ with the right-hand-side of Equation (1), where ***b is intrinsic absorption parameter***.\n", "\n", "\n", - "The sum of squared residuals (SSR) between *Eobs and Esyn* was calculated, and the optimal values of the parameters $l$ and $b$ that minimize SSR through the grid search were determined. *Here i is the number of time windows and j is the number of station pairs.* from Hirose et al. (2022)
\n", + "The optimal value of the parameters $l$ and $b$ is determined by minimizing the sum of squared residuals (SSR) between *Eobs* and *Esyn* through the grid search. (please refer to Hirose et al. 2019, 2022)
\n", "\\begin{gather*}\n", - "SSR=\\sum\\limits_{i=1}^{6} \\sum\\limits_{j=1}^{N} (log_{10}E_{obs}(t_i,j) - log_{10}E_{syn}(t_i,j))\n", - "\\end{gather*}" + "SSR=\\sum\\limits_{i=1}^{6} \\sum\\limits_{j=1}^{N} (log_{10}E_{obs}(t_i,j) - log_{10}E_{syn}(t_i,j))^2\n", + "\\end{gather*}\n", + "*Here i is the number of time windows and j is the number of station pairs*." ] }, {