From f0821c0f864fb59daa26a8b37380d778fd278692 Mon Sep 17 00:00:00 2001 From: thetazero Date: Wed, 22 Jun 2022 12:40:03 -0700 Subject: [PATCH 01/20] Initial work on KF --- Simulator/kf.jl | 11 +++++++++++ Simulator/simulator.jl | 35 +++++++++++++++++++++++++++++++++++ Simulator/tests/tests.jl | 22 +++++++++++++++++----- 3 files changed, 63 insertions(+), 5 deletions(-) create mode 100644 Simulator/kf.jl diff --git a/Simulator/kf.jl b/Simulator/kf.jl new file mode 100644 index 0000000..09504b0 --- /dev/null +++ b/Simulator/kf.jl @@ -0,0 +1,11 @@ +using Quaternions + +struct EKF + q::Quaternion + P::Matrix{Float64} +end + +function step(e::EKF, qn::Tuple{Float64,Float64,Float64}) + e.q = qn + return e.q +end \ No newline at end of file diff --git a/Simulator/simulator.jl b/Simulator/simulator.jl index 71400ed..0c79e7d 100644 --- a/Simulator/simulator.jl +++ b/Simulator/simulator.jl @@ -9,6 +9,7 @@ module Simulator using SatelliteDynamics, LinearAlgebra include("mag_field.jl") include("quaternions.jl") + include("kf.jl") export initialize_orbit, intialize_orbit_oe # Generate valid initial conditions for a satellite @@ -255,4 +256,38 @@ module Simulator return q_hist end + + function my_sim_kf(control_fn) + x₀ = initialize_orbit() + println("intialized orbit!") + # x₀[11:13] .=0 + # x₀[11:13] *= x₀[11:13].*10.0 # Spinning very fast + + J = [0.3 0 0; 0 0.3 0; 0 0 0.3] # Arbitrary inertia matrix for the Satellite + t = Epoch(2020, 11, 30) # Starting time is Nov 30, 2020 + dt = 0.5 # Time step, in seconds + + N = 1000000 + q_hist = zeros(N, 8) + q_hist[1, 1:3] .= x₀[11:13] + q_hist[1, 4] = norm(x₀[11:13]) + x = x₀ + for i = 1:N - 1 + r, v, q, ω = x[1:3], x[4:6], x[7:10], x[11:13] + b = IGRF13(r, t) + x = rk4(x, J, control_fn(ω, b), t, dt) + t += dt # Don't forget to update time (not that it really matters...) + # q_hist[i + 1, :] .= x[7:10] + q_hist[i + 1, 1:3] .= x[11:13] + ω = norm(x[11:13]) + if ω < 0.001 + q_hist = q_hist[1:i, :] + break + end + q_hist[i + 1, 4] = ω + end + + return q_hist + + end end \ No newline at end of file diff --git a/Simulator/tests/tests.jl b/Simulator/tests/tests.jl index bb26472..9abac74 100644 --- a/Simulator/tests/tests.jl +++ b/Simulator/tests/tests.jl @@ -106,16 +106,28 @@ end display(plot(e_hist, title ="Energy")) end +@testset "MEKF/DeTumbling" begin + function control_law(ω, b) + b̂ = b / norm(b) + k = 7e-4 + M = -k * (I(3) - b̂*b̂')*ω + m = 1 / (dot(b, b)) * cross(b, M) + return cross(m, b) + end + + data = my_sim(control_law) + display(plot(data, title="MEKF/DeTumbling", xlabel="Time (s)", ylabel="Angular Velocity (rad/s)", labels=["ω1" "ω2" "ω3" "ω"])) + +end + @testset "DeTumbling" begin function control_law(ω, b) b̂ = b / norm(b) k = 7e-4 - # M = -k * (identity(3) - dot(b̂,b̂))*ω M = -k * (I(3) - b̂*b̂')*ω - # print(M, ω, b) - # print(dot(ω,ω),'\n') - return M + m = 1 / (dot(b, b)) * cross(b, M) + return cross(m, b) end data = my_sim(control_law) @@ -162,6 +174,6 @@ end end data = my_sim(control_law) - display(plot(data, title="DeTumbleIO", xlabel="Time (s)", ylabel="Angular Velocity (rad/s)")) + display(plot(data, title="DeTumbleIO", xlabel="Time (s)", ylabel="Angular Velocity (rad/s)", labels=["ω1" "ω2" "ω3" "ω"])) end From 259d751e92fcd653d96f4eeea058d0979ca7493c Mon Sep 17 00:00:00 2001 From: thetazero Date: Wed, 22 Jun 2022 13:02:14 -0700 Subject: [PATCH 02/20] Minor kf integration progress --- Simulator/kf.jl | 6 ++++-- Simulator/simulator.jl | 22 ++++++++++++++-------- Simulator/tests/tests.jl | 2 +- 3 files changed, 19 insertions(+), 11 deletions(-) diff --git a/Simulator/kf.jl b/Simulator/kf.jl index 09504b0..ca1ddd6 100644 --- a/Simulator/kf.jl +++ b/Simulator/kf.jl @@ -1,11 +1,13 @@ +# [Simulator/kf.jl] + using Quaternions struct EKF - q::Quaternion + q::Vector{Float64} P::Matrix{Float64} end -function step(e::EKF, qn::Tuple{Float64,Float64,Float64}) +function step(e::EKF, qn::Vector{Float64}) e.q = qn return e.q end \ No newline at end of file diff --git a/Simulator/simulator.jl b/Simulator/simulator.jl index 0c79e7d..70f8999 100644 --- a/Simulator/simulator.jl +++ b/Simulator/simulator.jl @@ -6,7 +6,7 @@ """ module Simulator - using SatelliteDynamics, LinearAlgebra + using SatelliteDynamics, LinearAlgebra include("mag_field.jl") include("quaternions.jl") include("kf.jl") @@ -16,6 +16,7 @@ module Simulator export rk4 # Interface function for updating state export IGRF13 # Function call that provides magnetic field vector given position and time export my_sim + export my_sim_kf """ @@ -268,23 +269,28 @@ module Simulator dt = 0.5 # Time step, in seconds N = 1000000 - q_hist = zeros(N, 8) - q_hist[1, 1:3] .= x₀[11:13] - q_hist[1, 4] = norm(x₀[11:13]) + q_hist = zeros(N, 9) + q_hist[1, 1] = norm(x₀[11:13]) x = x₀ + + q₀ = x₀[7:10] + P = Diagonal([1, 1, 1, 1]) + kf = EKF(q₀, P) + for i = 1:N - 1 r, v, q, ω = x[1:3], x[4:6], x[7:10], x[11:13] b = IGRF13(r, t) x = rk4(x, J, control_fn(ω, b), t, dt) t += dt # Don't forget to update time (not that it really matters...) - # q_hist[i + 1, :] .= x[7:10] - q_hist[i + 1, 1:3] .= x[11:13] - ω = norm(x[11:13]) + q_hist[i + 1, 1] = norm(ω) + q_hist[i + 1, 2:5] .= q + + q_hist[i+1, 6:9] .= step(kf, q) + if ω < 0.001 q_hist = q_hist[1:i, :] break end - q_hist[i + 1, 4] = ω end return q_hist diff --git a/Simulator/tests/tests.jl b/Simulator/tests/tests.jl index 9abac74..94fe422 100644 --- a/Simulator/tests/tests.jl +++ b/Simulator/tests/tests.jl @@ -115,7 +115,7 @@ end return cross(m, b) end - data = my_sim(control_law) + data = my_sim_kf(control_law) display(plot(data, title="MEKF/DeTumbling", xlabel="Time (s)", ylabel="Angular Velocity (rad/s)", labels=["ω1" "ω2" "ω3" "ω"])) end From 332c2f1594acc86f4d7391a24e232d9622db72ab Mon Sep 17 00:00:00 2001 From: thetazero Date: Thu, 23 Jun 2022 11:28:11 -0700 Subject: [PATCH 03/20] finally implementing the kf --- Simulator/kf.jl | 52 ++++++++++++++++++++++++++++++++++++--- Simulator/mocksat.js | 35 -------------------------- Simulator/simulator.jl | 6 ++--- Simulator/tests/tests.jl | 2 +- summary/main.pdf | Bin 0 -> 83787 bytes summary/main.tex | 23 +++++++++++++++++ 6 files changed, 75 insertions(+), 43 deletions(-) delete mode 100644 Simulator/mocksat.js create mode 100644 summary/main.pdf create mode 100644 summary/main.tex diff --git a/Simulator/kf.jl b/Simulator/kf.jl index ca1ddd6..0e0f94e 100644 --- a/Simulator/kf.jl +++ b/Simulator/kf.jl @@ -1,13 +1,57 @@ # [Simulator/kf.jl] -using Quaternions +using LinearAlgebra +include("quaternions.jl") -struct EKF +mutable struct EKF q::Vector{Float64} + β::Vector{Float64} P::Matrix{Float64} end -function step(e::EKF, qn::Vector{Float64}) - e.q = qn +function hat(ω::Vector) + return [0 ω[3] -ω[2] + -ω[3] 0 ω[1] + ω[2] -ω[1] 0] +end + +function f( + q::Vector{Float64}, + β::Vector{Float64}, + ω::Vector{Float64}, + δt::Float64 +) + θ = norm(ω - β) * δt + r = (ω - β) / norm(ω - e.β) + return [L(q) * [r .* sin(θ / 2); cos(θ / 2)], β] +end + +function step( + e::EKF, + ω::Vector{Float64}, + δt::Float64, + ⁿr_mag::Vector{Float64}, + ⁿr_sun::Vector{Float64}, + ᵇr_mag::Vector{Float64}, + ᵇr_sun::Vector{Float64} +) + # Predict: + xₚ = f(e.q, e.β, ω, δt) + R = exp(-hat(ω - e.β) * δt) + A = [ + R (-δt*I(3)) + zeros(3, 3) I(3) + ] + W = zeros(6, 6) # related to some noise or something (ask Zac) + Pₚ = A * P * A' + W + # Innovation + Z = [ᵇr_mag; ᵇr_sun] - [Q zeros(3, 3); zeros(3, 3) Q] * [ⁿr_mag; ⁿr_sun] + C = [hat(ᵇr_mag) zeroes(3, 3); hat(ᵇr_sun) zeros(3, 3)] + V = zeros(3, 3) # Something else + S = C * Pₚ * C' + V + # Kalman Gain + L = Pₚ * C' * inv(S) + # Update + δx = L * Z return e.q end \ No newline at end of file diff --git a/Simulator/mocksat.js b/Simulator/mocksat.js deleted file mode 100644 index 14ecd39..0000000 --- a/Simulator/mocksat.js +++ /dev/null @@ -1,35 +0,0 @@ -const readline = require('readline') - -// while(true){ -// console.log('hello world') -// } -console.log('hello world, from mocksatjs') - -let ω = [] -let b = [] - -const rl = readline.createInterface({ - input: process.stdin, - output: process.stdout -}) -rl.on('line', input => { - if (input.length < 4 || input.slice(0, 3) != ">>>") return - input = input.slice(3) - let key = input[0] - let value = input.slice(1) - if (key == 'ω') { - ω = JSON.parse(value) - } else if (key == 'b') { - b = JSON.parse(value) - } else if (key == '?') { - let control = ω.map(v => -v * 0.001) - control = JSON.stringify(control) - console.log(`>>>M${control}`) - } - // console.log(ω, b) - // console.log(key, value) -}) - -rl.on('close', e => { - console.log(`eee${e}`) -}) diff --git a/Simulator/simulator.jl b/Simulator/simulator.jl index 70f8999..12c48d3 100644 --- a/Simulator/simulator.jl +++ b/Simulator/simulator.jl @@ -247,11 +247,11 @@ module Simulator # q_hist[i + 1, :] .= x[7:10] q_hist[i + 1, 1:3] .= x[11:13] ω = norm(x[11:13]) + q_hist[i + 1, 4] = ω if ω < 0.001 q_hist = q_hist[1:i, :] break end - q_hist[i + 1, 4] = ω end return q_hist @@ -285,9 +285,9 @@ module Simulator q_hist[i + 1, 1] = norm(ω) q_hist[i + 1, 2:5] .= q - q_hist[i+1, 6:9] .= step(kf, q) + q_hist[i+1, 6:9] .= step(kf, q, ω, dt) - if ω < 0.001 + if norm(ω) < 0.1 q_hist = q_hist[1:i, :] break end diff --git a/Simulator/tests/tests.jl b/Simulator/tests/tests.jl index 94fe422..7e620f3 100644 --- a/Simulator/tests/tests.jl +++ b/Simulator/tests/tests.jl @@ -116,7 +116,7 @@ end end data = my_sim_kf(control_law) - display(plot(data, title="MEKF/DeTumbling", xlabel="Time (s)", ylabel="Angular Velocity (rad/s)", labels=["ω1" "ω2" "ω3" "ω"])) + display(plot(data, title="MEKF/DeTumbling", xlabel="Time (s)", ylabel="Angular Velocity (rad/s)", labels=["ω" "i1" "i2" "i3" "r" "i1'" "i2'" "i3'" "r'"])) end diff --git a/summary/main.pdf b/summary/main.pdf new file mode 100644 index 0000000000000000000000000000000000000000..3e72e95a0939a0ab0d390839f0056f62a59f9d62 GIT binary patch literal 83787 zcma&MQ;aTv5^mYHd$(=dwr$(CcH6dXTf1%Bwr%V0doE5g$=qb7p1w-zb){03MW!e! zPRB&g217Q#JhBeMLdZzyU~CP;%L~IGYi4iZYDvh+$w>IW3x+}5%GTA)nUF!;*2vXN z)Xc=e)C`87AI8Pi+04ie#&a`QUCv>f9jWK8-q0Qesqy&!YhEOnT@v(<85 zL?wTYo*mrp7o3^ccKP!Fyw^#0VN5W+&CYnNM&W z54KQ!cCR{&_+g#n0(N^S34b`L+6(0cc|AC9ee|Jniw$721(cno2A>MEH9HCrl}U?ybQpT8U=ewhPl(Ab>W9aI7WTm6C6eYetfX73Iy za1QAlc;I&83)eWS2eW3}U-a|?CO${U2<-Mg=^=5TV6ywdJL>1h+WoG0`WV8ozQ}0> zVvf0D%LqU# zC&~TxK6^(BeiStp9a~3Hl*);74b=dx?B_v_=VtS@aCn#4mF>tM(pRRF>MT+v5HGiRW@ zMrNjBUgnsEL85H3V1-Dx#APu;hcS*c#k|I)C<0`{1z29jG_WK?oCI6vSigw19NJ() zvsbtC&)78F#sS*-mbQa{PA7F$pRLa{OuI97TIZwYfJ6%ock9=EiU;N$_g~f}@M#)XC4FTgk{dybi>s+CM%+`wx z9qoZRIlR|7Iop`7uAPb+n8yOnaNNET&9KYB^xq7ySQb#@FG}RF&c3XwE&-saJBEVv z_=FgJjOhW46ch>B$tvdpVFf?5eY%P1A22&(P@BD##{rntC0Ifbk2b@5(+Cn!@e z7T}JM5HQVOM@Z=nk(@$9B-{!J=hnb=0dufmLqR-9#y^LbF9t-IVs^TEc%fXJyu7^Y zLa+4K2(zn)vNeK;71i*HIENrk?%|ps{%vrokQ*WZ>sY*CgbEu_w-1%VBg2zRVCOI( z0&qr3Ob`z*SEpbd;hezujX-YcDnWP&sBZ<66#;n7Kz}@hgy`I>J4b)7Uj(qw-vc<1 zV7zNQAx8$vZEQ>MP@w?q$ExePyt*1er0vHxBB;~D`*TERNFgJ_Xo+9}X-MFc(w2aN zF@)dCIq+*ZH&-_!kHA5{xA6y_=_9780Xj3>>l-3R4UCw*mkHrqqWv+Kx5M`<=bL>5 ze8>D>LK@2L#Y+#IH2XN*zm65jIHFRo4A~OMqXNt}5i|+$$;m1K22fB=Ai}%@?6`fX z=YP3`0gBE)=mM|%XGiDzP%VFLK@TBYgLwT6-`fMWi3Bw)I3n`r`p|x38Qa?f(i9tL z=0S|%LruO>f0W=FKT!pC`&gDyEvG*u1?UOS?$-7c%HrXgT?2-?0ZhL~QCxpwj#X7p zv;fv04y7$E!4&i1vGI`OqcdaR$A^c=ATak%AU^<}IC7{@4ZNJ+cj|#!+yG(!v@cKj zU#g8C1_VX`8=-9fZ(Ql^{kSev!UI83w&SBy5W!wMApZgoKJ`m5 zBVQ6DzJzF63+oNA`!Aj-bnM;<1q&S@$w6YA->ra{9vvQk6YLzh2HxPWg-?-C{3;6J zhoAMSOb!enr1h1@WXBUgIlDu689DMRM`p((o)5in#-$<6324F{8;0|=^LGLfWIusw zb8(t{_bU1E6NWE4EaXcfjT3!wpzK39QvzbaxB_vrehVaWUj-r{#|xi0P*8xJ0{}w8 zf(Z9}5$SKpeDb#fpL=>PngCd5U@pSnp!mod05}1Q)-No7Sn~(Py9o=sejIze%oRbj z1lF^`eWS>F-@b{{>#tzlm@?qxi?b3J;aAkZ#drJUS6LwK^qt{dGy3U2BgX0eLVkqd zSAF|OSRFs09Rw_1Kfu3N?gRF-v++E)@gLsZpTGM8=+t0`(7+}>nJz%lA)BkA?z+9{ ze{To=@sYF3wb>=i4xc`abS&LqL$ak-lJj`3eCSfURUkT8sl+|VIwwBPTGY1!vlY>4 zRsW2>t@0{)J9puk7zJ4xe(&U|`U-$^b7j{Z$KB;7*ZU>lDbp0_WYW(~$RSBI?{Acj zNr2CN0v_v04qUdhdDyK#jFdY$i9nR+DDjc1~ z2`ExiG#S4C#GFt)GHD;hb!ENi*Das8Ul5g&u@M&vat{$?mb@hb(#m9c?^LvL>v8c_ zY{iV&*JyVmVVw_aOD_S$+Wfq|{IZWbypmT)(UI6gLpg=#OU#hUM zl5n7jVQ~9sAB`h@f&_|JKa8UO6NbJACvZ;fxQEKRgWmRAFpv)<_;Z$Yc(ppH$?6~6 z;DJEvZ%7`e{->sODN(eo#7(smdxxBG#7i#keWXY+mW5mzDca#yFNz(zwt4U=XOnSL z-Ll#SXNK>IhDkYA_%`_g@1v9+K-(hvI|@mvb%0^`Ja=H0GSa!8wkj~FmQ2Wt)kG+H z=)LBU>4qxlsmBzy)fDNTROw@6Fryetyd>N9Wt0wKBaF#TtB=vG2AIfF1fnC_LX*fI9)_cwlCJb9^@#U7!8sE>( zc$Otf5+*l;h?u0#~By78WTr$Q~Sn#rO-wr!-ri@77LR( znIGv4UCKBg9}A|AZK=ub-|1olgL~E6&wK}JbT6cEPX`-ll39epe=3L5wed~mvM51( zTa%G6J5swT$xP-DN72}k;!g~wkCi^T6E%>s&7OPwYUk5uV2m$doH**_60mb?NN|C$ zK%$H5W*;EZpUpxZ<1SlX5Qw!}n`P^<8^Q#Mu`x}Q*S-V82!+V4&_W^u(&VoPmrw_Q z`aFg3Z38ZdEb4KS12_Hg$68h0Ndd`omb zDF((lLOfYERJ_38sIzA4yA#_LP7XcT2QGE$u>#2Sq#ITK9G*SuD{Pm2?ufHP`5NA# z`km-6F_Y5cy1Gviu>@J)T}@BaPRF6S&lZ((7ZqEU-4l3!f_9Na`B>4-l6;S8l)K7H zG>Z;vq5v!AXu&{y57@uu?=%2eg^-A)Kh7Ia(;`}<7Q~9Xjt1_e>lVTub)&n^M4W47 zX+OL-7#vlcm|BUL@i=JP^Q<82{xX7vjp5_gdGh9I;*V{73bGw2)b1F~iC0Yi6P)|0 zNOgXhag&sy>P<4X6C4!U4Ub7}xE;B7SZd$Mt1_9iDR0gUDRW2EommgE=;pEpr>F@n zr9GAvvz?o0Pp1C(dUG5$Xk>l-UtW~_c)W)N## zN=9)7HU(^GGy@bdvt)hiKDf0^X(Xyd_`U>&8dd7Fw)P5zosCwYIzWNaAt=cEC&Nv> zDi5&U(j7!>g4T#w;J#x&gK0xv#_mi!!x6%$@sV#mE4a~NM>0>68vcX^Sx(*mkZUXC z@}{FC{!*OOaVR#jOqE5u4kTB{RRiH#_0jnRtb!tN49g z!mtOJq-Kk}9T1-x)KTbCCgorGq}Pin!zQjAv<+A%!caT<9QKKj>tH!SEdw9XG!i|C zOf>se%9In;7S@yfmT*#d`DqG58=LXOC^>Aag{cZOV2!dJR&XxIGgc^ndz9WwdLg8x zj>e^G&R92zQA(mo-n~!uxosx#Mdz>}AQg=*iUao{<5{uU_5%AhAG$jr357*R)2VcS zo5CJg5lx`d7Rxq{dv#09@_$;AHVZguiHfGKp#Ddy1?l%v=f9Ls&ffn`<Yc z#W%I=ACl&-7+pf~B3o>vt7GGqUpHDW?RhD1%mWv{yysc=86fIHjG)VwlZ#6Mke;iG zwyCNdH0W@Sg(s+<7tlG~QfZcs=wKGPwN}DrhLBGKPkfwnDxDsmoV-(xDb zjGcE0MMaNmHR3+q>-LXZjgqO0?cgX{vUKD7P~H;PoO8Edv|Xm8sZ_;Oe0LwwVWvWv z#M^6#TGD_LXb``lq{qKGC(7z{(kOtSWuI*qZ~iB-1nj0HsZC5~B>W&^El( zvkMgKc&tb%?I>N9*^&|mOhcJP>Scg_;CD>#t{lwVN@Sz>pGgpjHH9(qRAxyD7>r!ZUW%I`tK>8 z?V?mpsBTA6(M;mO^q2kJ7JlNUn-%r~t?+7hw?zLwQAj+|>@#Pu_6obAMOyIvaIR;u;*?96I|K%7GCpryQELK#BnEc6X zlW~(NYcG^`ly*}8jg!K%!fZv7tD+qKnwRd{gHDcNNVFgF3Kfsi>W`Ez?R?fo?Hw~W z>f2LHj)G`;ZqtR+mj?~)N3Ycv)Q9a{s3z3#J8B2?O0>PWwf*V+Q(ztH_OX%&=@)Tu z%oH1}Y=kW;elt0Ijld4p)HKvphIJz^Z5E&{A}{W(ZRHsVxVm6>2iV)fOYPu(G&Xea z)YJLZKn3-hWf@#vKq!HYsQv5`7z`)a-|FEoymj~13y@-u)qShu$wj&2q^r6yh!zjli2wMyDoQ&vVdwl8$;L@g zuEP0g#i?!KIyoSVCA$}q^|w&f!#J%GV_4x2^dW4a%+v_GXL0n=*}56S5PPkvxPMxL zFx$e#8GBY+`7wy!gI_zd!?1(ciNXYIBm$)#rmRP(Sq10-Znuo)Y<^mfplM}eE*;{9 zNzG2$FR$}r3qf6r;gZ%{>EkyHdI`77ITPr9<14J+K0MWbWqx8WDVousWzf+C6{xA5 zzFCOTtHBWc=E}+iHyQl0&UO$9tCB4n$K|pHEI2$eL@Gt>#fF*C!Rc@(s$q6sWLO2o zW?++I+~|_1g~-U;l_$xKDsaS&4R!>qB?Wf&0=_Thmm&qA)MNJk~4h%wH%|S zwe?gc@T3QmYfD1^U)k7HaJ?*QBMPRzz0%#*<0)t$y_AO(HLKhL9$A;NDnKxrond{Qyh zvm(h%#YB5cnObKc0w?fz?(OfvIFzhvkKdeLWQRGk$z2fWX*Gl?N>p%;yJ%~Z28?g| zCGp1`*J(abXT+**k}-W=`(dK#4+>6Mb{jC61ygb`XeN{*1P7awLD{%CNZ5wa>uTPQ z!GILiLr&x9H914LUPkQs6&r$gr2Q2FZtWwN$hp(Gz>Yidl60tjh;oirH$I15{zmUj zLwG!%+=Sokj$)f|b8r4*@I&SdY?u1ie%B32EaHIBfj<-YACJmhR%Ay>NFQ~(dp<>u z+2O40PgY`!m`rp0j6d6ji+_Z>PGgtSjw$h-U)kGI07x?rl6!#AS){G)4r6c|f%(*( zK?p-Imq@#4Zi=ZnWO60X^Prth5lLUUX;win2iEkYg^eR2QpJIpry4jM^7Qp*J|>sj zZ91K@a4O7TkY%)sS(saM2A{~}mz=|-$5MtoRbdBE08PD;dKUj+Oy!TeY65^&_^@bH z&RNm$MGG>S4dl5O?IIR;&RRI?k0Cl5XO_q_v#vX{vtk(FJjH|GVRNu*~JhqY9gJ_qA(#BBc3sO+70VVdg7v$9jp@ z9Iobq6Y~2M*8twf7F|3XA)K>_oW(z8{dV>0U2f-&M#9rXzB5E(j@?eIer^M$IK6K9 z%iE&zlT#{_&{APzo_PeyyZVZcvswh?H}=fkA_ zDEITG+2A>xp8{WdlvaFUQVLw24wJU_Dk&p8`_*w(O^2;LClw`)&%E*Pw^=sw~dVL?`r^9^o#>F1Phw? zU$5PRv4E$bjJZpPA=p{1)Bw~3P|8nwV%hEs_cm+1ZY|8+;ORl&wt9+qKT?isyKpDOBdyOQ-4k<7VRqkO7&M;pE040tA~4?S;^vJ^ z)Iik~&UI2=l^K(-QcL!-J0$lvGry^YC`f1yxi1T&y@N;uKoOC-*>7o z^6_ERwStA8S*M2Bie|$rmc`#RNXS~mL#!|<_VRkJyC%>Liqe$T$B9 zbJ@f^oOtg=eM~x>Ku5Ax9F?o z;=MEWgphUdOD`1@ zG!GAb_1sf##ASnFhxJ;6++?XInrq{2_*<#+-gdkeAFqUYNk{axPvu#}LUwNfT3blt zfIWKv9R|b_wfCrhIEWi1JiK%gMy&$jQRYNDn{cPWu3|Zf_JT0e6DElkLgdGOvKWKL zi*|dzJX{3M%g;x>tnC$|;ai4Yd${3FJB)^4}nm6!YJZEPs$tnieqX5DmmHSFMR#^Co3~F&!N`T**6FZrg_v; zWzMk%Te_PR6kH$s%;;bQNU=}n5Qp);jYcMQmokUS1=o1-b4qboiOu%?#K-}HUrD8M zEXA4~g@2CZ=5j<{*7cqQvodXM(f~Pb!AmY>yT+e5?rcSpYsY?1g9Z*li6kT8$G^V%Od_i<`SDv1XJyz#r2)gzl#;oc|irAx_TA* z0|qhRtgFzIXYF~24i>ikp8e3;ljY?2=IYmRN4;Xsbud;iNPau^` z--{rqu=4%&mP;`+(ZOF;o40&}#m_O`E3_VOZ!Mk7w;DCrMT!*Id?D%SR$)%(a($Jp zk)})>PZW7O_8FBga$Y->n*Z`DlhW^1Ru*~`wVZTiOlP)3cwSL@6kj+J%pTl6Y^@6V z4fM`NLamjA0zQ*4w8RQ4;!ckO`;Y18XtP<6@%4!r*l@jNspUS7VF4~F%Wl%b$nO4v z-8GZ6hMp5d>lL2DCG$grPbDVtlp&YwEpq`w{H=E`Z~adk&Ng_ub>X65%)_^3n7mBZ z%@VZc-x#n*80=g>daCf;z|uBVopd9^18dv*1GemzbBa?I*6UaKqFLIiT|CjQbx&fs ztM)BkV-&m&Uv-VCJL}f*YfOz}$ydYUMM49yT7N&eSxo>bRt!@q`{a}0Jal0UKhQww zQY#?nDV65ta3r?SxNves36pE1|@KLFyfk@My!<(W3DuBcbvW6|7xqrk6{5gPk9(2+O(%Z({5V3>NF^ zTp6gP1{L3Agp5$)j$G9L+CGug)+1-AOA;VT&^m*_|=zV3FpkK!)O<}REIbTEpgq{brW^LS-Aqh zy<7rIPr{w1!bn;uVuX0ZUQ#1;xCbYd;qh zEaG{-+t-V3e zf#vd47he7JWYz zQ8vTspK;QU)3J7TreZy5EX*4cxmzOC?W1KW$7pk+VZl>Z-;^j!th?@C;Y)~>pSo-* z|5F)ezOO(8*(lHCa-BCz1f>fZnExP!v80)u&+XAu;mgeZEB-cQ{|haom<_O9Ul`#N zXW&aCM0ksov7CPG$B3^6dUO#nkY)~jev^)vdAmo(z?G$CP058S~E;?-h?(Z9A`%Ut#R5@rLBlM>HD~+ zp)64x=R!fZ2m#wd$Ab`U1zT;3WyF^fY$A$tO?qO@=r*q>%eHoPJ~N{7<&iH% zy&yXe7Ut|BXT+e(?bWXO^Qc2<6NW;Ab%_Q!l+75G9x6f|&lhE4j?tYhC1Z(Zz4-aT z5>q>3HMlCPayTCOT}{`{)Pxcm%Ddt{h_z3fXIneqlvyA9T1{u2eZbPRm85p2BX9Gi z^F@%Q1WY>10ah<{wny`_ohb0zzJ$D-$V)qP(ZJKW?eAsnVez1i%%nc(}WkyO(vVjtK#_2N@Mzuk~eFl);0STnm3(zPO^I#3;tTcGgKtRR)tHT6D zFo)2Z8u_Tw@Jt1XzMvy)@_HpJYFs3)Mt3oTC39)Z{ZbkI0k6Inup5w~qvdSAu(vxH zO0^g%g5gQKL=6o6^zWS?;bWe%z4&J_k3Ilp_1{XwJrx5VbQbZGYobwx83sf6qnA=RYM2M3M8Fv%L5wU;ow_NQR96upAG&|{CtDA1d>#cG>q;Z$;v~)0Mg7TbSu;1)p&E6~ zBvCpiH{@rHS+rK%FW*`fsY$uIcVY>SnD60miCH|!j*W)WtmUzU7{dBa>K#&Wro+Ln z5yx0Fmg5G_ITZ>r4)oOZP4Ur`AZD;UcgEw-WI6O(opLEQlzLaY(HEQt;xsUjJGfa!ce$?&&JwfF4`D=SBBw3=M^kV_ zbppyzB}dZ~z9bbM|y z{=2rQ*#T9$8Ckb4m)#NPUEJG6e1mb2#E@XK`Mh8aVmA?=$OS^z{P6i&@VPuFK)5?W zf6jaLTt@nx$tC$?N^{^L>W3zT26+w1tC6gr*fpfFe2BndVo2^`q%NcVv*#Bf-MH%< zdgoLEF_&VQ$u5TKjV#6iefclOvPcPKIQ9KJXPW#WKU(QyL}`(`3a1HMxn|xosJ+yn ziE|=WdH-VyizF95_@P*`Ma+~df=tXt>P@+$&v#-7Up%13Fu6To*2d%bJ#1YHkzo9b zL0N(NYHLi>2l?KZi@>zYh6k%X?6}&-U#^)W!6zDSuXhHL9%PoI{wn{H<7aU|eS`^B zSH(mDo3AfQRr{*RlbfiW+GLPsEtj3@%?Ys@k3C&5jMcd@K>g{UBTuDj#Wu9iKa#al z{^@pL{W9j8VsP2lXDKP&I99>Md7t#3@l%27be?t1`c8&M5R0eWM6^~-u>A52K%p;y z7!j=F+6PZM#fEsd#;E0u(r(L~NLD664dP+{ZsxW0TQ))>ljwo$i?UCU@#!gB9P|4I zQ}cGiQ7^E4T1`uhvgj!KEKs9l0Q^iMPzOt>OM+M;^_un|BKOep{BE=ZJg*Q7oN+I# z%78#Kc5Q>S`Q&hmXqn$Zjt$JY9E(4f-R$g2!eim?*rRqQHi25nzkuPiGNw$*aJC$o z^Yl%5lc9b&eS!b*lRJ2JN*q_s$@{Bwl=GxkNS0}+myA@jyqXg}Qv6!84l+Ssd_086 zPQLPN^(^@NDY7F+W+Es)&-6rpcZcVpW+)FA^`GGu8ZJl04N`bsbWb1XXSpByO~Pk+ z1Cg^zE#cT73{Rg65Xivytxaxvd3{v7;jaTi2#Pw*TBYsKdTZ36*g7uDkA1JLAIL)EyTMz!04e#?@6dvn`xaawa>PBk}H@7>2 z^u`JAs)+RaNO-8mZ%91#~X=@edw$Cm!N+O;E>#ZJz>p?f?BNBGAG zYt>XB$15^S(jGO6E)$%8P92pa9;>}SPnna_HdHjNYba-SzM9xC!ecN*Ws5(BPWT#r zZEdb%$SPirA<->Vv}F{w>}40ABF#Zs8agFuK-YC|a$#@x1zigMNqJOfvZLnFXD`8D zMpg=KEp?-u_0vZlyUwI*bq~{k`8wl4bd7JR7{J_01*MWgG4pQ=eI0_7+psing_pj;~pN(^KPvxLa2#pcdZ*5v$cC50da* zN-xIOhh*CD=Uhv?9td{!YtkKJSPqh^^r%wr%Tw)4(w#KU^}I~NYnG$J1|x|tq%WSn zMclmT-}ys+4Rb#x7Fh%)Z;s{Jm6J$wM|QmgT4VuNM|#!vs)i$h5fc1gCB>i0Pag_u z`u3P-gjCG~Ab^veW;*caJgZ5t|UbOn-RMSqi)78h-EU+FD0|FzbJI?y}?YC-Es&>r~sj z!wRTu@A6-v2#cBu(R_O$zi%4x}Mg`ElXVz+hTbb2Yvs1HwjTi~bv5Up7F&Mki{@;~l z*^Ow67mm?Ns&lv7PB(^ErW}lQ<1_b@?8jy{N9h@hsKt56V}H@5p7V{UL#|eaT0>E+ z8mGdgbZK>2;63QW_HW{XCD1% zc5T!wJujc_3+!~DE1Ho%sF2qB8`YCUwo|ON2jg|}^Hwj#C?Eu>{hW*c?D+-g%QM8S z6Ns}|cJu6XxfZXnrhSJ+!*j<6Am3Kk|M?hwDd*7k(|WJB$DIsHMLh9wL0t@232>tj z0Ntbt1_A9-U#kW*>gUh*9vWCcHBPA@Ns^h6L=UsF3R!Hc7`9Z<)e z=}X!#^*$D(dGW;Y=J5#T>uvILVpsj9TDs%PG(Lba9;6dt1YwE#hkb@V(%cf1P9L zXnCl2HhCKX2ES^opCJCRp?~;vJA>ECzkQQ>@Z>l;jRYOT5W%o9-oC3@aN3reehwsy z?YWv`6x-=m2o5mtE?6&9Z)spz8Y}8(4}Ft7cE+^|Yb$Smr}_e)QxFL{b`TC7 zMT^JA=|Gc@qdH2&Tf7FvBJE&|bEgU3V$BMC{L{E{HZZk^N@IfL@AcvPTH1!>c_>TN zY^w3JYF~VJoM1Y1J-D8$t55_j%6R-A0q{DQFzL_2h1F)GhWuJgPH9TZ?%0LWP0To3 zK(J&kQk-yGt3vWe+6VJ&bCQCq@#WLPnznu#A}+{BOA8?Uo$7hdPaF7I#q(T(sJVCQ zXR?LMb`Y0F$0?bYd#`P;yg}LhS%s|cCc|8d&)q089-IoyD7c#Wf{me1 zO<%Oqxd-R6g;XG^?Q|@Q`+9aGd8w+?r3qF++@OY(bzV$Ooi088r4=8vzF&4NNWh*G zcMHys;S)d8PD$(o82pM6c!5>YkXXPLd;1<V^f| ztYC3?{hVtRr_lF{xo({I==6fMn4S`j?i(%|JspB)8#z1^eug`M&8SDYvjYB}FDYLg zMsQ@*lk@&0oee9Z~e0Uj|B$m zv?sly(hKKiUbDFigEj+ts7W>m<}ZjnB(^f03tI&*vS=2IZm~}Nk)M=!-4EuBRna9E z{5BHJn&4IEnaGu?-#qoJW7LJV`IumwqvD~rK~lDOi@uX<$57L0ON^U!#KrSu-=IyT z-xjb;=FyNBWSGe{zd~ZK-3bq%^LANXAvTU00X2`Yr~vuEV0!9E+e#U@pVdt0 zqwG2nr^IiCAHDHl+bAT?a2XO?)*<41jk3O0;M+Ee&;|tCKRR=u4@oU(gW2T`Wn<*& z^5RncekqhAev#ccQuAk24U}+~DO%)76nF6hSVTAPt=&uvn@2lzDZ3B?A0UhdIdZ85 ze%e%`1XvhJORLW*G|RlUEndb&u&8VNakp%lmQ*Scak<$C*#PQ%$Ycany~V0$vgXgV zp7`G`G~{tC|3@0)?bjuB2=$I9{$Bx#@5Jl>iv-U8f0Mvj8M*$edNUKUGchy$xBGuO z;2eyc|62hk{Qu#AyF#+&+^qiF4LT>=B>aol^MAa~CWKTZiU?iwMe4BYGV2ZSuf5Z0US*lh}f1y_6d+jJk&t9L&*El985csR6FjFgNNNRa_(nc$lEnI=#w zfZ~}0>f5aMf(Y{_@*a$Hrv(^3j7!J?5HY_?_#b^I^f?fs_!ss4ED{O|OwiEqPe1Y{ zTrlbn&u1*w<@4HyIY2T)pkMo4+)pszzJ9N7*Vf?=gi!wuC_wI)&zMMs#*WU6uUCNX z7gkkC=pA$(gxDAe7}*Ni2uM*@7mkpM3IqyZK?sU;tw8j9Py_BE3{>)4^>$hWsM~mp zhid<|g+koVZ+u4LSPTVe|E+b!F>*4xp8eRu&lE5d>Ht3p_}*MR!y5HFkl z??NO{OqlD(bYnnMHei#@0~6Oe;6Ci!Ufmo zS1y4NOQG;{VEw`Qq88|h&eUhotxulsbyMOZN-)i`MAhe?P&&`x5gQ8&M=>PcFI|7@ zRmjiNQx}- z?O;wT z;O5NlRY>ml_eOdBRy$&8!?&@*t};`l#UU%~)@VWMd?FpgmWyb`z)7X;JMUfHA1foK zrwD@G!Jubjb+DiCwcX*sJ$xRXBGcNsSpiqFHZA>6o8+1UpKaW|t3l(@?hm;P-vLG= zt%YBu1CxTdrqB_V`3D$M=v(Ayz2p8-)1ry3<0%V0-Q~khOEY~`Q2@pPyf1Gw0>Adi z@#MwP`)mXr!Z(TM=%$E!94t3JUT~*2a#nkK+{r7C>%UA#7QjSveo9PYZry6~@>`w) zdzRQ$Stqx2YTSAFY)C;|jp}NiW4>*FnoMpUjB5yo4i1ua+P!Od?)cFnz*wVO5Qwb#WkvFJ}$f`;-22^JhR;rkF=!`>-5bbC57RppHbjh|Ba{BL6iFjTjW| zw#H*O4N%DeSk^GS?Q(cPnG21SYQ$zx=}MD#^ozzY`ALsbhyxIy`QfvCBtUqSEj!y! z#@J7FZOJp;nZ}ct>UUd4?fmDbl%ojxrP^!rqqsqT^mD=_wH@GAHh;tBR*byS1}UB` z)nrv-`-8IWJBQth^+!!h;*V*ro~OPXvQ3b7AF}iie10QscUrV!k=|bX#7UWHV)jss z%%l<2Ve|oX9`iXtTq(+d&nl47?ENE*W64-utHxBQ6$#%RhvzPionrLO z$Nd75?&ON`@oe@Ud_J&CMjJ{A(m0a}7o2GAIK+&cKbAZmK z#*A*T+hvQpl=meMTs_l<3%in2UaO0B^kO+RQ++R*+^e*z_rQk@$+U&P$tja+;K93U zJsbCl$tSd!8`hLUJ?t zw6+d?nvn(2ugY|>?1hrA&P`Yy{ zyz!r-m`$}%mw=K5$RNn#&6T9cDNZc{f-!whl4R&0ux;+P!oBqOVZczqAJ1`2iab=6xL6-|!y42O8W7xhjUv%^JPKm+f60 zS9xcV%R#W4w{P^?e`q<`lY24SLo}sKU`^*_lFfA!(mTpoUpaKBuf16=GW$er+9eJ5 z(LYhX6dY>gKGawD-*okCRkPVxx ztf2RS?zpO_cp-ZKvbE0En^o}=tRll6CXGYNUyk3rHw*Qe+Mb1~;XzmUu-=QeTUo=Rq0_Kw6c&UzL-q|LivM;i;L*vd)6B`E z9iM^Eu}WEccH%YnLmAKC%uVa;uZ+S;hn*?aO0QwR%z9#nY0-wFw098peI#&qJ?ZHC4WE_j$yaOG*|Kzdi2Ak6wl#P{Q~wAsSrOW|$*X19vZ3 zVYRMabJ@;DGyxQ|LTw>ljYBr(pYeQNI^;I$v3~kh67l_n^ZM1kUkFs9V*c3A;h5Ff zB*XN`nV3pn)A7FkW8kb4yV!CWyu^grmZy(e#SYlE4 z&+$x+yoPz`+Gv^st8aElW49!C>M&5gX-Guy4R$63U8CQefY)cpGJi=p>?0hS14pO* zf*iC_WMVC(6V8IeWp|6c1~V#Bhl=uwJd>hLO8W6IP4Kdr>aj+GrWsM~(0}>+t~oO< z!jnvO#Erza%mFzu%J_)gT;BW=f8sSV*2v%qJ4a@dzMEce#6?9^v(5Y};fBIdtcg#oQ>) zZ5CBFpy)}KnMs^sP!T~_JK*Ctz0(-nPMSYB5us1*AsMwxtCYAYx+cE5-0Lkz!P4Sz zmR3z?%)Xgfc@B@_aQ6e|AT=+g5l(Y}dd7O_xZ#6I;1RHda*`LFiVeXu%1_$ORxIb4 zZ|BjIyI#(TG59;aWmq>DIp(E6#SS8)V!(4xnM-Nw_TxrV=#;qQx49rnEfBqsqB1Ez zmHIbcPX5M?2(CGA<}K6ew4r8tlfz~Rv>4H5ktN0+GeLu?+^>Eq+_rDhyOVTmYsI#Wj%{l{d=7LoW6* zGUxo_e9O!0C1Zc178aZAFV65@-AXyf1~FH{F3=#Z3fAjes2RZ4 zR`?}S4Z>pK{1A%#!cSBDXsgVITAW#o9<1nL*R2w3(Aia|B~Jh63XA*SQ>z;+-7y_U z*P49CyOe3g)uPo#F}zGjMjMURCOHZA8k%#aURNto1g??9+KI**dMp6m+IJ6>Dq%3pXmgoZJ7a329brap13^|}fbAhk{7 zLqF%8Ip_5Bl$IUOh+c20g4zj}^NUP@Clzu9 zI@!wmYqXEXNorvv88uww$h%dWVJEU&?EOfL#cAc=k@}d& zpz;SFA2=^#1T^*GOZLd^AcZiHBQ;{K#3o$7?hyv*7W>*SGZS@__cb8WrEENRJ6UP6^=4lx3z3}o6}Q!TA^2XTweRI&;xe<&JH_BYZEat8+luysAaKz!%6#I# z1j1Xiu&qDLVF2&+#HiA{+rq{fG-t@{)>V4sO%;wZvvj9ZXT(&+MWqC?3Iq`B-WE|O z!Sbniq87*av!{xbo;v#_XXS!|;}@bgsx(p2>Ny2>V+&5!FI`)|rs&5#ZHRnT+>T{km&kc-84G)#KIb}f26Vz4xOj(u( z|9s{NC#fD+ft3?WF6fji?BN7=Om<-YZ|IdgrQIc3ofs~^jq+Gz6GtoS_8kNCn+3P( z05|RWC@`aYDSpkVFFYBVLPR9#pzZIqu{a3 zQ4HgvwlAU(lq6*<*(0@GfUemC%b%J2KfR)=DV+M~1yAT>+ydDFGE_0gL0xt|Xd_cI zleLyy6T5csX+}jJ?wV#Vrfztd+9jPwSc*C>%S#;Wx3lEFfsR%xo)2PLwx6C9L5FNb zOY!AL%ml|;%Eo^#c)qTSq<{G?W0$2CtIm}JJx+8g0Rvu~OhXCmkH$#)bq~o9ywh<8 zcq&XRXhf-ogpIR%r}dKs4AQ?;y+F#!_geG@J4&moqL-*x)@Z0;W1R__K_sF8F0SKA zv!sR4MhmvHUP^$EmER2b-SZMURJfmzDrT7F^To;o-7eFktkKpp-DC3#iJ!{Hu+|&t z^dlR%RX2e9fd+yeEG9Z7jvE{m6+Z$_hU@GR28+NFC^A;G7_G1xyEv|%4wOmk@WZW2 zW(D}E39ey8?Mi^_rKgo~h2GvpG~r@J>A_!jy@X1AZ59tmyje+4l$%`se)=aOlvZzP z#PCbL#QOnCjLEKs8=T?i#vGB%UL*>FE!wRr+{cBKv+or(I;}ckg(Pcz?;Fyk#>$pD z%RN|9QG2UvS*vkyhf8nT zay$Ofv-RQM2+* ze+QT=uCJ$bCi_l;;oV< zmc?H&VHv9fJ8(JDl333dj!E6bk6lOm_O8x}){;GqxDKQ0ej{13PqXH~S)%la7{e_G zaU5zg&y%zB(DNlHU2p*@>fp4{4}C-=AiJGy%Pr%$c*Z*xdiYk%btoV4;F5f(O1iYo z3?D%p_^hucVFKgb?(cjiDI^uEB$- zA**%(Jpqg8-q3*(nIA%Vv1l)Pk_Z_}QFg!4=WW=epfU@}m#fDU$@q|B&4Xl@_OTlV zCMbjrh_}^3y{5R33k92Lk%$>hHm|QwhMS=aU2gM?Y)LExJ&9h)pz~~Utl}OQ-}&%< ztL@Zi;f>ibd&LKS*SPJYHoE|&o_62L7YKxOxOMbh-8<1*nU-sLdKd@-#~zMv!tk1B z_Kjxe0Ui0Y;EN&6YsteRSodS77|bW@Y=0C@p>~z~klU9pSBrx`Ig=PpOFnMt5#Vo| zdUK38(zc@mUDgCiI!IN=1CHLia28R5NN=tt7W&mf)9aoC_*un-j0duVzOYjMNJm9H zHF_f6qH~?3Reh-FTTPzEQ6i5}+UWXD)BB10JhP*Gg*#C7bkZwzQoW<77*y3K9u~%n zwK|(oxhW}EikwKlXdk2T#c2d^=*Tr_@?H|84?&psyi7>A_}Y6Q%d^0%c>UEd{nf1A z)n=Wril}O)52}OmP8IX|fxO`b>viTm%p4auy}h+g9*CUICAlYdsBk(kF~Efy~@XY8gv3ZhfPPBBxSCoa-Grv|J6M(P0pg z$6grf#@kz)SdpvjA}y3L8qyt}s%c?1`^2hfaZLd%&oP{?T0G^T_b{Qf_=7GbV-Vht zR-bkL7>ih{C8wkwno&b|Qb_xhU5)fqq-2-x({7y63&3l#Byj`V=?kR!J)vztN`-fbdJ4wxn2U>VN`tZ8GWoSK1nb3E!2?_&Env zMwc5KA+X43==M1ah0kNhHE@?_ltj*r5#-ingb1>5KWVRvvhNvlAcVSbp{bqu@kCukjA3!R42BiPLs6h?A2@`+-oQi| z`tf)_9L7%X_B*y~ug~b}s;iI%j%U1_s%Xe0EgBKkz}$|B(+ev;1!&`2X-huK!mI^?&d|W&p?ktqU$1 z2(6m9G0z@IBvuqANl8PL9lY&zL(Wzv1%o1umN46)N-cdO*BVF!5!j?WmTFxp-mNReJg{vWD7@`K=kydnt5uFg&fBD+7J7{f|2 znaH~#y#Gf|Mm?<0j3`lYzW-91qrm28W)V{74~N`G{g|a0$&{@ z(%}nGuv^3OYOF4{PFt|KxhsjE@~PoR1A~D7R{%)ZCy;s(DC$A{_1%n*BIy;XZ@lXD z4^l1;28i5RJWhzPE%t4nONtMT0x2B^FB1l~dkxacF)4Ms0%btR_6bCJoQTdP?nRCk z2FJfO?gd9da`;4UBgEQ=h2CW%xCQC?LHJxCferw9hQUDs3E?4$JN&7R0S^5BuD^SU z6>#@j?}Z#i-2-ZK;o{^S|s z`sWeH8#$*21e%?HKZ!4&(XFQ7uHVotO;Fn|A?Y7pjvqLoTKBO%U;RTyycp~yfe%Ho ze#PgWMM^RnD1QSc6hf!ZqXT*+e?k~W6|dG8dpB?h$T%u|en|gZE|6XS{8x}l*$GB83DII4+ExGbzc>To#CBl`?-O=Z3mpULkKvE-^V)W%bu!5 zqlZI1EOa9OEqNvkg917xA`ur#zngSr*`~*G=6#FV*ASH(9It0poCL5kurY za!wvL%On7vJn|DfnaeG02|PpDmkGa(g!Fx4HB0AM8>P?EKKqv}KF&5J+4G0_+X_W0 zEj6Ho51uw7XXqHT;`i41RStAX;RW~_7ue``p8h)Col%VEq%!&;mY!HdKDdMt*|cdS zOLtX{%cM?qq3)t+%}!irusZi*^4n3EP{&^{!llrr$0Jz7vX@|CbsvXI)-bg({P`)m zk**sG$gLUm=(F?cKs!nDGDTpCp;!jyxeW_6oQ{oQN#4r+rYqo9Ut5{5z3c6n_nOlx znL$*;zwK{Xl6&0fmpGTGqoM}UzFvrqIj!aN3`hwD>d7*?N?bLhpbu}LC!HM$2&}>9 zd2K{mk)u$&(7}iEtv&K>kWE6I#ur?8hv)vpIKEI1LBhSVb3smTS>vAz+-?u^)YIuYxYoe;?9<># zT~!3U*OUkIrty;O!A^B^%GlpHI7$AmgswKq%hAEsOnR^Esrs#yD=^Q1`mVOBK^rYL zzxfY?Tj_zAwcpM52*DGGT$(l5L_8OzA69?9Ag*f4znFTtARbY zNR>FFdD1CYxA<06KymHVz_`+N$j^fYHglL;^98+bH1P~23kS98duRsR_0e$!Opk2< zjv@B>C1MUa=Cg64@CM*H!-K#qW=}o~@Q~2iK}OGl={=AWE5d5d-mh`zEr5AJL+nV^ zikcH5rx&>yaYlbpERyi|eHT;Xi!cm7EeL2oukVagz4~S zx%VP;3W~={|GAwD)oG#@plY=(T?!!Ee`xu{_CRFkY z7z^D@)2aqP@jI9B;kX~i*G?jdHJJ+Ju#bnoj^ne@G&b}VQ#f&cuR8M-{Nqbz=Q`d= z)OZx*ejCHXVAHop_Ak-Ql`Or3MtHxzZZ6FuVul_!-Nq(6&N-r6?!w3d9#i6EX@r7{ zeh=Q)>A*4Y?I>9xc}wwuw%0J`-kAvu^0>;qtI0L09$4m_C^^DqjXZu?%Kc0uN(eR?8`S8%oNnqy_qCngv z(ml9W`7|XJoffe~ux`Bmy;THS9Dh0Hd9ePAqp}bdk|GkAEMUyWKzCibWMQ z+;Ppw?oD z%K)QsjVv;3N`xyE5VtF^>rUKTr6H=ulGoZsV)2!@;B7swfY)+1 z9o3a>wWmHdrU)sOe0%fFd%Lre3NQwKsGeXTU^J2$E&d>?>uo(9-drq)j%Q>}a=?5E zV4_2=UkEKyj!=_c8$w->XKjyJ7~7r%OSX5BWtF1Cdl>l{Daz-Ot@=u7WUtNiYRC{!dQCZdv}+TzIf zusq1luxIg-@~R!96tWd2S&2n5+P;Q=qch<@1n7Y-A+xbTSbK2F7w1tkNF@uMX1q7dH&%aG0a20AmFBJN? z%pRsUxfsNeFBavY6C4~{eBCt$(u{p+iUsjwJBDKZd~MVrRJLVm=^o;4PT~v%y?SDa z%mwumG4ZKV3M$(DSqwR*{mQG%(+K!(K{!M||A4YPqyDc9 zeL>(YXyen$PM+$N-Tty=xH!I)OxVv$n0=GsdQ@s7SvgdD6ZLj;xNWBznL~&ZVaIK* zlRGb!(2PjZ!XdjfZ(?sDbya_Q@-C`yVs_Z=c2KK0k-*P9YZw?InZWqe_UGzEyZK*Y zf1MyB$k4=8$>gy2;Qham`@BMo4)llvz6btl9pn@j20+h1*Q7b=I!k!T+qY|UW#V32 zW`uHNjqAP|ooVtgw`cMf(s(#tw(O3y|7&!>=3^STw3R=F7LK+T0IRb}vWg&BI)6 zx6)kD^ZOQBy`E(K!H*sc&d^x`0K8IMiEfYnLJHW?llpX>i$%xTPlpo18V*6kAF>P^ zaVmCMzf^9?SR7li3F*NH98+w;MA8GNzoNDrwK%>NTN1Z7mM090_9#ro+_OwgZ_2Qux+kWy@ z1iG%#9jW9>Gff!tvbILw!XpE#jIvGAs{g_Hy1!W%9CRWYHFmyW7c_rvG_M*s**hE| zM8}gc`X=UJe$%dj&5(ty#2#8oHSfuuv+?STKWnDE$=VvN+h8%~h=fc<)Z!Vvpt#*y z4{(^wAeOAIZv}Jb5~(Jp{c<$(HCAb0ssQqSu@p~PweLLzabV{vMDJ4Z%nS}; zMSamXey#uyFq7Ai2$MqX6I2DVE!au<1tp?x+T$>SkzQzz879x4Z~7kGoj#^p%!YQK0mnAE1M3?VWptun#OA57>ung-etXQb;=K z$7B9n*v_}f$XY+dwU;k&o@9JZ%{mC~MTfvqbR=bwLyo0woKew2W|P5iPOZa8+$GPq z6x*zRt)kUy_AX+{968omj3y7P6eIjKvmPm;MHYDS@kYbi zHB0H5{wdG{DmkuVTZe`~4W10k=eQbXVp!NF)>Z`uTeitqDg8D~dJiwJ$t`K~jZJs1 zS8fT?X>+(?7jDRr5t#M{}?T34g!e)#NfAT74sdCSTmlIQ2UF)>q& zH?4G3Q~ifPyJWXDLRLkrin^^%db6&i0bWp~$xq2d#lP8sNJ*cyy*m;>;Dcs@aVg_J zE?P)7mDQF7SN%GbSQ$%JQ2+PAKA6vUL#FYu>5p)B&^xhb>(i8(eZe(ySzQwP;hMEu zn&T}+;^~!ncpH6h(&>x}++K?c9$2c{=^_g^no?}#IwS2~#nE*{V(M(fKCO@{Hra%H?WN}P zbfeIGX5SJ+W>w$}*splKgG6o|lqyq$sE14dxIm?9K;~5yfG#g)L`p=(n41bE6avR4Mm&A{X}q9l53irUboe!PRqEUZ#sV65zWG}< z4?)#L&9^T`iaa5;4^-jFt{5`mx1Qm*YuK_?O+SQX3L>h|tC&F?jLdtRQf@Mk<$FwF zr&v@syVUVLxDLDvMiMG%qCZ7QCJ!b{)&{8tO#`S1e#DgpZsX40HsWH!n10lqgYLW+ znhC01ciQ3+IRs{!C3HPGb2{?~)0t?%lB3RvhehGEEwd(G@%_}N|H8maSXzSDX|3I< zLk%nytl0yiQo}Hw9$INgYVGmMxmwNS=LTpNaDl;`;R|nba~SK|0cyy3GG6^HrhKYE z#FfMa$u@R1lS1q4igsvRg=0#9ModU;BON5-LrH{d4C=X*{ZrIOzG0W}fq*u@kp%>a zS?Ztwl|a3c`nZ0SW<`E-_RwJLYSHzqSx0-%E?Rxg)#z*x>ybqbIk>fE_jqbdzED}H z0F}-DJg|CUyRexRhOOHmaVe)OPJ>Tx$!vUlzw&6y=TBIj5jKf&%O?))DPI7V9p&n$ zjp2RaCSOu(PUtKHA-KC|spK8Wm}q`W4K(IV(6F8?Q|wPMEZL$02yCH`euXWx@osf8o`G^{sCWleV;3ls5gbf=Vo0C6S%tr}~#`rT4d zSaDUH5%tT2OV^h?|9XS#iPby{7PiODfFG3G+YPUM7A0OucO5j=pi5SAt(qORn+S*2 z2i%2xs_q1rj?asKdYs^8Dj6oYDIcnB#ox7THIJnZfvC8GP--AMmhScR@fA4@&Pnd8 zwUUGWrR5)8Gm<>ST^WKWws?*BJtP9zt$K!WZrb;pdncZ{&xg(S_YeMa4XHth=yZc; zfZSG)ST0VWEG##wkS^9w`EX2!N4=X-RM%keAZ)8O$6sW)EyeS_~q2X5kZ%r0=<&${4O3v@D~BO zogj5zv#Ecr5ewl-IAu`9h@x*jcpDZg=0`a|zaHVDlt5CMr5RGXq^D|XyEEdL^!Den zLANG}XOi&qD2^&8RedLNVP=D{?{Lha{8&Wt+{zZ)%Zd|jzHs_+uS|r^@K{ow=8|?swUD z4Ew+F{I?E9I%Dse``-c_fw!TVXM?;NZSrR43F*d@3o<8ZhMIj=<4E#xxSWDZ?-CWa ziEqlrNM`a|6uGU!@JIP3ET?zMDpL52>G>;iSHkxL>>v-3!8Y~6YIUYq8mzK~2;rk| z@z?(J1prPO(p$ZDbMlF@&OGF{2jSO`&HSOme;OEP-=Atszsn~KyzA# z!8S0axya*?4vs|Xmp*LY;7_VVROQ58ug2Mu9PQ!5fkq0X!XQLM^Qt!ZMU%vVN~=Kj5VneTHOUf9xjN+nsM zDM;{#W8%lqa(q?(jiwLW5y;8e5LT@Ur6n66eqe=rB3|KO!0OWhb3A8@3hHC;b|ZSj zxF|sIJDr31Fcq0w+M;4%-4}X!8!27z>MbQL#10ver0`vGQp*F#P-3~#?Dm-jqW}DU zk8wc&%U}=mB>J1azB(thpTyD-c5OE$2|~xkPg}lHmm2AVBUS1{ev#2%b_-o`Gwi91 zrT}+#&84Tum`(|Z&pps1;0Sbz%977B!V6khIm{obrgwdTbSg73b9%0;!6|%)Dp#gD&HOdk zrt#9cC(a^C`Lir3*K)Zota1AmN<2~|fp;n&TIL~4HpeB@VZj(S-`!ctwU_Dk zLJ?hr)TWTgqnWL24(^rO&-b?!wC-ucEEE24;EuAK1J4tX*)kaODi7$ZStUJ~4)v+0 zAnl>lotG3N&CpcIkwcm8lXDNj z67ESVy2y4}-T{uNQ(lSYyfyJ5&V;KTLO95bYHem-SVzKSe8)|XJpFu0-6 zH2jnZ=A<)$Sry4?qy_ zuylf!0*srZHD~fB3;bwh zXDe*3N|X1pJu+c+QR$z`=J3H7nfookX`4}*)+p`$t4e78XY|{9ZnehIg0#$J*|6 zQST8Y*bbe1k^yh4s_2=9B{08{ar?u_l{-R}e7mT#v<_z2%R5y)F;^6OoQUt9g0a(% zh-)#oUA6mjy?M0Xh{I>%wleRLowiX!X4QDxn3bktnJP?Be@R9kd1Qn~Bbye`5Uur< zz@uY8oDLRJp3x|inqS&Y2b3i5k>t+i8=G~TZb-v=O>Zo+7upB9P-#{jp#Dz9Yt&=I zuR}%r&zCCWhi{)uXRy&;CqSp4ICOu=16_1du*_OkE6bvCzCzo#puoR#Q!WI1LlmmN|BO-8#K`K0Nq|5@YylvJmHSi=#yUx%YV!!TJXl z7nf6kK%gCh2SC913xacw!#-nEfZ+CRZth725B|U@Is^$RTF9aX3XfQlpn?F9DM9|A zfMCUgL`wS$f`Fn(e(DEEa=`kN45HhC&mMx85~DH!!Ty2G4~-z&na2tfKGbM{KqsMt zU}0k)y=5cCIsa|z%R_+>L=UVF`e}7f_T(f4>pu z(DpuE#cS_thZ23#`yZO1Zd3B$Xt4RUUw;ML&_B-=P_Lna@@uw2e`d^nqe?Oc@8}O; z#nje)0fM7L;d_X|uFpYL)IRY7|H6C+t@SB_7{;ZQ{!0S^IfnE15V?BuOXj$B4ES;j z^y@o?{ru92V*t^9L}Bu4rAPRZ*gZkI3kFrRBVh9D`MG*65d(n&x7k(1=*2aImO%L| z;l};D_FXqO0Sxj2pA@vX4gme9-{aTAC%!NN74C5Vrto_6b$?M&R+^dN^ppO&!|&ly z$k7{%jj)S}0tMP@pu{AI@wNXmL&J^oy2Bv&Bc&#`N&;#yn`<|>@>mnY_k-qtP{YXy z`a@S7BHDljwfl-ZHYLlowE!5hJ0PF@JeTJaSrCsBU(T8wL&Dz$1b^>308Re` zuj0@#f$KkclZ2bMsf$^=Flk8hzrQaw0;jRRb*Fw?2A{yU{ct8zNK?bJ~_eK2m_}0;&UIXU9ngTq2%-H>Nkpo2YR6tJCzd)NM zIRf>qvCzS8@je56wF3h|d<$78k?!%{MuZfRIXGM>d)@-B_w}Fcn+6QFnW*|U+0bZy z?JMHVyJ^u&oD6mF(}Z-Sx>=F*@?_Wt#Gb2U)apnr&WvYiFskKTK+${M+`6W2@&ZwUbtJ{hUmR zVyi9__0)4tQpv~i_LfZnwpz*giCcDV$thNMpIiWfSUwq-U=3BUhG}D(giE^0VcS>2 z;>L56ho_M4q5 zznTgzI4xBkv(8`n!$|E*SQ^x=jF~=u|Gfq(LKM|MTGSSYVA3h_`{x$}B7@)4$v1Fn zO@%U{%Fn7QR{Px@l1JEJ5Grco7;C!{Ja!zSdA4v$&{3$cgl-P0ROYr3Ta9Q_YI6S7 zR!Vxr_*(!A8Kek&1oLlR-i-(8;XNC`T>D|jPIEY}-wxk@AU`V*7R&&I;q_Uw9rMaPVf+R3TJBtxHkoq+q@WZmLT6AYE zB_bNSqaN;(`f+IE=o{m{b&_jb)_%$8zx9SSa}$-!CR(o%S8d#3hNrlhv-nD!3xZHJ z9O+Mp!=rTBBy-<6u@ZN_R0p9b&rcN7u8+y?l!AyB8Y*x``}V#pnxrJbs2})c%(#XW zvZxF+*$+yRCo#PM&rf#B!~aC6@)cNuN1*B*jdD8e$a&Kxrdc>iE{~ks1RN^Bq=MiZ zyjke-v72}}OZv-slzF^2{3RC{G6U*^9AEJE&&+mi0R1f1nyrt1H#;6e_Yh!>?DTdY zHnKmPy7hRwkm&)4aA%%XuxsK2QxGpajVBYo9}{+1q(CNa0`(mDYfXT!Dk>&8u!LT* zY9z{Z-odwPb|=wXC5wG>=A%aAK{+4w*Um6v68b$2Bs#j@7P z7`=}%pgw4e0z&(OHQ|xhzY{l0bmV$;dkd4Ew8Y_cO%OZrv*h3yCV~hRdzj9w5gnt1 z4^BU2WEHtS6u{skT5}2EF?-{@R(=;1S6FKM%{fQM^EKUVEfu>(ssMRlP7)X;e=x7+ z*S?i@<%;>~2HnFa&%oxz#^$fDu?IBFD`HR&vqKh~RI+phsBH406t=gMk}D!u9sg5}KwAcjjqi z!t~9s3s0l9Rz3R&nzUpXL=!VS0ki1_c`yl6;Ujx_FYz$`luwQRXSMDJ1@kQRV68!;v%UuIwV6|7Xras&60lnXj$rVY68VYC4Gnq-V*z7|d>GyMyWisJ4c+HT0g&F(#y{ zszzr+k3mB{dp_eSs35i_C+uG(xF;=TcsB=U47XCcy^#Q9#xbW}8tSy~-ld>;=6xLK zt{Ux4GG;}(lPJo-IY@r9t;|=yEq$40-baY{^|zh>Stm{qEyktkQS82_xstdbV1Me4 zOVOZn-H4!)4qN#MX?{oNU$+A!tL{p+!{iRkC}MKTpDp=THzsZYn=iO(&ilK+R0YZCZM2|0>`t;|E&Eh zQg?xYWMnlkZ?x~wRbj*s5JXv4EhOc?jZd;@ZZ_Pf)TkRSKPyQcv|c-@(osk4Zb-lI z3s3fBG05*lADfZJ8X0Z9k=*hH%IavkH)2atwrn@ye{yxQ3O1IyI(XEV4To6hvIC3g z@Lioly;z#5*y}GXPf+Ku;feb>cvUS6CW+}RVT%7g6G9F^RhR#UtJDD)iRNqAl|&8D zH_9cEEar*7)lYqAq)ro++5WX(2XS5H_CoFK7x0W#j}}st;1T1|{?F4B|8$*M*K~o{ z!W92PNk=F$&IPpOx+k~gr@I(Fmh_;&lZ(i4bA9@CsZ;FPb;kP8>D5WOzmfDS0U?h( z?W-MPP>;S36<08nrVLa264o4pUSU0o?5o?mb~-4pKlFw&20uUQwiuO(g0$TJOy4+4tC2e+m? z8D&B!ba1OU!^065O5z3aK%hE4VeiJeYTuq@Ssb=lxz{W*v0=Eh#0#$e<|6*4tu_yA06}B#Nqao_`ap9NiWH+R}egxLQxFtE9L`I*GzmgUaW)#NM;=?+dJ%&L)le_JZ`Gk~m@mj$Bp2|=oawI5e z7;~@54-!~AMetF#czlq4nGmWa!*9T_n_0MULrz|6t|c=i8g2t5DF4-T-Dmi+VmrH% zx>KR@d0b~dotvhcz`J-nbxy_Qz~$*f3H{J2dszz_E#f9uL8w30@Jc+*F@1>1#K^A( zD^5<&33W4=L?7KGm}W&c@IINDs&Edy+aFnu84hf{AN~TaRy2#}8u=XqVY={mF~HIw zT!uBAmwh;QhZ_!te&rMpnw>u@#vSd7%wHM9o}zde(dlAh9iD0JjdMEv=C%7j`c0QV_O6!* z+XzKl6{*-dD!-n6bL41TxD>e3O@%OOAL3#r(ax*LW;ck?Utf;$ zpODgKx#rqM?<5vCn0m?z@DzyOTkRM0cZF&!5iC||tF$N8q4UEz+9nTQJ16PGc$I6c zkyr<*BA!YnNq{Z6@fZw94E~ombhkFF6t5Y^4^VAt@xrDCs_<0@#tSkz zp!U*D5olhD-vvx>R1HYQL<1bw9>ES4uIYEA_@KFJRAQyWaNwn_OaF6vAOo{XsW@|= z2#!w<%@1MU*V73bsqeeTC{$g+7$-Ab5Td5O*0e*yw|zRGPjfX!@1`Pxq^=wrz4K4|i-5q% ziJu7G-U@7UIIiqu;uoY}3c6}m?sE<`{gvnq<^bl8LQT00+jK5v{?}dhCjQ4G*+u55 zFXA5%!=Zn!=3)a$>>p$FPaw4vB?lc}n%?r+=-xz(J?LZTi|lsJ>ww0uKpG;;q21~t zpOpk3#ABhT14Ao!^9>^~5#!f?!XnY=yY9chMhu-zpqp-r0EO3^(Hj-dY|(<;5e4+}%FyH@{a>9dzvFtekN320?yfWhU+TqAh3el#Lt>8mK`O zh-wnl*q$UOcv_QB{GvYR{MYY_cUXQPiEunfLHvQrC~*fXb5xdGtS!7kpO@U~S@*&M z!P_Lc!$s{k&17O?ulL|l&s-RM_a1B`A9tOQH>~>WZM^)WE|NzG!_b7K_R8I13n&y{ zD+0mF`(ylOus(esDe*`5*te-_BchXGJxICkpeipf_j5s@bYTbKwEoxjls+ELzs@1Q z7E|eYN`bd0wHwZ;WerWzDd^XI_lLIO?R4FuW$)!ouyS|^1^}eIVF&jU>_kgUeZzs} zggk6LOnxX9Agi_FovhFZ+*Ty^kM#M`C6{%#)SP;Tf{NIJO1wNo?@Q{K1BNas+_yDgFtEFLs$&$<+_^oZ9a&PhCVvq#5@XrrxsI{$BP>pP^x% zH#WEnUC)yKoO3+s6NiEs*^{YEbNVMhLiN~S_)GX5hvHr8x}cKrPGP~KZfFT$%j>UO zakVY+MMtJ%-Mjsn%6aN@w8UtUO$qzny6p|H2wry481M{{_E9)g@rSelhm8M9L(&(v zH5sWk0!XHaMxk2m8!4qgf3uq;bT7(maja9mSkNuvlVcZa7W&plUYv%*kxR6+l&K!6 ziPmkKhb&N5$zsW_NA*g<|2tNvOAHU-q18>a8htgwEorv;{+Q4CyPJQAa*zeBc|~A- z^x1!z>-#xP=xPR%U~}s0#&L&MbCzgwq|{afbAn4DVqz>Tg>;? zm7??B|LJ=fYo{IEb?28DPeHcYTl(4DED>jn77Fc|Bw8}gt3{Fc@D6IC@5WrYeY_Qi z*5xN2C11#pN`QSaN!IvzfAxt^EY)K5da1kpx*gTKMVfv0GTBOX+@j3NHD*Q6KbCp( zO2@VUGhK;Q9z$t*!0$Xz^>AUKyCb^6{r-=Q23WZzXBFlwl$Y{Gg8>b`S@&{Kc2kVa zf%*bY&6Y19z}0heeyA3Q?nWx*aoRTlj{g=HBy%P@QW<%OuwN!{W<b%E3 zwndwqS~oy1(T-LF-SBJl#urk~&tCZxK~GmL^>pCC*I~?&7w?XK@bHAMzy2M!kFlZi z*ul`Or2kxRRJPBVv!lig_;-@wJKZ>!Yg9M!5f;3k*3-^iSxWr4QD~BpA?QNILLmv( z&{mz;<2({k%YkVoLF>g(D)_lg&*x1^Rz*mgt2FK~Xpoo>^(hZRX@u5a?R%LXT9&kr zu`RlI8k&xPj$R2@JyLExP3*tu^9U0x*k4qx5xgp5;_|Iwv?|nbpL2DI(&mdZ=i^cX z=5`Wx`K34e@cQthnjCBtPG?>!c62l+TQ~#I48 z)E}J$iU4jvM?88m+L#S}e;%gSfxYn{lt7+#9;6d!$i%4r!)=v4lqfbb9yy7vf^eR9 z!kfK-w|z{z2r}0}Y_$_8k_{9mW>33j>(o)7Y}MQLw=7w-`Bi-!LxN7>jdMVCW6)I{ zfz%wbR#3|hfF1Y|aU6a=M)7sK6^`>!8-c&i#JO#kltiusfCFiH7WC=!W`k89(P zUrg=_1HIwTXZCLe21mncNZ#(w%Wr5ma1~F6bgs~8oET??Bjq3{m-~YObGeg7DrF|9 z6FKyPr+RECko4K(56dJnp?Y^NBs_X>4I{t79{GAD!u1|Tn%8P%$C9i?3>PqL>B4pUmDS}gtIM`++qP}n zwr$(CZQIpl8@Er~!M%fja55q%?~WYiVDGh`^=O|?H_5Ph%}Dj`6hNBd{sw+0xBAe> zS$|4Cmuu}SB^HNuU{0qlu6tE#`jxzBk~jcALg3*^n1IHw@fUOJQ--=o*D2JMLW@&s z%oNh%g`q2ERFo>JwcCUCMba4sPto*|D;dRe;x6FreSRoeA2s-PGxdBb4S65}d{q?m zc8LjIgSCW?)EGC^ySxH`Y@w}b8%%3GPJ-kb!{CpX4V(zE+Xlxvs~vmu$Ki$;anAd~ z>^cZrV%_OkBkHQH;b0!2m1T1=?P<8AIGlvHpLaCRDh4d$#5)^Ci+t_5-LuGh9nVZC zgxLe#id!ChCwtI`>iXxM0Q=?7o=jRF)UlCzmp``p@I^q&e12o+i=!Pv4+AE-TgQUF zD;{^`*~03U##N9PFE{c<;f-RFyW~kL7nl76(Hoz1 zs7qpAjh}4*GRng;2oGw(EXSjB-2v<1B z6|~7&Pz0ZeeH2daiwP^`;GMB({ZRX)^<(op9G0IOstsu~-&3nI=?IqR&%l+;&yD$J zXKRh*Z@h;*#?2n8QVlZkOAC(7v&e-PLp}5c<%Ps*2vlRen=q!e+H99ip2%1(%hF`$ z>&geebmbQW;Q)|Vsfx=LLZcU5YaeQtCFS%v0!E9?dF~r8IC}I3)IqIPs~WinSA9FE zv)&hz%*@OlKj@y~-j(tjI~q!!Y@6@#m`^b$d}ZnbFIkU_^I|3(l{LbFNc}JauZ56OcJGeohZ+y2Es?mu+G!WtTWwg?}$%4 zHDKkVV@R+;GNU3Xh7HANhwc*34NYzt6pyzMc~=Wb(#rDN6qHNpmneDdWWub3=^M>G zm%Y#BjgA&dMM~VcOCw{yCVGA^5L`%$54AS%GjMm_24&M_(MaU}71h4{UW)GT$1`#D zh9Es9@vsMuXm`nE-W_8}H2pAtrk+1Liza38Y{jl|B1Q7$+PU2fxe(VqeZ{EHVJ+G1 zntQE_jn6?K%^Kx%!=Ct@3OgBn%H;gS32>&LdHfekC?b|*^5!7maiFKJ(hMR$hf_xh zOYsbW%q%0L{kGn&)bR7|N5U_4Q4Fx5vhgypzZQ{xu-RK`6Us-%`)HV_*>&i2{;Fzh zaDH1V$*#x$vgz*L@-FCse2M8n@Gpfxq!~9qg<0@vD)w9(HHB6;mBqWpwV6;jZgQ9_(I`e+?n>SvqV?13 zEA)KvzLYF*g4QL}%YN%FJC@OvT+N`l(9mY6(Dqfm>RWNDk^QxhDXP$41CsR9WCPt& zG|il+`)X4riU0|#gYZG@S=4zy?u!!mx+ckQZ9A=yMHQLBVoVe=eTVT6r37 zUb=Ep3M>60Ti;LC2XawmFW!&%n{aOsR@JC%_VWWEv3`p94=~B_Utp5)|HU5vg&^4& z{u@lPvNHa!!KCXyF!{B}qOG8?K=r4%?|R`6KKb@SGgW8Z5DgXeHZkM^aUK(7M3En) zn3!v85|bGDek9M6@At0j%tsEl=}S(!0l)v?LSp%W^YWnt2el*Fj}>Nh&qHc*wq^!?dw|*r zx&!p=?d5lN{{qxW(J_I7{{UekuyUp?PKV2d5%R%<3HNPv{U`@$F3e0%s{=zmK0Y3W zvbEnA&b2IXbpd39GI99OPM`%n`>X)HA>d{LS$e+HBeFw4bFG4H--z>q{o7plX6Sn( z`(gjGY&ERx*J4=cn0Ag_|ISl!@{WHAegmoe{PWzaJ_z`031E`6eeOPH(b^fIwAO58vorJfK1iAh|iiSPG?z!F|_$5@q0UQ_) zazQ!-<{?d4Ree;04r}rKdndNG8VMFeaF35b1v}c2zRqQV{6qL{X?Cz)SAB5e$YhIV zcS*4QoBV$;heozBUeUwt>_f}RKA}475x))_gWG|E_zemI4uAn^Fn|-VW{__Vf!$4v zXHLlP4n1zhtY{=fqc zg6tl@ksSneU0bop&GKwO{XcEX@Gs$kCiAhazd1_|z}pdns1}|>mb~=ved6Zue1a;G z24GAIqwNAWzW6hFKM_{3p4z^2#%Qrz3S)iIl!OH2S$czr{!jt?<;^457Gx6ds0{${ zZle#Z^mKa7p#k@V@;0pGewiy@0K4BAXR9MZ|17`Jx;48gZRh~AuU4?0b)dPwwu*8H z03r5%EBs#WpL5{z7?Itkw@e6D>Tlzp6_#HoV1Y=|vUz#%0#}s+nKQd%3~u+;u#k}P zPPoE&$mlf`3+`5Dpp%k5KP0adCKAtwYaUZC#1Yp12f5C0B20fEB87J!(V;T_o|9;` zJ)!CB;zO6?ZPuJ$17BEBQx`q_V6SeRn@Y2TDQi| zDt8XbMlmB`1R%*OA#%q^Z7)xk6y3Aur z7EsE}SX^EPJShlcnvM}q$ z6NQgGFL9@!48eZ79}rzGMv_frav!WDPY_x z$LnQf5-<6=*R>A%gCH48pq=a2`}}-0?aM(d$gdkI8>@2vrn{N%jD2Q@HWnzJqjNp} z3ihGw#+IRzPrZi@&nlm)bh)d&sttnb9uRr5G$0(W8A~QA>;>f>H#LTuZGcoXf6(_x z*mAU^GwJ6swm*Er_^q5aEkMnebi$Ut)yd?~@j3ARm076?o zOZ>5^vK6@Ax#mh|YPS6QUpR{)%8~6guqoo61b5#|igLwf`fSNyfh0WcN(V-(a?HITT;yeyT83d@4M19f-L!b zwuHQ8zvf!IjwAh3KAK9Bu`?DLB3&*#4}(&y_BK|FnA^!mqD>l}GV&99lDoUs>pPOM z@yDxJpQ6CJ@UTogW78z9z)x9Y&Q=uE4(Fo`fXw;>C@T~pk)OHeHoEe&U@o{4okF`H zYNezHkuHR?7U-U0;X~iB4sysJ?}RI~HllaRohJynIjaX^?6TJDli;ids@X(1S6gL- zPjp`!BqK&*M!NwUZS4=7$IE*kwD;|Y)+tV?I4rRcYl9pNJv*h~Hp?URyl3{ugj|qo zy?gCVR5vZWV_R`XhQ#BQ@cv3wE0O)yz@k?$>#2nq?4!7-XcVqH)ByFG*L&eU|4L&DhMsu5Izh}C&Ra^c*i!DF2Ssz}6Kqi%8l}Xg zi?K$?mH9FVCBU9z#g~GSh41Y@JK`5sz+JOf&g{T+iOzSlGkZ7<0$S+d37$rasqa}M zqN(9^G$!`K#WOH%gK>{95?6A>!9rVcyJ^xvT77A)-(Dp+sMg2ne+(i-?e$d*>!c-# zJEe(`5;2n5((C72fKL5w;kB-X2get++XW+o!97u6Zxe^+BwaeAw=rC&pwS;vCxtFQ zLAx`vLXh&Z4-=WfuEykIu3o>6lj`ZbJLi)@=l`xyA=|w&<0KVAGdw<~x1ZKC_%M%d z)|skk04#VSi_ebHKH+=ypHZS9OJjT~l)OjQEnIbD&a0ur-Z^2PII$WAsZcv@vB(0; z&JOQ?TL-1sl$Rj-OMG3$j6D&8uum^5copFW_pCj>t8Xd0965v0)I1`|3%(n{_Zd@* zfhFw3!FbEjubnx!T9*ksr@^S-~621!BK97u z{V~qr^pizONV?0AVW=ZzWyf5HT;%}Ixt+gpPGuf+v;OSuITc?CNn8)F+AW-@!M~$O zr;LYM1Gub!YUJ)%90qqxX=B_K;W_B7>!cEfX5wjUs+KS>l;APPX|a zwru9~TlUQ{jguJ$bymYgirmq*``BZR(vv=C&a3iOb5t)^VSh6MY>M*2Tt-iXRa4}- z*>%+Xii(ye>#T-#t6J_)M(goKvX=*f9?qbHdVNw)VktMH*3QHBBL0o1c-O3-8a4c= zkwiZBIDEw}lbT_o#v{0YW-R0F7~Ff3ndv?}MTj(c)Zkj-4ilV_whyX-jv|kc1<>Gx zZ0p)aC-4K^z{F^-!-$&HIC};E8PvVy9^$rccURL_M($WvmoszhG*i{_w8??mG`&!T!H10Ta838un2IB%Fom+F+*szgSXX$z0e|s^ z(*?C&_pD>BW$3PRRbyFBOwG(+6hG~VYZ_%~4gdFPnEcaj3?;wrg;_NDZuh$6a$pbn z#Ji!L(H9Y?E|iZ$_@LopSJqHaPm$%0HxG0o-&?CBEgBO^ukCj&pS5sR4uBZ)&Q!m* zIOB3n`@|!;&v7aB?3N{Ibn7n)KD8&)H*H_>Yp>6(kjgb25ro691L2E+CE?|L8n&ay|L1acg#Bw@zU`m-tZv~t_T8kU<5DDClvcEo5A?N(sit`@ zDmBGu<#1I8>nSRNw9pcFWd+%kO7KPl+HjTqgY-IPEV=5|efguw=8mUGM68nOAQX%# z>C1C>?#7rd(+80Dn(3FWJ1z~Q#cU(IW5B+F_ZR9a`v!Av}6>NHq zbrewB0!LS6+XAjv1xD^5kQsxpmYrooLQF%kz}I^whZPSh_V6+@6Sd!86>j1pjG2+% zFWlxXDtw*R!Y-G&$^PnR&a~B{QrF01H$+i)taTABL06LLt1yhdu@Hc7~v z2HF|L#BXsNG6+8xhM`?m<-gk0L_|Yal=kiHZ)f6RvlYd12uL+K8|l6&-nXKAp;@_c z<+zKHgNJ*uw}j&b<;x0{fwY1h)e|diq-@~)Gk}NximG%fRXeFWNYnB)cN2XeHUTl% zMBfB?IxsI{qz~ot8%5zlVVr;9+@gzumFqWp4J+kHud^+HI8i}~F8Rs%PG@SU_--$} zj7ymG!EU2#iZQ=02>!YaQ@lrZlq=|HDP6Um>7QXJ)$q6ZaDU|%kXC>@omiEs6+Wmn zjn3@jw94-1~NrK-3b}ge_xF4gY z`!85ID9y_yDDcds-v?jqU5O&qYpCuxu+XS+@4lm^owoYD#vcpOGzfAVP!fq9q6w%* z=i(P$Gr~Ww(~N5DnzY#iJGQL^n~t@*TC&TkFe<4K>r|V4-t!Y+8P5jSpEMR6HLZzJ z_^L*i=7?VUC#kZRxn;Sbar-+t4f|oMNxf069M=snKupzDD|;v zr%_phP`}+B32;%mhyWy#Mf1P)g|Jc^pUT#iXWDN?1tN@6U@sqr4yN% zIJyM)n-^&<$7}fG#P9A;CBS4hMY~AILzm zZzJRFpn0v?WrmLo>~I#L%VqzP{z=sExGC|{DDCMc{z0&n6}qfef1?Uv#~N*?wu*{t zlc#(#%2qhCRPe{>DpI&(!IHwomQ5qA!DqtW zae8A|E;W?l0&+gmHYnXM>pVpav`gx^t`~nDjC#;~Gml|P4nfCcs_J^) zcAU?!kcZ+xS;HG~YG9~pbV<4+P-NcuxD7D$4RU96H6~pNnO!Ov(FlK4 z>2^j@Amu@LEdXaTs6**|`ojau0owklU7u6dsR^1fBBy?B-qwWV995_!{C zNgLD^8b%6B>i-RVt4U(mNOUZ8FNKJ=q+DKEjrtPS$pm%VDn7l#K z)+Oy4z)|QrlJkR>^W62JxI!VmXVgvtc{cnKB_@Ieol)eni<|6K4i4o1J96%hyAvPd zD2X(6_|5f zKc3dX)|=ju9n&pWw}8eVT5v#q8*g92#*EreCQ!4aD$LN|)iM#I9NO^sLX(zzGB^E&^#a4lN|MlaQee zIl}~(rrL!0D7b5`PlPo%n2V-(RD#v6V=V-J8s7a)bmIrM@&+s?&7+n#1W zB>~YnOO|2dlq^;i!rm86E`+{$Mr&5q*wxX~+UTo6O@qBtT_GQ4DIrShPPuv@0Q?Is zY{BV0%ZJK+yq+c}X|;!fgg_Q=Q&}ul6QKz<6r=E1TfAu7zXf{23|Cv}dMq)*Q+-x1Nmiw44MbsOQd82GnLzgRG zGAz`y{iH5;-0o)Gs`-G@c$%xBZcWLSu zSD0ui|Bt;>J zMdmPh;oK&-wT*ljb}w;r)L8WeCe_xTLx3&C2pp4g&#C=HP6vzeD)m9}H}ffzTMJ8g z6_(?2-ajBRcZit9_(Qy!*In?5dviIr55eu#ykeggYA}a65-~qx{~;?*sd?6>+o-6f zF#9H@;HsL5Ji9_#n9C3X4Rg+f^JvY%8O@;wnjqgcup~cQ{k?<^o&)U<{QuJLba0G=dgi`nN8|2)}Kc^c<68pd?P6C zWQb+3n9}2-;xzuX2S!o;fO7T(d%AF%v$;UINw(Qw^^#bNv4bgkS0qFzWUVMMXTv8ya##Cy|seS0Fr%>^DzHV;&$?YTazu0``f^o()Y5L_$Xl1LF}L5 z+Div|Gv>(0Dc9DREVqaNS*?xrD1Cdk4oGwy3$4l=fsm7GFJ+`-&wOC`X+O=;&Fh?M ze+cFWDm(kQqO9wg4|GFFWdoa*TI>P0(|L}y`j}l9b9!xVMYT_c6r z*`2l+q{UJsVww7g$-IPamD-kVrC|*XhC|*_ciqwk5AE6s1x#}6LdYpJ2!0??G^|L~ zsKJ)l%$eXa88d06T{Va-#i%~l)QnPg6JMed4q=9E8~S5wjP$wvQk0Xs?ac!!+>KI{ z3(B21_Un`@=us9lMGm@4sq$#W0VQ77UtqO5-@Eh4e+*B;vpLlJq0|>&62(mIRs=`0 zmO~_ew~q2zTDbj^&13VXrOSp4Rhgd`Pue#uiEFtMrQaoiGks9pks^GknVhA~H7V^- z0iIKQ8XtW@QK7waDDd)n1?2mF}k~@}lA5NGCV&#!MAM zFtY+#`J+PFBg@Ec1e{En*q)yp5Oev^UJ=D-p4aDbp;4-ob@ryxsn#0Q@fx`eqhzd$ z8)I)7`=6es#}&_W0;qYMe z26bvCc6y0x@k-Je%W@+l#iQ5Iu;AhM>)r8L2O`I}Jb5wh*RgMdQLW!nYPlG-8i9DR zT@&6?U*l*J+S#{B44N969s!_M(B-6Yl&PiCvHwcMm!UFGrCoBZ((+C>>KN+@I#y3}4>M)dZpOF_=WwZwnt zKj|M?QPj5CQe-d3YJ%y~Vm_9J5{hu-Tx_d85jb6J@Ty;gy*l)er#+R3F9|B#i7#Rm zyU8F$&;`R*Lw7-YM1N`GdX488Pqv)$-P14Ydv#h7=;$Dvis0Fc{8?^YiGc!cnOR8Y zTU}X^u&^i>BbUQ~;0Y4m9W6Io!q8dh07k}LvS1#pK?+i##6-^#?v&e&ZEfI=aJYp~ zJya!8fLM`kyT3`T&GMxhaD}(l(via*D16;g{TXl3l+0K>9tE+7WcjyIJ~aDB+<%(~ zvn{Qu+hgb44(@2d)7JJuKl;24;)pT|?YrbVT`2p~VsAGygtxL@f>=zj7WU6Q?l`l| zvZkwPOIYG!UDiZ6@@k{W$h^d`T{fxw4wjBm`BmjGSH;BXgLhozpdb&@x15(j-N{O;~}* zYFpg|SNDjFL1BnBk(dmQ=<7pu4PCAyN&V~64^UMw%<6wIZKnTX+H4H}U-2hBBO@Ej zf2D2dnV1>>uiE#2(frvACYQ6>KpoKu>iTbg@js6bRhy-qt7`!GDxY0Ir=Xqd%btU} z!?VUmm#VX^JlHNeb))C=85Dbc0!4RSdN-8Zz|Kr?PP%_M#2<>wAcm;l_t{i3>Y<}hNkmUj9V$0r6Bjt|p*-azHj6+qJ3+u8qq|AK(r z1kSFKioOEekIZN?%?M_U6bo=a5t&}M%m_=ipZ9vt)yj-6d!&*)iQOzGVi;qysD z*tI&;gO&lZ{bg+U82j*GL4HSB`};bN07T*B8k-!xBUPJOUF}(2fPr}+u_zVuvL8f) zDzh=>pl?u^^3DKFACv0Wv=AJ9yRIO5CVJnp9leD=Vn)U$RzHZgmnK%mhF16cx2DG<04RN9zdO-D-_}#B zbED(EgQ-0$;|I~mDBq&DHe-8EMtM(7&AjYL#x1<3aCLSdY2KRc#Jk#YpO!|~)`xGb zG|lz&^q={l^h%Jzj_l+@N;!>7+f;byYrGlYGhjVa1A~Jj6EFY{5Ffmvu93kRC=Ebe zb-DLF{>=|E{WE|!=-O}(IAicOi^b@B9G&&;dE`AZHnBB5e2e|aeK1a1kXKbjEc~iG z{8l_fDCf=_B4&};*?&uC*GE+kqNjmqZ_NI(MI6BAN1*;2Jfnrzo)UKgKvTfz8QXa znGd_75u527nP2q=RF~KQLsC*h*|u1G-#;#XtPi6!89?0Mzv=zr`)78J2|;h?b7=dv zcA&0f!yW2@(UpJlzJ<2303wIJ2xxxuM|eZ9dQ4vs_1=FdNflK!9W!H*ze7zv0>9TL zd=ozBe;2Oz{$cb(e8n5PgK%*ALjU)Jr~#r*RDSM%$mR}oZw~>JzWgqL!`^IOif`2l z_-Oi#ZxG&=fcfixYgQk=6hF6LZpA;w-(*Akcn_jggQJq>rXWnNZ=uHDnO~B}Kl$uD zXYkgWKb1GJKgdxrKLL2&Y+o=wkNcmrukn4yy0?Z@HUEa%9sh<(*Wa}tA7$eYvah-V z^v=%jWsAqN%HMMHKYP4*CV%EY%%EGDY;7m7thDjAzU3&RPGx#o`*|gMO4DVeZ@Q0` z_TPg*vjwyhaoYBsNV6EW9?YFZv(8)R80ITq#3ezT)GQMmHXCOn^#djzb+EP$exLLi z&4=H<_yDVUc`)aD=c(x5p(x&*5%;5OMdsq}E0S&)ksTtlEtY z^d`;Djh2aCfF(`EpJpC?4v_denuhrO!fSh$7ILtwZF^Xq>BZ{Xs*sPY5ipTrsy# zR?J@~w-|SKfgHm&Qb@n2w|8b?V=M;=K^|Z7&8y)Dr3pE=doqXvz5-rRucw#}Qaa)* zDY?`^&Et+~4)wG%_{<+=ohEI_wP@o;@V$2YRS9IIfd7inlC?Woj;qZ7kApl1gea4JURleRQ7_8T{$m)MX%%=F(?aV@9BnmF<>)F2fJP&5Z z=T=0vCv35f5L{6@^xeR?d4E2TQKg|4kx}m=+NJdzLOW?J>LpXzid4Ly4Z^tV9M3K8 z5xnz<2F3qa)Sz;I{)T2N*B`Y@99hPyfy}Ac5Cke;B@w$Bn0ke~J9~snZf;~;|oA<_fG`T(AiUf zkx>**>7-JY2>a&BPBia4Pat{_L|C5Gx~x>{ivCQaDfHK_7Q;%k;FF-&`26B*{uzkB zcODUB&|SBw1Eys}GF`0*@#0}mb<@X`i2(zxlwoB>K6WWxi_Ku$;1Z%l7h6XnZ zy+I>>m-E&-PlnDws&wPAr{k8zfhL|f*-WBK_~if#a-^W!T{Z<87A46JDO}4+s>CIV zgKh3NqI9FQI$rBuuInbwDj#J`D24=sdQF6MIc|TgTCs1$!+ZnuJmG~5foy(zyC|aT zb9;=yz?l5jXVmr3jI*x8#KcWwQVXxCQ*#PttAz75$i%X;Vgdzshz_-5SCXSQfE9Dh ze_rCWW9LqGS5N;#(6;ZI?iB~DE+TStJZ5J8?400i%wym^ql5 z{V+t5<@;eEAy(76@3vpNy3sLPGIkwIfIinEG3JXh#A?3HQV8<`ASE;|5WGdOUNk~^sblhsG;n75JJO3-60c>o!_q)RXAK)d+5~Hwb?6&+DB`c{buMx- zR)BYK!+2WB(pq8Si>5~^s;+hg8R(>o0qKWL?p5BB%obE$Z zo048D`aji}z0UH#;>U7vH=m)H&A*|m2WTo9w3+cuH6%id+!-dg{I4fvH^dB>ezr3W z7x`<5W)fxv%kZVDDF{Z}W>sU{XJS-nGAOm$>~jqVHLtigz+{CtI?BMzNr95ZfRVux z;E<+{s~gO6!xpwnLxh_x%iN(%2tH3pFAQnnZNpf2GfTo7?m}|s{Lk7ceP3txb>me{ z=j?4S6C^g?nv z4ce{YL(pH{+yU-jH|1OfBnh#>EN@R#*lWS%rT>1lYRc^pV%r7L%EKiGm07RkDzSMH@alG3NGY{(O3I$#`fkb4M>kG@p7dD zZzhq+ahV(H=S^DmtoHiFk|&Xl@*pP|=t|$~&EeW!^Nh)zD98o_o zutsQ;FYa@}$O85oA{u_f69zU(68TTYq@pMPZ*4Fspawry>?kP_O0~*j3{|y6tu0)6 zwwit(TK)0^U0rqleJ+d|=Las0mQ}5^VC~8YM+EF3QVgecXk|2YsjWnFX&@S7;ZtZn zEh|Z{`uca~SD%Dm1_QL6C|ylia2~?fxNKz`n+bN{mX`<1hn1n?=w+|Eu&SC0#8@}2 zG_jxMkim#n8FMAG>Cr0fS{OE3=CiITVmUUTZnS?tyF#!^#`iBbEm-ck`a;tdvO(le zLha<3(}}rAW>4W3SpcZEEl~8)xR^P;SgGyOPk)UAOSltk_cVPbWOy~`G}zl5S6OX(M8?^Kj1 zr(^bc9DW@aL4_N57RU5lnJ4|>ZMc7hOE6bidG;i}2D?_wMmN2tlanaSrVtoXEH$1M z#SK0>Sj0bF-whm>77~y<;fGuY2hFg=|4n6eU1rRpbNgZBNR@$R_H=ED4^kSjMRF7D zNffc$QxF|x`Me;OcV0{H!+1{>t5OpPizl<9DD$^JF^^ogSz)j%f}KeYu$7L1I)Md- z5?ql`>8)y(d>mEjpRYWg#=HE_BJ6>&@Vy=DQgQ|PXkDv0Bue(r1F{6C=+P0J83#~C zqh;%MkL<}xhCP1|rb3Zc-z5W76jCwyQ*}4}70_w>$jR6))|8lD7oC4gFlEjpYAq0ckU~hQ5 zu;-XE?h(6vhc?AX*q$NaorZEq@Z>d@;&YQ~r__1{(;mIN@&(*dq^;!LlL|_Lh=i!!asGvYden+i@?pKBn2d@k*ZXR*LNomWggvo5BaDiPG_Jg<3Sd|kv=1MPO#C9jqE-}ot&GXZY8QW>MH^ae^U zyHrldI=hLc@U9&u4#t37(eJvv4(f+Kxe`GK#JTF76HM6zr)yidJ$wcH!c?`i=l8QES`qz_C>1bgbPJP`ag1r^eu1h2Y1_aw)~u-r8(L?Mc@X8KVD#81W1_F zE!EKNh_khi>GK_AMvcDYxofW5WN)8BQ3A-jzGYsnI_j0dza%UQj{xpd!X0H5I4ibi zXoCqX-^p_YI3mG{=cI$1ik-`TxadPN!>KrE`e0%n>_VMOWla2HgU(VQgK;P@k3%E0 z?P(zPUGySB&5EYQ4FviZsn~rzc-xgi%uv&Iz~iATyvaC70``>9pC;v<0=`FFEU3DS zBew#T2XX}HO`j+0&i--bji?+VFJ-T|L^8uf>;s(h5%mR^kf@q0X`kxR4aW>G$hI`P zU4%K?2Q`*{kj90FEIFa=u*kX z*g4GtB+5SSL~O$$=N2vTze4Iz*s8il%`Wo#XeOa5-k_r%^RKnlwOC`4B&s)-%+ss& z9i!!6zSehL=R7r`gKOZ@NN}Z>W!SmF-HCEQl0(uLcPO*JA_ygY%rvNU7QG~1vxaFK zm`JK)Pivuvz2N_Jr&|lx3t=uQ!zDLTt(NSkhO?HEQxGc>8@*kYleSIkDeE7nc?50X)HV>nNTn`8U_i ztY4Dl6i{5?tV-;S19`bULh5qOjn3kuyGhIfxt5)=laR#}q5ek94s1n zVevQ1Qx-Q#?^N!Z;N0m9N`Bv{7w5z)-W&mw-DNz&d_zGT3WB`)2~Q#q&+dt^OP{K9 z!okZzDr?o`R2s>wmRqut*Vvu(l*Dv5U2f;-hUu<8`W|!bb30`dlrJj#yL24M-)P>6 zX`eb%ggjQuDPsYnfW$x&F>&l-m*&xdQE0=new-1z546_2cLg{IHWWEMjvX$|ZlQMa zK_ZrRx%{eQ6%{t%8vjCBYRUAiyEjy|%(jAV*({hO;?Y}}Z2;dzR=c7SgQ^~ZNRtKc z_OVgXXB&V2Iv4rNYc}a8!Cup`NVvk%Vs!W>5Hy`0r0$pZ<*D=NsoZ-dU)UduDX!`w zI8JpIIzoJ6J3B@bHSc(9n@MlDnQE75h3Gt`?6My+YXIB>C)f0Ga*Vo_VgxE%{FV1E zfOa1%5(&;oozURgb%`P;d>u$?jnKIu6N5HTTy-+wkF`OjQj`46ZQ0qqinEAT%2>7` zlN$p$G)$F=C{X!3uofsGNLvm;=4|rt?_7DSX7prah$eL_X}j~Ogx0RU5`0~| zF|vJt{QHv8WK{E*%QkgLq7NDPUzd8r;!}(G)F6rL>$-zR)y&3Jel-WX{8BB zizL1WpjC``knSy}VYz+A*WjvnKwnvYQr9S1^Pbn$X`N%cUgf6Qrl0?gB&amO$YGuZ#_NlFn>UHIvrT4_ zyG=u~k_x=A^w_qDuF_8lb7NIy4g}8BmdGo8dw2pmYPElmg@FxDusa%HsSeXa% zYl|HJWLL!XijM57)J9a6L?iKc1ZY#IKF_eoJNW!b17i?(clFr)I3x8pF2N^fl@;~$ z`x_-s{-)Ao-uJIEdTkDd-V4G4V>?#s?X^r$?Jo+I$_YVGp9j3q^whftSuYupp7eC7 zxp@=iE-uY&j~w0Sw=qncdkH?t!z5{Svog%OYlqMhxdWsnpS76V8Nv?s2MSd;9liQ^ zDjGwAQH%DoO^Oo%KXomX_`$dM4dCTzM~?<73xI~5!wS+ed8!8Xr(_i^jO>_W^sZ*> zha*`gy>N=dXyAf>16y*)L$HWrMA61Nz@bH8=aL{1Fil36Q2XkWBfVEtk2@E5n6JAD zyiCK_?S0qWDZjf*a*CL>-Sd!m&HtSQdPB+q7NH=)RM{Zb%>j07^VhPx{ce-IfQFwr zhcRC`9Ht5+*wxKUx%!&9gri5|oMlpa+P%yVj^zOrs-9%B4JBmw#zJoiydrM#F$mw~ zx~JCZcAo=!`T$+AAg`7A<#ywio+7SKq&BC$Q*z}QgS!nclE&25c0==0#9j?a_g5oV zhhae9vZHqaSu+wA#9#n%aNP1Kuf5pd1)?!1Ueyj%*X=J$be=$1n9uhSX~OBFNcPP# z(SzHtZ34*gCh@apQZ0$q;2wwGtR@LW?Fu-6?)q^?sVe=Q&;-Um_!& zX|9pJRPv)xH>8d3c73K=dy2=>T?lV|7E^FZ;_=E1bX~hPmJAo6rCW6`-=qa8zpkd0 zluoXF)5ebR$Q9CLuS_k1s1HEM*v6wTYcDR_vKcbia62^!@#NL5!r~{AV_w|dt=s*j zRiEgG?Wnx)I^WjddPVS`qzou&ay`m)t}}(Ut|X|q+{%aN-`?dnoNyHPIQ?ko)x@sA zX|BD(OiMC<9^|W*TV9LO?J2S*BE@MPpHWduKXmBdtYZGFmoKAB(PJ6LPf<`{y3 zFkc34rmEs7N24M)PaWQ++C=NItaLFlbA)*m~tsWXW1&65>zJAFrG^Bz)#{;{d{Uf;ULo87xg1wj`kZI2OYBn0e<; zOOu~lvvt~*GboA;qTvdw)@}&alU_;I5To))atVLPkHR0yuw~WU)>f=N`9MbAPNEt! z-C{;tm}(3vAjh>i+bdu(YbBi#TgEtrS8&G?X~BQxXFz&DXkF)fHew*J90|&rIPMV0 z$r0M3v4|zKm@A=l%YN||Xx47u7(}oW$vh|Y1DYVD`|YKfN9Mxhr35*o;Y4$y#9hTp zC(KKv0XTecVZGVU7dtm$WNIji+JdCfRoM)S6E#?|=*O1jGCe!E)CD9nv z3A&xUKh9d@f8!Y!TPver=>*&fv56$a5pww6&NU%w*N9>tBOyU(&nQ=E97yvRB(U zx7Pv)t)rm(c9YyR-g`NP>}(w2AqXcnzNC^72(q7S>3YHL>5wn?q0NndrC##)ajTx& zYCD1FlZB}-@~CuS!5Fx0J~y5m=a{7Hc;)-?lvjJrj2!qG0_Ksxp(E--CYSf!wNcbX4QYlK+ ztu{uveLmU>t@nFE$?1t9PlfO;)>xwFAQeAA0$MY|?9bkSBy~~gnZpCq7Xdj1xm{NI zWxe`$hQk)CrcmvUy$+?CvF=t&1NdvLybQmRr5vVwGQyd2jR?--%rvK;Q;}n2oJgs% zBuq9cO4Y32PmYim$o!{dz(vI(t`?Nk&Nq1_+5Kmu=X{}AmLB<1ZRfR3fjdZMe{57+ zDyMA3(+dWuoQ0~X!A*cnx*juBDtbOkMowSzeqV_rrNlNr3+YY3=3=Yhdc)PpbxnI? z(1Eotn3nin=8bZ8V7lgs?&c%%PHK4@{{ag=2*y0ax&n$E$lpT%nENEt6vRSCBmqf!kIG)nZ1}9!+p7JYu8J(x2o??i10ExIl-S+Ix{SVfO9#cX1Z2SWwuZW>enZF*V!*-)?nAsRifIj{<+kj?ln6^cK2Nho{khg0DI|M-DDbtX`TH5Z zdxCbdxy-6s(`1l7G%5mA`Qp6C!~Q}FW~XK2wF+Tq9i?Opj+L5T7gINlCBJHPeVOER zv52nt?;%vbGET{dEHWvd!QLL6c9&U)n6s_}YzJ?PNYU(MggPs;{b~O+VLB)F^~&;W zr{Y7g^|x*HFcQ{n3K1x%m6NIZfD?K}Yd8_n@chLFID#4msRF<S2|z%^9(=!X>QBQp#spI|7ihO!C2HD** zRslBL?T}wPm3P<@SqpD^!7$gE-=P z1koT@ywkEeQ1?A2Ey#0&r7L&5wR>2Vt4UN?_iKenSkC5DiRcr zO|1Dc|1b&h70lMR6uL?f?3p$^SsY^OPkBURx%_44vD2U^s%uIi0oEcExf6mvk-?FC zJ!bK}FQHD$Gjctzca-AB&D1n3t-cT)pzk*MRO2LvuF2ZLL*Lhn14B$Q@H*YLiAU=*4YAnQ=LE+YI7* zYI;UO&SaD;B)^v$6#`GJdCfCk1!+CA!kmTf?y|;W;d-5yf7c@o8-e`V{fQmS#gwmx>3cC{Nb?zK zT#Ly_HisI(pJKSwZ(5iW3L3(cniS%e5C&5rk%EAKM<3vXc;1jpi7TpLnrj?A;?FR` zb|~G+-g@e`ToFTPW(p4fh+pw}4UgB?f3T>D4@L=)WpOc1L?s~AA28^@#umHV5BB)W z54(~b-R_)-sB389{3sAyNXAi_D1h#b=$C?tX$bUD1*lT$i(J8wi^ z>@&=!!+#${`Qj*vHitJ9D+Y`n=VUn7YC@-wgF7ls*gmMD&UBPL3tzQ$cRW0*B0uO4 z?>d4Jed5@)E3g!Efncfpfp@@TRZlz&F%H(MK%6H6Wnk?qUqNe2 z0_3d|1Y_D1kNEbbrnpV-W#mbkE-5xsuVgI+^l-k0pUwVcIUcAp`jPg>zPwUsA)Dj61r!Vrp7wk`!cY*}A2g4c(qAsZy6#M@x5`lUDu6 zr}Jh+`Q7KH^Co7Ze$m13q{f#h_niH~$<`Z)wNt^$=z-gRQ;tW)1me6tL z?kWz?MYJ6!_NffB6@0!{r?enCRgCdusd16WV}cVo3@&p)aej^EhR8jLQeT%c?h$s#x_nf=%N7F(JY^pJVs;S zcP9o?;8oc4tgPkoWY7iY**JASLQ448IUpC#LUeYQ7s2Xa^47ix<>PNoWF27Le7ETo zrI?5h^4TM?+Jn^_UZ;LlSs_t^RwvJ}!VY!8NF(KQYx7WH&%3ie-n~;#Sw|UOe4f6L zf-|iKo3&|#BvTlg-8fPrWs4gUs7($Z>}q$>AXV;Dyo*Xtw&rL>=u{c9t(553#NF=? zF?MvN5}6!S+858+;-IvKAld8U;p?0<`l0QB49}w;1Z;&q<5;rbjmA4|&J8A^kNTPB zgZwVa=$P+YVWxMCR7l=(XnYk3O@BClu6VV&zF$#K{B@uZ`zT5>v|9+c2u#2Q)#SJagg z7U#O^u$~o`A+-0QwGTNSd34MIc6|+HDg{fjKXbpl(38xO+tQ8%&R!YkQOw1LrDE`ep6W1`z+$9gmpWinAYXN~+|L@! zYkLcdVv~>dAv<@>Tl9$`XP}A+-IX_3#|`13>FVfp5^VL5OUBoUR*gGtUoP<>Im95u zvTHx`@Okw4g1|GDUn-K^*;$w5L>re7Bw(oXz9(Ll0pY~M$HMH1oO5Dgg6diybWD#UEGIfu0t z({#ZRW#3XM})*Kz7x3hZCB}c#hH=xf|K%tF+>u+hoLroKDi6-ed!w zF?WsK0Dp0W3O_I(953W*?z#evyxx)wOJ7+bff1Gl^VD?b6n6dTtzTwN?^B+Y!>iUn z%A&^PxhO3=u-xuSR7yB52H=UXg1gUO8PAkNU-L0lsAW)b^5;-nqgZlj1!W}VK>wP} za-wm@-4*p)jE#UikM;j3#ZT=FAfTG8Lq;QmOQ#4Dh3S^PGOE4zahw!dYpi++6zBto z2k9(R^0ThRHT;Hp}j=YxRqnh?)E@1}}apz8R}M z8nGB&v$am4qRo(hK$!kP=30H7@`3s|y?yQ^7X7i}pgU!Z_o!9zA+@kQ%5YQOq92Y~ z%4_)a$Q)s!*NamV@~F`ABI%fFT(Hd8*nC1WCzw1Jog};$^tSsWy+Ar@Yp=6+xce!~ zUBGO!J<7T+Zi^+4N{+X0`NZE&|J?M?-vLvMSFF|;vlN@a(O2dehRZJvP7glamABLh zNr)x1dc3(P2wh}ppzs($VZ-1mC1)z;Z{NR#ufI= zn`2PFHX@=9^x_BV_*j?VG(-p)lkZo@^D#S@)$;2g3xnn97C6x?DQk<@tY}JFtvrTf za5nAp*uUcD1pSP6*TaE0d7fDpLy8zG#}jf|BK#^nKYjH8eI%+i3HqwTe8#R~Ru^JS z4k^*M770;;SE1!_+FTa@Ejx}Fg|t(qD(W|I#WZjT3Lld;h+Bau)lUy_E0Sg#wx?3T zEvJ`3uKPCUo2u>vTl-ez{B0j%y`(UgRgDjP)0lC?B#=*C3$iLvCU?Af z-^Mm=8fVD2S{Yhzv5sU|3Zo?)V3tXVO7Qc{UTCiA=RH1Q!vokU`3#s`@; zbLj`~F#R{84EkrM0h_*;!~lvoOdA$6in0qyZLiQhLeK$YfRYvnRQSm&MFzGVoVu{a znn3-!CMF`u`zd~_x+EJTV`Jmwi&=f^lE$qzJzKs_6)jL3k)1Ji3HUp0RV3gxU(|~) z7^6a}TK(|LF8JQ|;_}s4X%v_cUCOk&(gY5Kiq3?j?&C6~LN&B}?EHEIQ5>gwJtH~o zGCJ2$!G;5^xwnh~AhT1@uimy#&mzcM7qI(5c}~g)+7q-+k*7?$Yp=TRDe+>w&@I!A z-^m!U+uQRT0{t`YMQIicZU&IgkmiB`pDmkt`4&0%aVP6Z*z`uQjjBQ*gOEn>N(dBj zXw9J1=E!afFNBzXAW}9$=R zfK2nfYxoq|B+zfWQ|+XFs7JzCsRWFQ5!8UlbyfA^aF(THc5PmAlQlrB2UVO-dfPba zeDgnI^qetVH}CWvZfCV+1BECyae_z;3bhj*Rvtipg&L&fu$fk)ZBD9wyWtan*NI~`6^K<1o&QHXZD&!WX4h#GtM+W4ZWp*^VL z9t5ASO0k=s1H{``+d$zjXupkvO_K26u=mi0NchXCEb~DEx!E+t%H6vm%7?V*vzRt_ zPPisq&K203k{KUj;Ydn1=Q1U5gc&HunxMclS;SvnUH9R*qukX|7@p=B$TXTAZTsw{ z>k}CXM&&jSM`@@1Y#7pobmS@Wqh7YKPi=NbWGs9n>-cviGXXAZKykz4v&b`lW8uxB z-(xaYgdJWI78pqLxSy_B=q3Ie-lrV3xg51SLz}$dY`LeFIs@53RB9C03Q39^Rfpn> z{3_~QVVHeMsIYSe3#U{#Qrd!H`{xQo_w&hAP+VR0rTQspwVbE3w}&U($_ zQ}ayA#_*`YmIka_cH}l9K2D2agI~>m6pa%)j!SJ=10@AX<>FLZJ6@T=x|{Jp@u1b_ z24=SOsjZXxqW?7faU=05M$p3!w%`@PdC$SI%BvfE3^QeL zSv4R5=AVouD}k^@((cE)Kq(Tk|2YXpGt}(S+~w|)DAjW%N5jcl)(xwujy&Ks&+jty zbRCU)s>w~w+pwcs$Y_Dcxr@3Gt%JCIxckF!L;Ic$TLtI3&ab21i{hz4-ngOU@M5*2 zdWVQHmu5YSHWI&1UJIxajw-^YGc3Au9B5!me_K_7AEg8{C*QeU1b0BFnL9JB%7x%z zkucR5G?gr4pNo6$)h$Zy&}fvcQNDNVo99j4INBfxqfhB(2X;h=Q9}jUAel`?Bh!gq zEpUDg(Sw4&b9O>Nl?owYYmR&Q;c_a6bPDh6Ge;txmK)9P+nBoHLmIk@O=5g4e>3?c zzlF&IyYaIuZ;Y5+^C;=*-x~DIDUSvzZtj527^#O1VAfbS_;+H96bM#L3~V*ZQ+`!@ z?^uE=THa-tC4}S}V=^R~p^O8Na-*CAp4Mp>i2QT(&f*J9Ap0W%MNC1%();3}x1 zHPce@u!Sg2H&P@IQU#DSi05WcUDM)5l|O6-Z$m#%(PsNXA?>E`@6OKcE_rml@%v{8!~-kk#2RbuWkX?-0we1%)@R=WmT6s~j|np2O~m&?rg z;pW0~^eEx5Krdrd28|U&0G?k5pU%We* zV+fL=6uo%7RAqaM`FawYt-vAN!TjrQ>x)M& zFwL7o6M<)Mdfe5M=M|}`m#Fx3>nd6No?3ru?tL^|WX_2pbdWjEm02CZWY-WJ5Xc+A z$&ag*nLjG~FP7|NIHV{LTH&qvgcJbUfZ~O6q60 zyr6#7y|E-(@rwd!!;i=Xs=4l0kbgL@n6Uwh8N0^^S2}z6N?>a2^*1Ld|mZ+s&u61-m6bO+mL%pi&U}0y*AWC}dpcs+C{k98& zVZv4GCPtDeVnbF21Tc6uV-^elg|{Y0ZxJu@{_S}E?NrH56essugV~ekD4rNVWpXKtTW4=mgEIg*k7+SJV#fx%mm9Gq(?HKrp zT0`IT=}xtxQ$9&y>zsBbx!6x-{`c|kSeL1`nU_)zIK*yfHgYp~sWF&z^GJ92wr@@| z;r=(2t=23amE}q;F*E%gxKP2OT<|52=Rzv!A#?`z?AradstF{aWt9v(DL^x>XFbG{ zg%#X_L=XrRs22JL-ekz(=#hOAA3j}+zIA^#6|rF{>4-syT+e__v!9r!!G%gm@N;13 zlD5=7C1S2$ai48dg8)nR&ngDXKUma_IrO9p?YboVRB)(0X)U;pRIAwvs#s~})Q04%`9fJvP$c75^`IuOt{>DK+34!`MxKYokzyUx& zKq3Va`id+gvWfy282)`7z_R;2_R;)Xpl_AB`GhEtzuV~){_y5Ev3K9VSp~ZFfC?}; z;ow~e3>*{a_a&BtR^dk_;cKa1D zU}By^3j+}-*z>6epy5~+me7a6j6wnwk|*iK%hKqv@Zq1sh`RD;+{1Q5gab0E$OI^Q zdV4z_LktmiH1rkqgI8uEF zS`idp|OrppISd}_XC@3mQ|lsLwwo4?2;<} zSxLi@5tG0o!zBX)5W>?j6tH?BzRl2#DGq9>`us*O@^$6~*7=7I?9TN;KYv+$etLrm z0lv)WgjsUy?FoK&L0W;O0}cPx2Y<6pdWC=3U4JJY{1CqQFXx{gA3ki3|1|^g)hLkm zQJ>OqXVv*>azDY$yMPpa&n-cJiPeGIpa!Bi4Yjx_^sVl`+4s1VfA_N{j&Dl=f55$U{PSMAEb4S z5|sofz<^_Zj+@K(JgNYep|S?ck|Z*MLRTz);2B z)vB04MC_H*$$E_~7gPe$G!Y-ufm4=tTsPz=vA#Yh7In3hGxG9YwX;URP-(tPNrQRk_+2p&)Q?mm7{7btFkTP zv(_6md?-o99W-49Nd$^^ZRZ>?_VApR9@1=)Xducl34Uw~WY&KNDSl)ZK+|oL5DJ|A z33Ml%aTUDkpLf5za<aRa4Vb(?*jk^!fNnxa3%VX3;%~YfBtx4`k%-E!4P|uLgP1 zGDMa>T6ho)C&^Hk;uTp{+eY)VZBbZyE(mTmqAFR;t6DRCY`=G1LuvxbV8F8XD*nenEtL z_upgv^QQ$q@akx{GQ@~yTUs2FHyUUqi|^eUDe=}*4_!{rMFq-l%Wnr8xN{;WT8WDS z3#T2)qTH$R?ZAsEcTyR}tNNPwCCJiNd6=DR=1-dJk# zIwg@hDWZA@V)ic=d0fu@_N-iW#|>9Esq^aH~ZZ5sU@geswLE?;GdF5=GyK^^E1dluujni?v2}+ zt9Y_IE5$EebaM8Nw;#to>G#sh20G~4J=6gMUxqQzBsZ5;D@oXGaATutrfGuDbp%#2 zP8FUZ2{uL%Vzp-q2jn3!Ccb=OiDdEVrdGI(0=mqbi?#S29sG~8x3+|CSQ@(68j|#j zPU(CIz>OT!!z-DW#IFKraC}(|Jac|exx|4~8g>i*$nTL>`FE;@D>vsP?eT+Y9gu!0 z6|DeGV!Q2X1nL?5qbI2@)pGezhPUvQ?3}X_Ua_O7z@^I@lUi&f;^D-~FRnG9M=@O9 zk|jJw*;h2#PxAZZ8`ud5Xq4gpt$tqI8yW&7AWSj3J2L8?`Fh8|uNKl)rCyl|GlFWQ ziD}wGXl3=yut(d!$qDwPjQR>V->&sU!~=h1Afl0(POQ&jbj1*1%y z1-jDR7^r`IG+reQ_hu{ev_RmBs~L-80cTXg))St~~L#r=dEbR<1Q z+cR&sL{&DBv7QhTK|cTS`05@|QQmg8bqnniW2)JBEgs-npgWKgv5dmcNefEtG=yIz ziFgcBQGP%fr8YV>4m#szKLhPeT;3v}qP|GfQvB#@2(VNfuO2lu zRp&cFvh{|E=r^^*8M2O)O9vx@CDp*@IaPnhgapE1y6hh62NQ2H6`|#1u8eZU0&%EP z#poO6+RyT682#Ka(*|Cv&ewEo@pwXuIY8V)`bT!UNdIYAxyo;@hOw9ac4#->U~?4- zqYDx>ZnF4P!WXUj&H-7iIb@xzMLirpcEVdoRS-9~yNvU(MZ}q1IYWE2Wb7QvF(&H^ z*un~HRQ#B#wgYnB+@Rk{j4jtdmb0>EDRI_A`E9bFOn9w&p(L9#yDn3r)z5LKKfCvC z3Yf|96e9ltBfUNTb8Nx+x6B-6j6A1yY7@E)3(lyt!Rfov3!jHdC*WR9$Uf)!5>;TO+QDE%_V8|K88#=6BmL9gT(t>e^GWH^$BE)k7g^a38BOaRY}BQ*%ijb(zCyzY$H0k9?u5-O^1r)QlP)6-n;#4=)$DX+q2_g6a$;Od zI~suZ0791ebld&P4Win+qS4j*QYtYOdoZu*xj;tzq|Lhs_H z-is$d8c(7fH@QPC)slZp#MoJ5hz&U((C@K%JI7vT3ZL!cGrbb z`DCRIVuRDKr80^WPsrXXQp?FyzW+A0S(9_JCnt;WYR_q$auD?j`w)DD_gBQCC@VD% zXSM3W1B)>HjkH!d7$?F)=!iCDY)m95cJW6Q5<2?sOL;U+xY_9p*bevT*;DAV@4*4B z+LJD6&c6X5Xosn3Ql#UwJflqDq0=p`9zJK`O(#(?IfL)$Lc7P}t$O+6%A6=IDpR2P zEMd{PXM^2p0i#QrKRsQq&QW!8R2o+%27wTAgD=IN))Ar55bJc47M{HgJ^T}F)|!i4 zj}R-XzjOgFIKWmze%ZMO7km>Ts`)+k7kN9L!WV(EWBaU+o!}pFTGORi-iOL$RIS|d zvgTYBs^|AU$rTJI=(Ag;>$V8rcJje;hgS$^YTxy1!Vv1fekH z{Ohw^&dcA@)YN&J;9@Sv__g{>Bkg#oszgEGpJUKk?!ydF*$NA9rlYXPdtGSPAi{2? zeyi3aJIx#@{K;XRA&_F)eh?i5Bnr=jooRRYG;=71)bJ9XBX($wKU2PFphdKGf;hlU z`(sU6cmk)hhUuzdihL!<3+}#mu1f$&VUj6uy%dKH9=QZ-f;!trIJ~t)pleyYF1{v7+j$*A7MRmnbfAB`` zRAo3~j-=Fi!%I9;EGf^VKjEM4M@5(h!lmh=h$bX{O9zEMA)PtNRQyjP>gvuvDC2p} z?hSI-wC5-Lf)S4*y;h-H9x%|qcCFsjjl9hq29ZKs>0@mB2K^A zCz}xJQS#Wj@Qii$+&x}3-@Nhf1xFXJ*S^B?x?kz@PZq2h5JcTd*p|d-DC^Hi+zj#a zzj1DzWDOYb(C+D2+tX*VsugUg!(68lMB~vTDP87mJ|}M$VoSYPqN2B-9yiF69HEZW z&SsCvN=PVNm@Pv%L9mmO#OBzz%WApu^HOxXtpZ3Rss9i|wf~z#nI1NyZtq+xzC&%Y*CT5bEDQyAiC)g4m z+(jnL}R`Nl)|SB+8bjTpzo0$UnriJ?XAI zL+u&Gb!m(S;@^QDFFQ)OEqzzVX9X@Gc~b69`j6s!lIAKq$lZErVj}bN8ce{jVK>rC z%hzj(MmhvXZ+mh9o;7eQq6%dt@_8-Q5k1=66SLCL@=_L(c3h}pcRC;yY8sTXu0PK) zss9vp3uGmh4R;oBFSqljX4|S%mytV%MG8weorF<<=IxYLe?`TqEwU^v>Yr3u(<_&J z^mE9VKgy5D;xnL;)@1MpnW*yv&rsbMtPsO9Ew&sS*jd8qFK1HxxE@wzDBrlKTLeHa zo}gZy5vqjGZ2Ei#BnArL&tX9a(|k1bqbr(>U0&s&{q%RL?F@L}Pc3GtR8+Uud#Sh- zeJ+w>pM=qSg@#0ZK0tBRcEy>q({~Uqu3hGQOyQqpX~7V+5pR|Y6grC0NVB;b3r(-1 ziBM@f<3)AMT&}Ca8KyttN@_^QYQdnScS)U2YXIZ+{%KvFV7=P$F88$~*oGODV3`aq zAtHZm%=z~%_eI_7I^7SBDNYDpAF&>3bh}$Zo2FXv%rXmt0|7Vo_1cnM7e(?g8Zxj6 zPJL#&W8xK3;tjns*`D;D3&tXGNWR$aW@lr=#x)a&COq!56>nrWdZ^f@p)uHje@7tt zJ3DaN2Byy3tAqBBca|c+Gi`-(L5>4R*QyUL2rfQOwR*r!=j__F-uk|&m|m+3dzY|h zQd33fX-s5la)lDCx;lIK573Q$%)KY#{CTcImGSwB0@ri0)M7E)bzzIm<0i{)KvPbp zs@yIymJ)SWz~TuSM@g4%Y9JY5cXqFDHv3T!wpPpe1SCR@Q;hsCj6p;#L z6U4GX;gw|Q69O&H6-9aZ6_+}50=U+$f8S?uD%7C3sqEA01D%&(3~SC3{@Q4K=pSa< z^&c5e+PklUXg87hYV!+!dKcipuj8$l)|#}>|@vhqUj z?U8l(o0~tdZ83Oo1$t?z79jGI+`V`-VM<>RP`QzMM(BlpF~4z6h&)fu)?j}r6SCWS zaEle_pK-n`j{b4Hi_8v}P_t`{fn&!b*syI3tge;#BBTg|Fl~se`_Q26IeM;tI57Rr z64c+~X3ksZ7?uN>!vr)EIiSc@0w?Y9J`$jAC_Iz^V?4_S0p70^o2=bJQgg?;bGM9! zg1QtwmGsTd2T^dpA(A2%gL&Kno;O2-f=HRwQzxHLDEpSzh(=AR1j&nxm^jB06 z*5y-$lZ4>ZtzofgJN)?8+D4e=O&wTzhN(I)zi3dSH;=Fl@=G6w0K#}kv#4D5wTmDV)ukV<5TAeI|wMfKyMu@BTX zVHR3x(Y;?Ty$-gosW`8uzmbd`n}rTab%~K<;0%AXC{}{hYpZ?%!%vcO=y4Rb^36V^hWTcHXp836U-@n|nwXU!-+yA!W@f zS%Hwc3U~y+Ve1ay9`KUL47-GAqo%A_TUoW8NVLa z=eWEyNXbi|*hI*b`hKQqTM&-_C{89U8pWdpgrTgPy*lb-{tXu`X2s@NRuhB!s zOq)_gb5B(ptz0v1k;n^=Z2GXT5e52%aC&)PdAA4TYm@!El1!!a>{qY^v5yxJ0a%Wzbp3a`j1T2LYz5znc#Y*4XGAD4j})_B0rr)i~+ zT!2ZsNU`?%q=CEf7A)f7$ZzyZp?T#1fD{Gq)pO~f2Keu8tN)Z)D>7uq~-*i@>k5!Z0rf554^=qe-7kFr69&=ydnu#)_!cGBs6AL}p& zbQLwUB(koC6tJ^(E*(~}i4t`%!|h#RJ<(+fM#q}f#q@>UGyWR0tKA=IP3L+%q$Cbg z>-SnfjQ2*8P@gxJ1x``lbYE-fSE?(lG;0kOZLmC`071!1L!WZ-UDb9Jwb`wwtym#r ziYp`j{(i(ddn({*V?*fD!9q&u4WuR zyRXT3j*aQlHeT=K;6ulD9%f~$Dc@PSuxp!+#GP?A-gVHxy|tYy2kQpfd3)Le6gsJ{ zGsbYTT?7@G?F}9ubKdA>#pBmhj;DHiSBPE64}RWjn zrb6k~h68!iam5y#dT2zze*x|+G7W>b ziAmb8f6YKlIY+5#o!7eJ(c_oSI4#@$53WV)dHL8fOpO=EfQLmI58I2A@ICa_MRI#z zs*&4s%|=nXL>OJ%+^0@rZ1PC+TA6?Fd8l2ofMe0dF>Cbq$+T(52MYBe09$6zg;9Ne z;V@#-DPyB`GJ{wnGys0;-87c>bYjN}*2+8*sr%i)gat*GN5add$u;9J;j4kTBhEPI ztYDp(J@on_CWO+H+)z{YV$jyhO8M?R5~?=@t^Cn4MF=FZau?S?P4FrhnMmju0I(maxhe00B_m||BUURUQ=4ngJitp|g# zuOM|`E#ANkW;b2=8IJhNfz{H2`yH2G|DQmL(Mled%5V)|FSS20acogCSo@+Zw{Ans zJ(luh6cR80-V^GasR^AG?%LFQ>xZq5c5MlPE`IEvG2svnu1F(G)~rE4A3D_GE9(kz z$s(G!?`Je6Tb%;(jsUq*S-UL@mg~ZAkc0gR8+sZp?)f|`wJmzhKEl-|H;I!ml%I^g z&98L#XNxNN66aHp?I!z3pM*fnjfY8+PLt$-{q}Y5zOu10vs5T)+MvmYlVs1k*@JRKR`3#9Pj@U6LS2wn2?3>KVa4WjR~3mb7Cf7 z_`k%EfSr-y|L+cU2UStI-e6lGit&H)cZMUN-QK2UgaL+uVIG5lU7*~KA;1M%sQ&BU z)D#%J5aSQ^5!2Iin*I9qtF_zur<(V({pPafwEc#+kqx^smoTr3cL1(1(7)&nWIm}*W9UyS?mnhIc0%QaK2r_U$z}0-VtYUD zC*HuV3>#`ZE_lHAnc|Kb#xpBr|0NPaT^#}^gyc;%4`~S(NbI&Vt}k!UC4jhV&)2UO zC*R&HOy6eBs|^xkupo|?;6loe!#sR%H=P&+6ai2wJTejjh#(FS!LBZ;uNFj%&o+K{ zWS5~M41k|K7`qQ}8;k_t0zl9Y@r_V07Xe^m8yIQtH80e!1_A^W@I`?E$^hIUcmVPT z7Z+9-=SQn(b{F>owhpB5m>dG&=i}>g2CACD3k2un(eC4x7;TP4ZAoe6-=h4`8w~{o zVdsyZ7Z2huFN6YtfP#Vo8Xge=0Q7??4iECJJnHvWF^G!*5Yam!db_ww>h&`ZaGD3W z5A?;H3I|lJ1L6O;E93`)0vI@t-v8UU<=6T3tMWrV@k{meTO%DYga*eR6D;Mx+cu4^SfCQ3sg`lh<6jN#*()O?7J7#YfqsB zdI4L%8dTUP=U0A0zxB4-)t_%4jt2z!{TwXl1sL#aDh!Xg2KtR6SXlhN2_iT<$Nxw* z*bnikKS)MG1O?FE8RQKyXBMLu0q_MfuWSK-|3mx7AAlbxK&TU7ju#972jK>!7fua< z0t%rHxcyhi4_Lt8pBm*R;0v7~7>vh&q70kskp2q)ydUr8=1lXe8rj*#uhtKq599zA z5SeC%uusWwhJKK%_TuQg>y0mFzoy1ToIDd=BqD$9H1D0$V&Akv!L}K z>4^1Ie$>aE_2&D`#r96BLkkc}*uy}0KfTue?Q-2p9HNLfaK%B0hS)Rcib66kU}yzt zGm&`M`G!lUO5Oc(`1Byl_9!)E2(J%S!BcP|TP`dzY4&z!(A#omsv{JO(N6?R{>C=r z-0-h0k3-Nq3YK5_9H#W+ve_ULpViS|7<+pVFy=t1Xr@ndAuA#Zi<=K?`lEPKm^np~ zBdbku301)cVcKTZzhflO3>{YpU39aygx!l(#)+c zXL0NL6ypb#fO9P?=SD;;U65)(^$Omh5y@kj8iv|jnCvh_pUs)xPdRT81$8gW>VS@A z?64OSmG!sB!%m8gj%#9mylN(kjQ8206R09DM$u+j(D}miAPWE}Ep|vQS1MfNW9I0Y zkKPacEO@@^9~sKK$LkgL9OCYTQHT{`u!-|Mk;a;s<*eXksNZHcsRA@A;TYMo?*kP@(A3 z%Y}A%l3`6!r4itOc{K8PyGjeY@&whpZAH7EsH0r-Gv68l5W8B9`I`gdEiz%kz|RDT z8(nGc-PHzsxZX5utxEZKgB)|Agaw0)lV%b&UnfSPqyQX|Vngg9_hYcG30xuhtp_NQ1V{_%4~|6g z$4d|jI63%3lZ2W3$z;@=A}Y*h6Q3=tU^EBWRa+0E_&(rKuAiTPu!MW|q^lAUQ15SZ zJqK?2q0-MId+rm&pO*0X);AlNd}K;d3b#YCO@9oUuV@NWs>XN>Uu+UGav$nfV&+0) z6gG1*TcCoqGVuclnqOI_xA%2pbLHT%e&ImomxtMvamm1ftdEx#7lNuJr0n#UMmvqAJ3DvmaDl2>Fv?nDN{QEICIib*HCv)>G#a^V`o74X6$yxKRs zz-3?o(-<^88D6DOCJkd;@&%vGDy3G#Wf$D$8W^Ufw)fB7|2%Wo)0^q}?smKm`(K5@ zT8^ktRi0i;p2UaYoPJB8ylyf@W^j+pI(~n#_ju3413>yIQHB{7yW~Z98LKbuqskX0 zF3TKGdReU$W2L9bsFC?0P`1Qv!Lkann9z8&{Btl*zOemek}N)1k&f$0+lR&8qNu#^ z9a4I3=A_N?oY(^xT@+xnDFESC*B_7Y4|Nxt31#;Sw`58CI~kHqH8TU)LR#xvY!;)r z9)H}NE~{u+%U5-BIRM#%0{8rrJBGoG(`CnDR8_jOw^z!k#~At-mYfHNx%Qfd2Vw2i zJBPtgbM3_^pk6gfwpGEFUYeGLDCBED5>LEq-7(ww-oIy7 zRY_r$VVQZ-{TvVXkh40tWx=+8I%iPs*|8l4@0g^Vn80y)BO)&!PDtO#;#;G+86DA3 z6~|uBfO&ebuw%rh-0RA5JE*PEF=f6pI0#BdO_t}$$CIIcNG47=xTPI0I)$KXnVEy& zE7?1U03P1I)NjMTLmZJ+E#hm{>loc3L+ka8FmCRPX!U(AH}=x^Jox^HxX&ZAK8W{W zcDWceg@cAxzZPtH01xd!gTUTxBi*>pUpGEcj&FJx3E-?Qpik<@J%Em}h0e;#L8HnLKXID? z!2aQNS^xrz<@ytW-j-#;EB(HdWTqIho5R6C$ftaba{@xIZoKjsMKQ?Ar-V&?1qUDc z99Y>%GsAXt;s(h8FMTBkTAwmSChv~~pK~scr}c;(oXZBAI0+U=S@j)!^8GrUZKq6L ztCb^q`k+Jxt;zt?JNKp7Q2Q)cH{Ug7Z_5#`L~)M59I z8dKij?IxxYBw_ftgCvmJ)Lc>oq7GGE5jRzrObVl|2+81`dp3{)IJPzwi3~qymj@}h z>6$8s`<#y)2e8I>YZ{}%h*=~|CB$GVS88n~o_3tPn8Fmapiqv)c+F+8h9bHsg&6Gc zjl}!lTV%&d_-KVW7P0WnbmZK)(Hm|z5BhOicdC)g(mdS-mkKjM6=%3+8oQ9-JFFIV ze-ck{%eHdh$)viL?+89UfPqzn7{c@g7&a%0T$nD2Wug%8%)XbQP#dfH>J_g&I2&a~ z?D#~P4JLQnN4XCc-ts=_l^Ks9FcxSi%}UDck@g>mYRZAK5D+g4swEwIYK9&->v2^OB@FYM z`nPat&4t%B*4<^v4MnOL`4^yPmp&8GypDVxT}=Bvj6y{--uteE!!|bcRBOunx+khk zywPjUre}lO23OquMzd=0NI7fuPOm0$x-E@Hzqo2in+H|uFmMo78hG*>m!(|$2Fnv=b0g5^+*X3YoJcJS~ z)aTGKhD3&+zsP1xGZhbha|v0LFgp#`7zq+SBM;tedCr*pMsufwC`hi6if7QCZ5lJq zsoFRz6e%v8tEK(!;c+FKwyXVo%)edZ>^+i!2PPVSOI|jJG$L}Y6gn~UO}5+!Sb5dm zM~9=HjfAKtiIgfhikTM<;RdCqf!*<8f3$RHu?4R)odb_lU&Fhe={_5Z{ug5nhy9J* zRZwPcyK7_Vq}w4Z1hnuZJD2@YRx^5%I+?fIx%q2sRc-NyC(gaQZ1!V5)#)TiOI6%D z)Q4YB!!@_9qUK2PrmRzGgmHer=}B$8VBY0?n~PpSM}Hshit=zYb=7{*D!w(r#sp0)7T#YJlD=XP@h<}LOj%g}Jh8b?If(jEC|UFe^c2j$EsZuFC7lJsDJ zYBLv}aJEw0dfZI~@W)Iyu~dSLp@c7r)Dfw9X?RW371>G* zht^tIp!$P@&8ge39Fi!@Dje|Pmv+mC-};}yrMsa0m`Upqigggqdwj0w)b0x1eXIiG z2ryYsc1;(=_$HW_-z1wDSDg%hAG^|dFYi&z56fHgX!H&z5Dwr_;Ct>{&5OhHA+yC# zQ`{wLD{(=_Ei^-S&KpCJ8-(0?N?DsuQ1;ePOtin(dAg0|uZcZDSK~PbUE?7-p_MmG zOD1IR)o-{ME8_?+i6n_Pn|_YTa;7-s{$$&?+=DHkoePCS5slU>sH4i3gFmG8qK3aC&ZVQKHCrkD@uSI#ERG(gQvD{Apu`z8-oyeTV z&gZ#MNkJNHcxI%Q{k;adnLFhh{8&4r%)Al3u6nXp2%0{?NJ}!rbyyL>2JJ$!=l2{QHgSSbQ{@ z`ySH~%~~QyshGC(3&;Cm*`}KAkInpd45lpVR!e$v z*aM<^q>4Og@{NrAtbV=S1$epK%iQr5%Hbr1I^T}pq0C>CPl+g!yDdH!@CQZwt|N!W zi8A8P$1QTus-On^uyx(-sn~AmdXKg2O^{aUQ|238exNPPJC15_;TBL7!kFF_kbO8< ziIwc1m+j%g_`0Y=;Pk}22%RKUMy9XA>xnuaig+DusbPYHIe)Q?<{FoB*f*B<(6~%hY*`}(ga@Q<10@O{blTG9*K>g?3+E`4rc}A-6&uiL zVNPP=ms8n-Bp3ZVA~_`o0hx-u$Qh!eDpvR3)7Fa!MBWUEsWubMgI6^h@h9>M2Xr7d zYrg8BDyyKS&u$&HslRqvVC&goCy{)!rM!}NM_g`+iU}9`egwOnLr9j6F5a^0cGBFi z>WEj=Y;}@U%3DiB!<~hwn&WM8uRN0=dP0{%u4uqvP8PM?f+CX)qIh$D#Kau-&8MB! znQ|G#BzK1b=raS2WFD{(R6TzXK4o-FS(}XSY7-}drQ83l8Nxi*$Scp1xQRXB9HqJ? zOAXmM2mEMn{S-zYXIvFD-+a-u5idiR4h44B63p;xp3IZ{s_+~+l4Nwek`9H`D1_e= z1v?n|8az#EE?9XT?d(#Lm-lB_yd58Zp7I15#E9#J2-Jf_=8`+%DQbhl-6w8!6n)8H z#e12U6j@2z(Y*rJ)W7PUh?=RDkE=3u-nXF5zjxtZvG76Gvtv;Rx@Z~I05qn%6s1*W z6tr(oQUz!Vx-l^0-VG`Yrc6&~#eJbOjMuGSpXFsrx+AJ3-dqU*$AUWf^J^Y$~{r)qJcOH2(p%&-Bi(K-ZGp{i&R9& z4_Qkv+wo(5$WHUIhx$lsu8E@M){q5bkj-{Lh^E=D5;sTTmfBe>la@6>x^VN+)Haq> zXh2VFi3zbZE8~X8BXL1faH^+yR?DX*Lq0Qj8PCrv`UB%QSN5CuZl6$yQkf3Q_)=)M zCh3-IbVmcy^U!#Yr<%b$4_1G{n`y4cXMfgXHf`2gm4q$!&->D2LqYW9lRC9_|KIH6fnJ?-L#* zNT-FYG@O;ViYNv*27?f2_wc)u`6y-g$WGJ>ft?R_)N9{uAn*NL^{N}nN)V7PKWM6^ z`-9%9yv4c_hF>h^*+~3k0~B<#JNM#UKuIWcUAFJ1t}dLkXK6)eybGmiN|O~@cQaSn zxp7;S*;)k=3M&R?Al$1Y9_8=a6PAMM7PkJBo04I3)DH%+9F}4)T%s-NTRkOjod> znc)+~>AHR6q$?Tvsop=`z7?6g+M+A_U^OguHM$nZ{H~!oIxkX;I?HN{QvI#%^4<7H7CcJu&i1STcV=T0MT1ozE!mZ zI3#aptu`a=JeE41-=ST3lpsMeBTWhBK;49WsdAG(fN_vQ_EKE}!SkJj3IPT`t(89V zOkII0YxUZyeIc?+!*AZhcMtKNYaNl*vsip3H*dHoThA?EfAVBi=(L41bGH6&T(FS`3d*a~o5@1^nGIU&&2I7HLu7(>6h}8J zbWdHyDn#$TX0%NFSG)gP2kjw$Q@ZXSmZ<1rsY&~}g0aTD6yKmb{Ju6fy_hd0Vj04|Fd6NG}luQ!m2k``$}& z&>}CeBSkW`y{%0SAbNGx({7nCIMvrYmBj65%Edzn6flvf?}ZimR>P~C*pCQ=`!;rz zpGh7J(za5(C)FJeI{#w00Xg)C;U|JCr_UsZ13&&T^@7ai11-HN#^nlQ*wawKPFI}Z zdhNzeKqfeAG_x`X30=|Z)s3}97$4mqH(NIlE~r?})EgjwYpb9p?r|r*niX`i;g#~k zqpSS%wLJrW&Eo)8rku&NQ^&iria{W@s*}2%Tr|_dhraZB=VbhvuCfi%OW0ePluVjH z_2I=@Y5^-26>@m9sJH!orl=j5fi(8Xu}ocMCB#1Ufq1mxl*+Gq<$8hnaYRw}6@hD` z=tEkDTqi%VExom!ZNGgr^P3@n{0_x3rGbY#yr2fs%Ck7T1pRLO>>0zhWbbzF)k=Hl zq;g^#ptHfr5}X5?LB#uRpU>=dS)hla<_vJ>W9s~_gBiE`mD#c`;iGHy`thoy`?P!F zw~r)?9wwto9LkO>8ra>8i9X$9{!YUciCHC@NCSnFOIoEJdd?wCeYmcr1s4I}rgz_5 z5-o6#L((Yj6Q?O`{AQa1#jEe69FNe}|GB}}jgzJ1C!*hwX`_`CIUxXOWebSpG!dD^ z?Z3Y7pBTcZF9TD5=Em`gG({h5R@b2N6Fd>I4=*t=p!W9|WmJ)^^nb||b5YpChh0U; z_fBFDYh~v&3xP-#>vspYfhAwcyAU5QWtxrw!!MGto3gtwBvmB+ra@8y-{z~nYwKjN zygA-X=1)QDO{OwdaLwXwL_pj5?WOE(j`%2O(2RyA#qlF{ZpbH$V|L*yNg0?pLO_8u z?E5F;i-k6aGyH0F^4jjkCf_d)3Q_5&Igv!r34dDga64KvzrE~}cYI~+`1t?3TzI-cVyw*X*0 zI?G`eSs(ZyvBUj;j~cQt{JW?jGsAz!37H5O89AB%6)OCnQ9~vMW`_Sbxlj>GNl_K6 zf&rc->^vktM7cu6$Xr|pDWGNu<2Mj=h(2->Wponu4T(EF2|WpAKo1Zx0a`c=!LE~* zy}yclHyFC8eQpAD5hSJv8mk3rU>F0$A?a2C^Wj)}`mMv|r|+ik=l(0MNw8QGbe|Cq zO5C^o5q-i;0+t^3gr7j)wyO~XNA2uF`Y18I3{5$Nh4P`vp(Y;@eSU-})Q#ikNASl1 zh0il}Z39QOlACsr2LRCpxU~Vt{gB;?CfhYY903WJ04c$M_sJ71>_g5=Ab0cpCX>v( z<>oeKntH@f5c~s22nsPwL4l%Bw@(f;gn$u8p~apsy)7Rl|H`PA6$=g~u~aROse3R5 ziiLv|jHoDJ#E1fXAWi^-D%dc;8Uvm{A*@ss_j_L8O!f~*v!Aeo1UY{wEd*f`{d{nR zzXki{LyJcQ2#EqqF)ZM8>DzDo!ZySP-3|PbjFxfWz@_MfQ9^-Y zjDfc|d(yt6j?*Q>UTG3>xY(NmvTmsAKDA$J<$7Y@-QuvoU!(QSkasg^fff9suQ>2v znt;fFi~94dc0psgz*smrMVYmQEW!emfS|O&rTR%@J`D4Uri>sd^05jE+AcAC4edDi zBMebOf}uFN-Ow}GgockYGCPqL0}>%=`QZ5@IvWe;>w5Y~2xr@WriqS4g~Do+X?)D9 z4clgPbmXDsHQhXS%|}X}CL9ic;@uHCd!^f<;E(R|+At6jcE3$K9Ru1WINhZ{Wej0< zUl8%% zv_D4>fEn$a^o4RTR7d&lz!fRlK|zQF9vZxxeuP-&aNra&92ZGRxrmqg^051LiC5c8 zH+i}GQ6&~d`A=zV&45a7PrZ`-(ujuzP7{eikp_llO@t=Ni_^y$_QKPHz?P1%*JrF( zkgyFlceqQly(mx)l!IXuJ(Whf9no@RxgFbz4#31(ya2CAQ`XkC=`@y7?Yi587-f2E z({0~6pJ=7!*j$T$JlU8VT&r1B4jr|Q>Em11)=WP%pY=E!6PEBIceLJ}3lfS}CbY1b zjB`9x>AQxs4;Zhtm`0WU{6ls)nIh4?PMaSM(tfC{+D4NEsI}v{4rE;}BjlUxap?xx zFsdEH5-^tzy8JGnc)PsST+$%|CktZ=jo(hy;gm7jsq1=3;0UrI&(o_{r9g8`&-#SI zTjt!TiL?lriqm2NY84EV&h#>ARKVR|dP0l_{*kKWezDD`F;6QbeJqi2UG6Dr#ry!@ zN7q1wkcQeYH4vOMmvD6CviisRxrH?d!gf7b#kWj9adi*$+~0mR&emGP^jyrO2Fdo> zzGNrd6=#Exc2$eg>)n&r@o7f&cYhgAhw;<*-&5S#8&Qo-&eAqwyYm9WW3-yQnYL(v zxyDiwSFxG%vGr@)8llYV)sYSx7SlO0@cIhAg3)4cdB^PCx&4K3RTxuTdz1y;&3Fk; z-^E%M+M|i}O@o~=N{5%{Ze?R`jwJk$L0D*T}fJhfK7o7Bx0uJ_Et+NTWV4$44%@!er8oeO1> zFzaL|eYbb|fE4N>Z&t;Fn>0Yp-NS~96N8zd`OV|Xc_I6CT;86JPo&f8kj49S#^5~W zL>swckK}m!FPUGyT;y(>$KjScCu0u2(}U{y520FtoTCAql({-%A*(Lsdk|;cMlKII znTtOa{RN%6Ti8xhA0x{OO}Tv|8Z0K+O)rPxzFvS=HsP#Z2BdV5D=&X`d?wKgC&GYS zZH$6;8)pwZtD(o4zY61DdM??kMJ3~&6fGWUT(GN1aeW#~$V>%pU<20tLLGI&%E;dx z7ne8aF_K96d~tM|Kt}1RdbOxb57Lk4rR<{1pSado&lgq-xe6}9ybbhK?|(2)>zvZs zgf=20UQ%~@DiMn>svenM5M?HBlNqEx}vGdpepo)0sp+bYLRFo7}SB|2$7;y{$Z{I2B@2 zI5ITSm9y;VG+0wsiU$5o${+pODYJiT{SV+1uFy$dJ8ZqNedu69)UkZmmRro@J=5}j zI{&4W3)C~WE);vWBJEvqbJGU;A;!9w{IdNm+|bEVw?Xr6RP*!7dh%NGS7BCbIGNAg z`W;K?ET>d${z`B3X8d>);8u0;q4NUXVpG?^*15`CN9FN=j@QG(TvVyDyf$!}aJXqg zk4p2;VSmffjE+X4B%itFvI2ur${J7_X@!Zsw93I?pm?v|qATE}7RMWLy`I z$wtgO@6?Uf4M?&^$iBXPav#7OTW?Gn7gj!>km?be&>3%Tcgj!6+iWb{x}dWh9i3)* zM=+#PjI2E*f0(J?AGkXMXq%oET;S8VuO|dc8jCJ0l2~t?4k%V3wE8l4lI`0+iVY24 zDB!=eq}au5o_y#@=uR+e8a0X>puMxKH`}C9ify!wBO zR6|QCzOyGyg{}MKk^LSye1cZ4gZ|A_r#3)g7j6w4pHtAYhSt-w69SC}F{-APPN7=G zp+=T0rNFWS!i(Zj)bEU5YN-@6CML3fY0YkSxD7JOfh9R?s3+k}u|=p4$3vVBDpm&y z5e3x>=$){kb^1Z5Ez9(M@JTQI+6_psei#}O=d!b0Ed$iCYX+5Jg>6BJU4p}&jZ zyv`vUE`s+Fe1r8ynmN8uu%8LNRq8numB_+IeOW6?nT!ncvQav7%tgw=&iArGF%}Vs z#%kmCIk?TcS>-9-VUtbWJdg;Vbf39#%yc@s%`*kaa^1MkR8PjscxB$L9ZJx9W!8P| zSUHp^{5o{&KGQrIPv@1XaB-UeomHsx!joTnn$fJ=@tVJD-6QGY{V*cAxQsr{+jx=F zi>%6i@IsWef~?XU?Kv{|2_KOezD4h0Q)2LMK=b8m6WAJ9K=bhY%ljr^;bi&GK@%@+ zo6CUs>+%iyl2E}TJ2~1*88Wm_U_^hlBvhrH?14lC8XPl*?bC-BPkde(ULe|Z1G5t} zez@|nJ*)KPo6p-vz*O2$goa_K~2Bg>_W@J7@Y4c@S=^upd?}4&Ougv`O&}sRdZFKe)e!h zQW4E-$!Ikt?e=L8%vRK&@;*q~SEIJrdy_?#A3MTReE=ohT8#Viv>F7BC;hsDa}7%` zmOFS)ashhffTA1c9;_|lHX>gn);)XPw>f9W2`|q;Ki|L{pJ&3f!eFf5wfv0>i#}3& zao;6F9Ud3maW7;?b1!*^dM_}Ucrrx$JpB&rCSZs2CXPc`=euo;;|_+q^x4GDnai$CGUDtmgJR=r<{x;D7VAe$VdG@`&z~Npzh1eB?f-h^;*Y`Ue2c zKGr`n5~LMj0L!{3yK;@l@Hd{kIJn1yjtBe?s748QgAa!-mRZSed+FSI>niGtdJ&E+ z?LYKq<=1j-NriHcOhc502f)+Yk2|shQTRmXV5@!`X_u z_lozb>9|`|O|*Wh8teOWDTk~YzKFFs`)3@4{M|%Q*d112Q4ap>)DYenApWo`N-~Hl zplqli4xl;2W(X9&o_+;40%t$r7A-^zCLr;oUqB@}Q!4)5%^(ay9#Yy*ussGtXhZn6*7i*K@v;U2%D1-^FJ7@!vd$lxtTmSXDnoYVa)+c1${60aH8 zR8k%jch+{ncIpR1Wt~apv+R9UkFXW#Nf|x0o-o`5aR~hi=Oii2#EX-ZrhXq}{%rkG z3_s~VS?L0@xOUpjYx(_chUqC^`ATv!K$i>N=a}Rve|9Trp_@JfqB8>4Cq?9n8kxQG zX!H*~)Th;1&gu*Wz4zq<6fLy|vKpeZrS0q4Tt{R@>VZR+Vkkoew_5?+P3=f$*>>*t@b3@iG=*=a$r-D!({Do;DG^|Q&=ftq7k&<%>qEb>v3c9_bNGb&HO`H$7h z6$ETrPTgu_`7)wamCXr~N6OdUA!*#C$!j?-9nLNS9b~;g&~GK7OQ29dd5~EK4nQzK zseTcl{BW_J{J<1iG(ivw5aM~GU^!%__K8>Ui=ICPL_gW{q z>rzhe>JaF7ravi(DZFDox(0PH-YX#04hBX0llOfGar@C881O5>Ljug&{DQNoG~wh5 z=jqLY+-|%bZ><)so1XTkt_E1=9nDkz?pBh_7aHJv%yHjl+^le=y-U{;2xg=a;5j-~ zS5;R5zBi9e8K?R4b75Y`2mQ(QeEM^e*LH3Z1O#_g9jsOO$_14Mtj$c_CC|-7P422e zt9aWL4>yNBX-z!(L|k>g_=xwCRek~dhbQw-nUR9x#{k)fq!*64^T40)>Ih05d1*ye3VsD=;n)o}KwwjvK zv*$W$>rTr~tGa(U;9*vp>+SJtg=YxF?9N(%!K9h%5n%A^R3Ks;0y*hEJaHfbktdXR zP&J^bZ$>yk79ne((EQs5lsqw9{6yMxkgS+Md6Ggv^KxdTU0Yhf$%MQGbY9^4j7y9y zLnQuD-y-Iikr-H7jzlnI4C#n_>_mZ+zl^|4Bvdb_GyG7m9XKl!-mg&-AZP3eYj3Sf zD?r76jDrRWl6e|PdWJ><06U7Otz+NJBdOL)2j>Ej_a=^W-lZnu#cNBbcKS;)M-io= z&_djzJ$T2BiTIJ~z@VL`@rG)cV>EimA|5$lK32(h;Jspsy+G%O!9&v(<5eJsAdqZv zlH62`0)25d?Km&Q|Jxj7){Er)8j&PIw8Hym|sbrZY53p)D;Eq93y&4 z3{BI$8~vk;^s(F|+*m5F^Rj`9!lmj&W=C+cutB#{%uN^$aDpDZ4KHjYdO;{54yyP! zIAH&~Ibdhz{Ljt7-8^XWeV>ZrZX+>(=8W zUVCTb=`D`BDH0Ekj<6ROHx*p?QrtY)ifJC{yu4p-QR#xoms>=oTR(oKH(c7Vl4Z$U z(?M>2J6)F<;Jx(T`pjDl-@X-p-n@(;{g*jl@7-|dNHb{|u(24UeKvkkoUNCY;&GU; z&M`YJIEOxAI#IxezWtaJAL-rEp^2046EfJI0AsEvE3uXp`&%aYt%vhFO9T63J1 z!2|bRsQ82^rTemSC;i%eVVkg=a{uhMhFzq`qeiE0?b48Q|Ix;Y!S+6}m3Z^v&LY)P zO`xo-nl@{%nBp_*lT5JhH9S+J1zu48nI>u0@4~KtEV=LK6-1>0UJOh7n-O}=K(ouJ zsFZfp(Lz{2DU4+q^r4 zi*iMj&i=i-lR!6rH@;+E z6oex_gyTDhyns&*A=NOt$3XhiRLq_d1RtSK6aebj$+ei{xX+zx+c7t9G7{R`W^U|I zR0P0I_0TLIIvJ?|U=OSLlP#wG#X~P{^TNJ(=*>TQXyq3Vwfx_B=wB4|(S0W9i-&q= zGJWw-sH{S%7arsP!9%${yzfgKH&-F1c^m%Xp^J3~FL?N^#DDS7p}|hL2%V8F`s>>w zxqkzhzvq>Gi8ud@{QQ3sZ(#oknbi`d(yvO#EwKclf*g07m7AVWbgOv$yNB1e5xmK* z`1GqO6)Wu-LDA{$D~_JWlH8KM7cKt~Z-&`==`E=>rTlHNEcCIkL2DhOg?&nCZ{kzTb7nzkCxC!JS$ z7_WyoU4z<+XCA(RPjq{sbRWc*39T!=7X_XRs1lHro;YBLJ!xmINp$U;54u+sN;{hU zB+hZfzdLg1>di=K)Ij5`R;9k(4gKQhF?p0hGGS?KCru~DWS@WRu1iWkot6#L|@nUqKsf_ zwSk=K1Rn4&$UwV^F9GQwP5}QU-QXY(0Cp#u(7vRbYVwp2!>2Be`(4xsW z^qmH99)CPxNK?YDpWG5)YDnmJPUwblnJtU@`t{p8{l|^wZEjBDSA?ELvXe4;G2Z^b6rryrjY3|6GMKP zW)VQ?62j?_9P%GMSkM(dWij{-Q2G79DivyzcvqsLn_^wzAz(obZdkLkR`iu zXZq4^lK$G1i*09nyj9?acaBZVyT_-CG!wZmZ)fGp*s`;VYWMJn&ing08UtEA{o zC%3dM)OLW!e*WdUXNW8TQb|@NU?3)C~CnoF`Q1^KUT2`FAtJ z$?$)NnN=>two}!rV4&k}Rz4C1I#|#=BT&veQ0c3yim1Lls1ReJ288Fw=vvD(>sX

