From 31d4dd5f24f48e4dcade65051b67925b6cac0080 Mon Sep 17 00:00:00 2001 From: Thomas Dehaeze Date: Tue, 5 Nov 2024 23:34:34 +0100 Subject: [PATCH] Start to modify the way disturbances are configured --- matlab/src/initializeDisturbances.m | 302 ++++++++++++++++------------ matlab/ustation_simscape.slx | Bin 160671 -> 160679 bytes simscape-micro-station.org | 299 ++++++++++++++------------- 3 files changed, 338 insertions(+), 263 deletions(-) diff --git a/matlab/src/initializeDisturbances.m b/matlab/src/initializeDisturbances.m index 3e7cbc1..37a4add 100644 --- a/matlab/src/initializeDisturbances.m +++ b/matlab/src/initializeDisturbances.m @@ -1,183 +1,229 @@ - function [] = initializeDisturbances(args) - % initializeDisturbances - Initialize the disturbances - % - % Syntax: [] = initializeDisturbances(args) - % - % Inputs: - % - args - +function [] = initializeDisturbances(args) +% initializeDisturbances - Initialize the disturbances +% +% Syntax: [] = initializeDisturbances(args) +% +% Inputs: +% - args - - arguments - % Global parameter to enable or disable the disturbances - args.enable logical {mustBeNumericOrLogical} = true - % Ground Motion - X direction - args.Dwx logical {mustBeNumericOrLogical} = true - % Ground Motion - Y direction - args.Dwy logical {mustBeNumericOrLogical} = true - % Ground Motion - Z direction - args.Dwz logical {mustBeNumericOrLogical} = true - % Translation Stage - X direction - args.Fty_x logical {mustBeNumericOrLogical} = true - % Translation Stage - Z direction - args.Fty_z logical {mustBeNumericOrLogical} = true - % Spindle - X direction - args.Frz_x logical {mustBeNumericOrLogical} = true - % Spindle - Y direction - args.Frz_y logical {mustBeNumericOrLogical} = true - % Spindle - Z direction - args.Frz_z logical {mustBeNumericOrLogical} = true - end +arguments + % Global parameter to enable or disable the disturbances + args.enable logical {mustBeNumericOrLogical} = true + % Ground Motion - X direction + args.Dwx logical {mustBeNumericOrLogical} = true + % Ground Motion - Y direction + args.Dwy logical {mustBeNumericOrLogical} = true + % Ground Motion - Z direction + args.Dwz logical {mustBeNumericOrLogical} = true + % Translation Stage - X direction + args.Fty_x logical {mustBeNumericOrLogical} = true + % Translation Stage - Z direction + args.Fty_z logical {mustBeNumericOrLogical} = true + % Spindle - X direction + args.Frz_x logical {mustBeNumericOrLogical} = true + % Spindle - Y direction + args.Frz_y logical {mustBeNumericOrLogical} = true + % Spindle - Z direction + args.Frz_z logical {mustBeNumericOrLogical} = true +end - load('dist_psd.mat', 'dist_f'); +% Initialization of random numbers +rng("shuffle"); - dist_f.f = dist_f.f(2:end); - dist_f.psd_gm = dist_f.psd_gm(2:end); - dist_f.psd_ty = dist_f.psd_ty(2:end); - dist_f.psd_rz = dist_f.psd_rz(2:end); +%% Ground Motion +load('dist_psd.mat', 'dist_f'); - Fs = 2*dist_f.f(end); % Sampling Frequency of data is twice the maximum frequency of the PSD vector [Hz] - N = 2*length(dist_f.f); % Number of Samples match the one of the wanted PSD - T0 = N/Fs; % Signal Duration [s] - df = 1/T0; % Frequency resolution of the DFT [Hz] - % Also equal to (dist_f.f(2)-dist_f.f(1)) - t = linspace(0, T0, N+1)'; % Time Vector [s] - Ts = 1/Fs; % Sampling Time [s] +% Frequency Data +Dw.f = dist_f.f(2:end); +Dw.psd_x = dist_f.psd_gm(2:end); +Dw.psd_y = dist_f.psd_gm(2:end); +Dw.psd_z = dist_f.psd_gm(2:end); - phi = dist_f.psd_gm; - C = zeros(N/2,1); - for i = 1:N/2 - C(i) = sqrt(phi(i)*df); - end +% Time data +Fs = 2*Dw.f(end); % Sampling Frequency of data is twice the maximum frequency of the PSD vector [Hz] +N = 2*length(Dw.f); % Number of Samples match the one of the wanted PSD +T0 = N/Fs; % Signal Duration [s] +Dw.t = linspace(0, T0, N+1)'; % Time Vector [s] - if args.Dwx && args.enable - rng(111); +C = zeros(N/2,1); +for i = 1:N/2 + C(i) = sqrt(Dw.psd_x(i)/T0); +end + +if args.Dwx && args.enable theta = 2*pi*rand(N/2,1); % Generate random phase [rad] Cx = [0 ; C.*complex(cos(theta),sin(theta))]; Cx = [Cx; flipud(conj(Cx(2:end)))];; - Dwx = N/sqrt(2)*ifft(Cx); % Ground Motion - x direction [m] - else - Dwx = zeros(length(t), 1); - end + Dw.x = N/sqrt(2)*ifft(Cx); % Ground Motion - x direction [m] +else + Dw.x = zeros(length(Dw.t), 1); +end - if args.Dwy && args.enable - rng(112); +if args.Dwy && args.enable theta = 2*pi*rand(N/2,1); % Generate random phase [rad] Cx = [0 ; C.*complex(cos(theta),sin(theta))]; Cx = [Cx; flipud(conj(Cx(2:end)))];; - Dwy = N/sqrt(2)*ifft(Cx); % Ground Motion - y direction [m] - else - Dwy = zeros(length(t), 1); - end + Dw.y = N/sqrt(2)*ifft(Cx); % Ground Motion - y direction [m] +else + Dw.y = zeros(length(Dw.t), 1); +end - if args.Dwy && args.enable - rng(113); +if args.Dwy && args.enable theta = 2*pi*rand(N/2,1); % Generate random phase [rad] Cx = [0 ; C.*complex(cos(theta),sin(theta))]; Cx = [Cx; flipud(conj(Cx(2:end)))];; - Dwz = N/sqrt(2)*ifft(Cx); % Ground Motion - z direction [m] - else - Dwz = zeros(length(t), 1); - end + Dw.z = N/sqrt(2)*ifft(Cx); % Ground Motion - z direction [m] +else + Dw.z = zeros(length(Dw.t), 1); +end - if args.Fty_x && args.enable - phi = dist_f.psd_ty; % TODO - we take here the vertical direction which is wrong but approximate +load('dist_psd.mat', 'dist_f'); +dist_f.f = dist_f.f(2:end); +dist_f.psd_gm = dist_f.psd_gm(2:end); +dist_f.psd_ty = dist_f.psd_ty(2:end); +dist_f.psd_rz = dist_f.psd_rz(2:end); + +%% Translation Stage +load('dist_psd.mat', 'dist_f'); + +% Frequency Data +Ty.f = dist_f.f(2:end); +Ty.psd_x = dist_f.psd_ty(2:end); % TODO - we take here the vertical direction which is wrong but approximate +Ty.psd_z = dist_f.psd_ty(2:end); + +% Time data +Fs = 2*Ty.f(end); % Sampling Frequency of data is twice the maximum frequency of the PSD vector [Hz] +N = 2*length(Ty.f); % Number of Samples match the one of the wanted PSD +T0 = N/Fs; % Signal Duration [s] +Ty.t = linspace(0, T0, N+1)'; % Time Vector [s] + +C = zeros(N/2,1); +for i = 1:N/2 + C(i) = sqrt(Ty.psd_x(i)/T0); +end + +% Translation Stage - X +if args.Fty_x && args.enable + phi = Ty.psd_x; C = zeros(N/2,1); for i = 1:N/2 - C(i) = sqrt(phi(i)*df); + C(i) = sqrt(phi(i)/T0); end - rng(121); theta = 2*pi*rand(N/2,1); % Generate random phase [rad] Cx = [0 ; C.*complex(cos(theta),sin(theta))]; Cx = [Cx; flipud(conj(Cx(2:end)))];; u = N/sqrt(2)*ifft(Cx); % Disturbance Force Ty x [N] - Fty_x = u; - else - Fty_x = zeros(length(t), 1); - end + Ty.x = u; +else + Ty.x = zeros(length(Ty.t), 1); +end - if args.Fty_z && args.enable - phi = dist_f.psd_ty; +% Translation Stage - Z +if args.Fty_z && args.enable + phi = Ty.psd_z; C = zeros(N/2,1); for i = 1:N/2 - C(i) = sqrt(phi(i)*df); + C(i) = sqrt(phi(i)/T0); end - rng(122); theta = 2*pi*rand(N/2,1); % Generate random phase [rad] Cx = [0 ; C.*complex(cos(theta),sin(theta))]; Cx = [Cx; flipud(conj(Cx(2:end)))];; u = N/sqrt(2)*ifft(Cx); % Disturbance Force Ty z [N] - Fty_z = u; - else - Fty_z = zeros(length(t), 1); - end + Ty.z = u; +else + Ty.z = zeros(length(Ty.t), 1); +end - % if args.Frz_x && args.enable - % phi = dist_f.psd_rz; - % C = zeros(N/2,1); - % for i = 1:N/2 - % C(i) = sqrt(phi(i)*df); - % end - % rng(131); - % theta = 2*pi*rand(N/2,1); % Generate random phase [rad] - % Cx = [0 ; C.*complex(cos(theta),sin(theta))]; - % Cx = [Cx; flipud(conj(Cx(2:end)))];; - % u = N/sqrt(2)*ifft(Cx); % Disturbance Force Rz z [N] - % Frz_x = u; - % else - Frz_x = zeros(length(t), 1); - % end +%% Translation Stage +load('dist_psd.mat', 'dist_f'); - % if args.Frz_y && args.enable - % phi = dist_f.psd_rz; - % C = zeros(N/2,1); - % for i = 1:N/2 - % C(i) = sqrt(phi(i)*df); - % end - % rng(131); - % theta = 2*pi*rand(N/2,1); % Generate random phase [rad] - % Cx = [0 ; C.*complex(cos(theta),sin(theta))]; - % Cx = [Cx; flipud(conj(Cx(2:end)))];; - % u = N/sqrt(2)*ifft(Cx); % Disturbance Force Rz z [N] - % Frz_z = u; - % else - Frz_y = zeros(length(t), 1); - % end +% Frequency Data +Rz.f = dist_f.f(2:end); +Rz.psd_x = dist_f.psd_rz(2:end); % TODO - we take here the vertical direction which is wrong but approximate +Rz.psd_y = dist_f.psd_rz(2:end); % TODO - we take here the vertical direction which is wrong but approximate +Rz.psd_z = dist_f.psd_rz(2:end); - if args.Frz_z && args.enable - phi = dist_f.psd_rz; +% Time data +Fs = 2*Rz.f(end); % Sampling Frequency of data is twice the maximum frequency of the PSD vector [Hz] +N = 2*length(Rz.f); % Number of Samples match the one of the wanted PSD +T0 = N/Fs; % Signal Duration [s] +Rz.t = linspace(0, T0, N+1)'; % Time Vector [s] + +C = zeros(N/2,1); +for i = 1:N/2 + C(i) = sqrt(Rz.psd_x(i)/T0); +end + +% Translation Stage - X +if args.Frz_x && args.enable + phi = Rz.psd_x; C = zeros(N/2,1); for i = 1:N/2 - C(i) = sqrt(phi(i)*df); + C(i) = sqrt(phi(i)/T0); + end + theta = 2*pi*rand(N/2,1); % Generate random phase [rad] + Cx = [0 ; C.*complex(cos(theta),sin(theta))]; + Cx = [Cx; flipud(conj(Cx(2:end)))];; + u = N/sqrt(2)*ifft(Cx); % Disturbance Force Rz x [N] + Rz.x = u; +else + Rz.x = zeros(length(Rz.t), 1); +end + +% Translation Stage - Y +if args.Frz_y && args.enable + phi = Rz.psd_y; + C = zeros(N/2,1); + for i = 1:N/2 + C(i) = sqrt(phi(i)/T0); + end + theta = 2*pi*rand(N/2,1); % Generate random phase [rad] + Cx = [0 ; C.*complex(cos(theta),sin(theta))]; + Cx = [Cx; flipud(conj(Cx(2:end)))];; + u = N/sqrt(2)*ifft(Cx); % Disturbance Force Rz y [N] + Rz.y = u; +else + Rz.y = zeros(length(Rz.t), 1); +end + +% Translation Stage - Z +if args.Frz_z && args.enable + phi = Rz.psd_z; + C = zeros(N/2,1); + for i = 1:N/2 + C(i) = sqrt(phi(i)/T0); end - rng(131); theta = 2*pi*rand(N/2,1); % Generate random phase [rad] Cx = [0 ; C.*complex(cos(theta),sin(theta))]; Cx = [Cx; flipud(conj(Cx(2:end)))];; u = N/sqrt(2)*ifft(Cx); % Disturbance Force Rz z [N] - Frz_z = u; - else - Frz_z = zeros(length(t), 1); - end + Rz.z = u; +else + Rz.z = zeros(length(Rz.t), 1); +end - u = zeros(length(t), 6); - Fd = u; +u = zeros(100, 6); +Fd = u; - Dwx = Dwx - Dwx(1); - Dwy = Dwy - Dwy(1); - Dwz = Dwz - Dwz(1); - Fty_x = Fty_x - Fty_x(1); - Fty_z = Fty_z - Fty_z(1); - Frz_z = Frz_z - Frz_z(1); +Dw.x = Dw.x - Dw.x(1); +Dw.y = Dw.y - Dw.y(1); +Dw.z = Dw.z - Dw.z(1); +Ty.x = Ty.x - Ty.x(1); +Ty.z = Ty.z - Ty.z(1); +Rz.x = Rz.x - Rz.x(1); +Rz.y = Rz.y - Rz.y(1); +Rz.z = Rz.z - Rz.z(1); if exist('./mat', 'dir') if exist('./mat/nass_disturbances.mat', 'file') - save('mat/nass_disturbances.mat', 'Dwx', 'Dwy', 'Dwz', 'Fty_x', 'Fty_z', 'Frz_x', 'Frz_y', 'Frz_z', 'Fd', 'Ts', 't', 'args', '-append'); + save('mat/nass_disturbances.mat', 'Dw', 'Ty', 'Rz', 'Fd', 'args', '-append'); else - save('mat/nass_disturbances.mat', 'Dwx', 'Dwy', 'Dwz', 'Fty_x', 'Fty_z', 'Frz_x', 'Frz_y', 'Frz_z', 'Fd', 'Ts', 't', 'args'); + save('mat/nass_disturbances.mat', 'Dw', 'Ty', 'Rz', 'Fd', 'args'); end elseif exist('./matlab', 'dir') if exist('./matlab/mat/nass_disturbances.mat', 'file') - save('matlab/mat/nass_disturbances.mat', 'Dwx', 'Dwy', 'Dwz', 'Fty_x', 'Fty_z', 'Frz_x', 'Frz_y', 'Frz_z', 'Fd', 'Ts', 't', 'args', '-append'); + save('matlab/mat/nass_disturbances.mat', 'Dw', 'Ty', 'Rz', 'Fd', 'args', '-append'); else - save('matlab/mat/nass_disturbances.mat', 'Dwx', 'Dwy', 'Dwz', 'Fty_x', 'Fty_z', 'Frz_x', 'Frz_y', 'Frz_z', 'Fd', 'Ts', 't', 'args'); + save('matlab/mat/nass_disturbances.mat', 'Dw', 'Ty', 'Rz', 'Fd', 'args'); end end diff --git a/matlab/ustation_simscape.slx b/matlab/ustation_simscape.slx index bbabff0790b02d1cf5a4ebbc2d82b939db833461..eb5a532c1b0c359193ec3c558232d916c48387c2 100644 GIT binary patch delta 20940 zcmV)9K*hhG=Lx6h39x({1?aEfiho&E`eO}OZ}i18+JjZ@D%@5E^2?P2 zbvPW7!?@)_{uIBhAD-K<#-(n0B>=He6tykl9=66FnKBZp1wNM=qzgcWg&@l)JVQ_M z`^xM1?!H(dmXUNGlS!PCg3M`3D7j&IU*L=To#*LjEGs1JIDA6`4q_Lr;?A$bpf@GR z+avQk6K2&Wyr`V#RGs{9&@7}(V=^tWjAk=BnVkjg4f+%o_GPnsY?@4`?dyZm5#4gK zr(m4%7U`wYsiA+?kd1x-vlAW9Uj;M%CxhRU=X@rAM{nCO6u#$I7`oHqCRNx9(5B;1 z6h+24;EaV>10vlf{q?16%Oi1{B52WzL_U6e<9ma%^M`sl`R4jiwaq+3QAcWv1fpaOVAsm6st+|#XJNq-ft z9vidh+ot;AY#IBNo8?c0BeqkGTfX4R*Bxf`M+xJgmpHLwETD=DPzcWe^hI(=AR?^K z@@=^JLA+a0-?s5uh5J)B2l;3RVWK=rcC2Te z8$^i^{y=F`t5tA=hHZ82wuhL1m5)Dv%E_h^Y>86o+hbw4uIhun*z67w>jVxi9_R8n zpAltdH`0bJnJ(Zvw4No~Yop|Mejd1;OqPl{F3w5f@PL9%h~ zKyZvS?XvAl%-Ag7ERT;DD!g5`=00pW*zb>VLnK$7GrtOEID6vcz>rI<*b+`+2ZA>1 zW3~CH_lZ2js!Vr=-`P`mcm>7YK}IkKMrjJDOY>4nM}S1eLnOFBTKtCM3-ltvfD(?O zARPTio~IC^R}hU}CVG~yhsSP4<>QbV2m3h2*H3PKs=?EXpYNai@+Y(Lf5rp^C_LJZ zle>f|f8#a`{@!0<2wQPCAhq+_c6EHnCG9Oh-*9f%0mD!z7L#a|r9hIe4YvLFQMS`K zsnhoELjZ>+`8*W)&}?6h%2K1F-sg<>Gaa;kJn zExLNYS3{yb&nsNq6Fy&imo$Yw$!)i{;}MCDf1yp=l(ieZZonFyNwd?yU(%ScR-3h^ z#{oN>hW_9ZPPvp(Hit-dH*VQw7?qVPg`yfLEP_n!V2UoEjU`Z@>mRnR5$? zf6X)!o(cIxd!4oatpR;aHO5%%L}9F++k9zwLPg#H?B#F{$%DyK5;cQXyBBo2LGwA; zD{;>d>jbUipxHU6@^1QG(Nka5|3ZU zU@fM6SuKL`DF^m>d6o|-iGJP}_GFM*f3-ieCv`kmAtP$x-9cY(?pAR(JcZcObU~tC zw{_SJf}p(v2Sh2B7dT^?Xh=cYU%Ln0?gnb;zfyX5%D7d#&m%@&0vyGm5=!IYD_KV_ zo~y?Tsvk0}LQSW&{4-&>ook>B7h2X zxi8c83688L(N7bmm9iCvl?YYxaoN{_8JD)*B2jR=88y(FE-QuMP8gVY0ZBQ-IY)ze z+^*`AAqYXd5jInn8~2xHg|-t3e|qd->t0o${T%zdxy12!^rM8B$|eFk{O4b9unGO5 znXA#-HAUg6#O}e)(JFM4X^OCMHZe0u(N&NgC3hS2_+47-s>5cZY z=7`*y3Kzi8V-*De0J1HU0k|KR>mvaPf0}J?n>rN7-}h4>D&1FIf_VtZZklOrmfcG2 z8kK63HeFQ?aa+`2XU42aKK->J1QMJ;N}}EuRBeUe`1sE`=LZKv_tVUuD1Q(QJc7?v z(9l!`VV8_O{B*8{?D0ce{nR_@I`e=bzjvZ2-J$Xg`siG}CN5-rl-m0Rz?st0f6f%$ zGCOAqF!j#A-R5x3I(Be}T*mFSh8F#tDZ93EK5nOo&qi?0H)3CA$PI;>KJlR^tR5dX z%{FZA=ToHi`e+=wXw&p3A^rfmj8L_wYg%qX&krXaesRK4GGDSQF};`&_eF{Ri;2$F zi!ivP2nv@=FMQ8mr?bz0FLwh~e{rt*JO`EqY*jeP%s3>J1w3+X9^ND_Ylj1GnTrpG zW}PWE|KnRcnwDuv2K>CGW5iDmSVlX4)OjYaVk;ZQ%;76M4TT!X$i5d$Cvbi{jn#}G zl?{JJ)MEom{KQk5(`FPhc?j>&Bci-IQ4+!a>nszr_O_sD#X&PoZqyb*e=~tF2fKr2 z8u_Ew-iSgv{~oi%>!u+wa!5#B6rGdW3zjw0jQ0l3J2d5`z|lK#(C<*3={65>@viI@ z21k>v#zRVYk+TaJpi465b91o_VK$`b zkB45rV3C+{roz{ZDVjb6==TG^DE40)Y?sMAGMS}`Yvwc86&0WJAa;HW{ZUl9vMFNJ zIZPSniKl&moM(XSR0R^4DUo&ok!DR~P2@WeIX@Np5n)zBS_MAXKKk zt8E!osWj8x)hy~=TQ!w6mG4015qsB&r-M+r$WUokrP50KSF5OhwQDMCD&K+1BlfQm zPt~bp^R*6Eul!Wme@dS!I|Wo~HI+4$2c`0UUPvW3Ic5hdcsdA`*|na^9ju|Rs5Gns zDnU(UP31dKdCYsI;AvkfZ$mcacB|D_FJK`Z2~|M`O-UP@vZosYIL+Yf>EZ>h0qor` za68MRSz*^6Mq5TmmVg}U3QrcWPf5r=N?`5>EV9teQxC5Uf4a?N9@&J_Fm9@my-0;4 zG!YjnAiq#0YIUPXEk7Tw*4jhu+LHUGR@p12l3GEhoK}6$$j-K_l`X%QjjUcyt+`iZ zcT|CF0i>)}t(;b4pUAE+DQ&8iJs6f%3(Bc&?G@P_S|#6Tl&;x!IjwtL4c36xwF?AHx`cw_dw%{r=^);geM`6G+0QmS(rjD#k7ft{DXJ|005*F zmtlkf9+&YH0wOPG*ue&L7z#$B9U-#lQdFDl*O#PTEK9UxJHv(n%Lj)NkG$vZx#Jx* zdGokQv~QO}6#_7Se*kcyImU$sZOgyVfCc?OC*%2^bg{t#hm6Y`mxlbi(E3K*JurL! zh{IfHoskWSxo)yrwZ&R8w@+KF2QwVyA#RI)CS-GkLPn?_Kw!0}7kEo36JtHd;!U2! z>BA_(E0ibfl5SZC;Pvd6czypnVwi3aeb577yN=xvE4Nx#=697D)vHSMZ9>8aP5uh} zd`*5=xKJ-QnZ+5)=^aW#oaq|Z%-OSTVCiCgL)n~LkO`l2Vk(dBO=oe*qGQ`!#^4OP zu8{Us@}GaD5$4_kj$w6Y7x-(JPZk0ie=Ww}->2vqcAYLE~VX(`IGFHysLmm9>HY>Do)p@ggj zMoJQ`;Nk-%`L^Z(XQ6!Y8=`u(OqfT>X3K(~UM1|cz%Q1*TC>;Os$ym*s?u@Ge|^_f zA)fa&D0B!x8gGRG2ULysbk$ZKcs*$R98vzCM|+P4YMHE!N2zu)-mzUCJK6%Trg|pF zUy2}?6i>!7LDUIz#Z^@>i*b=7z_&m*N9rW-JRbTZVDVJ&p&SCcFAZ$p=+1s1jNGEN z9z1QTggyl-sPh@P_Apl8sIfXme}BrFeOS+d>Q+zyJ0>GsB|sFovfYI~-nzaOWgU|)h*5LQ-l=qjRF z{O*c9l7u~F4;XrvKK9p)!z^K(jMrW zqAQYPV8DKxyzVUTPw<*o$259J647T;9CfA6X`FbujZw<>Tcy@ySz7g=d{&&QNAzE+ zGAtcCT!uRweT+j&M$0k*e+!nN>-dLAf;^J?dGaJT4!?r{ej;#j*xR_8mf~1A#s3O# zNCVDr2Hyiyyx+n@MIAkr_*Xl4TD$K5{{27Jzw-$1GU zze59ZqscWq>U0i=B0T8ZvqRBxyunaZ2g@;=oBg2({)SLg zTyyph3k8P-0MU$}EES20PhllP{04tqRG2NDlX3OU^#@Q(0|XQR000O8001EXMNG4L z{yG2vDVL8K0wI5GcO18|=J)&x9)H*mxgrcyp{f9F;(I8|Z+q+;D~^**_Us2)qHJEv zqF17=ID7WLZxwn-4u@A7LB0Ip-~aj3%e`L@U%xpze)Y!> zg7Y8j9lm;Y{QT(E_kaB0my;L!>4Sg&=y!kk+tI7%$8Ue%eDUhV@khVg+xx?xzWdL^ zXD54Kzdru?Kc4<{_{R?(JzqQIgT3kZUp&?!KiK>D<dcsM^_bJqYXj<#U++x^eevet*~!tbhpTV?@q?Gg-yc1D`tpO1!XG~T`qJTF z9zT0}qPKtK?C|HuzkK)d@Pm*3SI~R(V2{1M17Whi19Eo$f3FX8=~BKret!7!(QCb$ z!`GYFgI8bArLUZhzILL2J$iFwCwA$w{My;99sl_F*^gfzzd5>aanF8u`uc;9zC3(! za`l(~cJ%z@hqI%PUY#I(@Da_-@cH4<_dlG}9{qpm(aV#=*PBP@s}KI_=+)6rNB=r} zzH#Q?zIt=Ge(HbtaJskitMSP{e|~xV`tbFmrc>Vh6Z|kS_ z{OIKP^$PQU`R(TmCsBSk{p#`2PlvBokL~Xv<@GaKUGMot(b{hpKNTBa`ODMqzFS*@ zcYlB7uSbV(^;GOyHWu#m!oKHWSm1UNtUtv|hZp;){;>3)jI1@uWMN=m__ zJlNYO$DkkR3t51DOtYu+&9lQ-hkv%^d}V+4vwa4=a6q&hKq)x7IN>Hbufh^9#nSA} zx?@-Z-T{_$36^7t#bqYW4}U#;xfoWnbKa=l3n)=IQYO_eB(fA@sA3{UQG=OSSCv-=v^Q$jz~2oPVXw(BzT1=dX;HO4_2Mp067>Ea@l|TU6^6P z1@@8^@4~u%=yAE0^lEg2h^WoYStKbS_93b9#NO+mIWu1pF3YrJ-3*Gu08zntF4ENg z58kVe^^p1WzXr_VkeE|37t#(u=iThA^E?HRV*&y_XHXAMz$~ReJ3$3}I2sVAT zWpKJ`JXLnrTwE1WZi9T^1yb0*H)5#xZ2F+*(^p|Jb8@8S;JIw*^$LIAJ;MzOq%o=Z z?9*;gwp<>mZR{XY$RPytyFxK+R3*hEK5|z3KX29zJqPpjqfKhufLXC z26m;SNl>m4&Ai z*I$e21_&-E97K@omE>8Q)QALgp+YY$8@ztVNy&A+r%t)lqpf?(^SwYN-fgZ$Fjc=kXCAA{Y?R>BhdT}MrMjtDCPYrMzlantu zbSGOpk~j)^M$CVzmOL&Oa)DK{=ous+%tDlH+fV&EmjQAP88}ervGr?bc(8{-O*Xav zqXIA6RT&VYVzXuY2gao6)j!$=Bi0a{x{oSC6rL(5a*-oYNCmHYtj83a9dxTfs4ymN zpE0~Id};_5Kow&`@74A$ku7S7t5(Wcl~)WQ6K;mKER%mrbu(~h`s5B73m(}0o5K?9D)WZ#2&IlFNK!6Hiq5?i7|j>`i~5wi0shM3js3Znls zal+XpP%(c(9wmbB+D@#TLA9Uv?2MuyL^ZU0iL{A)0hPYN#KJK-_&|)i@da~45${FU zQ-l*PCm(_uLXH#!93*=NQOUwdZEFQXlpM^MPCXf~brUG&IDoo`kz%40fJ*2DAUK$g zDX#{FiDOdrin*|hl*J~lK4qJ9-?s63ws$U!)6IXqRiG(Y5}@j*;y#e>X~OE6ue#~=~cf~bU;blxpJ@`ByCveAXCyr_=uc~vz>@9J| z@<)?oPN^%i{@Q`JU@eVFmIqb*yTMLFB_bfFnnIH06QG)TWV>vR)}q8ox*~sDswUWO z_-7oHMfZ6}(JZt)WD26f5z!1{1`Av@hf96>5Y--50~t%EaaWi^-3%533NBii?o6Pr zAK2z{)SH2NX97~BWT6WUFdLOS=bh9@4i1nrC0({ipt@)IW=*uJFJ4y}L(alLh1s@j z&_pa|^Ao6b-D98E&7k)o>ydx7=52^Yhk6G-lVXw7e5aNT<~S&g?vvk5HuQbntSezQ zC=c%?zb}1u&H_1O-cCTMjnf?vN~3b+>e}pt;Z${?D}tE;&_kA%oCqrrqz)@5086=6 zxxhJF4s<}(sS z-rpX@kP{ZwLnTr`=eK|NR9m=eymLa~ifgOF$)JR}=N(n;Xi!>|ToMT?%!1V`n_7y* z>*&)8`7=9cKLmO=&g?0l9w(T!WJ1K=TDvZd7pr#q~eSsd5z(*%CHw@75yUOp*X*7exgO zlNw_&_CTv{0ux(4s{N#x-})U6WO`(Y}ndCV@Uj+$&6 z>e76d2kSSI1|`SUvuG!f?h&$%kU@EE=;{4;kC1hQ3`#;9fpLG2kgXtOQ2tmwf_7K@ z9wF<^*9Rq%L(fK52SPH1QdtKgA>n*ywfSAGATIbE`E=vRyy8r5ZEs0G>j}cHoWV$A zb))nm&bP2^$39%#=xXqecu*oa$Zl-kECA$#ZQgp~`QirGWLH{A=Si#t>u#Nle#xk& zTI|}gxPD+mSc!k)RwtV=69NQ97_mbYZlXW4m8xsls^CE>XlsD77-1n~xOKfJlFBwG zle4Icvq0*esE^q(jdz~C#wxswNk2d_Bo1^5G>L#nlQf@Q9){0H451FQrX+Z8vu)Tj`I77+ol~ti+OoCWqZV%tM zWxvega4Z;_F^0$`*)@SHiR3THy8g|x9&@XOA<&ndDMBv>LH)0t(r|B7%fY=$7 zapidZ0%uG%In#U-lbeJIP|Pt}4ki0AZ)D8}Xkd0JrbO1eB!ZO4Qk3jsAOZ#}aVwii z7mHw2IfH+NNGEa@+l08tI2z5CKV6``DxK>>kRzm#+#*D8!J} z__a2M#h3C_Zcw({ z$N6y78`Xv8EP1x_+cDC7al>o6D5BG^M`ZS1)J1>Ky3lUDFt_nlJyd4Z1g1PK=Yyj< z#EQj;R(B#3I9bm_Z$O+3%AP}%9oX^SC7*)@tck9W`sT14a-~L=diRi{HG$8TiLnnr zjwmE+L(SHlwXaBc)eKg7B8ub09;`^XU5@co!K#GEVC+_LkC z?8JY|shKViI`{Ri7_ZwkltD>wtHl3);3Tn{(W+Rc+H9&!yHIk)X#E%&llayL$}bNx z>fR`cr2zzq?a-?k%=5WS*Hi*=Qx0Bz=)it?9)HP*RzJ!}q6BILmn?fL36Sk$2woLY4qboH zs+z%|6@&VA>aL+WcRpH~9IG4XBZlAw>2lG+D#(pVzMGY`gD_0A#&8KC>#8MPn1%`^ zh~j*yPy;c7^)lIKa4A^ncb_B+#;$TpRpU>jC!RG*RPCvXGKZ!rl`<&X?typFA<8l+ z!2@C7NW@kapn)Y&V_U@y>NQwAoArNhf)re2>w~~af#JK#lhjR6@+Ym#lKpgcU;xOD z6^N_DHoMDKVQy5u9lF|Nw09+I=p2l}u}WcUSz(4#K-)ETL{`MYGg;IQKv&(8x*5)U zyb8ATcuk`+=@QzTLUmYW-9>iQGjUneiBsbUy@0huiOH#37Xpirs^GX8{Jnn&SalOz zU1FC)yb!LA*nqe}X><_S0T%&I5cS(=RpY#GDH~_`*#as~qKA-7Rq5iBLk_VNtwdqk ziEj>`s?hZG5L4BTTRjFTR=+CXIc#`HIsunQ4JP6uuE*;$0VeF zriqRy!KojiYKc6kX6oAVshq*Ug(k7e&yu1yQ(v#jB4^l~V4uvbYyf}VpzJ+`4omBH z#gjQ_t4eQbF*|b(`U=Ng9&QW5X;5Zf#lnt+Q+H&Gshw>`wv_R@}iC)wPi1c;P+qR6X;FJTAKt@TT>nu3g>J5MGL%6x~|(5$c~Kv{rW+MYS%ndY zU_l^}a$3$}qw1g1Qc2}TN=dAOg*s5y6o6HKd5$r3^%tlaUxg9|an^jEgGltM70agg zrOIYiARCjn4&;9sALj5duy1H|vIU9MKtzV@t-L!RoL2U-sHHgNpP%itc5&KuGuSo> zrHkdr9_%SFy9gR)g}~Yg%vKyt!jq50mKEIH-RUmd|t&cyh9?e!60T>t?t? ztv4uHehN2f$n@*S8Z=-nutSLu(keQyaFkhGvSjcch+R9dh0vwo_`#l#xQ6}PJSahK z0NRNw7$S$11yw3o=wm}blaOPI$?{gn2kVVvp_4v`=#hXtr`b=3lMGq+@Os6N@Beywv#`%d3wCQ3*&U#PQedKmltd_lxywA2DrRj zl_j=P1n8}1heXSGU!^!5n^v(bNpACzv#mi(99=sTLWO(6AA1#81> zLE8;VkQ*!>@_drWO`6}C!D>0PH;*aLH?*5aimoklMDfLfr_w+JC8^n~x(1NZ8Ztl+ zi0Wpr;+;@hbOg1@qhPYs&eGLa*vgP$mei2I68xY}hT zQPF<`2_)17PO+;u$GgVfg7g}c1xKcxa83;XQf(BBX;{SufQgG9&|Ed%7JxD+BaRaN zZq!_G%L^7zQmx~YmsP-#6s?2)!Uh@lq@q;LU}>2VFKogo$RuSK)((>wFaYy|TA*Cii?-|) zi!5idU~@`MX)7TN%7h1QXt1yaUVyWY-pavwg%m?t45x_TvRX5NGFjOu%oAmeq-KAh zRaVV`u!9z;8~?!`a01K8wCeqk^k~cDRoC}opi>Lt24%z%csERxrD!GFDGeb3%-7tX z&x}E*t1Y6)1-2*0EF4ED6xGaDa|na)KwRDUrK1PAXiaZ((EA}JU&cPmRtMI$Y!t1j z1#!dD;_*GnR9NG)yLN%shUqt?!l-|=c-}=sN)s%gq*Q|VK3J>`vjte(RoeRGm*+>v z2hX4Wd~*2u(aYm!KQ4xbdbJIuR{d1a-U_PS8Z?TX+hswYj&T0gf`z-vJ#PgI4fQGi zD`R0ydO5YEJ8lC94Wu6+sozT+QcOTrXS-x7W|x@MzU1hG{nb_7uWANsNhg136kxXI zVJ&DHltK>BKmn@fvMN5px45e{Sqs<~n7B5IS|Ib+Q3r}Kd131y|1hXe*HGxHThM?b z420$pOKZiWYY$&i)6++)Z9-eSzGFv*1- zik4uFa93I7ZCE-DZt{O+O9y`j<%DDQ@G!JA%h?4O_L)jpL{P*&0tra$*q5e5B6BS=?nQWze8uClGA(s%wC{h2+C+n(i;Z`B2@uFvFUC5-Tk7MNwmY1w{R8EFKmG?_^G%0u6m*ani8K*pYumSU0})dmxvyMac+J3I?~}^(Y_MqbgazHn@#Xj4uyD-I3j%#07ua z?xWjozdTyM-HroO<2yi!nk~ywqDscnWmE{LKr_u_ZOfeZYH33Nw$!n%a>I2KSYl3? z>c?6%yfmR{Wb(YBGTv^s_i4&v*pI$50&RVqEo|@O$W(uMb7A2fk>yb>5_P7CIj>RK zoVG`nSI%&QWI2vZmN)lr*b(4>vWqb>`)O4R90HVOhSMzeO`>`c;3{XhNx%(Fna?MH z-YC{^aa=7%a4~^bPa&n4awtB3vLe(f6F7UezQP?$R^czO5WW<)@ zv63AIvZSDModv?%Gs7xpxIwKL2c_N5;r704BuB@m1X3dFUYs}G{AFV5isAYxsW>KI zzGk(*CxN~130y|zyA3Ysy$QhE)HJ32Ia&draH@YnuuGTKO>mQ{Y+$;)vFQG;45wrj zI4!Yp-Ks1g04A`q8fr|_)SDR8nlTz1 zyk1@G*`}$Sz=I5oG~7aTnAtONfzC ztknz6hX&3Y^cWq)$VBL)nLI{*tHZ~cbMk+XNd$dY_;KY7CfK}w>N(~Vd@>e3S#xd< zR!&quR?*mGdJRgudu*tPyba8)^wXJjycBaPL!4wE&f2J1`7IGg`nIAN1Ue8`H@=k} zRh7H$ELxj)J#M{i8?@Bkp+TAVRimUlUvg;{cPYvWc6e(MS%GZ3_%u%)uvG3p$8LW| z0p#lkR=G8QwRHd*l#4$F5bOXz(W&+9^?~*|X0>AGx4UO6xm%wa6&j>}d`Mu%xaDGj zsV|Pn6Xu@gTXo}`p3GLr9Vjs6q)Jrf8z4I?HU`$cE|B;Z&wvKy^v__o+&9>Q+*<5x zLdHqFoQf1(bBU$P(OKuCj(#{*!&`qcSiIPP39Z>u+W4#mWrLFWrzpc6pbXu7MB=PC zQ!`1eBn)aVso==gwN~AUp4jEtA+40A3Sbpu>T;%x>Idj=ym^|KyaaEOY?D5zL(`yS z{VAGSxccB@FQ4H7lByP6nxHz97K2N_|16yO}WZnTQR7++CZBoB=s2h~dKSNzZ%;^?TM@Bfsgro*qmb4soRw#yY zuwFN!qQO*~C7pVZ&CauBE1Ektj9?t1z)M^s>faE5#Hok4!R+*!d+-- zXK)9koi3);_E3{t$H=PxCUJknLCtv*RLPTi8tzmLf0NQ~hT1WMf+`5Kj)(pvpc$NI`wU4yLwl0^1HKPaI2IRdgD8BS4?N_rl4<#4f28QYS8{ zt(36{5EUw*+TY%c+3ncRtQ-F(-ONF0{|RLc8IQX_*}k6clpH21LIi&pOaBd0>)j{_ z@`x;IzxPyZtDE5lk;Y@w{UN-k&Hz|mVotg~WE~hx+5m`COE1|~G5Or(gI0}ygU$et zyh|W0RVlks*S|@C3`WGqNQwa=qA-#8p2PwT8Dr`o*lI>sbU|=lP5vklI3^=UR7eI) ziLAUAL*Dd8>hLz^_CSAD{m=om+bl65l-#Nlc?{PXg<=d$by5#eWUl)VsXNx%7WEw zR-MAdl#ULFb*a%x^hwVSckK)ouY<{@2#PN7%poZ7vgILJ=Kp``r!6oWl(<)z?ZQfL zYJ)WSq5I;{(Z7n1B1v!*~` zv1_-l4wjSXJzAtoM2x1FQe>yY1%Twm2bD7&fU2B9(G^MWg9<41`T#Mb?K@g1Pg5y3 z9+qlv=uX@Zwf=u#)qU_mokCA3VJ;D&X6=;R*VEKj1xZ)^q^b$*E~pK*8LqmGd`BLX zkM}q)S#B0{0-&C@8p%vqjhE$^%V@bkkt_*A5gex)ZiZ&BVlzaFF>KmD%RpB*gK_qW zoi%(7YMkkT#wM`(jxmK~F^XPcy9xg~)D6ndgWj2)pss&`PB7`Q3S&-cgm>`?n`7z4 zTys%X5Roy7b*pQ>1B6(ZLqeY-c~aLMX?d>#tAz}i(PE5+lRo+u058d9kF18Ahp3{W zM;TVmpw3{-j;z3&BvR|b@;OLI)h2@oqslAd7G)@7lIpo7L2WluiFCD-52QMXBIQU{ zVc_rK$T=fKBNZoDj>MbDg;$n zzn~*buyO)3-^pe47E!2%VhM_l9^D0(CAPqBR3^Uc4CY-(GSr2-;7Klz zxB+7{$x4)FOGRV_0D`)7=+z8QWF0EHB5tb2x9EQ@Jx*PtBdmEe{6WSChs?5H`okx(nKd@xE=ngtYdjYxl5j2CJo%t>G)@yn0`S#3K^h26uUH z2@6eQvh-(QYv^jQ1<|#{toFOjLIwe4h|zBWtQ)j@Ai8egHeJUDd=pj&L0p=voq@T-V*%7u5LQ;Vr5CiR^| zTcNnwE3Sjjpmh97a&g*$hP21uha7)JO`f2}&wSeLyyOya0vsXgat=YcCm1&Z1+A8SxKQhIvu#5i%toc_6W+t@4#8|xvOa&oJ=EdzNOJ#tUY3m=~kp#9(5x?a)VU6Gc4Iz zbXHez;k$maE3M<^{po{!EIt0cBe}Q>%lV0b=s%=jpdY3(5v1)K@55uFv|-pzat43w;R; z2#iTzITXf-UA_=l(35|XSUH0VFO_9k;QW;>@m1w-35}RuO*vzEP=fl%Zql2&n^Ig2M&M8LdmkcM4a{HBwbh^ z0%Pc5Kr3gkEX636yCEb^0rq)NHhscVL#c=za8-gXMoEM*AgNhy0m=K|omf`0XI+26 z&4$5C&{Q{mFLwogRMxyesbN2`7fOD){@cf->napWx*unel7QNoA2vUzk>xo#Tzgq^ ztC|cKx)C&N!l-}SIWU}B+FtZaCQO}XyjnFrt24kz*1kzfy&3CU5fVm1K(V~ip3K(D z8CK!&YEn^TVe2PlW_(siCPwrb2`Fr~kpcMRsyLJofc6V>0Q_Dv}+W!8{^&|Id z!9nTx0-t+s^=@#h;0I;!S17yZR_`Xa3Vu*df5P28$8~=fI4@IoPqk%Od!Z zDg0~N8{RDk^L-2t1V1QSZ)LCVW4&%7)(gV0tUNDIm)L=%*|+kcF)^8G<8?W5=@b@p zHOZO=3F?0WqjWV6uA8A-`h-y#`OuJgLsRSc4k*EL=|jN@my$)L=@e322?|Kwc4m_0 zjD?qK!Ccu2~9USeCsZKQ$xWwpb$7D(O58hNpIOa%!w#$9D5Ohkmd(f5CO& z8I)s3NxKk44XMo-9mOoFaU>)wY&?~OK_iX8O3Qx+@S(>`tDND%o>lU`;HHHwxoKfg ziX9U#x3zbo?NsYTwjQb0*u>M+xiE-#IeN;d4m^5FJ=(5vhJI}~C~Xed(A<6tIm=Oe z!i3rC=E*g9z6M&>Y950yM32N1354iaFUe<`-g}Jl${8%NnWBqP{l8=>@ufDyb_Cd< zG&+AGG$e*?p*90p*8yFN2prRVJlH(wlAZ!v$#vueK2C0a))QJ&U(k`;TswnBXeF{V zXuTFS%_7FsA>0O~&9SNLW-rdHr*MqXin&k*l?OZv_W8!t7Co%XdDdR41bnJmGaHYT5J|n#PV;eMjm_qdIlNQ7c%CVbwxf8@K_bx$P^f`O8 z%7ZoYQmJB=QzTO3JQ-2mCB*6`C^cs$?j2JTPXS|c?^$qp53v>H)03y;E!G@E@=1R( zs@1$4SYDDkq5>SObglXhCRteBJyc^Us0z!JHj7d%Lss1k1~Aj)3Fu;GH^gUp1*j4O?TSIe1OZp1}rNwGj$GQsxsWh6s zm^xq#O3i}`=Ha>ohDqq-o_0dfmy| zUG#W=FeV3|1REL+SoO?4W??yHWg&4c2N7$tlzqsh+$XU^uWr7swOV62bXtE*H8ZHi zBr*lWSevmX(KIRvpPC%IPBaAj-68oF%BmY@4I78RF0GmGRw%(*FR==kxO(ql z*(-)G>{RL^4+f>xJ+}0r>sNn%NgsEnr1y1dP_w+YHhXkDvEJ5Cw`>gWLlDR7O)q=o zRpks;_CxKh6}|2~it!9NHUe!>irphEPqo>Cyb3;2NXAy#4t2x!4XQ^7r(sN%U1G`) zgITf5kw;UYcY?;g<>X8iA;pBrO1sEnf;8yNr1nn`3{lN-Z(Dcb)bfUujZ{*P4m!4W z>c*ck3zJ@h!Vtv=ix@NxHAcNf9)ST6TXw2(Og?@EI+$-4aSys3f^JY|UOm2+`plhp zYwK6xv#xj4#r29)Rd0U^0U&fi9YBjA^Qr}*C)lfQ0%OEBx>M+9s%YK%N1H7`>p(Ut zA0KhKCe0J=8q}7%)x1@VC})i%tc~wJ#4W}}RL-EjsM>57qX#YZM>euqTTnJCxgX&k zWp@W9D{}=(eZCysg>Z?U4WY^ztnwuJ3tbIOD5(ItXkfMAY)}GU_1tno zEQi%Yo5SrzH&Dqo+7DA$@p;*ro*iv|+h^l|*lCBI6NnI?`b+-y;9UFsn z<4^s>^>?-j-O$1r;+T|p#$B9wJXHJt$2oDQvR;y`Bi#~}tt9KTsg#{aj3GoLTatYa zCS29!*3cn}Qj~4%$~LL&6=GuOHnL_ygleaZN!@6#?(BKq@2^A6_R znEKwfOh*|^u`5`YDD2vL)(ur*=GN@3Tp?rODfR8Q{NXlRd$v?)Y4na^TX3iQ{%r@o zvf9sQ?Uhm6<$L&>7v&84)9TKZ2k-Jv$S}O0m54idGAM-GonnU`)~Gy{w@^EGRiemV zMn_s9?i}_AQFba#Jg^iS-nHEIvr|vO3a>lj_NgmT>@|?57;ryBC{=UfRZrY@>E-+- z@+lp(Nv<1twuLO{ZL!^m#$pqeZ9snp=yGMTP3>M^ZX`L|P}niTtJB#+X0M zRH~hOB71?Cd%Rs@bNJ$)XQ0Tl6RGlTYtU>Czvn>+OF&+Mm|p++W3X zeB|&?=L72}7lo6(Gp^R0n=1vT;fwR?*5&;M@`Qn!hF=rstpca|kvkQ}_GA?zZ;zk9 zYA0S?mwDYux#9J5kz(c7m%)$Q>Sv?|ENwc?U%Xls)XOjG%l&Xv?}+M9ygSzS+epbt zLc+na?p*N*k6iq6o@m{8@-f86m}-ea_|`t`)AE4FH5<3A@HY==wN1ls|M)zaR^hplA`g&A)8c!~v6KUdrmFEv#7bExLqdMG+*@^=dHR5iEVDA=*({3*Rl0toE6yFnIdDv z=ch`F*#U&UF?*&j+{W#kDC`t`qlSnW@GTc_)l8ts@`{8Z!MUIa}&uvI*_ zY>!^XE_0~lJg^5`5NLnfDY+V*0&F=^qlv^Cm1vKglK80MGpO5hW@1#`ShVYfGs*Jn zxJ&%YxlO9)pJLy|=JvLt z?*FaFaF>ztdRVnyZR&kDkRkPkSnWa<;y#y6?P$%aTn$2iMhS;5Nn!5vCr3iL;NF8` zyD2rv4N3N9gQBzyKxIGW%0OtF^>Nolpaj4t!Vz9ede!O!JS zBo8{fxYC3M!T5ZgZEMe_-!W7@XY%J%w(Q%LwSDafHTij3K#QIK$U)4D0$-mAQLM%9 zZrJ4EPRGTJd+f9Z6|rkNMo&?_*oO$fj_lc9|j<*ZkGbn4X$J?mzY9A-($W zuBNxna%Q%>A_F6w9vmWR$%P3&z~i0@hv6U0CmLnz<4VilsIn-o*mvyM-}55p3PGu# z6Jjbr#Za%~9v7J`|A}V2Ne?sro~@Hn95hI#)&~GZ-0h?Loppu19CP-E$+7Si3k8Nu z4|;|I79`!eYv`M+(*G^xU1m&wddtVJZc#`RS58$@I1#?@RgTa|$Lsn#?k&WVGhEr* z9fRMDWfWld)mo%J)VShh&)S;Lv_04nxMElsyW6wkXTjS`pF;x#cV&(x4A^DAR3kd3 z$iB@kE|E!N$ez=Db^6?lquJ4d5~Y*r(#C^#&m03K3vpLz1=$DVo3E7raf0#uONac4 zs@wpJ$e4t>#H!Lmxur&1wR>7rj4)r;0HeznCYk7c7yl7;R9U!E5>hnw##trR<$8o1 zAsHq7{RWbu99}EbpkvOKp}&67KWRb_5Dxcbun&xM+;IBWN7=Bu>nK>}O5T^|LNoxS zb9;VcXOw1D?YraW4_|YckR9T$xJ=H3vf83}VcbmKT-o07nlAD%#M)} z(>z60UJ4x>(hIps`=xkY-N_LC1yZ5?t(L^MIt0Zs_PUK`Q>>`G)Z6(Jx3;J#q=jjz zWNFm}NmRH9Q*srf1fAapteCV10A8MqY3X+Pqmm0JZb_>sWMm>p3ZWE9`<$EV`1@9A zRxaeUM|m~}J)>y#colgXTRCdCzCgh`g3RfA<^SFKa}eke&0E+;g6#>^dwo zkde(1(hPWs43{XKi#s6EIVBtP;vc#@Ys64^AlCwC?zYsGw=v&bZ5i}e$A}9$E2%!# z@j{BK*G;K}Q!5in=RK09ipS@it1eg0y9 z@BEQo7@fJJMaQhy30H;~`CG0yaZAKX|3&YL0KhLQ?$SCYDBS1X2c)U-W=WF|>S4}G z6E23syy=*ln!{aPT(m|ZNfk!imLl0Uf0`n7n$#&)SM^AkJ}cq_SAn`{&V9} z_vko%tUwEcV=Ylq6_^x)m?OgT1vwoOf~4fQEzG4O1+HSQR@wYNxGC<-+6QCsKGr`SlVqa?|`S=lDT z$k3&}UeebOpKo4<%e*Sfzl1V^5*J8kPMdJX*Y2OkCxNEX!K>?Za^?HU9OXvi>~v?y z=GL#qMnBRqks4?7Fr~rL#8$HQ`qaj++6`vcVh0EYFEPgH%larq5;;UuqAhD$R-%Ap z<+gDHJjSIp%V4x^-J(`=+!{G^i@h`gwNvR0y+qa;)8rg8x4{jV@2nMB)RMq_Gz~Y} z0*C44MqSPGaKt`3hwArU9VKVwNLWr@%v)c|idx{z+)9{|a-}AQb1Dz99&hQce=(MS z2Zfq&>Q1HNLY|3~>HLgqGL80INUdx{rh(uv&6B;>XYJO&9mGX0Ei15TfruPqgVp)D zyba^k*-y-F&kGsP1o7YITFg;NIJy{ zVGW1pAac0j^i=#=qPQu0Wn`^6E+-|L^R)3ZY3#Sn3j0c zvG8#T(~Z#)r1ecVlX?=VI+)vBsf(-X2Fsai4F6%)Z0JkSWW5sIotK9v&K1wDBV$1N zt0xHv{Cz+OpzJgGks0vh2ANtjvFX)0;k^KZ@OAB`#br@A!5OTbz52#X8?bfHDgUhMP8YFHo*A$d5*%mwSY5T&Ht81oPE+* z6%rl5vP$jqd(8-MoZNUn-&;wZjO^w%`%LF%7d(F0D2j{(rMthnfxliW($@l!$~!D# zT(y@j=`5}EIw5j!QS?bQM4w4ba6!}x-oq9&)nqY}KQn_3SG%NLhSgS>3h(EsbIfj2 zgHwZT6UL@%!#6QF^m{7B1h-UM8O50qF|H@U@@PT`$-_J#A`Z!CtH{&{`oy?T$d16ijZFTARLnhGG^GeRMUi%}Mm$4|+>j zGSjp}WGHoiLk`29!V+DLkJW9JEsE<-;=WJnOWw(v&7*Qbn+t+~fWW^cx3x76Kr@2B z>lzg83b+9~ln7c%NT|#ekc33gfIOh63>y@nlW0I6*oNe30tzVsy+8vRf|5`V8cS&7U4b>C za8qC8X{|C??}3W0z&&MDz*-4fx&mke4%M46Cq3bqvD!_2pQmYWHuVTkk2Gv*EiXU^ z?vIBTU?u1RmBDYoPcl^Q1t>vE-f&agK|=6^BIM@{r_2v;>N=jbo`CfpNb)LtYRf6O zDFw*sDxeFTr13NXq-}*tc|P6wO%9yj3UPU?wzwI4+6Ru^wX~@To0>(iT~~s}ec*Fy zt;2c`WZ(9g=li~hwlb5dYEpK`53)grc4~Yf9lYoPc2LP6UR6KMd0A4;+yrA&r z=pm=u|1*e-w|N1I&{_aIr$-X`;GQV>8)O;CoALX6OVEATlEt&UB|zoyzgo~(Ae;ca z;)CWZ*sUJ~yQQo67BC3V7VLn^g5XZDT6s1FNIn?uWl9epe0i|7j}J9GC?Ds8&8>=Lw(Z39x({1!(9c0&J6k8#I5Qu4a-A6E%BsiP0L3Zv?JeV}S;6fBf{) zuB~%1K637l!{OdLzuxOQ*@E+xH48*Df)X&&s>-Yv=(${`GnDusj1sjqus{d!D7og7 zxonuU4j!Csz=aBY5}k}^vRR-_2o1y7%MEDZvq+d&e78=E5RYz+8zDc$8gNd?CDwls zgc3pE)-Y{GN7P#?Ia)TKuI>z!#886<6FjCF#XqdB!l8jbZw$o~I)KyuPq?d$?|D30$4dnp0bkI9qu58Qg%7JZ==TcL z&4GuFjkD?!o>fkB>R$dgXdYA2loVx7*|cEA^d#tDFr;v>t=iq=MOI9^*GHu%`sL)0 z(Kr_!GDwqSLw~IypL_#QO9KR#zFYztvoamAUj@lqt@P%T=zJ!B%Z}SN6y48P7~0d~ zn^IsaKr%_YD2gI8SqFSvj)k!v5b1WNUth|$WP9Rt+M+==ws?7XdG5InuHHU1+r{rN zgt}|j8Cl|N0j=#^-9D_d&-Wiuk-e>MR$&|aFAsz6H|JoIEqZMq^#gn!>afl>G4^lr zJlG91dfDi>`J)?u9>dahP2PB%@4Yll{)u8r%$M;OOCqpjIIC_Ji~Tjrrb!0|SibMt z2v2blpLjGI9pM_M52x*aiC5vVQ$myoJCS%Ll1P_a zAx}6rw6sKee&m(>*WvEDwVR=9>p#I&aoE8se_@=ky>4Oq3YFhx%*0143)mMF8D$;f zS}Ig3FA(-+BBTgY&S&{C-25!wqiN{6c(23bwL3@qbdQyWP!A4i&?Vwbl__qe5ygUO z2TUjpr?~BZ|C{H0b_`GbwzhSwyxxYS<-DE%t!|*Y1+RD8l=uB=7Wh~M6{3VIUNFgv zQ)$ewurr6QxI+G?abu>A1$9Pw1QamhKsC~sQKSThf~IQFM$D2V$1%wnUeg&1Zg43) z;#dmAiH9TwbqEv5oVP-j9L}63w?+pTJshiRUY=urzE?i}QBDq>;3|?z-(4ERechb( z#bI@h*e7tLc$v$~d~*DC8V|*t>inwWpUyf+SWRUT7r9i#Np^yf!OS8kxqyOLW4OFf zJUP~F+odKN2g!llnc)&?yKOgARESl6_`KZS=i|1k@%LpDsY-el48Xq`GGnwF(&^(X!5$yt9(B^<{6bwLuwlA(-_}Bx%ssQPcMFc zdh)Bk0JH0V#smZ@zC%QlwS*~u+cpsXo?jslThI=0`y2x#$a zbD>Csq?#CO|ND-T9LtWJr2Sw(KHc5(9PgeJgHKP1@V9F%Tm;RXU*-UA;fcA>p3qY?O@!2Jnc>Fp7Ycb`^Y7vZ2Ik3~^c{ZRV{B>X0lYVM{)&9($RPkJe45@{8 z2Ynget>VT!fymNiLBejQdDw2X>xVmVKonwmfise^h6JSjy>rlMZ=eSLE2Rf#j9az) zJR;;Jz)>71p)?-8l2zoQxq7^$`XR+C)O1?OKNXhSxdzIBVb?ZN1%4HgAac6Q*2Hdx z&o1tONGLZ2iJsSglYepeJ)ECK9)g{s7iu=ceAl?JG{-nE_Zkgx4-s!5P^zW;ZZ)id zT>L@Sjv|X&XJmy};+8m{`yySB;K*tc{xV@&DO+M#ickd~mwg?Wa%o#l64q}wqxxFY zWvMXS34IeSATDM&<7hCC+f_X>03nDs#AeE3$J9KSEj^q`gC8c3xn-Sj$&KZ) zXA1!DLqeYrKf(8cEOOEl4AK|0C?*5XKF&zJ8L?7!mx!Z zHXPWwVsBtyDGBBEwrbZ`DZ+u>2CzP9!Snwp1=k|~N-Ufb&1I+Y;GR)UA3?2&vu+`DZ_ z{OleI>lYW<;-Pfu4fnI)h}@bA@c#;Z0|fv8ge{kFW&t37H8C_Ucx`O$T5WIQI1v8c zUlFC!ed?AtF9hx=t8H2KklHP3s#|GSr^*DkjY@3Rakjv(KRY2LB!&c7(0xeNR*^HF zc|0?ZGk)oOocR;w7ovej@RbTKHB~{_C1Vf2T&W>@{?JlCcF#J_JYdN0o+(OasC9e63mNZ!rFMS+aG^A{3q`lg_JsmWz5Q>eF)n-h)_dr4Cw zdOacTn-cvOhOX4>Fu0)z3YAMOyycJc*{8pj%|KOuoI#KK+OmMH3MCnfLqb`=J=EfE zOvb)}Zl#raTom#m*e` zGZZ7;q#;h;mAt~>XtGgwNC{7Jb`1k`L&m(HbDN3|aKIq9)PoGhO_}}=kKgW)?;y7y zVt0T&_K2vjc9#)Uo_PorpR1g#^lj=KlMn5G+*iQ5ArVFq%RPDzFhdl>$zGj@kn%>Z zpt5ZamLbfB6#eqh`)4cyGtOlAnlMGvrvUwY;*%fS_I}PMs`ptMv_>&%9g*z?1)}k5 z=#Qebm05^szOFpcv@03s88WsHB_l9XHrgd@G;222Y&flg58fmR+mes0y&BvOLC*@<5?#E+f zUj6=@OeF`F5#qBbF${-X*ns<{fL~*so$TI+5Me^Yga>G zani6#I0xlA*Q5&O_`?~f^M3@nWu}R zS_9bIk80b?w^?Af9!48_NH%^P>IzL3u+K@z{*+L;7qCb|H_tu166iD%J2HcRQa9$& z$c{|m2u;K^1?1CIf>t-m(DLzcwAK!2*EZSnSY-!F6|sVLHLUuMo}Fz*E8B3d>RG)S zT5~7QZpi}K#z$4GS~aZ3E}mWAO3I^^EeWfl1=Y|tck=8OtdiRpm2@mvCwV8n+-R0k!K2 z967Tre*^#kycL(>wgDcO{S*QsFKgIA2Xq(;Mxq@dvglD%o9x$*q+Tpbv}8NOh5^e5 zhZ2vx=kB@V9W{CLxJk6{mrfM|Fn_-RaG^QIg$8ZQztDgM{og0!`JQyK!2*Yj%Nv)5 z{JYTlMx8xAZ!yrH+rZMr`i8POw;&T- zaKb82?@ecM$*yDDT*lzcx~`D+Rq|gy(+G2K0mrbqvkUx3jx)*9JCZW|$gb0SOk*Zs zHpyzNk%B3kQ~W)~Ki;S4E`PyMgMKdh5oHg>kUwz4w+VXE7I}J|#_WpF>yV@i98vUR zs&f_kxRw-{PZeeuLXXw^wpZZ=2jImJ5iO6TkgB23h}(JL4Tn`2-0{f3^<@_ zyr-+S^1$msni9|4P}f)C{o*nMeW`$l*617YMAt@Yq(Qzi5%P(hu~z_o|5 z`bLe_G5S;1?8ABvRDZXEYVd<@5GwR)sL<|(YI~puDm@3P*MRB{LpAFNfJT1=c(4a* zbPQA>tv6QINKLmlRz(#ZiBsD%CH?(49R>Rm#DcK0!b4XP&Ej|01W;b!*mJpQ8w?s2 zn?n|RmB9cG7;OI&2A{<+Y^1oNWY-1*Ao?vJf)PMrJaBq#c7Ke%U_2_t+fbdx$!QEi z9o0w&$npGVpPQYY=!UGdWxP&NlI@b+QzU!Q5?THIkT30lzA3sQIR*yox5?|y^8N&` zd38*qcO(&gCdE-#>YT=jm)jVnY`;}%O_rrq56Wl7sd_~Jr7FYHvBPD!!_mh$q-3-# z6R=wsh=lLa^vv33E(FJ7l*x#t7$2Yg;V^m@P;(t3}^5?Fva^VJXF-t zlUn$Z3Cb4U9_Zko>4<-|gQvCY{_o%aWBogi@J_DW#3?=ysL=Ekfx4ujDBq6%lH}}X z#kn@i_|Yhf*2_4p1t#N4Q?E+PJ)~0JKMv`flQ=C(i+}jO(MZ0xl;_WV?YORq?(R0- zwm}-nO=wj1adW>THa8D$E<3KKgZ6hn+Fgn^s&JRM81z#lD#=+CDR!)Ny2t5+r|ML+ zSA|Y-uRKntb*fJ1a45oqzCAk>Eyo)SMRl+ov$@$Hir{YuMa4B||FBSSSO5^s_*s#t z_!L$$#3gU=$3=zN(m5Gd-(vp)P)h>@6aWAK2mk;8App0PN96lD001zTloiFcvA3p%+KR7sg@#N&`@r&>O z_`xry&koZE|NhbM{_wZs7f(-Kf4}$&ogv(YzB z^smRSj_t&*T$bNBn~mcio;>;S>yuZzH}1h zZ_}S19{+UoV*R!KJ*K>QM(gXnxTv`C*Tvt8t&jZW@ps>CEWw*UfAZJkqt|*Wb}d^A zcXnZaIr;g^qhF7n&yDcq$>XQrJpT3Qi-%hS__nTubu-*bRmd}jmVUQ`QE`HsbC?W zbBuZVOP)QqZ=M{ze>nQHt?UcCwXNgZJ&5Eua|Xupmw4|WhTt|2<8b*rw()x%+IT72 zW|#3+(DuzQ-@W?n)ncqcjL3Vy8@RX>V26ptWnH`+eNIub*VR7Ml@(41n25m(U~Gq2 z?F@S1Y>2vdAPP#_e8_TO$sGV|_9M3jEELXrz?FV%DXsv5e+rBgqfg=igiNx48bQhH zy8%)L@iDZ+s&ayR2MN)xUZtDk%|~7isx&01z{fPJm~K@&E=JWM5xATIQ57mANhG%F z!Y)M;LGLq4!VZO4JA+<064?!iNyYAREafq=6gbZ6ptnr(g&3SFq!>I$)kPtLodG#j zNsuszPae2OfAQ5#uv@(ceoS1ckixy-iq&jqLdKLB1!qrU7G@~f1;_y?+p}t>dk2`+ z{6U{tF`*_=^3(ydx(W0yL@zEzHScM9R>>yd6{6%-q$NFA9e@VN!H|&4-tWU_1PAOT zE8c~5{m9dDCF#`)LJ?6Th_gsiKpbLHPmjIVL33uleNshESb2he#pJL|mTh|zJxhsx53Th6eY9;- zF_mdgf7F&N#}owbSaf+KWi`Xq2Fm1E%_22;Lexd>F}vP8xCtbm9R$^kvVOjIz#QE9 zOM7lHxfZkZ?*z6W!aPNOtKeGBM1UuD;>A#i1&gxUw zgHWQ#Dvo>PSLF=X$S=AozqrL%ACrz(Z>Qb=)Bv|JIr(x!_p)s$iGxJXL|N66f2ZX_ z4p=3Nosp*He0q7aZHL{ z{iA&_Vhs_fhomAz;i-Zm7de83RPd_DdQ7p|L3bKyieu9Dnc!XFQvqx%ID|yrtL$jjI9KLS+dj~nOX~4kI zVFQlFWZyx(oc*|gXwk}oi7itvhvh-bgxUENAY}ErK=fA=C!AeE6(iA;e+2Yh+liGk zsP^-oolz8#NDVDtB5eX+LZxpov1m*V9wO6ze8C(@#Cy^86yZe6$p=tF$VCMKgJjPr zDp@$GZLMGklEIAWlzL3AstFWx45;p5h#|!U7)$5`5;~lYDUSw)Nn=v=iaFRvUSbnh zpR!GQXxn%_+dG%W>E_;Pe@yOOO@-7G^pM(t;WNkN)vL;7hj%ef^h>Xf80#0FRDmY2 zEO^yV$z>f)RDO7}otU1yow^weTPV(}grbBJfNj$-7JFZa0|m9{MH+UYX-x9nKsT#< zZV`0kQecR5{@g%@x`-@^QE?Gc5DQa}s*3bQCQ&F0bwRFbd{v3Te|E1duR>gGpxJ3m znmva-g*PmZIfbZ;imIb1=OoFq8Pbwm%8-dw;;6N=y|T;v!kY1KFW!^}W!IrgUEU=> zDuOOWwfsu7x+*27PMc4%W(9hP(Kij<&xq2P^txr=?}uL}(NR)yn6j7opy~M-Q+2M0 zeD*P7AkQTb0o3sie<8+%UQIc|4op>zulGQ*zGunakb=<@dYtM>y}lIGO7rT~VhdAB zV^Z+o1MVkb2Z*9q&phk;>C!Z$uLkc_RRvb$>Vlt@zCssJ?QQ+fbMUbTa+Nb!aE+)Z zV4KYz^~T^u;B1R0IVdNuk9TG3FvaXp7r8`ts+;G#h=yxMf7^iq9F%cai*i5ms5{BR zQn<_@J~NCXUX#MbJ&b$h47ZnL5e8-4jlp?GL=&(Drv$>B(Sy1hrjwFS&gykoEj|(I zsh(9ee!ut%gYxJg{H~P$U{ypEN_?~EdB;)##!4P%j46d;Jf|M6s&0ZE6#odLa_{>p zV6+omm{rG|e{b)W8hl79S=2K-Q$s!As+{2#b>Xnwdu*3r%^8Ocs!&9x{W~u9&l$Po zse=`rAM--xU2~09co~yY2cPHT1l}^+iG!NT>*^dKdrO?L{Lv(tQ|ijBzj5FlSW9D) z<)DgxKiFxgM1;|)rjTU$1f*tOv|TnAH=@L&bVasQe@(F4@Xs(Pi|+HDqFLhdkSPob zME3zzIa__f}Djh7G~SBK~rQgn?Isj*FE-m-3)pk zvK|>Xf4mK`=uq#_=crhWYQD#o4dyT?jqa1*Pd4=Xx>?u4Y)~HFOWsG}(XN2dkarUh zYU6YVgvL?1a&>L?!f>iO&=sM~h~zO#OHPCp2vUcY6QZSDt6bomEe9G%b*cm%xGT4` zY67c3BVbdSZzb(E086 zo@xu%jdwvPTyt$zI2n{M_q?a79Suq=CYMA&g&A1AvZ+f=xQRZkkUz77_H$5iP`(=} zwsQr|vPIR|rATv_5Br;sxU=*FHg&V!ubFBt>KSonPC71fl&-3dbu*YkQ&XCaLLhDe zf7hVoIL7&cDz~b*Fow;4%c*h|64??qZSU41-%OGqWfw#R3`aG_Lg;~3-2^7Kd{FyI zO+WnwW3)19dVpIW!RigimI97(P)gk@fxR21lA6ctlIp0*wuxPu@9JRvM$(|{1=)pE`edPJ#1~+6^T1n?gtN?(cWb{iWHPu4bmc`8j8^THyf44f> zOgR!FD8eE;QsE~0nypk_!&U_kN^Pb8IXPL9r^F3tm~d!jyO$28t~ z_8P13GA8{jA-yAX4F_@=JOX-&qF?{Aw~$qOQ1&W7R9{K#s_j)dL%%^aDA^n&w^TP- ze=X%M^J}fxr23F)3YkXYz69HAe}*%Z%wAa)O2#CJCFJ(-jXUm$(tj*f~qjumr@;pp0wB>lZjdY1&05?M-;U4R&o(MsIPcGATn7*)<-e<9K-Itf+{ zk)joNqur!`@TlxGHdPwwha8BYB&&VwNqBY-W!}qIg(MLIWHo-RjbZVnyc6`m$fS00 zUrW2X39MwZ)n}$aUTrHzN&}HteU`-r{#QbP={Et7|7~Q5u=L{lC_~`YtGtNB)n<{ zt2`0Jabk~FB-}2?cq*_e;Q@@jvM*Ib?7s|3pGTkClS1y;c_en?f91?f7lh7zy=%tn zb`51v65J~BzZ*D7tY)+-mZ>(ID$_2MTr*lfM#dz*^$X>f2N`v5l!Vd%B8u(Ms~ODm zxlA`y0&rldx5-vImktzax~QY3Zk%TfDW~mjCd))k&G=TVS=1*?s-|pLQMH%Vq4a9r ztAb8RKdjo!w-`)gf3o8l!aK6jlCn4Z7;Gy)M5ilm*B5K@od(=e+d(~z}5$WlLEtcl_#m2pyW?lnI(I5ejp%b#|p&tVcXqht1vez-;P~v zGTOV6HFORpaI8|;T2`3h1Z2C$E=DV2;W=8=4q{i`lDZi#db|p?^mt99GU)>CO`$rh zvhE_g>Y2DK>cpvWL|)KZqJ-qstpi|@s46&Y2Y)XDe^%WD*O%Ch5H5wQ6Ez@iP#PTt z_rOJ<6C(XIS=Bi2Tgt{+ezu^Bli(4OsVZIk=!hYdqLm~Zw>hU}Za~!pb`z{Yr0w9S zj%<--qZ024_QF<>ZL;NrtmeKM{mb{EtsNbcL{HI5Cs~-UT2I;8_OJ!RSXbg6?xJo2 z>)^wpe+NQl6>>3o+dH(vU>ps4J|J?W2J1hKN#m<>*iSXW#w|g{WNo8udq%}eSJd&f zfs0IHmDiGjH&b7)$|7gjoM4~Kt!x0@f1vDrgbqvVZpD*1XX{FDYB4)=4*Ckm zT^()*!f8-uUd6(mgj089i>aMmMz)mkc7T*Y33X5{(r#hXO>1={+1apZgGn5QoIGY` zYr|woE$#xtR=X93Gb{u3QQ47`mGn1#8iT}!AUmvC?Qcy zJyGOU!OBo_)yBpxKhC7*#X^9+Amf zEFGrol5PK6*9Z~m!dvYOwsEI8B($Kye|3)@g(3OiV%FOm(}cMO$PG$;TNeL5Bw9ng zgPa`av#wrb9lm=FseQRbgaM5*uzP%5Nf9E<0eRQ&^e!60T>t?t`tv4uHeg-#b$n@*S8a7}J z*pUQ;w2saz9Ay@lEE&89VmA(KA#@2G-#d^fu3`T+4@!_5fcD}FLSRT)NTq^>J~jk2 zNpu`Tvb>e(!FuCZ=%mjOJh3rL)CH)j@r|F#YX4bWKxXouiWDP^PQ6^Kf5ySy0)!g? zHz?!2mRRqI?c^_Po*u97(m36=6Zk>t@`8l7zS z7jK82?cd57)LlchuQIb0$Irnsu`O`dBD2cNh%u;>*MhV`S@LtcX)n8fEjBs?Cv5d{ z_<}bA$vd>s>t?t;8y)tQQZtyyc0`{4p5ss2qjDH zEM0wttqd7uNeu}s;0Gnj3+y1sd|1HSC6gj}XS;6eBmirUN*<$NenD+lKlQjBRaoC2ZCYRy2*(aJ_qo+xXIY6e9^R2mg$+KtYZrK9 zn0`YFMy18`e?B5onqYw?r4q&W!D3^W9l+vkrL9kYd3t=Zdiwb1)1#LUo}WDVaWOR1 zt8FN?>ZgMCR#5HEpi%7HE(`i>go~#ZEWEAU^G>kPP@nR@G8V?9ms3l+<1TQ}K>8s? z^?M^j3JIgt*)Exi*(Hu@UvltZKf0>>Rn1^6>4b~|f6UfAtOZSjQpkV}6rg%8>*6DP zhr3#nwLpD=i5sJ+1u}mfb)XoN7q$-a_d$KShC)}}f(Bjm2qBar3@8lCLp3mpO93o2 zf&+v1B@s-QoHzFcy2MW14EGL_AuFM|$Cwel#h6)Pl1n=jEx{V_wzA5*uyh*SAD1mzhkNc4J0g;A=N*Gja z zZ6#@UhLwiaqkF;1A%Ge;Q{M1d2yj<@murUGe*s?JRw8#-cxkCzz5}bar9@D)ss$;Y z;ou6Tl5Lw@NlBvV6zDWJ|K zi~?5HH+V}yR}JJHT)frC7gx@3??9~o(3D0hn3Dea-T_%Mb0$PjQ$A`#6-FIsV5pg7 ze{0Av-H#%qOpZa;iVNYUsVf2iua33@H8>{0yaASdSWQiPSpkDf+vQd$de>EQV|LK~ z9V>iW8RA`FXTG=HE#ijUe6!^na8pc{9D}NonAN+|C61~|$b01EGe-#9SuG1eR5gQ2 zA@*h}WuH<|L8VwD>m|yG2}}LqAT7ode+^1MTT;MYR2e`>Rxj2E^pNNGXK#*hC&r+n zMEx_LtgE_(TLGx?qGxDb$fT!_W91B%m#lVFPKH30_etMZ97b>n)P&OO^J)+^1&Pv_ ztag_9ro*=iI!u%T+tT@@rtsoUjwDenYYS!!7IqF;N}03E$x_#Kl`}APq!HGQe{cOB zqRZK$WF*D}1~>3}ln3^tN*1XNZsQZ<%Y#t&WVa`A;LqB9bldG$N9(uSX<%x64=7Qy zWjRVz$ymCK3IP>prg^Mwne$#PZ3JXX9qTGLTsMIw=D<`x)}rC12~8uD=M9zd_Orbo zr!0oU;5#GGHpkh)_CAeFl{Xg_f8G;W9@HXHXDTx14Jw<{?#S}W8E%m*r;*9>=Kc+P z0vyKdLP*ShTGavvz_QGEmc_nJR4)QtZc~*FOqVwn-QSntl&k`$B{puFHR_k`&Ty)m;Wp7W zHl4n=Q1q6bkS=USKZ#)>Nb$_zynbg_Y)st@Q=)E~a&ud@hr{3`{9eMkpkvIb|;-=XMEz^l#lRLPWPdgmlj*S<~N4E;1ze;AZ4pSp+S9v1Uz zzw%l^6xats&w^fZSRRtcYK9C8AmjD=6uoSVVQiR#-b8kyYCU`XLi-%DS~2t6-LsY4txt^#ji~=TBs62(ao|!%vQ)9 zBXG<~m8i-$Bz9J846S<|MDZ=2K^&CRKZo6N-(UxFYq7J5nI`daCQ@|GC6q2l=beu_ z`r%X!Z^>Zce{usRvSv$Z>%A6~4NB&pp^Wx`GIH|~iL>HN%_Oyw2-IFu!I75757O^E4@W0dJCQn?9*S)1YMi8Jb$Sdhnr_&+v#zRf{eS zsLqUw!KI&anPNnB@z^Xz3uG5Z!YQB-NKiS@VS&{Qe_!}l?|c@$5K7{Ktub>l?|>Dm zB{L&#Q@?en877+Lk- zB#s!=f1Hm(DtS^*!=0+(Z&TWhykj5@EZPpdf`-@&@AM0z|FUE$+1OSK#79SJP-UMG z5~wfO!PM4GVA}!Z8^_XC6+Mo;5uo3`cf!e}#4f28QYS8{t(2(<5EUw@+TWgx+3ncR ztQ-F}-ONF0{|RLc8ISuw*`c29lpH21LIjygfB%iB*1J&<(G#(#{oYftt!{=}L>iAx z_XE76&H!3oVotg~(K;}gv_TZ7mR_=};^=di4_Y<;Ejj}{@-Bh2RHf`jUH>KlF&Ysc zLR1V$j1eXh-$$`PL#B{A2)3Hh6?!O3})DG3IS=qz-RmZVzPD ze-9l{yUP+25|dkZA`fsA#4BR@w}vu?cK(v;@G&Z}&!I;H-X%bcQkOI>D>?R_gU761 znqDYUw(n@6JWZwCcvz~vp*wLu)cS)}e-F`v zI)y%_M7czSnzd7MU(Zrs6(n8tld2}LyP!7McDU*;@*R0lKHk&3WVv0;38H%1Y9w>a zYP>ARTqer}3el1<6v1Jt;bv&|DmEd+5MbN>Sq8eg8H}?}?5yD{sBxwT8k@lCJB9?w zVidi?b{qb6s2h}@N4+z9L0to#e_+yM9mbs02=C$(w#U+ox#6O!AR=QD>sHr%4+ya^ zheSSw=%c#s6qolZuv*BNnJmUwIO$j40^lXN?1|No^B7cA^eDs18PpjJ*+nbxCMl|Q zVfh?UqN+_s5hj&akz16Zj7h5JmISr^NTsN&ojgR_UQ#NtsVuG>P_MUZHwmDji5AS&SZ1gLwsnF0%?j6&5+6_vXf-^P3Os2otQF zz|8mPGI@(AR70@@MF&srlFJfWU^glgUv>ubJ|r3GLS6DCmq*;7F`Cgzlx9msWCZ}B zx^(2#3{TNIRCGn$RE=-Ze_MK-x<*OJE9lJfSB>k=1M^idZiBxLbb~VZldNlqh~AI- zRMW?{h;hAH1Gc%hqrAsAkG%_(IVeM~->M#U*^#cDb5e0YWv4EuS>bu(rlJ>QK?)X6 zU_l3};=JR;9;HNGsnF=$U96a)L=buhK2Z#TJc*@97B->@-IeFd<3!|1x_u<7?guQ zM@ED0+XFHXv%^bPf4!Qp(owm2+M{uS4RABHIRcd2R^{O=Gh3=U{8-YOn** zwYS=8|*=LXK zQg(IF87@Ts6M!WUKuj zF@*#k&|+?)$}!ABnI!A6o}0|FLFst2XZHcOhJ<~q6DWe~Mg{^VOAd3qocSkR8;x}{ z+&d5oTY-XB%N{P(y4-HtPzSS7>H0+PVfGfmY*ey7f8l$Wy+tq^m7rJa;yujX9+(Zv z_Um11*m~rgI$?&Yt-+`4nbpvq57n68dL(d(IWT%Ju>`pvqE)Cz3ocE1ijsW~=ISO; zS5u9CCxls@2CxbkB}xn(WsKgsx)bH?3U+mn8+1zmvu~g<`XSfSJ1ww_p|m$p*Ey1b zz`6wme_@Ei0)$;TkX7Sb8tvuUBj%XyM5^UcHv%L#NVT_yB|DGK>IyD=)6o4aSss+% zpMqBB<;5(7L~^;Gc!w<$bP3No8^mdE=bQ=c*`~aCT(z`jFSQOjgYxqv+#q?>NGLb+ zYI(+5+$`5^C7+P3awMD=S33l&N7|dDcL6#tki7bN`e|rE*`S2_8fDA%S-%<79phl3FThC1l=P9IFh=b1rNDxo zf0V?^8B}|^WLaJ!&U$i2U04sm1U(FBWTQltPB)hxGwQN}XoBS~WhaGr$zBeUp@WGuF2P5`{=eVtJ)KnXQ#Gti$2eq>6!st)G;c z@mV36n8;^}7-72|M;-14rQj#rwOd$syL%IyT0R14_xm?DkKC^X2c_c+eBN`bZwI%E zeozK~jk5RL>f6b!q92sgpK$k{f8%-^I4<;q68ICy-V^w5DS?lEPy&Cl{P!T-T7oIm zva%3Q;nETmy$)L=@iAV78IDg?aU<08H1j7DwQ)>qK!Cg zu33sOEX&@IpPDJ&wOER(RnkQo3{UOmh^W+9R-vBLZ zHBTUn!Bgak1R`>*m*jJto_mb)${8%NIR+Pk`hUq%;!ACY-3YKje`$0gY)A~-L2U-K zt^>Lj0Xd}kc(8fUB|QbUlIz3?eVE+*tS7XgzOW;=xpoGN&`QzLp!HhRGz%0`hj1H| zHm9bpo4q)*p28snE9Me2syyIXu+KN9w&-D9&hz$CCE!!ln%UrtJq<42G_*CTW;NZC z%H^P@hK1S~iG^)#e}oZ6CCf?nQeY;M{>$pi2kRO^4Lwv)1qZ5h5hI-JUx@{c4mR~w zSM`s&8P?q`^gF`qf3`uRhbctQ-)TYIpd7n-mwQ3na_bxobr<;#8CE%g?HgJM>RG!C>oOM1 zJd>arb$V`Se=dwTD1DEV8+H}D;nvU`(USf^u(Vk1>R5MaAeBb57gGm}L8*CE!Mq$8 zW;JysfGv)EKV&B)L>zqv^&lCNIcO0bM`2b=#ySm44{2ICgI;&Cb{9R~ACAevC!vN$ z16Dn=k6BnwSy@P&%OGNHma>Oj%6&$52C$ z@TtkM>qHZ{z=9k9miH7B%hpr-#EWNS3Dz;9i!ygMe}LRod=tYPC2*`*Eh-3ldG>m^nJ6V}f?EPKW9g`G-W*_|r&~6L4-v%iX4A_ac~v=smHkkAYfZ0v zPhvbn4vjz?lw$W3m#5n7Kwg0l5t6Z0wnN>peM4#yB@W>A{{Rt7Jg15imq1i=sg-R( zF+C=HjmgKaL5K70BHn{;kDwcrnOBdmr9N{f-rD+A_^j(4baB1n ze^k|*LI5KMn>ct?m^1o+?_m`POy|&^nNf%Et#a`d33{>WA~YYWOoCHDiqN7>thvN1`20Ck_1 zqb#KfW$6x+qa&35H;0(66ct@&d2Nn!fA8RiYpDaqpaeenhD~RLl-GCjE={Qb-f~UZ zugw)K_4#sem%=4_HiRl?u*#FsU+QXTLP-VCWdo}PXM+;>s^^vyVj0$7+6;FW-4IK* z(IHG>#TR93dIBA0$5g>Ie~op38I-D5z%*#Ae!#TIC2@?|lDRIro~qQ~9{qL0f1nm= z9NeZnimcKn`K|X_z%(Wue@;>@FF@LjC~UB_1)YRIo+L})R5D}^fvx$es$2*ZV$NNg zJ(V-uJD@4~FbFs0=g^o`J6PR@eMqOpSb|t8>Rmi6D44=!kE)O2teQXt6cqMNa$7B# zFeracEhJMvkd~)J1FO*=*;06de~7B*lCF%7qPkd7*WD_Zb!-gQjX(7h*T2~&bVCbW zoOwJ{@7u>YYV=E%s1%u@s7NANDooO5DGHG?Lza%IhHk9>TW5{r_%>S%Bgo>nTTYS*6BFM7By?EFo?1Dm$6M#fcD+Kv?LueYyS z3lAP zTG?LwusV%-U|OZA+;}}ljd)2kSBo=kdV+FioB6RMEAdR3JtcQZIUAFA%ir2g@)%YB zrm8S{l{EX54t-7W(tUw&YphFy#O3|hpWZQ(lzP%qgOc;@)<&6y%z2`TmLvK#hBVzm z5_YrLc#ICRsUuQ)ruk1^>lRs^uy9>q9Z4}DI|bmSzxnOOXf&N8&908QTVaN*O;FG4 z8YiV}4dR^^8Kt|iW;IV}Xxx}I=x{P5 z!&_{@3cGcPqVB(M!Ob3iRFKo6wV1B@sN?!wioDN=hv&rr&yERo&%yBh=W>1B+V^(n z#(p|vp?xJI)b7?i{5Kb=ytuDQe@~%L_-bd$i_vd6Xmy>|GrG3um$(#iW?amxkFa1aDpcpZ*UsCw`s*bOPxT5UajMu!j zEioDCa_VnPbTUrqXQbxL<80c;kGnkzd2UC@FMJa+vaXN&^5v4nC7hNDs}B*o@8~Ht zhub{tYi`KUY+l@Oszvn!`N~3x3nsU{B#RPbwzk2KDn^hUu}@c}cR?h%M5p$;T|=Cz z&UxI(rr*k;>)(IUG7TxAGZ?o6h5AzNFQ#XXIly{Ethb|4Z>WPyqlx4fHQW!aZ>K3^yNB;mdft=wc(D##=v};@T$igF z!c+IOJL?l7v8C}^D2YuUQ-yDc^hwS^*!6pnDgW5%>wH2m?k-A?6W*NknrKo}esj5p z*3V5oP~UF;v%b&QR%RZLG1w#T*EP0k2u^P)cQTDsy||a=_-5^_jogW1ADQq(+=Er-h%oE~ue(YlUdx=yrH6=yqC!Q~g@5^fECJ@xEyS{9n84M|S;VfIzWw1xd z&#eiPKNtVa>#&aZF_SHA8J+o0cHP~XdFc1ff8Lc=J2mxMVX-sSJLWsAU2#4E`7%zR zSX!&{L_%g9&+~QUGP~v^+E7E72o|=T3imz zShHdeXS7V;UQ^DyEU{TN`<~TSH}O0P)vcKOy$cE8=j-m^_@pyI~-ukIc2MD2_5TEYCL>)d^mLB2ovA@e{=>*Pg%PHr{ zI*Oxd`AWU(@K!U#>T$6dgvyBcjl+H`!%ji#XG*Gic-h~kJaT!J z`)pXt&>{EUspr*aG#&rBW*S}6e_{Qm9e5qINNZq0kZ9xQof?Ku$4(lCmj%`;Qe1uV zi;^2&Pd+h=&+xT%3K|?#vN_UJ7v~fYnkbXHy;qa*r!?)=pDJ4zfWt`Yz<@!BqHl=1 z$-Zf&dmebXyoeB;gD#3by!C8_g2YW2+K7+}%io^JUAEC(*Gbz3yi6h`&$kE7{Aen8 zaXjkF4_iqagd6ohx`oY(MwZ#;PDnfBw??w~(D9!K_rH+lz446Kn0j1r?;iYq=#^>D zx9Z8I3%|ByKWzVMp2}!g&6GF|-XFYyRLbKZ6rX8j$9UXNOWNLIv90I(eUIg?nA@O; zpDj^l@!f(=m!@0NWm0|P`9r4*TP5)%mY4BBs?DmaUdtWb(^k3o%M(|tdu&!0)w8w- zeh34w0|A%g)6K{2uYBetzNg%GMf!0Uvs>xk# ztv8%}-dg7M`={>ek>Lvk!E(t2y`*$3vT_0GS|MC)T(txo@;)mRAz6Gy2)>Mr5FUsm z-exP-SSN{he2`7jnRagw9q@_qwlInUE;l2_B9cZBpAqiAiHtLjx}SRYI+|(rsAiU2 z(s}ng`@q6^)-Cv&<*fUEhQJXE>2ekS4C5RwP;jSb7o(}n-2b6zm*r1Sk9W(+t2OKT zQ(vsC)9}qkEx{yWhfF3O{AsrQt7$ZjzUgY*A&yEo@nGUwgk*tZtd|9XOdA{~iIk^KCsaXcZ~DBSr5H$k@hypevhnIqNnhZJ<czLR0#i_w>G`RnId-jMeB#>Q?J(VEX+g?6QW<# zgDVsd?w?T|l;#CY#1a8x!29!y+4-H^p-?rSgG(edIDpf@TOA-}fGV9G%fQ+mc|1Mm_7wJH0TlS8{ZTII@wbQduv?T*6(MsgN-6d-`Zd zZEY+WYhW0@`ZY(Ru`hLhROHg&E0lW#h)E&;VScDvo5$k{(dBk?DEYGa~v;>k! z19v5Xgho!)-V6mnpK&e;e<-?$4Sty8utJ8Y(*YT1e4^<{Lg#%?#!`nH0U>1Ti(wbtha1SoMziR{k#UU>JCgcmlhN9y1x&U zz3-awta=K4r)Ify7_QN2tUB&@wj^?(YO<`hJSc)yH&TPLBPU19r;LDLBz*c@#ED ztZnVyee5n{qXu>;jaO|{k=R)k<)^M&}Eo|X3> z6xnd_V;lE_(prJOCKgo(uQA&EsWnNW`9>l^A)0N$a=s}?kT(0-CcncCqA)o>C6O3y%_d?ALz>#WuY1f{@mF&) zsxuvzf!pkoIqPjW+foEB`5z^Yy^5Nz;Yn6}umQS61-00q_^SX0kWs*EDxjdstAG?< z4Jow=PQw)v=6?yvBVk1zGIIpQ?5$fRNC?>AQbR zjlq^Za5nFbzfS(wS#Kns87U9xW8jUMApBbfFBh{BY6$=~Lr*Zk4&iy|J*>&|Yc_3w z%M2V8`~MztASMPl4NTs{uie`LZM_CaLD#T=oUlC10_f$#H0dwk75X+n``zF%i3|wTKsQ*O9}++lA8aWCkaCB?dJ@vY z!9#a;2Q-1RbO98>AgmAla);wM0BPdjIFR#vX(!%lUi9yd!|!2BK`XAnMgYm?yMqX%RtqR~iA6$eWC2v3;>i2-Vj}Q^SdlBIckj8a*+{p;OWhWGR z9d4^EDS#F}RLBZoizgtv^WXYA{G4yWM=<9~cMF=ic*5H(?h`;UA0`h7K-de0S6TwF z;)9!x00=M$Vq`qENRF93T2R-y#QX^@YO^jD&9a!NI!u!u}Ak0x02wX}kc0{b1O4 zM*x;Ez{e%V51x{xc!(JQ??##*+|6RP0QCJ~5G@tJ9X`~R|AnU&aQ{R8aR2mHfvV94 zTT}yJOYAp+#XkUUztSTB^mkY;cZXL9DF9v}@1WTLxYOBLf$9(lH*J^`KsF2_jiLZP zplk=U5D08n__xr9kdWek9h{mif@W~(j}4G{5FjQ0f9d)qfYA2;Y8wWH1_ApKIX3Xw HK> #+begin_src matlab - initializeReferences('Dy_type', 'triangular', 'Dy_amplitude', 10e-3, 'Dy_period', 1); - sim(mdl); - ty_scan_triangle = simout; - % save('./mat/experiment_tomography.mat', 'ty_scan_triangle', '-append'); +initializeReferences('Dy_type', 'triangular', 'Dy_amplitude', 10e-3, 'Dy_period', 1); +sim(mdl); +ty_scan_triangle = simout; +% save('./mat/experiment_tomography.mat', 'ty_scan_triangle', '-append'); #+end_src * Conclusion <> - * Bibliography :ignore: #+latex: \printbibliography[heading=bibintoc,title={Bibliography}] @@ -2979,200 +2980,228 @@ arguments end #+end_src -*** Load Data #+begin_src matlab -load('dist_psd.mat', 'dist_f'); +% Initialization of random numbers +rng("shuffle"); #+end_src +*** Ground Motion +#+begin_src matlab +%% Ground Motion +load('dist_psd.mat', 'dist_f'); + +% Frequency Data +Dw.f = dist_f.f(2:end); +Dw.psd_x = dist_f.psd_gm(2:end); +Dw.psd_y = dist_f.psd_gm(2:end); +Dw.psd_z = dist_f.psd_gm(2:end); + +% Time data +Fs = 2*Dw.f(end); % Sampling Frequency of data is twice the maximum frequency of the PSD vector [Hz] +N = 2*length(Dw.f); % Number of Samples match the one of the wanted PSD +T0 = N/Fs; % Signal Duration [s] +Dw.t = linspace(0, T0, N+1)'; % Time Vector [s] + +C = zeros(N/2,1); +for i = 1:N/2 + C(i) = sqrt(Dw.psd_x(i)/T0); +end + +if args.Dwx && args.enable + theta = 2*pi*rand(N/2,1); % Generate random phase [rad] + Cx = [0 ; C.*complex(cos(theta),sin(theta))]; + Cx = [Cx; flipud(conj(Cx(2:end)))];; + Dw.x = N/sqrt(2)*ifft(Cx); % Ground Motion - x direction [m] +else + Dw.x = zeros(length(Dw.t), 1); +end + +if args.Dwy && args.enable + theta = 2*pi*rand(N/2,1); % Generate random phase [rad] + Cx = [0 ; C.*complex(cos(theta),sin(theta))]; + Cx = [Cx; flipud(conj(Cx(2:end)))];; + Dw.y = N/sqrt(2)*ifft(Cx); % Ground Motion - y direction [m] +else + Dw.y = zeros(length(Dw.t), 1); +end + +if args.Dwy && args.enable + theta = 2*pi*rand(N/2,1); % Generate random phase [rad] + Cx = [0 ; C.*complex(cos(theta),sin(theta))]; + Cx = [Cx; flipud(conj(Cx(2:end)))];; + Dw.z = N/sqrt(2)*ifft(Cx); % Ground Motion - z direction [m] +else + Dw.z = zeros(length(Dw.t), 1); +end +#+end_src + +*** Translation Stage + + We remove the first frequency point that usually is very large. #+begin_src matlab :exports none +load('dist_psd.mat', 'dist_f'); dist_f.f = dist_f.f(2:end); dist_f.psd_gm = dist_f.psd_gm(2:end); dist_f.psd_ty = dist_f.psd_ty(2:end); dist_f.psd_rz = dist_f.psd_rz(2:end); #+end_src -*** Parameters -We define some parameters that will be used in the algorithm. #+begin_src matlab -Fs = 2*dist_f.f(end); % Sampling Frequency of data is twice the maximum frequency of the PSD vector [Hz] -N = 2*length(dist_f.f); % Number of Samples match the one of the wanted PSD -T0 = N/Fs; % Signal Duration [s] -df = 1/T0; % Frequency resolution of the DFT [Hz] - % Also equal to (dist_f.f(2)-dist_f.f(1)) -t = linspace(0, T0, N+1)'; % Time Vector [s] -Ts = 1/Fs; % Sampling Time [s] -#+end_src +%% Translation Stage +load('dist_psd.mat', 'dist_f'); + +% Frequency Data +Ty.f = dist_f.f(2:end); +Ty.psd_x = dist_f.psd_ty(2:end); % TODO - we take here the vertical direction which is wrong but approximate +Ty.psd_z = dist_f.psd_ty(2:end); + +% Time data +Fs = 2*Ty.f(end); % Sampling Frequency of data is twice the maximum frequency of the PSD vector [Hz] +N = 2*length(Ty.f); % Number of Samples match the one of the wanted PSD +T0 = N/Fs; % Signal Duration [s] +Ty.t = linspace(0, T0, N+1)'; % Time Vector [s] -*** Ground Motion -#+begin_src matlab -phi = dist_f.psd_gm; C = zeros(N/2,1); for i = 1:N/2 - C(i) = sqrt(phi(i)*df); + C(i) = sqrt(Ty.psd_x(i)/T0); end -#+end_src -#+begin_src matlab -if args.Dwx && args.enable - rng(111); - theta = 2*pi*rand(N/2,1); % Generate random phase [rad] - Cx = [0 ; C.*complex(cos(theta),sin(theta))]; - Cx = [Cx; flipud(conj(Cx(2:end)))];; - Dwx = N/sqrt(2)*ifft(Cx); % Ground Motion - x direction [m] -else - Dwx = zeros(length(t), 1); -end -#+end_src - -#+begin_src matlab -if args.Dwy && args.enable - rng(112); - theta = 2*pi*rand(N/2,1); % Generate random phase [rad] - Cx = [0 ; C.*complex(cos(theta),sin(theta))]; - Cx = [Cx; flipud(conj(Cx(2:end)))];; - Dwy = N/sqrt(2)*ifft(Cx); % Ground Motion - y direction [m] -else - Dwy = zeros(length(t), 1); -end -#+end_src - -#+begin_src matlab -if args.Dwy && args.enable - rng(113); - theta = 2*pi*rand(N/2,1); % Generate random phase [rad] - Cx = [0 ; C.*complex(cos(theta),sin(theta))]; - Cx = [Cx; flipud(conj(Cx(2:end)))];; - Dwz = N/sqrt(2)*ifft(Cx); % Ground Motion - z direction [m] -else - Dwz = zeros(length(t), 1); -end -#+end_src - -*** Translation Stage - X direction -#+begin_src matlab +% Translation Stage - X if args.Fty_x && args.enable - phi = dist_f.psd_ty; % TODO - we take here the vertical direction which is wrong but approximate + phi = Ty.psd_x; C = zeros(N/2,1); for i = 1:N/2 - C(i) = sqrt(phi(i)*df); + C(i) = sqrt(phi(i)/T0); end - rng(121); theta = 2*pi*rand(N/2,1); % Generate random phase [rad] Cx = [0 ; C.*complex(cos(theta),sin(theta))]; Cx = [Cx; flipud(conj(Cx(2:end)))];; u = N/sqrt(2)*ifft(Cx); % Disturbance Force Ty x [N] - Fty_x = u; + Ty.x = u; else - Fty_x = zeros(length(t), 1); + Ty.x = zeros(length(Ty.t), 1); end -#+end_src -*** Translation Stage - Z direction -#+begin_src matlab +% Translation Stage - Z if args.Fty_z && args.enable - phi = dist_f.psd_ty; + phi = Ty.psd_z; C = zeros(N/2,1); for i = 1:N/2 - C(i) = sqrt(phi(i)*df); + C(i) = sqrt(phi(i)/T0); end - rng(122); theta = 2*pi*rand(N/2,1); % Generate random phase [rad] Cx = [0 ; C.*complex(cos(theta),sin(theta))]; Cx = [Cx; flipud(conj(Cx(2:end)))];; u = N/sqrt(2)*ifft(Cx); % Disturbance Force Ty z [N] - Fty_z = u; + Ty.z = u; else - Fty_z = zeros(length(t), 1); + Ty.z = zeros(length(Ty.t), 1); end #+end_src -*** Spindle - X direction +*** Spindle #+begin_src matlab -% if args.Frz_x && args.enable -% phi = dist_f.psd_rz; -% C = zeros(N/2,1); -% for i = 1:N/2 -% C(i) = sqrt(phi(i)*df); -% end -% rng(131); -% theta = 2*pi*rand(N/2,1); % Generate random phase [rad] -% Cx = [0 ; C.*complex(cos(theta),sin(theta))]; -% Cx = [Cx; flipud(conj(Cx(2:end)))];; -% u = N/sqrt(2)*ifft(Cx); % Disturbance Force Rz z [N] -% Frz_x = u; -% else -Frz_x = zeros(length(t), 1); -% end -#+end_src +%% Translation Stage +load('dist_psd.mat', 'dist_f'); -*** Spindle - Y direction -#+begin_src matlab -% if args.Frz_y && args.enable -% phi = dist_f.psd_rz; -% C = zeros(N/2,1); -% for i = 1:N/2 -% C(i) = sqrt(phi(i)*df); -% end -% rng(131); -% theta = 2*pi*rand(N/2,1); % Generate random phase [rad] -% Cx = [0 ; C.*complex(cos(theta),sin(theta))]; -% Cx = [Cx; flipud(conj(Cx(2:end)))];; -% u = N/sqrt(2)*ifft(Cx); % Disturbance Force Rz z [N] -% Frz_z = u; -% else -Frz_y = zeros(length(t), 1); -% end -#+end_src +% Frequency Data +Rz.f = dist_f.f(2:end); +Rz.psd_x = dist_f.psd_rz(2:end); % TODO - we take here the vertical direction which is wrong but approximate +Rz.psd_y = dist_f.psd_rz(2:end); % TODO - we take here the vertical direction which is wrong but approximate +Rz.psd_z = dist_f.psd_rz(2:end); -*** Spindle - Z direction -#+begin_src matlab -if args.Frz_z && args.enable - phi = dist_f.psd_rz; +% Time data +Fs = 2*Rz.f(end); % Sampling Frequency of data is twice the maximum frequency of the PSD vector [Hz] +N = 2*length(Rz.f); % Number of Samples match the one of the wanted PSD +T0 = N/Fs; % Signal Duration [s] +Rz.t = linspace(0, T0, N+1)'; % Time Vector [s] + +C = zeros(N/2,1); +for i = 1:N/2 + C(i) = sqrt(Rz.psd_x(i)/T0); +end + +% Translation Stage - X +if args.Frz_x && args.enable + phi = Rz.psd_x; C = zeros(N/2,1); for i = 1:N/2 - C(i) = sqrt(phi(i)*df); + C(i) = sqrt(phi(i)/T0); + end + theta = 2*pi*rand(N/2,1); % Generate random phase [rad] + Cx = [0 ; C.*complex(cos(theta),sin(theta))]; + Cx = [Cx; flipud(conj(Cx(2:end)))];; + u = N/sqrt(2)*ifft(Cx); % Disturbance Force Rz x [N] + Rz.x = u; +else + Rz.x = zeros(length(Rz.t), 1); +end + +% Translation Stage - Y +if args.Frz_y && args.enable + phi = Rz.psd_y; + C = zeros(N/2,1); + for i = 1:N/2 + C(i) = sqrt(phi(i)/T0); + end + theta = 2*pi*rand(N/2,1); % Generate random phase [rad] + Cx = [0 ; C.*complex(cos(theta),sin(theta))]; + Cx = [Cx; flipud(conj(Cx(2:end)))];; + u = N/sqrt(2)*ifft(Cx); % Disturbance Force Rz y [N] + Rz.y = u; +else + Rz.y = zeros(length(Rz.t), 1); +end + +% Translation Stage - Z +if args.Frz_z && args.enable + phi = Rz.psd_z; + C = zeros(N/2,1); + for i = 1:N/2 + C(i) = sqrt(phi(i)/T0); end - rng(131); theta = 2*pi*rand(N/2,1); % Generate random phase [rad] Cx = [0 ; C.*complex(cos(theta),sin(theta))]; Cx = [Cx; flipud(conj(Cx(2:end)))];; u = N/sqrt(2)*ifft(Cx); % Disturbance Force Rz z [N] - Frz_z = u; + Rz.z = u; else - Frz_z = zeros(length(t), 1); + Rz.z = zeros(length(Rz.t), 1); end #+end_src *** Direct Forces #+begin_src matlab -u = zeros(length(t), 6); +u = zeros(100, 6); Fd = u; #+end_src *** Set initial value to zero #+begin_src matlab -Dwx = Dwx - Dwx(1); -Dwy = Dwy - Dwy(1); -Dwz = Dwz - Dwz(1); -Fty_x = Fty_x - Fty_x(1); -Fty_z = Fty_z - Fty_z(1); -Frz_z = Frz_z - Frz_z(1); +Dw.x = Dw.x - Dw.x(1); +Dw.y = Dw.y - Dw.y(1); +Dw.z = Dw.z - Dw.z(1); +Ty.x = Ty.x - Ty.x(1); +Ty.z = Ty.z - Ty.z(1); +Rz.x = Rz.x - Rz.x(1); +Rz.y = Rz.y - Rz.y(1); +Rz.z = Rz.z - Rz.z(1); #+end_src *** Save the Structure #+begin_src matlab if exist('./mat', 'dir') if exist('./mat/nass_disturbances.mat', 'file') - save('mat/nass_disturbances.mat', 'Dwx', 'Dwy', 'Dwz', 'Fty_x', 'Fty_z', 'Frz_x', 'Frz_y', 'Frz_z', 'Fd', 'Ts', 't', 'args', '-append'); + save('mat/nass_disturbances.mat', 'Dw', 'Ty', 'Rz', 'Fd', 'args', '-append'); else - save('mat/nass_disturbances.mat', 'Dwx', 'Dwy', 'Dwz', 'Fty_x', 'Fty_z', 'Frz_x', 'Frz_y', 'Frz_z', 'Fd', 'Ts', 't', 'args'); + save('mat/nass_disturbances.mat', 'Dw', 'Ty', 'Rz', 'Fd', 'args'); end elseif exist('./matlab', 'dir') if exist('./matlab/mat/nass_disturbances.mat', 'file') - save('matlab/mat/nass_disturbances.mat', 'Dwx', 'Dwy', 'Dwz', 'Fty_x', 'Fty_z', 'Frz_x', 'Frz_y', 'Frz_z', 'Fd', 'Ts', 't', 'args', '-append'); + save('matlab/mat/nass_disturbances.mat', 'Dw', 'Ty', 'Rz', 'Fd', 'args', '-append'); else - save('matlab/mat/nass_disturbances.mat', 'Dwx', 'Dwy', 'Dwz', 'Fty_x', 'Fty_z', 'Frz_x', 'Frz_y', 'Frz_z', 'Fd', 'Ts', 't', 'args'); + save('matlab/mat/nass_disturbances.mat', 'Dw', 'Ty', 'Rz', 'Fd', 'args'); end end #+end_src