From 9843ea8b31437b8d44d526be8f0213af784a80db Mon Sep 17 00:00:00 2001 From: Thomas Hiller Date: Thu, 30 Jun 2022 12:14:33 +0200 Subject: [PATCH] see CHANGELOG.md --- AUTHORS.md | 3 + BLOCHUS/BLOCHUS.m | 6 +- BLOCHUS/BLOCHUS_loadDefaults.m | 2 +- CHANGELOG.md | 12 ++ README.md | 11 +- doc/blochus/BLOCHUS/BLOCHUS.html | 6 +- doc/blochus/BLOCHUS/BLOCHUS_loadDefaults.html | 16 +- .../functions/blochsim/fcn_BLOCHUS_ode.html | 2 +- .../functions/blochsim/getMIDI_Tx.html | 4 +- .../blochsim/getPulseTimeSeries.html | 89 ++++---- .../functions/blochsim/getRampAmplitude.html | 137 ++++++------ .../getRotationMatrixFromAngleandAxis.html | 124 +++++++---- .../getRotationMatrixFromVectors.html | 78 +++---- .../interface/getRampParameters.html | 4 +- functions/blochsim/fcn_BLOCHUS_ode.m | 2 +- functions/blochsim/getPulseTimeSeries.m | 7 +- functions/blochsim/getRampAmplitude.m | 3 + .../getRotationMatrixFromAngleandAxis.m | 70 +++++-- .../blochsim/getRotationMatrixFromVectors.m | 10 +- pyBLOCHUS/README.md | 4 +- pyBLOCHUS/pyBLOCHUS/blochus_pulse.py | 7 +- scripts/example3.m | 2 +- scripts/example4.m | 2 +- scripts/example6.m | 198 ++++++++++++++++++ scripts/example6_ref.png | Bin 0 -> 54675 bytes scripts/ramp_2ms.mat | Bin 0 -> 2885 bytes scripts/ramp_3ms.mat | Bin 0 -> 4232 bytes scripts/ramp_3ms_raw.mat | Bin 0 -> 13517 bytes 28 files changed, 552 insertions(+), 247 deletions(-) create mode 100644 AUTHORS.md create mode 100644 scripts/example6.m create mode 100644 scripts/example6_ref.png create mode 100644 scripts/ramp_2ms.mat create mode 100644 scripts/ramp_3ms.mat create mode 100644 scripts/ramp_3ms_raw.mat diff --git a/AUTHORS.md b/AUTHORS.md new file mode 100644 index 0000000..545bd28 --- /dev/null +++ b/AUTHORS.md @@ -0,0 +1,3 @@ +**BLOCHUS** is developed by + +* Thomas Hiller (since 2019), thomas.hiller@bgr.de \ No newline at end of file diff --git a/BLOCHUS/BLOCHUS.m b/BLOCHUS/BLOCHUS.m index 85c84d3..c34c9d1 100644 --- a/BLOCHUS/BLOCHUS.m +++ b/BLOCHUS/BLOCHUS.m @@ -33,10 +33,10 @@ %------------- BEGIN CODE -------------- %% GUI 'header' info and default GUI settings -myui.version = '0.1.4'; -myui.date = '25.09.2020'; +myui.version = '0.1.5'; +myui.date = '30.06.2022'; myui.author = 'Thomas Hiller'; -myui.email = 'thomas.hiller[at]leibniz-liag.de'; +myui.email = 'thomas.hiller[at]bgr.de'; myui.fontsize = 9; myui.axfontsize = 11; diff --git a/BLOCHUS/BLOCHUS_loadDefaults.m b/BLOCHUS/BLOCHUS_loadDefaults.m index d7e96b0..29d460b 100644 --- a/BLOCHUS/BLOCHUS_loadDefaults.m +++ b/BLOCHUS/BLOCHUS_loadDefaults.m @@ -118,7 +118,7 @@ % pre-polarization B-field amplitude in units of [B0] init.PrePolFactor = [100 1e-6 1e6]; % polar angle of pre-polarization field [deg] -init.PrePolTheta = [90 0 360]; +init.PrePolTheta = [90 1e-3 360]; % azimuthal angle of pre-polarization field [deg] init.PrePolPhi = [0 0 360]; % pre-polarization B-field switch amplitude in units of [B0] diff --git a/CHANGELOG.md b/CHANGELOG.md index 99495ae..9f29fdf 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,16 @@ # Changelog +## [0.1.5] - 2022-06-30 + +### Added +- custom pre-polarization current switch-off ramps are now supported via an interpolation routine (so far only via script and in MATLABTM) + +### Changed +- now multiple angle-axis-pairs are allowed as input for function `getRotationMatrixFromAngleandAxis` + +### Fixed +- fixed an issue in `getRotationMatrixFromVectors` for very small angles between vectors A and B (also edge case like parallel and anti-parallel vectors are now considered) + ## [0.1.4] - 2020-09-25 ### Added @@ -43,6 +54,7 @@ Initial Version +[0.1.5]: https://github.com/ThoHiller/nmr-blochus/compare/v0.1.4...v0.1.5 [0.1.4]: https://github.com/ThoHiller/nmr-blochus/compare/v0.1.3...v0.1.4 [0.1.3]: https://github.com/ThoHiller/nmr-blochus/compare/v0.1.2...v0.1.3 [0.1.2]: https://github.com/ThoHiller/nmr-blochus/compare/v.0.1.1...v0.1.2 diff --git a/README.md b/README.md index 2f909b4..0f3618c 100644 --- a/README.md +++ b/README.md @@ -48,7 +48,7 @@ In order to work properly you need to meet the following requirements: #### Operating System -I tested it successfully under Windows 7 (64bit) and 10 (64bit) with Matlab R2017a and newer. Always with the latest version of the GUI Layout Toolbox (current version is v2.3.4) +I tested it successfully under Windows 7 (64bit) and 10 (64bit) with Matlab R2017a and newer. Always with the latest version of the GUI Layout Toolbox (current version is v2.3.5) **NOTE:** So far I did not test anything on Linux or a Mac. If you get it to work on either of the two systems (which it basically should I guess) please let me know. @@ -63,7 +63,7 @@ I tested it successfully under Windows 7 (64bit) and 10 (64bit) with Matlab R201 ### Usage -1. By executing the start scripts (see above) +1. By executing the start script (see above) 2. Simply type `BLOCHUS` on the MATLABTM prompt (make sure the `blochus` folder is on the MATLABTM path) 3. Check the example scripts for the usage of the core functions without the GUI (inside the `scripts` folder) @@ -87,14 +87,15 @@ In no particular order and without guarantee that it will ever happen :-) : ### Cite as If you use BLOCHUS for your research, please cite it as: -Thomas Hiller. (2020, September 25). ThoHiller/nmr-blochus: v0.1.4 (Version v0.1.4). Zenodo. [https://doi.org/10.5281/zenodo.3829464] +Thomas Hiller. (2020, September 25). ThoHiller/nmr-blochus: v0.1.5 (Version v0.1.5). Zenodo. [https://doi.org/10.5281/zenodo.3829464] Note: Even though the version number might change due to updates, this DOI is permanent (represents all versions) and always links to the latest version. ### References - -1. Hiller, T., Dlugosch, R. and Müller-Petke, M., "Utilizing pre-polarization to enhance SNMR signals - effect of imperfect switch-off", Geophysical Journal International **222**(2), p.815-826, 2020, [DOI](https://doi.org/10.1093/gji/ggaa216) +1. Hiller, T., Costabel, S., Dlugosch, R. and Müller-Petke, M. "First Measurements of Surface Nuclear Magnetic Resonance Signals Without an Oscillating Excitation Pulse – Exploiting Non-Adiabatic Prepolarization Switch-Off", Geophysical Research Letters, **48**(23), e2021GL095371, 2021, [DOI](https://doi.org/10.1029/2021GL095371) +1. Hiller, T., Costabel, S., Radic, T., Dlugosch, R. and Müller-Petke, M. "Feasibility study on prepolarized surface nuclear magnetic resonance for soil moisture measurements", Vadose Zone Journal, **20**(5), e20138, 2021, [DOI](https://doi.org/10.1002/vzj2.20138) +2. Hiller, T., Dlugosch, R. and Müller-Petke, M., "Utilizing pre-polarization to enhance SNMR signals - effect of imperfect switch-off", Geophysical Journal International **222**(2), p.815-826, 2020, [DOI](https://doi.org/10.1093/gji/ggaa216) - - -

MATLAB is a registered trademark of The Mathworks, Inc.

diff --git a/doc/blochus/BLOCHUS/BLOCHUS.html b/doc/blochus/BLOCHUS/BLOCHUS.html index d568020..ace6cdc 100644 --- a/doc/blochus/BLOCHUS/BLOCHUS.html +++ b/doc/blochus/BLOCHUS/BLOCHUS.html @@ -105,10 +105,10 @@

SOURCE CODE ^%------------- BEGIN CODE -------------- 0034 0035 %% GUI 'header' info and default GUI settings -0036 myui.version = '0.1.4'; -0037 myui.date = '25.09.2020'; +0036 myui.version = '0.1.5'; +0037 myui.date = '30.06.2022'; 0038 myui.author = 'Thomas Hiller'; -0039 myui.email = 'thomas.hiller[at]leibniz-liag.de'; +0039 myui.email = 'thomas.hiller[at]bgr.de'; 0040 0041 myui.fontsize = 9; 0042 myui.axfontsize = 11; diff --git a/doc/blochus/BLOCHUS/BLOCHUS_loadDefaults.html b/doc/blochus/BLOCHUS/BLOCHUS_loadDefaults.html index 5ee46e9..260cef1 100644 --- a/doc/blochus/BLOCHUS/BLOCHUS_loadDefaults.html +++ b/doc/blochus/BLOCHUS/BLOCHUS_loadDefaults.html @@ -170,15 +170,15 @@

SOURCE CODE ^% initial magnetization [A/m] 0101 init.Minit = [1 0 0]; -0102 % Earth's magnetic field [T] default: 48µT min: 1fT max: 1T +0102 % Earth's magnetic field [T] default: 48µT min: 1fT max: 1T 0103 init.B0 = [48/1e6 1e-12 1]; 0104 % corresponding Larmor frequencies [Hz] 0105 init.Omega0 = [getOmega0(init.gamma,init.B0(1))/2/pi -1e9 1e9]; -0106 % T1 relaxation time [ms] default: 100ms min: 1µs max: 1000s +0106 % T1 relaxation time [ms] default: 100ms min: 1µs max: 1000s 0107 init.T1relax = [100 1e-3 1e6]; -0108 % T2 relaxation time [ms] default: 50ms min: 1µs max: 1000s +0108 % T2 relaxation time [ms] default: 50ms min: 1µs max: 1000s 0109 init.T2relax = [50 1e-3 1e6]; -0110 % simulation time [ms] default: 50ms min: 1µs max: 1000s +0110 % simulation time [ms] default: 50ms min: 1µs max: 1000s 0111 init.Tsim = [50 1e-3 1e6]; 0112 0113 % --- PRE-POLARIZATION --- @@ -189,14 +189,14 @@

SOURCE CODE ^% pre-polarization B-field amplitude in units of [B0] 0119 init.PrePolFactor = [100 1e-6 1e6]; 0120 % polar angle of pre-polarization field [deg] -0121 init.PrePolTheta = [90 0 360]; +0121 init.PrePolTheta = [90 1e-3 360]; 0122 % azimuthal angle of pre-polarization field [deg] 0123 init.PrePolPhi = [0 0 360]; 0124 % pre-polarization B-field switch amplitude in units of [B0] 0125 init.PrePolSwitchFactor = [1 1e-6 1e6]; -0126 % switch-off ramp time [ms] default: 1ms min: 1µs max: 1000s +0126 % switch-off ramp time [ms] default: 1ms min: 1µs max: 1000s 0127 init.PrePolTramp = [1 1e-6 1e6]; -0128 % switch-off ramp slope time [ms] default: 0.1ms min: 1µs max: 1000s +0128 % switch-off ramp slope time [ms] default: 0.1ms min: 1µs max: 1000s 0129 init.PrePolTslope = [0.1 1e-6 1e6]; 0130 0131 % --- PULSE --- @@ -208,7 +208,7 @@

SOURCE CODE ^'circular'; 0138 % relaxation during pulse [0/1] 0139 init.PulseRDP = 0; -0140 % pulse length [ms] default: 20ms min: 1µs max: 1000s +0140 % pulse length [ms] default: 20ms min: 1µs max: 1000s 0141 init.PulseTtau = [20 1e-6 1e6]; 0142 % pulse amplitude in units of [B0] default: 0.0061163 min: 1e-12 max: 1e6 0143 init.PulseB1Factor = [0.0061163 1e-12 1e6]; diff --git a/doc/blochus/functions/blochsim/fcn_BLOCHUS_ode.html b/doc/blochus/functions/blochsim/fcn_BLOCHUS_ode.html index 20ecc7b..b295981 100644 --- a/doc/blochus/functions/blochsim/fcn_BLOCHUS_ode.html +++ b/doc/blochus/functions/blochsim/fcn_BLOCHUS_ode.html @@ -140,7 +140,7 @@

SOURCE CODE ^switch param.type 0059 case 'std' -0060 B = B0*zunit; +0060 B = B0*zunit; 0061 % dM/dt 0062 dM = gamma*cross(m,B) - ( (m(1)*xunit+m(2)*yunit) / T2 ) - ( ((m(3)-M0(3))*zunit) / T1 ); 0063 diff --git a/doc/blochus/functions/blochsim/getMIDI_Tx.html b/doc/blochus/functions/blochsim/getMIDI_Tx.html index fdf113d..2f98cdc 100644 --- a/doc/blochus/functions/blochsim/getMIDI_Tx.html +++ b/doc/blochus/functions/blochsim/getMIDI_Tx.html @@ -317,7 +317,7 @@

SOURCE CODE ^% now create an artificial y-component by shifting the original signal -0249 % (x-component) 90° (quarter of a full period) +0249 % (x-component) 90° (quarter of a full period) 0250 % get the length of a quarter period 0251 t90 = dtL/4; 0252 % cut out these first samples and put them to the end of the original signal @@ -330,7 +330,7 @@

SOURCE CODE ^% pulse this is already taken care of during assembling) 0260 if strcmp(param.Tx,'MIDI_OR') && isfield(param,'PulseAxis') 0261 % depending on the frequency, there is a different amount of time -0262 % steps for the "90° phase shift" +0262 % steps for the "90° phase shift" 0263 shiftind = sum(t<t90); 0264 % now shift both components by the necessary amount of samples 0265 switch param.PulseAxis diff --git a/doc/blochus/functions/blochsim/getPulseTimeSeries.html b/doc/blochus/functions/blochsim/getPulseTimeSeries.html index 81ffcac..dfc8ef1 100644 --- a/doc/blochus/functions/blochsim/getPulseTimeSeries.html +++ b/doc/blochus/functions/blochsim/getPulseTimeSeries.html @@ -214,50 +214,53 @@

SOURCE CODE ^if Imod.useQ && Imod.Q > 0 0125 % get line (band) width -> f_L/Q (simple bandwidth for bandpass) 0126 Lwidth = abs(param.omega0/2/pi) / Imod.Q; -0127 % apply "Breit-Wigner" formula (here already normalized to 1 by -0128 % multiplying with pi*Lwidth) -0129 L = 1 ./ ( ((df+Imod.Qdf).^2 ./ Lwidth.^2) + 1 ); -0130 I = I.*L; -0131 end -0132 -0133 % get instantaneous phase due to frequency modulation -0134 theta = getPulsePhase(t,df,param,0); -0135 -0136 switch PulsePolarization -0137 case 'circular' -0138 Bout = I.*Amp.*[cos(theta + phi + phi_ax + phi_ref) ... -0139 sin(theta + phi + phi_ax + phi_ref)]; -0140 case 'linear' -0141 Bout = 2.*I.*Amp.*[cos(theta + phi + phi_ax + phi_ref) ... -0142 0.*t]; -0143 end -0144 end -0145 % if B-field values are NaN, set them to zero -0146 % this can happen due to the interpolation at the end of the MIDI-pulses -0147 Bout(isnan(Bout)) = 0; -0148 -0149 end -0150 -0151 %------------- END OF CODE -------------- -0152 -0153 %% License: -0154 % GNU GPLv3 -0155 % -0156 % BLOCHUS -0157 % Copyright (C) 2019 Thomas Hiller +0127 % apply Cauchy-Lorentz type formula (here already normalized to +0128 % 1 by multiplying with pi*Lwidth) +0129 % this is basically a Cauchy distribution PDF of the form: +0130 % PDF = 1/pi * ( bw / ((f-f0)^2 + bw^2) ) +0131 % tweaked with some algebra +0132 L = 1 ./ ( ((df+Imod.Qdf).^2 ./ Lwidth.^2) + 1 ); +0133 I = I.*L; +0134 end +0135 +0136 % get instantaneous phase due to frequency modulation +0137 theta = getPulsePhase(t,df,param,0); +0138 +0139 switch PulsePolarization +0140 case 'circular' +0141 Bout = I.*Amp.*[cos(theta + phi + phi_ax + phi_ref) ... +0142 sin(theta + phi + phi_ax + phi_ref)]; +0143 case 'linear' +0144 Bout = 2.*I.*Amp.*[cos(theta + phi + phi_ax + phi_ref) ... +0145 0.*t]; +0146 end +0147 end +0148 % if B-field values are NaN, set them to zero +0149 % this can happen due to the interpolation at the end of the MIDI-pulses +0150 Bout(isnan(Bout)) = 0; +0151 +0152 end +0153 +0154 %------------- END OF CODE -------------- +0155 +0156 %% License: +0157 % GNU GPLv3 0158 % -0159 % This program is free software: you can redistribute it and/or modify -0160 % it under the terms of the GNU General Public License as published by -0161 % the Free Software Foundation, either version 3 of the License, or -0162 % (at your option) any later version. -0163 % -0164 % This program is distributed in the hope that it will be useful, -0165 % but WITHOUT ANY WARRANTY; without even the implied warranty of -0166 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -0167 % GNU General Public License for more details. -0168 % -0169 % You should have received a copy of the GNU General Public License -0170 % along with this program. If not, see <https://www.gnu.org/licenses/>. +0159 % BLOCHUS +0160 % Copyright (C) 2019 Thomas Hiller +0161 % +0162 % This program is free software: you can redistribute it and/or modify +0163 % it under the terms of the GNU General Public License as published by +0164 % the Free Software Foundation, either version 3 of the License, or +0165 % (at your option) any later version. +0166 % +0167 % This program is distributed in the hope that it will be useful, +0168 % but WITHOUT ANY WARRANTY; without even the implied warranty of +0169 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +0170 % GNU General Public License for more details. +0171 % +0172 % You should have received a copy of the GNU General Public License +0173 % along with this program. If not, see <https://www.gnu.org/licenses/>.
Generated by m2html © 2005
\ No newline at end of file diff --git a/doc/blochus/functions/blochsim/getRampAmplitude.html b/doc/blochus/functions/blochsim/getRampAmplitude.html index c337653..f629ec8 100644 --- a/doc/blochus/functions/blochsim/getRampAmplitude.html +++ b/doc/blochus/functions/blochsim/getRampAmplitude.html @@ -140,77 +140,80 @@

