From 83277301f60313bd9512f9e827257e2ebb4d562a Mon Sep 17 00:00:00 2001 From: "Chimitt, Nicholas M" Date: Tue, 27 Jul 2021 15:00:22 -0400 Subject: [PATCH] Add files via upload --- Turbulence_Sim_v1_python/Example_1_full.py | 47 ++ .../Example_2_tiltonly.py | 39 ++ .../Integrals_Spatial_Corr.py | 94 ++++ Turbulence_Sim_v1_python/Motion_Compensate.py | 39 ++ Turbulence_Sim_v1_python/TurbSim_v1_main.py | 510 ++++++++++++++++++ Turbulence_Sim_v1_python/data/cameraman.png | Bin 0 -> 38783 bytes Turbulence_Sim_v1_python/data/chart.png | Bin 0 -> 974 bytes Turbulence_Sim_v1_python/main.py | 25 + 8 files changed, 754 insertions(+) create mode 100644 Turbulence_Sim_v1_python/Example_1_full.py create mode 100644 Turbulence_Sim_v1_python/Example_2_tiltonly.py create mode 100644 Turbulence_Sim_v1_python/Integrals_Spatial_Corr.py create mode 100644 Turbulence_Sim_v1_python/Motion_Compensate.py create mode 100644 Turbulence_Sim_v1_python/TurbSim_v1_main.py create mode 100644 Turbulence_Sim_v1_python/data/cameraman.png create mode 100644 Turbulence_Sim_v1_python/data/chart.png create mode 100644 Turbulence_Sim_v1_python/main.py diff --git a/Turbulence_Sim_v1_python/Example_1_full.py b/Turbulence_Sim_v1_python/Example_1_full.py new file mode 100644 index 0000000..c5d53e0 --- /dev/null +++ b/Turbulence_Sim_v1_python/Example_1_full.py @@ -0,0 +1,47 @@ +''' +Demo of Multi-Aperture Turbulence Simulation. + +Nicholas Chimitt and Stanley H. Chan "Simulating anisoplanatic turbulence +by sampling intermodal and spatially correlated Zernike coefficients," Optical Engineering 59(8), Aug. 2020 + +ArXiv: https://arxiv.org/abs/2004.11210 + +Nicholas Chimitt and Stanley Chan +Copyright 2021 +Purdue University, West Lafayette, In, USA. +''' + +from matplotlib import pyplot as plt +import numpy as np +import TurbSim_v1_main as util + + +def rgb2gray(rgb): + return np.dot(rgb[...,:3], [0.2989, 0.5870, 0.1140]) + + +img = rgb2gray(plt.imread('chart.png')) +N = 225 # size of the image -- assumed to be square (pixels) +D = 0.2 # length of aperture diameter (meters) +L = 7000 # length of propagation (meters) +r0 = 0.05 # the Fried parameter r0. The value of D/r0 is critically important! (See associated paper) +wvl = 0.525e-6 # the mean wavelength -- typically somewhere suitably in the middle of the spectrum will be sufficient +obj_size = 2.06*1 # the size of the object in the object plane (meters). Can be different the Nyquist sampling, scaling + # will be done automatically. + +param_obj = util.p_obj(N, D, L, r0, wvl, obj_size) # generating the parameter object, some other things are computed within this + # function, see the def for details +print(param_obj['N'] * param_obj['delta0'], param_obj['smax'], param_obj['scaling']) +#print('hi') +#print(param_obj['spacing']) +#print('hi') +S = util.gen_PSD(param_obj) # finding the PSD, see the def for details +param_obj['S'] = S # appending the PSD to the parameter object for convenience + +for i in range(100): + img_tilt, _ = util.genTiltImg(img, param_obj) # generating the tilt-only image + img_blur = util.genBlurImage(param_obj, img_tilt) + plt.imshow(img_tilt,cmap='gray',vmin=0,vmax=1) + plt.show() + plt.imshow(img_blur, cmap='gray', vmin=0, vmax=1) + plt.show() diff --git a/Turbulence_Sim_v1_python/Example_2_tiltonly.py b/Turbulence_Sim_v1_python/Example_2_tiltonly.py new file mode 100644 index 0000000..4b1b805 --- /dev/null +++ b/Turbulence_Sim_v1_python/Example_2_tiltonly.py @@ -0,0 +1,39 @@ +''' +Demo of Multi-Aperture Turbulence Simulation (Tilt-only). + +Nicholas Chimitt and Stanley H. Chan "Simulating anisoplanatic turbulence +by sampling intermodal and spatially correlated Zernike coefficients," Optical Engineering 59(8), Aug. 2020 + +ArXiv: https://arxiv.org/abs/2004.11210 + +Nicholas Chimitt and Stanley Chan +Copyright 2021 +Purdue University, West Lafayette, In, USA. +''' + +from matplotlib import pyplot as plt +import numpy as np +import TurbSim_v1_main as util + + +def rgb2gray(rgb): + return np.dot(rgb[...,:3], [0.2989, 0.5870, 0.1140]) + + +img = rgb2gray(plt.imread('chart.png')) +N = 225 # size of the image -- assumed to be square (pixels) +D = 0.1 # length of aperture diameter (meters) +L = 7000 # length of propagation (meters) +r0 = D/2 # the Fried parameter r0. Most importantly, the denominator is the desired D/r0 ratio +wvl = 0.5e-6 # the mean wavelength -- typically somewhere suitably in the middle of the spectrum will be sufficient + +param_obj = util.p_obj(N, D, L, r0, wvl) # generating the parameter object, some other things are computed within this + # function, see the def for details +S = util.gen_PSD(param_obj) # finding the PSD, see the def for details +param_obj['S'] = S # appending the PSD to the parameter object for convenience + +for i in range(100): + img_, _ = util.genTiltImg(img, param_obj) # generating the tilt-only image + + plt.imshow(img_,cmap='gray',vmin=0,vmax=1) + plt.show() \ No newline at end of file diff --git a/Turbulence_Sim_v1_python/Integrals_Spatial_Corr.py b/Turbulence_Sim_v1_python/Integrals_Spatial_Corr.py new file mode 100644 index 0000000..c7d986f --- /dev/null +++ b/Turbulence_Sim_v1_python/Integrals_Spatial_Corr.py @@ -0,0 +1,94 @@ +''' +Spatial Correlation Functions for Tilts + +Zhiyuan Mao and Stanley Chan +Copyright 2021 +Purdue University, West Lafayette, In, USA. +''' + +import numpy as np +import scipy.integrate as integrate +from scipy.special import jv +import os +from math import gamma + + +def genTiltPSD(s_max, spacing, N, **kwargs): + D = 0.1 + D_r0 = 20 + wavelength = kwargs.get('wavelength', 500 * (10 ** (-9))) + L = 7000 + delta0 = L * wavelength / (2 * D) + delta0 *= 1 # Adjusting factor + # Constants in [Channan, 1992] + c1 = 2 * ((24 / 5) * gamma(6 / 5)) ** (5 / 6) + c2 = 4 * c1 / np.pi * (gamma(11 / 6)) ** 2 + s_arr = np.arange(0, s_max, spacing) + I0_arr = np.float32(s_arr * 0) + I2_arr = np.float32(s_arr * 0) + for i in range(len(s_arr)): + I0_arr[i] = I0(s_arr[i]) + I2_arr[i] = I2(s_arr[i]) + i, j = np.int32(N / 2), np.int32(N / 2) + [x, y] = np.meshgrid(np.arange(1, N + 0.01, 1), np.arange(1, N + 0.01, 1)) + s = np.sqrt((x - i) ** 2 + (y - j) ** 2) + s *= spacing + C0 = (In_m(s, spacing*100, I0_arr) + In_m(s, spacing*100, I2_arr)) / I0(0) + C0[i, j] = 1 + C0_scaled = C0 * I0(0) * c2 * ((D_r0) ** (5 / 3)) / (2 ** (5 / 3)) * ( + (2 * wavelength / (np.pi * D)) ** 2) * 2 * np.pi + Cfft = np.fft.fft2(C0_scaled) + S_half = np.sqrt(Cfft) + S_half_max = np.max(np.max(np.abs(S_half))) + S_half[np.abs(S_half) < 0.0001 * S_half_max] = 0 + + return S_half + + +def I0(s): + # z = np.linspace(1e-6, 1e3, 1e5) + # f_z = (z**(-14/3))*jv(0,2*s*z)*(jv(2,z)**2) + # I0_s = np.trapz(f_z, z) + + I0_s, _ = integrate.quad( lambda z: (z**(-14/3))*jv(0,2*s*z)*(jv(2,z)**2), 1e-4, np.inf, limit = 100000) + # print('I0: ',I0_s) + return I0_s + + +def I2(s): + # z = np.linspace(1e-6, 1e3, 1e5) + # f_z = (z**(-14/3))*jv(2,2*s*z)*(jv(2,z)**2) + # I2_s = np.trapz(f_z, z) + + I2_s, _ = integrate.quad( lambda z: (z**(-14/3))*jv(2,2*s*z)*(jv(2,z)**2), 1e-4, np.inf, limit = 100000) + # print('I2: ',I2_s) + return I2_s + + +def save_In(s_max, spacing): + s_arr = np.arange(0,s_max,spacing) + I0_arr = np.float32(s_arr*0) + I2_arr = np.float32(s_arr*0) + for i in range(len(s_arr)): + I0_arr[i] = I0(s_arr[i]) + I2_arr[i] = I2(s_arr[i]) + os.makedirs('temp/In_arr', exist_ok=True) + + np.save('temp/In_arr/I0_%d_%d.npy'%(s_max,spacing*1000),I0_arr) + np.save('temp/In_arr/I2_%d_%d.npy'%(s_max,spacing*1000),I2_arr) + + +def In_m(s, spacing, In_arr): + idx = np.int32(np.floor(s.flatten()/spacing)) + M,N = np.shape(s)[0], np.shape(s)[1] + In = np.reshape(np.take(In_arr, idx), [M,N]) + + return In + + +def In_arr(s, spacing, In_arr): + idx = np.int32(np.floor(s/spacing)) + In = np.take(In_arr, idx) + + return In + diff --git a/Turbulence_Sim_v1_python/Motion_Compensate.py b/Turbulence_Sim_v1_python/Motion_Compensate.py new file mode 100644 index 0000000..0c8ba90 --- /dev/null +++ b/Turbulence_Sim_v1_python/Motion_Compensate.py @@ -0,0 +1,39 @@ +''' +Warping Function for Turbulence Simulator + +Python version +of original code by Stanley Chan + +Zhiyuan Mao and Stanley Chan +Copyright 2021 +Purdue University, West Lafayette, In, USA. +''' + +import numpy as np +from skimage.transform import resize + + +def motion_compensate(img, Mvx, Mvy, pel): + m, n = np.shape(img)[0], np.shape(img)[1] + img = resize(img, (np.int32(m/pel), np.int32(n/pel)), mode = 'reflect' ) + Blocksize = np.floor(np.shape(img)[0]/np.shape(Mvx)[0]) + m, n = np.shape(img)[0], np.shape(img)[1] + M, N = np.int32(np.ceil(m/Blocksize)*Blocksize), np.int32(np.ceil(n/Blocksize)*Blocksize) + + f = img[0:M, 0:N] + + + Mvxmap = resize(Mvy, (N,M)) + Mvymap = resize(Mvx, (N,M)) + + + xgrid, ygrid = np.meshgrid(np.arange(0,N-0.99), np.arange(0,M-0.99)) + X = np.clip(xgrid+np.round(Mvxmap/pel),0,N-1) + Y = np.clip(ygrid+np.round(Mvymap/pel),0,M-1) + + idx = np.int32(Y.flatten()*N + X.flatten()) + f_vec = f.flatten() + g = np.reshape(f_vec[idx],[N,M]) + + g = resize(g, (np.shape(g)[0]*pel,np.shape(g)[1]*pel)) + return g \ No newline at end of file diff --git a/Turbulence_Sim_v1_python/TurbSim_v1_main.py b/Turbulence_Sim_v1_python/TurbSim_v1_main.py new file mode 100644 index 0000000..9518300 --- /dev/null +++ b/Turbulence_Sim_v1_python/TurbSim_v1_main.py @@ -0,0 +1,510 @@ +''' +Collection of Optical Functions for Turbulence Simulator + +Nicholas Chimitt and Stanley Chan +Copyright 2021 +Purdue University, West Lafayette, In, USA. +''' + +import numpy as np +import math +import scipy.signal +import scipy.interpolate +from scipy.ndimage import gaussian_filter +from skimage.transform import resize +from scipy.interpolate import RectBivariateSpline +from scipy.special import jv +import scipy.integrate as integrate +from integrals_spatial_corr import save_In, In_m, I0, I2 +from math import gamma +from motion_compensate import motion_compensate + + +def p_obj(N, D, L, r0, wvl, obj_size): + """ + The parameter "object". It is really just a list of some useful parameters. + + :param N: size of pixels of one image dimension (assumed to be square image N x N). + :param D: size of aperture diameter (meters) + :param L: length of propagation (meters) + :param r0: Fried parameter (meters) + :param wvl: wavelength (meters) + :return: returns the parameter object + """ + a = {} + a['N'] = N + a['D'] = D + a['L'] = L + a['wvl'] = wvl + a['r0'] = r0 + a['Dr0'] = D/r0 + a['delta0'] = L*wvl/(2*D) + a['k'] = 2*np.pi/wvl + a['smax'] = a['delta0']/D*N + a['spacing'] = a['delta0']/D + a['ob_s'] = obj_size + a['scaling'] = obj_size / (N * a['delta0']) + + a['smax'] *= a['scaling'] + a['spacing'] *= a['scaling'] + a['ob_s'] *= a['scaling'] + a['delta0'] *= a['scaling'] + + return a + + +def gen_PSD(p_obj): + """ + This function generates the PSD necessary for the tilt values (both x and y pixel shifts). The PSD is **4 times** + the size of the image, this is to simplify the generation of the random vector using a property of Toeplitz + matrices. This is further highlighted in the genTiltImg() function, where only 1/4 of the entire grid is used + (this is because of symmetry about the origin -- hence why the PSD is quadruple the size). + + All that is required is the parameter list, p_obj. + + :param p_obj: + :return: PSD + """ + N = 2 * p_obj['N'] + smax = p_obj['delta0'] / p_obj['D'] * N + c1 = 2 * ((24 / 5) * gamma(6 / 5)) ** (5 / 6) + c2 = 4 * c1 / np.pi * (gamma(11 / 6)) ** 2 + s_arr = np.linspace(0, smax, N) + I0_arr = np.float32(s_arr * 0) + I2_arr = np.float32(s_arr * 0) + for i in range(len(s_arr)): + I0_arr[i] = I0(s_arr[i]) + I2_arr[i] = I2(s_arr[i]) + i, j = np.int32(N / 2), np.int32(N / 2) + [x, y] = np.meshgrid(np.arange(1, N + 0.01, 1), np.arange(1, N + 0.01, 1)) + s = np.sqrt((x - i) ** 2 + (y - j) ** 2) + C = (In_m(s, p_obj['delta0'] / p_obj['D'] * N , I0_arr) + In_m(s, p_obj['delta0'] / p_obj['D'] * N, I2_arr)) / I0(0) + C[round(N / 2), round(N / 2)] = 1 + C = C * I0(0) * c2 * (p_obj['Dr0']) ** (5 / 3) / (2 ** (5 / 3)) * (2 * p_obj['wvl'] / (np.pi * p_obj['D'])) ** 2 * 2 * np.pi + Cfft = np.fft.fft2(C) + S_half = np.sqrt(Cfft) + S_half_max = np.max(np.max(np.abs(S_half))) + S_half[np.abs(S_half) < 0.0001 * S_half_max] = 0 + + return S_half + + +def genTiltImg(img, p_obj): + """ + This function takes the p_obj (with the PSD!) and applies it to the image. If no PSD is found, one will be + generated. However, it is **significantly** faster to generate the PSD once and then use it to draw the values from. + This is also done automatically, because it is significantly faster. + + :param img: The input image (assumed to be N x N pixels) + :param p_obj: The parameter object -- with the PSD is preferred + :return: The output, tilted image + """ + flag_noPSD = 0 + if (p_obj.get('S') == None).any(): + S = gen_PSD(p_obj) + p_obj['S'] = S + flag_noPSD = 1 + MVx = np.real(np.fft.ifft2(p_obj['S'] * np.random.randn(2 * p_obj['N'], 2 * p_obj['N']))) * np.sqrt(2) * 2 * p_obj['N'] * (p_obj['L'] / p_obj['delta0']) + MVx = MVx[round(p_obj['N'] / 2) :2 * p_obj['N'] - round(p_obj['N'] / 2), 0: p_obj['N']] + #MVx = 1 / p_obj['scaling'] * MVx[round(p_obj['N'] / 2):2 * p_obj['N'] - round(p_obj['N'] / 2), 0: p_obj['N']] + MVy = np.real(np.fft.ifft2(p_obj['S'] * np.random.randn(2 * p_obj['N'], 2 * p_obj['N']))) * np.sqrt(2) * 2 * p_obj['N'] * (p_obj['L'] / p_obj['delta0']) + MVy = MVy[0:p_obj['N'], round(p_obj['N'] / 2): 2 * p_obj['N'] - round(p_obj['N'] / 2)] + #MVy = 1 / p_obj['scaling'] * MVy[0:p_obj['N'], round(p_obj['N'] / 2): 2 * p_obj['N'] - round(p_obj['N'] / 2)] + img_ = motion_compensate(img, MVx - np.mean(MVx), MVy - np.mean(MVy), 0.5) + #plt.quiver(MVx[::10,::10], MVy[::10,::10], scale=60) + #plt.show() + + if flag_noPSD == 1: + return img_, p_obj + else: + return img_, p_obj + + +def genBlurImage(p_obj, img): + smax = p_obj['delta0'] / p_obj['D'] * p_obj['N'] + temp = np.arange(1,101) + patchN = temp[np.argmin((smax*np.ones(100)/temp - 2)**2)] + patch_size = round(p_obj['N'] / patchN) + xtemp = np.round_(p_obj['N']/(2*patchN) + np.linspace(0, p_obj['N'] - p_obj['N']/patchN + 0.001, patchN)) + xx, yy = np.meshgrid(xtemp, xtemp) + xx_flat, yy_flat = xx.flatten(), yy.flatten() + NN = 32 # For extreme scenarios, this may need to be increased + img_patches = np.zeros((p_obj['N'], p_obj['N'], int(patchN**2))) + den = np.zeros((p_obj['N'], p_obj['N'])) + patch_indx, patch_indy = np.meshgrid(np.linspace(-patch_size, patch_size+0.001, num=2*patch_size+1), np.linspace(-patch_size, patch_size+0.001, num=2*patch_size+1)) + + for i in range(int(patchN**2)): + aa = genZernikeCoeff(36, p_obj['Dr0']) + temp, x, y, nothing, nothing2 = psfGen(NN, coeff=aa, L=p_obj['L'], D=p_obj['D'], z_i=1.2, wavelength=p_obj['wvl']) + psf = np.abs(temp) ** 2 + psf = psf / np.sum(psf.ravel()) + # focus_psf, _, _ = centroidPsf(psf, 0.95) : Depending on the size of your PSFs, you may want to use this + psf = resize(psf, (round(NN/p_obj['scaling']), round(NN/p_obj['scaling']))) + patch_mask = np.zeros((p_obj['N'], p_obj['N'])) + patch_mask[round(xx_flat[i]), round(yy_flat[i])] = 1 + patch_mask = scipy.signal.fftconvolve(patch_mask, np.exp(-patch_indx**2/patch_size**2)*np.exp(-patch_indy**2/patch_size**2)*np.ones((patch_size*2+1, patch_size*2+1)), mode='same') + den += scipy.signal.fftconvolve(patch_mask, psf, mode='same') + img_patches[:,:,i] = scipy.signal.fftconvolve(img * patch_mask, psf, mode='same') + + out_img = np.sum(img_patches, axis=2) / (den + 0.000001) + return out_img + + +def genZernCorrHighOrder_v2(): + """ + This function generates the correlation functions found in [Chimitt-Chan 2020] and [Takato et al. 1992]. + There are no inputs to this function, and is provided primarily for (i) convenience (ii) further modification, + specifically with respect to the "A" parameter below, which should be modified simply to include (n+1)*(npr+1) upon + calculation of each row. + + :return: There is no return of this function in its current form. Instead, the functions are saved as it takes a + good amount of time to generate. The functions evaluated at sufficiently long distance and resolution are provided + in the download for this simulator package [LOCATION]. + """ + lamb = 0.525e-6 + Cn2 = 1e-15 + L = 7000 + A = 0.00969 * (2 * np.pi / lamb ) ** 2 * Cn2 * L + D = 0.2034 + L0 = np.inf + k0 = 1 / L0 + + MM = 65 + NN = 150 + + svalarr = np.hstack((-np.logspace(-2,2,num=int((MM-1)/2))[::-1], 0, np.logspace(-2,2,num=int((MM-1)/2)))) + fij = np.zeros((MM,MM)) + outxx, outyy = np.meshgrid(np.linspace(-5, 5, NN), np.linspace(-5, 5, NN)) + smat = np.sqrt(outxx**2 + outyy**2) + thetamat = np.arctan2(outyy, outxx + 0.0001) + xx = np.linspace(-np.max(np.max(smat)),np.max(np.max(smat)),NN) + x_grid, y_grid = np.meshgrid(np.linspace(-1, 1, NN, endpoint=True), np.linspace(-1, 1, NN, endpoint=True)) + mask = np.sqrt(x_grid ** 2 + y_grid ** 2) <= 1 + + for i in np.arange(27,37): + #print(i) + n,m = nollToZernInd(i) + for j in range(4, 37): + print(i,j) + npr, mpr = nollToZernInd(j) + line1int = np.zeros((len(xx),len(xx))) + line2int = np.zeros((len(xx),len(xx))) + line3int = np.zeros((len(xx),len(xx))) + line4int = np.zeros((len(xx),len(xx))) + line5int = np.zeros((len(xx),len(xx))) + line6int = np.zeros((len(xx),len(xx))) + line7int = np.zeros((len(xx),len(xx))) + line1rot = np.zeros((len(xx),len(xx))) + line2rot = np.zeros((len(xx),len(xx))) + line3rot = np.zeros((len(xx),len(xx))) + line4rot = np.zeros((len(xx),len(xx))) + line5rot = np.zeros((len(xx),len(xx))) + line6rot = np.zeros((len(xx),len(xx))) + line7rot = np.zeros((len(xx),len(xx))) + for xxi in range(len(xx)): + #print(xxi/len(xx)) + x0, kap, mu, nu, ss = D * np.pi * k0, m + mpr, n + 1, npr + 1, xx[xxi] + line1int[int(NN/2 - 1), xxi], _ = integrate.quad(lambda x: jv(kap, 2 * ss * x) * jv(mu, x) * jv(nu, x) / (x * (x ** 2 + x0 ** 2) ** (11 / 6)), 1E-20, np.inf, limit=10**8) + + x0, kap, mu, nu, ss = D * np.pi * k0, abs(m - mpr), n + 1, npr + 1, xx[xxi] + line2int[int(NN/2 - 1), xxi], _ = integrate.quad( + lambda x: jv(kap, 2 * ss * x) * jv(mu, x) * jv(nu, x) / (x * (x ** 2 + x0 ** 2) ** (11 / 6)), 1E-20, + np.inf, limit=10**8) + + line3int[int(NN/2 - 1), xxi] = line1int[int(NN/2 - 1), xxi] + + line4int[int(NN/2 - 1), xxi] = line2int[int(NN/2 - 1), xxi] + + x0, kap, mu, nu, ss = D * np.pi * k0, m, n + 1, npr + 1, xx[xxi] + line5int[int(100/2 - 1), xxi], _ = integrate.quad( + lambda x: jv(kap, 2 * ss * x) * jv(mu, x) * jv(nu, x) / (x * (x ** 2 + x0 ** 2) ** (11 / 6)), 1E-20, + np.inf, limit=10**8) + + line6int[int(NN/2 - 1), xxi] = line5int[int(NN/2 - 1), xxi] + + x0, kap, mu, nu, ss = D * np.pi * k0, 0, n + 1, npr + 1, xx[xxi] + line7int[int(NN/2 - 1), xxi], _ = integrate.quad( + lambda x: jv(kap, 2 * ss * x) * jv(mu, x) * jv(nu, x) / (x * (x ** 2 + x0 ** 2) ** (11 / 6)), 1E-20, + np.inf, limit=10**8) + + base_1s = np.zeros((NN,NN)) + base_1s[int(NN/2 - 1), :] = 1 + normalize_1s = np.zeros((NN,NN)) + for anglee in np.arange(0, 360, 0.5): + base_rot = scipy.ndimage.rotate(base_1s, anglee, reshape=False) + line1rot += base_rot * scipy.ndimage.rotate(line1int, anglee, reshape=False) + line2rot += base_rot * scipy.ndimage.rotate(line2int, anglee, reshape=False) + line3rot += base_rot * scipy.ndimage.rotate(line3int, anglee, reshape=False) + line4rot += base_rot * scipy.ndimage.rotate(line4int, anglee, reshape=False) + line5rot += base_rot * scipy.ndimage.rotate(line5int, anglee, reshape=False) + line6rot += base_rot * scipy.ndimage.rotate(line6int, anglee, reshape=False) + line7rot += base_rot * scipy.ndimage.rotate(line7int, anglee, reshape=False) + normalize_1s += base_rot + + line1rot = mask*line1rot/normalize_1s + line2rot = mask * line2rot / normalize_1s + line3rot = mask * line3rot / normalize_1s + line4rot = mask * line4rot / normalize_1s + line5rot = mask * line5rot / normalize_1s + line6rot = mask * line6rot / normalize_1s + line7rot = mask * line7rot / normalize_1s + #plt.imshow(line1rot) + #plt.show() + + if m != 0 and mpr != 0 and np.mod(i, 2) == 0 and np.mod(j, 2) == 0: + fij = (-1) ** ((n + npr - m + mpr) / 2) * np.cos((m + mpr) * thetamat) * line1int + \ + (-1) ** ((n + npr + 2 * m + abs(m - mpr)) / 2) * np.cos((m - mpr) * thetamat) * line2int + if m != 0 and mpr != 0 and np.mod(i, 2) != 0 and np.mod(j, 2) != 0: + fij = -(-1) ** ((n + npr - m + mpr) / 2) * np.cos((m + mpr) * thetamat) * line1int + \ + (-1) ** ((n + npr + 2 * m + abs(m - mpr)) / 2) * np.cos((m - mpr) * thetamat) * line2int + if m != 0 and mpr != 0 and np.mod(i + j, 2) != 0: + fij = (-1) ** ((n + npr - m + mpr) / 2) * np.sin((m + mpr) * thetamat) * line3int + \ + (-1) ** ((n + npr + 2 * m + abs(m - mpr)) / 2) * np.sin((m - mpr) * thetamat) * line4int + if mpr == 0 and np.mod(i, 2) == 0 and (not m == 0): + fij = (-1) ** ((n + npr - m) / 2) * np.sqrt(2) * np.cos(m * thetamat) * line5int + if mpr == 0 and np.mod(i, 2) != 0 and (not m == 0): + fij = (-1) ** ((n + npr - m) / 2) * np.sqrt(2) * np.sin(m * thetamat) * line6int + if m == 0 and np.mod(j, 2) == 0 and (not mpr == 0): + fij = (-1) ** ((n + npr - mpr) / 2) * np.sqrt(2) * np.cos(mpr * thetamat) * line5int + if m == 0 and np.mod(j, 2) != 0 and (not mpr == 0): + fij = (-1) ** ((n + npr - mpr) / 2) * np.sqrt(2) * np.sin(mpr * thetamat) * line6int + if m == 0 and mpr == 0: + fij = (-1) ** ((n + npr) / 2) * np.sqrt(2) * line7int + + np.save('./corr_zerns_v3/Zern{}_Zern{}_v2_-5_5_150points'.format(i,j), fij) + + +def focusPsf(psf, thresh): + x = np.linspace(0, 1, psf.shape[0]) + y = np.linspace(0, 1, psf.shape[1]) + col, row = np.meshgrid(x, y) + # Don't have to normalize as psf should already be normalized + cen_row = np.uint8(np.sum(row * psf)) + cen_col = np.uint8(np.sum(col * psf)) + temp_sum = 0 + radius = 0 + psf /= np.sum(psf) + while temp_sum < thresh: + radius += 1 + return_psf = psf[cen_row-radius:cen_row+radius+1, cen_col-radius:cen_col+radius+1] + temp_sum = np.sum(return_psf) + #print(temp_sum) + #print(radius, temp_sum) + + return return_psf, cen_row, cen_col + + +def centroidPsf(psf, thresh): + x = np.linspace(0, psf.shape[0], psf.shape[0]) + y = np.linspace(0, psf.shape[1], psf.shape[1]) + psf /= np.sum(psf) + col, row = np.meshgrid(x, y) + cen_row = np.uint8(np.sum(row * psf)) + cen_col = np.uint8(np.sum(col * psf)) + temp_sum = 0 + radius = 0 + while temp_sum < thresh: + radius += 1 + return_psf = psf[cen_row-radius:cen_row+radius+1, cen_col-radius:cen_col+radius+1] + temp_sum = np.sum(return_psf) + #print(radius, temp_sum) + + return return_psf, cen_row, cen_col + + +def genZernikeCoeff(num_zern, D_r0): + ''' + Just a simple function to generate random coefficients as needed, conforms to Zernike's Theory. The nollCovMat() + function is at the heart of this function. + + A note about the function call of nollCovMat in this function. The input (..., 1, 1) is done for the sake of + flexibility. One can call the function in the typical way as is stated in its description. However, for + generality, the D/r0 weighting is pushed to the "b" random vector, as the covariance matrix is merely scaled by + such value. + + :param num_zern: This is the number of Zernike basis functions/coefficients used. Should be numbers that the pyramid + rows end at. For example [1, 3, 6, 10, 15, 21, 28, 36] + :param D_r0: + :return: + ''' + C = nollCovMat(num_zern, 1, 1) + e_val, e_vec = np.linalg.eig(C) + R = np.real(e_vec * np.sqrt(e_val)) + + b = np.random.randn(int(num_zern), 1) * D_r0 ** (3.0/10.0) + a = np.matmul(R, b) + + return a + +def genZernikeCoeff_chrom(num_zern, D_r0, wvls): + ''' + Just a simple function to generate random coefficients as needed, conforms to Zernike's Theory. + Uses smallest wavelength as the reference (D/r0 is scaled accordingly) + + :param num_zern: + :param D_r0: + :return: + ''' + C = nollCovMat(num_zern, 1, 1) # The way I wrote it, has to be done :( + e_val, e_vec = np.linalg.eig(C) + R = np.real(e_vec * np.sqrt(e_val)) + wvl_rel = np.sort(wvls) + out_scale = np.zeros(len(wvls)) + + b = np.random.randn(int(num_zern), 1) * D_r0 ** (3.0/10.0) + for i in range(len(wvls)): + scale_Dr0 = ((2*np.pi/wvl_rel[0])**(6.0/5.0)) / ((2*np.pi/wvl_rel[i])**(6.0/5.0)) + out_scale[i] = ((D_r0 * scale_Dr0) ** (3.0/10.0))/(D_r0 ** (3.0/10.0)) + + a = np.matmul(R, b) + + return a, out_scale + + +def nollCovMat(Z, D, fried): + """ + This function generates the covariance matrix for a single point source. See the associated paper for details on + the matrix itself. + + :param Z: Number of Zernike basis functions/coefficients, determines the size of the matrix. + :param D: The diameter of the aperture (meters) + :param fried: The Fried parameter value + :return: + """ + C = np.zeros((Z,Z)) + for i in range(Z): + for j in range(Z): + ni, mi = nollToZernInd(i+1) + nj, mj = nollToZernInd(j+1) + if (abs(mi) == abs(mj)) and (np.mod(i - j, 2) == 0): + num = math.gamma(14.0/3.0) * math.gamma((ni + nj - 5.0/3.0)/2.0) + den = math.gamma((-ni + nj + 17.0/3.0)/2.0) * math.gamma((ni - nj + 17.0/3.0)/2.0) * \ + math.gamma((ni + nj + 23.0/3.0)/2.0) + coef1 = 0.0072 * (np.pi ** (8.0/3.0)) * ((D/fried) ** (5.0/3.0)) * np.sqrt((ni + 1) * (nj + 1)) * \ + ((-1) ** ((ni + nj - 2*abs(mi))/2.0)) + C[i, j] = coef1*num/den + else: + C[i, j] = 0 + C[0,0] = 1 + return C + + +def psfGen(N, **kwargs): + ''' + EXAMPLE USAGE - TESTING IN MAIN + temp, xx, yy, pupil = util.psfGen(256, pad_size=1024) + psf = np.abs(temp) ** 2 + print(np.min(xx.ravel()), np.max(xx.ravel())) + plt.imshow(psf / np.max(psf.ravel()), extent=[np.min(xx.ravel()), np.max(xx.ravel()), + np.min(yy.ravel()), np.max(yy.ravel())]) + plt.show() + + :param N: + :param kwargs: + :return: + ''' + wavelength = kwargs.get('wavelength', 500 * (10 ** (-9))) + pad_size = kwargs.get('pad_size', 0) + D = kwargs.get('D', 0.1) + L = kwargs.get('L', -1) + z_i = kwargs.get('z_i', 1.2) + vec = kwargs.get('coeff', np.asarray([[1], [0], [0], [0], [0], [0], [0], [0], [0]])) + x_grid, y_grid = np.meshgrid(np.linspace(-1, 1, N, endpoint=True), np.linspace(-1, 1, N, endpoint=True)) + mask = np.sqrt(x_grid ** 2 + y_grid ** 2) <= 1 + zernike_stack = zernikeGen(N, vec) + phase = np.sum(zernike_stack, axis=2) + wave = np.exp((1j * 2 * np.pi * phase)) * mask + + pad_wave = np.pad(wave, int(pad_size/2), 'constant', constant_values=0) + #c_psf = np.fft.fftshift(np.fft.fft2(pad_wave)) + h = np.fft.fftshift(np.fft.ifft2(pad_wave)) + #pad_wave = np.abs(pad_wave) ** 2 + # numpy.correlate(x, x, mode='same') + + #plt.imshow(phase * mask) + #plt.show() + + M = pad_size + N + fs = N * wavelength * z_i / D + temp = np.linspace(-fs/2, fs/2, M) + x_samp_grid, y_samp_grid = np.meshgrid(temp, -temp) + + if L == -1: + return h, x_samp_grid, y_samp_grid, phase * mask, + else: + return h, (L/z_i)*x_samp_grid, (L/z_i)*y_samp_grid, phase * mask, wave + + +def zernikeGen(N, coeff, **kwargs): + # Generating the Zernike Phase representation. + # + # This implementation uses Noll's indices. 1 -> (0,0), 2 -> (1,1), 3 -> (1, -1), 4 -> (2,0), 5 -> (2, -2), etc. + + num_coeff = coeff.size + #print(num_coeff) + + # Setting up 2D grid + x_grid, y_grid = np.meshgrid(np.linspace(-1, 1, N, endpoint=True), np.linspace(-1, 1, N, endpoint=True)) + #mask = np.sqrt(x_grid **2 + y_grid ** 2) <= 1 + #x_grid = x_grid * mask + #y_grid = y_grid * mask + + zern_out = np.zeros((N,N,num_coeff)) + for i in range(num_coeff): + zern_out[:,:,i] = coeff[i]*genZernPoly(i+1, x_grid, y_grid) + + return zern_out + + +def nollToZernInd(j): + """ + This function maps the input "j" to the (row, column) of the Zernike pyramid using the Noll numbering scheme. + + Authors: Tim van Werkhoven, Jason Saredy + See: https://github.com/tvwerkhoven/libtim-py/blob/master/libtim/zern.py + """ + if (j == 0): + raise ValueError("Noll indices start at 1, 0 is invalid.") + n = 0 + j1 = j-1 + while (j1 > n): + n += 1 + j1 -= n + m = (-1)**j * ((n % 2) + 2 * int((j1+((n+1)%2)) / 2.0 )) + + return n, m + + +def genZernPoly(index, x_grid, y_grid): + """ + This function simply + + :param index: + :param x_grid: + :param y_grid: + :return: + """ + n,m = nollToZernInd(index) + radial = radialZernike(x_grid, y_grid, (n,m)) + #print(n,m) + if m < 0: + return np.multiply(radial, np.sin(-m * np.arctan2(y_grid, x_grid))) + else: + return np.multiply(radial, np.cos(m * np.arctan2(y_grid, x_grid))) + + +def radialZernike(x_grid, y_grid, z_ind): + rho = np.sqrt(x_grid ** 2 + y_grid ** 2) + radial = np.zeros(rho.shape) + n = z_ind[0] + m = np.abs(z_ind[1]) + + for k in range(int((n - m)/2 + 1)): + #print(k) + temp = (-1) ** k * np.math.factorial(n - k) / (np.math.factorial(k) * np.math.factorial((n + m)/2 - k) + * np.math.factorial((n - m)/2 - k)) + radial += temp * rho ** (n - 2*k) + + # radial = rho ** np.reshape(np.asarray([range(int((n - m)/2 + 1))]), (int((n - m)/2 + 1), 1, 1)) + + return radial diff --git a/Turbulence_Sim_v1_python/data/cameraman.png b/Turbulence_Sim_v1_python/data/cameraman.png new file mode 100644 index 0000000000000000000000000000000000000000..7f772d36f02bb0d566384a6c740bac489d321595 GIT binary patch literal 38783 zcmV(yKOcK`6ak>-!+Z{^AuR%XQa-un#a?CJ1tg)2I2f&_YLc8;Z69V{m-TWok8 zL^LZp47FNz!dA4>MgSw41bQvfY|`O^Xga{MfMN_qM}dt_cq{|aMLH{oWdp4dEXXXe zTO+t9Jb>&N7*QH934j%YKr2e^FpFq&U$70Y^h5mI6vf70T#7vMD#(1$h|lg6UZ8ec74r zB`8Np6^a&`B}1h4mVk=xj#vw$nk8f|5hSb*iMi03q!9L^P7Psr34`nmE|E(Cu1o{~ z1*r(eS0zd<7!nBxMATxa#Tr2;4VuV;WfX^CIxIkehMGYz!3qOFEv8IM2$)LCq)Z9T zGy`D*9zl(zTNs9;vK(!*X{t%15ookDf-F}$!H(wWvH`H8V3cW>X%<=5kcPyXpKRFO zsMC9nw&14S{{OFl_^L*5pbfhfJwuosq#;@a637XfN=H7x}qb}|`2Bq;_#RS?lZCBgxN$<7v+5@C{7HtljZaa>!V zv*qZ53a}dlKqZENifF>IBYSDU3Q$oC6{5Mg7cf#yC#f+45G8;?IY0!^a?>D3hiR}aqzTF-P{E=~ z;gJD@Y@rAz*cmhlxm^HFFhIK-=4Oz^2F-w@s%?VkW+QHILd4g0xJ=7bi>#nz)~=RY zfr2^_)gc4KNUA|*r)6b8w8a1t7S$FVLfT;>9A;1y#U#@7vPCcoX&EpOq{=YiphTFk zM%o4numzi_g32I-Dvh8G*yVx|7AzMA-O++%4^lQ8Z5rtq=!Hh0bp};Nln|ZVDedbn z)`Hv?9%`Js#~Mz%v5-J#<5F7=$enwVKpGaQ7cTd9X`Bwkv0#t3OnXmQ9c{R3G0d@I zEKO~yKxGl?LISZbC3p0Qq82v=nbVTZbh47D|fM3>z2KRVZ+<%kpn0ocaO89 zwuQw#SryxQu^e02GDcS+zP4g8a${7w93r`o^IG;mfa$R27DTmy*omQbz!l}gK3|eq zOQVk>;L;v$6kQmmd!Wk2a-rwu9prAK(d{-|gK*d8?i#E;+}vY0E@Bu1xD~9b0AXXm zOOtA#)ofM)CW64oupQ}AR743X7zRk%M%pF~(5$5> zXM-ZR3MF>a3T10k_W(^qse*>e+8hJ|luKTb%?NmRn;=GrLJUR@4V#7voi;?nY#U8mHS^_vI({#E?sk>+oG~|Y{ zOQSo6j4Fb16b-WwO_LoEqUFkf5eJ(jqoh7-Drb4ZlI4J1BE?!kF)Bu&Wk3+k=s+sI za;Pgn6jT;V%_3m7gk@->Roy*U=)j1m(k`^j(T9ek%>{<%y8_&4idh*UIV}YsFAM}z z9pwyXK#X3py>3JE0=|lq>_O}0kv`jk`4!8gG+;mQ~k444uGqvv8Wa+Qew?zyO?OW2+_fo+n}YLB5)`%Mp$msoQ0uC%qp=#dU2B$A_Pz) zAdM`6+nEMZ5?b2VCh2C2c61AgVo?KJI=(`}e4;K22G|LyV+b~Fwj5@H0-_gUem|x_ zCAd>9NJJPy4hTUy-6k7hDmz$ZNOU6rqUQHy(z5Mzqcwz<8*J0q)W%>Js?cu8YM|h* zGTk=28ywa0?U*Mj$0j$*U6jl&fC#acHi!`hJ!o%^uNRVRDJW7=%J5>#OL_sUSUV^i zNfK~Vfp!#x9D^iH1yNE{Jxdijn@y0>vP*UlagjTOv9#x7kK$-F=%@e;*oj6FMVA`d z0^Mk+P%*l#==o$@1Bq_+SUfA(2H5Rzx62ajAz1Jr)b?Lk#b8w8>Z&-aui8+mx_SK4mgg4OcND5+_s8l z4Gva!AGNdgYI=`Rh@BO+n{30v4zYSa+f*GRhDdavv)LLM(uzt3fK|#O6;23HS`6WA zB$C9zQX`i5ss92Mw>%-zL-?2FkI1bwGYDAH{1`9d@W%Y^>&}ayb0!T%Rvx8w7qdbozX4-12O}L4i&p0FkgyHCbhR|4O zt_y&sTPcP`$W4keNE4~{txC{nQ41)!jg*RTt7PxJj^K!HVC$$_cI`@6R2NNSODj;a zV_)c2H=~6#xGw;*Sa4g&>R1a&mlSNVGq{2`>cSX zb}9>6Z!1cjU9lN~WeYDATLH7+5-Qm94{xbiP-#KYBT7bztw?C7V&uqZ#8^wtGxhFtW(6Ob9gAI&B zqf8=g5Hdi)up}A=HpXs=f>CX>Utshu%LTx-4h}-QP_dx_Rud&RR0Gw`?wz(>M>e`l z)0ML%r>o;@f`mgR!vHOJcEQd@4x)`N*ecCe>4adTD_e9)FlIq(u>mv#cI;UavMg!{ zwC452CW@6uz>RtNl^O+>htai-o$k&7%R9Plpc`D}28~jw$TmSWS=p8hV}XNkSF{Ku zNEgUr!fI9;1OoxqT#>L=1}GbXh-hX5Tm_mcYRSbW0Fz{~Mj2QJ>j-mmv%4tEd$-zl z*+!AgAvN4$H6%c-aLE{|Y6=AxBQcopyLvRrJ>)ib zw);4`sw7^ty<>E_1Po+baM-3{*bxh~$8Rw>(8YEI2ZfTM4rrhPCTMqc7E0|2#Tn>g z87O*(AoMzh8!=P?jvR|=aIk6jvY~7^nshbYpovZMLMI3TD1~j;0uB&#FC?Nnm1AwS z%1q}-GX<Y0JxU}!uelLs88F=r!Vwp$>UK+n1}zflD6id8iC7T|g0e%56r+X603)<8pszy)fl1zg%yTStTO zHL9?g1kB)IvjNgj&YPz?7J*$k1T>?JlsX!U!TI`yfG5;TrW?fExMmSGv^hvOX$S@&Lu*75+ zb%x$2S%T=SiJ_s`*#U=WCjfLr1`H1s&MM5<)6otWOO-9>&@Mn7JJBs;q_OM;wY;?U z3bbB>7Q-DKz4DP5NvKe5(hT1^KS2va5Vj@ux8ZjiRv$0e( z44AMNNdiD4y%4=vwn?yK?i@HKAg-~nngbx^BBSV3mSiZcmF9)Jj;o8~ToOfeYIw2TffY;ZvyMbb`&5xp2| z2b36&uhCNIGaHf@sPniex+NFO*wYr;pf=OvZ0s>06umSnsz`X0R8fT(*!M=9TKTyW zCkr*IU{l&q>b7H`_uYVtMAK+DS@>#T}}f)PomX*recjusIF4FzP-i?JV~G>eTe zxVAYQWM_mpmYqQx^u9J&4sWaCE4!gCSuRs$a;|#NRv=}%R)8as7N|a0Be4`9AR~KO zkt4Yz#9WZk%?vqAIYMIgk{J;Sd1X*ibgBdCga|N2A%OW8tKAY7XbDyrG#5v-m<0wKYV6cN?-H6P z4Bd_r?1rUITZNrrRSu&7^<=QYai_ zK{Pqm3Ro5bWB?uV0jcVR0A!_MhA9QW0E%RaA}1#tiCPT@QVmDbBG8>B_HakTAR0B- zDQN&;cXqo0VQ{J7sF7`p5?d4J9RRzd(P}vxO*JBh)W`-LN32yr0+Jk zqQEhhtYtexk%^VTlrau2nfMprb8TfTxuFE`P5uE^C*-&~#dGfuWKxntO#=IXjf#GUm zR0G)E=!6&mWj4a0E(qKOaDix4G=lIbLg;-7c1Tin$JYuNPQfaG0kPmPP^pG`MHCFq zJy|Nt#EES+(t@$elMU>7{nHL79B4A;;`Op>yU1L$0Eu*p=~*NQCWPjQb~=nGjRdT4 zq?S%ZfkIOg4l`6MIhM{!NsZ7D80r8;H@$2D333ovi4KxXGy-lZfE?IDxh=iyStxQy z-1fy8*wy3+Mr_Fuog-pf8}0~O2&$xRB;ZOgBy4Bp;Nn`K1lQ_X(yPXk+k5`mb?5o( z-u$xrUw-M_;qEY+jA#jA=2%;hx?6JYOVdN%3Dg7&o2?;C>;`-6#n^$GWMqgnAOI~y zmxUom!;ZQDqN(miz&+SR)zZzs+qP^un*?&G z1XEFga4&nB#OWKSXjO1@5**p4NOx*D=W|xLMYu!luxQ{APyOAW*>JeC8JLANfGdCE zEn_rXF~B})*g!XG9bxyzJG}_^ExMM%C&#i(s z+2H~)6eQsgmZT$Uw$d zgYJ3EV?%eb4Z`d;+IBMt(=}Oju}bD>KO9c zOic9I96xj4OP)M_@ZJN+Km|Jy0vC%g)!{(`U7f+6mjz->z@;aC#HQdB?Wu^u)V#(n zz=A-?a&Sr+tZE3q+0c-Y)4%2<;Fe)8`trpO*(E?R=vpP#2%M>;ms=#JP zN14Drze_RK=YPsa{-77K>Ys(>vi>O$FOr;}`in#S77;W1y88)Ihxqj^LeC&Mt_n&dd zMllj?N=P@D0_z~KOiL=wN}H{AqQf~CuGCQMqFU^5HMP9?#sfmakY1BSt9iuU_>;~y#2Uy*MP46mT7>!o9vqlv>hIW_R zM8{4ydXs5Sbr#Jw13;loVcBiR*JTdI!bG2q&?a`0*_dsxXhtv%5`zbVvjn1owU_}2 zX;Hxr)Rgb4A#!3J`r{vJ%x;4)LNt^+WGrpq>F03co>zcy+1-J_pv?*k4HvQ{1E}N* zMpT+!8fpiDbSp@Rz#_qEAe5!M$CesMAkZ`+S_mzpq)^lzpk*xsjZ3X;!=V(ppj$#4 z&al^jVzxHsUP6I;C)lAxHMC^3`z$MLO9360xBZId+Xq3iDnlBKRdz_g@#fQSA1t-e z>}@ZHiHh!W+0xD`qPur*t{O4yxG;8i??CTB4Ypa{s$+2puG=-*;b?b`u{awFCk#0I@`bfnc?bFx|ks9gRb7 zyMt{;m;*8_RVpW|*F3JaVFYnn2(-YtVqroPcr)L7A z1e}*!byi-7w#k9Zl|2Yypx`3A0Yinz9ad!*r3+|^sncMW?UG^DIE);JBZuJTF*zHH z?3VW47kNh|aYp-OJRmY!P$f$3K-Xebg^3YX4+GxT^ic%%QmG-M5^y1Ek}fSNHdz_x z?>0=FJ})rq=eteO)P(_l<=MrhJ%@433T{(;%29%Q%A?s}3}PUwMfC37<%ZD!&;T23 z%j7zWPn}TLs$)A(r5to!uA=_N=3-;RjcC zqpJ&h+t8Tghl&Fw*gKAwW(DB+4BXXaXgV6Sp3UEl&pn)Dg@qVq2gL%*B3Y34sZcLk zk(OPva#=_cXtw9OYY2o}6xk{6NYM-g6#`d#6%}cVuF3(<9(lx4VnUmy>U%rb(T5;Q05OfQ#Sed(t+Hu&=6_utZ` zSHA86sB_kAY%hD%(kh1ylnUV<(Qc?4AWEn(W!c*jQ@%T;8FkobI(ZCIql4{b!flOa zUDJ9FQnh1%0qD4RTb2!gMynv?(a=!) z;NSg(eOl`DnKzt{L@xI{a^+>8dgWK1_KJA+17G;|=O26XmH+UWOK4xaqY>2$qY*4b zM*-O&7VW{E$c{azI#{D4JxE8@qJ4=DkEwAT%f7glsm$PU8KLlW^96Tl=ih@2#&M;f zZLvneDETh#l+_H=gxp}5og;If(#S@(Fr5)NK`WLZ7KzXge%p(80QlxN-TblZb2IS# z(Jy(^d0ZMl{!d==(hofK(igOxT^j3e8J}{$F$7Lxbznl??d2r;t6Fz;3|b`^#c&a| zYzfoGpooG!70e?Hg8Nu{Ds>FNQbAS0EzM1uFbbXYPP1WWtu-WIi?uY+2s`QWopuFpnB zU2uH&%Pv0m^0xLKmf;FkhfQW_1cxJ7^O>-)Pfr-5fT2OdO;maDzHqeUMz_I=9Lq4H zx1~ik+*{-j<7=HMLINtAh@_oZ%g{<95J(F+4;b^KN$k~Df=~pDltkDaQ~A6aX{Gs# zAE?uDk9&Uao9}zuvjYb4>Ob+R=ZxD=|NKQ>z4h<7hL%fN^8L>}|KuCFDKN`G4nYw; z&rGx^(W33#l`x~SFf_6CMzYy(gmIL*s!I*Qqg74Mp2ip^Q7EO*(T*Ym09$m<+Ap|xia6>_YAbW-&sz&t?00atc@I^%gq?ALLmzq6hM`V zQ(nDv4M#6*wg@kSC`MDUrPWSQr2?`>f-$lNNS|$l+dnJCP>wB)5Ud^H0a_b_I9@7M z8~f4?FL23dKt=C_i$-HtbTFQJ>pmehP&c>df8gP_y!xN~DEOxPZhEhyJAT8SPV8wwg#|+7*aNkVaA@gz1le#@>oXYhAd7slZV&griT{{BHhfU5j!cz zzBDVEAUjF`$##qX`lh1iHGt=j*YA1ub8dXbe?0ajVtL6eedzZ6zy34%Z*KjtAO7(5 z_wH1IHv5ZT@}2hj_go}fm8~wSiRrY#;FM+=O=dTjR2&2JHZ_5DTi!7rCs~P*Uc`1~3wi%<@**4j%WE;mb?#x}#$aYPgWv~KK zTHoo;(=D*>TU>wqm*U#qhx0AhOP}%JohN?pOK{BQyo!LZR=g)pyb>%Z8Mca>}AgeG>8UACJbLJFw)FU zwH;HY7E^m3h2Vl(;2as-LMS$-gAVc%7z;w|W!LJF7(r}j`^t6sq4%ANej~L%|LpJk zqfcD^FwP#i@A>r*>uc8AuQ*%R|L(7S{pZGFp$_;_ATUFx~xdiC(x-?s;U)5o5F)%!0lUU|O1;-&xo;T6Al#Rk+bvB$hNJ49Wl%iy_6uJI(i6@pC7Mv%fkkuNRWtoUh({m zrTSm5FWL7?WG;i+?me&PAAil^_L^M!zGvU@$yfi)-nS zx${^4$CrP7(3KuYX_v7ld*|T@^>jbL%yOWhJ+fP4I+mlu9o?+6t%G44LaizG}A^uqgYJbcE(rSn_;$aCwZClA|0@BZ??{@%a;36>Z9o9k}_ z-&xGzURqhPK?s6i`XV&4iKUGWlWH6~-Es>{L0dRR7!iGpDL1c@+_7|Ekcd34@eQyv zrI1mIQVao%5C|0mVuQ`_5W=Ol=-?{!&|H=#P&%g7lmh$F!(}T2# zCc6VWa$e0TYYI=Ssjcd?i2{3m$B;Wr&_jZN=v@xLjbXy(&Lz8K#xPq>$bp_ zEO42mUdkkVMl@*WSSzR|?kc^4b@^{Ua{~QK(=%)jP5%fq06zZzdf#>)mp}FThhP4& z{b^P2t@7hnZyfRD;p1QUWk3GFQ@eiOKl`0e(pH=<%^-q2I+3R+l<(*G^^!(2~=)s1Wx zh}nTnspKwn2tbUq7sHy8h?#5QD9KW8-+Cud_H@5D)iemQPcjH>hmU>wm;PM(=AG?G z_3^3Waom33T>YL$`rlrB@>>cDzW)b5g6HJt9{EX<{FO`i;H#G9RUiM@pE%3q_SGx5 zHFRkUL$=(+wA>bQ3%l7MI?M*6rm@try2a)m1)8G=vgzqs1WOCq%9s}(my;@+Tcq1^ zC)_l|7>Q0{P_wvhQxXEwcRjI`Rr9G#jU^^eYZ@freBJkd&mZOWfA~A^d*k0fe$`Lh{F#3byy{O! zGg{1}b7PLM*B3GlTz?Ix268z4NM!})LdUX7bKHj!{WzoZ5p zcfmGFl$4etV^PRu?iJ$(eFUn6mJaWx(bW~m3j99!AQ$QdN7o(d~y~u10 z!Cq^mM>rs$fq}Bk*g4AB@|4k3Zn8=_3@47^Wy?W0xi55cbR)!e>y#NY{i|-YI8qL8 zE9>whf9uW(o|=o^CQ7tn3&j8%@Xvkh%JS?huU)$GvCB_APaAdpqR;PLKIO45DOVge z4h$d2M{e8#FP?hcyY3wS=^wu7KY01-ANvdZ>Yso1qwg}#zJfANSsP$!kBaS3)0hT2 z#$v$MAPgeV5E`8oP1*spWY7YOYjdQLBc#%hh}NPJAuG6ZDsz%UokI+nV$n;9@*%wb zyKlen3G(L*%`m~Ksbnkg?8E}MpS$(1Klc?c|ImuF2VR(J(&Mi*Z>?v|CcKJSsIiYk z;Dc}c&?mp+qx#Tyee%a{T>JhH`8%(T>mbl+?J;t)kuY~{+-zH8DV8y8pb}-PMp$T~ zu`!HkwPJPm-HE1HZ>(UjBxB}hXW3t-J2)S{@Ntn5-RfD)F8 z9KuLBDQ^OzCTB;QLDbQ$yTbVr$I8u7VGgw@to)w;?;SKGH*li53hbF81DG!|V7&cx z@A((J{DHH5xjIF@1pCXczVcHx)4ktRL^sHbtDpArmp!}v=$&QjFTMOr|JU0+*1kBa zjcq~u4MEK=%y1ba?PbirhbrVQqgb}LKs3aiw6i*<&Yh;QG_8(WmeI_fc8wUTVM|5t zSuxdDQ9$m?4(_SgR~7|%>9zjh_m9c*sBvQXku@=YbN;($xK;S0U;46-zVEp@>WPaf zhl8EV=hEN%UE4P!&{&xgQ0?Hhj6m#~J)l|a zEoAO;bo2+wJJFE^9Z83l53pTo9|2Hf3nC=>S5|NLEy%a=a-(eul75((HpoN1`(N1uEp zJ=o(l1gYalU}wT_?OSmOu7>Rq9#v$ut~Q- zPX54wKEL?LPyf2*1J`bU$@LG!ewP;h*wOSzHDuC9U=UyD}VOv&P_UbgSd6? zSgOJm!^MeJHFAXaDTIiMHaJ?42`{@U#KG23RcKJ&14UHaD`{JQD)6W3R7=X!0gD}(wczkR;_3m50`CmU( z7!QxCj-f8e|JsmPuCh5zg==;lS|fAwwO@-V9( zxc8U0EJp26Ki`k`k6r%UqfxquhKd{LE(c0Ip98QPK0$V1L3?1(?B1JYj^r?CwoVY5 zhtOn=!@!VR@wGPTaGq>t#2^JS%;NL_MaM*TDIMPPmp@ojL_06O?lxx&aDqYRy$FQ? zx687>=}Uh8y<0@CWtspoZhiI@fbM?P&0pC^-MV%EuiUt}{=%ny=B0PvxK-Y+-Foti zAy^2=n9E#ZAkmSXR(er(OcNWYWu3teQgk95ibh%%2?il&u$Wc{9f+?L!^OaJK z&a_YW3jv2=h|r=R`)ePY29GK2I4%5>24I#&5CG@>ufN{K0-P1U=R1CM%V6VDPoO;Z z{ont(D=xAQpLylRy)SuS@#Qyu^Ubqwsvo<7M3?L4=U;L+Xa`gfjg$X8(e9|89%K@* ztA$vGno(hh0+#lm=268?3k5KBp{>wKHe3k|dR{$~G@F``LaFj`?5>A2E)Yv4wn`%|E#TXFP1uorI$#59^6l#+_ z*1aHQ8jWg+{LO&=_&)|MAX` zidd>v4Djk#ef_;p-N+TuI9~ho<^O*B`F+oE_Ui58BJO$RcijdDue{+C23#r!msYEk zLrp?uPHE%Lkg)@S!gQOC;F0tUvdANCw_VhcgS3bCu%lQpjLZrT3Jqw=XNv4>&yC{~ zGTaf!+yDE|H0Ckhf?##K5vFnXc<0@%AlFO)F<@P~@huj?R3Gp3>;L5beee9>kt<(z z@_$LxTmHv;MHj?=?oEFQp+4O6NbW_#VK;V&L&6O1nqf?axlgA-&A{1RrY-KB zFj$@CKwFi2*z$xdjNM@8j*+YPWIJcVjz$)kZed?-!R)2%q;mY#p96qA?Rjbz;Pxqm z4u*^gM$wq5JlPF`-~Z?!xP#p=7&~N4qZzUZ7Re?D}HG<9|63mOg`Hq?81M>!$WTc0FGVe z!{7QNfA0muUf!Tzee}V@r*Hkzw|>EIJ(M#$8vnQ3pL9uq4?XuyZfO~GG-d*Cj9`{? z?DD`#d^cZ4bKTGqfnqOcFWXG&*r)slXss43l<1);)y+0rd(?4hu?ZX9odYne)+x+l zsc-+Op04`SVsFh#$9J`k#6*oxTcH7H*LHJAwNPkTf|IV=vk(=J8H2Iw(^cg47IA0vO>FK*V6(ss{VYUwrT1dp0=r4Y=<$^#%?c|L`XZ|G}S>*7$c{@XBY_!_J$p zh&E!_MLK5yQVW8LHD5!Amz_mu3|25W-9H}7JSK{F@8y`O3hLami0LUJm=2X((-A=m zXmuKB2ZGbSb=d#r+e$c@dI3;?X#$THU_JDpIPYV|Lc=~=V|F~ z>u~LGN7nw)-zLU4@m%J|?tJpwZr=j6)G*xZaovLnU@cIUe#wROld>A*JE_H#liKd@hWs?@;{W*mXT0y*Klq>j z{3rd%TL7Uie#?$34t2S`_l#((Q7|M|Pa?mYVUG=B8s0HOWeS=}mkKRLp%lx4gxjQ| zT68$Ag(bo$IDnnk~6@be{7k;O9hyOt%|A7H2}@e3*?;Yg_)y&;Qn+K>(&}z9-|E$3Nx1 zJ8$_HPyD)%`cfSAIUoCd|4ysytAES|%Y8!JXAak{nrHRCbkkNh7c|tNp3y4}4^e}A zPkw|}&#GPB9?3y~!|Lb?#8<~jOQWePAeFPKrof0?{=vH%0NeAdaf6D`UCj5C8lmY4 z14)6ny9am*Q3LVWzwIylTnBqr4yW_Nb>-ju+R3{pufidpuHm`**); zYXq#z)v|idY-zX=Q<4vHnw`L!I$xu#8C^_yMh@7W0}J34MmAcQc#>cQklCDb9Hd&% zskY2$uNrDokN<>a^ptzdEa=FF>s0=9o19|{PK|=+$y29z$rD%j<^TDg`-SP?0`nv! zA^5Y;+&El$L2tbCGr#Ne?tT8-w-`?!*8K}UdrOV3z(tenb(;rWU=-#dHOFwYN`Yp_ zuxA`&!5l2QMIECU)c`lh1Jr1hwwO6DBXR`0sGTvzr8_lfk7Zj*Y34uLW^}`ZfFOa` z$}>9bMJDJ(3**IeaSF{e#h*4m_jmuLUueNIjxopG0`lX}8hG_<|0#eUS7#rz-{z%7dI|j8r3E+BlU71-99@e(a6DN2dLVV#ab}hy(gdr z0}%q-)$Ku-X{$J0-Zy^Ph7C!vEX$G-Nh)%gOVFkpbJBzF{?MKn(u_DT_yvFVJAS4$ zm#k5B*F-QM*_Ag32{#M#&7OH7TMU-T~u7ciM&(nWq*IE~H>DO9Eto;>j914}m$e z2{1*8ls^4``Kv!#bLfsS<6G~3JqqF%{z{+FwXQ+#?|pY7Pl4JEopKsU76v%?j#0CE zi_(~6(G(YF5=6&Jfc{FjJ>pe_H=1Y;-ck22gFo{={$jJ3m-SxlVUQ-rWp-;Lpfr5spO3J-T~*HvFOeOIF9sCIcPKJsA}bUD=CEV(QKPszV0#2@$6 zjnhqlZu`E!@pJ$1ruXAEs=GaN1ZR*2KKU~*U~mm%JKEhzfY|>r>Bc_jn#RB;*bN~7 z2gIO;Ck=dYj4Ss&`e=Bds8jvtW@x3AMn_VSA`%f&Fo`}SLUQ{c&-j2fDF6|4hPjui ztDiYlh_+G9$Xu5U5t0N-T5`f?cfm0F#+_NN1Ko|Qf8~SU{;Xm$%{JOSMsYR`=lZq= z<4nR?U|IzM^e?=C1{DQ5(OomM!#0dI+KoZ)Yznvlt@D-C*7hh@xomnPP#yHFkl3vn zr(!1CZNxVMY->g+bW|cv={Hs+0xF8pr?X>dHS6=`;%w? z>yt7`h!iO>$9! zlY08rnSA5Mp75M5G7-W#3Bu9ye5QJtqaQlMG0*7uMu;NVLX4SW*8`Zo?gqx?vvHYWonO{2lXZ3+^$VXrKv{_tsm7lkq!`fxI z7zF^Kp2l#yVs5$d@tfCP{?zftoxHT#W;9TQdfge>P^l;sijxNKh#fTC&0adkx~~`O zNnV;m-O-lOj6tgEbH&-F(1@? z>@$f8rwhX@#rgH^_$6Qe?zcU#(FVt|jMOUx^;9{4*$xh=ItnZep)}Y~>8kI2ES|rMU$Zm6K+hMv!Hi#JEJ|6PC)Q5@=SF&snOU0e@;SsSPc3S)cf9>rj z@9TVVJfR{qg|<#NZBesj0wKKe^M3jl?j&vw1YqQ{_pFpKr;3>XuD+>h0&Y?w@e&=RbOb4cjcZLtN!>xk3hg)0U_m&Z7FWq=x@3)oq%;M$DH6N;(zD%4V zQ)TP8baU6xI+}lfC|#WNo^vcS&3#Tk=tL`^{V%_3-dZp2@skcp1l4v2!6K+p^OSgy z*T)=zcM#y_))r<2{yc^OGZtSyKfiWdzUX87Ghh1y*RH`00}uZGPbvaEjRzW3Pk#jF z02M)V#MjMNp8-PC<%<`Lsf~)7v2Sf-WNA5a&vpkerhtJi0F8>T3t9@toJKWqO8Q*M7I2}$P(F?1oU|R|n*d6=fEx&jjBzY?PaW@pejm4k3-@$2Q z5UGrpK6vSTy|O7#wj=^NwE!n8`LCsFEN^;dp03Bf^0FrmBg;7Z}*CM98$Z1k;D8@wdieP-zjI5PgU5imwgk5G6 z4Xf$4&Ohxv&-;{yGoQqKGj=9<82Dpf!9hri<&`mXqacV}tqwr!yM3{D_?smYid3MQeEe||WJYSaaYe(lCn9?T&DFEuf`PpCSiJ!@_ zMVvo<<<|3`Wp)g<5Y3u%T?hhl3ZNMb2N&Faa6xE?gRI+B8k&(vaMC$DcP6YM+1YB6 zMyvB~d{{$j=2%B`34zXlEj7*nIOPH3@OvKkSGP|P_+6F6R0(?v@szwOles`e45Xd>5rKH+dq2ALm;R~W zbn92OTsqtuH*d6V-t*crT%4zDm)%-mvNi96Ac!)@3lA@H$$*!>4XY#QIWY#9>p_qi zoBiCCMCsBkr8f+1TJdcR6WcbP!E1 zg`LhlpgG$Hwf4~1yMx-jS2M?`YPYGf7nI=Lrx=O~FYoS&2-#CAIKTQox!nzny>D&D zW(K!gH51J+Lsr`E37*XL`JfR4t+);Vgwvz_0&a_$k`4cY&-dYYu(^pv@riC3;C z$$S3bZ~D)KZ9o4vjK%4kp}Cf6j)}^6R|KQQNWki>X=_6V>Y`)5APjg8l0(}eSBr{0 zAwuP-_JnPNvRg%~5kzL9Q1g(Z+l^(qQA6WnCu6{g9Xfqc*NOMx*2$n@fRW6W#o44=Q zfwPv8l&RAQyh#|?=wsdEbshUz+-#r%jZwl8um5w8Y~#{r z+P4qOKO$W{HklqXq0mj)HL7=%izwAaZPh!4?o&aS*rKZog?4M-xv>l;(6eQM5uKS7 zMC9bDXApCC$B3E|pc#<|izBiS?Zk&R-z@}il9Wt=FlQ?CoCdn#(kSkGUwZB$?#PEz zHrA;F;E5C>fQLTh)8vr2_`KhGaPRT8KfC#T5551{tEk{M8{oU>gNk`~#=RUsFD-If ztC_Q?4XPRY9JtMJ7NP>`8I0#401mL zRg)#xoR(Z|>^pAXDqLLLdy(sYyx17wI}W3We#L9( z=kDWoW`QLJ)t(uMQ#V{$bj=~$cDeR_+lGm1ONW>vZ2;dD{ugGdC@b88Q|=l2Gs*n` z`({JPQn46aWEMM6A(p{?S9`emKD@X9<}rTvCYbq#_T*5v(Z0BS=~-TR&108XoBJdm zuz@f6{hyO}tycV;r+gqIp$H9(-mAVWIBs73>09tFsa&GDGz)C*QmqrR5(3F3 zhJf18H8T}OD%Qi