Ez+|IZFg{=w9LW z-S0Sn7CID)Z$JRsJ96P}6@Dl5y~vLgzxyTjSKI*(v|AWd#?ClT%mJU4tM8L~suZ*@ zIo#SEM~$nbyj9QF?`~nO!PmyRG)V>@)U{dd6OVQ!OOg*hJJv-QKkWuICuUy>zXo80xjVlC@DNt%WLB?P(IhEV4VstB!Z;p$ z%gqNhGj%JPoktR-U#30e%fs?#C{=#>t$m#9F)2xE>xtM{i^v4U1kV*(}JZ&W6#iK!xe!a(IV zus)h!wR5ZDS5?`1)Zi9CQ@wORO{_o=G85JOM>|zVQ!0|Q$Xmc}S)_N#Ys?6#q}c*h z>a%_xITNYibC3W_EIHz`+5+HvtQ1yj2Gv~K)PmlQOpEf0!nOCt47y)a8aYu)Bx*qP z#gsOb2gIr{#NWG#Y3PcC^YF3)ZTB>_HvK5s6*Ls?Eul z1(+k!b3n9M4(s^!XxM+1nQ)9~PvuN7%1$j8e1?gl+;o(Bro9AXZ$EJJPoDhYpO(`x zqgSb)IowWZl0QuHt3D2xF{KIGXgbWWL`VZu=$R+?8+$MFj}g;bk1?C44{?m4s&mxs zF?E06;Z%%0rSwYups>rz0UA`@R;+Z|tRzs`Qx0(UQE%LnRY}~ry0s)Jt7cot>(m?S z5!3P?V>@&mAP{02?x43V;VM_rhUF|AO>vjs+V+-Iou0MFx#WEGqPHV(uB(4LXhhl5 z*^D!455f&@DSjSc(ig*T$|b&EIzM!>FMoE7eLN~ZbT%%3ZWw>Qa=r2moqcB4eEyvH zKp=bNOIZFSH2MrCd+lU&|6R~FF8-yTv3@k~Qns{Xue+zt#0y@cwxhTa|6<(8tn2;o z)4mh4-M+u>UyPQ%=BrRU!TUR=S^qv#f`EhhpLx^&u^s$xMJ@b1kV6S_cp`YN*T+n= z`~f<2HBt!)Woyc}(}y{*s{|tEo#*z=(Lh=GbjWqw*RuK`;r!IOOle#8ltmMQA2bAt0f`?(ulqZxABrC#L2d|a z*G(8Pgd9*+3Y1?DU6?E7_nfYrklq$BAiE+`P)Q&(KSOwsA3&&n#$fPp3gdYxU_kk) z2`g*BLliGmAbBf{Kv5W61{L|j!Vw9{3Uf&O0dQzT@k|l9%3(VSR5Iw4K#-CW$rCwq ze!UFALV;kx^NK%dUM%B_htW!`xmcG4`~7h<(o@g)PvVM);VCC%q+N>E z_7&rCN9^|J_HceUdYe1gO_?PIVU{a?uETD}q%5H@MIBE)mq~F>q%5fHw(WTDyH}{0 zr$*h@_K$IpXQX;mG_6C8%%vi@*+7lR(Rx|QP`1Goxr06O@Y)_JNE7#xBlML6sQs!_vBSqY8XiG;WVz^#9X0#qmFX?_ ziL;Fz4hIXWbglP^rDVGM4mZkpvK}{kf8uR`IPoeeH$Fy%w6lH@wr$_}XBT2dJFGsB zvKwt)#)|O0JVoeokNGG~^I2n0qx~B|{?DkQoP({Av!SsA0fn8BiITAzEh8NZ1053+ zB{aQ|gR#Do?bng7-T&G&x3v-0cQPiR5awcHU}9omWoBe#WngAypk-hoXJ8=zYA0=L z^uIMxa?rQ4Gd3ci7uC0NG=`>^QxZ|56Lq$-GSIiN`KvEVX6B9rU!Q;UiH1PM*unAZ z)CA0Qj4X_-oJ?%2%(P5&Z2#))e|3>S%G|}6;ICfM8yUOM3mY4l>)TKgSla13nK>Hk z6S(~IKiKG)=~(`1LJ183{B?EyiWM>b&!CZ#yPYusy{v(SqLVc=y$k^}<6ogx2S+CY zCMM4RjwvxQu>3PD@V^IQ%hfFvH93*qFLZOFTUmqdeU$SHsQA@t0+(TwE=7QXfo07-C50vQ<9y8uX zK1dYDL)I`mWHe;hk)Bl3_*2cFER{h~CWDEqpW$)@i46-VuNubEZ?}W?T%t&KmK*7D zUvk!H!cd%qN-{*zfT}c?+?HklK}^2^`D}wBg2acc17=8$=s7eCTp%(#3;6`T^2GCw zSG=4r&E$D~nS*@xv0J>`+jz2H;F&1s33!)~&19LqzTGQG=|SbOTU&#NXx*jX`Jqps zmzF&4j4(>BCiWn4ewej<16&rWS{cYw_2UYb)xXL*Kq#3VG`m3>#>6^)DY2lTueC0N zT=3kDtYU7U;O5=P%W-9_4;AtpAJTonD<#yj#}TaYu3gW4y|sq&o*A&6(ea&TI-Kiq zn^BS}Z*7vulY=`AAVMOyMIDRHC%nLhJhkW2W&=_Ff@Pj80S;3J zWtBd(Fb!EH@;&k=f&$fUN3>iKzhx}tGV*4!RTfiO5k+tg%2)A%Say=dyzD9#DmX*J z79S!pDpG?8ZxT|C4Jr-LLy@d()cRcKDt~!?Mj{RvHOA?OVGD5wyKEA&kf@?x!moQn zeF>!8-Fm@W6z3Q}O?joAydAo5kjx1l88dMsNrWJIR{ty$a6C~fQ2Uu|J7uGEr4pp) z#Bh(IMFyav51G?_SHO|jHuU@D9Jp-0eggJ`^dU&U`A)L&o;FVZ`QdOi$@`}~8h!M{ zqh^%S$S5G}hi>DIh%&NwDdte62-$;stW!bzp^}k0OL7UL(^gm+|L=8!f+BGIr9*!Q z;+g5YU#zV8t&{1UIJlEcDBMM%9M2>bCNeh;Z z3VhVP{$MFICu^)Qi8bRj%9iu|-l$%grCaUatY)Lth2+ax#HYfoZUastsDQVo`DuRK zLuQ8b3;PrDVh>$+?1aXIjcuE78%HcQEYRkQ*t=P7qPtocqpiL^>#HH1H&CHhRkL9! zw}#`U3*pno`ww>tE@$_Su z0t$t0C;AXU+aU}AR%}(RgQ|Fj+z9yszT*hggz=9!bKN~%yoCwr8T{@ulg9(XfU|-9JNTsS9SH^!`q(I=sGmD zhRuV^*4ygoq5GFG?`QBOdHFQ?_$i(#XMk{+~ zz~u2@q%&(cE?sG=1J?^F@zks5V8Y2YxJ(u}@-u^%SLR6Y;k4Bgiw4|+3hGhs@^xiV z@nB)=VA_y7j+!ZocW&@uPI1q^Ai0*igw1Y-{4I3D!F@x)HX91wlKThCE*)GCW9v`z zUfs0vb}$vSS=F#N@3n3_ri|%|0IL30fmH$ClVOD3)#Y~-rSOyGr#E{y<_<8NTHf@F zpf!KwRy(Jn?LIIp+{LM-^JypBN%+(?5oD)GQDwWEuSo=W(%rhH1O= zSQ&gg87v=FAfnzCIo@#6O{3I%ksAhQn{Bzw7H*Q`-uHPHEX{T{1r*W068Yus!;u{I`xI<#od{SNzf=;1znYua60rW|&Hc3@ zOrXglD8eQrBq+!s!Yaxsz%D4nCd~P@#>BuNB+4wrz|IID_}`m+okhyn#?;A-fRX(_ zA8Lj>$_+1|j}f?yIhe8Q(}_kDD^3!p7&9f#f)u2X--wtvgQgY?1{nmV7AX2>E34A| zMurjS3I-Se25FQ-6`EQUUQiTLLJ|JkEatYDcZG@2HI5(V=2V+FI8Y`ku@a z$rfy>vz9K)nEJL26G Date: Fri, 24 Jun 2022 08:59:31 -0700 Subject: [PATCH 04/20] incorrect mekf, incorrect noise adding --- Simulator/kf.jl | 38 ++++++++++++++++++++++++++++---------- Simulator/simulator.jl | 17 ++++++++++++----- Simulator/tests/tests.jl | 11 +++++++++-- 3 files changed, 49 insertions(+), 17 deletions(-) diff --git a/Simulator/kf.jl b/Simulator/kf.jl index 0e0f94e..be00fdc 100644 --- a/Simulator/kf.jl +++ b/Simulator/kf.jl @@ -22,8 +22,13 @@ function f( δt::Float64 ) θ = norm(ω - β) * δt - r = (ω - β) / norm(ω - e.β) - return [L(q) * [r .* sin(θ / 2); cos(θ / 2)], β] + r = (ω - β) / norm(ω - β) + return L(q) * [r .* sin(θ / 2); cos(θ / 2)] +end + +function quaternionToMatrix(q::Vector{Float64}) + s, v = q[1], q[2:4] + return I(3) + 2 * hat(v) * (s * I(3) + hat(v)) end function step( @@ -36,22 +41,35 @@ function step( ᵇr_sun::Vector{Float64} ) # Predict: - xₚ = f(e.q, e.β, ω, δt) + qₚ = f(e.q, e.β, ω, δt) # β remains constant R = exp(-hat(ω - e.β) * δt) A = [ R (-δt*I(3)) zeros(3, 3) I(3) ] - W = zeros(6, 6) # related to some noise or something (ask Zac) - Pₚ = A * P * A' + W + W = I(6) * 0.01 # related to some noise or something (ask Zac) + Pₚ = A * e.P * A' + W # Innovation + Q = quaternionToMatrix(qₚ) Z = [ᵇr_mag; ᵇr_sun] - [Q zeros(3, 3); zeros(3, 3) Q] * [ⁿr_mag; ⁿr_sun] - C = [hat(ᵇr_mag) zeroes(3, 3); hat(ᵇr_sun) zeros(3, 3)] - V = zeros(3, 3) # Something else + C = [hat(ᵇr_mag) zeros(3, 3); hat(ᵇr_sun) zeros(3, 3)] + V = I(6) * 0.01 # Something else S = C * Pₚ * C' + V # Kalman Gain - L = Pₚ * C' * inv(S) + print(S) + Lk = Pₚ * C' * inv(S) # Update - δx = L * Z - return e.q + δx = Lk * Z + ϕ = δx[1:3] + δβ = δx[4:6] + θ = norm(ϕ) + r = ϕ / θ + qᵤ = ⊙(qₚ, [r * sin(θ / 2); cos(θ / 2)]) + βᵤ = e.β + δβ + Pᵤ = (I(6) - Lk * C) * Pₚ * (I(6) - Lk * C)' + Lk * V * Lk' + + e.q = qᵤ + e.β = βᵤ + e.P = Pᵤ + return e end \ No newline at end of file diff --git a/Simulator/simulator.jl b/Simulator/simulator.jl index 12c48d3..b3dadea 100644 --- a/Simulator/simulator.jl +++ b/Simulator/simulator.jl @@ -258,7 +258,7 @@ module Simulator end - function my_sim_kf(control_fn) + function my_sim_kf(control_fn, N) x₀ = initialize_orbit() println("intialized orbit!") # x₀[11:13] .=0 @@ -268,24 +268,31 @@ module Simulator t = Epoch(2020, 11, 30) # Starting time is Nov 30, 2020 dt = 0.5 # Time step, in seconds - N = 1000000 q_hist = zeros(N, 9) q_hist[1, 1] = norm(x₀[11:13]) x = x₀ q₀ = x₀[7:10] - P = Diagonal([1, 1, 1, 1]) - kf = EKF(q₀, P) + β₀ = [0; 0; 0] + P₀ = I(6) + kf = EKF(q₀, β₀, P₀) for i = 1:N - 1 r, v, q, ω = x[1:3], x[4:6], x[7:10], x[11:13] + b = IGRF13(r, t) + r_sun = sun_position(t) + to_sun = r_sun - r + x = rk4(x, J, control_fn(ω, b), t, dt) t += dt # Don't forget to update time (not that it really matters...) q_hist[i + 1, 1] = norm(ω) q_hist[i + 1, 2:5] .= q - q_hist[i+1, 6:9] .= step(kf, q, ω, dt) + rsun = (rand(3) * .05 .- 0.025) + rb = (rand(3) * .05 .- 0.025) + step(kf, ω, dt, to_sun, b, to_sun + rsun, b + rb) + q_hist[i+1, 6:9] .= kf.q if norm(ω) < 0.1 q_hist = q_hist[1:i, :] diff --git a/Simulator/tests/tests.jl b/Simulator/tests/tests.jl index 7e620f3..660710e 100644 --- a/Simulator/tests/tests.jl +++ b/Simulator/tests/tests.jl @@ -115,8 +115,15 @@ end return cross(m, b) end - data = my_sim_kf(control_law) - display(plot(data, title="MEKF/DeTumbling", xlabel="Time (s)", ylabel="Angular Velocity (rad/s)", labels=["ω" "i1" "i2" "i3" "r" "i1'" "i2'" "i3'" "r'"])) + data = my_sim_kf(control_law, 5) + display(plot( + data, + title="MEKF/DeTumbling", + xlabel="Time (s)", + ylabel="Angular Velocity (rad/s)", + labels=["ω" "s" "v1" "v2" "v3" "s'" "v1'" "v2'" "v3'"], + linecolor=[:red :blue :green :purple :orange :blue :green :purple :orange] + )) end From 0420b848df98c4e09cfa0179c5649f2f5efa09fc Mon Sep 17 00:00:00 2001 From: thetazero Date: Fri, 24 Jun 2022 09:50:18 -0700 Subject: [PATCH 05/20] Still wrong --- Simulator/kf.jl | 22 +++++++++++----------- Simulator/quaternions.jl | 12 ++++++++++++ Simulator/simulator.jl | 21 +++++++++++++++------ Simulator/tests/tests.jl | 2 +- 4 files changed, 39 insertions(+), 18 deletions(-) diff --git a/Simulator/kf.jl b/Simulator/kf.jl index be00fdc..35401e2 100644 --- a/Simulator/kf.jl +++ b/Simulator/kf.jl @@ -23,12 +23,7 @@ function f( ) θ = norm(ω - β) * δt r = (ω - β) / norm(ω - β) - return L(q) * [r .* sin(θ / 2); cos(θ / 2)] -end - -function quaternionToMatrix(q::Vector{Float64}) - s, v = q[1], q[2:4] - return I(3) + 2 * hat(v) * (s * I(3) + hat(v)) + return L(q) * [cos(θ / 2); r * sin(θ / 2)] end function step( @@ -49,24 +44,29 @@ function step( ] W = I(6) * 0.01 # related to some noise or something (ask Zac) Pₚ = A * e.P * A' + W + # Innovation Q = quaternionToMatrix(qₚ) Z = [ᵇr_mag; ᵇr_sun] - [Q zeros(3, 3); zeros(3, 3) Q] * [ⁿr_mag; ⁿr_sun] + println(Z) + println("q's should be about the same", qₚ, e.q) + println("should be about equal for sun vector", Q * ⁿr_sun, ᵇr_sun) C = [hat(ᵇr_mag) zeros(3, 3); hat(ᵇr_sun) zeros(3, 3)] V = I(6) * 0.01 # Something else S = C * Pₚ * C' + V + # Kalman Gain - print(S) - Lk = Pₚ * C' * inv(S) + L = Pₚ * C' * inv(S) + # Update - δx = Lk * Z + δx = L * Z ϕ = δx[1:3] δβ = δx[4:6] θ = norm(ϕ) r = ϕ / θ - qᵤ = ⊙(qₚ, [r * sin(θ / 2); cos(θ / 2)]) + qᵤ = ⊙(qₚ, [cos(θ / 2); r * sin(θ / 2)]) βᵤ = e.β + δβ - Pᵤ = (I(6) - Lk * C) * Pₚ * (I(6) - Lk * C)' + Lk * V * Lk' + Pᵤ = (I(6) - L * C) * Pₚ * (I(6) - L * C)' + L * V * L' e.q = qᵤ e.β = βᵤ diff --git a/Simulator/quaternions.jl b/Simulator/quaternions.jl index e12ea3b..362d1ee 100644 --- a/Simulator/quaternions.jl +++ b/Simulator/quaternions.jl @@ -89,3 +89,15 @@ const T = [1 0 0 0; 0 -1 0 0; 0 0 -1 0; 0 0 0 -1] # Forms the conjugate of q, ⊙(q₁, q₂) = L(q₂) * q₁ # Hamiltonian product +""" quaternionToMatrix (q) + Arguments: + - q: Scalar-first unit quaternion | [4,] + + Returns: + - Q: Rotation matrix representing the same rotation +""" +function quaternionToMatrix(q::Vector{Float64}) + s, v = q[1], q[2:4] + return I(3) + 2 * hat(v) * (s * I(3) + hat(v)) +end + diff --git a/Simulator/simulator.jl b/Simulator/simulator.jl index b3dadea..3fc2ee2 100644 --- a/Simulator/simulator.jl +++ b/Simulator/simulator.jl @@ -266,32 +266,41 @@ module Simulator J = [0.3 0 0; 0 0.3 0; 0 0 0.3] # Arbitrary inertia matrix for the Satellite t = Epoch(2020, 11, 30) # Starting time is Nov 30, 2020 - dt = 0.5 # Time step, in seconds + dt = 0.01 # Time step, in seconds q_hist = zeros(N, 9) q_hist[1, 1] = norm(x₀[11:13]) + q_hist[1, 2:5] .= x₀[7:10] + q_hist[1, 6:9] .= x₀[7:10] x = x₀ q₀ = x₀[7:10] β₀ = [0; 0; 0] P₀ = I(6) + println("q0:", q₀) kf = EKF(q₀, β₀, P₀) for i = 1:N - 1 r, v, q, ω = x[1:3], x[4:6], x[7:10], x[11:13] b = IGRF13(r, t) - r_sun = sun_position(t) - to_sun = r_sun - r x = rk4(x, J, control_fn(ω, b), t, dt) t += dt # Don't forget to update time (not that it really matters...) q_hist[i + 1, 1] = norm(ω) q_hist[i + 1, 2:5] .= q - rsun = (rand(3) * .05 .- 0.025) - rb = (rand(3) * .05 .- 0.025) - step(kf, ω, dt, to_sun, b, to_sun + rsun, b + rb) + r_sun = sun_position(t) + inertial_sun = normalize(r_sun - r) + inertial_mag = normalize(b) + + rsun = (rand(3) * .01 .- 0.005) + rmag = (rand(3) * .01 .- 0.005) + Q = quaternionToMatrix(q) + body_sun = Q * (normalize(inertial_sun + rsun)) + body_mag = Q * (normalize(inertial_mag + rmag)) + println("body_sun: ", body_sun, "body_sun(predicted): ", Q * inertial_sun) + step(kf, ω, dt, inertial_mag, inertial_sun, body_mag, body_sun) q_hist[i+1, 6:9] .= kf.q if norm(ω) < 0.1 diff --git a/Simulator/tests/tests.jl b/Simulator/tests/tests.jl index 660710e..ef434e8 100644 --- a/Simulator/tests/tests.jl +++ b/Simulator/tests/tests.jl @@ -115,7 +115,7 @@ end return cross(m, b) end - data = my_sim_kf(control_law, 5) + data = my_sim_kf(control_law, 20) display(plot( data, title="MEKF/DeTumbling", From 255da731370ccea8a6f46286107621437f079565 Mon Sep 17 00:00:00 2001 From: thetazero Date: Fri, 24 Jun 2022 10:40:23 -0700 Subject: [PATCH 06/20] Upload --- Simulator/kf.jl | 7 ++++--- Simulator/simulator.jl | 9 +++++---- Simulator/tests/tests.jl | 5 +++-- 3 files changed, 12 insertions(+), 9 deletions(-) diff --git a/Simulator/kf.jl b/Simulator/kf.jl index 35401e2..36857a4 100644 --- a/Simulator/kf.jl +++ b/Simulator/kf.jl @@ -48,9 +48,6 @@ function step( # Innovation Q = quaternionToMatrix(qₚ) Z = [ᵇr_mag; ᵇr_sun] - [Q zeros(3, 3); zeros(3, 3) Q] * [ⁿr_mag; ⁿr_sun] - println(Z) - println("q's should be about the same", qₚ, e.q) - println("should be about equal for sun vector", Q * ⁿr_sun, ᵇr_sun) C = [hat(ᵇr_mag) zeros(3, 3); hat(ᵇr_sun) zeros(3, 3)] V = I(6) * 0.01 # Something else S = C * Pₚ * C' + V @@ -59,7 +56,10 @@ function step( L = Pₚ * C' * inv(S) # Update + println("Z: ", Z) δx = L * Z + println(size(L), size(Z)) + println(δx, size(δx)) ϕ = δx[1:3] δβ = δx[4:6] θ = norm(ϕ) @@ -71,5 +71,6 @@ function step( e.q = qᵤ e.β = βᵤ e.P = Pᵤ + println("β: ", e.β, "δβ: ", δβ) return e end \ No newline at end of file diff --git a/Simulator/simulator.jl b/Simulator/simulator.jl index 3fc2ee2..0192745 100644 --- a/Simulator/simulator.jl +++ b/Simulator/simulator.jl @@ -266,7 +266,7 @@ module Simulator J = [0.3 0 0; 0 0.3 0; 0 0 0.3] # Arbitrary inertia matrix for the Satellite t = Epoch(2020, 11, 30) # Starting time is Nov 30, 2020 - dt = 0.01 # Time step, in seconds + dt = 0.1 # Time step, in seconds q_hist = zeros(N, 9) q_hist[1, 1] = norm(x₀[11:13]) @@ -294,12 +294,13 @@ module Simulator inertial_sun = normalize(r_sun - r) inertial_mag = normalize(b) - rsun = (rand(3) * .01 .- 0.005) - rmag = (rand(3) * .01 .- 0.005) + # rsun = (rand(3) * .01 .- 0.005) + # rmag = (rand(3) * .01 .- 0.005) + rsun = [0;0;0] + rmag = [0;0;0] Q = quaternionToMatrix(q) body_sun = Q * (normalize(inertial_sun + rsun)) body_mag = Q * (normalize(inertial_mag + rmag)) - println("body_sun: ", body_sun, "body_sun(predicted): ", Q * inertial_sun) step(kf, ω, dt, inertial_mag, inertial_sun, body_mag, body_sun) q_hist[i+1, 6:9] .= kf.q diff --git a/Simulator/tests/tests.jl b/Simulator/tests/tests.jl index ef434e8..30af980 100644 --- a/Simulator/tests/tests.jl +++ b/Simulator/tests/tests.jl @@ -115,14 +115,15 @@ end return cross(m, b) end - data = my_sim_kf(control_law, 20) + data = my_sim_kf(control_law, 15) display(plot( data, title="MEKF/DeTumbling", xlabel="Time (s)", ylabel="Angular Velocity (rad/s)", labels=["ω" "s" "v1" "v2" "v3" "s'" "v1'" "v2'" "v3'"], - linecolor=[:red :blue :green :purple :orange :blue :green :purple :orange] + linecolor=[:red :blue :green :purple :orange :blue :green :purple :orange], + linewidth=[1 1 1 1 1 3 3 3 3], )) end From 49f2d0c2fde802cffeb58828516590e6648d708f Mon Sep 17 00:00:00 2001 From: thetazero Date: Mon, 27 Jun 2022 05:45:21 -0700 Subject: [PATCH 07/20] debuging --- Simulator/kf.jl | 10 ++++------ Simulator/simulator.jl | 2 +- Simulator/tests/tests.jl | 2 +- 3 files changed, 6 insertions(+), 8 deletions(-) diff --git a/Simulator/kf.jl b/Simulator/kf.jl index 36857a4..c54e6fc 100644 --- a/Simulator/kf.jl +++ b/Simulator/kf.jl @@ -42,24 +42,22 @@ function step( R (-δt*I(3)) zeros(3, 3) I(3) ] - W = I(6) * 0.01 # related to some noise or something (ask Zac) + W = I(6) * 1e-5 # related to some noise or something (ask Zac) Pₚ = A * e.P * A' + W # Innovation Q = quaternionToMatrix(qₚ) Z = [ᵇr_mag; ᵇr_sun] - [Q zeros(3, 3); zeros(3, 3) Q] * [ⁿr_mag; ⁿr_sun] C = [hat(ᵇr_mag) zeros(3, 3); hat(ᵇr_sun) zeros(3, 3)] - V = I(6) * 0.01 # Something else + V = I(6) * 1e-5 # Something else S = C * Pₚ * C' + V # Kalman Gain L = Pₚ * C' * inv(S) # Update - println("Z: ", Z) + # println("Z: ", Z) δx = L * Z - println(size(L), size(Z)) - println(δx, size(δx)) ϕ = δx[1:3] δβ = δx[4:6] θ = norm(ϕ) @@ -71,6 +69,6 @@ function step( e.q = qᵤ e.β = βᵤ e.P = Pᵤ - println("β: ", e.β, "δβ: ", δβ) + # println("β: ", e.β, "δβ: ", δβ) return e end \ No newline at end of file diff --git a/Simulator/simulator.jl b/Simulator/simulator.jl index 0192745..29457a7 100644 --- a/Simulator/simulator.jl +++ b/Simulator/simulator.jl @@ -266,7 +266,7 @@ module Simulator J = [0.3 0 0; 0 0.3 0; 0 0 0.3] # Arbitrary inertia matrix for the Satellite t = Epoch(2020, 11, 30) # Starting time is Nov 30, 2020 - dt = 0.1 # Time step, in seconds + dt = 0.01 # Time step, in seconds q_hist = zeros(N, 9) q_hist[1, 1] = norm(x₀[11:13]) diff --git a/Simulator/tests/tests.jl b/Simulator/tests/tests.jl index 30af980..535fd73 100644 --- a/Simulator/tests/tests.jl +++ b/Simulator/tests/tests.jl @@ -115,7 +115,7 @@ end return cross(m, b) end - data = my_sim_kf(control_law, 15) + data = my_sim_kf(control_law, 5000) display(plot( data, title="MEKF/DeTumbling", From c07c817f4b17a586d09c825815f0b73aece34464 Mon Sep 17 00:00:00 2001 From: thetazero Date: Mon, 27 Jun 2022 09:45:08 -0400 Subject: [PATCH 08/20] Simple MEKF works, full MEKF almost works --- Simulator/kf.jl | 54 +++++++++++++++++++++++++++++++++------- Simulator/simulator.jl | 12 +++++++-- Simulator/tests/tests.jl | 11 +++++--- 3 files changed, 63 insertions(+), 14 deletions(-) diff --git a/Simulator/kf.jl b/Simulator/kf.jl index c54e6fc..7c6c3ef 100644 --- a/Simulator/kf.jl +++ b/Simulator/kf.jl @@ -35,33 +35,35 @@ function step( ᵇr_mag::Vector{Float64}, ᵇr_sun::Vector{Float64} ) + q = e.q + β = e.β + P = e.P + W = I(6) * 1e-6 # related to some noise or something (ask Zac) + V = I(6) * 1e-6 # some other noise # Predict: - qₚ = f(e.q, e.β, ω, δt) # β remains constant - R = exp(-hat(ω - e.β) * δt) + qₚ = f(q, β, ω, δt) # β remains constant + R = exp(-hat(ω - β) * δt) A = [ R (-δt*I(3)) zeros(3, 3) I(3) ] - W = I(6) * 1e-5 # related to some noise or something (ask Zac) - Pₚ = A * e.P * A' + W + Pₚ = A * P * A' + W # Innovation Q = quaternionToMatrix(qₚ) Z = [ᵇr_mag; ᵇr_sun] - [Q zeros(3, 3); zeros(3, 3) Q] * [ⁿr_mag; ⁿr_sun] C = [hat(ᵇr_mag) zeros(3, 3); hat(ᵇr_sun) zeros(3, 3)] - V = I(6) * 1e-5 # Something else S = C * Pₚ * C' + V # Kalman Gain L = Pₚ * C' * inv(S) # Update - # println("Z: ", Z) δx = L * Z ϕ = δx[1:3] δβ = δx[4:6] - θ = norm(ϕ) - r = ϕ / θ + θ = -norm(ϕ) + r = normalize(ϕ) qᵤ = ⊙(qₚ, [cos(θ / 2); r * sin(θ / 2)]) βᵤ = e.β + δβ Pᵤ = (I(6) - L * C) * Pₚ * (I(6) - L * C)' + L * V * L' @@ -69,6 +71,40 @@ function step( e.q = qᵤ e.β = βᵤ e.P = Pᵤ - # println("β: ", e.β, "δβ: ", δβ) + println("β: ", e.β, "δβ: ", δβ) return e +end + +# for the "simplest MEKF" +function simplest_step( + e::EKF, + ω::Vector{Float64}, + δt::Float64, + ⁿr_mag::Vector{Float64}, + ⁿr_sun::Vector{Float64}, + ᵇr_mag::Vector{Float64}, + ᵇr_sun::Vector{Float64} +) + q = e.q + P = e.P + W = I(3) * 1e-6 + V = I(6) * 1e-6 + # Predict + r = normalize(ω) + θ = norm(ω) * δt + qₚ = ⊙(q, [cos(θ / 2); r * sin(θ / 2)]) + Pₚ = P + W + # Innovation + Q = quaternionToMatrix(qₚ) + Z = [ᵇr_mag; ᵇr_sun] - [Q zeros(3, 3); zeros(3, 3) Q] * [ⁿr_mag; ⁿr_sun] + C = [hat(Q * ⁿr_mag); hat(Q * ⁿr_sun)] + S = C * Pₚ * C' + V + # Kalman Gain + L = Pₚ * C' * inv(S) + # Update + ϕ = L * Z + θ = -norm(ϕ) + r = normalize(ϕ) + e.q = ⊙(qₚ, [cos(θ / 2); r * sin(θ / 2)]) + e.P = (I(3) - L * C) * Pₚ * (I(3) - L * C)' + L * V * L' end \ No newline at end of file diff --git a/Simulator/simulator.jl b/Simulator/simulator.jl index 29457a7..778057e 100644 --- a/Simulator/simulator.jl +++ b/Simulator/simulator.jl @@ -258,7 +258,7 @@ module Simulator end - function my_sim_kf(control_fn, N) + function my_sim_kf(control_fn, N, type) x₀ = initialize_orbit() println("intialized orbit!") # x₀[11:13] .=0 @@ -276,7 +276,11 @@ module Simulator q₀ = x₀[7:10] β₀ = [0; 0; 0] + P₀ = I(6) + if type == "simple" + P₀ = I(3) + end println("q0:", q₀) kf = EKF(q₀, β₀, P₀) @@ -301,7 +305,11 @@ module Simulator Q = quaternionToMatrix(q) body_sun = Q * (normalize(inertial_sun + rsun)) body_mag = Q * (normalize(inertial_mag + rmag)) - step(kf, ω, dt, inertial_mag, inertial_sun, body_mag, body_sun) + if type == "simple" + simplest_step(kf, ω, dt, inertial_mag, inertial_sun, body_mag, body_sun) + else + step(kf, ω, dt, inertial_mag, inertial_sun, body_mag, body_sun) + end q_hist[i+1, 6:9] .= kf.q if norm(ω) < 0.1 diff --git a/Simulator/tests/tests.jl b/Simulator/tests/tests.jl index 535fd73..8d1e0f5 100644 --- a/Simulator/tests/tests.jl +++ b/Simulator/tests/tests.jl @@ -115,15 +115,20 @@ end return cross(m, b) end - data = my_sim_kf(control_law, 5000) + type = "smple" + data = my_sim_kf(control_law, 5000, type) + name = "MEKF/DeTumbling" + if type == "simple" + name = "MEKF(simplest)/DeTumbling" + end display(plot( data, - title="MEKF/DeTumbling", + title=name, xlabel="Time (s)", ylabel="Angular Velocity (rad/s)", labels=["ω" "s" "v1" "v2" "v3" "s'" "v1'" "v2'" "v3'"], linecolor=[:red :blue :green :purple :orange :blue :green :purple :orange], - linewidth=[1 1 1 1 1 3 3 3 3], + linewidth=[1 1 1 1 1 1 1 1 1], )) end From 96f95a156840e32fe98e7d29211b5f5bb8e95d28 Mon Sep 17 00:00:00 2001 From: thetazero Date: Mon, 27 Jun 2022 11:39:59 -0400 Subject: [PATCH 09/20] Working MEKF --- Simulator/kf.jl | 11 +++++------ Simulator/quaternions.jl | 5 ++++- Simulator/simulator.jl | 4 ++-- Simulator/tests/tests.jl | 5 +++-- 4 files changed, 14 insertions(+), 11 deletions(-) diff --git a/Simulator/kf.jl b/Simulator/kf.jl index 7c6c3ef..54fc55e 100644 --- a/Simulator/kf.jl +++ b/Simulator/kf.jl @@ -22,7 +22,7 @@ function f( δt::Float64 ) θ = norm(ω - β) * δt - r = (ω - β) / norm(ω - β) + r = normalize(ω - β) return L(q) * [cos(θ / 2); r * sin(θ / 2)] end @@ -50,7 +50,7 @@ function step( Pₚ = A * P * A' + W # Innovation - Q = quaternionToMatrix(qₚ) + Q = quaternionToMatrix(qₚ)' Z = [ᵇr_mag; ᵇr_sun] - [Q zeros(3, 3); zeros(3, 3) Q] * [ⁿr_mag; ⁿr_sun] C = [hat(ᵇr_mag) zeros(3, 3); hat(ᵇr_sun) zeros(3, 3)] S = C * Pₚ * C' + V @@ -62,7 +62,7 @@ function step( δx = L * Z ϕ = δx[1:3] δβ = δx[4:6] - θ = -norm(ϕ) + θ = norm(ϕ) r = normalize(ϕ) qᵤ = ⊙(qₚ, [cos(θ / 2); r * sin(θ / 2)]) βᵤ = e.β + δβ @@ -71,7 +71,6 @@ function step( e.q = qᵤ e.β = βᵤ e.P = Pᵤ - println("β: ", e.β, "δβ: ", δβ) return e end @@ -95,7 +94,7 @@ function simplest_step( qₚ = ⊙(q, [cos(θ / 2); r * sin(θ / 2)]) Pₚ = P + W # Innovation - Q = quaternionToMatrix(qₚ) + Q = quaternionToMatrix(qₚ)' Z = [ᵇr_mag; ᵇr_sun] - [Q zeros(3, 3); zeros(3, 3) Q] * [ⁿr_mag; ⁿr_sun] C = [hat(Q * ⁿr_mag); hat(Q * ⁿr_sun)] S = C * Pₚ * C' + V @@ -103,7 +102,7 @@ function simplest_step( L = Pₚ * C' * inv(S) # Update ϕ = L * Z - θ = -norm(ϕ) + θ = norm(ϕ) r = normalize(ϕ) e.q = ⊙(qₚ, [cos(θ / 2); r * sin(θ / 2)]) e.P = (I(3) - L * C) * Pₚ * (I(3) - L * C)' + L * V * L' diff --git a/Simulator/quaternions.jl b/Simulator/quaternions.jl index 362d1ee..24a754f 100644 --- a/Simulator/quaternions.jl +++ b/Simulator/quaternions.jl @@ -87,7 +87,7 @@ const H = [zeros(1, 3); I(3)] # Converts from a 3-element vector to a 4-elem const T = [1 0 0 0; 0 -1 0 0; 0 0 -1 0; 0 0 0 -1] # Forms the conjugate of q, i.e. q† = Tq -⊙(q₁, q₂) = L(q₂) * q₁ # Hamiltonian product +⊙(q₁, q₂) = L(q₁) * q₂ # Hamiltonian product """ quaternionToMatrix (q) Arguments: @@ -101,3 +101,6 @@ function quaternionToMatrix(q::Vector{Float64}) return I(3) + 2 * hat(v) * (s * I(3) + hat(v)) end +function qErr(q₁, q₂) + return norm((L(q₁)' * q₂)[2:4]) +end \ No newline at end of file diff --git a/Simulator/simulator.jl b/Simulator/simulator.jl index 778057e..24e49bb 100644 --- a/Simulator/simulator.jl +++ b/Simulator/simulator.jl @@ -262,7 +262,7 @@ module Simulator x₀ = initialize_orbit() println("intialized orbit!") # x₀[11:13] .=0 - # x₀[11:13] *= x₀[11:13].*10.0 # Spinning very fast + x₀[11:13] /= 4.0 # Spinning very fast J = [0.3 0 0; 0 0.3 0; 0 0 0.3] # Arbitrary inertia matrix for the Satellite t = Epoch(2020, 11, 30) # Starting time is Nov 30, 2020 @@ -302,7 +302,7 @@ module Simulator # rmag = (rand(3) * .01 .- 0.005) rsun = [0;0;0] rmag = [0;0;0] - Q = quaternionToMatrix(q) + Q = quaternionToMatrix(q)' body_sun = Q * (normalize(inertial_sun + rsun)) body_mag = Q * (normalize(inertial_mag + rmag)) if type == "simple" diff --git a/Simulator/tests/tests.jl b/Simulator/tests/tests.jl index 8d1e0f5..c4b9057 100644 --- a/Simulator/tests/tests.jl +++ b/Simulator/tests/tests.jl @@ -1,6 +1,6 @@ # [Simulator/test/tests.jl] -using Test, LinearAlgebra, SatelliteDynamics, Plots, JSON +using Test, LinearAlgebra, SatelliteDynamics, Plots, JSON, Random include("../simulator.jl"); using .Simulator @@ -107,6 +107,7 @@ end end @testset "MEKF/DeTumbling" begin + Random.seed!() function control_law(ω, b) b̂ = b / norm(b) k = 7e-4 @@ -128,7 +129,7 @@ end ylabel="Angular Velocity (rad/s)", labels=["ω" "s" "v1" "v2" "v3" "s'" "v1'" "v2'" "v3'"], linecolor=[:red :blue :green :purple :orange :blue :green :purple :orange], - linewidth=[1 1 1 1 1 1 1 1 1], + linewidth=[1 1 1 1 1 3 3 3 3], )) end From b5b8f7fd780b392ddc3ce09a4f59add489c39b93 Mon Sep 17 00:00:00 2001 From: thetazero Date: Mon, 27 Jun 2022 12:20:50 -0400 Subject: [PATCH 10/20] Plot cleanup --- Simulator/simulator.jl | 19 +++++++++++-------- Simulator/tests/tests.jl | 20 ++++++++++++++------ 2 files changed, 25 insertions(+), 14 deletions(-) diff --git a/Simulator/simulator.jl b/Simulator/simulator.jl index 24e49bb..dfc484c 100644 --- a/Simulator/simulator.jl +++ b/Simulator/simulator.jl @@ -262,16 +262,13 @@ module Simulator x₀ = initialize_orbit() println("intialized orbit!") # x₀[11:13] .=0 - x₀[11:13] /= 4.0 # Spinning very fast + # x₀[11:13] /= 4.0 # Spinning very fast J = [0.3 0 0; 0 0.3 0; 0 0 0.3] # Arbitrary inertia matrix for the Satellite t = Epoch(2020, 11, 30) # Starting time is Nov 30, 2020 dt = 0.01 # Time step, in seconds - q_hist = zeros(N, 9) - q_hist[1, 1] = norm(x₀[11:13]) - q_hist[1, 2:5] .= x₀[7:10] - q_hist[1, 6:9] .= x₀[7:10] + q_hist = zeros(N, 10) x = x₀ q₀ = x₀[7:10] @@ -283,16 +280,19 @@ module Simulator end println("q0:", q₀) kf = EKF(q₀, β₀, P₀) + q_hist[1, 1:4] .= x₀[7:10] + q_hist[1, 5:8] .= x₀[7:10] + q_hist[1, 9] = norm(x₀[11:13]) + q_hist[1, 10] = qErr(kf.q, q₀) for i = 1:N - 1 + println(i/N) r, v, q, ω = x[1:3], x[4:6], x[7:10], x[11:13] b = IGRF13(r, t) x = rk4(x, J, control_fn(ω, b), t, dt) t += dt # Don't forget to update time (not that it really matters...) - q_hist[i + 1, 1] = norm(ω) - q_hist[i + 1, 2:5] .= q r_sun = sun_position(t) inertial_sun = normalize(r_sun - r) @@ -310,7 +310,10 @@ module Simulator else step(kf, ω, dt, inertial_mag, inertial_sun, body_mag, body_sun) end - q_hist[i+1, 6:9] .= kf.q + q_hist[i + 1, 1:4] .= q + q_hist[i + 1, 5:8] .= kf.q + q_hist[i + 1, 9] = norm(ω) + q_hist[i + 1, 10] = qErr(kf.q, q) if norm(ω) < 0.1 q_hist = q_hist[1:i, :] diff --git a/Simulator/tests/tests.jl b/Simulator/tests/tests.jl index c4b9057..8fb5170 100644 --- a/Simulator/tests/tests.jl +++ b/Simulator/tests/tests.jl @@ -117,20 +117,28 @@ end end type = "smple" - data = my_sim_kf(control_law, 5000, type) + data = my_sim_kf(control_law, 100000, type) name = "MEKF/DeTumbling" if type == "simple" name = "MEKF(simplest)/DeTumbling" end - display(plot( - data, + p1 = (plot( + data[:, 1:9], title=name, xlabel="Time (s)", ylabel="Angular Velocity (rad/s)", - labels=["ω" "s" "v1" "v2" "v3" "s'" "v1'" "v2'" "v3'"], - linecolor=[:red :blue :green :purple :orange :blue :green :purple :orange], - linewidth=[1 1 1 1 1 3 3 3 3], + labels=["s" "v1" "v2" "v3" "s'" "v1'" "v2'" "v3'" "ω"], + linecolor=[:blue :green :purple :orange :blue :green :purple :orange :red], + linewidth=[1 1 1 1 3 3 3 3 1], )) + p2 = plot( + data[:, 10], + title="quaternion error", + xlabel="Time (s)", + ylabel="Radians", + ) + display(p1) + display(p2) end From cdbc42634f60553cfaf5f39a6d28c6fa9ab7ca4f Mon Sep 17 00:00:00 2001 From: thetazero Date: Tue, 28 Jun 2022 10:16:25 -0400 Subject: [PATCH 11/20] Fully working MEKF --- Simulator/gyro.jl | 20 +++++++++++++ Simulator/quaternions.jl | 30 ++++++++++--------- Simulator/simulator.jl | 63 +++++++++++++++++++++++++++------------- Simulator/tests/tests.jl | 33 ++++++++++++++------- Simulator/utils.jl | 10 +++++++ 5 files changed, 112 insertions(+), 44 deletions(-) create mode 100644 Simulator/gyro.jl create mode 100644 Simulator/utils.jl diff --git a/Simulator/gyro.jl b/Simulator/gyro.jl new file mode 100644 index 0000000..105a870 --- /dev/null +++ b/Simulator/gyro.jl @@ -0,0 +1,20 @@ +using Random, Distributions + +mutable struct Gyro + β::Vector{Float64} + Vβ::Matrix{Float64} + Vω::Matrix{Float64} +end + +""" update(g, ω) + Arguments: + - g: The gyro | Gyro + - ω: True angular velocity | [3,] + Returns: + - ω: Measured angular velocity | [3,] + +""" +function update(g::Gyro, ω::Vector{Float64}) + g.β += rand(MvNormal([0; 0; 0], g.Vβ)) + return ω + g.β + rand(MvNormal([0; 0; 0], g.Vω)) +end \ No newline at end of file diff --git a/Simulator/quaternions.jl b/Simulator/quaternions.jl index 24a754f..2e9eeb5 100644 --- a/Simulator/quaternions.jl +++ b/Simulator/quaternions.jl @@ -17,12 +17,12 @@ using LinearAlgebra # For I, norm Returns: - M: A [3 × 3] skew-symmetric matrix | [3, 3] """ -function hat(v) - M = [0.0 -v[3] v[2]; - v[3] 0.0 -v[1]; - -v[2] v[1] 0.0] +function hat(v) + M = [0.0 -v[3] v[2] + v[3] 0.0 -v[1] + -v[2] v[1] 0.0] - return M + return M end """ L(q) @@ -36,11 +36,11 @@ end Returns: - M: Left-side matrix representing the given quaternion | [4, 4] """ -function L(q) +function L(q) qₛ, qᵥ = q[1], q[2:4] - M = [qₛ -qᵥ'; - qᵥ qₛ*I(3) + hat(qᵥ)] + M = [qₛ -qᵥ' + qᵥ qₛ*I(3)+hat(qᵥ)] return M end @@ -57,11 +57,11 @@ end Returns: - M: Right-side matrix representing the given quaternion | [4, 4] """ -function R(q) +function R(q) qₛ, qᵥ = q[1], q[2:4] - M = [qₛ -qᵥ'; - qᵥ qₛ*I(3) - hat(qᵥ)] + M = [qₛ -qᵥ' + qᵥ qₛ*I(3)-hat(qᵥ)] return M end @@ -77,7 +77,7 @@ end Returns: - q̇: Derivative of attitude, parameterized as a quaternion | [4,] """ -function qdot(q, ω) +function qdot(q, ω) q̇ = 0.5 * L(q) * H * ω return q̇ end @@ -102,5 +102,9 @@ function quaternionToMatrix(q::Vector{Float64}) end function qErr(q₁, q₂) - return norm((L(q₁)' * q₂)[2:4]) + return norm((L(q₁)'*q₂)[2:4]) +end + +function eulerError(e1, e2) + return acos(dot(e1, e2) / (norm(e1) * norm(e2))) end \ No newline at end of file diff --git a/Simulator/simulator.jl b/Simulator/simulator.jl index dfc484c..72ca82b 100644 --- a/Simulator/simulator.jl +++ b/Simulator/simulator.jl @@ -10,6 +10,8 @@ module Simulator include("mag_field.jl") include("quaternions.jl") include("kf.jl") + include("gyro.jl") + include("utils.jl") export initialize_orbit, intialize_orbit_oe # Generate valid initial conditions for a satellite @@ -266,9 +268,14 @@ module Simulator J = [0.3 0 0; 0 0.3 0; 0 0 0.3] # Arbitrary inertia matrix for the Satellite t = Epoch(2020, 11, 30) # Starting time is Nov 30, 2020 - dt = 0.01 # Time step, in seconds + dt = 0.1 # Time step, in seconds + + q_true = zeros(N, 4) + q_predicted = zeros(N, 4) + ω_true = zeros(N, 1) + q_error = zeros(N, 1) + β_error = zeros(N, 1) - q_hist = zeros(N, 10) x = x₀ q₀ = x₀[7:10] @@ -280,48 +287,64 @@ module Simulator end println("q0:", q₀) kf = EKF(q₀, β₀, P₀) - q_hist[1, 1:4] .= x₀[7:10] - q_hist[1, 5:8] .= x₀[7:10] - q_hist[1, 9] = norm(x₀[11:13]) - q_hist[1, 10] = qErr(kf.q, q₀) + + q_true[1, :] .= q₀ + q_predicted[1, :] .= q₀ + ω_true[1] = norm(x₀[11:13]) + q_error[1] = qErr(kf.q, q₀) + β_error[1] = 0 + + Vβ = I(3) * 0.01 + Vω = I(3) * 0.1 + gyro = Gyro(β₀, Vβ, Vω) + + progress_step = 0.01 + cur = progress_step for i = 1:N - 1 - println(i/N) + if i/N >= cur + println(string(cur * 100) * "%") + cur += progress_step + end r, v, q, ω = x[1:3], x[4:6], x[7:10], x[11:13] + ω_read = update(gyro, ω) b = IGRF13(r, t) x = rk4(x, J, control_fn(ω, b), t, dt) - t += dt # Don't forget to update time (not that it really matters...) + t += dt # Don't forget to update time r_sun = sun_position(t) inertial_sun = normalize(r_sun - r) inertial_mag = normalize(b) - # rsun = (rand(3) * .01 .- 0.005) - # rmag = (rand(3) * .01 .- 0.005) rsun = [0;0;0] rmag = [0;0;0] Q = quaternionToMatrix(q)' - body_sun = Q * (normalize(inertial_sun + rsun)) - body_mag = Q * (normalize(inertial_mag + rmag)) + body_sun = randomMatrix(0.001) * Q * (normalize(inertial_sun + rsun)) + body_mag = randomMatrix(0.001) * Q * (normalize(inertial_mag + rmag)) if type == "simple" - simplest_step(kf, ω, dt, inertial_mag, inertial_sun, body_mag, body_sun) + simplest_step(kf, ω_read, dt, inertial_mag, inertial_sun, body_mag, body_sun) else - step(kf, ω, dt, inertial_mag, inertial_sun, body_mag, body_sun) + step(kf, ω_read, dt, inertial_mag, inertial_sun, body_mag, body_sun) end - q_hist[i + 1, 1:4] .= q - q_hist[i + 1, 5:8] .= kf.q - q_hist[i + 1, 9] = norm(ω) - q_hist[i + 1, 10] = qErr(kf.q, q) + q_true[i + 1,:] .= q + q_predicted[i + 1,:] .= kf.q + ω_true[i + 1] = norm(ω) + q_error[i + 1] = qErr(kf.q, q) + β_error[i + 1] = eulerError(gyro.β, kf.β) if norm(ω) < 0.1 - q_hist = q_hist[1:i, :] + q_true = q_true[1:i, :] + q_predicted = q_predicted[1:i, :] + ω_true = ω_true[1:i] + q_error = q_error[1:i] + β_error = β_error[1:i] break end end - return q_hist + return q_true, q_predicted, ω_true, q_error, β_error end end \ No newline at end of file diff --git a/Simulator/tests/tests.jl b/Simulator/tests/tests.jl index 8fb5170..7a95533 100644 --- a/Simulator/tests/tests.jl +++ b/Simulator/tests/tests.jl @@ -117,28 +117,39 @@ end end type = "smple" - data = my_sim_kf(control_law, 100000, type) + q_true, q_predicted, ω_true, q_error, β_error = my_sim_kf(control_law, 100000, type) name = "MEKF/DeTumbling" if type == "simple" name = "MEKF(simplest)/DeTumbling" end - p1 = (plot( - data[:, 1:9], + q_plot = (plot( + hcat(q_true, q_predicted), title=name, xlabel="Time (s)", ylabel="Angular Velocity (rad/s)", - labels=["s" "v1" "v2" "v3" "s'" "v1'" "v2'" "v3'" "ω"], - linecolor=[:blue :green :purple :orange :blue :green :purple :orange :red], - linewidth=[1 1 1 1 3 3 3 3 1], + labels=["s" "v1" "v2" "v3" "s'" "v1'" "v2'" "v3'"], + linecolor=[:blue :green :purple :orange :blue :green :purple :orange], + linewidth=[1 1 1 1 1 1 1 1], )) - p2 = plot( - data[:, 10], - title="quaternion error", + q_error_plot = plot( + q_error, + title="Quaternion Error", xlabel="Time (s)", ylabel="Radians", ) - display(p1) - display(p2) + ω_plot = plot( + ω_true, + title="Angular Velocity", + xlabel="Time (s)", + ylabel="Radians/s" + ) + β_plot = plot( + β_error, + title="Gyro Bias Error", + xlabel="Time (s)", + ylabel="Radians" + ) + display(plot(q_plot, q_error_plot, ω_plot, β_plot, layout=(2,2), dpi=300)) end diff --git a/Simulator/utils.jl b/Simulator/utils.jl new file mode 100644 index 0000000..a06857e --- /dev/null +++ b/Simulator/utils.jl @@ -0,0 +1,10 @@ +function hat(ω::Vector) + return [0 ω[3] -ω[2] + -ω[3] 0 ω[1] + ω[2] -ω[1] 0] +end + +function randomMatrix(covariance) + ϕ = √covariance * randn(3) + return exp(hat(ϕ)) +end \ No newline at end of file From d0a2c20cc52ff16175c76d2c615e3faa8438b2b6 Mon Sep 17 00:00:00 2001 From: thetazero Date: Tue, 28 Jun 2022 13:42:01 -0400 Subject: [PATCH 12/20] Fixed ordering & changed time step --- Simulator/simulator.jl | 26 +++++++++++--------------- 1 file changed, 11 insertions(+), 15 deletions(-) diff --git a/Simulator/simulator.jl b/Simulator/simulator.jl index 72ca82b..7827096 100644 --- a/Simulator/simulator.jl +++ b/Simulator/simulator.jl @@ -268,7 +268,7 @@ module Simulator J = [0.3 0 0; 0 0.3 0; 0 0 0.3] # Arbitrary inertia matrix for the Satellite t = Epoch(2020, 11, 30) # Starting time is Nov 30, 2020 - dt = 0.1 # Time step, in seconds + dt = 0.01 # Time step, in seconds q_true = zeros(N, 4) q_predicted = zeros(N, 4) @@ -288,12 +288,6 @@ module Simulator println("q0:", q₀) kf = EKF(q₀, β₀, P₀) - q_true[1, :] .= q₀ - q_predicted[1, :] .= q₀ - ω_true[1] = norm(x₀[11:13]) - q_error[1] = qErr(kf.q, q₀) - β_error[1] = 0 - Vβ = I(3) * 0.01 Vω = I(3) * 0.1 gyro = Gyro(β₀, Vβ, Vω) @@ -306,14 +300,12 @@ module Simulator println(string(cur * 100) * "%") cur += progress_step end + r, v, q, ω = x[1:3], x[4:6], x[7:10], x[11:13] ω_read = update(gyro, ω) b = IGRF13(r, t) - x = rk4(x, J, control_fn(ω, b), t, dt) - t += dt # Don't forget to update time - r_sun = sun_position(t) inertial_sun = normalize(r_sun - r) inertial_mag = normalize(b) @@ -328,11 +320,12 @@ module Simulator else step(kf, ω_read, dt, inertial_mag, inertial_sun, body_mag, body_sun) end - q_true[i + 1,:] .= q - q_predicted[i + 1,:] .= kf.q - ω_true[i + 1] = norm(ω) - q_error[i + 1] = qErr(kf.q, q) - β_error[i + 1] = eulerError(gyro.β, kf.β) + q_true[i,:] .= q + q_predicted[i,:] .= kf.q + ω_true[i] = norm(ω) + q_error[i] = qErr(kf.q, q) + β_error[i] = eulerError(gyro.β, kf.β) + if norm(ω) < 0.1 q_true = q_true[1:i, :] @@ -342,6 +335,9 @@ module Simulator β_error = β_error[1:i] break end + # propogate everything forward + x = rk4(x, J, control_fn(ω_read - kf.β, b), t, dt) + t += dt # Don't forget to update time end return q_true, q_predicted, ω_true, q_error, β_error From 844e6eaf3e6d6341cf2063aac3953bd44dab3371 Mon Sep 17 00:00:00 2001 From: thetazero Date: Wed, 6 Jul 2022 16:17:05 -0400 Subject: [PATCH 13/20] DeTumblio + MEKF work --- .gitignore | 6 ++++++ Simulator/simulator.jl | 24 +++++++++------------- Simulator/tests/tests.jl | 43 +++++++++++++++++++++++++--------------- 3 files changed, 42 insertions(+), 31 deletions(-) diff --git a/.gitignore b/.gitignore index b6e4761..3743b8b 100644 --- a/.gitignore +++ b/.gitignore @@ -127,3 +127,9 @@ dmypy.json # Pyre type checker .pyre/ + +# LaTex Build Files +*.aux +*.fdb_latexmk +*.fls +*.gz \ No newline at end of file diff --git a/Simulator/simulator.jl b/Simulator/simulator.jl index 7827096..b813bfd 100644 --- a/Simulator/simulator.jl +++ b/Simulator/simulator.jl @@ -47,7 +47,7 @@ module Simulator r, v, q, ω = x[1:3], x[4:6], x[7:10], x[11:13] if norm(r) < Rₑ - error("Error: Satellite impacted Earth at time $t!") + error("Error: Satellite impacted Earth at time $t !") end @@ -244,7 +244,8 @@ module Simulator for i = 1:N - 1 r, v, q, ω = x[1:3], x[4:6], x[7:10], x[11:13] b = IGRF13(r, t) - x = rk4(x, J, control_fn(ω, b), t, dt) + m = control_fn(ω, b) + x = rk4(x, J, cross(m, b), t, dt) t += dt # Don't forget to update time (not that it really matters...) # q_hist[i + 1, :] .= x[7:10] q_hist[i + 1, 1:3] .= x[11:13] @@ -260,11 +261,9 @@ module Simulator end - function my_sim_kf(control_fn, N, type) + function my_sim_kf(control_fn, N) x₀ = initialize_orbit() println("intialized orbit!") - # x₀[11:13] .=0 - # x₀[11:13] /= 4.0 # Spinning very fast J = [0.3 0 0; 0 0.3 0; 0 0 0.3] # Arbitrary inertia matrix for the Satellite t = Epoch(2020, 11, 30) # Starting time is Nov 30, 2020 @@ -282,9 +281,6 @@ module Simulator β₀ = [0; 0; 0] P₀ = I(6) - if type == "simple" - P₀ = I(3) - end println("q0:", q₀) kf = EKF(q₀, β₀, P₀) @@ -292,7 +288,7 @@ module Simulator Vω = I(3) * 0.1 gyro = Gyro(β₀, Vβ, Vω) - progress_step = 0.01 + progress_step = 0.001 cur = progress_step for i = 1:N - 1 @@ -315,11 +311,7 @@ module Simulator Q = quaternionToMatrix(q)' body_sun = randomMatrix(0.001) * Q * (normalize(inertial_sun + rsun)) body_mag = randomMatrix(0.001) * Q * (normalize(inertial_mag + rmag)) - if type == "simple" - simplest_step(kf, ω_read, dt, inertial_mag, inertial_sun, body_mag, body_sun) - else - step(kf, ω_read, dt, inertial_mag, inertial_sun, body_mag, body_sun) - end + step(kf, ω_read, dt, inertial_mag, inertial_sun, body_mag, body_sun) q_true[i,:] .= q q_predicted[i,:] .= kf.q ω_true[i] = norm(ω) @@ -336,7 +328,9 @@ module Simulator break end # propogate everything forward - x = rk4(x, J, control_fn(ω_read - kf.β, b), t, dt) + control, dt = control_fn(ω_read - kf.β, b) + dt /= 1e9 # convert from nanoseconds to seconds + x = rk4(x, J, cross(control, b), t, dt) t += dt # Don't forget to update time end diff --git a/Simulator/tests/tests.jl b/Simulator/tests/tests.jl index 7a95533..32d874f 100644 --- a/Simulator/tests/tests.jl +++ b/Simulator/tests/tests.jl @@ -113,15 +113,11 @@ end k = 7e-4 M = -k * (I(3) - b̂*b̂')*ω m = 1 / (dot(b, b)) * cross(b, M) - return cross(m, b) + return m, 1e8 # half a second in nano seconds end - type = "smple" - q_true, q_predicted, ω_true, q_error, β_error = my_sim_kf(control_law, 100000, type) + q_true, q_predicted, ω_true, q_error, β_error = my_sim_kf(control_law, 100000) name = "MEKF/DeTumbling" - if type == "simple" - name = "MEKF(simplest)/DeTumbling" - end q_plot = (plot( hcat(q_true, q_predicted), title=name, @@ -155,12 +151,13 @@ end @testset "DeTumbling" begin + println("DeTumbling") function control_law(ω, b) b̂ = b / norm(b) k = 7e-4 M = -k * (I(3) - b̂*b̂')*ω m = 1 / (dot(b, b)) * cross(b, M) - return cross(m, b) + return m end data = my_sim(control_law) @@ -174,7 +171,7 @@ end inp = Pipe() out = Pipe() err = Pipe() - proc = run(Cmd(`python3 -u main.py`, dir = "/home/thetazero/Documents/pycubed/software_example_beepsat/state_machine/build/" ), inp, out, err, wait = false) + proc = run(Cmd(`faketime -f '-0 x1' python3 -u main.py`, dir = "/home/thetazero/Documents/pycubed/software_example_beepsat/state_machine/build/" ), inp, out, err, wait = false) close(out.in) close(err.in) Base.start_reading(out.out) @@ -185,28 +182,42 @@ end i=0 function control_law(ω, b) - println("wrote to sensors") control = [0,0,0] - while true + got_control = false + got_time = false + + cur_time = 0 + last_time = time() * 1e9 - 1 # Time in nanoseconds + while !(got_control && got_time) write(inp, ">>>ω"*string(ω)*"\n") write(inp, ">>>b"*string(b)*"\n") - write(inp, ">>>?\n") input = readline(proc.out) - println(input) + # println(input) if length(input) >= 4 && input[1:3] == ">>>" - if input[4] == 'M' + if input[4] == 'm' control = (input[5:length(input)]) control = JSON.parse(control) + got_control = true + elseif input[4] == 't' + cur_time = parse(Int64, input[5:length(input)]) + got_time = true end + end + if got_control && got_time break end end i+=1 - println(i) - return control + + println("I: ", i) + + dt = norm(cur_time - last_time) + last_time = cur_time + println("control: ", control, " dt: ", dt) + return control, dt end - data = my_sim(control_law) + data = my_sim_kf(control_law, 1000000) display(plot(data, title="DeTumbleIO", xlabel="Time (s)", ylabel="Angular Velocity (rad/s)", labels=["ω1" "ω2" "ω3" "ω"])) end From c7a8c78223a610d835ee0cc201057b21762d5aad Mon Sep 17 00:00:00 2001 From: thetazero Date: Thu, 7 Jul 2022 10:18:56 -0400 Subject: [PATCH 14/20] Works if detumble task runs at same rate as MEKF --- Simulator/simulator.jl | 2 +- Simulator/tests/tests.jl | 49 ++++++++++++++++++++++++++++++++-------- 2 files changed, 40 insertions(+), 11 deletions(-) diff --git a/Simulator/simulator.jl b/Simulator/simulator.jl index b813bfd..5248d6d 100644 --- a/Simulator/simulator.jl +++ b/Simulator/simulator.jl @@ -319,7 +319,7 @@ module Simulator β_error[i] = eulerError(gyro.β, kf.β) - if norm(ω) < 0.1 + if norm(ω) < 0.01 q_true = q_true[1:i, :] q_predicted = q_predicted[1:i, :] ω_true = ω_true[1:i] diff --git a/Simulator/tests/tests.jl b/Simulator/tests/tests.jl index 32d874f..c973390 100644 --- a/Simulator/tests/tests.jl +++ b/Simulator/tests/tests.jl @@ -1,6 +1,6 @@ # [Simulator/test/tests.jl] -using Test, LinearAlgebra, SatelliteDynamics, Plots, JSON, Random +using Test, LinearAlgebra, SatelliteDynamics, Plots, JSON, Random, Dates include("../simulator.jl"); using .Simulator @@ -167,11 +167,12 @@ end @testset "DeTumbleIO" begin + Random.seed!() println("DeTumbleIO") inp = Pipe() out = Pipe() err = Pipe() - proc = run(Cmd(`faketime -f '-0 x1' python3 -u main.py`, dir = "/home/thetazero/Documents/pycubed/software_example_beepsat/state_machine/build/" ), inp, out, err, wait = false) + proc = run(Cmd(`faketime -f '-0 x100' python3 -u main.py`, dir = "/home/thetazero/Documents/pycubed/software_example_beepsat/state_machine/build/" ), inp, out, err, wait = false) close(out.in) close(err.in) Base.start_reading(out.out) @@ -181,16 +182,16 @@ end i=0 + last_time = time() * Int(1e9) # Time in nanoseconds function control_law(ω, b) control = [0,0,0] got_control = false got_time = false cur_time = 0 - last_time = time() * 1e9 - 1 # Time in nanoseconds + write(inp, ">>>ω"*string(ω)*"\n") + write(inp, ">>>b"*string(b)*"\n") while !(got_control && got_time) - write(inp, ">>>ω"*string(ω)*"\n") - write(inp, ">>>b"*string(b)*"\n") input = readline(proc.out) # println(input) if length(input) >= 4 && input[1:3] == ">>>" @@ -209,15 +210,43 @@ end end i+=1 - println("I: ", i) - + # println("I: ", i) + # println("==========last_time: ", last_time, " cur_time: ", cur_time, " dt: ", cur_time-last_time) dt = norm(cur_time - last_time) last_time = cur_time - println("control: ", control, " dt: ", dt) + # println("control: ", control, " dt: ", dt/1e9) return control, dt end - data = my_sim_kf(control_law, 1000000) - display(plot(data, title="DeTumbleIO", xlabel="Time (s)", ylabel="Angular Velocity (rad/s)", labels=["ω1" "ω2" "ω3" "ω"])) + q_true, q_predicted, ω_true, q_error, β_error = my_sim_kf(control_law, 10000) + name = "MEKF/DeTumbling(IO)" + q_plot = (plot( + hcat(q_true, q_predicted), + title=name, + xlabel="Time (s)", + ylabel="Angular Velocity (rad/s)", + labels=["s" "v1" "v2" "v3" "s'" "v1'" "v2'" "v3'"], + linecolor=[:blue :green :purple :orange :blue :green :purple :orange], + linewidth=[1 1 1 1 1 1 1 1], + )) + q_error_plot = plot( + q_error, + title="Quaternion Error", + xlabel="Time (s)", + ylabel="Radians", + ) + ω_plot = plot( + ω_true, + title="Angular Velocity", + xlabel="Time (s)", + ylabel="Radians/s" + ) + β_plot = plot( + β_error, + title="Gyro Bias Error", + xlabel="Time (s)", + ylabel="Radians" + ) + display(plot(q_plot, q_error_plot, ω_plot, β_plot, layout=(2,2), dpi=300)) end From c6b5da251d2560a181a7b17da76cb002a6edebfe Mon Sep 17 00:00:00 2001 From: thetazero Date: Tue, 12 Jul 2022 11:25:33 -0400 Subject: [PATCH 15/20] Tranpose hat opperator to be inline with Zac's standard --- Simulator/quaternions.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Simulator/quaternions.jl b/Simulator/quaternions.jl index 2e9eeb5..ad4e9b6 100644 --- a/Simulator/quaternions.jl +++ b/Simulator/quaternions.jl @@ -18,9 +18,9 @@ using LinearAlgebra # For I, norm - M: A [3 × 3] skew-symmetric matrix | [3, 3] """ function hat(v) - M = [0.0 -v[3] v[2] - v[3] 0.0 -v[1] - -v[2] v[1] 0.0] + M = [0.0 v[3] -v[2] + -v[3] 0.0 v[1] + v[2] -v[1] 0.0] return M end From 4c1ea35756b1e77122e3155dda5bbdcbf0a5cf40 Mon Sep 17 00:00:00 2001 From: thetazero Date: Tue, 12 Jul 2022 11:49:41 -0400 Subject: [PATCH 16/20] Fix hat functions --- Simulator/kf.jl | 6 +++--- Simulator/quaternions.jl | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/Simulator/kf.jl b/Simulator/kf.jl index 54fc55e..a63bfe0 100644 --- a/Simulator/kf.jl +++ b/Simulator/kf.jl @@ -10,9 +10,9 @@ mutable struct EKF end function hat(ω::Vector) - return [0 ω[3] -ω[2] - -ω[3] 0 ω[1] - ω[2] -ω[1] 0] + return [0 -ω[3] ω[2] + ω[3] 0 -ω[1] + -ω[2] ω[1] 0] end function f( diff --git a/Simulator/quaternions.jl b/Simulator/quaternions.jl index ad4e9b6..2e9eeb5 100644 --- a/Simulator/quaternions.jl +++ b/Simulator/quaternions.jl @@ -18,9 +18,9 @@ using LinearAlgebra # For I, norm - M: A [3 × 3] skew-symmetric matrix | [3, 3] """ function hat(v) - M = [0.0 v[3] -v[2] - -v[3] 0.0 v[1] - v[2] -v[1] 0.0] + M = [0.0 -v[3] v[2] + v[3] 0.0 -v[1] + -v[2] v[1] 0.0] return M end From 4d2e1acf9ced0f90d5e45171c18577f2853133a2 Mon Sep 17 00:00:00 2001 From: thetazero Date: Thu, 14 Jul 2022 16:02:49 -0400 Subject: [PATCH 17/20] Update KF to use rodgriguez formula instead of exponenentiation --- Simulator/kf.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/Simulator/kf.jl b/Simulator/kf.jl index a63bfe0..70dd4bf 100644 --- a/Simulator/kf.jl +++ b/Simulator/kf.jl @@ -42,7 +42,11 @@ function step( V = I(6) * 1e-6 # some other noise # Predict: qₚ = f(q, β, ω, δt) # β remains constant - R = exp(-hat(ω - β) * δt) + # R = exp(-hat(ω - β) * δt) + v = - (ω - β) + mag = norm(v) + v̂ = hat(v/mag) + R = I(3) + (v̂) * sin(mag * δt) + (v̂^2) * (1 - cos(mag * δt)) # equivalent to e^{-hat(ω - β) * δt} A = [ R (-δt*I(3)) zeros(3, 3) I(3) From fca9a19957f625dd1794823b058ffc95877f3c06 Mon Sep 17 00:00:00 2001 From: thetazero Date: Thu, 14 Jul 2022 16:09:40 -0400 Subject: [PATCH 18/20] Docs --- Simulator/kf.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Simulator/kf.jl b/Simulator/kf.jl index 70dd4bf..00e93d8 100644 --- a/Simulator/kf.jl +++ b/Simulator/kf.jl @@ -42,7 +42,9 @@ function step( V = I(6) * 1e-6 # some other noise # Predict: qₚ = f(q, β, ω, δt) # β remains constant + # The following is equivalent to: # R = exp(-hat(ω - β) * δt) + # By the https://en.wikipedia.org/wiki/Rodrigues'_rotation_formula v = - (ω - β) mag = norm(v) v̂ = hat(v/mag) From 8c246a22922fff86a451bfa66859b19ccd10bacb Mon Sep 17 00:00:00 2001 From: thetazero Date: Wed, 20 Jul 2022 10:43:56 -0400 Subject: [PATCH 19/20] Fix MEKF to use appriopriate hat opperator --- Simulator/kf.jl | 33 +++++++++++++++++++++++++++------ Simulator/quaternions.jl | 3 +++ Simulator/utils.jl | 6 ------ 3 files changed, 30 insertions(+), 12 deletions(-) diff --git a/Simulator/kf.jl b/Simulator/kf.jl index 00e93d8..6852554 100644 --- a/Simulator/kf.jl +++ b/Simulator/kf.jl @@ -9,11 +9,11 @@ mutable struct EKF P::Matrix{Float64} end -function hat(ω::Vector) - return [0 -ω[3] ω[2] - ω[3] 0 -ω[1] - -ω[2] ω[1] 0] -end +# function hat(ω::Vector) +# return [0 -ω[3] ω[2] +# ω[3] 0 -ω[1] +# -ω[2] ω[1] 0] +# end function f( q::Vector{Float64}, @@ -23,6 +23,8 @@ function f( ) θ = norm(ω - β) * δt r = normalize(ω - β) + # println("q: ", q) + # println("L(q): ", L(q)) return L(q) * [cos(θ / 2); r * sin(θ / 2)] end @@ -33,8 +35,21 @@ function step( ⁿr_mag::Vector{Float64}, ⁿr_sun::Vector{Float64}, ᵇr_mag::Vector{Float64}, - ᵇr_sun::Vector{Float64} + ᵇr_sun::Vector{Float64}, + logging::Bool = false, ) + if logging + println("Before: ") + println("e.q: ", e.q) + println("e.β: ", e.β) + println("e.P: ", e.P) + println("ω: ", ω) + println("δt: ", δt) + println("ⁿr_mag: ", ⁿr_mag) + println("ⁿr_sun: ", ⁿr_sun) + println("ᵇr_mag: ", ᵇr_mag) + println("ᵇr_sun: ", ᵇr_sun) + end q = e.q β = e.β P = e.P @@ -77,6 +92,12 @@ function step( e.q = qᵤ e.β = βᵤ e.P = Pᵤ + if logging + println("After: ") + println("e.q: ", e.q) + println("e.β: ", e.β) + println("e.P: ", e.P) + end return e end diff --git a/Simulator/quaternions.jl b/Simulator/quaternions.jl index 2e9eeb5..d5cd415 100644 --- a/Simulator/quaternions.jl +++ b/Simulator/quaternions.jl @@ -39,6 +39,9 @@ end function L(q) qₛ, qᵥ = q[1], q[2:4] + # println("qᵥ: ", qᵥ) + # println(qᵥ[1]," ",qᵥ[2]," ",qᵥ[3]) + # println("hat(qᵥ): ", hat(qᵥ)) M = [qₛ -qᵥ' qᵥ qₛ*I(3)+hat(qᵥ)] diff --git a/Simulator/utils.jl b/Simulator/utils.jl index a06857e..524e74f 100644 --- a/Simulator/utils.jl +++ b/Simulator/utils.jl @@ -1,9 +1,3 @@ -function hat(ω::Vector) - return [0 ω[3] -ω[2] - -ω[3] 0 ω[1] - ω[2] -ω[1] 0] -end - function randomMatrix(covariance) ϕ = √covariance * randn(3) return exp(hat(ϕ)) From 66c766f35e3f1aca872e4abcd0cf201378fe9f3d Mon Sep 17 00:00:00 2001 From: thetazero Date: Tue, 3 Jan 2023 11:15:41 -0800 Subject: [PATCH 20/20] not sure what changed --- Simulator/kf.jl | 12 ++---------- Simulator/quaternions.jl | 3 --- Simulator/tests/tests.jl | 2 +- 3 files changed, 3 insertions(+), 14 deletions(-) diff --git a/Simulator/kf.jl b/Simulator/kf.jl index 6852554..993fca2 100644 --- a/Simulator/kf.jl +++ b/Simulator/kf.jl @@ -9,12 +9,6 @@ mutable struct EKF P::Matrix{Float64} end -# function hat(ω::Vector) -# return [0 -ω[3] ω[2] -# ω[3] 0 -ω[1] -# -ω[2] ω[1] 0] -# end - function f( q::Vector{Float64}, β::Vector{Float64}, @@ -23,8 +17,6 @@ function f( ) θ = norm(ω - β) * δt r = normalize(ω - β) - # println("q: ", q) - # println("L(q): ", L(q)) return L(q) * [cos(θ / 2); r * sin(θ / 2)] end @@ -63,9 +55,9 @@ function step( v = - (ω - β) mag = norm(v) v̂ = hat(v/mag) - R = I(3) + (v̂) * sin(mag * δt) + (v̂^2) * (1 - cos(mag * δt)) # equivalent to e^{-hat(ω - β) * δt} + R = I(3) + (v̂) * sin(mag * δt) + (v̂*v̂) * (1 - cos(mag * δt)) # equivalent to e^{-hat(ω - β) * δt} A = [ - R (-δt*I(3)) + R (-δt*I(3)) zeros(3, 3) I(3) ] Pₚ = A * P * A' + W diff --git a/Simulator/quaternions.jl b/Simulator/quaternions.jl index d5cd415..2e9eeb5 100644 --- a/Simulator/quaternions.jl +++ b/Simulator/quaternions.jl @@ -39,9 +39,6 @@ end function L(q) qₛ, qᵥ = q[1], q[2:4] - # println("qᵥ: ", qᵥ) - # println(qᵥ[1]," ",qᵥ[2]," ",qᵥ[3]) - # println("hat(qᵥ): ", hat(qᵥ)) M = [qₛ -qᵥ' qᵥ qₛ*I(3)+hat(qᵥ)] diff --git a/Simulator/tests/tests.jl b/Simulator/tests/tests.jl index c973390..48c4474 100644 --- a/Simulator/tests/tests.jl +++ b/Simulator/tests/tests.jl @@ -218,7 +218,7 @@ end return control, dt end - q_true, q_predicted, ω_true, q_error, β_error = my_sim_kf(control_law, 10000) + q_true, q_predicted, ω_true, q_error, β_error = my_sim_kf(control_law, 50000) name = "MEKF/DeTumbling(IO)" q_plot = (plot( hcat(q_true, q_predicted),