SOURCE CODE ^case 'lin' % linear 0060 Bp = Bmax.*(1-t./Tramp); -0061 end -0062 -0063 return -0064 % --- +0061 case 'custom' +0062 idata = param.interp; +0063 Bp = Bmax.*interp1(idata.t,idata.I,t); +0064 end 0065 -0066 % --- -0067 function Bp = getLinExpAmp(Bmax,Bstar,T,t) -0068 % linear + exponential ramp after: -0069 % Conradi et al., 2017, Journal of Magnetic Resonance 281, p.241-245 -0070 % https://doi.org/10.1016/j.jmr.2017.06.001 -0071 -0072 if numel(t)>1 -0073 % linear part -0074 Bplin = (-Bmax/T)*t + Bmax; -0075 % exponential part -0076 Bpexp = exp(-t /(Bstar*T/Bmax)); -0077 % find change -0078 index = find(abs(Bplin-Bstar)==min(abs(Bplin-Bstar)),1,'first'); -0079 % merge the lin- and exp-part and scale the amplitude of the exp-part -0080 % to that of the lin-part at the switch-over time t(index) -0081 scale_point = Bplin(index)/Bpexp(index); -0082 % in case something goes south due to very small numbers set the -0083 % amplitude to 0 -0084 if isinf(scale_point) || isnan(scale_point) -0085 scale_point = 0; -0086 end -0087 % the final amplitude vector -0088 Bp = [Bplin(1:index-1); scale_point * Bpexp(index:end)]; -0089 else -0090 % linear part -0091 Bplin = (-Bmax/T)*t + Bmax; -0092 % exponential part -0093 Bpexp = exp(-t /(Bstar*T/Bmax)); -0094 % Bstar time tstar -0095 tstar = (Bstar-Bmax)/(-Bmax/T); -0096 % amplitude at tstar for scaling -0097 Btstar = exp(-tstar /(Bstar*T/Bmax)); -0098 % apply -0099 if t<tstar -0100 Bp = Bplin; -0101 else -0102 Bp = (Bstar/Btstar) * Bpexp; -0103 end -0104 end -0105 -0106 % if due to division by "0" the value is NaN ... set it to 0 -0107 Bp(isnan(Bp)) = 0; +0066 return +0067 % --- +0068 +0069 % --- +0070 function Bp = getLinExpAmp(Bmax,Bstar,T,t) +0071 % linear + exponential ramp after: +0072 % Conradi et al., 2017, Journal of Magnetic Resonance 281, p.241-245 +0073 % https://doi.org/10.1016/j.jmr.2017.06.001 +0074 +0075 if numel(t)>1 +0076 % linear part +0077 Bplin = (-Bmax/T)*t + Bmax; +0078 % exponential part +0079 Bpexp = exp(-t /(Bstar*T/Bmax)); +0080 % find change +0081 index = find(abs(Bplin-Bstar)==min(abs(Bplin-Bstar)),1,'first'); +0082 % merge the lin- and exp-part and scale the amplitude of the exp-part +0083 % to that of the lin-part at the switch-over time t(index) +0084 scale_point = Bplin(index)/Bpexp(index); +0085 % in case something goes south due to very small numbers set the +0086 % amplitude to 0 +0087 if isinf(scale_point) || isnan(scale_point) +0088 scale_point = 0; +0089 end +0090 % the final amplitude vector +0091 Bp = [Bplin(1:index-1); scale_point * Bpexp(index:end)]; +0092 else +0093 % linear part +0094 Bplin = (-Bmax/T)*t + Bmax; +0095 % exponential part +0096 Bpexp = exp(-t /(Bstar*T/Bmax)); +0097 % Bstar time tstar +0098 tstar = (Bstar-Bmax)/(-Bmax/T); +0099 % amplitude at tstar for scaling +0100 Btstar = exp(-tstar /(Bstar*T/Bmax)); +0101 % apply +0102 if t<tstar +0103 Bp = Bplin; +0104 else +0105 Bp = (Bstar/Btstar) * Bpexp; +0106 end +0107 end 0108 -0109 -0110 return +0109 % if due to division by "0" the value is NaN ... set it to 0 +0110 Bp(isnan(Bp)) = 0; 0111 -0112 %------------- END OF CODE -------------- -0113 -0114 %% License: -0115 % GNU GPLv3 -0116 % -0117 % BLOCHUS -0118 % Copyright (C) 2019 Thomas Hiller +0112 +0113 return +0114 +0115 %------------- END OF CODE -------------- +0116 +0117 %% License: +0118 % GNU GPLv3 0119 % -0120 % This program is free software: you can redistribute it and/or modify -0121 % it under the terms of the GNU General Public License as published by -0122 % the Free Software Foundation, either version 3 of the License, or -0123 % (at your option) any later version. -0124 % -0125 % This program is distributed in the hope that it will be useful, -0126 % but WITHOUT ANY WARRANTY; without even the implied warranty of -0127 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -0128 % GNU General Public License for more details. -0129 % -0130 % You should have received a copy of the GNU General Public License -0131 % along with this program. If not, see <https://www.gnu.org/licenses/>. +0120 % BLOCHUS +0121 % Copyright (C) 2019 Thomas Hiller +0122 % +0123 % This program is free software: you can redistribute it and/or modify +0124 % it under the terms of the GNU General Public License as published by +0125 % the Free Software Foundation, either version 3 of the License, or +0126 % (at your option) any later version. +0127 % +0128 % This program is distributed in the hope that it will be useful, +0129 % but WITHOUT ANY WARRANTY; without even the implied warranty of +0130 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +0131 % GNU General Public License for more details. +0132 % +0133 % You should have received a copy of the GNU General Public License +0134 % along with this program. If not, see <https://www.gnu.org/licenses/>.
Generated by m2html © 2005
\ No newline at end of file diff --git a/doc/blochus/functions/blochsim/getRotationMatrixFromAngleandAxis.html b/doc/blochus/functions/blochsim/getRotationMatrixFromAngleandAxis.html index bb6b44e..518e565 100644 --- a/doc/blochus/functions/blochsim/getRotationMatrixFromAngleandAxis.html +++ b/doc/blochus/functions/blochsim/getRotationMatrixFromAngleandAxis.html @@ -33,11 +33,11 @@