Ztj4d+M<#~jyT2GWMS)j+`2r9!c)nCuuB#O6fVK}@3VBfr+i zo&b_5nCJ#zblHZ>$hvd5dFSI#fB5#LN3H~5_k@sr_No5;U$ovO1>bnb*USTr1KAaX z{TYwIzxvAKk0!&l^+^q`naGe3YFd2&r=p79M;TS86zMo2b~YY+*~5=KdKl5o8qMD4 ze8c_peOIXkDt+JytqugkRa2fHp>l|bn)Bvg_%$%&?o*hQsB2D*nsuVE9d!)U@z$Le zjt^eBKj(uZgx5UK{X2g9)d_yBj}JBaGJr=8MxexzUqX(X=RZ=x;sQ+eZqyiVjGC52 z6J6fDx9VuQy2Y|x^j766^m?AlVOg(UKD)9y$PyF{AYdOO~DE*^Qs*%T4m_=@T=zW9%SVhF%X z-rZeLDRA~c!BUT|FFvC`@Ve*TV~T7H0JOw3IxGTJWJ)!U#+><8W5nLv%i(r}&mVsD zQa^oTz5iv8J@#I?e~Jp8xShLq znr&tV@Ecab?YjQ)dmmWaC$20u?8{F-8=$q(KSY)_Z}*a;)Sw{BSJGZ9?o^;#rr?*b&ch~Xt`6( zoltj)$^<)W%I=CZVqdzET^@wZ(P?>3$*w-cJN{GRV1MGn6n||`b?|#+$uRjonvq-Y);^MeZ zNLn$%O*G_eE=p&w6+tKEp3Zny@ZuSt*MsqdjXL0zd1+)rR83qohL(A7NLrF2QX~-> zVK0+`hf#K(ko>!{sS#5i*+cM#!Z0{zVX=rxvgWvWgZYEk6nU!t zBO)d`JT;kQ#N^2cMl9=+aNWFiqagG}%`M;u9)&tR#E;!H&aaKLhmJs?SQDoV!ejgD z1@)l^f8u5s<2Ww1@_PTH56KW$oD>Pm(hzuORNG0RyWMmF*$;Cf9+1nTPTF$x z4ly{>kg-TL7&F3tJ*j@lP|sd#pR?o5hPZuywqU10m)dMmGYS5*Ib$_1C$`CY$eS@Y8=0pJ2l*MpQTH?(qKf zGp~QV_Z#>8$TzJ(22Qi6?mZh%J#=rlun87#&AD+qNpC`N3duEHOWg_>VYzNyj@@!Q z%!-oS>}a++6Dpyn1#BoW1teu8hAF$~4(itT`>q-eJyX&N&N=e{`c&ODaF;#>W-9*O z%E5T;Bm0yPrhgga8-BuXoN>(4Wz=D(Z++yceB;}m$A@13%RlycwxXjJB0T!OifXr9 zHZ=!@cS{=0y>xiM!)r=1pN3 z=br%^JZtjFqj70J#)IAvN$`+nM)$H^)ZzTz-~96z{-p;WzxkRy79d!)#NH7;-r8jB z8cO7GLG4tdvca6Mqyr8}jw7tR$dxezg=$W2l&8-}Q9=$WYHzeKQZ&ll0xd_!)VX#l z|LKKKhUw|uH=+4?e+}z4*poCrUHjsR6q%;_H(cK7bYRx-&*KTC822z3x4AdB$>_12vFFkx1IL zPX*V&j170QfRP==vN~_}IegrUc{b7P(?(%5mmOv%~vP)1#v0-j%ENp1L_ybWh99#K&{5s%8V*MQz)8-m{)hmvB~P zlWeo|#*g;g?lU6TB|aX7uhu7k9hEEH^l6R@QO&q3eLVszwtM<@SpU zo`53fB#VjqPtoAyoln!dV18xe-m91Id+5H)m)s8`>Ar+02^H!r=F*q%*xa1jj-}Bi&=nJm@{PXI)Z@=&U zx_rh3anV3<8O&Gat%ji)TggZx1EzzWmBU2`7Nv0ahF+2$FkP|S%=+#$Z%&FrA8-`l-$SOeUFIsER~Q|V+8~* zx8wXfYl`CT9`8n>7?T<{X2!dZA3K_-iCS*Y5sD|pEPVf$on3voIv84(>j1v~ultK% z#vW*galo06`}xf~u~Z+SRS&=Bmwf$mKl;6Y=;gO`&+X2A;~k5m>T=&aUv=;8nG3KC zGZj*^E81u?cGr2K!+onSCGT9Ew{jl|V^34P=5Yw_tV4x4+*WUgg@;Sg=5}~UaqISj zXT8R?JMVbMPaC0`8Qv`v5cJxCe>AU=58(C7R~ymlPt z2q|GdEVq1i&y73#@shT~==;vT>fP>l-2b*~aqEBr)?0Z{T$GOoRT=l*ckAMKV56Hu z1FQ|`7&QBWU3boM8PcH+dCwiK2V*RYvw9JV%pBri;;7Nf?R&!a6hV*Uao1Voy_fmK zGglTO?!WpM`))CR_X)_=Mt3P^gkn$bOMAwPcfi~)audYiYp^-y5L~{u=;U#E=^+^M z!OmTA<$>~jFWJwpUfGv={0=_p(I@fbfA$}K^y^}(K$g&Ap7#}}(6xux^}dVSgt|v9 z*1h|kDD9*7i&gvZGLVO}H3xZqd8}qF8j~3_p|fSU(>;RwMV(+yU3=XG2lEMkJa;1Ea#5+2>LeeC9)*yjx3;sGOfI}XRrM}7V5o(t+VpMUusAARti z;{z{z?B%|C{ryjW^fsP8J5(G$^27aNg?riqr{a#25pwTyj3of#EIbM;>9(x48>YUT zBAYQ5^7tg7&--4y{&t^t>SQ=Ozrqw4*7F}_-?;P4&0DX2r*G%_kN&`&ao$w`j?oA2 zVxevwuO4)MeD0Z>hl|6KhsuQs^E_k~6=m_T4EM@)+;e|jeEjN~iv9TfnwQR^JkHjO zJ&&JT@4wdUt9s&@*S_ZX%meE)AHTrJtMSyeXCooqy6gU!7r?KDB%M$ZqxoWjj=A5F z%pIaJFhT9b{f1Pw$)&Zt|CWI_{w=>rH`Z$}ZAL_ZSGR1C8dxu(aL*xcz2V)PKXBhO zdex%tkxO@b*r^Z`%?A!CqW|Q3%2Mv}YTas-+a$q@d5+^VmAaNWdo)R0f4G)rS z70W7Ac8TOC-e}*lc$9z1?@LxlYxY_xMgLAI_=U%-@izFLgei zU3>3$t!F*YZ-K7unm*aym`rek&W7UQZfH_pD!sq&@cxOuZx%B*?MLOt<$Fg{<P! zKT&`$zW>S|a!~;+s{sNa&;XusX2>BUAoRJPDvByOLWE$bdZ3mxBhJbB2&7zeNPy zp=G~rd{n2N$N}T_-B~K8skN*?$C(&DPUH!`O`seH4I)Ci-kKeGb zy!hrj;msbq^$o(C-o;#L%? zAPS#)3w#ybZnv4?>v835#>k|f%ikvfSRH(D;!&jSL8rg&xn#_uK95TP>1=&H>BD%jH6XXS zS9kV{=9aK2c7CU|{=t(syk@bF1e{x0 zHW0r470^!ot50*}k_Ulii6%R{Z{~XM-+td(WIz`6{P?}Kr=H%)IC=)atZ*(sj}8rS zpPvqg!~_UWj4S5twCr|BRE>~Gf&n3SoH#K+5LPBvqmGY$1}RFhb#+5#%}8B zm9JNH+TZ)M6WS%6#JeB~nAK|b{onk)1rF#rAoEZk9DnQPi&s06TdxE}Gllsvx99H{9mUsS= zqzo76Nbtl!M6Y_2;#xaHN4IX3s;3Ns2O@_a7zbA`h;wHu>%}!1Mzn+dy+=) zdbZl`xl9UO9EXo=p)1Sq%2h=O_(|EYlu!}a6C9xv-OH{Q2aMr25#YrS++*AD&IJ9_ zUwrD@d;YgCoxT3qFMTDSSg07%)OlkU(D%0fvEMUv*4|5rMS{Fxf#LJU;po_d1LVj6 z7$Rfn3F>$6-G1##i!7nz#;NAy)(1yC<{TI}?id5&)^kC=RNVQBhzn;qS|h!8Tvs(Y zm^Q4`d0Z*ZbhWrkZbi!c|4hOkeZXZ+v_jDk;-iP-5pH}R8FPa7`OVtzoA9-$#_S0+DcqxRP?!R-V5LzqmeA3hmd^M~%f@H>+Gcor8*5ZU` zC6%g(=B?a0g}|D|de%4|CxhUiT#k7jfkF{~;5GlKCSH!jV08?_Cnq5hc;wF{!g1AX zw+wdD_$vg&KJ{1WWoJUcirataqu2)~m>24of5b&x>brxi0w`pgsUtb|x-NhLPD`oH^L$(|Y z;RR2O@IoMYR4)U-YxkW(<{fLD9z6yR2vA}U>5yG9;{H>c|NO8$nXiRA+rM__PydUL zy4zuVwqKN0!O>SfKY2_ww0~_cL*lNsrgm%|r)0YA6C~h)^+Z6GBShkY7&PBIdhN3r z>fyT|QWOaz(P`;Vi8Jt%nrX(lw!Cbrh6+0sBN-UT%Av7l0%erRnQw~4!))GiX&%3K zG4df=0g^6Ee*WTVzbKAZH^HL;$tQfZ=;jf53Em^JNYw=*tUa*dXtEerdiI46KWg0^ zm!rq_|84o@=;#00%;qOKwA2z-xIT_2Un^JW2n4^56WzZwlAmwZ-=JZgcDQ>{3;}aEIt*ZFKUK zlI@+wRNK$?`0=nYCh#6v=dsZPW~FeJJOBz-w?*~v{JFeXPbb7!017!sUR&0T7MgqK z8H;S=8vx@#gmJqHIFqcMIkXuluVow!BEe)7{+3_9`@sX_AN!~D$HlK46j*$TzXj4=O1BmVb|UbX{@v)9o+j zx~*Z5?w4dcB$IFv^hUWWCbjp-9D!>rSjUYKBFGdHW4#cXJl2iJrfH|JPKqnNc(5$w z+0q2MT=%=0f?u$6R+iJ1|Nr88tXy$GNDx0G$pc%S55j%hA;3#_=*CPD*ZmjYIqfEx zuA{JEeVYE_U-&a$>m3cN%$XDBo~~q(M?d~gHRbJt*PqwXOxwBDDdcj~`Y_D)i8q%v z_l1>VZ&b7rD^)BhW!Tq4uM%0(jr~A*%B8i|-g zIFod(fEJ>?pTjEY2yg)jTRrb^@i=_)v#Ouz2dV@@K-k4G)?g6~c#wnshF*VpYp!vpR|HXwc^OB5p#>&`~ z)?mm?xyLMLLP+eS%qnh#VVN1>pgOjk+#QCNksc}1ek~IMnGK~t=3Rr#bCD%1n;KM^3p0!^&81mlbCwBuv1svny#?D{=t3Usbix_pY(?tUihGzE6k@KVI z?Po(LT!%Le%TmsG2Lrh+0UV3gtp_2OR#{6787%oYu491#cW4+}BO=MH=fHp*XnpkK;SO2{aI1jrRf-ffnkHS=|E{On; z`OBd~!J;h>gCxZ7I6(N~Z@aGPf7<+`KXmhJFW*|NJ%24^f^n&WzxLPv+h48_eG)n9 zdVvtEe|t)^8xN0PHV(j%2V|%fIN%P!BcgDG49EiWj7i53nI}b@Avmy}bBPQHi2G>S zsCfqx35koczIs`oQ#Z&>6h*=W zDMOBgW5CG+L(ykjts_wR7@fwP86P-P?NK-r;=cdEcmVpc2LRdQ_vngo^!3-D`c@*AwKFFPz`#3BPltcy@BZwU9Y@P} zYaxledHw6!pI$k&bH69et700Dv=`C2-^ z(z!@v-V4brF~2l@A_GB=a~8c*t-jjI>JAr{Fj&QS0m$KMqsqz?;u&2O1xwH}L0V`I z1jZ$xyaOlu4EiR6ebmS2f+C+nA-V5mzfMP~)YY`pe zU>G{#3~@3NA!sB(&B21VARFnvOv|Y2^4PyHLoJ(pSS>wT*@9L)4Me`Yu=J=B1lZ}+ zoVW`#h=3(uAz<+IJFb_qH)lW8{*4Wac4_CiYw*H3c-OFHK(2XG|HZ%e>em6eCn68r zZ(RGx&$rKi_2-Os$dqrKcOEn+0>Kfj5nfhO0(xXQNc2dCq4&TXX~w7#=s{bimn3S> zJXnpIKs3ObhFaV-=NU?0vpV)hGUF!G)mL0K0dN6O`9v4GXn0m(P`^apU3kQnP8t8C z2p}Ntqecm09}p0{e9K?r?#sK+4g+VoN6EPeEJqF&5Ez<(KJ$Bj?(hB3n?R0{&~xGW zXYP16ir)Dt0~&B)krgOt$wRI{Xtbccu$~OmNX~c^U_2QXdQPag|Wjv}#5K#6SP5^q6~d^dnb)Q#oImC6>=LBt%$i6$(h|vYM@EuEPL{0q4F?Op@c{!Xyp@hw&%s(GiJB>F zv{OPiCWyJQgu-#jxQ|FP5Inr^SGpdfM`iu;f4v0jFUP(Jgex})UHC&+*}&yx`)OSm zAQKQNkIAtC0M$99%Zp#{f69ylB&LCgEH;E#g9ip6EdzS$dw=O~nBV5IS?(Qo&wkz6 zCVcxjYrtuQ0JSqsWUY5bbKuUPmBdH`YKh!52w*vLf{+}5Rgx+xJZV7=ETKUyEk`XB zMgYJ`aLo-vWWWSt9%jQOla(z;Tg%daRq1^q6E2HC1}ma8;)N;T;>`SgFB{{5E33Fe z!S4|u^vWF9J~;o(EAK^1SLX8ZA*XV*o-3mz_r84cyMOhs%zwv0$Sg=?Py8;?Xy@ea zwzuRYT1T$r1V(d_6GO*w^OSLP9<}FATLkNg@j3XAd*uazQw?wjk{Q8-3cLm`omf6d z6(D$9OKB_zbk2o$G_Afet%BhzMM8IEm#y6a#_UJfw z9|8ir@PSR(Z)LAVM;>j!2s}}M+=F$@qpygOw@iso|IXk2v*jJl5j}7gJa)@!%ie#k zHsBythml zN2~(vx^k}`u#bm^>?kbq0Z6qCS#oq9m{ z;dg2G*7;+fdq@BhM1qR3Kr(4;2n+?+!dHmW!K>f?KmLdRa97BhJR*GVBcSuvg99#% zYEjZ)>2E@s65a^U5@`zMIP&J{`O-YSI@0B1&&Ns_X_~SnF2z`6c7m$y#6)(Q~&;9Z5V=;V(r8?@#4*wnS}*; zE*xoYwEzr-X8^2L4NsRQ{>VZ=1l$65H897ycYqBd7)!pcNph%|WZVe@PHE>9A~oI{ z<_$CfRtNR!f4OqI!6OUwvb>=Sbqtq&(A6r9<@H(G<%v+Nj~R2;ve1`ZXY_4@0HU+! ztKrQbUgR-IPTVkW(RnXCdS$)!!f8ofY6|2R{$2LR|JjKr&a4Pt{WAG^c>4<$WSfMo zc@}HSMb#+Bxlud^F4DRVnBk}dcqC6wa|K*^Y?we$Rb;GmhLa%*IFEn?8FWD8C1Hhq`IJarZt+>{n&}_4~^LWG27{vf> zjXE&aLlb&L;=&t1qOidlC3=s6ut-g0JchhMu-00%1{LteLIaI(f}ya?1KCOf_X;id z-q*7KncubO7&*MS(=P;nTyd>dVjy|8Od4I(MoR+p@~T}Ij6exXEaiYUIPpmE^7|I1 z8=vxrq72cS7652OAZsZkDJ>BbLb@Q-B{uBUetu?!7R=YR14hZxr{= zc%jid;IB?UaD2IOUf89?gRn}8T$F}z!AN=p#=1a*tO}>6GlL`u30YbNK<5ZZzwet4 z(koBiRt!mUE)-cHMhu0lNJ=Nk2@BGo(yYeLi~lhEga4*L07*~3Y5lA>Js30;Bh*q- zi(7`4g*Bn)%xN3@Cg2Ho(uLZSbVLY=1!#=yde1Eez?k)lF~b>oEd@!i3N1Ru{1U|z zhEyttKPlf`LTgr8dRnbdxrmG6MPqO=6$e;}2TMT!kZ7Q51wlccH6H7D+7J?Z`JYo5 ze5Uu+CLwQ2WJu@`tT&Bk+!Ru12!yR^;jJ@C<*)qBzxD_JVZmAItmdB5t^M$^!bDmjU-&k>jfIAO%*wt08WjW5A+VG+=$>l%s$( zOvL0ZA}T-({&4awyFB=mT^ar^7D}%8imT@Mg0;JJ1Q9&?jfD4QJSX8<2au<~=Lk>@ z4k+Gy<#cX|d8kz}Ux7~@2M;a~M%u;!hBfnO16N*BzxY@FF%zPe+{`JgE2fQ9ffov;69%Fs z1&SGwXflX8G57`bhFKcFRyVhQsoNeA zwoK)N=PrM>ue`-TC|OV~IRJ*p4W%AHGJ)Lt#MM#+I`zHuul?eG@nzyRXm`$C-Ob