DESCRIPTION ^SOURCE CODE ^% getRotationMatrixFromAngleandAxis(phi,n) 0007 % 0008 % Inputs: -0009 % phi - rotation angle [rad] -0010 % n - rotation axis vector [x y z] +0009 % phi - rotation angle [rad]; size Nx1 +0010 % n - rotation axis vector [x y z]; size Nx3 0011 % 0012 % Outputs: -0013 % R - 3x3 rotation matrix +0013 % R - 3x3xN rotation matrix 0014 % 0015 % Example: 0016 % R = getRotationMatrixFromAngleandAxis(pi,[0 0 1]') @@ -110,44 +110,80 @@

SOURCE CODE ^%------------- BEGIN CODE -------------- 0037 -0038 % make "n" a unit vector -0039 n = n./norm(n); -0040 % get the individual components -0041 nx = n(1); -0042 ny = n(2); -0043 nz = n(3); -0044 % matrix terms needed -0045 omcos = 1-cos(phi); -0046 cosp = cos(phi); -0047 sinp = sin(phi); -0048 -0049 % assemble rotation matrix R -0050 R = [nx*nx*omcos + cosp nx*ny*omcos - nz*sinp nx*nz*omcos + ny*sinp; ... -0051 ny*nx*omcos + nz*sinp ny*ny*omcos + cosp ny*nz*omcos - nx*sinp; ... -0052 nz*nx*omcos - ny*sinp nz*ny*omcos + nx*sinp nz*nz*omcos + cosp]; -0053 -0054 return -0055 -0056 %------------- END OF CODE -------------- -0057 -0058 %% License: -0059 % GNU GPLv3 -0060 % -0061 % BLOCHUS -0062 % Copyright (C) 2019 Thomas Hiller -0063 % -0064 % This program is free software: you can redistribute it and/or modify -0065 % it under the terms of the GNU General Public License as published by -0066 % the Free Software Foundation, either version 3 of the License, or -0067 % (at your option) any later version. -0068 % -0069 % This program is distributed in the hope that it will be useful, -0070 % but WITHOUT ANY WARRANTY; without even the implied warranty of -0071 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -0072 % GNU General Public License for more details. -0073 % -0074 % You should have received a copy of the GNU General Public License -0075 % along with this program. If not, see <https://www.gnu.org/licenses/>. +0038 % for only one axis-angle pair +0039 if numel(phi) == 1 +0040 % make "n" a unit vector +0041 n = n./norm(n); +0042 % get the individual components +0043 nx = n(1); +0044 ny = n(2); +0045 nz = n(3); +0046 % matrix terms needed +0047 omcos = 1-cos(phi); +0048 cosp = cos(phi); +0049 sinp = sin(phi); +0050 +0051 % assemble rotation matrix R +0052 R(1,1) = nx*nx*omcos + cosp; +0053 R(1,2) = nx*ny*omcos - nz*sinp; +0054 R(1,3) = nx*nz*omcos + ny*sinp; +0055 +0056 R(2,1) = ny*nx*omcos + nz*sinp; +0057 R(2,2) = ny*ny*omcos + cosp; +0058 R(2,3) = ny*nz*omcos - nx*sinp; +0059 +0060 R(3,1) = nz*nx*omcos - ny*sinp; +0061 R(3,2) = nz*ny*omcos + nx*sinp; +0062 R(3,3) = nz*nz*omcos + cosp; +0063 +0064 else % for multiple axes and angles +0065 +0066 % n should contain only unit vectors! +0067 % get the individual components +0068 nx = n(:,1); +0069 ny = n(:,2); +0070 nz = n(:,3); +0071 % matrix terms needed +0072 omcos = 1-cos(phi); +0073 cosp = cos(phi); +0074 sinp = sin(phi); +0075 +0076 % assemble rotation matrix R +0077 R(1,1,:) = nx.*nx.*omcos + cosp; +0078 R(1,2,:) = nx.*ny.*omcos - nz.*sinp; +0079 R(1,3,:) = nx.*nz.*omcos + ny.*sinp; +0080 +0081 R(2,1,:) = ny.*nx.*omcos + nz.*sinp; +0082 R(2,2,:) = ny.*ny.*omcos + cosp; +0083 R(2,3,:) = ny.*nz.*omcos - nx.*sinp; +0084 +0085 R(3,1,:) = nz.*nx.*omcos - ny.*sinp; +0086 R(3,2,:) = nz.*ny.*omcos + nx.*sinp; +0087 R(3,3,:) = nz.*nz.*omcos + cosp; +0088 end +0089 +0090 return +0091 +0092 %------------- END OF CODE -------------- +0093 +0094 %% License: +0095 % GNU GPLv3 +0096 % +0097 % BLOCHUS +0098 % Copyright (C) 2019 Thomas Hiller +0099 % +0100 % This program is free software: you can redistribute it and/or modify +0101 % it under the terms of the GNU General Public License as published by +0102 % the Free Software Foundation, either version 3 of the License, or +0103 % (at your option) any later version. +0104 % +0105 % This program is distributed in the hope that it will be useful, +0106 % but WITHOUT ANY WARRANTY; without even the implied warranty of +0107 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +0108 % GNU General Public License for more details. +0109 % +0110 % You should have received a copy of the GNU General Public License +0111 % along with this program. If not, see <https://www.gnu.org/licenses/>.
Generated by m2html © 2005
\ No newline at end of file diff --git a/doc/blochus/functions/blochsim/getRotationMatrixFromVectors.html b/doc/blochus/functions/blochsim/getRotationMatrixFromVectors.html index 87c87e3..ca3cb50 100644 --- a/doc/blochus/functions/blochsim/getRotationMatrixFromVectors.html +++ b/doc/blochus/functions/blochsim/getRotationMatrixFromVectors.html @@ -113,46 +113,48 @@

SOURCE CODE ^% cross product of both vectors 0041 c = cross(A,B); -0042 % check if A and B are parallel / antiparallel -0043 if sum(c)==0 -0044 % check if A == B (parallel) -0045 if all(A==B) -0046 % in that case the rotation matrix is obviously identity -0047 R = eye(3); -0048 else % A == -B (antiparallel) -0049 R = -eye(3); -0050 end -0051 else -0052 % skew-symmetric cross-product -0053 ssc = [ 0 -c(3) c(2); -0054 c(3) 0 -c(1); -0055 -c(2) c(1) 0 ]; -0056 % rotation matrix R -0057 R = eye(3) + ssc + ssc^2*(1-dot(A,B))/(norm(c))^2; -0058 end -0059 -0060 return +0042 % angle between A and B +0043 alpha = acos(dot(A,B)); +0044 % check if A and B are parallel / antiparallel +0045 if abs(sum(c))<1e-128 && (alpha==0 || alpha==pi) +0046 % check if A == B (alpha=0 -> parallel) +0047 if alpha==0 +0048 % in that case the rotation matrix is obviously identity +0049 R = eye(3); +0050 else % A == -B (antiparallel) +0051 R = -eye(3); +0052 end +0053 else +0054 % skew-symmetric cross-product +0055 ssc = [ 0 -c(3) c(2); +0056 c(3) 0 -c(1); +0057 -c(2) c(1) 0 ]; +0058 % rotation matrix R +0059 R = eye(3) + ssc + (ssc/norm(c))^2*(1-dot(A,B)); +0060 end 0061 -0062 %------------- END OF CODE -------------- +0062 return 0063 -0064 %% License: -0065 % GNU GPLv3 -0066 % -0067 % BLOCHUS -0068 % Copyright (C) 2019 Thomas Hiller -0069 % -0070 % This program is free software: you can redistribute it and/or modify -0071 % it under the terms of the GNU General Public License as published by -0072 % the Free Software Foundation, either version 3 of the License, or -0073 % (at your option) any later version. -0074 % -0075 % This program is distributed in the hope that it will be useful, -0076 % but WITHOUT ANY WARRANTY; without even the implied warranty of -0077 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -0078 % GNU General Public License for more details. -0079 % -0080 % You should have received a copy of the GNU General Public License -0081 % along with this program. If not, see <https://www.gnu.org/licenses/>. +0064 %------------- END OF CODE -------------- +0065 +0066 %% License: +0067 % GNU GPLv3 +0068 % +0069 % BLOCHUS +0070 % Copyright (C) 2019 Thomas Hiller +0071 % +0072 % This program is free software: you can redistribute it and/or modify +0073 % it under the terms of the GNU General Public License as published by +0074 % the Free Software Foundation, either version 3 of the License, or +0075 % (at your option) any later version. +0076 % +0077 % This program is distributed in the hope that it will be useful, +0078 % but WITHOUT ANY WARRANTY; without even the implied warranty of +0079 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +0080 % GNU General Public License for more details. +0081 % +0082 % You should have received a copy of the GNU General Public License +0083 % along with this program. If not, see <https://www.gnu.org/licenses/>.
Generated by m2html © 2005
\ No newline at end of file diff --git a/doc/blochus/functions/interface/getRampParameters.html b/doc/blochus/functions/interface/getRampParameters.html index 45e3f83..10984c0 100644 --- a/doc/blochus/functions/interface/getRampParameters.html +++ b/doc/blochus/functions/interface/getRampParameters.html @@ -132,13 +132,13 @@

SOURCE CODE ^% now adjust the direction of the pre-polarization field 0063 % 1.) rotation by angle theta [deg] around y-axis for the z-unit vector -0064 % this means an angle of 90° around the y-axis will turn the z-unit into +0064 % this means an angle of 90° around the y-axis will turn the z-unit into 0065 % the x-unit vector 0066 RM = getRotationMatrixFromAngleandAxis(deg2rad(data.prepol.Theta),[0 1 0]); 0067 orient = RM * zunit; 0068 % now this new orientation vector gets rotated by angle phi [deg] around 0069 % the z-axis for the x-unit vector -0070 % this means an angle of 90° around the z-axis will turn the x-unit into +0070 % this means an angle of 90° around the z-axis will turn the x-unit into 0071 % the y-unit vector 0072 RM = getRotationMatrixFromAngleandAxis(deg2rad(data.prepol.Phi),[0 0 1]); 0073 orient = RM * orient; diff --git a/functions/blochsim/fcn_BLOCHUS_ode.m b/functions/blochsim/fcn_BLOCHUS_ode.m index 04e557e..b39f6f4 100644 --- a/functions/blochsim/fcn_BLOCHUS_ode.m +++ b/functions/blochsim/fcn_BLOCHUS_ode.m @@ -57,7 +57,7 @@ switch param.type case 'std' - B = B0*zunit; + B = B0*zunit; % dM/dt dM = gamma*cross(m,B) - ( (m(1)*xunit+m(2)*yunit) / T2 ) - ( ((m(3)-M0(3))*zunit) / T1 ); diff --git a/functions/blochsim/getPulseTimeSeries.m b/functions/blochsim/getPulseTimeSeries.m index 1739567..e63441d 100644 --- a/functions/blochsim/getPulseTimeSeries.m +++ b/functions/blochsim/getPulseTimeSeries.m @@ -124,8 +124,11 @@ if Imod.useQ && Imod.Q > 0 % get line (band) width -> f_L/Q (simple bandwidth for bandpass) Lwidth = abs(param.omega0/2/pi) / Imod.Q; - % apply "Breit-Wigner" formula (here already normalized to 1 by - % multiplying with pi*Lwidth) + % apply Cauchy-Lorentz type formula (here already normalized to + % 1 by multiplying with pi*Lwidth) + % this is basically a Cauchy distribution PDF of the form: + % PDF = 1/pi * ( bw / ((f-f0)^2 + bw^2) ) + % tweaked with some algebra L = 1 ./ ( ((df+Imod.Qdf).^2 ./ Lwidth.^2) + 1 ); I = I.*L; end diff --git a/functions/blochsim/getRampAmplitude.m b/functions/blochsim/getRampAmplitude.m index 7ba0ecf..29acc85 100644 --- a/functions/blochsim/getRampAmplitude.m +++ b/functions/blochsim/getRampAmplitude.m @@ -58,6 +58,9 @@ Bp = Bmax .* (0.5+(cos(pi*t./Tramp)./2)); case 'lin' % linear Bp = Bmax.*(1-t./Tramp); + case 'custom' + idata = param.interp; + Bp = Bmax.*interp1(idata.t,idata.I,t); end return diff --git a/functions/blochsim/getRotationMatrixFromAngleandAxis.m b/functions/blochsim/getRotationMatrixFromAngleandAxis.m index 793a597..97700d2 100644 --- a/functions/blochsim/getRotationMatrixFromAngleandAxis.m +++ b/functions/blochsim/getRotationMatrixFromAngleandAxis.m @@ -6,11 +6,11 @@ % getRotationMatrixFromAngleandAxis(phi,n) % % Inputs: -% phi - rotation angle [rad] -% n - rotation axis vector [x y z] +% phi - rotation angle [rad]; size Nx1 +% n - rotation axis vector [x y z]; size Nx3 % % Outputs: -% R - 3x3 rotation matrix +% R - 3x3xN rotation matrix % % Example: % R = getRotationMatrixFromAngleandAxis(pi,[0 0 1]') @@ -35,21 +35,57 @@ %------------- BEGIN CODE -------------- -% make "n" a unit vector -n = n./norm(n); -% get the individual components -nx = n(1); -ny = n(2); -nz = n(3); -% matrix terms needed -omcos = 1-cos(phi); -cosp = cos(phi); -sinp = sin(phi); +% for only one axis-angle pair +if numel(phi) == 1 + % make "n" a unit vector + n = n./norm(n); + % get the individual components + nx = n(1); + ny = n(2); + nz = n(3); + % matrix terms needed + omcos = 1-cos(phi); + cosp = cos(phi); + sinp = sin(phi); + + % assemble rotation matrix R + R(1,1) = nx*nx*omcos + cosp; + R(1,2) = nx*ny*omcos - nz*sinp; + R(1,3) = nx*nz*omcos + ny*sinp; + + R(2,1) = ny*nx*omcos + nz*sinp; + R(2,2) = ny*ny*omcos + cosp; + R(2,3) = ny*nz*omcos - nx*sinp; + + R(3,1) = nz*nx*omcos - ny*sinp; + R(3,2) = nz*ny*omcos + nx*sinp; + R(3,3) = nz*nz*omcos + cosp; -% assemble rotation matrix R -R = [nx*nx*omcos + cosp nx*ny*omcos - nz*sinp nx*nz*omcos + ny*sinp; ... - ny*nx*omcos + nz*sinp ny*ny*omcos + cosp ny*nz*omcos - nx*sinp; ... - nz*nx*omcos - ny*sinp nz*ny*omcos + nx*sinp nz*nz*omcos + cosp]; +else % for multiple axes and angles + + % n should contain only unit vectors! + % get the individual components + nx = n(:,1); + ny = n(:,2); + nz = n(:,3); + % matrix terms needed + omcos = 1-cos(phi); + cosp = cos(phi); + sinp = sin(phi); + + % assemble rotation matrix R + R(1,1,:) = nx.*nx.*omcos + cosp; + R(1,2,:) = nx.*ny.*omcos - nz.*sinp; + R(1,3,:) = nx.*nz.*omcos + ny.*sinp; + + R(2,1,:) = ny.*nx.*omcos + nz.*sinp; + R(2,2,:) = ny.*ny.*omcos + cosp; + R(2,3,:) = ny.*nz.*omcos - nx.*sinp; + + R(3,1,:) = nz.*nx.*omcos - ny.*sinp; + R(3,2,:) = nz.*ny.*omcos + nx.*sinp; + R(3,3,:) = nz.*nz.*omcos + cosp; +end return diff --git a/functions/blochsim/getRotationMatrixFromVectors.m b/functions/blochsim/getRotationMatrixFromVectors.m index d55f399..abb7d2d 100644 --- a/functions/blochsim/getRotationMatrixFromVectors.m +++ b/functions/blochsim/getRotationMatrixFromVectors.m @@ -39,10 +39,12 @@ B = B(:)./norm(B); % cross product of both vectors c = cross(A,B); +% angle between A and B +alpha = acos(dot(A,B)); % check if A and B are parallel / antiparallel -if sum(c)==0 - % check if A == B (parallel) - if all(A==B) +if abs(sum(c))<1e-128 && (alpha==0 || alpha==pi) + % check if A == B (alpha=0 -> parallel) + if alpha==0 % in that case the rotation matrix is obviously identity R = eye(3); else % A == -B (antiparallel) @@ -54,7 +56,7 @@ c(3) 0 -c(1); -c(2) c(1) 0 ]; % rotation matrix R - R = eye(3) + ssc + ssc^2*(1-dot(A,B))/(norm(c))^2; + R = eye(3) + ssc + (ssc/norm(c))^2*(1-dot(A,B)); end return diff --git a/pyBLOCHUS/README.md b/pyBLOCHUS/README.md index 3646d1f..b0d82cb 100644 --- a/pyBLOCHUS/README.md +++ b/pyBLOCHUS/README.md @@ -36,11 +36,11 @@ I tested it successfully under Windows 7 (64bit) and 10 (64bit) with [Anaconda]( 1. Makes sure that all required packages are installed and work properly 2. Download or `clone` the **BLOCHUS** repository and put the folder somewhere on your hard drive. -3. As I recommend uying an IDE like [Spyder](https://www.spyder-ide.org/), make sure that the **pyBLOCHUS** folder is on your PYTHONPATH (e.g. in [Spyder](https://www.spyder-ide.org/): Tools -> PYTHONPATH manager -> Add path) +3. As I recommend using an IDE like [Spyder](https://www.spyder-ide.org/), make sure that the **pyBLOCHUS** folder is on your PYTHONPATH (e.g. in [Spyder](https://www.spyder-ide.org/): Tools -> PYTHONPATH manager -> Add path) 4. Sometimes a [Spyder](https://www.spyder-ide.org/) restart is necessary tu update/sync the PYTHONPATH - - - ### Usage -Have a look in the [example](pyBLOCHUS/examples) scripts for the usage of the different modules. +Have a look in the [example](examples) scripts for the usage of the different modules. diff --git a/pyBLOCHUS/pyBLOCHUS/blochus_pulse.py b/pyBLOCHUS/pyBLOCHUS/blochus_pulse.py index e7b3a31..56b0443 100644 --- a/pyBLOCHUS/pyBLOCHUS/blochus_pulse.py +++ b/pyBLOCHUS/pyBLOCHUS/blochus_pulse.py @@ -484,8 +484,11 @@ def get_pulse_amplitude(self, time): # get line (band) width # -> f_L/Q (simple bandwidth for bandpass) lwidth = abs(self.larmor_f)/self.imod.qual_fac - # apply "Breit-Wigner" formula (here already normalized to 1 by - # multiplying with pi*lwidth) + # apply Cauchy-Lorentz type formula (here already normalized to + # 1 by multiplying with pi*lwidth) + # this is basically a Cauchy distribution PDF of the form: + # PDF = 1/pi * ( bw / ((f-f0)^2 + bw^2) ) + # tweaked with some algebra ltmp = 1.0 / (((freq_offset+self.imod.qual_fac_df)**2.0 / lwidth**2)+1.0) current = current*ltmp diff --git a/scripts/example3.m b/scripts/example3.m index 2ce93a4..8a5e4cd 100644 --- a/scripts/example3.m +++ b/scripts/example3.m @@ -245,7 +245,7 @@ for b = 1:numel(Bstar) for r = 1:numel(theta) - % initial orientation is towards y-axis + % initial orientation is towards x-axis orient = getRotationMatrixFromAngleandAxis(deg2rad(theta(r)),[0 1 0])*zunit; % initial magnetization in the direction of B0+Bp Morient = orient*rampparam.Bmax + odeparam.B0*zunit; diff --git a/scripts/example4.m b/scripts/example4.m index 3bf7306..5754d42 100644 --- a/scripts/example4.m +++ b/scripts/example4.m @@ -140,7 +140,7 @@ % switch-over amplitude for the "linexp" ramp (factor*B0) [T] rampparam.Bstar = SwitchFactor*odeparam.B0; - % initial orientation is towards y-axis + % initial orientation is towards x-axis orient = getRotationMatrixFromAngleandAxis(deg2rad(theta(tt)),[0 1 0])*zunit; % initial magnetization in the direction of B0+Bp Morient = orient*rampparam.Bmax + odeparam.B0*zunit; diff --git a/scripts/example6.m b/scripts/example6.m new file mode 100644 index 0000000..8eb30a5 --- /dev/null +++ b/scripts/example6.m @@ -0,0 +1,198 @@ +% script: example6.m +% This script demonstrates how to calculate a lookup-table for an arbitrarily +% shaped pre-polarization switch-off ramp with BLOCHUS. +% +% The switch-off ramp is imported as 'time' and 'current' data and used +% inside BLOCHUS via interpolation +% +% See also: BLOCHUS +% Author: Thomas Hiller +% email: thomas.hiller[at]leibniz-liag.de +% License: GNU GPLv3 (at end) + +%------------- BEGIN CODE -------------- + +% clear data +clear variables; clc; + +% because this script takes approx. 30min, ask the user before continuing +answer = questdlg({'This script needs about 30 min. to run.',... + 'Do you want to continue?'},'Continue?','Yes','No','Yes'); +% handle response +switch answer + case 'Yes' + % nothing to do + case {'No',''} + return +end + +% output switches +save_data = true; +plot_data = true; + +% --- LOOKUP TABLE SETTINGS ----------------------------------------------- +% global lookup table settings (change as you like): +% three custom ramps +lookup_name = {'ramp_2ms','ramp_3ms','ramp_3ms_raw'}; + +% IMPORTANT NOTE: +% This very coarse (and basically useless) Np=10 x Nt=10 grid takes for the +% 3 ramps roughly 30min! +% local lookup table settings (careful with the discretization) +Np = 10; +Nt = 10; +% pre-polarization factor [B0] +PPfac = logspace(-1,4,Np); +% angle between Bp and B0 [deg] +theta = linspace(0,180,Nt); +% ------------------------------------------------------------------------- + +% basic settings +zunit = [0 0 1]'; % z unit vector +% ODE solver error tolerance +tol = 1e-9; +% ODE solver options +options = odeset('RelTol',tol,'AbsTol',[tol tol tol]); + +% hydrogen proton +nucleus = '1H'; + +% parameter needed for the ODE solver +% simulation type [string] +odeparam.type = 'prepol'; +% equilibrium magnetization [A/m] +odeparam.M0 = zunit; +% primary (Earth's) magnetic field B0 [T] +odeparam.B0 = 4.8e-5; +% relaxation time T1 [s] +odeparam.T1 = 2; +% relaxation time T2 [s] +odeparam.T2 = 1; +% gyromagnetic ratio [rad/s/T] +odeparam.gamma = getGyroRatio(nucleus); + +% pre-polarization switch-off ramp parameter +% relaxation during switch-off [0/1] +odeparam.RDS = 0; + +% because Minit is not normalized to 1, M0 needs to be adjusted +odeparam.M0 = odeparam.B0*zunit; + +data = cell(1,numel(lookup_name)); +for i1 = 1:numel(lookup_name) + % load ramp file + dataR = load([lookup_name{i1},'.mat']); + dataR = dataR.data; + Tramp = dataR.t(end); + + % create the output data matrix + % adiabatic quality p + p = zeros(numel(theta),numel(PPfac)); + % final magnetization vector + Mfinal = zeros(numel(theta),numel(PPfac),3); + % loop over the local parameters + for tt = 1:numel(theta) + for pp = 1:numel(PPfac) + + % switch-off ramp type [string] + rampparam.ramp = 'custom'; + % gyromagnetic ratio [rad/s/T] + rampparam.gamma = odeparam.gamma; + % primary (Earth's) magnetic field B0 [T] + rampparam.B0 = odeparam.B0; + % amplitude of pre-polarization field [T] + rampparam.Bmax = PPfac(pp)*odeparam.B0; + % switch-off ramp time [s] + rampparam.Tramp = Tramp; + % ramp interpolation data + rampparam.interp.t = dataR.t; + rampparam.interp.I = dataR.I; + + % switch-off slope time is NOT NEEDED [s] + rampparam.Tslope = Tramp; + % switch-over amplitude is NOT NEEDED [T] + rampparam.Bstar = odeparam.B0; + + % initial orientation is towards x-axis + orient = getRotationMatrixFromAngleandAxis(deg2rad(theta(tt)),[0 1 0])*zunit; + % initial magnetization in the direction of B0+Bp + Morient = orient*rampparam.Bmax + odeparam.B0*zunit; + % Minit does not get normalized in this case + Minit = Morient; + % orientation of the pre-polarization pulse axis (x-axis) + odeparam.orient = orient; + + % add the ramp parameter to the ode parameter struct + odeparam.rampparam = rampparam; + + % ODE solver call + % OUTPUT: time T and magnetization M + [T,M] = ode45(@(t,m) fcn_BLOCHUS_ode(t,m,odeparam),rampparam.interp.t,Minit,options); + + % normalized magnetization vector at end of switch-off ramp + MMn = M(end,:)./norm(M(end,:)); + % "adiabatic quality" p of the switch-off ramp + % -> orientation of M with respect to z_unit + p(tt,pp) = dot(MMn(end,:),zunit)./norm(zunit); + % final orientation of M + Mfinal(tt,pp,1:3) = M(end,:); + + end + infostr1 = lookup_name{i1}; + infostr2 = ['theta=',num2str(theta(tt)),'°']; + disp([infostr1,' | ',infostr2]); + end + % save the local lookup table for this particular combination of + % ramp shape and ramp time + data{i1}.ramp = lookup_name{i1}; + data{i1}.PPfac = PPfac; + data{i1}.theta = theta; + data{i1}.odeparam = odeparam; + data{i1}.M = Mfinal; + data{i1}.p = p; +end + +% save lookup table +if save_data + save('example6_lookup_table.mat','data'); +end + +% plot data +if plot_data + [xx,yy] = meshgrid(linspace(0,180,Nt+1),logspace(-1,4,Np+1)); + figure; + for i1 = 1:numel(lookup_name) + p = data{i1}.p; + subplot(1,numel(lookup_name),i1); + surf(xx,yy,ones(size(xx)),p'); view(2); + shading flat; + set(gca,'XLim',[0 180],'XTick',linspace(0,180,5),'YScale','log',... + 'YLim',[0.1 1e4],'YTick',logspace(-1,4,6),'FontSize',12); + xlabel('\theta [deg]'); + ylabel('Bp [B0]'); + tstr = lookup_name{i1}; + tstr = strrep(tstr,'_',' '); + title(tstr); + end +end + +%------------- END OF CODE -------------- + +%% License: +% GNU GPLv3 +% +% BLOCHUS +% Copyright (C) 2019 Thomas Hiller +% +% This program is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% This program is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program. If not, see . diff --git a/scripts/example6_ref.png b/scripts/example6_ref.png new file mode 100644 index 0000000000000000000000000000000000000000..6daa3ff1d10f9975ad80b247c5926d5d70008149 GIT binary patch literal 54675 zcmeFa2{e}N_cnZ$3JsKK4wVo^MW&2RlCgv+q>`j4LS~PoQijYTR7x2uN#>$5Mv*B~ z=FIaveaGn;^#0cGTi^Tt*1O*IzU%+H*Lt2O-M9O?uJhdI-p4-nvG;RFMM-wq;&qED zidrUj_`oszhn1oh(9zAu-{>rQWq|+CUpcI4Mp1lQ$UpP0Kb5k;U(#NaQ#?r9vS>Y> z7{$*!`W}D9aqZxVYpN#auUYDu>Qjm`#}#DuZD&7zL7!bw@7e_o6H_Cz?d)>KXZgfW z)U3hhHW@uu}e%$b~KDAqT#}RYm9lLk#-p#&ir`R45v7I}8It2VEik*@> zAaneBaA%E`BXjdZsp&IHr=_>@EMCmCIgIJx0fsH6^z@r~4l^G-@I1*25>MRvMj}Dvta7>opx3eQx#NKTy<*gGL+w`RTw{E!uy6Sn;N8-rqk+ zFNiV!`_I&Jng8CuvWcs5Aw?~(^4oE;)_$UE<6evHy6LI!-iG>aLY2Zjp6Yu)R!C;rya2sppC3o12?+W9Lzh>~v1W#l<$=Nd|d7 z+YTN)m=IW(W-}5IYM*S`Q|qP;I~I(5{ra_>H|JEt>J08}+eCjh$=p~ZF;F0|QS91b zX+3=99(?5mc47$I>hJQWuMstEDt;>8@kY(qxI4Z`a<4^a#i2+X!F5&LH7P}+t-)gh zt>yOKZ!cG+5C1T-n;5?Q<}pPbwlU?TzFyJlpeQesk?xwQ@u8DRrutJtp^|pv&yPN7 zILU-F_&LxZMAW1A7=2Pr6If0aPmj2ch*rQWNY)Wjcow>|% zg@B!UHOb~zqhY^HjUQ;)qO)f#7sE4QUo`-d2=X>3g+!1@OoT9#l z$DDZ8@W3ech4E$VkIS<47L$G1tTA5d`cEI|s*14|tq?C@L`VwtMy1zY< z?WxOrzef(9au;(`AY2r8)ZZDbMrddOFn!3>16r zGlf7Vv&oz-9+Qo1Go#J^x@1)Z&KI0`b>>r~T2%VL!=8`9lkHJ)_WikBZkn7F^|7+Y zek7%%tYN@%kB+nxi9syS_&Vomze&#pRwk8qnq2Y zt*`*6KWJ$APQ81XciMwhVimtmb*j|>Ph-SGiysZ2?mn@6eXXsyHqAy@{SJ&_6O7^Z zYRecJ>hQBb0fR#$r3+`Q}>eVjcZM?iyu+oX1jG5_) zo}TI?vzh5BHw_(4&C21fDyh+*H}Gzmsdqy|!?N{4mnvUBRgyaXu(Kw`a;Uu&54xR~ z*C;&k#M?`n1%2{-N#V(gN+R{nPdp^WtZlR1HxBV_VR4_q6D69|O_!ZCvFrP^+HL&S zGU~!bj(Y3E(o8x@rp=FTv1;kRG;1w{CAo=aB=;?1bv`aK?X{B%-zFN*+ElPk@P2S$ z1KS4b`Q}-_uxzN-rbl9GutaUq#Bi5Q+f#p-;pM8gHa%%0xm=;u2}T;9B8_X3-8SyU z{>Mo)+}`qdWEk(HjyBvi$(S0(qGtv<>l+wY{`&ezH^YAaEe5!lOHt}E7h}6B4MN@0 zeSfHFXqy+yaj9)LfwyS?{-Pv8S@+eMgWSxKYqoQ9i?jxt#ikAa5GWMca{i`jgp%D* zsfI?377HseP<#dX87}2zig4o1LAVPu#^KP3qG3yDVJ`n}>_knUJ8N6rDLXgqtySS0$Q!npu}VHC)wE z9;f@}!Z+MwS*&(SayJ(j*XOl`%;M9&a+#E8LsNltp^rcZ-$E%7JeM|}*?jA#H`YV# zx-;YD-x@7ztqhjQZf{vX)LA)@#~T-Iypk$DBuaU<=-Bps`sL8Tr4XQ^qmynwHSWG~ z@8<%Fmd*Bx@6@H*RIAC+zRJqV1cenbClk|q(kCv)Ue0u;|F&K!R5I;Mm)%H>hKhiqO) zU)=J=H#wCHSK6{JQo{_hrlagD&ddShOs%zQfX- zYhC(efBvKWx601XNLUX(7b))@YYiR2{s{%W!GS7KmWqiDU#B#v5G1riN-9IhWY>kS zVR;79O#CgYhwkzvKhJxEZAE}OnP?o-92@zo_-X$5P{;bdm-DH%Lrk%om4{0}?$@SW5Gu(NHRSMz3DBjA{!{O@HzS}FdTuh4MDT~$=IQ`xc zZlUz-Nh(b4zHiU zc*Jg^-WZl}GB3V4vJ+Qc#OaDB&W_&Evd0#85z{dpm7-z046$}qe0V^)3Qlp_aAF?g zUbvEV&BE5~rP@M?@c%pls=cFi)*UhCHEJwRWIc>0&n~6HY2`$oDICdsr1I;FaGf}Q z<&59A2<6bcj4ye|dPn7Pa1S*V1T+=&`-o4!@D@-VxFKNG-}K9{bV@ENY5&u+UC z+sk534mQ!UUT5c%WpND8H~cB31T+-cE!oQQ*@G?AcgR6F?%}Sphs^y?n#1M6Gg?_$ z>FevKSadzuYiV2%^mL)j&dWdEAc&Q{ep@6pBc2)@6twru$J;j;CEvY!cO8e^S3~VB z&VC+0?b0=0U!BbcCc@g`0PbqxByMi2|EXj7Uh1)@H1)$3gM`wohk)6-ve_I?)u>~6 z#{qMQM=$HwZfJ;enCgg%OJ65pkH7{XlOkkuxS_8OaEFFbe7Qu@*=+aLhnBU55f96f z%{!u_qY?huoC`FBBFQ6W5mlF~2>izyg)uG=z6j4l= zdGV*FGRV|Em57xX9?r`Pvy|n^Y8$R-soG%K|R##s&buxCW z?8cJSd>arid_Uk~v04c{)~4nfIW(U3UNu#;l;H6K`42osb?H+FS8EH}PmWGcO=#R> zUbV_an2}Mv?D4_7{<`)_dFyZPdz8|XzPIJ%V0$U=`OCZF>X)QUk9Syh{a^_>b^Q3X z_a;FAZWFy(jJzhMpI9_L1_uTOtsltCgzcyc;d|Ah)Ni&rjJHLU>pIw(d&9x&)(qWs z)5Bf^ubeA5!fH48?bp}~i^&MsoP^oeFd4woT@JOkSgpBR3bWnTulC5Ir@|%VL|O_8 z4ri85U1M=!*fr=V9Ng8aA-Q#h`oNu#fMplq z%wTu6=0TI)mKjmqCq9+8hB~kgFiS6A0QY;xDt%k2`^-qj%+xQ#Zw7~W)(ab-x_;im z!U9p^v)g)&nmGDd?6FiS6GzK^>=}tpmy1vT^R(Lx)VWIC9gZ{?_)Pt*HqGDhO@KGU z=BAgIS9NvuI(uiOhhhhN#9PkZtG9$1}YN}0nrmJq6;E>2{N zZ9aK6gv+`m&;W492Cys?OHI4+nu>1%SEy}20-KI@0!5YHcQ(>g_KcK!m*=3V_ST>% zguDyfmEQHX>SC<;1CFJXpl7gVzV&?3ftmXPH_ka+83+1A&^87x>Hqm)j}iQ9d6HRc zZ_Y`RH1qeHTZ=*+;L3pa`aZA!R&yQMQ&hL`5@|8haxp6B1XsxQ20%g;k7XwlE~?c; zGGZ}2945~esKMnXINbGCmP|@kdnjh^Dk_P1(k`@(oBOt?PQ*%yCx_g`2EXn%3iygd z7Pyh$IyF4}9HL=@!;|FkhP(WW7&e^H*40KD%m!qV+%2G)ouLt@v&&%Bs#Qt>JL`>f zCw{#DgQh2Spb#7Rz0H!AQn|?y+%0z<$;=~_=kT7yu?Flq>n=4t`b!jHr#q!5cCt<@ zVG}(KwQLg%q4d6hk>_Hz;DE&YK>`L!&Hjk;QM(_Rca-fhVHz;Z*<}(zD?R?_qn2t{ zs%!A7N1sTy9YbEVYKPYG^7Ydi5e_rcr7@A+`wC3#2ZBm3d-?m9_j)Z_w(i<{(JhDW z1t6av90g!e-rM!YrF?l8V#9jz>%Ca%GqIQ9V52qT&*%i4h(EV+#`P{k)XK#-S0TB& zQlH6A_c}hFD&v$GKN~;id5$ig0v!hmquC9fWzrgae}2?Ma{Pkt2jpWqfcPx=@K1JF zsz`-*NHro90*_gRADPe|RMo%1E;3SW)^U72wM_|R&beq-3G44BS^VjvXQn5|aB#E{ z`l?Jaw2vQWuFJ&j8vR-n)qP;Uxuqqe6FZ;5JTZ&T>+eMoZSRPt>*xEqtX_V(>qi0& zEdx;O*@6VP_$H@?OWZ{oZYOsO-N~M*IP0-d*tiea{?SCf(XfSrte z*WgM5cfT9u$f?8>>&>i?`OTypZX#o<%Y@`!>7D0_!{r*YJ)~L!&N93^NrEA}XuMKW z!D7?%JhoSmSQI19@+h-~NZ446e#?qCC}Oy|&$oOPGqXalNFOrDIPr8~4V+spjnp~hZHstwxoS9b34gGdc{?y( zqrD&3^XcRJzVr_aH0sFdCnSxW1Ta7A9;cHgR@=M)4nRs-YHF~3ydBszje z-|EJh4&9mLWG?pOIKC1DDTc}1P!etSS-51?$EQmv_INpL|3$fXVpvuy*v6A=WbTi9 zq*i$Ub-KAsfew`l6>c&rsjI6GmfXVoB=o=9FeO-l6e0pa7fThhn^p3R8Y{!eO)cIA z^9~z$`9Sc8W|Dv-+KlR*uq5#wOp@OzLkb_ z8}4a(a-isJY!~pbo}Qk$!9xb_r0f|G20ONIkFMK7f&LN#q7R+wG`RX=0hM#w>q_<3 z?eyiHgC)v9k2Zs&ym6_=g1-Fapx=25tmhjtJMU_H<&03;&8kKa5GVCZ0jqDEDB8nN zG+7cl&7776dZql4I#v8f8Sz31ljn7`w8VhnVD@M0((U#X!MASVPyY3Iz45z8n@s#) zgLr9|;7TvIj0KTW@lp@b$}QP>5i!b!PO#v1GM}}fN*w;>NABTuxkxE zER?4lN4=QbLylAHVeqc3lZX^|@7`T`N>1W1P-{2h%6);AoS$ZI`;hwvfSi;1S6T`J zWtM|(&@YWtQ&U$HEL9V%P*eY=m576rSZf`X7YW-5x%zqSp4E)fKfZxeT^r@+pN9Zx zj%XWcd>QdRHmT*oYF*o_Kv9T$*7Y}6KjtsaTZ;RZ_!*3}BYQ(!$FeS?mq=~bjkm+i zjDVr|)Ys*;m6x|eR%-kk$SPRF9=%Lw;@AxbK;3QazrUQ4{y~z7gI-s5sIA^uXVv;} zfyOZ@LRajWJuC-~GzSN6W(&sa{uWvFwWn}L?WIwO2D3AgJ;G~zDR#mY0L%}bss?op2FPUked zDtteWLytNvJDWGj-FFfdFqahB6-vJsY@%@O^(T3!-&-9h!O<^fNEOza ziSuCvE@F1{)Tzn$G>p~ABEMCsk36TNGM)9F2Xbuii{=nV%1pA^wZTZ~ zzRWbQO_N)`ASo$HcBA?6G+sGzxiw)&0sz*sp&osgxaXxA|R@ww-+s4kRv^9O1ry^yn#4+ndDI5U^DFP zUg!V8%ra1Nbb#d#fCZyBrd}_?!<)eWMe1cr&5RWtyu0e4N9eU;xn&uF)vC%8Lx3`$ zPn9Cf)Xhwfto^EyXwmfpUnUMpjf8A)tmMvssBv=+;1*{g#K=JCOlrU8y_fB*Un%K@ zdTSUCGcm}kGZ9{S@U)&@y(tBnE*i-NT-7V9ftF>q zcaQg0n-#gSN{nY?34v$@EW4|bNi38BT}9&R17WfNPiRX9WtLKRHZ9$sCYSkb9Y4Ls zh?FvN%P}O%PqGs)RhoX;Q~uphW7|HaMTghO1^8?=I7A{uquKehwpYM5rqUGNNms0y zXYwQGl#W+H=oxy}3+4ELN&XugfVa@?fn^1^WgyO`&syKudE0E z$+wr@JI({yF4O#BtGwS}M>&5>%mJ-LHG zeJIAH^pwK2Ej9D1S`K-zKAfim05!nZe@mcz}E#co7jI zHuPQf?C?gNtDjih)VQS|taN$g$5GG26zY6~{oO=&%2{`;Ks0=zLGabj`s5Vo`#&_z z+eWwogx0f4Bwrz?Vv%|=S&fgL&T=bE&x5CI%*ATtFe-f!-cnFl$SP@j-S-1cZ9+GU zEyfCjpnQj*r{P=U8Zi(a=QQRSHUgZ?;4E#BAY?9%V**eWnCwf2hHcoT$dvipFQMv| zy}?YO(I@Z%VO*>}ont}8Y_Kb_uJ--N12OZ9+{_n$DiwbMI2al4hyy*HmUmrFF1HB= zrW{fYG>Jg=7U!PLepGMp%}bD$zfkHiZaOos*bAG_MW5}H!?9t(2hc|SSf9mCJu*D9~MFcPhjC_uJ2^@m$ z)gK>hqn&)cZh=g}`uK?ZAOcDa?Lj-khO0s9ElTf^=TH0oULyQHoKzZ!Tlmyc#R_0y zBx#G6EHM(vFL|LFQ8(El?54^t{R102%u)aH+}8TVMG_;YV888&U8vs>v)RO_FDL{@ zjHe1`APz~v@m_A=oyniqkOE-uOJm%uW>VSdK-Xy1aMn!^CxR=mf1gPmR@LfoytfiO+G1%Xr9xB`E3j}sOiV$>RA0^B_^lBE@{X>yTrW^t*>tth-Bmlz z7X)1WVHmNeR}l$g8t~7eZjQ928^s0@l%9{!zs4ctM{EX2bu&?;0%g{?h~0P2^73-d zVnD?eKs)$)!yWUQm&&1%OIB?^{Nx-rJ&hnMyw7F>7d*oKwW0za*n)SR{d_zzk~PMB zPv@CQ9S%Sr*6u5A0`@@CgkK?Q~_A80@aG0pXM|)?inmcS~Ysa)424 zvKK7x?feP>sK}0Kg&39I^ok;_COx?Io*K z2YXysc;L7JN)jkUbQ98djoQ5My?SffWe_ZEnPNL{Z*f17e60=Cg+~4u>K2KxG?Lek z*5Arb^!m(+$66#8z|8T_W8G&9hMyB|xsFvt0=~wqw6{JhK0f|QHvT2+w$41vYu@3r z*q^0G)zvMq1)V5z(mBEPB)jYtrTK>6!FRLi0K=Bn2&W^2p!jhLpqa2e3uF+c0GEF8 znC*M?)da=sLvacB6(0y1!q-ksjCdpg5Wls=gCUP4k?37&4SCXcom^elx~s>?DL`fn zwr=s_#TRE!UJr_u`TqQ*{@OSsPUpCE%WNdkqgTpQ6Evxv6wVm?(%%g8p6p*2a=sLn zV+2S|m|iSde}DhFC$p>Y(Vs&h&r2Yyso-PIy;*llvNVg%H4>752+`3M8=1cZ342&+ zd;&mw_I=j^_U;G1{{CDkiyTGXAq@(Jy|B}vNIp#rpUw$X8({-jjNpE2NvO$#Z$zE#^e6(gH_l*FInH0&bCzTz>Rh#9j)#G$f z?lM7mzK#IvMZ)$E%Kw-ssrIN_>Xrr=N zTll$hI?r}zy$d38B9pt@3I;}Ssa0u?fjVvsDM3p{ozQal# zUY)s9MIM>N+Cxk(%$p46Q)?BJg*}HYx~r}7cO)L)K7-m6}oLCi<#0t=g_H);` zAONr+aN8G<`^zk6M#0E@fEcgTbfXjj<~)tO8nLHPa0VU)ZQy_vUCJgY-qM0% zq1DenXnFj4kh20z!UB2BcszW3e4KU98yJ+^S}dZdi=1;S`+Kg2iFG{7Z{1rBcnHS2~a;A#U+&_$`@tb z#{g)%vRCRZrP>a!Syb%9u~eEyE>l_ok`T`i_-=Dm)D{}ugbRa&biy0OUVjE^y-0H6QilE1DHXbp`JfY|kcI&vNNX{<{P6~k zP{qV+rwR?^3OT306#J)?%}Y=_abnavbTT|| z<3tAs^J0mg1)|;^Bo0yOEq*A z!HwC}P3Y`b3R@u43q^ujZ5F{ASxwO6IP-Jk@Z4Ojb0Oy>4JAEAhbItJgb$z`%f!Ms zQ9e)>@B8@m(SBVNLh~9y{jS@2T3UNV2w85gDRH%y-KPd7m++_VOicxq0fMS@XyZI8 ztai?*!klO`=kG~E1k|-|dz!Ze{y1Z(>?BUw#0|E|Fr^e8Ufz@MuIP<6V{;kb$>9K1 zy)lf3CAs=LjOFVBFsVlfU%rhTEy_7fBh%-mDWmwd8lSErUPgVc6jzb_%0T}#;Szhz zgIo2hEd4q0S5@~*8h7A5_hcPoVPJL-AheF7m;qB%(Bl+Yi zHclDcOZ>K`EN+j^cMq%00VT-@ma^#^TwO!?4DfRm|euI|RGwHXc*279wA* zWIS=sf4Rfc5991P?;V`vWEF+zX=T1nKUZnNFZM6J6tg7nbzc5%zCC`zOS%*-mbtNW z@i$J7E|jsdvQY6c)O)8UV;nW9uq5lz!^`=myLT*YWo~Y5I{!m-vfsH4TtKLfv|FIA)a+3TD*AD0XeiD>B6Sk!>-UM8|Iy4ZG{C>n_TmFmIBvBI zSD*6qEe*R>o(jXSZ$_!Hx>`%3^EDEc`exWRveECRS+j*D6e!gqHsjlN>#d<2pUE-7 z;nU&E_K=FK2PpP;ofmOh>YD)!9{HUF(|g=L$#$$CNzop@A;iotza6%HH~56+@775r zr~M_gf4%*iZf!W09|JO*H-m_iQri(e9}rjd=RucVh$nWM1b4ldA4DwrlzLr5+DmaymM{^irve^*y(w6BTvNhly51U zPCkV96*J+Q6bWBWg>NMH4a{`YKOYovqW2Hoa%lKo?4JV(Z(6G8JlRp1rRUd3`t@=V z(v;ZA*PuE2qUixA_WkzBX8M>e3{h`Fu&5wZyM7xe++Ej5keF*u2`0duw&-rDuAE6@hddEE*uZbJ1%xV8G`KZQ~}h5GX- z>Bg1+`lYDn*EdVp3?ue`KOY6V;ZLnbiXXk*!3f&M-w?-PqjWclz8LDbJpLs}?&3=` zribEsp(_mC(E4xp0(X+7cWdLSS;L_*W?Fo6dU_<|mhLcI-BBJI%JI-|2xE4e3*n8k z@=c)Ro7W0CVmwgqq$#)$w_~edI*4q4n`*H?HtH%O1B;Fh04~Zx7ei){M!Wb=AW(v- z_#BDtS{FEn1>59}HhPLggjin$tCPh~)S7@wuak2rS5#6F~KgWyc`kG=W0frsGb?>uM372W2pXm|uyR6)&y5D%lx zxP~AYCvpd1mS=@FEUQ)#v5z;W;^HEdIsCo0Ot)p}o+mF)ZcB?~CW&zt8C4=KE|f75zPS)0<0`{h)zufDV%_y4IQPLl#Y#2vpT=J5JUexdI*V*XO^j z&t6u#CErU>Q?7k4ac=|UyE*qsUJJ!)t;5U|@-&p_-ArtsKv{3d3l$@A8h z*^VBPz8~Aw>!^ImjuY66_v@!-R35(XWvUp|}&2UkD z_{+VC#zcXfwz7UUK(A?6hHf8RId1#F^keGS{`m8N8_%YqHVX0`QQC?`n8gqh15d~i^IWuJ_`@$_;3@b@$7G?aUKwKTNjk)X>HGf_B~Sh zln==(Xx9vFP%(rNf;PwkO$KLdA&mt~qoDII;hw_lkDWxSpP$7j zm5P)Q`99QgQMF*b*y#v2|Lx#vDDMRfo99!~lXGN&=gw|I#{o4=BaetYF+BpEWiP$Y zOJoR=(29wQ&XgaT1V^{^jU+`~5S#-X9QObo-Jo(#Ir9jR3I#-~0#T62qmJ`vc>JJR z2=z<`+xG3ZykRFG${bK7v|zIW$*s=Ws- zLtr{UcmGH>H(FNW>gvj~3N$j)(SJd+Ooy;-!wjNLp1liie*@bM@jb7fGGsglSkM)A zhkXvv&A-mR-Lt+ML@`nUwe21o_X;I_LOQm2?q-*-#=@M7Mj2+hT_fXxl!N`Je3*~h znq8Y;uwXZx=GJhztX9Xg)Nr@Z;mH2N^}T`JZF=X~CD{*UsVq3S_TZVdG>>1=@m^fI zURGb#h zaC*U_Z)?wY28Og|i7T71F<&@8ezV(5JG6V^N@1JSrcD{5lgaD3d0!lFcUf{!dnu;N1k+i!4u>(}E>7e!&O97rU-y@s!C+QpVylkHS)AaOdRO81< z_mu9;kS-K^qqQUH8mg)*7kq%;h@=TnzOAy6rl_3NB>%WDO9hJ`jiJfQgS|k%Qwj0> zfz7)>-x|fG4Wn}YJSwVVyp)AkoydwlJzP#v{RfctvHzLNDf8Twpj0wh?MLgO71%-N z)XfmM7IC7cxW z_UsiCb7Fsotc^%*2dd^%%dZkU|2E5#e|w}wd#wbWP7XMeLc+q}Ey~I)p8NOwNH7BL zNVMsa*=UUc&pD2}{a#x71=KXeb`osnIKlX(TB1K(AanL?BdG&HwKh2U77Oh+w~FDT z`%Z-I^~!BRAsf<;lC*r`%RjE6Fl;|Q@5r;ypFbOxJR>K~hX0iZgnCkGwKxh&-SpT6 zVzGbj*#S4@7@Fj>D29WWZmUgz-0I z$v0R!jx;3>rkeGvl;dsU9dFMv|F=1H>6UD8LTi-*=-(H0&wd%K7U##@#$0jPp&N59k2mXRIdqQqa1Tc{`uDv?JgDj`IGct zVQWy#J#T-;K^O(yjNH^L1P#V3{b;WwQiG?7*q50>(YBQjf`WS9HM&pLTm%OL)kK>X z;b30{WQ^ks&6y5m?@2L)V}ro#6H@~r?Pq0+@+4|RqTf>>kuYl;LeUc%Ow*qOM6nrB z4mV*tQe!zHTF*9J&xWl}E(2Aag8Kil7sp%ZQ3NDb5^6u*7TX0XArAHI{TtWg(9S(j zn7sZgT8hC4YhKy9Wh5#piqwZd$e)6K`qs6m<2*y?@d{6bL`mY;R`e-q;dnfF@F4ju zgsiB0g#=jbyL`C|^t>cGM#{>{*y+&KIEnNctv=o1-bzO}z5nX!!0pt=dK2}_^}YtM z1RIMh^q)Q6Y~JgAL~(mt<8#8GSbBwOla{R7%cu&0_cL5e%|?0io!LI-?Wm=d zSsrzGm5sCb^JqG;{BFjm$laQ%#!_!NZbZtN%qYK#$Vp={SNz1fX{2FKZseIFN}EIP zh>XF(g9fh*auNz%+0`vFUa;TZnzut2;+wmPt8JXCJDw{%liZfr|Gt-zXOXRnv6TDN zi>}-Lw4Hpo{Z&TBFSZRYr6yF3XL{(G$|QS+*YCObowV#gzhYe7jw1R6lbO0uK_!Wy z@0ViN;Z4|QRe}yQfDy{0c^c{c&~2Tv=}OS9+xE4hV)U-=WJ5_rX?Zz{?%Y#4hw_eW zS6kszCDnRuk<7tLiT!Q&%q{%1>0=&mQP*@2oo9ROqVg-*oV4@=2g5k-#EaTzFYC)* zPHG-=bvF4`zjO1_{a1bJr0y@0Sv#F#64)M17CUL+_S@qD~De>ao_oo zNS>gepTQb;-}-H!`1ZX`(r^&3wnH{vLH{b;!2g zK*fU^oO#c@Bg+?{WR+hXV)4VU@{~Em22ab09TPAhWN?iWb{mU# z2nj>OOi{_QbHSk@uS=EBF&M(s#%D?ot#h~ff3zOnv>OpdDGs&Vhle>{CtZYZZW}^y z^3H8C(5XOu3`IFT=tt2?@dI!O(oMx{n;O8O%+3WMqHxk@S6&`^C6iuaEsg{CoX28U ziJ@WWq?H=|uq%6mxA#dq6y)26fXRZJ)j;(XR>#5_gVUMVE4u0~6k5qPjq|AYk_13L zo+XWcH7wElgC@6>#s?K%Ev@1Xkca23)KI3`ort3G(byJHzch0f47+=AqY67%QKdMX zz4+9~^mIwHAO%h{6w@eX=efutry4`Em}9_urm^q3pGIj}*?Hrts;YngvG6BPwxktK zuRVg|QA=6X%6Sy`!MQlNRWX5q+v&u1kiIDU42Xe;We%9TP%dJ*EBiocX3pc=UaR=}&rkfe21g!*9&((&R_xjZ>$IXE2FBA^^S~-> z^Ql6R*$WAo{U?l3DL&reI3JNf*}CS(j|a_aZf;9@xE>qMCjS&9irHfS$;%T>&uD)YT2Gg<{4A7=2=|6(m&X8aix3iqew|sO^t&%xz3EHpXeZ zp3xydM_5P*8U=GyKSCwDNQo3<0*b3vY&m4xWu(jd(H6@p{F`wdfpI-Y4&=UGm@lp7 zrQZJi`?JujN7}R5>5jeljy8PKK1yT*ZD!U`^rGIexd0;#M7Pde{xJcC%SDOe23Q04 zTa|G6;zeM0bWRLN4KzW#0ihAJde0qckOrmC;a#82p<1-eS_nzdsa)?)YGoZPKrG*( zlUiXgO+%>&{bp46VN?r%G0@GIO4_13I-~}i@zs!d!#S%OJTzDT!vMT|>Rhen-f&xU##md(D1N(F5kX{ck*hCk9UgWydzukX6}>a)1N?z!DK;EIc$BQ4s{y7VQJ&j~^RA~lI+rqi3faVX%*Pmo|P&R4m`RWBS;&e>8XBh9(rrnj0X0<(5qSLv;)TB7b#PaQ@ySRqkNI*vS zv{W|kX+y_=RE(V~<{T_H)E^NzAHDgxtM5rB2B)yJ$j;gnhzZw9!h;@>J-F7TKF`%P z@AZyO{T4bbH0+{e9_6VDEHQ5` zx-vV0nK)8KR1Z)e4O)O4S$*ot2jrG3?|dX}xsXXWlg64~zv!HPqHx*)DI(N%uEN`a zd7$V1n{Cp=nco`y^l5WTkV&^z`SJc8y&iL&qCzb8!RD%*RFQClxV^h;W*ty9n!JAc`~ z;G3P`&**Ky%OtM?IGjw$m`Is=Xz_Wi!M!!R9ve*+%k#6G262=lJh!$^Slc5NXzQy? z5CM#1v~9tiCa;3_vGdojoSuJymkU)jJLgf!+X5BHAPhE#$*ZcWs^cqir;n>o_vl)_ zIL^cB(P|}zOKh3D#6FN1N7Odo4?ts%AJ+_YEL~d1)zr2!-AP>r;`1gUKL0|_hXE6~ z9iY{-9&3iYjg--^-Xc(94y*hU1@N2Pfnh^wX=z4wXisOhgdCHbWA1_>NGkHGi4K@# zZ80ELMhR{F=^>VDF_C)yMbdEr4w?>Q6QIM{kw&O)raK2#wBmYBEB@g|a#(P^OP4EO zH$$&zbn#;Q4`_GmHcD7uyDAF-MDA=s<&ZlJ689VS9#_#G~ zfYxKXAe2rP#u%^-j&oLf4OUw-2n{RIpnm!CMP~V%goKUl&DaLgLqhtsi~4Zv2!{w3 z$=uvrz@TsgxZF)O7NGw#MnCZvSAhZ4Hv6*YmKgc+-}0_-^5ypxI5De%G@o~Cr{j}7 zwyk@yKqBu&ks3iK5hQZ80YetzlMGIttb()yde+K>nCR$KtO(ABiG>No0u^j?JN*1; z%nerDttCnaM?1F-tEs9EpgtfL$G;2fkB;?gJwG7eKnsgLLs$YsE9#b!z1`ij6D%D~u=vIA&rR$)PPTPHoG)cp2pu{G739VaLGeNa z%(e~-amA^$PE0Z zzQAYU^gg%i;vrKNFz^$aoU@8p_9gmjUABTIapa3M_moH0^~t(TT!J|OCq$AhNmpCj z4w_^#ngAoP#KKQLJ@%r&fbZY3j%Qi~e0+H|FhgTc+Ek*$MB;WX&G_?*D`vu}Bhp~* zxN_2T`c8lwla4!l?QO&lf*7CjsD5Jo!IViy!2#^syHz?rZ-Yr{^#& zg{(0?3EOj&=%vf85kdnq22<*$V%ML=Ufekmy?GTcSVBM7h`hMlExId3qj5J!k00%v zXBxA4)q>{HqT{ApntVByE|+0pdH7I&waX3H;{dKFB9~@e8Jn_XUxnCJ z6wQ(zV9&Ui?tqS*CMRINV0ZhS|B4YFK=A~GH-hv?ogwpEMkg=+Q6A7IqPMpiKD}M zCx}S%SMiggzC}5;Enk94IzTaPkBf|pE2`%7&i(gfQ|@Hb3~aMm4Hg(!1(xOCe+>19 z!*THrF2{}@y`MBRd;l55{kb?G7qxwuQfZMyiwJ0q9h$G3@CX;7R5UWGO0hKA9jNg1 zm^RnPGp`~eBXLg$HuohJ*>$6%Fcrfk-9Ze26kKK5pX&PQb=Q!(`nb_zDvU!ebAlXM z9|^rtf!Q{U+`cUOrcV&>f&FyHfW${HUbKaA;p}!-#zz>Xp@vBmZ!*h)dV*1^K+7Q; zgyA}}J6sDAc~~nXdt?iy?cBhS-Q=)UA+9-z$%RpYEE=KknogjM1^6zGk_cq(IJE;o!pYCC1G&jZRSQsMt&N5 zgoMr;ue@`B4{b7G*&Owk2qx5*TO-N+=QjZ%nKsJnju#uSf?M4?{_6=Lpg{%=L2!wH z#Ok2|78A`755J=4ptd$dffFO9s)~zWBm9!KIfrq~O901bJ8n#RR~)8W9azo3AGMl% z*rW&_gg*Q8R;7~4C*386D>k? zGZ$xN6S;EQJsVF<%Td3btX;G0nnD+%#`O_IwM}9O(!@&`wl$CgstoWz6mI_#YuJYO z0O_G{=w$lWz)u|b)qJNaysK8p+u7M!*Hl*rD{#VfVW?x9xwf`8(5I_xk)Qel}Jsgwpd^|1DUSC&?AN33jqJE136Fxm|mZ8S__3|iv=b(RGvzJ zqF`_w+Dub-vQ4<@fBnW%v0wN7{QNM2Ga?Z)7d)V1^&QYa=i*sJrU$ie-@Y}*8EcL6 zGWq`^Rc|RVNpcWF7hT~J&qW%e2+3>R1c+y3t*xqR2B_THc9XQqlTNfmRJsm@=y1ou zW`lI3&{s$+F@sX4p^1sI33~Ls{a;a($ThFtM+8&nx$Lw(%z*$#7+4t(!99DPmP5>1 z*kPoNN?DtEOA3p84#H8-$ol=8 z$sInGCwGYQwR{}%i>k_P@AsFb4zGf^`h+mq<&VvtMMtT|IK9aq54P&pjv0 z)V+v&r8vIw$m~}he=)sg(XBgZy$>(PNX4H1aV<<>I!sp8ySzLv^#2Fp`+iz3vxXoy zA6mUpm06E>CzANpACi~^5K)fa`MeIJDQGs1!~r+4wGJipxBHzmu7vYM#XjY&qM{;v zUMNUmgOF32nn|cVVVo|@6=_5)C8lmd--j)w%c{Y|MyzvA1c^2qM0Tci|N4dgP?-=N zGSUc6TggOzJdhMI|1D#@!G)u~n*`B3$DP&7W)HRN@s2RrGp{O0S*RA{EqYM~%^8al z#}GoC@1l|_w00$oj5gXXQ6JUO)tw3CIV;s zoT7*b%Wa-{WdtKpDCyk3n=ikrlB5@7dy1>rrhds{=Gi2w!hwf5OigD_8YeaTXH4q1 z+($Yp%4OO8=)Xotd-C5`ob8TUZgqcI7I~)S+ak@ss-ly-;NzW@Xjco7>__vrLz}!V z!l~}d(;qNg!v@a5tqN$Cw5*^f>41%mO-o;W)~%T@b(4z9wkS|mrK1jETJM;||F4@Y0d!<~AWa$FC; z&k9NI+xKY8&c1dBqVZfa^?DyP3D8!M?xZ5yM+i?ab@93kjM&a0PD=V5aaokUasq?t zS>XQiQ)XY#BRtxsdvW%)*y%KInS|9Hp#C_mgh(z+z@%VrAd*=<0Tq@?@&_LejL1N- zWUsPtsGy-2$pAKlkmJ8u9$Q6mpCu?3IVq{XlLMiyq0!nNQaAcpm^7X8&-C*rC2t~G z!`vM;qcmb_f?*g$5Z-0tuoplP(a5}!L3YE{j~FyIT?L2@iN#c9k$o95m0(7c=0u4j zTdqkeEYf2y$oY7S;=yaQq_fjOXF{?FQj{5QkK#v|!N@W%1y1lpJie}<&b&IY{r+LV z15(RJGKvBw`d}SUi*%KM2J9LU2>rEM%_uJP^Wnp#v|>W{0L*3t6D?k0+#BX#BakEGU$g7Xfy{HYzG^3w4Tov-a!!@?616EV7uj0$htf`LXxevB%~(P&Qn z;IvTUA*c;g#kuF&jm1d13dB#mq1Z1&=&F8@Ie>p^vQYiAvrkA)w(kA;VQ6Rwr(qI9 zC;Hm*5#iqf$^$89I1Yc_JN-EY`U)Ri(u+P#K*mVc`1y_E@*+2+#tWD(J%AGTPdUvliK| zD@LBSUH3T>ItPB#U~m5BYz~n0jKe#(;IdH9p6J+Er8uHI;P=Br=jcZa%sJTom!M@h@g9cuQ?laq^*e5R_FrCfIbiEHKg&v_3! zH8zIB{MYM$bwSI=7s@-xP4_3CgDAyn9US=LLvu@M3X(BN5F+3GVh)fhkin$(=fQ1e z-v9mn_MQ0q_9!|~!GHhYxV#iB?q4suSqI&Zn6!A(^;N(zu->13TT>-sBBLn9-31V&k4KC>^IhH1o9<^-Q+p5_vwuat+rL!-LdiqYZm)zb}(dp?|<{K;_&%{E}P!_ z(JtK2+kKp-17T&skz>mfs-|y=MqH{>b<{Jeb-X#gR}$#K%DB; zRH2x)U(WghWyI-4D(9S)SN<^%qUZfrDe1?j^d&^7cMDr_Ug>_PD5JPR;(mW7_5lki z$+?!2l2W_A`h-c7CXF3CcHB7U{H!F>VMI7RPzEH&3HHB803&)_k56~hZHGiaEjg@L zg*vJj^A>>Noyho;XOC~-o=4v^y2!3H2^`2dCX{B>@26i_>QJ1pi+EX(cM#lbr+Rnd z*|$|a`>@N+;|D(g$C#P(gkL3uBE);%z7!rG?G=nncyeyxljWekM_S7UZr3pM_8mmL z^oMwu8I>{dPsB@8iA^RKethmwgaHvcwHb5=TvZwZB6 z;e%^s1mo85nwg@ge#wuEeo`4co%Q6|*1RO=n3Rfm!}I6Q&ujIvzayai{+QJSH4V~* zV4XILfM&E9;uFd!w@G~hnjkW8r_~e)zveGOZpbOR94sg`B?h?As~U+48;F@#0Wgnd`Bj|? z<>q{Z4LK+}fM`urd6ryCprAa3>dW~q$Y~uve)FLd%h&sd<*UZ>U4S582&(b(IVm|* zknNOQ-`Vcj`LCiLM40s)VZu)%iq#d0TXTQ-(Y3Q;Yp-d$XXFuB{AdgbaOP1Q?2&5} z4UaC8ngx}J^HS>iyw#ChKI-Jj9y>?W5FY=7r%&Aq-rfa}`g{&@G+F8WyssdxGO2`9 zF?aPd#CB$7PUkK@YcIy7@-8mMmc1|L#*~OwTThQv7`KRoWT8A@_AT!#iN^9br%BRK ziu|MG9KFcDD`U6LU<-=CevtkD<(C!2V^uK$u~yYTjCk#;&~wh9hWW+#CsM+-|5)Jp z9Q-2Oqaxq5>KMFy@uQQ&SN;8YZGYs1we`!ZrZ|}`y&7cFLVgvnpeO)0?5srXl{u^S zZ`5yJo+cMa+|s;n-@erIUc0i_*E{kaYgQ^#^Bd!i|Mn$O+I~dGN2ZRHW)4>Y*7H1j`|IDd6$*d!r=?rN^iad~Jx8$;ij&xe z9+4+bd~RLA|K50%2a`F5hl~2340obd_{~pOO0gW3Ma@@)MMhSTAGtutqgA5sfG^E_ zF< zU67kt%1{U}k+P6pNm|*{aM4}So+R;bVKv-TX*3z&?etuhp{Q-WH2 zsoRM7A~8%^JG=Q0m&>waA=acN_F0LXoC-KgeOqANI_m)H#Ru;?xO>kYRn7jVNeUq= z3q>Pm5B+cVT=Oi>p;RPIa}`AgD;}MvbnJ#E5A8BL(hn&VF9)mL5;}?YcF#;@HoivV z6V9rdB%Bn_V^-w@IGv~RfM(X=*5*2D^0@579``7=iJ0^$GYE>WaEey4T9-7UuLCGh z%DYwAJKgpaw0|OFNEy0u@Jg<()$v~;`we@uNUo&95KGH{EFh^(+0rlyN2r@5HN&7`QyirF%e$1Wp0M~_USrsbz=bpg$g65 z^y)?(w*Z0Woh@4oY%K=d)Cw5ak@|EI%7CLydQcTcUnIEoDpZ?foJgcwb|>-~O8n=} zwKj1#D$NlovMwLvn9w03B}k$&cLK|zjMuB@{pMbAqK=1xJIz^1GwZvb4sd7pitCJD>64;zRecJO0O2x}pmopg2!Z5!Gv?e^# zdxXwC%MyeEZ_WONTa7O4zkJRNiJh~_RxL8%+$>#P-qF&WzM}i^9zI?iQ|yo<$tM>M z_gfhnY*dk7C3J!hoR;*$qD@2FmN|?-TC?x;5Ipx>m$Cc?dGqSJ6wI!HQo}0F=!3}V9?72}k@LlYtlzT7 z;j?F#gYx*bY_3?dVDm(^1`?IryQ*JFRtRtMDZ|lJ3v+Yxz<_U9MB|u?hYn>z(J0PQ z9dnNw7q)DLRAODp^scK_8SLoKq03AyfyY)5lN!fJ>;sYB#O=FN@fF6{%>al6Sfi}O zgff-walB#DEa%L4uS@ETChVbXvZ!)7RFVL`k6cqWcj^01zT~nj`Vz4?ZmrC3E`i{5fs>-8w~GkXwe$k5y`v=J{;X{O zTMUTi|4f?9m=7R@438Nesu#gcb5?{LKmOrm7Hx1IsQA@8hd$pp$7#eSE=mfNa+TH) znJ8{a#~mT#B&Z*J_xOuFU)<)tu>GcOZk9Xk4i?NWjxPBA=yoT)>?BKqgi;7RYhK>D8(F8_ zgKOD|yFF2U4M}Z2-8mM-M6ZFak7i=?^-o4l+0wP1ZpxGj{fwRp|4ni>2^Y~;HJp0V zqoUlPzE2}2#2IF{`$mc(V<2Fum%-^U44$%V>C(d{OZ>n*CUq5UX3%}~(jT8+P2_Df zq}p{a@3!F&4!nP|{A8W5We1x~Kl5?P?BhjRYqj)al>~|=UOt<)@d8BwZNF?QEYZnm zH*4C+GZDWwoi3$)o-kz`PiWG4VaM2Rw$Gmi~_zQE`aup-XEl%i;Q*$(sV7&RQz2 z$`bETxh4Mn7Pqw(O4NSqLb@zhyYl+xH|26QT~TQ|ZiPoznxi{rVuhJ9RzYgAE3aeS z1d{lKWlU`1H-pz@oIP`9HmZLLmyv1%b;Z+mZ?rx}AVfvo$rc%(Rj*De)~$f|v1r;v zXb!Op9UX5oOA&ClYK>}AK)t?mf_~G8n2LHwXnOBS1Xkmr>-oz35(< zFbpmiiG-B4^5z8w7R|45eJ9;5Uhl>~U;5As!scN|j&x6Hk;+oGNL8ANL9iK}4|DW0 zFAym17lOc1M}(<*{Ks$9DQlTbnh>)pt)t}>T2gmk*~wn^73=&`*~(RyrlTG0z*I9Y z96EG}kL%sr3oo>oF0@%I4KT(^0jA36S-A)!a4FFs&-3Qbe|Kx5@*t6~bnIL>(4(Xr4cdKmr3cY;O|w53bkb2u!cckN+Z z0b^!<6+k>$lGoY)*sue1uy{!wI-syr^!P7bGJa~|->-WoNe91HS?MOU25GOLvraui zSYJq`u2@}84^r!hYuj6koci4Pl5HwA(JJ3F9;h_y7Q#=*!om_hP>kM&#?et9kx!dG zz4fBGH~0BI%VZh%|FDddqj!v_QmdK{hY@4b4{5iY;Nu9IOWF@n?#WI^z^w2-~)3pPR^xDpf$JMJk zT|fEW_@*#@FO^wo8a^>|_?xc;ZZyf_7az?Hfb_s5%tu6}8Vxv@kA66~8{K&9`K0Vp zRueVP2R-?SXqUoy2a8k{rgxQWR^|QI8;bSOnlNDk;R2e{BV|QL*Bf1)tj&T-{p`0V`BN=e~^Ux-*WB$F3JCXNtC#RqF>yu z4v)VW4}QI@H0`e2_8qHK&lsD$TX|4|p2FFC;xx^#e`}&JHTT`+Z$PBrcJlPhhlku$ zcdn{%%Tj)S+-UPhjcxbkbyTd`tKad1o=X3G<6%eZ+vk3Nx{2?zJ7_Yt`1a_;X6sA( zoU71_?wc`Q;nqj6m!OdruU*#}=#9#DU3a6Et#kCn!*7o*E7DRZHlCR=I$n3nr;Ogw zSx<DGGE+Y_Z<5!}vjCyz@Htc>i};7^{N-lf=_@n6inPS*2smFE zN@=?ijm_9G@(Atse!?gvFev*@o{D6HOw5&}Z|8`t`HZr^uIzrUbovV@SeRa!zi{CJ zs@bq{Evb=*0%cX!@E0sxc(ZH)&U<2&^EP}jp8nx)U$INY~X4bmT)VfyF2w z8Z;@mmJk^5#4@Fl{uTX5AF7e-tD=)EX*R0+gfS%|0duR64uJBqjEb^Wuxu#tmuGkN zC=$jDfbnVR##Vd?u(MzET!sYWlVYtKV9=Z`0~^@ZFi0+m6?iu7BShorvR_jWgwWV`o*$h-1ZhHBR#P}h=*$x}OYVTB7ib?P&~mi_zojU1)a=qKfPw)FmJ29aO~%yIfiD1K{cJRKkj!fuXnPb_Hs%kM4jKBGqWLaFuj7YOo-iy1(F z*+=^7T^ZVT$rYlRJM4gh->hvHC(_aeYDjF%)Lu0$8D2>o3$e*E&T`VEMRrmJl)_Iyviy*Xs9GwH24Ve1;?%l-zAi*Z4ehBq~F?K;a z=GH=e-(u!{oT0RJE3XOtf>0W>!HJj92om6ajxmClieB#@)UXsd@?-HbGL*FQEg&=8 z?m#A8bx>C3dqrI~l%xU#(zR6dgk+KuwArV($Pq}`Go{-{ZPQgWncS&YGY?YFC>kEo z1bs;D@wM9ue@0cy=>8TlDGv7bZ*65v^|b_&*ttJ!1F6r_bovY&&u8z}@_0cOLLvqy zGV?2gZ4QDe&J7S*+QoE!Wu65N$3_h-JB^lAN85^&xY88AX?;t_w%r-cOWj_E!K}(4 zHuj`5<2|!I{t(jP-tO(n8ke}f5aIZ^Az*mTb+@-RDc_$_IwyKHJ7Qrnmt!WU(pKgm zGBx3L&KvlUZ|K0uv$YH`TE#yb@74GVt=(7&Qcx&b^JOxbDP2XR9dcC*UZFy3j2w$J z04}-kgd{k=P?97LPGM0-@q=^1GQL$vby?_5VvKezL4B`+L@{_=D2zCES1C;N<|4IS zspSBMuG8K*;>JOf3jbQQ0yp}TfB9KIKWD7h@vKMsokkw#^K;ir_}BX50O#=T8AI;0 zV}%og9u7tc^I%M1Z;+(tQH?4F1imeY@2puqrbQQ1;F82chjsd`I{7x<$*lYOOv@=- z@=iE4tl@ZRf1b)3bHQJP7dC#NbKJ#?rHKEYr(P|sP)jaM1XU>RxraRRK%qA!Kw;sN z0Ume)=k8s>Ts`FV_A$KlO?nM5sTa)u6Q!+6Zk;$D4&em-v{kP;?QRc(`zV@(c?X3B z)*R&Y==vW;T5Y~WTyz`n^zEp;QQV~T6`AqLZkB^Af#{oT2&$K0JleMYTe*MVA6pc#Zix&{FtFIYE6df>yp_|Slr=w3Cba!ib6>yAM{c)wtm<14S8eOo zSTnV{0rNEP_=3n5e1KifcyuF9y4KTN859cvH~JEz#A24D@oyP!UvAziMLJ?{ok!LY3;?y=kj$ zKcL5xunczw;Q@56_q++>Cg_|QaNi3YU@sqpLm>Bt257h#Ts+da(Zl6b-?3{KkcPMq17V}Jn#AF$*MxK2+j5;ym#Ik3v*j#UxC+tjeW zt*Q;{I$ku%C@u|+)}nW&Z)7eE&*HBvrr*B*H(aMZ9FScz0_Z7bmGEhc&r+c9p^wGX zLW(k#$aj-oU@F7!nG_L8E&Xx~)>*+q_NfEi4{Q9ZuDqT-u5}w=x6`w-xq8ez%$eed zI(uNlDm$wDB-N3lVoFkU(WhAkzw3JUzu$k&Izi3`6;Ejzvpl(^QA|bbFLW!I{P5!% z@*tFD8b=_D=d#@AV_YO0!dtDYOj^yG)7zPD`(3@azU{&MqX}KO0@|ymOq&*vnz*M1v zcQ_dup|IrAvAuj8{?Zm8lE_6T;hz|lay)`P8Ljqo1J^HL0uFk#WRbKcA!LewkGV zB0z&Gm{+-DAIS*i8e30-raTJjx7+6jd^CXe!IjNNCsRS*Zh)vG-L=w&E5+x`K zhB5EHP?;@qVwC-N2(SI7_y{%QhwV~gam4c0Q>SV!YTkMBQYyDS4_rjj4v`|K0fgu8 zMWyGFQl+O7M~r3}{YB9-dJ$0o)(r&R!bGM76YHvqW*b!c!ChYednl_nhNEX8lcRvj zGtb7x#vVQTdag;yBd$C*4&^;OG@Yx3^wfL3paOeyn_hprB8GH0J?VbGp6jv+@t0- z1PzQ&vu@L#AHqlRa?;oVvs9hCfeoo!a2-ng-0fg6uppI`0lJkQt<%-E?T1hk7;Y_>`CL%{Y~j%*<@sv}vnv!>7yLC4_`Z0f;$(!hww#P8OTFzSYQ)BXeEb z!o^J&CsHW`+2hvT{(&btcIl#pdI!1xyMKPe$t!rNt;qz#9*|?B*Wz7kh>sbnq)*J+ zbnZd5t=EX!-KoOG{PSP8%fDXX4L(A#L~MJ<;&%+f!v|K301(+=ruqWJf25Q*mSBgR;pY?efg+7pI4*(`9)i* z)=<;wbOZZOu*LJzC;~9W?Z+qCp8v0F{=Yxg%KhEUG@Vx4uZ`~|BVrHGbF0X$k{!C- z<~rL_mxMWh$+Qh?(7!8chV#0`vz8>b(`{^`i5sX)W@lwF#W&;ay(1}r9-?vX+vRfs)7q;LvJoY!f~EY2&-yXxy-0J;M@ zo092bCtsN|Efxm2Xm>kZsp;p}1sG6?ZYj^0!uN>F;5r)v0b)o(n>KHr9$vyh+EB0M zxT0DkPXc284yo(jVYA*P@X&X(m#Dr{GiO556x4c%B+r4Y>?@?g{Dl+>huu^HlmM$+P0$iBiR)MIWU>b(OzJNP>dSPv#2C`09$fH3oGfbmcuoO{D4iTyVjg zui4z*Kj@RxrQ@OwE5r1RaXgW-sx5q+9=0#uAYYOdl4U8v!fM7qNZOEM~A;Z|iP z!c@QZ>a%*qIE5m_xM~l^EjitwPBS8C?oRr5uZlqpufFI9CBAW>Y_wa z2PWOXWhKHdq)&B5&_w`3k%^%sALu?)^I`s`glh`LR-bPEk9T)AWTQ?0&{46TpQ3xB z^M{mNA@W{66o0q>#t_!SdS$jdK+K{`g%1c}MMt=HZ|N<_H{Qe4`ryZ1j9m)s(~hN) z9yd-C3}j5qhLf{L-C}yvq>NK!eVsfrQ+Ib+D3$PDBfXu@EQ&LZj16YPGzJ^3Z_sF> zmDR^9G0;{n(g#V59CBq0)pRM!+Zbo>8h|f@}BeM=Muy*47~# zy4kKxRH0?#;^K`vJUgB0hd1!)8<^>Qx2M{=M?k5s8mM3vrjbnkg{8eE#B zGzKby@^xzWCG><679Q30c!nXjzP4v(G8TW0-VV<&Ny-$p$A!)HUbDW__`ZM&G{~&5 z5#68r@6!MtvvN5*exXI!$$3>*Rig=Z^6_J)(AzA6R<`dFvw6M4Ag}oO16ZUg;qXdB zt|xZocvjyEMDTmLp0U5%_WT8yF6Otu;d6hG32xGFa_v+rBQ4dGoQxqwkJDp!+ZL|G zIo1C`zC^xf);J4`3o#=>UaUYS6w8#=1f?Q>;0YGw`r)kMO~;3QGMZ6*TlhGK%mnJM z`uC9I2@M7q7%)9FZNZL$gY%`^dqb6NZ}R=lPpMXWx=f9C_S(>3{CDaNK4IMElq{S; zpb{;q^f^P<{;|w-p(f)0cB??)9lk#1Fgl+BuKS->z*P^AI^vYIS4?Mghb%F5r@ysG>D7`yT zqo+-KPGK+CE&?Yuw(4eVow8`PrNZ6_1>XS4`~d(q>r%H4f@@|k^o7K(m-?wlD8h|k?eA%)R>ke9Hw z5&r)Ej`(~yel2Mdf)TmI`QqzSJExqOk=Sm#k+UIgMscE?C(61^Lra4nvK>>wO{T#WN2sWVGV0wb6w$T(%cwDz(TA>z{1u? z@IM;!Wd6)EHckC^2W!KtDha*p184o)?W!0Q9W5?bV zX9lcXy|;ZvY#D5GpSZ_sewo8P;jPQ*-d0P4W+#J?bCzoKSztrtbEW(OEMf*py!_R5 zXL~SQzdj2e<0oALYte08M7`%Za9qK)(_y`|38h+r@Gf4wc+4MhF07sV+d~`1JuQ4o z!NU@To`Sr?4<|OLdB&Dc439=HZhYwg8CNX4wR%jdf%?HmSa_6YJEpsz$>j+-XVGx%1YuYT#d9$9!`IA|5lqB z-21o%OjUje58vF^ ziJ8VDXSLGHDLF|0h`L*9WaqykT!t!YH_CQUdthq&X2|77a!UTNgcg_UsTvF7E7 zdr^d}T}ngT#$XJ&4xPdeQi(diCBZe65Fato*}K%Sj24fNRB5@mp&QPGY1Ex}QGcEG zUSrqgUaj;_*t};LIPbgANNYMYY@)Ab<@V1mEaO}Xo+#>CO=doxsNThH@A|L>t;+o|BjWGs)N5QPVZyK$wZ`F5H%8DcAUz_6yP9?$zw}>R0D3x>5_W>4 zj+*f5^LK{O_|%Er)QPx?`he7xrl{U;KVTq=LFy(+CkW4@G`1v8KD)rr&v2~nBAO-Z zwqAp{qzfG96>smfxhQHK+~Fm>=BC{6Nv9YD3s$;1u=@x7`qwGJ!nTx}8WBE5$vfQ= zrn4jUk}QvF$rTSs!bF`5Wv&A=Y+^>19m_su6SI7w*YGFOu>kR*eMZ|!+$91bZqcp~ zce+Oqzcl;p%}WAzpD)2+f9RRD_@la(=&L7pm+GOx+dy1 zCtli(w?wk&7#fzCySlGZXt=m-_{M+8hmxN|PB0mS*!k3))q6TmX|@|iWIt*|^+sTP&ebWKRu~NOy-Qh>Qea7$#%SWB9$tcHKJ~^?kv}h zA)5x_%Zc*vD@_BsjVhhOB`86Kq>jnc=crSJ<{f%CF?X!e`r!f`^Y0(QrC)>psDC<# zDv`)Pa3`4|+95*J+@Ze9U;$sw!gc+he>ZKhO~Ly7`xEj4b#$unNo3{$g>i8J`c(nHiZjSE2G+xpSHh&4Xrdjm$4)J6aB_o4@58(LBI~lyhG@f zR%f|{YP-o<15WtuobBc1g&AxRhFNZD=saulS0rVM7Cxtu;44jDHs}4OX5aai*nIj^ zq&YbVX61YO+mIA8y0pIZY8u{obxcKg1v2;6s=0y9j-V4$k9?dxYa(;dQODA?7M-o2 z#qRFqq}lPgbd64&-7}`{q}glR)En-jgRK@jHgtE*_52wprjjrLdP?oYyZyC;{zLj? zw~D322i|;eBcw55wCIU(%=U!HG_)1m-v;wM$(%i({ifaT!|D4y3{6hvU`R(tJcs9z zwmb)`Lmx!^y!_#BcTCM5a~9i1-YA(X`b%#4g=Dz_%N|Us-8+@G^>6PVHMSC6rqmBJ zPuQTU6Q2v4sX^+!h;i3?L}RiZ;v;Qt-Mo4HYcm_^49+Lv2A-I5e(O}P@ROI@TC~_^ zVAgM4<(cgFv%yQ9Oh%Ek8`0t@v0;dyF?tN>7**bI8G^_Yo*cnu5USJ-tL;LTF@pdG zA)!&bWU3FbV(a&HsN!cNIHQ@;py0VeUv2ulvefkE zd;MhwzlF7}vyK-xG6=q3wz1u|5fN6?S}bpRGWn+=KfY0qpX-_ROZLpHqqaVj^kSa& zxmaDr7**k3let8>FhA??{V<+RTE>GUWemUZd;b#@^5uGr{Ro*QcSXn6*G0R>GrO0G z=aP-zaa{?{`3H8Uk|@RJvt4Kv7&QAik`Q1q=PoIY5Lg-#G2mpTKphfKH_on`^1{@+ zwb_l7lmTOU#*EZVd2y^xvd8!Y;(Sp~!c-a!E5WI>x279&@}aeV?@y1v_os()e;97{ z_x?CyHtpGhq_nK&F2Q%WLs1v$52<&1(}i)iwiM(0>L%l^IOdW5w)aoKVTf9kiQq9WgPZ{-6AnYY;1k0t+?ku%^=o^qsQ`R~19o^g{) z@HCW%f((m8!6~-s(y?R5&Yee}YPVsR4CR0iKPjkG+`I09rQ>e?X;xlGbdfvM71BJQ z5hM}0Q9eHR2{7?p-;nKM<<9u`X}PP^d|dy7Yr{?xk7gaiGNiL6!nt5(S*@r9Y&%uY>N^n_7)d3*QT6tTz+#R>u&GPZ8^rr`}sF})K_|HckU#< zsi5&T{82wj?n$i=NhBoz4Q0$yLRJ0KZp^WsdEz0~$i3i?*2`#4m=m|{Mt|~n_(#5a zd5=!%-kfH)M?Z{U?CEIqrEheQQoTK22mCAJYdrduGhYXDQ#0?S)@} zVOPsW^h{Kh<6u2)V-{ZBn-Zs&(QB>}#{OzSg~+~J`Kri$c{|lm2C0o_H?tvjuMzL0 zQ)&Eq!|Y$HCmka;5gLmc9F8x^6iq7gTbtfLfyyaa)zI@CEK<*~&tq_dDk+5RkV4aN zt*iuS|b-P{!Gdz z=&{L{8$B}Oju1YKW+y#&blenr*2p%s3^8QBSAvLC9?T{%_y`YtFIpjEe9 zns@rzeaX)k73%DxcGhje45FAqrlONM!rL>{T&+VgfA{t@Yxs}rMHFmz9(7-buljm= zBO(ySnB#0QD}%$K@#h0#k<`~G0tjaN8PGWE0sS*=`6J?2Sp%%h;qMKHJwuT#w$W~p zEb)@|ynBTT?}i&!rQ!NIF_ea-#_^r?k>D(d4m^~hov zDTC>!Ehn6hiH`19XC0?CIC^>o2+H`)UXh+7qL~qQtSs_*#65f4<&S1d(>OgYQ8eK-QDTP>X@c#4Rc;gJ+l-e| zQwGv&qb|1{q1HjCo)3WY+HUo>9w4I#(lcCr^x&x2uV*6LTWc0Mn%4dKXL{;(>H)4g z<5$jTI1x%GrBN`DMq2_}4aghh;gW;1F%r7*N0D2#K(Qq2Px7AW`AbHmB4u2B#xY}H zxZtoPb1s{|CMG6Arf~JzlJsN*&sT>$``E>`)bu&r(y7Bg4j8r~X5_1e+pmc@jVO49 z3U=(ohi-4#Q!L7ze~+05%*$AR}AP7KOn@a8qd zcUPd+FB`LG>T+RSHW=_e>Qk^tmE+s`!duWUA;sB%-Q}<@h2sOEvtC!&NvmOUGQ3O` z+2vJyL{@gTm&s(3141;nz?-^$%44c-Y`3?rj$IE6wO`bZu045K-Pmj|NM;X<7QeJu zulFgVMZ;xH1GW--xW)~tE+NShNc znF~4};>YH%;1^}EC|B`qUfa7|!uxWxR*|mFH*7p_0hAWVV$)B~G7A#)l7XKeC=L<8 zbP`e4KeqL!V=_D{>JtP!QAnIJZX7$y7zuj3u0`0M?{T+%V(q3eBYTNmj<3Fv2v6oi z!UInsveJpT$mzo4B0dvwCq3d`ohHqjPhXta6c%WM$b$EezHqDG(v9hSYNZj=Aa5Ie zfAQNUj~HTfx<3+?6ZjL z;ygX3u5WLLdlCz9(w67f-zi!a6WwdBnR?_<3`mHNW1ua@1I#4CNlGeg8V#+Yz`>bq z(6#(!^SF@GrV3@`wxxmHUVRXQ{W=;Qy}Q-lQnQ(+cemPY>-j82`m0&E(cPsH-fppDN61e|iH6xAhq(e#rQJ_2E00{GHuZegDPN%#_p#QspSIL7 z+nefXm%F~()3+)K4K-_*4qUnB7X!mRo(E5aO>A>zDH7Yijc@y^88VN6k5QA6vS!Y8 zu>HZ(&ojAw!G=McXIugfq1oS`5?l1C#PGbU+*GQ3sX2i?|otYBs|xHj9YIygUZQ)D0!`>Ue%~re!@Sf4TH&ZdfL(Ro^Unv5%*NlR>v!FOw%Q181URr>^ny?~gyt?X=I+PE852wCOJqx}J$i z94-Qut>px18eZ6C(>0l3tQyvITJg;`Z^8Bi;J@$a!Mq}h-tZ6rgm0~g^)Z^@oFWi0 z?D|d4Y~cOU(_fHNH#|?Mol)@SmpQcR67wwpFDLPpi7YtjGnSv+n8qRcKsLtIB=mbl zg$=#t-F}pqswjVw1jj{@PO6C!_wGo6kuu+=e3SCe+k3g<4p5HLRrM|%uY@Cpt?g#+ zt7A8Bh}zK}mACUJv%S?Dtq-klC%$uI1G4pv{rMFEkUWPLG3NEVsBK>;qsH(pFQc{9 zZSVr4dW|2Q`8AW_us+kg)t{ zo66|m#}~~wS#M~nG@9gNCb2!dBeR#9?JlAH0>ZC?^tapl zxS?dj>)#6a8$#4CS2vAhnA4*f@>xP@^pT?aLlu2*v2Dl`N<5-B9#srH8%N_t$H%N3 zwx>%tVLY0*ZMJ-ASD)zF-NqdaC$OL|O}0vZ2&ya{l1aO2Kd$y(ztQ5q%0PC71U9bN z=<%m3xMBLGy3>mi%a{AE{|LxpQQl_*>(C)?ZQMw@V;7=UX>lGt@Lp%dcVGH zHG~q*OeH4PYiKj&e`%*(Z#<$F#3l_Xlvlf3$n5rh29BH~%H^5X>AL#~jj#jp^QdH$ ze9N2O{L8C`Qhl%vxpmMaE_pYKUzq@ef@gk`&W9(o7=>GGpSK#}te;azlZD~&{gP^u zF?PnR;YUN>Y-D zwk1Fei9`EMa6A20l#WNEo?;zWJaOd_*c_rqHHnm2u)nc0%NO)U8ud&YmGNsSVf^)9 z=!rVhrrvn*=Z*IFau~PzcC|qFrMam&#Pooes*Te~%_=_A8A07!G!ffv>lYLr{+{gW z9q||&+jKcmREcD6jul@W3wzBz`4(BJhZKK0FM9yE({UHTJ7H;SiZOZ4tY%&weC-!Ij$L1|t zu$I0dLoUR^`EX1McZ^tzsPW^QvtK>j1 znz?Q4crAb-{`FT&UjgjoGvZ!q3|p}c1r)fLr?IalcDIZ7 zB=|a7vBi1G63{hw3d%?AXG&(zZ z_Cx($mm}2L6YOQsSHgz%-zuZiW$QaG3aqi*_ZK^+Jed?FSxT0-RQ=7HH{8?#h4N#A zcTf5?b+34R9rMmgI7p%<<$pQQ%^bQCX!wvn-T(-N%)fUyNt7+|DgYuY*=*ldS+8HC z&yl~jw$@vq(r1?2`!5;eSZ?`Ut_QVyQiRCp=#^kBK&QS>m{%k7hdB!93Xl47s@bxm z7)B{w!-6B`J63v|<6oZ!Dn5|5o|VNFZrH9PwenV~(sWoNy||+K2iZx}zV=&)e+1|< z&EwXac`?atg8GELWJ*yh93B?6&!VR`$Xv|c4=6AdgCHqCPWxWV`^sQgr%uAyb}Eh8 znqg`k#u;=|(-9b1#nbK>#H89H09<|$27fpf9A*Gq1eJ9!JR4iQeR$#(l!^9-#~sH% zkOxXJ42TxNc8riA=L3xOHmN4g=fSu4A=6%(RB&HO?6dfK^JWD*LHRSs$lK;a6g%Tc z+$ieje+jIvG}79V7+YQ#Nt~NHN97Toh}anQy@h1*oL0;QisWSFA`}ERN)2$HQMfXW z{Qvbxza-s7x8ScYJ>6ii^1r_}G#mcz@4p0U{P!>a+cofiS{wl5>Q6y+$wx=u4Wx!b zL=l$l2s%0W;>yJ~#+X+N;;G5oCgn*w#D-QcE0A%L6o_)2nUkQ=#4P$x9f)4*Em}lK z1JFRVKcHqj~Y-MRWuOelavWv-;pi+T=X|kdN4Z67Z*7E13)A?U6#< z)ar`yE_BRhQxU3NoCCP?bO>YY>X>~W(+&~D)RqiGff{Lg#;uKkVoWfYjxi#_l7FJ2 z59HUit0O0kZk^jAkJzMX)3k5zsasFs0V`)ysOwRfP~_eMix!v{i^NoLQMDA&poTR7 z&?KsWg+Z(B7o<2-%i9d2b1e051?*n7#1V>Rv+y#Nynt)im0nZ_Ke>=O%RGc_w8w(d zU>Ku`33?W`8bH;b@6GB2a%Fq8mg6i`8Pz1Pq@v0!OXFL$6@lG%daB|)5MESC*dXMm z%4KYhwO&r)Ruop*B7rTv4mKeCYqo_LoJx?Ji08WuK0#$f7}>w9pI0EX=d4_?{AYhH zT|-ezy|E!ob6(uwPH?jk*8dvbn_p#Y^>oc}RI;XYkz-B6^;rNKw zQG|ac$&47K-t?V1kmH^nkI1C2?02)FgC;0uJp)<<5^r4N%DLFW<31MsLzDY3+;4Eu zCj>83pWR^fs|FXSQ?^sAc|Vz%FBsE4~TL?v$#%sx9Ab0GR*n zw)@vpLyHJ$jqjY=lmvo$%BUkOvHqfocpYBWnqLwy5ZYU`T;3ictmrjXW(E(3SXnY#PZ$z_TRH%b;D-Op2D} zho|RzUVMn(;e{(a$|5V4Nlj#B7Qd_L%bTvuAB}$jWFqF7Iz2ZaX4yNdQ&*Ty(U0qc zB5fY<%f*L8A_NE}F&SlxsKKgP%Vk6lP$!uajQ;10(G}@X#5U{r%LucLP3-!KWJ3J| ze=LLhRSJ{8UhDocUd`Eg+@&rn@L+aDXee4QpKB$}s?sPu{oJCd*{v7SyGP0(T96`g z#p`f8_P*SGZh3RO<#xq_xHMunChY|{Fh zu1>xY!`!wudt zz(tGdy6=j%gZWP}B)}qq*E;n&JUyAHi;DlG?E1UvO4ce^BI4K%4a%=Jwwy5FgLa)} z1~yJ;l7&x=>2uSrnj)mWS-(NDMA~-!2DPZgxtYQl3*SrI_xA-ea7h=a?FW6P^yd@C z(%yV|L%)augfWSBayoASZB6En0%v$Nr-T0+?)->hJ~_@nI>k9T*3mdY`VFLGk|23v z^f~3cG&GYbO(eqK?%D4HaX47}(LXZ`u)f7QnIdew5xxZ!$-u>nPd89pzleW^T-ZHa zCclP0(ao3*Z31$^l%8R;5UFhXap-A^2f)00rnxu}JF8hs+$?npp!TOq4KV7YW^=P7 zR@5ADHopEYBF5h#X|XSv-?Sy|hlL=DB_JSbePXCRQrd>FvbX#+ZyAAbr6imSpY&Ng zo9iNVd}vVt(S9e{X?CFN0l(^sF*8vbl%4eLCUYb2ycNY~?Jv?lLKwb2t;x5a)Vy0; z3ft-I(TQe+;_&DL+(gnI#4KMG;qyv^Zhi)y<$(yIXj3B28k@}r(&t{03B^OoW#*Y4*q0j=?1&#)!>twaqD zLQrx_%BZmSmKIM60{X*B*cNwMY;%4u#VfbBmIMN#m(`Q`l$OZ*hvpaPn;f#SvigeY z=u(snth&eLF*kgcS0v~$Krs2{e+YR4lzUUJ` zg|0Q-t-D`e#|UDrr|sK{QK)7n<1!hpWFsa)VomJ?s|eecyZ18#m=HGggU|F54+vEg zOFm6M+5{GrwV+JB23VE!dm-jAD|e8w(<%-S*OUR8A-Ld}rs@W?2A_9<)^*s>nKYIh z%g6B}W_pshDCXi(qcGB2y~Rj(KwmRrdfdQpoVz`7cE`s=+Qvm)xUbbu_B~g3H&e%K z1?AbNNsq~a3s>lk@(F|9kaJsKSNDx5;u5v1rFnaM<5__3b&d(+qvG4^>}(MFu{ql% zb9X7sEzm_S^B+cMr*hX1<7{-TD_)tDw_8>szQtPXL6M>6Wej?w-7FTS?cV!`fCWTP z^b=@=T)`ZAhavwr*U??z zcXf*}eaUE2&|M#tznwY5%f<|$@?g|+2@@3YEiyJIXk-f!61G3Ayh{#*gvvTOR>|US zG?)k3X8o`|w$5f27Mfl@P&NlO^o)F9d7s+ip7IxGb{S3f+IipOuo|RqAJ0yz(?-dI z4{1}+TI=F&KDXacrKR%$cu2?xSlf9)9ER`=qFks-&YnLx&G2E2R_$exZ=&b~xSq4i z=ViGC_HR&qIz+7Y(`&zD9qd{;1dZaIV52E|HIT5$=5ujz|k{dQS=9;^r#CqjCCvrPKmu*q;6-U^cfgv9`7&ou?qOHVStpP#Qnh>@1 z)rln4J4#FTta16(okH_A3HO%S6LrIok{R7Bc%qp>vr+Ezobb1`3*S!5-A{$Ub%FCT z(%#wwE>(BH&*y*G``e`r`G6l4H&e~|c+-D;yhLdW0EI?xJ{SwPFf+Y|FoGd`EIMGG zTfsdUAMb3Wl@Vn(WS=I3YlSKjGzjnG6UM1L?(gs4*rH#X@Q=XWe1^l<&Ye0mQ!Zrp5WrY@Vi9a8;@hpaBzm`Elwen| zy6^~_6Eg@V(i+5Ij0+sqswu`)|91yD!cZLKJy`m?^$xFGQ6>R~xRJ5a*a-@3Q@e4P z0oG?gMv+CnuT(h`PdAp%HbCon!zBdwi7aL4z++@w&y1Je0_b$!k-~U#J|(g`Nof_u zkrkp?vDq$;vI=rQN%Je1*)Z}y95M4;V49#P9Afon+gHMKCmwBHV-RqVFxfBcw^a@^ zWVl-+D<-=ZF{nakSczbXc#K*7+f4_K!n)K814uzwC6$nU+){zBAr0(t3x0+BJnpnd zH>yh!z$SwtZHHZHK8NX;6)? zER!VHA4mxu3U(jVu8vP3L%U?5Yn>$FbxP-RI@A3CLJ0^7vaAkeW9&qhp00(XhR%{m zmVY^@I~~WHLV}Uc_>Du!U!}X_$ewAFGaLPRkaLGdMd;IhEI;I!+c$64k8}fM$E%5c znm9-NfOZ|Y^Flq6%#%ju!SvWt9;R-j!L*AjUtC_xb^pk{C=}>e zLIv`FiaX>$e)e~G3myj7l-dp@sc4d8QdiYT2naHeK&DmTS1Ttn=I9$wS9#+nirKasS}c8%GS)TQI)O z;EfS&%{JQFbY6F;?PU9kkYSzgMQqJzyQ5=X?@KF=7A*~!yJFiLjp@(6=fApPsr_fK z%)DM0Se?W5Jc8?0zpU|!d1MV|B7uzSswl`A&kQ)U2CtSP7ed$JhV|Qb@g06vkgQ>T4*EBolE(W!*M^mb0wKRA&MkaL zK)~~-Pi?KN`jjOwI=suYi_uiI9U@e{29y6FBkzAp=2al)%bl=T{elb^QL~;sdIY`s z?aGT@?W!rdjo#A;wrk`m=sodzq+p-4#N) z+S@e;?Dq5Xd%MQt6{#>~CnRL3EoiiAwg2L;Tnp3@iFtY&pGWlsXCYN?4#1Y?UCSLF zacLnPBZK2};DB*i0&`lQKYy;(U<3L>GNPhizwc%DH}m0Kf+?6H0FqNDPqNxNYMY)% z<2C%swH1m|?+zNH8`n~d&`);iHM@o4buDYWTXaZ7#PhQYr>t19qMwGv=+T|Kbm@YZ zigvOKS~9vlEMf%BDDah)m2t~8)zvo-ZA%fhlfJ$_z(P`zpG^@j*Ybe>&Yfe?=$SHQ zfRRxtxQSE0QfwH39@^ov=FWv@cLA5h^_Z^P)xe-M`|I7icMrBLdq>abDK>d~ zzXhVbVagD%f~dGSkNz>WYS%6)F7B}Y(3lXohX4Hjmt3gJ3%sc&?b*|}yIaBAu`5=5 zOP%tZ)RtuDtsma*=;60(`N$3|io$GE%UWxOem>V&vG-Klo!6(EumPJIG_8)Q->mf@ zEp2M7+p48Ym*yukAQ=Bh97sWL&J1U&I!vZ0a>%xf3xR6)XFEDq<2G~%m-Xp9+lCKq zQ9RSl&FF4R1i@YjX06&fV8$TM7}wh|rLpexShvvda1AZ3Z2{(%mR&g$a5^w|oN1^R z37MZRUzWq|{@9&t2k`V@2Wk0BT;hY7WAIqVUl{DU(Rj3toVdx!Q7qMGpMA=uGRH(r znev3p>5qfkfB*gWg0J#raiY7i@ukpEhy3bI6pFn~|6XC0%+~QgDQ*npb{sfx0C?+u z#!0vxq_wQzwXQV5U)_&UFugZ!eYF?g&^=K6`ms)}-wutj-&`{^bV+b22>)$ zWmEL8<@#>PM)39oT+?FRiWmqnLgyqz9QfIwVhbNewP{hWLw|1ERj_4ZaN@MZ*Cw)Wnv&PH6SrIIx#RhF*zVIFfuV9ARr(hARr(hARr(hARr(hARr(h zARr(hARr(hARr(hARr(hARr(B00000000000ZB~{0002J3IG6joTZg{G?m>R$9WA= z5tUM@P$8rcDe^6|44s2xJmxZm6hdi`GL=F}gD680LUYQLA%q5FC^Ka`T=G`J)oH)? z{&WAkd#%0pdY)(R-|+o@_Ih|(SXg%Pu&{8>tz2{SFP1-pgFLhTXR@>Gv-h#*m|NM% zyWeMgSXgv_@8_MH{+$1f!}&krRFd`g;GUa0|7@BYV=mockH7b*{#hqf-B!~K_UJ-jram7=OMm#iGX9?bUJ1$-b4)%TbE|(I7z%9Vewe=E{X_sl+4d6HKi=?3v_49ZrLe=D}1kyJQ# zOM%p(h>E;tk({H7q@Has8`J8Gh7zeu34P)Xs;e(6KHPmL<6AFvRr{pH5Hi}bMzx3&qSZ+Wtd^s^gZjOesU z>)oWU-Pr7z%=(f;jr6Amzwm*qpe!}gr#cE!WG%I6>ZD(Fa2EM<@=U9fzBORbTKM2Y ziU#Rl1F@5wwpaBvi4K}*<2PJjGoeZJ(8Bsnu}{~*`cgUbkjxn95gK5xfN?Wf6Obj(f0TD8%EO zoF36r4|elvS5RB@h_3ob6BQ_5c3Pk4s}IIb`r)>X21I8AY_mvSll0zz=xvA-J_A+G z0ft0(LzG{X_1h$}m*~G2JQ@D$#w+#`A1ILQZg-|DQHUQD)OxY0;yeY3OkXHCJR`PL zw1Wb5rax2&oZmCRzLko8rcYEnU*TEQA4(;DQ9-{OeWj|E3K^zvMqp8R?4GsR2p^dK z86o9WvszTJ5sol@G=dK0m!Nl(5%JR)L!Z_%Xi~<+S7SsS_?_99_-l-x;&(KGa#Or!@?wf> zK~>-XUSUSKF@wsYk=hmyGw?F`F+)c@e`G6OTMfH_Q=Jefn{nK?Ds)tqo; zj%v+7ai@H9!k0O`>Jq1Q_$>%$7Fc3mQYW?F0=-P$Ebw%sFP-w#DR)O6E5fN2ai;2BYxHd4%H-b`$~{icSQBl@J+MXlp@<7T|FVTS^FG+2 zc-C~z8ar|?>~N{vW$vUMxgU1;|NcLt&~kgLpgY4r_QhAdwtq8V!^;0~`yc}!?!?vI zYGI&hQP)! z7$Z=cO^$Re7=iqciJYwwBalvgAkbkn0@?iA7mf>#VDNaR?#9+(`2RJ_Dt==amlqE= z+nW!geTa5nYsoO^MIG_i-weUF7Ix!-O0B)_0yU1-I@>VEKehqlUk_T&Aqwvqm#erUy%)DCC& zVfwVhqh*vnjBZO{?PK&JdBF{p{TF)?uDmchO{NzIxV+3uYkClOF2m!Edk<=^Mb*4n z(1UMmvHnf>x=~TD()vrc8wIQ4h4+5#LXp7krHx@-FpA9|O_u0Fao)gGK}9D{H}AJn zaqR?q%&mJO^E*+mA=ZB)qXU(fdZla)I$(0PmNmM)9VrK^S0`?0$E4~(?>6_(cr?{K z82b1V*b>)WK04b5?xxi=UcEMKd!NL`ezg@(f;&c@_O-z37hk}HVhhryG$v-kKEgk6 z%h-#yW|X#4L#cAj$a7ZtP#oTbnq@Eb`MVq8Z=_)EpxOw00D?Ubs(*W)><9SW}@)BW0Okw!VTt%+|GGb@ATFLm`;&r(G8Eu6X@V_L?c?|I`R#6{+> zXli!LBZuSW;B?k}8)7Ho?%!cr&R+B(X7Wx`u zk6ManLM`&e$Iv(TU~ER&Vt6VY3So-gRhnsN$zUIpUwIcD*@30wzf#cv&d1w_Q?NLa zqmF+d8R=`Zl@rF3&=DqgPMRwT#rrdoQYCJ~!8ZAAq0ucE`1+pb;!DI16lS77r zk?61Y%G-JG9F&AL!`-waP`g(-S8O~CjrW3+Te3oNL~wGkqjv~+1G#r6=$(ad*wc&S zl4r15V$b#`!ofIxdg%%Im8apj$`YcRf^asgI!#RD6!cx&ewK8e0K4a*5>d_LPz{ec zzPIogDAgr71DZ$SERd&mqSYTx^*rfGL4Fu)6w0_McLY_oMeFsx`#?D=WhS%K8!|_1 z_jacrMx_vUn)0PX2;SOSdEkl{T-&XcZm1r_m(Pw}mB&0_X`kr0y!ZfeT|L{zR=T0k zQ!e#whzp#P9^T1ecZTt+ym_&g_k&WaT@k3b54RPQJDdj`FmF0cnmxxJYdM1H*}-<$ z=P*foS!ILlB`$secj>req`c(lO-mS_NTwHLnxj_(Qd$kB_}Cj0!@ryc`Ju3a&mP9` z3o|guYNBH1tTfN05d~50x7|%k4B*mkHv22er}FRaa6^LsA8F6>JB_o*(hS=u87(dsfzA-a%kzOUsm8I1NU5?c8$m4 j@SWw^KV~6<*NdND_YPSGLze^UrGlLJ=MVC4C->)8GtQdN literal 0 HcmV?d00001 diff --git a/scripts/ramp_3ms.mat b/scripts/ramp_3ms.mat new file mode 100644 index 0000000000000000000000000000000000000000..5efb23f32b8e9507b34142c594b3c0a176a7cc3e GIT binary patch literal 4232 zcmV;35O?oQK~zjZLLfCRFd$7qR4ry{Y-KDUP;6mzW^ZzBIv`L(S4mDbG%O%Pa%Ew3 zWn>_4ZaN@MZ*Cw)Wnv&PH6SrIIx#RhGc+JFFfuV9ARr(hARr(hARr(hARr(hARr(h zARr(hARr(hARr(hARr(hARr(B00000000000ZB~{000005C8ypoTXNII8@*NAF^dB zSxTh^ZBn8IMfYXTl6~LCjA2F+V=YTksU$5@XtTAcRHB3rvP7k9sYpeZghbL_>Zkdf z=llET```CI&wcK5&z$o<=l$B|PNC6gvC=e}(8N<@;^wFQdkK|;G;d}oQ*h!r ziH!YS6H23*{+&N%;!d0V|8<1^OP%F>|33d-6H{rvlWG6CMtytYZ++{f{(XL;{N^}V zNQV6(?*dqEApfpgLlpmxPesem+SYUzDg67p45mK&v-VNy3^H#P%=5HPe}5`Z=Ffv+ zjn)onQ3bM25nd}D)RtaSBI_3*HY2D@*JKf?vl!{(hmJbET|(+DMQOxkgMhs%r0z1b z_tk!iTc}FvufVTMqmM1T)W|;SIPiLh%JUoQWWSY|HDgMgH{U9M3%5z(u~{m){g!(uXcGHDxcU73z|H z^w4kX?(oG!kMyO76#-}YN(J>vfBM+EfAH$gtNNr*1C*arcirG%K>9U+P{lSOiE#tc zw;?QE+{u4@&XDwPh~r;`JkvHB5gd&0a+;+W{gV;F!x#&eAHS`1+L+*CjLjx5SpxEVE+PmezN^^pzIT*qXGi}<<39cKFqbOTDeb+{U??#NQcMpE4Y(a3gfTqjYx!KJY z1aC{kgbvnmA}k5+mZ&|fzIC~R6~W&M5?8hRV+M4iVgXxK8g958<2v0UBo*7W_cApL5%LX17lTzwl*bv@q!MEm?|J4wuDbR4858+;-F$jc(uck%?1-Y+Yx^4u=+du)a({J z!m}OP_9g`VUSv=BwnvoAwGCa7_JntP%*$?`bGN~s@NbXm6`#j0C^`@wIDk_oSGPah zf#|^jqUJi8{&fyS7Y?}A<{54^*OBPM5w`5z8&?K95}i2W%S1KVT1TQ6N1WJM_wDOU zC!!lCte-mE(7wfq=*J0Nr>C7St8^kda>Aa!+sZlW%=1$sQ4xj>p(-J;^-f;T+9x#0G2-=!(VE<|@OaElg=9Tj)QGM@fi z@q4nW;}eD}ns_>NMRuyqvfDRZ5z5n}D~uP0@t+iNgB(wnZWzo@Ym0Vw!#_NIy5XR1 z={3(hH+b@N>IU_O&NzJmcTD2x)g4dWCb=zicE?4YZru?vut;4b*PZCs9r9bvgx~#f z$1qRF9;grqcvWxbfulSQbFED4`}do?E#VbxArIh^dS1CS)kC*W44Ec>4E*b}ze#FT)e!JRf+XBPiie z??+EK^ZdX-knpqW)S}{k(*PwT6`W8vOFI#5yK;tiMKpIF=4IoMsy^YNqohG;`o=a z*Qc4pUrcavJmSy&ujGF8`)MSRP`H_Jzn>c?2bzbw?rcDzw6XAvK>@Yl66{VWLb z{LBKc4JmK5@_g+LT|eo?U$wo7zr7))zVK+jk2mqTH+sKnj6Xl@P5kbS!g3SehsEB+ z_uh!#GUCa;RB=3&ng@zvq^5SF{89!#bgzmmCON`JiK#fh%04!fczMD`qU78y@;BJ{%*z!v=&dW-Bwny_kC!iOsHN)}{1M=g zoZ*1q$L`t(RSwA;4xXs^%#Qx z#sBtw_5o=$mnhTUBhW3#GhlulLDSTpah-P~u%cJqxzIiWO6%4||C$lpuey3qCU*ok zU_bxHieQZ)*hUqnW%j?r2a&CMwQ28=j0^RK8{WaVDd(v0r$J0nFmX4_ z7=%~riG&%ZgBU%iGv(mu09Ya+lk3w4&|d#yu!9;vzQuWqPd)wkEiN%xFuosmjf@N| zRr)bU#9()6{afs7&f0S@@GYKSvlhB7`4$Jt`Z72-`Vj8B>sE$MAHp)3dNHGK(7J8g zJ@KSB=x^D{^;COpBc3kpMs9t~?LT#0_|C0go4C0P7x`JIV#T{aEz)lty80SN zKc@)2w0eyKp?8d&-cHO|eEw6*(oRJAeo0ad?ttMT5xK7=uV7!Oa8gL{6@Gqd`dDcF z5?!1z!3Sqv;AeoG;EADjqzr!9G^*2%r!wclXB>Wx-uKHYM!MUuDn!URQ>zU-xdARJ zM_MsyUKDkF;2A__&1KPyo}pXf&g+fYPZ2PseC!V2Q@lCe@mPxW1k2LpP4a46AjAzT zU$v|S3;N$XOgZ)#FD~Uqy%}!?7}9Gb{F?Dd*u2KFtqHNOw%N3rG(mVzikQmQM$}Jl zbDzU(z%d(E@pY?uZ0VMbIjR2$C-23jIjwt$g1(4&bB#K91*cxu)qH@eed65?I`=Uo z!yS$?zK4}t$HkiMYY})E%ul(Y%Y$=w>Jhlt3D1mMNPAztH zF$BD`xV6SbNSFB#>NNj0ZVg>=ixn<}PLP59@1Hl}o^Cv3FkXO6k+7Ruq6L_jaR1

jHaQ2LziISqoQpV|HOc>!*#)?k1p13>pGVE8^=S#MbJ)=-CS9p_7Teb7h>VMOG%3v!7P7u4X{s`XgGXayqm!L}`)5r;*0jNXv3Lg>R-s zY#XtYcwKb&;_rqu)E%{a;&>(%;%EQa@I5L8p%rWE&uuym)vT>1&fdrHiIe(q+%FkA zKb!gAMjnO!-bi1i6G=$0&N_Ux{0K6mUo|K89>%p|Zm+KMA@Cn}Qu7;3M9Mw8b&dND z;In;&W#y^_#M$XKO^(bA#lqxA zk^S-^kTl5-A6*!X6TAK6at(tZ9wX-x&f$W86JwgFav*|dl)YJed<%@;O34{Z`{Ocq zL8JFsKLiwvipCjj!kF{j>a*WDP-%DX4!Xm}frAB?Hzxa_E%RsFk}z+)zSvTm&thRw zoOjSmx))NP8GJm?W@4+d;hL#Y3@lOD7_X+~iMqUy>3WAf@Wk9h%=)_<;tCz<&q7_X zT_{eaD_NPLAlxRFir+&jG7{q@<(_+d(ey=hVqfw&-xtj8!SM zfk3YI`@Y-O;E&YnwA*fl5=w!7@tg&=l&)3v?KX!SxBQ*O8Z%s*roOj0-30C-p-j`6 z8?Z@sz(?=05qf@!D!(@|gl=TR+?(V2sQHw%m@3nQ{7%K%;G??Om#Xj6+(|(^DN6P0 zT^&SAh0y4Q>rhp>_gPVs7KHAmEcs_l1IeTM1!mM56y{dPMIKrOU3(|lGJ%y?y1iI& zdW;$@omb?%)l{AMzP+;RWS9y99B+-seO-dtAM?|@wlBiDzLL%`?FIOCW1vh|L=nY3 zQS-M9&BLp#xW&1A@+es=V%TUc2f+-@k)Q?{EVZZgE!jN{KeRrsnc5@{{|a|I=9mz~ zOJZl8jpT#W=JWC=n!iy6bE+HVbBC#QsU3$mKWnE7`7-VJW@eP9c{YS8UM@d z@+iUZY%BPRZn2Icpr@S3M=h2gQQ&L`^wTebr)yI+S>Qlg`#;UR45j!KWZ(V`%#C~)pSzqf=_IRl|N zD~@m}t<^DExvl|J`o}xhzq^D`FPus@vC@O7e8cCaYRX*d?o`2+TW&#Aqlf>m?_R;w zGR^U&Mmj;%Gp%oPln1y}W*J+^H8qHmo0RS+8XZJwoH1G3C>cbxXfpzPCib%~t$3_# z7(}Tq8vS88H;DR?U0CiiDTu028pyMsm^U|ZURsE75Oqj8GCO~SOC_bmej4xRQjTKt zGlkx8sqh~)5u01N)c&7RN8UV~7=JY&H&D%`+Uw%(Q6*gJ^3(Fo+w!DM#4#>4YuZ6U>HwD-OuJC9Vi%VxnP0(Fh~QE^D|KG!a=Fx*vg;8n z&ct|QPtsL7m-?q&CWOzCODReWkA+!qDe>76)jtfmRM*_bQzvw|)Xfz)U947fsUHIS e_RB5fQiXE0YyK#4snoGoDnqlm)c*nf+pUbthfp~H literal 0 HcmV?d00001 diff --git a/scripts/ramp_3ms_raw.mat b/scripts/ramp_3ms_raw.mat new file mode 100644 index 0000000000000000000000000000000000000000..9ceaf2b0c527eb906f8ae22cf89dab09ed280edd GIT binary patch literal 13517 zcma)?LwF?&vxH;Y6Wg|Jn-kl%ZD(S0V%y2Ywr!sH#QDF?-QQaEwl`HzcfATC>hdC@ z#2k#w#0nzn3>MaQ=JdqMcE)ZNjxP4R#L8k?GK!q+^u%H==EiR3X2gyTyu=EQ4#bk? zCd4cp#4Mb=%pANdY{ab0tSrR;Pl5cOKxGt=|C=S)K|s86mp8oGwf6d%@6=e+q|2yN z>cgyb2?WwYrK|#uicMc;D>8^tCPvO>L@c71YVbx>Cc(TH}Fj32WSwQTo+I` zEnOZM2=t62in#c>nTvMFA^1U!Jn8>8Vg8<|9a!IbjDzSh|(4cd=h%*>!WRYTn-$lUBQA4gqof=|+6u8Rh z+~7V@1wMDZK&d-bDqug zf5VKutXiCS3l27xk98k?6&&IKpVjPgrkv=%=8PU;&7OV3^P$Yb*nkGRF8EXHRyKf$ zW3-ZgbV-!Js*_c%{Y0ML+ytW<7k3lWLcGXFCzjv9Tu9F}HVp@sw=k->VgL(=*&J zU0G_(WsqYCdlx4a?T4A9r>3);5vKnf>L7g0Or%ai{m8X))V`TYm76N2XWO~4T#pef zB6>cU1?QoDDS6nt2}exUpSg?1SaDHulM#TAsQ*p|uNx7CP~`+jKt&@Q$Q&6H2RI5m+oa1hklRW&mZ+rB}dxV#DylM z&dR%xjN8m99yyuKw%PE*c_)17nINQir`{&Q2d=g5rZZt$F95V$PXImwweii#k+PRq zu}Kk31(<>oHmlNqzJwDbnG1He6HIllZt1;A16`TEHdcd83BMMQ&yuidFU7Hw1a@h! zGLoq_DR{L!57&FH<9P9!6of3Iz$dUhpYAkEA#XZngvOs$6FslRpcD@vb;r|S&y4BJ zqGd+Vm5$5c71YU7dS?qJ#T?vJYi)&Dx8>~uNCKtQ8Jj(eJiN2Sppz06>3Z3vMj07k`7vikB7-g${!;K?!T&@VoBdS@w6VeMpu_SWRH$!bZV_nYE=7u_v z-@X_Xmk|RL>9c#Grduz0zUc42R!k8-@oz_1W~alJuLb$t9IskR zp0W8`oO90M9<|i9-o6=?l0En+TXx>>=V1w}E%HI}mzb;+Y9;btNef(`@CYyKI|CnA z0bb`B(BVBw{YUFU!wy$7f%`=trv}pOWDFORfPc^neNUe2`FQ`tGy>e-{6t(}*B}VL zq9YyHzGu7Mu?gSh95Kh#j-6+Z7-Yn}tTrfoEIiNvcGYG`+v z6EPz0=u(~d#*~q?D6Q}YWFWvm*CMGaqy|&wXIzlkHSMdReLQVi>+F|h`Ej$;!ejP# zU}xwb2gUi!dA92#t%ucBb6<>$jk0x-ovC>%H)9H@(Kw?(RYcVJ8VbQUM~`TL;b9gI zcfXEyJ3;Q-3oyC5s)I)CMC~ZiV-#5{{UP=+MCG0BolVNq3cIvgm9UneV7~W5wK^?_ z==of#i9Gy6OeMgGIi(dse-?eY%}Vh%psDj_L2emfGG|fANpZJ{ZJqy#B;eb_aQs#Z zA;&iiU*`MA=x%1W{JvgRxDL_c22UC64?XZY%Tab7yW-(x)sy0GORy6)fLg?Am95ps zY)^RCw@ri{(4GzB#%K0{NA;Rz@ie@{&Jg>OZi(E5g` z7e&{yadIeOt{+!B7%^1si_&qWWNU^$onv>>ZaJNwZx^q_alhUNv^qJDiFhfsg9o5S zmT&Q(86a{efsJ54YlFTz8&}GugUG!&LuU|ywB%!T@3hJwx)_Y8UvsH~2j(Z2|CQ!o zF0>LTpj!#*KmLyT%*Bl}-loGYT`n43FL;cs;QeRr1}~hz%l)Ow>r^R4P>NYM2&lzv zpu9v@SQ7{D+S1o&owF5Q5ONo6=D>A2{VR=QF9R;0FF|!OWgWu&F&{6qo6T(MdLHZ^ zOsRghQ{4%aKxb66x3fH10*i<+>~9XoqNear5KLGLc`m?Qb?RLXnM2d+++nfza&pn> z&{e@BHQ=4PCtQQdVc^&lY;P6GMQO8nEk{%vQ9&5Y2DU#zWy-<$JbmxekQ)yet%I5*Lc+VXz0+*8tb;%(xtZ;(U& zH%YKY!E-){ruY0(Q0X&oM?qe+TMWIygF2b9@igdo*nQB^C3s$W!FOJ;8hL~5;Jrh= z61k9wKln)YuG~9DwX41iVgj~nK%P3f{^q(KAM+syawslj*;2k6~pbl zG$$#|KJVh+!u@T54~EWmGxh^HQo!I@ewW2Yg({Ulj%nzE8|LA z?d(Y5^FjD;jn^m1qy%%8$d#F4SZU#}O~;vKD?PS5 zrSX{?IZp5+!c{3a1q)*81D%3#JEW`WI?igCtYLK{wKZY~qyn}KA;T7~$oKazp=Gd{ zBCmgMz&i=oj@dl+oPa(2;yX^Vb|Kt$b>>i`xF}Hq$L(w?&koiKW2b{}1+z{P-yFoVhz8;{r0JiZiBCU?e zLgffGghpteFd6==i-_kYv;jzbwOKdJe#qttKim9XoW~m@2Z7-#_6#q={Q;z6mJ5V} z^D9JYvD`t8)Ka0PS>Vga?!aRxK4ELpc6$H=^msbk#s1emNHmd&wS`>*3dX~Odrc3S z7aZF~4>7^u0p_TBh z$U(#8lAnR%esi<6Iv@mnTD`~fnaLxBN>xfR&%_z3pP&9DE|g9BLB9v{d*Iipp+mWC zoipf_la%8PIOCl%#ZOvY9=u&1LF{6p8)*NBEKS4+cU+m~1oQyxe zl6=|k10d!auPd|`cT0POh0iP|R({{<8j}a)*~8rJVqTpvwvg+X=OiXnH(^??aq)0+ zvrzu1UGeZS==Wty;JMR7s*HadIe}&|SWq3G~f9AW{WP)oN3g7}dK>UWFR= z;x!=Ff;B~B@6Zm;McL$MTf%Is+TL<0>0NDYc`+G|b}pK3yt;I1*H*f@+2O~Gq(X3$ z_mtSQ<|U9`I3n?Rw>NS5whK)wyvLk>U?pz*UV+RTMMl8F&qPQflF!Zq+b!8sNdDS4 z=r{1-;7$~VfE*1ma#eg`(7w^<2SCp2*L7#0%NL78DT@A}1z0>78FrsM!lYJ=R{MPl zp9YBf0p4*@=S!S+t#k4^2`SEG07{2#Ol7~%eR=V372YrRLSJ;2)T65tk=|08i&G6R z?tHywZQk(+)TgLc-4PdqE!*_kjbUzhHbnYTn02JY-af346zHJ#&e-a5+Ayz;WkA|BUa+pH zj7(l=P(8|ib4e%EB=tFUwK;`>MqLzo-=-p9eJ=2~%~+&CBLJft zbYrkO|MHct7&PX_3FdXN+@Q4UiH5Ncn2w2x^kD`G;nXj#u1tDRD3xm)I5>#hR?VD` zI?(~5zIjkj;{9RTD;`%4FoO@mg0f{;6z1_(8!^r*C!c$3lzkggc4;vzj)?GF@p4xW zcTQ6OL{Ra3`77AI)lV)>j}2Z+17iy_;0OsADPu5&O zn0~;g1OlKOfhUw!OLtkj8b8f7Y?&KJxJaOF!~uKXtPhYD{&u@o*1?twDn5)GmoZE3 zZ@F&c1{0t2#akr-!qq9o#z{T9m0j~3fq`~T^-?jC7xD-4>NZgoJ!9RfV9AgOQyQw5 z)BOQSd_3FxqfaS(u#hAy0tNcXXb03~T?j~Y>Gm)ydF8@tEqVN2E$ePqL9I&lSo>7) zocGba2$#AoeoJUEgWV+?m{egF>&*;`Mb>bCNQ(V8DGa|DkBaq$sMAI&D3~x^xi{S& z7Wb$Cn5!?Ue*+J+m+9Y=)1PGq)Ro3$Y^A_xQxYs#GrcLF*w7$_uc7mr#^>W$g zaH&4D&G$YQ+EMs800o^Ya~m6zUgOlrSw;gJa)rhmlY^zc@EY~Hc`KSZ^8tF~zVKO&)50&2|KqP2x&Sr?fqi? z=$b&LEulhyB#zSO7a^F{vuoOHM!)l{)ixpV!1S>Vu$~A4K})p}P7H$iq=xUsJS`|$ ztzxZN9&I13gM-W|A_%<>$p*an{8N52;m>T*s|H+S(IZ4p?-PxnQzblz`?kgaT_r8Rl-{fL-b z=V|iKniZ^`QCU&Zq@(s6TbU?1zzTsG70C}W-;azUT9xJWGN&E6h}+%eaMTR>(|hh| z&5mry<_<{ysXaM9;`aR7bmXBaXvYK6lL4EfRwW25BRVX&dk8tqB+WYk7UhAKNZY#% z%W|F94fBry+A>fHfW$2rIR?>#Fi zi8XKka8SfydhGtH1-bY%B}{HRc~NGIRNY`h_@HV%PwS7B*Tg zR*!1a1v!FOeT;!va>K`(6I3c5?2Ts=W~adXcUC&WF9OMBOKTml8nSEspLD(;JyHg>@dg0 z3Ok~)B<1wusV$tyL&Snl%s;|lqr1zsNg)p%!}E-CWF#k};@J!4Xh>fk2}0Q-AugT$ z*_M>B?(ZRVXH7D%F7-;B+{p1cx*6MsVJ{07pEJ#T@{7Hr~RF;4G)5_CY#{U zZn82-2^q=eJ%4|>X2q@szU*ZEzI zD?hlkRh(*PjG@kMI&4d((sL7+o%v|hYj>ibaY*B(eiH6Wp;GuId)nlN3rR{Kzo*Hl zR_I{X-F1X(BXmPGm#+X?ee9{w&ktw!McXgLB>IuY)RPDxoNSj6$|y?%`72iKBAer!Guf8##WqYV&1Ph zj3XQ!1i8UED&))g)nJLf6n(Vyi&!|^0BvIgq&focD;f;}=E1Vbn|KO@B^dmnwF=>N zItz^@nr*;p`j!aZ9ZK3Z!8h`XZNHGDqCMvqILvIXES~*M zq;~u=GtxzsM1X|(t~36#^P@Gk|Kq9lG!p*v`S$fAzkmVX`{!r$LCLZw{U?<@dFUU- z@mk4ES?qYLyvfv~^h68s4NvKFe%tFWGwG4rIvO*=o3sJ)MR$PTG|of54(A>WYvrl< z6>$p3xtFofL+8xDbf$5mPSb;eWN5DV_yBR<6B~lhpg7eI{KfcfECb(1f(ON;@u|9> z;2$Uuke{Ej6;JlngDnuJdva@76;@K3>ZG&^T!S`?_+ibDg{Dn(YK!M@)3x1})%Ye{H=d(J{=?IR()Di)7c$;w`65~3CVWkFSFD-e*|d4a>&(xnA)DhG|+xREkeq` zHGXC!u0$pSd-PgVT2kEBs&aK9truv>ba^f+`VN~^RecYFzcWf3nam|*rRnB&8<|K6 zZUIa*#O_;D4>Serou8Yx{EO{2Z=5`(0Vb+M1*w9yU9DHYE!*LN1u2< z&3*cm_7?8ACyZu)-BTFN$TKlg)4$>GqsaUpN@~Z#)0FNy>L{&mxTgrMRhL3eFBD4Y zxlAVe>*dLS51=K$Qop*&8YXE^MWsLgQXJaVneR!S6W2bBC3Mr>>j3`(PYvVNyJDBe zi`~870qilf6MwLi^Qq|>dSoZt>-@&{gqCnn&Ym-5v-iL4JN>V& zwCIXO)p$~V-eo_Ndfs&dEC^MiyjRoN)nv?&>B*`-e!V>=N19kttaMc(1&p86YmlLP z$$OJ~N!>%t!e3%-zAF`$_6Lt#!J})h@7pv)T7spDu}Pst1`?kf=4RW?B4Z@cq35_$ z9G!;0qrL{W`V0eoKJ$BrO-b{!F;1$;(#i zOE?PEp15qJI`Y!3{8owB;8vYIcSF1-7z;%TPT~}d-~#hsu_%Z)KC>Te?d3G{>wihC z`-hW_O_2PYha<6cRHJ#xB-!X_J>I)xfluH}-+Q(}`5d6jEkzmn4~K&uI-;R#H&Tl`L}U=4TSgRGAf*Kc%qbB4~w&>U&R>BvNzN zB8flhz5L^ulaPumI;jH;y=`kf&M2Av-HF~LGu%0Uvu%e!T;*Z)fVxF;{Qm9s_iYG8FvshiS<)`*nhe9c7Qw{T z0g$!k9_mF-y=%kmgwoHF#SL4HWb#a?c*PeN`uZrT1m=0eRfDlZv~&81*hy~jRJD&6 zkQ?xG{Wu8{p_OJiu8Fw+m%{Pcw;oj9N{vHC9e2QUan+_=&X6VcxtgT^Ol#T~Mc?M_ z41RWf&%gNe3(`BzBfS0lYRu~3Pj&3~m-KoY5*IV@$UC7x23u6fG%t7|R=8~rXJ(!C z1Pcg@sYtzpdGrz1xR5QIxcKgLR=l^sgy!#MG#nm&e`pYAXw9t}*^9GFLvgP2gUV)k zHOFEbC>~MiDZe~O*mJipY27J?Y{2Kw>+4sBrzSBH~N_JCAxdde4nfy1bR`@etZFj0SnRt8;{eyyxRbs-`_sh1hgE6 zXD?b( z#+6Xv?3RRpm~b69ySm?aVWgp`YA9kA^Q|AMdGyxXk&XznLhmHNI*Bf=tFlIm!Nb;Bx89Cc@4})DDWPFBkM$vWi+zG)e&{LEQi+gK|X^ou7!Y91b*XbZg#eyG(O!6MKD-97c9k*j#?SvQ?U6$vfmg~K;TZMO;1t%*i zKFO8!9O`iE_Ku<1L8`Mm$wXbn-_!1^edZIIgk$*8Uvw75FY{4YGXmf@;ffFbHPTWP zKg|@QwoCWQZ;i!F+1A{=sf#T~I<$)6gv+?)rHA`t&Hj4R`9^EWXG8>-YRt@P=L=FB z82rS1)}3@WEVn#^%cK2YEUW1)VdpS{q8j>-#pryg3qSejt`vr7i>xtRd;mG|DZD=s z#WBc@`kk)4B3GbY8TS|;-8*A&V!??DwfEWpj_1dLXl}cmFfOVs)A_Hcq{ZF@YPTz; zZ9ju)_ri701mQ_MzM?39>vdv;Q*PNYEYuhwc+rWP} zmf?}N?2^9G+N6=KW11uNr|aC*ehcKI2#L+3b*8;GF|~HQbI8)!`?H>gLT(Mejvz5B zA|K@?|JtoFfcKd?kMrh!xv^BIwG z*b03=)rImF@}b(nn!4*wrj4J~=g9E=s16a^k$OVu5|(5+2Eai51&YymZtG7&$j$B7oa{)9j? zZkC5Z?Yff?Txxp_d#ic>&m-!$hQ(Eqd$pa~u)Br4NfjerPa75)Z59;H6OGrN=YjY$ zD%_|BZT<>I>HbSRi31-GO~d@Ie&A1POAYek5h;@I?62r(&dCoF zq%*k3l%Dy=QGTwfCh2ct($Qe_m~K0*~DyG#`rC zJ&wARphQrn#PZ+@Jvj{kRPI>>p#@j*dtB{1npc+Fur~r%0CQ$b1|py9t!Lk z!l0`&C5K)ZG%(2x+`irOo8$$9NPJU^a7R1mdb@vN-fRD8U(_XW{VOx{W`NY%LnHjC zR_CSNLm`-(WU}`6l8A$q`0%Z%rMz5=ybjFNlEr|7z8F`+MzQ4l?A~i1u)#`ijQh9JrA`P_Z);$kt^%iLQogC|PcZeQsP8@%zh$zeJZODYq$Xa7lp+Xb6o z`q$1Mie`T`xP(I9c(B4VZkE8lG`#ic(5br57}opn`ErL5eNiS99^Dn#okr69xdWogcUUF zyIr}_j{DiyoQ}hTV|ymOsjjw&t#hpx(D6ZZYwlH*c>IGYdU7h>lksrQUSXPFH-f8x zzQ~L^0#T>2pE)$ulBHwZ0AtHSR0Ogwgg**g(Q$1LeUr<|q9eaU?e%=wi;OazNsU#X-ji%e8IFwAIOOcy#SbP)L8#DcQ8rZ|QNz^TJ-?^jK2W={Z->`+_%=udKxsRj@ zMY{{a{&FoaKAHVj0_KgMF{5DcUW%?)q{Ay3l%VQ4RC^0ewvvLa;VfaTxJ4u4p>vO! zj2PRLPS&4igW{*3%o1ppi|F`u&KUgQG2@9x+dIJrNq2I-&p!}w7rcyCMRK8CsuJdI zPh&!_VbB_k+fs$_P?*zRxWJWiccisnmIiNyl4ey`EFQePARW!vTQ`W*^^bXJ(Msd_ zwKO&QLUe+&X>|ks1rvT~ckfm@?KiO6El(nbv^yrcZuLAxs<1^)4v>3tyPs?JS7)-R ze7nz&&me2pfwAXxaZ2h~n0C*(Bc> z{9v;i7WisY1~^)IeTFgAoKxyjYHWev?P>2T@ZhaCQr3!64&(OGHG6y@>J2*AT8<9V zYYMMB1xzG0go$tQdmHM&YnvQQ+iaU6nOaxhYf;J{RH;83mAAtB1>CHfotQviNphSL zIh>0fzNWRMCt_(^bP+HUY6J(F>v(PvQAH8>Dzw)l~ZkE983^<;(DyA(S|Za?_bw|A9ve zLopcMaK!~4EiLd~6!(e;A;t4E5kgC4B5J-z#HM6FB-v8brB#4T@80eaXxtCMr*z@u zqjwQ}r~`SXi%dohhgM~b&Kka_h6jE0I5Po^Y@`t!ffE(#7HIaega2fqrN7=)cOp00 z(|+a+Cik0I?IM5}zYr!~A1XW9pPs;RV<0ZgloScLI`6+) zP}Wwkg04b)dX6TkQZsV0g5B>Q+LYTKytOll>=3*VKCtz&NTU3zWo%U8 zha(*oSG#Zq7pX6$fvhVXG05-9IhR_A0&}k?FwKSvfnYn#$6fFi(rMx#JSfK1G2_Sj zte^4U^Wk$i)zlWAtrF=X0(_L&@%Ok-r#|?@Q#UsjfFLaVSXL$)gMzeU)WNfKnJ0O}_QHF7GFRHOhfaOVSU%28#@>IYROD_d?QtC=7Uh|$+af2Q z0D598xDtzPjaW_V@DOM&$F0vjdvB8hrxk90#yay`_-HrB)(@PJSlX7U7<)I&?FdGt z5<O;9Z7j zunz^yaUCTe$`vu5s(W_R$Gwa9<@aZ{L2R}`5R502B;I(g%p{S*^QdxV(u=GLy}p|? zh>vM5tk-Ylp?>0jf#&Ys`91DB;@))qv{`(nQI20s& z{~aSJh4P%5ytd*3IwxOmZg1_0qW4X`wsrygoP*>AJF^`Wo;j;svSYwJF&UUga|{)i zpj0bk@C?#z*(JtV>EI-tYr>{{1a`QHHamPhJ3vuPaA%Q2R?|QH8t=PGn9$#8)b$d> z0de{QncLnHgqPz`_NyUVY)Q6P`4iW5-{+OjDR&7>k(NcG4MKY`JJh?lPw$5KPAezo z10DKxdxQn0MLIkx^2ek5+9j-PT)5AiiK?)Qg?fO?lDQ2t8CBNUjT62a~j3$%6+qRyX_e82zO%?kw>k-NA1(&CVNoNfG<%t2V`6yJU=K5mM8Ij&D4i=oaGgM zUJ#wAZLhq$?TUX10!1vwe}MDm;t}nzMg#WSNKtQiz41-aP$|3X{n<%Ag}674Ie{cl zY>?sZ`E${W?nI|0p4)i^Az}K>RcGfCQkVn8)R!j^(fEoReS^Yv;(AXXwO%ajGCea6 zC8V_q=+g39h=J2NL|{rS`}-dk0kRSt2{-Ga)p_!wR#%AVSbOy?9MC%N)*fYaIS2A-!4e1<4gh<0p+wP=a*F>;Z^48id?wXcZF4cygE+v z>ASkjw-D3vzf!2P$~zZ9awgd%sM}xkgKHdCN%%8OAHS&)F`BzHNZG9KyuM~i+wgIv ztM%V$Gh{=#3v~L6jBuB8kC4gatr{p7NIx6E5aC2Az7m|424sG^4$Gd|*Rxm6TYHkB zz`oInJ$l$dtW_)}b(+3=ahW{vD}oaNzoYYM5jnkNDITYm_S@TD{zSU4cL**m`sl?# zogiR0;dm?(3(lX!pmi36Zs#APNaB)5`OCa277Z$)E`eX?SUGLdWwRyZLB;#v)dCwA z1~^@AHMeOoi^F*48fKMAgacWw54bVYc{cbLNB=g|4Sja)GGlHk`rTSxRD`%~%;Jfe z?Jd?|?~b;s4Ss_>?hz@dE2mBXa9@QX`_b!~}qREj1NsnV1};Ua%#S6FD( zj7hiF^}j%@bB#&agfL@dGsz6^R0j8LpZX(9PekAo{bU&2BrNe|+M&{La8kdqBjcr6 F{s&fBk{AF0 literal 0 HcmV?d00001