8I*R3j&G4^B|#F=^Ee>DBBpsOz5D}p37P=2^dJ(<$wMk{N%5Ig!G4h z+536$){i@qAma)mCynzatVC_NbV3?a19Mb5Uqu2qGT_VxI%%|wsWvksy)IY;0SDeu zwPdWPm~x5&o-)P*5gQ9MuRIV??RWn86JllZft4ce=>mX)uo46ru3)W8VQj@c17KI9 zP;@cbymqM*cRP`;00Ex=hLe|Hem`LXg!KsMJTwkV$x(Y0%+?}A!&MkYVQi})a&C6z zxj$omRD>3u$qV7V)AjnoS2`&Gk zPqqr+N*-n&dV(BZU|6FIjb-XhM@anDFH`?DwZDcrT0jQ1VazyenC2RE#DGn>WOK3R zQ4=-D8A*Pl{EL70m2W5`Ui!NAM!oSl-k8WGM8@{4OPOy}gc#;>CNifW0tcPB&;XQx zmUdUpe!wz8FFn&J;31b*n+ikvk0$`o`qTet@%Yc3?77<7l2=B1)mz8W z-n7NI4baElxd5pnk@>(ei6IjJdad}6|9NsGo2_>>_;m3^6mTwi6@?J+HuJRzxQYpw z;?f&twMjkasEt%arE6C%WKugv5hosTsRD^Ew$4eRMaUWysxT}GV$eXwYu6yx4w|6t zZ1~EH|JOaXvU~zAqrhx=H*i{+qTNMrOc!+%5pIv99r1R9wJbc=_MQg@gJMok@BGtp z%xtb|82GUJqz-44fD;|70}2fp8De82S2h^O#^s(jc=P=~^^co3Rr1(-T66v%{|M(G z1*>E)(3JSpsnT5d_S! zP^5`RmP89h4R_2s;d2+*c>qY<=GX19l+R)+*Awhn5Bj9&mFxL+Eh7Ni%9YK?3Lr{VXQR^_4sTh*l?I&&coYH-nK-uSb>=)!OIe(H~1KigxWEyG%d zrS3w`Q-WY2qaoLEdoIGFV5Bq2p<&S^R&y0Ktfop4Q@&p`+|@1N8*IBd|IwLVSZa*5 z4lJ{VYfiZ)lG;+27oj0j|E@=f+L*)A4ayiUcNqf!NAiHQtPhsCeYT8EiZ7@NOF^|a zxpg4UzgPBu{3m`Y{o?#fZ?+)Jf{KUyr1xB;{5xb`X0RC?chJt4@&ehy+qAb9P#jk$jQpdV!uX z-*keCbLP7VkoJ;sD-tfPB#Nd2h?J#pQ37(VyzO{mP*bHw9tdk%A6JhCMVtX6Vo1xU z005W+SQ4}SQaA0gh~`RI zMI{+`i@>4OMj8n%S64N66#mF}hJ&sC@Vno9=Z&Ac85BjRrjpNn8HL7&m79A_noK&+ zL*tB=3)zVAdaf{>)uvyGfR%wHqf~8LwVY_q1^=ltra*gZtqzqHi$<`n?f_$9D7|Hj z6~)~{PrzD2Odo+^y(i+Ft)fwE<(vQofNVJ?r7!FnlB|U`Wz)DeBO$S3!Wj?F z{w$mDh0S~|TToRyr;6G|HQ)&h?$-(cX#G#WGzo{;Bh8z%6jHkh0~9ih z2~$ptgD>Z{%pfGoP~^4_95!}7&A=zc!=&k{3i;VNP^;0;6-*9sLIyn*cs^@Js-4R% zW59}*2!&HTkEK~4b(L4Lt12IGMhsYOOu=ThXi6uJd(zgNkSWO{G|CtwoHBG4h&-w? zsjxXGWR7YBX_7C!_3UEnkn@0&d%9o^6JE51n=$!8UInZhs>oNni>XJX%`*2t!ZU7} zQM%?{pC?vmAI=~e^@+JMz{w!XVMCWjL@ETY<}9&w9NcZY-sQkdWV~pe`^fuAbo<=q zVLXZb;aYDvOh@bErg8PzEQp;kWl)Ft*0kG#2@mI$P2^f}UbN*gvfjd6sn_2xI-Rqm z@ghb8u}kJofLq8ItW#1~E=WxnCyu+)^QnfBw!3RB()0d+x3cFy5X`S&OL#PeT^B&QOcR4Dy5p zbLnudvD|9SMy`EUaD7l!+#fS$Cy^MCHz9f6IwiM=5@g0v+W3Z=S)^s7^DHfwl)y4V zyma6nK~Yyx079@1k>PT?pjSL+##b6@R`9vwKO`p=SDw9guytHTCK}A`Wzw%Y0dGO9lT`6ZZIAG_Y z+&<{^y-ztmlQ!f{TvTN#(_-qQ%~Nb-Y3d1;DDoM_UMJING`V95cEk*9wP;#w8b`$? z+dA_JmX51wh*_KG)_A#ERO5@5FdaHvl7%0+qnAnf6(qvDbH6VWrt5 zc-sR2fAW9#H}laN4oWu0Mo#oxY|R=*Y&M0R`O)`%>#`HmR&M5#qY!!ERC_TlR6YsF zGlYOi#Y9G*`PPEU&GYxKn7%DGiuEZJ3kft1l6ll%XU;-$ohpigtTcJYZ+e5p;5Z=f z`oVhZBW}wRnE3VcsW=doid41jqNz2qFwGI6=%4oItnu=+7hAt`>MMqV993n;Azho8 z)hP6$NA~cEFuL3ji!k`iXxN>29B+?w7#yFZI0gVn2M0bH-x7nd>>ph^itbKxvzjNu#DC-JON$et20};pUNJ=en#-a@-1I z;)c(Wa#^)_+l&)?L~~2T1bA16OsmxDe3c&fWLEHrQF8TYTz%H%YB;E>ia1avD71hG zkCz(M2NM{s`MdIA#N)Q`9e;AI9XBRegqiNLozR}N(t>)fw8>>$@I373K?60KNhH>0 zdT^?yy&w6l8mdbN*gg|sGxK2(HonBVPV;eBo2J%}JsiyET_4~(Pb4#xM)sx=EOa{F zHpdD#j%Tx0H6J8w!M~vVB7Fe$I#=C}SZk%b*HlC-3zQI@06Wk&Z!QhPOz_xqQK7Vi z7x0h^&hYXxgXKqdb9XAkbU@hit!~yz)A#QYGW-etW%W{OtE4B}S-M%MfGL-?TJ2Wq zLfuUljkQ!tw^j0TJ!*n8w9I##*v_c$FQ<}Ejn$oP{;N;g7eo-7r!K^<2q zH6#HO2I!WIBeL2=votP)EYIB!(q-WXdLaM+6oZ)0>NV;mjl^P+sNfy}5PSYRE|rsx zRlUY^!X}x%2qSH+HSfMOAGDciYi?N)<;RIPFm?fJf{F2?c#)m9whrS-(GBB8s53iy z{k?|f@t(i5kZvG(DLk!(2ieRMnS2oE)1b@;$|M<-P#Onzj_1PkW@J{ba6EWF=)ZYvCO2w~!O_#Do%p&D*kZGIWpDOkb>;`T^EB zc6K0$s|7n!)A#R}jQrMrssTa#FaP(!#cE$!&Pq*(*z^6^~bk){t4|v3NB#rElZnv|< zPOBiu)*ysz%|x56R>U@O;dH91B+C}QT}8d2N@NS9;{SEjxX8|+j2eYzk%ebYimFT@ z=etdrNA2VD_mFw?z`10MGZP8jSRo$i<#3T9cw{AuWe<<{r~KB3M^joK9Ahd#F_Pt~@N*M&7oDo7p?)r5s)l<3eXnwM?1&xO5o?3Vkj5=Rlqi%s6f&7dz*}+3BTXc@r(8JQ}^Dzr1(so!JNHrUrhh|B(mIT7VHQVXqRP*)e^ z!qqGAVqzLa&`Yn3W!`Peoh>uteKoXzNZX4E?Ibj>QZ+Z*Fb9(PdFE~;;(od|m^7VK z}e_ZLC{_3%pD1>7ioPho#T9KeCM+tEY?&oZ~5z^ia%t! zjU_9IO?*0YYN&#&4i9|mr0Bq=K6c)wBYmxGGS|wpbR;{1zZj-A%@bWzjpat$$_$gP z*K>e;9@o>U8J>K^kbwuzWyn1$C}k8X;Z^MI$~3cj6aWAjbQKq$N?TPUHsVMuAflY! z+9yQ#(T#uSeyV|7P;g`hezA+t*G;0A`uQ zs0ErT{6QlT005}Imz_pC>WZDH%7XecPj&8!Dlcb%U#5L$Fp+NGwAG3 z+w2TvGmWyvoUe~6t<<-f?=&JvBgvFC4+ypv!)rc;K^66TNg-QlbnZ{k+kS(A>eSwD zW@kax%Uy4t*gyswG>g(p$dP;BZS!u(MW??|duy#uW#p)g3nmv#T2TSIqKP3Ycv%81 zx4e4t9D$MMne&dj$3DJr5ie)s(WOvWjuzcbmDJ~fv?7lRnW@YB1OOmfxBOG*+eg#} z)(Rg#ibQrasp?OiOAVRzan70t{S%oAT@UgzlN8fytjudw4kNF2{o11S{8~Q2xU6a% z%jaJG-P%9Mjz`;U9v~hC5y^46{?zaa6)+3Z*d8X*&C^7Oa%ZnldADUA0^f1+5VP3= zu?{0sm+V~eMQ1HqfDc(xLA!rrXQSaNGSv`O?yY?20AUH8cI1$_t#gtbZ13?LKKVh6 zBT|ebgQmPZU8IAWXW`~komsnaA6@;!e>(x$ocd>XpM0vzsxeeP%j0A)G`r<-EA5|> zoycl*V<%%esOBdk7|H%_lyn}y?pf9da5|U!Q6cl_-b~sD#P1#51hDE){&)4mW+o4- z@z`j`&qkZRYklW(F)3Q>z3v7zo1G|+is))@bSP(Yq0Y)Bl4J6}uc|Xf8>{bQovkv)MowC)Sago@VKbK z=i+RA;%5aGb`;&81X?GL->O(r6rI}l^X;J;@U?Mwancjs2^-eYk+5~!?m;sTdv`61 z$HOgqz>0)d%mQ)C5P`rmK;+zMC+jD=wIq?j%rDW?#K{@u$p>6$r~EKqGhvykT2{sR|UW3WBel=qk@^XO7N%%s>n7TcT2sMs+3Yt?06a29nNhVINojy({5oMd3* zom+UCZ~1#&**or260^=TfAsMZtH<|5DL%S9=8t{Y@&`d!L`_!iilAq0(jXhyIGUtS z3?KS_vQ=GqG!*RHmb_8D=#7T#NmRD6XHQvD_BH!DjAdpp!!Tp)MY2?i!th3sHS1s| zBq6fIFtUtoY=f~)X6!z{?>pyr&iB{zoO55-b3f<)>v`_`T=)GuQ)2@nk;<)$7Ark= zh4s6EQ~*r^lMgvqSXiz+HP+X&DcM9!#7n28OsDm{oW+ZP|1KMwI%?u*+{}E@9i!w^ z0qecOz_Ku(Dcq35z?%tOJ-38&Yw&s>A7J>VQZM~?aSh>;pZfI+wWl2lwhNj& zcSV1LRVTBr;<6F{B1)XZu@VPi1V~qrPtdUOv~2i~gi!t5Vj=CvbN^C zfj{qroH4o55pqal`qJ&AzzU6^_o9-`0Qk-c280$5X$YNw0s;E$MqM;n|YX zkL(RA4E|$;(@3EVU&qnT)zoAXLNyR|d|0xg?9Xjz^et&;`l!NQG>?32w}hzBT?BV3 zm`av5N*UH54_Q{$yy*JDIs>-B3C&JqzXS`JwBF__q@zeW zp>8wB`HD;H+lM+x|05KUGlAnV3+EXxtt74cV~=t72M~l}ahnIfJw5WX&N)k7XKjw? zNdAKIEE0uR*x6ypVDAm3&+0d12=MC|`Q^^P-mDU7KKzp}WJ>ilQT#5Qg(UHU^%7)` z6FFiR2u10^;8PW#a1!*Lw0+_Y!OCrB^fiUthT0lLC76uiT z{7&+~fB*Z=!G-4=g?w*tx02QZF?LF-kI?Ph!;Gxm&&YYkagqRM)^go4gUj*T)xW_q z*J=0LzSZswIS)N#g<7D{`$%UTE^mx2KY>*X!e0+mjGh@bC zUyosUosGk+lQm1hb2j%Z%CfHbz)0NPVKF7tIo#;(y$}xQ2gp}pzlJhN`GGTWiK2<4 zNQG-YErS6Bca$GPe{%IXxZSiWn?HhURFlJ!^X8vkjmVH6tKQG&tS6R3XnEDOArk$r z8nW!i2zsHnKXN_&>%7v5+XljDAUED01(LKBAii1#^xb}4H%2=tZeGbpkYR67yz@c?v$k`s4(=+}8 zu#Ccc-eHy*+Q5^GYE7a~MGm)a64?JC#(mL(9DSBvT$c3y+@Dv<@)^2JDL~c(f|sg3DT`mS&uhcJck%%flRl zwx0xL4JaXup-%Xhy>v?E%|`kq2Ha;nB65%UWk%;kY(usfpgg2^@Xb$K$xhL8|6H(V zjYEEIW}ms|?6>#RR`;z78HWEDyVdxG=UQSbuliVV!v{c`R2| z`n_UhSXtw)Iz=|@-l0XW^Z=TXs|%`38~SjEb+5&!)51r30Y>8T-}F>^e{&}0t7fcm z>fb~_{ZJ!SXs#|Vt?m_ezDVd#L~$kr#*S;%1uSnXYpaAGZSN}cW!w~IL_!Dne8t;# z4u&G}K=_It+0%ndiJZr*SqwC~dNuRC9gmg)rpL`zXo0b{Cb0t9bx80=zacpHM~}ypXRJbL z2khiBwa;N&(v1_tV+Co?tOJ6n67To~dB#o}M+{WAMU}GOZ(M=%JS^m%2-gtJiPDCM zb{a>M8lxOM4tv=uvnXMU@3hEzW9J{|(5m<0vlHklZBC8af*6wZy#s7Hze%xSsQh%s zn=d#Q^lzGxl!UHMo_ztw*o)8DoB9oxsJPnaxI@l5@Bs&0;(rfTzxH^D{*a!?vs#|gc% zO?h{u;yg@GTpV&(+A1x2xmN#-4Qh20%C+i*V`B3b^J>n2^TWEMPW9JOd4N1s_fCKo7T?)2ej zvGxAM89#>-ld40V)cEt$Tmr{JvTV6DR~>1Q#g*@Mq(hCh5603Px2lGn>^sze9N4xO zaG2n>1Ed^nivJY&OwFo>mu-2_<&x2vO0l(}7VLd#zv!6ak=-yNCX}Y&aE#jcF&?4Z zB*PbCj&B?Ms2R`ZsHgJtz?kzNd1B;Ax8Kf8g!JxW?%y&Z$O&rf7}73G* zT-qOUa4HIY)Jxhm-$}C9ob#?&hD@Fxzm`O+AIXAOSc%qb0)+X#Qv4 zS>H~g+zr|kVYw10RC6_Onmc*VX7QoC$f*J&(E0N8iQKSnyUZeic-FAJBK^L~^n;zi zhKIhsf>6lh-zxf%knw~uQuadOaqqc>v(AM2u&7Albb#Q_3m1(kLEm;g7-xE(5N>MP z{~NI>-|FPnK}%zEg#uVMyWN+4$16lgt8y{1@&xc;i>rMguQPwDv*QoxYz`Lg-YOAfmo#Qy>!X{o z#F!{v`QQLRJu2^quJgM35H}@c9e`zq#Yfm;`OkaPcqA5Y^}3kG_;cvvD{trkkj;BT zORSY%Wb|9AakxT}{@A2tGjE0E@LDt<&X(We`KpZ&!~?AxT@W97)r4{?!W5PKnCUaD z@p}Iv4NVbb&mpaVGb{TfF}bo_PTq8fJISQ}QHacT^EK`BJUcG&D=Dn(vtL9-==rKk zFM+f#Ii`dFqO#VlW&W}E4JnCt))#C;oe~~&T3&k3AnouIbxFJTQhocHi85@RA1{3f zFR$qxbwB;~bU7`=yVllb-k3Gp3nF3mxBGra_bv2|Ku698{hBb@pz2+$-Mq`1Fo-Kb zTU&-mbrwC?mx}VegxfAudDT&mn=iD(@u#<=lqWKqtA0;hc)AD;{~%t*bNm*QY} zR~(%9@_Fr)3KqsoX(B&sjlQm0xj@*LxtxKWawnTbmYy@&s9LY4tU)l@*aLPx}>c=t%7o82jA{<*o`FELe}mDC}N{nb@G%#(z_om!n%QuPr& zjq*SGtz6>Qr;#XdIF@IwLpW}fX$_g86kz5^AVw@zOTlb@a*Pk8idW4wm)hw()xAo} zkByatM-IH=GQ<=WS$x{U`W%`4u6H`sQyS4a6SF*J+Qa%UKm34_b066kwu^Sf412y0 zRL;>iF?t?$SduJ}f8%)w1Sgx%PdFU%s~YVi)%T*qYMp;(Jem^S7Hyj1M15ajaOh?m ztwf#2urR+zFm_*_=W1#??*E5x61@ZxOW47KR!8OJd-(uxtgy~-MUqpHOmh1#CCSa# z=en`(KWlLUkJs`RN1}wHqG-2uj~qFA>M`ueu@?q%T0BE=Ev-6N3GG;F_W6Y%;Gf8E z<7B09}w@TQ3;0cU7Q)1a4frfXy`kosvPBrIx(+ z#Sw)y)q3;urXT!lJBn9EhhhYw`}OHEJIKa9gDtV*utyI@+@DPO4gZEUJAh(_HReA|Bz|j=Y!P0I9LcO-b!LKK3-&A4B@O!*BJjR`i0)bwWki|sjO_U7koymE~0 zaE^KNBh*1ji-PUmk;BQ55&?6fi5FB?kM$U8Ck9AWYAV_G$28hUS&n8oVHcZbKBQAZ zEV>FwMBuVR~JGe03x1|{I9I?D7zL3tRzjP#RclR78JM^RokM{Q}GA6LQH8MZr z{xo_K1~*3Nlr}$m@ac@iAmjsvHvB7!^H#queJX=~A;-EspNpCo=#ZwpJfYYC8>8QR zFs1bI^S+*AjZr`!H9)%lGm|7uzD z(NMq|Hv7PqF~_JMj=u2&-QX*#=PgO7eE7CSkuI^X$V|s}7QuIcmr5=iK7fs_j4$lc zjp8XqS6-`=?L8^sWGUZHdXUw2NoZ`+%}gQOz)y#4VYq{d_rBI-|6%j*PTz^6Ki}nU zqplLr-4%5|34McO`0?$dzA5N>ed>U24;f4s^!bNgoIr{SD$|ulHIxtDCa#x^fm=l% zMBjItGLBb6d_6l%y8N(BDE6}>cYkW)qjc_$12(pA9l`uXrS&8S_@8*BZcf05>W^XYCDjADE$<>!PkJG^N?Xv!2A@iNw=I zxt6PG-x-zba&zc<$ky2W(7x9A4v3=96tR6?MZ2duzZ1#(sgK)ZYk!<-8K$8?FjCcx z_Hi;#&r#AT=Hmj;odl({WM3T^JJ(|-^J%!=^>(?AvAGpd4B|fxY8ih_@-@J^uX~$b&{d{blvLyP% zTwCcK*?IxWq%@lxS)2$6Vt7o!?|CgBSMR_#n_rl(?H3Jx8t6_k1Bh0*C9C4fI>ZG` zH;#TVJc7I{bHg^wyVP8SdAal2U_s78T>%-lvt2%)!xG9dVJ-%e`mN|Yl;+}6CTZXm zy6YbkUmi>}ZI{}CU28)ZWEsttw}#h%vwToq;YCQ){Z5C**|3WF{FQ8_gUJt718X-% znRq`UJVH5I?=s&6fqMZTeeueZ7(?_W4nts^EVm%$czyxHICJL3+*@EkGt<+Ra&gqW z)*y1{U~fe4clJv9z<500zUJ;A+TK_O#XV_Jn@wU29O*^%7X&+l$H9mF;=+cDu#c96(#WAwBwlmupz5DX)-L%}UyAm(c1;SOy z))Tnl%D)aQm)7#SD``{9?V~z*PZ&Vfs8Db*dv`ta?1k%UBUEDmayp{C5mHRqS@t0G z8sCt@8b9e`%G*=bzh})AN;`qBME{z4w_$>YVO&!RXaBrG2}mGKHiPmjwJf?xE)=M9 zN#bV#-0{q%Jgpv5@8K1aan5OFBcmf?fKVk2b8&45vJyU*0cg!U%?^2Pp3w78`1Gpb z17X3Si>E2saD>5Qgcl3~0YMQgHef%vm>0d&j>LaFi6bQ z0|B;&g#^H5#EhX{a$-n71o%P_C;{Wp!0z<-OPvau{p9luAC@Esit_7STSd0zK^egYV G$Ndju&6<1w literal 0 HcmV?d00001 diff --git a/Turbulence_Sim_v1_python/data/chart.png b/Turbulence_Sim_v1_python/data/chart.png new file mode 100644 index 0000000000000000000000000000000000000000..47244dceb3553c797aef5071c0d80f37470d467f GIT binary patch literal 974 zcmeAS@N?(olHy`uVBq!ia0vp^4?&oN8Ax*GDtQ2@>Hwb**ZTT;AoKtK|4mIzd3kwi zYHFUIo}rej7W5fKsc^78A}ty{8WiMY7<-o1N;g@unDJ2qj$1S2CO8ylOM zGiSbg_ipRft^F%scQG(9^Le^BhE&XXd#5{Zk%NRoqRh@_jZhs0ttkE4|LQU8V{{B| z-(@-Tq{dVE&os^DX1Oc2&SC>P86CVobv5;m$Fb?_F7j{JG}y@7c_#Y!-D`1b(M!Z< zn&wq;JXes~(t0D1>`$r&umk8}R6{MWJd^Bw687kG~@D*HE0 zKz{4@4XqJ3t~S2sSYA=%rmlD_wMJo%yRTkL-Thx8Nwx_KB|CiI8qU$5_dIw`>5h}D znjgm<*`4HHK4-;^n|`i!H6^X`{rcyw-t=VCk&z8|t4)a3Z|?kB#JAOX+p@@KeLSxP zEly2)R=mSxyR+E0&O>G*K{w>NZU=ZI)@$TAE4Ka(+*%;?X6^53t6GB<_q~s{6{`z- z-+Nm7L7|~f*UcxT3Cpkl{k&?|l;WFREgbsh-H#sqJ1F%0{5|gDSF4wAZd+O>vMcrB zqiZRl(|7%QbeKDTVWU9Z-^DY}$208-O_codzcp2>RU~LdLE@*X6V{@;?tLk-GTm!y zZg$E4p#S@j71!%5%x7&~v)u3DduGig6P&n0wdcJ5zU)?O{N3}6+5s=-&t0dTz{-`r z@7`CL4NR@7Twp~ILS|FJ@`{Hkdm~rHR%muy)$4It@k{hW(ae>x5sQE(-QW7b`>WIS zxKlDI`la9a-^ZMY*9!1Z1iGeTm!D;cPt#a3;ZWBt;nTC19oqk|pZV$kkR2lP z{x7^>8g*~C$kOh-%r8|P@s=vG<~mcWj-K&|XKVjjI^p)!rJr8^WZyp9!guX`5w&@$ z!m37r&lLkN@A