From f9ddf04fdc5b4e3e34c3101683d36edd296fdb23 Mon Sep 17 00:00:00 2001 From: Thierry Zamofing Date: Wed, 10 Oct 2018 17:15:36 +0200 Subject: [PATCH] towards matlab --- matlab/SCRATCH.m | 123 +++++++++++++++++++++++ matlab/StateSpaceControlDesign.m | 166 +++++++++++++++++++++++++++++++ matlab/identifyFxFyStage.m | 92 +++++++++-------- matlab/observer.slx | Bin 0 -> 20958 bytes matlab/rscale.m | 43 ++++++++ matlab/simFxFyStage.m | 54 +++++----- matlab/stage_closed_loop.slx | Bin 37775 -> 37720 bytes python/MXTuning.py | 129 ++++++++++++++++++++++-- python/helicalscan.py | 32 +++--- python/shapepath.py | 1 - python/utilities.py | 91 +++++++++++++++++ 11 files changed, 640 insertions(+), 91 deletions(-) create mode 100644 matlab/SCRATCH.m create mode 100644 matlab/StateSpaceControlDesign.m create mode 100644 matlab/observer.slx create mode 100644 matlab/rscale.m create mode 100755 python/utilities.py diff --git a/matlab/SCRATCH.m b/matlab/SCRATCH.m new file mode 100644 index 0000000..6f11860 --- /dev/null +++ b/matlab/SCRATCH.m @@ -0,0 +1,123 @@ + +clear; +clear global; +[mot1,mot2]=identifyFxFyStage(); +[Kp,Kvfb,Ki,Kvff,Kaff,MaxInt,mot_num,mot_den,Ts,MaxDac,MaxPosErr,A,B,C,D]=simFxFyStage(mot1,1); + +[Nbar,A,B,C,Ao,Bo,Co,Do,L,K]=StateSpaceControlDesign(mot1,1); + + + + + + + + + + + + +function f=SCRATCH() +open('stage_closed_loop.slx') + + + [m1,m2]=identifyFxFyStage(); + controlSystemDesigner(1,m2.tf_py); % <<<<<<<<< This opens a transferfiûnction that can be edited + + %identification toolbox + systemIdentification + + + + + + %opt=tfestOptions('Display','off'); + %opt=tfestOptions('Display','on','initializeMethod','svf'); + %opt=tfestOptions('Display','on','initializeMethod','iv','WeightingFilter',[]); + %opt=tfestOptions('Display','on','initializeMethod','iv','WeightingFilter',[1,5;20,570]); + %tf1 = tfest(mot1frq, 6, 4, opt); + % Model refinement + % Options = tf1.Report.OptionsUsed; + % Options.WeightingFilter = 'prediction'; + % tf1_1 = pem(mot1frq, tf1, Options) + + + + + + bodeplot(mot1frq,tf1) + mag,phase=bode(tf1,frq) + figure(1) + subplot(211) + + + bodeplot(tf1) + + Opt = n4sidOptions('N4Horizon',[15 15 15]); + n4s3 = n4sid(mot1frq, 3, Opt) + + + + + + %tf([1 2],[1 0 10]) + %specifies the transfer function (s+2)/(s^2+10) while + sys=tf([1],[1,0,0]) + bode(sys) + step(sys) + + sys=tf([1],[1,-1,2]) %instable + + sys=tf([1],[1,1,2]) %stable + + + %0dB at 12 Hz=12*2*pi rad/s =75.4=k^2 -> k=8.6833 + + sys=tf([10],[1,0,0]) + %1/s^2 -> 0dB at 1Hz -40dB/decade + %10=+20dB + + + sys=tf([1],[1,0,2]) %not damped constant sine after step + + + sys=zpk([],[1,0,0],100) %stable + + sys=zpk([],[-10,-10],100) + + + %parker stage 1 + %!encoder_sim(enc=1,tbl=9,mot=9,posSf=13000./2048) + %!encoder_inc(enc=1,tbl=1,mot=1,posSf=13000./650000) + %!motor_servo(mot=1,ctrl='ServoCtrl',Kp=25,Kvfb=400,Ki=0.02,Kvff=350,Kaff=5000,MaxInt=1000) + %!motor(mot=1,dirCur=0,contCur=800,peakCur=2400,timeAtPeak=1,IiGain=5,IpfGain=8,IpbGain=8,JogSpeed=10.,numPhase=3,invDir=True,servo=None,PhasePosSf=1./81250,PhaseFindingDac=100,PhaseFindingTime=50,SlipGain=0,AdvGain=0,PwmSf=10000,FatalFeLimit=200,WarnFeLimit=100,InPosBand=2,homing='enc-index') + Ts=2E-4 % discrete sample time (servo period) + Kp=25,Kvfb=400,Ki=0.02,Kvff=350,Kaff=5000,MaxInt=1000 + Kp=25,Kvfb=0,Ki=0,Kvff=0,Kaff=0,MaxInt=0 + + + num=7.32 + den=[5.995e-04 4.897e-02 1.] + open('stage_closed_loop.slx') + %sim('stage_closed_loop.slx') + + sys=tf(num,den) + bode(sys) + G = tf(1.5,[1 14 40.02]); + controlSystemDesigner('bode',sys); + controlSystemDesigner + linearSystemAnalyzer + load ltiexamples + linearSystemAnalyzer(sys_dc) + + controlSystemDesigner('bode',sys); + controlSystemDesigner(1,sys); % <<<<<<<<< This opens a transferfiûnction that can be edited + + num=[8.32795069e-11, 1.04317228e-08, 6.68431323e-05, 3.31861324e-03, 7.32824533e+00]; + den=[5.26156641e-18, 1.12897840e-14, 7.67853031e-12, 1.03201301e-08, 2.05154780e-06, 1.34279894e-03, 7.19229912e-02, 1.00000000e+00]; + mot2=tf(num,den); + controlSystemDesigner('bode',mot2); + + + +end \ No newline at end of file diff --git a/matlab/StateSpaceControlDesign.m b/matlab/StateSpaceControlDesign.m new file mode 100644 index 0000000..a4cc4ab --- /dev/null +++ b/matlab/StateSpaceControlDesign.m @@ -0,0 +1,166 @@ +%http://ctms.engin.umich.edu/CTMS/index.php?example=Introduction§ion=ControlStateSpace +%zustandsregler: +% web(fullfile(docroot, 'simulink/examples.html')) +% web(fullfile(docroot, 'simulink/examples/inverted-pendulum-with-animation.html')) +% web(fullfile(docroot, 'simulink/examples/double-spring-mass-system.html')) +function [Nbar,A,B,C,D,Ao,Bo,Co,Do,L,K]=StateSpaceControlDesign(mot,motid) + sys=mot.ss; + sys=ss(sys.A,sys.B,sys.C(3,:),0); % $$$ only output position + + %sys=ss(tf(mot1.mdl.num1,mot1.mdl.den1)); + %A=sys.A; + %B=sys.B; + %C=sys.C; + %[A,B,C,D]=tf2ss(mot1.mdl.num1,mot1.mdl.den1) + %tf2ss(mot1.mdl.num1,mot1.mdl.den1) + + figure();h=bodeplot(sys); + setoptions(h,'IOGrouping','all') + A=sys.A; + B=sys.B; + C=sys.C; + D=sys.D; + + P=ctrb(A,B); + if rank(A)==rank(P) + disp('sys controlable') + else + disp('sys not controlable') + end + + Q=obsv(A,C); + if rank(A)==rank(Q) + disp('sys observable') + else + disp('sys not observable') + end + + + t = 0:1E-4:.5; + u = ones(size(t)); %1000um + x0 = zeros(1,length(sys.A)); + + [y,t,x] = lsim(sys,u,t,x0); + figure();plot(t,y) + + + poles = eig(A); + + if motid==1 + p1=-3300+2800i; + p2=-1500+500i; + p3=-1200+10i; + P=[p1 p1' p2 p2' p3 p3']; + else + end + K = place(A,B,P); + %K = acker(A,B,P); + %K = acker(A,B,[p1 p1' p2 p2' p3 p3']); + %K = place(A,B,[p1 p1']); + %Nbar = rscale(sys,K); + %Nbar=1; + Nbar=-1./(C*(A-B*K)^-1*B); %from my notes) + %Nbar(2)=1; %the voltage stuff is crap for now + if length(Nbar)>1 + Nbar=Nbar(3); % only the position scaling needed + end + sys_cl = ss(A-B*K,B,C,0); + [y,t,x]=lsim(sys_cl,Nbar*u,t,x0); + figure();plot(t,y) + + %observer poles-> 5 times farther left than system poles + if motid==1 + op1=(p1*5); + op2=(p2*5); + op3=(p3*5); + OP=[op1 op1' op2 op2' op3 op3']; + else + end + L=place(A',C',OP)'; + + At = [ A-B*K B*K + zeros(size(A)) A-L*C ]; + Bt = [ B*Nbar + zeros(size(B)) ]; + Ct = [ C zeros(size(C)) ]; + sys = ss(At,Bt,Ct,0); + lsim(sys,ones(size(t)),t,[x0 x0]); + + %https://www.youtube.com/watch?v=Lax3etc837U + Ao=A-L*C; + Bo=[B L]; + Co=K; + Do=zeros(size(Co,1),size(Bo,2)); + mdlName='observer'; + open(mdlName) +end + +function SCRATCH() + A = [ 0 1 0 + 980 0 -2.8 + 0 0 -100 ]; + + B = [ 0 + 0 + 100 ]; + + C = [ 1 0 0 ]; + + poles = eig(A) + + t = 0:0.01:2; + u = zeros(size(t)); + x0 = [0.01 0 0]; + + sys = ss(A,B,C,0); + + [y,t,x] = lsim(sys,u,t,x0); + plot(t,y) + title('Open-Loop Response to Non-Zero Initial Condition') + xlabel('Time (sec)') + ylabel('Ball Position (m)') + + p1 = -10 + 10i; + p2 = -10 - 10i; + p3 = -50; + + K = place(A,B,[p1 p2 p3]); + sys_cl = ss(A-B*K,B,C,0); + + lsim(sys_cl,u,t,x0); + xlabel('Time (sec)') + ylabel('Ball Position (m)') + + p1 = -20 + 20i; + p2 = -20 - 20i; + p3 = -100; + + K = place(A,B,[p1 p2 p3]); + sys_cl = ss(A-B*K,B,C,0); + + lsim(sys_cl,u,t,x0); + xlabel('Time (sec)') + ylabel('Ball Position (m)') + + + t = 0:0.01:2; + u = 0.001*ones(size(t)); + + sys_cl = ss(A-B*K,B,C,0); + + lsim(sys_cl,u,t); + xlabel('Time (sec)') + ylabel('Ball Position (m)') + axis([0 2 -4E-6 0]) + + Nbar = rscale(sys,K) + + lsim(sys_cl,Nbar*u,t) + title('Linear Simulation Results (with Nbar)') + xlabel('Time (sec)') + ylabel('Ball Position (m)') + axis([0 2 0 1.2*10^-3]) +end + + + diff --git a/matlab/identifyFxFyStage.m b/matlab/identifyFxFyStage.m index 55dc568..bfa22ce 100644 --- a/matlab/identifyFxFyStage.m +++ b/matlab/identifyFxFyStage.m @@ -73,18 +73,39 @@ function [mot1,mot2]=identifyFxFyStage() figure(); mot.tf2_0 = tfest(mot.meas, 2, 0, opt);disp(str2ndOrd(mot.tf2_0)); mot.tf_mdl=idtf(mot.mdl.num,mot.mdl.den); + %ss([g1 mot.tf_mdl],'minimal') this doesn't work as expected - g11=tf(mot.mdl.numc,mot.mdl.denc); % iqCmd->iqMeas - g12=tf([1 0],mot.mdl.denc*12); % iqCmd->iqVolts : iqVolts= i_meas*R+i_meas'*L - num=conv(conv(mot.mdl.num1,mot.mdl.num2),mot.mdl.numc); - den=conv(conv(mot.mdl.den1,mot.mdl.den2),mot.mdl.denc); - g13=tf(num,den); %iqCmd->ActPos - %sys=ss([g11;g12]) - sys=ss([g11;g12;g13]) + numc=myNorm(mot.mdl.numc); + denc=myNorm(mot.mdl.denc); + num1=myNorm(mot.mdl.num1); + den1=myNorm(mot.mdl.den1); + num2=myNorm(mot.mdl.num2); + den2=myNorm(mot.mdl.den2); + g1=tf(numc,denc); % iqCmd->iqMeas + s1=ss(g1); + s1.C=[s1.C; 1/s1.B(1) 0]; % add output iqVolts: iqVolts= i_meas*R+i_meas'*L + %tf(s1) % display all transfer functions + num=conv(num1,num2);%num=1; + den=conv(den1,den2);%den=[1 0 0]; + g2=tf(num,den); %iqMeas->ActPos + s2=ss(g2); + s3=append(s1,s2); + s3.A(3,2)=s3.C(1,2)*s3.B(3,2); + mot.ss=ss(s3.A,s3.B(:,1),s3.C,0); % single input, remove input iqMeas + mot.ss.InputName{1}='iqCmd'; + mot.ss.OutputName{1}='iqMeas'; + mot.ss.OutputName{2}='iqVolts'; + mot.ss.OutputName{3}='actPos'; + % u +-----------+ y + %iqCmd------->|1 1|-------> iqMeas + % | 2|-------> iqVolts + % | 3|-------> actPos + % +-----------+ %h=bodeplot(mot.meas,'r',mot.tf4_2,'b',mot.tf6_4,'g'); - h=bodeplot(mot.meas,'r',mot.tf2_0,'b',mot.tf_mdl,'g',mot.w); + %h=bodeplot(mot.meas,'r',mot.tf2_0,'b',mot.tf_mdl,'g',mot.w); + tmp=tf(mot.ss);h=bodeplot(mot.meas,'r',tmp(3,1),'g',mot.w); setoptions(h,'FreqUnits','Hz','Grid','on'); end @@ -109,10 +130,8 @@ function [mot1,mot2]=identifyFxFyStage() %create ss from tf MIMO: %https://ch.mathworks.com/matlabcentral/answers/37152-how-to-convert-tf2ss-for-mimo-system - %Gspm=[tf_iqCmd_actVolts;tf_iqCmd_iqMeas;tf_iqCmd_actPos]; - %sys=ss(Gspm); - - %normalize: [1E6 1E3 1].*mot.mdl.denc + %http://ch.mathworks.com/help/control/ug/conversion-between-model-types.html#f3-1039600 + %https://ch.mathworks.com/help/control/ref/append.html numc=myNorm(mot.mdl.numc); denc=myNorm(mot.mdl.denc); num1=myNorm(mot.mdl.num1); @@ -127,45 +146,36 @@ function [mot1,mot2]=identifyFxFyStage() den5=myNorm(mot.mdl.den5); num=myNorm(mot.mdl.num); den=myNorm(mot.mdl.den); - - %http://ch.mathworks.com/help/control/ug/conversion-between-model-types.html#f3-1039600 - %tf2ss MIMO - %https://ch.mathworks.com/help/control/ref/append.html g1=tf(numc,denc); % iqCmd->iqMeas s1=ss(g1); s1.C=[s1.C; 1/s1.B(1) 0]; % add output iqVolts: iqVolts= i_meas*R+i_meas'*L - tf(s1) % display all transfer functions - - num=conv(conv(conv(conv(num1,num2),num3),num4),num5); - den=conv(conv(conv(conv(den1,den2),den3),den4),den5); + %tf(s1) % display all transfer functions + num=conv(conv(conv(conv(num1,num2),num3),num4),num5);%num=1; + den=conv(conv(conv(conv(den1,den2),den3),den4),den5);%den=[1 0 0]; g2=tf(num,den); %iqMeas->ActPos s2=ss(g2); - s3=append(s1,s2); - tf(s3) + %t_=tf(s3); + %bode(g2);figure;bode(t_(3,2)); %connect iqMeas from s1 to iqMeas of s2 - s3.A(3,1)=1 %WHAT NUMBER ??? s3.B(3,2), s3.C2,:)? - %remove the direct iqMeas input - %s3.B(3,2)=0 - - t_=tf(s3) - t_(3,1) - figure; - bode(t_(3,1)) - %compare with tf iqCmd->ActPos - num=conv(conv(conv(conv(conv(mot.mdl.num1,mot.mdl.num2),mot.mdl.num3),mot.mdl.num4),mot.mdl.num5),mot.mdl.numc); - den=conv(conv(conv(conv(conv(mot.mdl.den1,mot.mdl.den2),mot.mdl.den3),mot.mdl.den4),mot.mdl.den5),mot.mdl.denc); - g4=tf(num,den); %iqCmd->ActPos - figure; - bode(g4) - - %sys=ss([g11;g12]) - sys=ss([g11;g12;g13]) - sys=ss(g13) + s3.A(3,2)=s3.C(1,2)*s3.B(3,2); + s3.A(3,2)=s3.C(1,2)*s3.B(3,2); + mot.ss=ss(s3.A,s3.B(:,1),s3.C,0); % single input, remove input iqMeas + mot.ss.InputName{1}='iqCmd'; + mot.ss.OutputName{1}='iqMeas'; + mot.ss.OutputName{2}='iqVolts'; + mot.ss.OutputName{3}='actPos' ; + % u +-----------+ y + %iqCmd------->|1 1|-------> iqMeas + % | 2|-------> iqVolts + % | 3|-------> actPos + % +-----------+ + %h=bodeplot(mot.meas,'r',mot.tf4_2,'b',mot.tf6_4,'g',mot.tf13_9,'m',mot.tf_py,'b'); - h=bodeplot(mot.meas,'r',mot.tf2_0,'b',mot.tf_mdl,'g',mot.w); + %h=bodeplot(mot.meas,'r',mot.tf2_0,'b',mot.tf_mdl,'g',mot.w); + tmp=tf(mot.ss);h=bodeplot(mot.meas,'r',tmp(3,1),'g',mot.w); setoptions(h,'FreqUnits','Hz','Grid','on'); %controlSystemDesigner('bode',1,mot.tf_py); % <<<<<<<<< This opens a transferfiûnction that can be edited diff --git a/matlab/observer.slx b/matlab/observer.slx new file mode 100644 index 0000000000000000000000000000000000000000..922620de259ad6a8ea75abd26025623de1238f5c GIT binary patch literal 20958 zcmaI7Q;;rPuq^nsZQHi3-L`wTZQHhO_io#^ZQHipJ^z`fbMHjV!+Ke78Cg+LnOUVE z4FZY^001BX-y)n+gdRa$IzRxx05AZ6^WU$HiL-&RfwKXFk)5N7f}@?iiKDZHi4(oM zjkS-Gj(iX!>Nic%tpazp95Ah7BSji!&itEUrsHUM6PXR(pq|}iS{6OIPyT6ku2=xkO>uV5~_%287yu+mbX*bcL_Ivz{;p3pdUGL?|m zio{v^GSvo_@+hWN3y7SHKH5N{RyilzO9iyH@{15mUn67J5H9R=dBGai$MaD*_+++h z>T_uYcTl!LnST%?c9Q=d0Ud*^*5En@iANE8ej7i{O}eDV+;Nd*Cb{`PIkQKUuZ|wS z3*Tji5Gyp0rPnZQ$|6e^dZ#jEZ&)(n^K_~@^_z8I2Lf+n-G7^%)6gc3@!Ic1_mXWy zinwC;UmtY+y*oUEKUx)6QI74+Al@8u9R2jlqRkqc=eE*gb{ikZA63r~P5k)c2`CB)^T9v?R6y3GX8 z-NWCPr<|k?qm(h$i|@Kc=y9iTXPoEJ=P&0aDQpG*q2NMYg(hyUEP`B`9#wjEEfK|S zVdjhS<_GKOHz<~|_ySpknnu{%MyfrtwgRrJ7l_8==EmG|qUKAbIUy=4PjkMMc4c)B z(vwnU$*Y<^ydKC5cHIW9Qr6c+YnTn;Yxjq`k*LD4{tW&p8tH#_gAPd5>3+ z{}cael1f+Ff4CR^d*uIr_(k2FO>CVk>}>z1Ky_m4;R6{_Mg59BI&Fwx%_(KLdCEQ; zy8r<#MYTv@UQBYFJK0-F%50TcLs7&C{WxPq;6gj0u6H>BpVy9*LiT@l4^Vlr+);we zFDB^NMLI?pm1RKc*LZW`o~WiNiq~F_;)@5D7>HtNZR{<@cnrm{)#kc>8g#x%Aoj^i2Eoy?K4#x;=FgsK4VO2Noxd^A{up zTEzJ~4fOf`x(5bwil8=R-FNrCQ|Ai{>`?oKuT7<=*ay~Ia;|#V*{?BK7`~a2)xe4l z%j8|LxK{MBaL5egtwHMgKtMNoX_Md^AGKI5haFLWiNCr{Z&FIhC|72q`&xdPW9QPL zvsyW&<8an?+zJ&IB9<`keRPJDn0T!S5Ak7TTgqyEu;Sw^g0YCk!rMBMJ+^Z?W zBH)8za_$tb&mVUcyMm_aBc#98oc#@W?Sp#`4-sMj8Z`2ZL8Cn|^MiW1_ zXtA1vK5y$J_yqx$l@64(T`f2GPr+^$)C0SZj$;9;27K7#*k&i&Ls?s$ev?s(>};PP zYKrC8d4Nm~lJ;>bdXcW^<3xR*77O$ZYFd~8?lrp!CjApT5M=t+TK9s55WeTT`OIx(sX{c8O zEsu!nRR`w@d3zE z|6aRD9k@Vc}oj(oHIN)@ZjRoj&a}TBb3~H z*xFLIu(0@hneoGtCt3vehqMfSpj19jZMw}zO{1Mm>(|1~&Q1!^)yOI7B0sZNh4W46 z+}F*|Zz#pQUy?fRiXeY?C&#j(3FPxC2n9iTV^g-Vi6E%J#cC@Mu{*u`zAotMuO#<4 z49rbY9KqpASZa_2=>o|bp}*J?A&X=DdPhr1q*&qp%S%I*FyDVoC`61hnYArGiZD>FmnAoh5D}NZtyy_=nL&cKY?G zcDF~T*T=pNcXs1YW?A+_Ga*gA|GPN(`+&zi$IXtJMQ@>`I*79y8xr&&H=FUxX3cl{ z22`$W_s%~DqUhH6&~>7Py>F?7NaS-r91+t^?n+|5+Bt=IqM!fIz0NO$$aA@d8v?@P z=Mcn-iFBU8XqU(1$k$18^&ZZ_ga9ZN(<2 zQ}HCBqlc}nrAb4@3oRB!o5$?q?n^A>FI6VD(hezEQOZU3L1;7^oUIcpQsX8-`>Rd z8}>%`E_O%vv(tYgs@|C17u=`P>KGd`#@4@&@G~?$PJWX%a%r22E&A_?Opw#HmmXGcz;(jR|=cYmTr# z3CsO36Jukjsi_;46ZjMRs(=3ySoZPSSnVk{o-gc9KsTsGA6?l9%^Og-*8<5q5Arbm ztC@KBY6cbl)grY3(e~v+pL~XxS+i<-k(Wrsd8XQ@waOAGmeq`>py$tR-Z5G0!d3BK zdEk9n;WV*$u_~!B{)d%IaY*RB*4Xs39IO%7UR;gafe)O2uJxD?`9PLtE)9NIwy%^q z=KZ8G4Aaq(VX~|QEjui(^8x7x;lvI>_k@D#pP1jF-PJO_diQOeC&c8Du!&k+U2TKT zj%bk^oHCp5wO)^CJxtMWc+Y_yWB9Wd_vVM~f$LR>@S=xDLs8YghLH3qlUO3=O&NjL zY{RN4FOPse7s!s>*{Su`9({BF^x(~A51p3Zq}CQsq?Ix@HE>vK!1zJXeKrct5rO^n zaP9x|oyy(U%;;m_g2pcI&DDVWynVmGG=hlvFSCgqdw$0$Olfv*^o zJnqY!g)~^z7!*kU@~Ubzz_H@89?0t&9A1vmO5WDjm_kfHhn9zE;n~1LbS~`+9piZ``r!5AA^pWv`{R z{r(gyF@2Aym_#a3nJkZoAu)tt94t~1+@R5fX?nKPKwjnr?c3-wsKy>WeUOlx=*S3JG%sUiQ-pDDZ(#-#UKe7}6TyN~D|Mr5Gut|8)Kt{e`0ig{<>h3D>({TshZdPJ z5*gZobtdK%rViWfxSE?2|9D!7oe&_|yE|?)z;O6$+<(4mM^m7<9_2q$%ZyO;E)J5A z_*Rv#tkzP>xAHFM9Lu%6eN+%<=y7vrzH)h~-po&{>$5pYOTGO{9aoCQI)2W0s5y*F z!f=x$1R5`0JR(IzX@|*yUkbg21lJ_PpQ!$g(wjw@f5M6~cWb{obu-B8hWGNx%KtG2 z#G!%5w+>qGX!#l!lHs|O=EwDEYkyJD*<$qU&l?X}u_m!}vS@z6TlmpnKJ;R{EX=L` zD(oKW9g|Y1Vc|@U&d%n=rS4gw=^-jd*;T8f9%#$-6glhzx6}gE{8KrBl)s$1|CCZ6 z>r@rHg1j9usqQNa`u;Kt!-3p;Y7R^9`IS`HHsC;G=@113f{=Nqx-y!Qe=Bb>uXN=dSTsVdw+^H~o}h+~5WjzX)A3&vIg zyxU&PFg`U4ydpwzEE^EIni}&_%WB>vva^tn#n(~*oy&%MOvWVu8 z?9cG2;3azWB`+|Ngk#7aB#2wjO8omXd;h*BBk!tRylkOY;$3&$d z<=WsPgfc7x5HV$4b?%W;h73C)%}cccELr3ehnrq}+TR??KH?T84@FECEz}D{nGGgY zE|O2Mt0%!2%oq-mfi>jm`TV^?8kmK)K)rE?MH1?5AxG7V!`2MZAq8o zqW8m!&sazEDtx~Q0fC?HoAZ{u$dw#Km#Gv&;FbYkxxdM@SPMcrgNe@PWqaJ{UrOBc z>1noH)8n@zHQT?8={UYI)$m5Uxidrw#+xuu6xMOv?KJfD_s|%)Zq`@=SEM1dU{fwK zU6IFs+WA-Li|dtU!_h!|khQl^Ot)-k76$H=Jhis~awtgc3@K3AG>VCiTo8#eT$^c= zFm5bVH63bh`HhHybmlkc*M*XhDCdwibYL_u_O06TPB`cn$rDiD?Shr|YoYx{&k?L* z8eQ}NmIyrAUX=Q%I^8l+eDIX_VvL$NjD#2ZQHv6lq0~L;35C@_mHA6)uWj6JHwt4% z3dsVsy$g+g6*On~%oeK(*UjZV59>^Ea`}T@Jn$Qf#}WH9Cm0&5Pru41-A?{1?4CpZ0z5MWh_8T)Vd{ViFwOoCo(i}Ou?O+8dh zHmz5v3TIkD2v^Ui4X59(d?}PMnC~2M_7Y=wprtGmS=3UWz1rOJI=TaukLX{9L9(6Z z#EuG`m4f&X5?=Mf?1^=hMet!=*mpire<<8%Vd%HRU{Ce8I=w=7aO$8|Hfd?R-ed_V~saw7~m8W+`nQUSh4~?SsUFOV#gt4;1 z)AV|QUQ%^|=!JrWaF1#Ueoel$Whx(je}^ak%BKlpjL|8fc-|r4H(s)6VIy;(IGV2NS%$wX?Us<U=4_U7YSjURQw-c2Cq2by z`vKrHh25Stl3?JJ@!w;T5iSGR?`Y6u5aEPyE-p+@qA0;+ijDLX6#7%994DNa(O&Js zO59UL@#Tk@Hq{l`S&SF`_At~Mz4I*0c#yo!F2-E}aHI&RN%8`?PKYxi=Mcd>CGc^` zr^D6oz^geYJhe^wwJ>k#mJ?z=!0oJ_NQ3BP2q85{pF0J-!EJ%}-*isKtH?Y!jg5Ev zk*@Qz)H%tV=~Teo zk{OLsZOWS4X!m0}4y~|rc~wjn7Y>2Pf;twO))`w@Xp)@>J@Zi3%NzABXIlXC6gM zBPo&^H?ErK`%8!A@a2Zk`rK@t*JTKPLu(^1->Bls49Mv{HDvb|ad9ug8L3jiwgj2*w2#(Sdc}~9Ld8fP(m`G;DYCmy4s*)Cr z=VOE7Nva6*w6`A+2mD7A0t~~xn=yw)KG2sb1W4NgJ&G9}BWq8^{K{2b%yUawpFB1j z@*y~D)PT2^x2AokuA#hc^q;0u*N5|i*O|D#lQnOixCgU5pHGF~byr`ed0{Wt4UKiP z0vma0c)LAs<_Kw}kd4UhFDH*A3%A}4wbwIAKKxDN4UbViKjwZdJX8F=^*Wm?WX89c z0ytB;AbvYCJQy$M3q9X!T#a5nopx@{pa1Oqejp!#-SSc;{+Jze4Nt`Af0w9Q*@VFByFU|;YN>v%Fc-Q|pH zyrOm3X6&`i;{clRr4EI#$>!fcSe&^X|J=Rf8~gRp(gnGX{@Q_8@R(g>>rIreuM%>4 z+|sL6s)19<_So?jpu?gW^txUv*Vo`%fzcIzyh3yToy4YRkQ%80L%5q;MGBB3%(B!N#Qy-lPf!RM8fFTMM5H3jg{s3yZvDV*R0##_0 zNH zhOs1Ywy!{7%3>u*s9h4zXmo0Tfb72<&)5y&IQec6fi#fuzIB6!VjLdUjpH%76CNNq41xZ4C>$9PH;~p9 zlBfejsTqJ8xphZqU#A)71*5${d9cR=VPI2fMzFvM`30}2$wOZ3pbi}WL{Y@QVjypG zWX9x#Ova+h&*H?6F;IC6 z5|*)hmxW_)$(%;e5x}%I68PuFL65`5D8S$!WSJPu)b>cFI8&#=?I|y~J>L8C4RdB( zH!GCG^BIGlW;~a>Hp+si5>@uwNV88=nZrK+)dWe;l{mSI%JvG^fOk2i98%;)wa@^9 z(k|hMBwF?EQGJWN$d|rlSX(uzNA(`lmy)BP;ZYVHM6QQ*j%SyD+I-or7oG5B=<+wd zZ3_&7syFcB+A z2g8P`AQO~Ia3btdIl~2ny@Xo&im9L%QyK*T5&8K*uob2Y2_(SN1YtpRuQM$vf;!r? zOj(id3QRg_UMkCr$i5yzm(^`@{3Q0@#5Ez+96eSDE8zwNTznt>O%!|HE!lMk$=XEM zGG`mvMG|7g%bXb0 z1hOuA@2u~&T*49XNoqgb`c$QAJYFA~>XFtfjlLXkNlDzF2un6e_ z9yT47Nuk+Ev?y`rE}NVNrc*_G`tB^>Nl3z}xgvAG%XW5^{ts7lkbs|F7fMA*(;%PI zfD-E(!#V2IuLjJ6%O?1@<(s+w2y>1C<{8U*vi4N{%Ta z2PVAV`Y(0;&evRr-jG`^UC;zp@ji?9f0;2e|O6Sn;c7@bKxOQq$+2HGF=G zZ|FtlD4?OLYXpQ~jJPOjqU%thavh7;_2@u`?;0eEvuJ$; z|EfS;=Bk+z1z)mPVb^3|TGHy&A4c8c0gSjmFKW%|FStR^wd9amL>JMXma;I`>s7tW zmUIg8&AOX~xtYc4G*`?EppVj%En92EVk@X0ve@ZqFyorKvtDGRwn(Zg_&^21(BrHV z!Os4{*mN*(Z6O#Ucje^+8qyLDM56ZTOK8mh#(d{a+I}ce{}#lGG~y8qe{>|ByY2p2 zg4=a7`AW=rAk|ADv199qJW=*d#t9O9)S3TG*6y3mvZ=<1eB6O>I|Om#bR%Ip3Sjo z7waS;#fycs@{fcO#KWRLuReMkSb8b>DPE$*O|x&&VX66n&D^#J1fom^%?-OdKdE%P z5apQTt{>NFZYNEP{BUjU7a=fsjVY-Noh9&vMA$VN+%S9brj~M_aITEQB@1@wS$0Krdil>XqG+|oJD8;KYwHZBX?r1jxf`@FDiu`sMC-5ud zKctGz4W!UT0UJ`lAm{1^BrxRUpq(6znIZXr`9M@Z$Hsk=df*)_DoQSaJAWh*;g|Q1 z3_e%)SF+*3xuJLa+s^HGbNd(CZ8hcQ=knXir;C-U@37lj%0Cn);y8(8`h+9{i#|j6 z+8h-u3b?c9;AtG@0P)tykb*)STknUMKyPd*9@&~L+;^T{T(Xcl%-}+dMKVSZ7^1(p zL4M+4NwK11ey*w{giS`4XZQc#W1nl$$)K1@MQA;(oq*cFVljc(meRdTl>QqcKlv)$OKqWjKs@@ec zjpl8s{ux=dayh90L0yQ|p9_s(c#?`y`TKsa1Jwt~02eJy{yC`y6%O)mHb>ONg?KhV z+BQY{eApv#R3j6@9pKt`ufWchwg*Ifjq>6m{$N)+Vfhe5J2KFtv~1n8d{JVqo$f^| zc1)Ynf`^LWLVX=;Sf*v+h9=xC)me&p!Lpb8mchfZiT{Dd`ua9V{n0TIl$<15-pJ9L zcoD5x!Clxc?~+O!UvRTjp=>2r-S%6!ldG4{I>i&o-a@>(9g)o@8oUoL(Qz+}C@w(d zExKin;#sq@Yq|NU9ENf-1^qXDX4|~nm@u*J;QL=VP`!BB?^F?nSY7G`2mYq+#65dj zSe*^kdMXOrtGu&g_OeDhFxn_T6^D!HuL^aYM>?=WrB4VjQT2x>Hf!bc>qciq!EK#V zmDD3mS%XQcoH4eZbjo?~vS&t1g~mwWL#S~vTvt1W;N%n`58C?g?m9LWH-QV;9$otSw#BE-6CB zZ;;^UwBP^oM*nB$#uoF>a`}I29p_2^mp5u;XKQLqiLq-T5v-Qd?swu@UcgvDf>{zs&-RF_KhU-*fc$kLl5y-q)W8i@rs+O;Fa` z@8ZRGR;}-xP><99{PPDALwLhP6F1U=UT2V0R3>SH9kp~F5wIflQ- z@xJySlv@4HdOn!&e-^qgK4QrBnCdJU%OpJDM%ju+UWMdDdV+J^yfEw-HY`B7cB$fD zuL0kxEk_Qm;*r@ngnL~b125?8=kiK7|eWoIw^O(0*` z^sRxWTKMAZ+=NP4a8cRU=e0qwS~_#;u?G(aVUv|{0&_a1S3)RQiuAD7_c^Qt{^397>8@abyBI47jY^396C;bQB*iNGIkE#aHF5ZxQ7 zp5ie8bfO8xUBtD*{GAZl+-LI_r!ih=uSMqH7 zeF&JT*r7sb0fxoIWcQrgXpAwv$F)2m1gQr@~UI0ZFMffudq_?4m=+T8nh&c#eDrZ`1>yDkLCg$O6ruVjNqj zHj>rs(c~fxflkx+%y~}W9ins)JNRA(^MrGa+bFnrcY0H_eRm{5w}vKaDu^;?O1uT!9#TCd$GExVkrV z;SO>{o9PID0fX3@^VF;U0Px-3_^m^^356>mtgQgU@P0|rO}!)fQWSeTLpjIr z2Pu)P14X>A#6XU`6lI4!(IFTuia7FF^z@>9WPcs6rpR(NXZA4kW)RZ{P^0I&*eC}j+A(k2>hE6zMDx+og2m-zv!F$RAs!>ynkB;1?CcOduL(B*$GF&jusaI*Mk^u7k zoACqB$y+O|f;`yuG_~)oYao;WkC1#kMY%W*-|FQ_8f8>vW@ibvcQr9o$nSY$1L1V< z)L{X$-HuICUfr^MmskP zf&DR7%;GbiWcQv^08k0#f`X(H}W@Tf%#xDiIK^ubwv7U!nZR z$cad4+QJ+2wr6bl%3PxADH7Uq9a`Y`MnxMO1 z^U=a(sK!Po6SWW5(Qf~eNxPLVb+A~iH)7FB72{N44WeL12(t#uenuE}?VT}v_;8<| zfC$#!*tOy=Jw7A0_MH$XUNV-nD+Q+<0at<&_8)gRAsE9W`SksTfBGXI&|={*WciRa zysD+%H~lxF90uHk4zT9s(wMB2BG^WEFT@qqv07APIS5bvrL@u#__jIPkX+2f;-A`? z%$2B1jFv027+~6P2z2^S&DVwEHA3M~;7SqA(rdvg!MMeM!*CQ8X)Qga7cU1$LOeR5 zw=&@DktNRbx;^a=ihLi8r`JG(^osGe>va+|1hso* zg{!)>7A4S}?Z1)@_fy9PtbN4}!7lMeDqdlo%1ppzS6Qr>8U4xC_gQe$c$zrT&Gje^LX1F$#TMW59G%sdy+Mo6jF);rf}-Sq~ar;JA`JaIq$M(||XE zW?7e;W!;y9ml+pn^<&lGc|J^P(jaQ(652!S~scY>eO$ ztAjPE7=lsl#9b(v3>9v2Jah(>u!tK`4^!E|{89kx`JYobr_K=+u)jeE@2_8OYIQ@0#%p)kaBW;mA zD|P%#tk6(fxScP8KL(~E_Kko)Z<~YUZNMeI&L+?A3opMwFVM#dhFHE`2cRb#E*rdw zqEja4_Sz)bAi|=vhU}8K;Xv|g8GzMqwK#!}ymlM39s}6- z%`Lw6P^%ow+5=O@AqoBoYoP86E3G6S$fajdmO!QP`#);$jiiFO2I7Qi?2K*~E7YFb zBA890h&vt5WPbyqnKI+PfjKx*aVwiiG~mr7NJvH?m)Iryk%^sCGS~X_b$Xf7kH9o83 zy9kZ~mj`L4N#M=(NW>CdPZq&nt5dottSXfca$`sb>Af5;ttLVs+LEe0NIvqGQZ zIxJVeK~(|s+o@7#JFOHC*w@N{abP%8xI-pDuatpK(RIs4GvH%Ri(zRS#R{A}4 z2DaGF#wDV~5T_8JaG>&0SUO18*~)2I;e>9Iy;YV<&&g1HmB$7=1Z2Wd!voLXW8)UD zdxozAE2i2eV5+HiFjZlwih&o2#dT&e=7>2U+Wwb3Z%`5XyGT}hT0^L&5U@WC`9)8tRONLmO>8FrXQjzFb0G%F#Aht4ugAVXDi( zah6xxAAFf+Nq=KlF-XuI8v|+RX#l#UA%T2;MNrtnD~7PRsH~$aE8QA5U!B?zOr#z& zL`zzdcws`6x9(5bU}Qt~G_dN8bvjrN6oMyRA)G#oO}S~xU{dw11~znLGZJZ!-adyI zul>wu#11orAAr2C1o#NLhe|0!YeEI+^>6CA4&|l=fXgKPQ$t~n&9PJk4bqP>hyCw} z3|ftf1o}4l^7kDT*`B7HMPMWkFG`oL*7yG4XCAK#QVDld4t;m!(zM|_T(YdJ{P$IP zDXL+0?<_sy*SIQC+3Gc&0^-^Wr;Y+DXl)}UH$==b8qqMx>ZuW5D_aQ0(Cmn|0EScY z8yh!-dUnym`6~~OD*IbZ%}cbe01U+PVDbUH%nNJ_F?+s23gbnznL^(*+3`Ye^bBN4 zsS1N`I7mAZYBa^;y@;*Tn}LuxLWW2GDs5`HwEPmZF#dT8eZVd%NBsOSYKY2$%^=(A zaYa-8L(Y=qc+A-DNj=e^(g|lR*Owt&D+LH&x|$Z0Q}Fkw6hu}CAwv3@tOT?T+1+UO zSZ`~P+5m`xmCM26*1CZ*aj{>Rh)R2!jlK*f~9 zNd|wrj=+JJgG&?G^OP2NZKlwfjNVhK@8$JkHhwk9pW^GKJtFdIGH-(+!dO*7dYW;6?8R;7fVzD4ENG&f4f2=;S z{is}aK;)}rU#5hLK@P{iRBw%N%%dMT{QYwgA`f;T?UC8x?*2*4+!)-JBcd+F&wpg) zdtp9!0?}r(rhgl}fEr6lQR&`!1&$cFtRxe$>@l9)XdigRk)&ixn&+j9{+Kl633L`7 z-IGm&e7wgT!jC5k%&FiR*$(oLV~-=?ori-*pH$y_X-rVKV7~d&(!&jwKQ8uhb6TRg zTZXPE2>k>*i2lA&FM?OCi)?Okwv@#J$&4RzW{ydlcXYp|nf5)Hs%AJeELH_wD9?;} z@F>+C?}{@e&9m%?!2XikR;|OY&x#7jdg$(qSlOK}m31)Mei`6|?u>i#@@PKn40v9z zUOXa`4Twi&V@;{FzFOX2+$=BQD5op{CJvUL&Jqc5GFfGjXM3skG}gQSa_p!VcPN*g zvVyZ$ZKAwa-fMIxpTOL!y`C~DY1Xi45W}skJB6l7v9sQZ%}e!e7LFG#HNG-J{cNd7 z=^1|R;gx$ej?f2*mTT05p&ZpWZYJU~4{V;{w>Qx@u-W4*8zF3pw!qsJdZ8k7HN_Dj zJ}l+BxD6Y`$MWG^U+(oIeU}$=v4XhkYlN)$!cpUCmo#>V7PB8`0iPc2I;;_M5KJXG z!NM1yr+zLRW=@ebnFcRj%VO^E5SW=}2fdzS25u0SVa&FHU(+|8XRDx&%~q*`b!b;8 zs@r>jR-C<-(Fw*QBB9Cw`>*8O`3>bU$sC@>x6j1CBTlFy9?yj z&?vhJeI3E0OZ7wmW$7jjL3$@igH<@q^IuZAQO%Bs$+dKOSn74Kc)K@Ean!=I^E~|c z0Pxz=-HwfDkF&}ML%fuDrr-zU+d3(#CFf(2jj-zHhf>zIG1bWMq+66EX}dqjOZ z)%ELJOb3^K&&G3{HtUEjfyut*UXVp6`I47-@An12KHMp#(FLCeS9kF7ZoB2x9(uL+ zyW~;q%GVg%_Jc3ws~=7qy{Mi1*#U;xYa_9~MjvbHfIzRR?wJO5{@S))$J%x20VOd@ zBNPM52~2|y6c%`Ud+h9=#mMoy*#$fFs~>KUi3xLLS=~R)Kb&}=$SxdI)c>1KJDK*M zAXR%02%lh1l-qO&bDd(&@921&M^-*hWpfAU8UN<&WzVkBW#M<0U(`SRvMTyps1F|O)Lvy5f>0>yK zcsUxn&t!4+D?d%$HVOB3NP~p{?G0UPeh4lm$G+KK{ri3q+`-QQ6^2(FAQ?xpx{G zbZkz@YNQhQz1rh4iQ9QS}6NsZsaT zK#ns^&?f>$X!vNt_^E$#Zk1cFr^oaA=Q-}bs}mn9qDK=A@wBo+@jc}+{eoez4#c)O zT|F2n{H4(0%j5ARNZU-WOhasG%Xud#sa3x}cVsxUCb3bZK1qAinrel`U0hAb?&`AF zkzJZy4em?+MTtsg#QQgcOcQ+RzJ_QsRPzv`SGq%o6%FHXyv1=9V(I)RuSD|51hC$eeP^;(jC0I})ysa9-wJT3w80m% zqrNMmAd_1S6je|6yavaTlNq&^upf^!X^ysS!e`vn5Qff=Mdk4{0T%R~uYWp5jk!0)=MN~vyC zu3icIjNtwQfwQuro8Eq%*>cvS!T?Y|tTYM)os^-^P2LFA^5R32gU}s#3_6QE?k?;{+8mTi zDE_>x5dOe07 zcPmduEdTjuQz$xudMD@Bm$??OGH;7bSQfZK1ukjjXG=%8g+V#4T-$9i!kj$?^eeg_ z>>ow73s=tbvAY743C21U_Gpzeqatc-+5E$6Hh?m^f~Puz`(fdspDE&|iXC=9(rz|i z9G69pf%_3BXi%$uPGkJl+3(^exx~%!Z2Ff5*>x<2K}lgh9kaER0V5Smg3e!;D;Zo0X0lpB8?0u&N?we`j{d4mxJd)gE)lYPCL2G|Mj&1ytoW}{ zj?{|L1Eil;jZ=1!9cwZiCv}S2;18jvRH_^ z+GyX~!b0m5P<9 zILZV(o4_)09m7z2+(Cjl0;;FqvMOeJ7*Y3gf^W92W|D?fa0;H_a2`#nFw4<#hX;$% zWRezcJM05@fa601HsO@wbCr@4)O!S`0dQSTTW_no1AC5^+yG_Y>~RzpJxiPQ`s9Aw z%S^o58=W(lrDRG&yYfE1hn7vxYG?33GmU)IO+4BtMdOb2U8$0nmQc@GU>W%rX6`&H zd2{;QTUxDq-u}EEsmW09)HBJ_uMD2e zd1usmTm~Li-g%Eey_K<*^JJYeW}n+H@*&-M6G*HXS$Q0)7V*WgILlQGnm|~6&mbAT z5?6(!a)I)P`&5KxkoGgETGA-tSp~Z{54E-Iur4 zk7vE*IO*rZ<_|?d!17y1t?$?E}ia*|Lmf zh_b|FnHb9qGj@eYAsH0nB7{Ur)?`Vx>|-mJP?CMhQdvU%&UA0hIO)D#_c!zU&Uwur z=lyz~=X<{2^Et2cJddrrvvK7Biu#R|Y`DYB)Yb9AVGiaC?gr_T)SHJDh4qYLQ+VF{ z+O<9AL>Qe-PPcN-`z(|&8X=nB)n5@zt2Cr`j3y$#&VvUtiJdb*eL8g@|9Yh^+Y%d< zZ%}o$-$><`gB!Rts^HDFhUZ4mx~sb{1e_V18by2s8{YT>Iog@P;HAvf{`q!S z4}Bg(dcj-r=wVBzi8xG3qN{^YVlLM&4>#5AgmoFWcEBpkC8YHF%HYKnyhempJsOXp&Yoz1L=qwHWj$BoKX&TZhB}h zB=)?W1CnI*^>}rc23a=k6{Ya*TJPdH4WH?HRE@KJ@7WU;4!$Q6dwW`LrPoGf!b8oU zu%2!Z`1JQy2K&XYyWFwWMW_|^m3uPx&#~rTXUpf* zQPF_nJ(U9q*Q)~F3KkSf3@N9N9%nUr#XJyyyE{ZlYrs=3wt)5a?im3K$b-B67g_vXWO}ybNd-OADAeT6*qt^h z*>^dh)^Ye&82c!L<&JSgx@n;tynto|c36K>qzq%)nSh2FS}Q@?kJke@xJ{#a(z>-XyT4HF zZ{|I068~xJwc-oYC(5bGc2Tp)lS}ZfxmN5sC~emLR4aGPz7*G*SDNOcP4{xYLFuy1 z(;nl;fQ(mkRG5JGPRG$A&mr8j=B-^3w##i1{Sn=2Vv4I~Ms$Wci>IDN7pI<=L@>5p zj(uCan!(K2zRy?DLfcG-@8;ZPGsL<>U0Ia1vTCzm=C#@D6yOSo>G}dKov4pV8Y!GrR8myCl6et-f351ghB?s5ti( ztF^5i7ioWg5VLi%x)vN3f&e#2Z6w6EbH@PV1ib4(+qwA%AhGc8KK_@5x!7t|8m?%~!X-cLJ%m|zYZ1)9*;C4P3VZ4yI|HJln z89gc1Uv1uT?cj_!zmMt97G+o|^4SXaiVg%LW zX3F8yU7@1$m%KK;>w52zi(mM~Y%o>(W8JN=O7c!!zQQ||YXf}B{Tj7ISFKW$62!X; zx_KVDF3ycTpRN)Y55v^a;ma(9BHy~m1zG51hbF1K$ZA)z)z{0{q&K?~kuK)0%PNZ#xWC3I}cIL31cL_QVB23j*qss98bz|CB`@*N+yhJCmne= z;4Yl(G4^j!ZFC46VgE#+2k`*^d1t{eV)#e(Gheya=lw_Rr|97$eP^&&tCW@Nik5nW z8@MM1hIni0s;R#T#Q_y!%0q>i2C7C|pTlNi5&W)1srekPJs-s8UGU9SOm02|1iWSf1XneeLD9?45&t@p97L_j8 zbg@NO?2gQb+`bi$%C+5wK@__YX$E%$JsUL#X0U+<0@(@7AOeeY!^udKesgXjjO%6h zG(FYt^L^=lsF_Dp#<@N$(qsdY^r2Y6`dhZtj@P3`LPdoi+<7M5R!Jx(Ob zrrotXRyRA*6!Q4At4U2JBCjll?enLCSTvpPX5-Nt+78|pbAGtAW%k3Rwi)ZR(`}!< zzjELhB9`=Q?2mI-SJ71AOOLF+CDf?+ezfry>QeV&CP=qQ_NETo4%-OKT39FIef($r zVnm;YOnOe4c+a;eSd1eiT}BokjlTN0z)0RaNzsR$rAm_nqHhR9Ul@q~&kEkZr@u;YWhuJfHt6Cc*wqlE{V)*d9t6c{#={KN9y5TYu*&m{Ym+s`~JLF)akP;lz6Bz z_(Ra;*q(PaM~Ix%C?}k12VH7r!h`P96C21)56||!+C{WHPKTzX*qv@HPVZ%D^tz(| zJ~Br%5o_B#J%3ry(w^m}68?yJxpY5UpPcxZ`sob#SysA%1N)ZWJXaa4Zdo%3D~`!Z zENEU2(USHkzmg$4IsE5I6f8?7A*T{L7(mR6&$)oh%rbal75jkij{E0nro)j6)yK|9 z&O>1ga0i2xM(tf*V>c65xD3`W*{ykJ%@G;XlGolH=!K4|FLvC2VSh{|XK0`+dWc7? zV5Cp~-d=R(p}6I0Lz4~SW^x(F$Oijo*11zTAu?q&@)dfuw2^_QDwc!Kyp&g5Z0wPd z{YtYn=INF5*Ghp{^8;}g`8DRhD9Gi|O4Ec}gY^&Qp0TYfYqN57-6bUX%J%LvMt2%;p@q+} zjxQJ%xCRK}B<=Uq+uNZFF9=56hAPcq_rL5E8F-fxqj1hM$n(-d+5C#d?Y%|Fl!vXM zgbR^lx`t1WRr2PkIHeUUJmvi2bo%wipK}M&2#g}Rw-uJMUXHwa8MSNQo1<@XU-;3( zVN`(H{CyiTz;_4`gMU8%(v%Dg0(Z_L5&i}x0*9ZV?fPheQQ#(26ez&<-=em(qyjU+ z-Jkw(5IA^JQ!uymf}&(>S6WZD%Yv0+nQkq#r-YDJ_Xenj^R-RdGJTm1%4F literal 0 HcmV?d00001 diff --git a/matlab/rscale.m b/matlab/rscale.m new file mode 100644 index 0000000..0ef5e32 --- /dev/null +++ b/matlab/rscale.m @@ -0,0 +1,43 @@ +function[Nbar]=rscale(a,b,c,d,k) +% Given the single-input linear system: +% . +% x = Ax + Bu +% y = Cx + Du +% and the feedback matrix K, +% +% the function rscale(sys,K) or rscale(A,B,C,D,K) +% finds the scale factor N which will +% eliminate the steady-state error to a step reference +% for a continuous-time, single-input system +% with full-state feedback using the schematic below: +% +% /---------\ +% R + u | . | +% ---> N --->() ---->| X=Ax+Bu |--> y=Cx ---> y +% -| \---------/ +% | | +% |<---- K <----| +% +%8/21/96 Yanjie Sun of the University of Michigan +% under the supervision of Prof. D. Tilbury +%6/12/98 John Yook, Dawn Tilbury revised + +error(nargchk(2,5,nargin)); + +% --- Determine which syntax is being used --- +nargin1 = nargin; +if (nargin1==2), % System form + [A,B,C,D] = ssdata(a); + K=b; +elseif (nargin1==5), % A,B,C,D matrices + A=a; B=b; C=c; D=d; K=k; +else error('Input must be of the form (sys,K) or (A,B,C,D,K)') +end; + +% compute Nbar +s = size(A,1); +Z = [zeros([1,s]) 1]; +N = inv([A,B;C,D])*Z'; +Nx = N(1:s); +Nu = N(1+s); +Nbar=Nu + K*Nx; diff --git a/matlab/simFxFyStage.m b/matlab/simFxFyStage.m index 1af6fd5..85e6f08 100644 --- a/matlab/simFxFyStage.m +++ b/matlab/simFxFyStage.m @@ -1,39 +1,41 @@ -function out=simFxFyStage() - global Kp Kvfb Ki Kvff Kaff MaxInt - global m1 m2 mot_num mot_den Ts MaxDac MaxPosErr - global A B C D - function ServoDeltaTau_z(motid) - Ts=2E-4; % 0.2ms=5kHz - MaxDac=2011.968; - MaxPosErr=10000; - if motid==1 - %!motor_servo(mot=1,ctrl='ServoCtrl',Kp=25,Kvfb=400,Ki=0.02,Kvff=350,Kaff=5000,MaxInt=1000) - %!motor(mot=1,dirCur=0,contCur=800,peakCur=2400,timeAtPeak=1,IiGain=5,IpfGain=8,IpbGain=8,JogSpeed=10.,numPhase=3,invDir=True,servo=None,PhasePosSf=1./81250,PhaseFindingDac=100,PhaseFindingTime=50,SlipGain=0,AdvGain=0,PwmSf=10000,FatalFeLimit=200,WarnFeLimit=100,InPosBand=2,homing='enc-index') - Kp=25;Kvfb=400;Ki=0.02;Kvff=350;Kaff=5000;MaxInt=1000; - mot_num=m1.tf_mdl.Numerator; - mot_den=m1.tf_mdl.Denominator; - else - %!motor_servo(mot=2,ctrl='ServoCtrl',Kp=22,Kvfb=350,Ki=0.02,Kvff=240,Kaff=1500,MaxInt=1000) - %!motor(mot=2,dirCur=0,contCur=800,peakCur=2400,timeAtPeak=1,IiGain=5,IpfGain=8,IpbGain=8,JogSpeed=10.,numPhase=3,invDir=True,servo=None,PhasePosSf=1./81250,PhaseFindingDac=100,PhaseFindingTime=50,SlipGain=0,AdvGain=0,PwmSf=10000,FatalFeLimit=200,WarnFeLimit=100,InPosBand=2,homing='enc-index') - Kp=22;Kvfb=350;Ki=0.02;Kvff=240;Kaff=1500;MaxInt=1000; - mot_num=m2.tf_mdl.Numerator; - mot_den=m2.tf_mdl.Denominator; - end +function [Kp,Kvfb,Ki,Kvff,Kaff,MaxInt,mot_num,mot_den,Ts,MaxDac,MaxPosErr,A,B,C,D]=simFxFyStage(mot,motid) + %global Kp Kvfb Ki Kvff Kaff MaxInt mot_num mot_den Ts MaxDac MaxPosErr A B C D + Ts=2E-4; % 0.2ms=5kHz + MaxDac=2011.968; + MaxPosErr=10000; + if motid==1 + %!motor_servo(mot=1,ctrl='ServoCtrl',Kp=25,Kvfb=400,Ki=0.02,Kvff=350,Kaff=5000,MaxInt=1000) + %!motor(mot=1,dirCur=0,contCur=800,peakCur=2400,timeAtPeak=1,IiGain=5,IpfGain=8,IpbGain=8,JogSpeed=10.,numPhase=3,invDir=True,servo=None,PhasePosSf=1./81250,PhaseFindingDac=100,PhaseFindingTime=50,SlipGain=0,AdvGain=0,PwmSf=10000,FatalFeLimit=200,WarnFeLimit=100,InPosBand=2,homing='enc-index') + Kp=25;Kvfb=400;Ki=0.02;Kvff=350;Kaff=5000;MaxInt=1000; + mot_num=mot.tf_mdl.Numerator; + mot_den=mot.tf_mdl.Denominator; + else + %!motor_servo(mot=2,ctrl='ServoCtrl',Kp=22,Kvfb=350,Ki=0.02,Kvff=240,Kaff=1500,MaxInt=1000) + %!motor(mot=2,dirCur=0,contCur=800,peakCur=2400,timeAtPeak=1,IiGain=5,IpfGain=8,IpbGain=8,JogSpeed=10.,numPhase=3,invDir=True,servo=None,PhasePosSf=1./81250,PhaseFindingDac=100,PhaseFindingTime=50,SlipGain=0,AdvGain=0,PwmSf=10000,FatalFeLimit=200,WarnFeLimit=100,InPosBand=2,homing='enc-index') + Kp=22;Kvfb=350;Ki=0.02;Kvff=240;Kaff=1500;MaxInt=1000; + mot_num=mot.tf_mdl.Numerator; + mot_den=mot.tf_mdl.Denominator; end + mdlName='stage_closed_loop'; open(mdlName) - ServoDeltaTau_z(2) + %ServoDeltaTau_z(motid) [A,B,C,D]=tf2ss(mot_num,mot_den); - - %mdlWks=get_param(mdlName,'ModelWorkspace') + + %mdlWks=get_param(mdlName,'ModelWorkspace'); %whos global %whos(mdlWks) + %for k=["Kp","Kvfb","Ki","Kvff","Kaff","MaxInt","mot_num","mot_den","Ts","MaxDac","MaxPosErr"] + % assignin(mdlWks,k,eval(k)) + %end + %assignin(mdlWks,'Ts',1234) %getVariable(mdlWks,'Ts') % in global space call: %global Kp Kvfb Ki Kvff Kaff MaxInt - %global m1 m2 mot_num mot_den Ts MaxDac MaxPosErr - %[m1,m2]=identifyFxFyStage(); + %global mot_num mot_den Ts MaxDac MaxPosErr + %[mot1,mot2]=identifyFxFyStage(); + %simFxFyStage(mot1,mot2); % end diff --git a/matlab/stage_closed_loop.slx b/matlab/stage_closed_loop.slx index e11bbe1e3ae6fc81572a2f1516355fe5070055eb..e594048ee4a54a04178a5195462d7c5081ebad1a 100644 GIT binary patch delta 15607 zcmY+rb8zQP@GTr$8{4++jj^$_F*f$b`NZ}&Ha6Ut8{4*R+j^hhy)W)NRWpBdO`YlK zuCD5;(Nn)z=sS_teQ}j zL6y?blhHyriKWylQkV$TDb8I#@G9w+{#)VLTe-O${98!d0r^XfHk}kbco=d-LvvVA zY-+8sOCqOM!J$!dzkVg<=e!(iJO9~Q9D9tKhR@wgenXpmA zD?Acxfd^Yzh*4>(WM_8AqlASS_Xjk)C~SjjYPQg*P}fLnDtC~4Yi|1V&83X&`X5wJ zrw6UHU9@UMAx<-KX-NhDPih_oTfc!uubnUolW_y4n9_k#rG_aH)XZ=|<`Pk?KD<*r zz9MO_?*8FKLWGr_$|*b0*~*Y-#w>a1^^13#_cK$cHZVku*S?8u|Mr0VHrtpm1HYxN zHvb=HYgqbs?3$s5q1Zj-u=|tIe|IuZp`oGOZtYgd-iq5r1UA)t3E(;r*Q-EBhKf!T!N^RKH^8J$Ejl?P? z<`fdPGKcBYKzl-9$P9pFb9PGq)RbK=b=#OBUfNENha?1!6G|@xW$|k#CDsIsWdYfl z^YaHA=jxt^Oq|z%hK7c5?r}EK6sLPFh|*7SCf!&_RcwiOhCM|S-)%H_p>AGf2LC-nSca{1G_27FdVmHik>10G{qbT4UJILiD>QN)t z(rQz82`OoPDd!`VPPJb<(yluCEKA_w6$w`U_!KH*IP=1);N<0Qdn-q7dn#^1GF3w( zQ)~CpeYh5QH8x&LYk5!1n##0<7Gta9GInZdX|KQYY$lHzX|m4>R;bngl*fTQ%?Zo+ zW4~K42%kTeJ=dRPQ^mO^_ZJqG*rV&h#lqx`*4Ta|h;*mIA={%%mZ&rRo?S$uDf}7V zZy}NX6A`6;JsT9Gn9UO1%QrSF_P+X?vL9K%eenpWK9I@})_eSlCw0QH-zB})bTI3| zDY?E*KrZTQOd7kgv{F;rK%hg#=M39DnJ^fWQX!w?2>ditB+H?A ziRA(GU{)F7UmN>NvQrvYr&2lnis8dj_V^8C?dH?#&+V$|NeKNk^s)>7VZ{Jz3tWnU zwFbjNF{wCAGWk3hg(1z19c*gf63H6g(N>-s8tXwz)4yCIukY^W?!G-iWWT2^HQn%c zC0@UBsihClV0G13z=O{zM+$L{4GFQU`eOskP@?HzrqBX`eunC{_V!f{kF0X})fQSJ z&jL0)IZ3@)m;Dz>HP8fa+}chcb?uhcn(N}M%WTx4Cxv}-CL(EvjT~Tg zv+!WM@aH|hVCLMt4;U^Fb%7=W3)jn1mI|VfQ{eITRj#B1J>N7}d6&Q03RLs8le?SZ zAbgk5z~i?Tm9Va{%_9_Ks;t`bZL7f(BnK0-!ATXjV_GU{pKX0p&>yC9@07Q6Yq-$O zXmHzio`=UB`}(U+GWPU6?;ydoWOHEdYcRmqrXM?@jbJ5;$pqRkF@i3uPs2-usH=MB zAtXzzaOP52_-Pv#W2`k-mV!j2F2@x;c3OU>ZRDc*P&vp3Z({iN{s#!}^-W9iYZPl0 z!Pw2B?L&EEMT&ZN$kxs$JFBbFb>tq@+j}|r!wZgat-EKT@sF8~$vL-^1|o1z1rcf5?M;?DOGQ+(rP^v)S}{V>m=%{~?~1qF(#)cS>~dzu_=z<_!9faR;o@*#-H71I+QtVDGr^1}R9%##7=E9XB}HMuYiMX65F2`x4*{jVxnr!vNg(F)FS;MN8AR z$t>ES7|(wUNXAhXP~=&e>`~VL+IZsPEZC0vJDnX$(hT%3@ zjmqr6zven+W6Q(Uya|-!NzVR6NRurf@}hu$`ZLBpe21m(E++35Ih_V|u8Ts}`%vVc zXGHF#t5#Tgw8m3axgh(a*CRB*Bd|kZbg6pq@?|_-h$~3y*l$fu>sF?U3wSzW*`mHh z`7xi%igu^vJ4r)n)@yV&$>vZ%KU? zOZ0|gs?>(74gkkLYcsPy3)OJkfy9?*S6x8}MJTcyvvVIGIFPX?=@tn`2P_5(g1q6L z3ti6xYLE~ZGFSR219~K!4PT~VW(kL#$S9ivz5@Jt48LTU+OTt6lE1EhE95VQRIA7d zXKdlv927aEi*J#&!$x|9vkmv_N+&EMRl4#r0gYxqn# zK_Wzo>3;0^_T%hUFx-cBGH%jcla}G)DXkJgdPinqbHY&Ts!HktU2_eNVK}8PNGMg% z5ql!p3@>xZw(`F76;ADk#$r-VLn)2)9DZXXCa}fvmnwM<&%1`1{LidW)La}g0mdPp;H~Nrh5(A&0u0HxROlKXN#vZ^qp68^GE3qTj;3ZS z4)%64`JsoX-Lw7f`QA`ULS^E&NC^M#O_lbo7^XJN>|L4sqBlv7it6gS0IS}7_^ZIf z7>!`wZQR$#k@4hYtim>sTlT`HS7J>ATcGGkOu;3Ec6Brf%lj5z7jpuyA?Dhn0MLJG zj`zK!ve9FrmNylL7$PAdAve0-TuYNZWP#1!EVnGgiS|v=CvW0$UBYh5cI#jz&gy88 zH@kDS8cuf3JDX<}zX-`GJmQpnuV5Cy@@l_WkC0g+d;F ze8a-5+m^@s-r=rbW~eVX4Jfy)2XvSi+x$K;OX9-LFE2LqCvj(~H&Gh?B1xNbGQmNg zoq-9zwo-F2%4d<$x{Y&T+|Kd*ab}d_ojg76X+oWyj1&~=b@)HBE z?XG6R!j#n3)zzhI!t(!Z#@8@o05ko4$EV533M}ugXykVSly(+HZkI9wneD@IsQ-?? zNf1$(U>K*)bZ9Kk9y-yBL3;LAOQ+h!;*Wyi>6K}ZNrD`853g);yaoQ5c6OC&BCMhm zGUpe~91%*UPqC?U-;Y!jZlHZTc7sHi_48W8Al>#zfo#^w$*Gq8N99!oxU;3EmKN$( z!M_m(TFT+U=~;2XL&&!I4@i5?9~Fi~Z-Ehz`_?hr8(Y@Du10OP+^zVQiy)j87mB*t zrOA$Ly16T{39wa_WIIy@W;W(9u8%CMMTuZC3W>4Xnh#ud`{Bc;l>q}xyo#Tf!7qi} zA$4b*Vh{mnTxq_w{JR+S_03Q4FfcHXHy|)Eo!z1GZ}s1kf0G+6&IXL#$LY-+jp`+y{2D#ma8pFg*qbNn4`~|}gUoe(k;NaYP&r8nJvW6yn)~71+yAFR}GPn8jZ^iAnRW` zDTSLsAOAZc+ON=Ztv#YDL<27bXeXrVpW>Fk=#8G={s?)JS0wwvMTvZ2(g-C=j&fDS zHx}UdYS3M=B@GkR-d!;(MM79z9I8UY-a&USHm99Zu}g>F#_1k3x5 zH?!J7TgyQU(p{aFlI9)QJd9fhMM}3#r0fwQry*Dg5p$nT%-j^4GQ3iH>hmZ98@9-K z`3@cv*JR4ev2IC~m6ZZg9ntZEDZ>vNMy1Kn9@7~a_&Ag62h$y|B{Q!uN1>M{ilN^t z7%Ke*A%O@2{5pPBr5%Ny3#VmE!uET@8T)cAitRrJV&R#PU&V6a5GIsr%uX!N06;mKFkjye>Z!3}U`_md#a zp45C@a}1Rq<{D~2Z)GD;unH-k?q7NI`&B3!D~K$BGFJVa~s_6s7v2xnIk8s@gb2Kf2BXfxL3L)qxcx4CTHiGS7D#- zOwEasqOsCuY65nj^!$e4)FAi&>L4qb54h1APXGM;bYBcFgWItVo^*sG1o`*vWY|#U z+sh{eRl1Q5Doj_mF2W!j{Z(i%@*US87a@=oJB^+6G0P1F2C z#Gx=Xk-p>AOPXyE+1j~ud=u!BMXq%anCZ=<0y_a4KxB*|ufyjE~7C06m-t8%3Wlelu=4UvS$ zG8e@*6w=~q9AftAUYqxU;&XD|W@0CjvGRcj(p%r9gLPlXxzDuwue_Gmn7-jyN=3hmTV%(Yf%2Y;)wT?AlXlKvQF9Kd9vcHFa zZbi&}94x9#=w-5NDrg8HV3z^}AKlVHJ!TG1y+XuK$##i6W$ZHma?LZ8U^{Zww04ot zQQz4&s{R1lsk9Tlqu*=o0_O9^@8@qGFk)La#o&Rnpy4MDJxp<2Kc@x70Dv3a1~Cpk zUWo8}zu>{j<8P)`V@n#i28N;Vv}=+!73{7>cRYbgyEC7FS*kKr~O!8(P-#YN|F+i2Yk$5Qcp zlG%~oQpozoM)UL&Sk_olHziz@uo`Mnu!7XO%0W(;4rrqC`Eema3oH7bV`5PWtMr8r zJ2$h}P+afF3MzU$|DoV0$oxM;9;JXB|L{fK(Z zP2M_MhFy8(5CVMP5qacn`4HcBry9|hUHfVBSg!z&g)AnJ6T9A+)vw^3vJz{BHBmuG z`4lI*rx?^gZMPkTls5OP;xZ0or)9Ge7sR91KQt(ki!F3+>Ey8>BR;7sAzS4OzCV0u zm=ZcCFNtGT-#4E zOM*aK5fzn=p`lxQ6vGWCA2ZnewXeT2>nnnC6!`K@FfXfq{TjNnqOvSGh zM9G8d7W3goXN@KN%-pkY{@n~PcFljQCw(h}>e#@Ca;00{e+H7(4+%O_ZQevQWSn`T zU~(35bTWLZ*X<$3&7?Duu6K~cv znX=(p&7`(RXXENI3yaXZpP)H0`?3$El|z-&f=b}+ z(W0%OprAg|l`DLhvSBa;B=&>FvOWSu?cY$cZC%(wVg&q>|8 zs;BPLt-~$BX1fk;aEix9*g0Uxd_OiMYGV}byI)x-^5p}nDN|1Vz}dBlCgP-|GL-2P z&u%tqiG20cd2@wSA~dh}&e;M<`ka61!YRD@_&nF@SXMAH(8ci2A@2T8R|vp8c+OQS zDFw8(w#<-;=W5EtWh6}ucms*x?~8owy-$n(3@{Fy@ zEn$rdxA8-}@O^t~-<^ztq*-%J&gNfu)(O9XRNEyhaZL@zVu!6CcSUZpJIwEIP36fp z79(?JS8tNujvUQPFCYB$8(XvAD~@|-x?FjjUB4hHdCx^;55D3bzPya6Q&4u&Z@lgN zZ6F~!kshMmoK{ibzaZ=BN$%iKwlC{4e-nX$A$)>?;eml-;|PF@pWg0k+;v|fxq;H& z4oV~m{ZEBXZm8xf5e^rtjpu)?lwwgO#o3oI%c$tqMm9e#-NA(Wab)B)2}<(Y{>tow z;UXRKbR!M-^n3|26CD!21YER;Jp2P+GmW?LO1WV6JFV>VG`#qjIv(grwR|=G!GFSRv7{MBQ@xh)a9t+2dDuLM;*o#RR{i5??c~4w3)>f!_U`kj zWEMIea#dO2>_A(Ewrfr@G4(w=zBN3Kc}4iF z^+REqPuV^J;~xD~N3l7n)@ndZ3&TU~gj14+%T#;YHH5&@5Zq-@z4P zB{Zxkkh`}aiW zG0aOiuI$>T)9f=uFaO`QR1DF;>}hf{zvRY}uMnIu1r%oE)OQSE!>pRiXL*VHiAnIC=Q?MPX0;vNMFiqfN$0 zks&D8mz-pVA7$@gf#Pcxd_3)#m2SLwKtQaS;y*Gw;-byhv)x@Y{J`T6?e z>nBOu>zU2QtQ&Vh7ShdD@xDY)?o!mvf4%Rq81Q!Xteu;6O!Jzyu$6V(^@VT&;eJKp z*im7h_$uk+8}q{FyIe{5($<|-n&y=0L&MI$T9Yzt9(%yvxfrMLAxF}+5Lfk3AHsY# z)N0S1x5>j_^+ZvH^TGdVxE}k$>eu}d`l$2Kh~WC_%Ow;&IFl``7cDZd^c7JJ1ig!3 zAGpbNn)(56IGs*SGA}XDlpFjSJX=b$22#F3WRpHy@51FRT3Yr?-qOM5uN@K1df7t) zN#1fKlUwuapr$g8pxZzP^$WU=o}5Ac`DeER7Y^h4b zp9UTBqsCn*jBC!&*eZwG8tCc?vnb9cM}k&i1b>*a{}sn|pCJR7^ZI|c%VrPU(3px# z#j_v`QFEz)0w` z4VT?kAOMe9wq&T}wew=UyF5!*ug(@1G4D3S4Y7;eHYu&kpk1R)DBVdiGGSNXWBJTX z9L(L&>gE8O19(J@w;gU8wEfmr;t~t!aAZMtpP{sz(t|!?KEaBa0fiuL_ANKrl5`{Q z8S_p+@56X_7yWzoq}7QFaR_j)1~f0U!4vcN~YT zmxAkAbv##s)9i+Gf4sHUmj*1J-OsRMOwPOUF#syavwaIj^8%NR{9e9)Y^u(BD}%)j0HC`7sNkD;+CrLQ0wFc(fFdm?%3aWO`bKnQp9=3g`eHNI zO#f8ZNxqsmlSByjK6-KQ8D*WB*L$Nz%RSL5dFAcU*H8v{6!#S(4qDmYp%DrbIIa)G zj??DZwcs=nrW-l-P$qcH?HM_FM8zC>&8Z6iFfzLwGo^;wJV?nuJMD#`CZ|2z*v1s` z3UP0YIF)@J@Iq?`Xd}+Sc~UFqYr@iVx5kZ0gOC2h*1dO&u(J0``WaGw{TJsh$Yu*u zf5@KGJDUa2SSNF84+d#ZUuGr9g=6@ zSX>9vs1xyn%Rs!OnC>sQZX--Pp8%eEaz#Qfv_U_lqrSiHx>Z`|qD9(Oljb>39(%!< zLCDzsnA+)GKN&UjPhv*o-V0wl zQ~Bk<$o}D+W2*u868N4z;{@G%9AR}6ZCj2;gK!0%Qzwen8Co@ zX@x2Iw=C37P4N%*^BcC1zi?2aTRPv%_$@L3 zXMalFrnpopHdJ?O%niRChxq=CJrhBvB0@9sJG-7CnF-O{MdxpPy0Cj=NofZTIdk{D3W{bAe&^_^%Ao>WKvIIkQ5MJr6|x$-lmeQmR0% z+K!0qe1H0VQizZrcH)lq@1f+5d#tIeLXCH4Len+{h8H0o#A&-VUot6C(NH)JVF(I; zOK(;<-U=%QR}pjmQ$OBg25c&qBM0K; zr{(WXy+O6N_jaSrmTN<}S>yqBg^haLUS5sVR4bzWX8mI>EZ80K z<~z-m)y~M2HF6qP9Tzj*4#doosk(Tg_tKIWWwd1{inCL9-d;eDc#7!aU$*iOnS?gb z5t)W!7RDvouzq*jk4NgOnRa5h+q96Ou@T{Chtgjfhv#|dCt$-po3^11?#Q`&$x5+g zQxhP62_YCsv~BE2aDeTWgcRaKDd=2XU*UT2LMw3;cAvRlM9eG=0MM5QA~@WD741fI z#U9P&IIT^+pqxhTV!yDsr5U0yGs%a|nfi!gUwsqdU9%!&1i_AaEsU$45T~ZwBI`AO zItX6tQUZBaxmzN0Gok)^YqkL^_|};n%~ZHd)g_60$`al=yI9UqDuRQ(D#czsmGk{9 z-N1kM!{i~ZfLI*@fPI^bdgd{>GBeSL>fX!w)3Hn|d4;9grw{q?xj%h;MBmnS8YIg- z>GiyNV-X5sNM@@zc|li^QxX>^&c`df1}FQkDdMV@DPMi(KNv<{6gFYZ0>`!B6o9h> z!SBJOIiwC!G;{^OfgN9!yO_^r%m;;X=QALUkNg6E{Ytn2ZqqyTIQ`YwuV!if&DSaA zwRhIu?y9%Zwm(#aHUT3vArtC&=>%C@c$Vn{*?0n(7_#(AD0EcQ>)Q?o%f&>!oW$r6 zNvdm!nD-BIo@hmhOd2`FLRi^cAxaruBXD)}oo1YpAIcGVSk7f587j-}Y!6Is0ghblu!pnOAM3l?w{6NoGP!Mk+PZ5WcQpQEt39JS^SvxwK%)tBH!CNC2W z_g%PbBF2bEB|?J-9T?0Js!z^R znDvG1R5ael%n_#`DZS4+<(_Mv~$Z zPv)j-I zf5|(5*k;Ni&xZXY>VpJgO}BF+nL#JS4;s|=p~y+AQ)0~Q)awNO=Xil<0Fer(99}Lg z6SGpevWX#MReY_|2fzXyjM7f?r#MBK|>Qjmw3_R9PM{Gr5v6w|A&NGB&b-+Tzptx<2tY&6vLD z6kB460RH>uTk_Aa@xT7`poR8O-}Za>LeE&<56)KF_3weuRMq}4`uc=MnmQz%--;%y zjhQw1Fu^cK=9U)B($6ojsHw(!V~c|$cI6`^RStX=v-c@7_&qyk%Hfhm>Q|z({nOoXuV4mi`#MDL3MuWrg;IHPQ6sKj|%^J+)vr;lW#euoq?X6kE>_gS(?kp zA%?KI`aLE!ZYKw4ju2*<8eKNFFAGg3GhcqED~+*| zUjh12gsIj7ztc#X?p_3+-kO&>?UIY0I$v)?_CmVik+>co4N9;bnLGQ1(K5vpR zLVo!UuXohlh}G*v22_l;lgj}W#dU#$cM4GqvN3QGkq?%8#kZ0A$T)rD_gn|5Ar@+YtH=XBM#ct~wm1XQ zH><3np`MxQ*H^&TkqFy=(;GMw0rULrT>3Z(7hQP@a3DfHf9B3UbNrFiaK-O+tI>YN z3)@Ea?S!SOl6{nI6tWasR0AgMG&G1qGwpw7P(FxW1@6B>jxEcMzxMj7l&gQzZ95>N_Q3gp|lbyIq0X@vk*M?c@^J8 zSNH<*X_%$XrIhuT$1lD8M2yaqdtQ)|ZiyZ82fjV7vVw2XE06+KjGH&zGtz{Cn3+LP z1@GyBPz}ct%u?IS^@S(Thv&`4?Zo~@O>z`#gOvj}=|8Ac^)xtTQxJf%)?H`eXn~+sxSt zk}=ZH-ROimcM1Zdud&^H-F{uhHOAv^k5_T$5*iSF;eijCs?R=ZLA%-eLKKk@DE1B99ONCijNdE%Df{D{*! z{+D9Vs;%D2Vu{_k;L*eQ#}o{*xhG;FiI%-3o#{Xg&QvhyUq9ClZ{-Y7*iNv}3!JT| zb|t^5Qrmu1cpSaBp5@qY-0|o?ae4~RsBnwz1u0`;a9JG=zsDi^~t}CFd%C%WD~KM=$ImMM zn0WIs1q4#1w((*hgK61){8Iqt+|Ydn)EV|0o8dD>jQ!xji7(yn0=$rP;H#z{|FII{e9^>k zj6w;17Etb)xsGHxn&}vGF4aE;ko-<1R~m1`+jx4p`~>0+dcV?STnqpLpH#I9y{Tpq zv9-`b0izi9@Yqkb+CL_ypBso-eD{MrkMvR@(yBEBRJhDmq2Y8f4C+bs8|5Y{IY89! zkmo<~_Mz$GuzlLpe0;_|34$8DCw3>_Hh3ol>dtH8=Rd$0ubr_Db9Q2?_D44HC-Ke8 z3J`U$L-~mDTi!eu0-0{Zg|bO}-bSMPr%0L;Rqp+rEB)a{=88<46Jmd|WFKz?mfSfz zOdZ)R_qa9^)_U(+EGA~I7pYCxmzxyzHNewbp#T1WcX0th?!^T}aoUA|mLB(B^Y9q&6i5r}(Kk!VUf zxlEK%SNz_s^ouu-XL(r>L{gcKhRtEeVOb1gl}d!=K({+kh+e>!ssJZaYU?zCWFSj@ za^*Kdb4Ws2!uTO(H&OK96yUBz`nw@Jvt5#L>knVTjXHMgJUCaba6-be{cX%)Rl5Xe zBna%Qz3t|>;9cP#Rr&=?syd+r`a!6MA6_>h5c)DSeCZFiH6LMQDB8M16 zWes0go4XM?gI_4N9D%*puy_+_1#Ue>Nn{y`62 zXZZ&5kMq_)eT56!FM!~04{TG|zHFRTBWb#-?a@9LLrd$fVsJTY1sJ*;-C1K*0-C)? zYQ3EmTuTI=Y$t!Pw}&l3=!^FOc8Ik zXrz})4XXSEau=508>Nf{qMF`+_{RKVFzDV=X^`D_D*QKaM+xMT{>+ibOvTY-sM2q~ z+|p>p(TG5`hSHbqy132q*?>!NK)3b$3ZT|zcGeuh(3GgLo7b$Ug|kaN+IClTKqkb$ z?r5~(V%h(?ve6z?@(=N*Xt+`{)1ic#p~)=cy?}i8i^{6tp2w}7V2v8#oTRHPaVohsFdjJ1G~=_!kRVfm5tAuE3S6ch%@f?1t3IXJRI~rpGHLvvNa)=Qb3?Mp z|6Uu)$zO--g?j0{02kxi$K_DQ78Hz^A3aqrB}(v$@D+LD!V5LrX>tuMOWT|oa|Ug^ zj(0=Sq_>c2j0LLX)6(Ut`gK7df{PtiJ>h{caL%Oxi2&F(5Gy}gQLO$JJ&25`RSi7N zW%D_pwJ=|Uuy0#0`iL>e>;WAV_pBZB8&XwlX#aJJ;SLr5HI>mWviL7UV_NE{d0TOr z!od!Sk=9Lb6Vp7dePNtGlkH$_EN1OL^SvTQb4J>3M1eiLLAMv*akHACHMcCcv;DdJ zXl4`a2T&1rq$Sv+w#T=9ByPfC&@8y7T_feD^xvd)6mlbW^v{H;4XJ2{BHtKns*BKJ zno}))qJ~zEgx8A_4snLN(gv-rQ7rX@<{;yp`NPV3hC|}(RWp|4r{wpH3S!wu)nko4 zd46vhzmI*}DCmg!WZzgT-zQLLl2@wWJ4XSH0YDkIH@`77zyQ%sfDojN%Z2&Ombg#% zO;@Sw{ylaMtx`3**L`=>g)wu`i?WKW#5C9AK{bgGs*}sY?GI7f`tHjNtvTHNUkq3;4ca4mIFJ5k2<@de-VM5n4R&wRU4=!4B#UHU2gWHM^VrjF@7qsmc*t+~d5RZvL@I z?@@b-H-vl1UR0jUJ;e;_#D%K7cC$an2d??DRS@Qf8yMHZPg&a6I!)RfSX5N?CIFY% zDe+TFB8}d4pIK8cJOjT`5!kQIyYeg5>^!+1nS>N*0cZTb?ay^*^5XX)23o=L@-Yna zux4rjw^#`QquX`G9VV|P&@_@Ky}Kpn?x}rkcKoC%1vVi8N%ZT-CTeR)c3$jueWmnC zWt>fs_t)TVDH>5AGGw)F{`PL&J);B*!DQ&`0#22+!4Y4l5q zgNxdQ@uu)q?<&YroAjDY#Z+lap3ty8@;!GW70)oJs}R15DlC)s6zd0Ln%E#>+Eg1? zep=}~XBW{eq0nkxT3Dq#g4Hpa1y#d-f7fi{@cIVuIUP|(Ay+_1v3!!$GfFawB7{En zf@A%sPptz7nS`(a(~8!8tZ2~NzCVl`9y#$pV!NG=%QI(ToFR1Z1z?7lf2gkJbQb)? z%Gl@@?1Zge0kTR?vzNyHjSB>zs;Uc?P{rNVU^LsHQa)Z3dm4nY2~9e^zVAM;ZiG77 z^Jpv(mZy3gW{}}eStlZ`QX)&iLaRtZn0pKvFqHCs#m| z;TxW(+co<0_>;&tc4Q1m#|g|75gn2|&v-tiBa(Tkws74uGOxnTNR=kQtGe4o9Eln4 zphjlX(itxhD-J!HMaGDJKtVz%U;4X))(x*2>!Q{K{c^_g#WMnln@LI`+wR0pTiLnlmLkfB8sDj=pY8k#?cY} z|B3YfeeZu`70~=gt4o3_f#%~Rz?ngK@e<%YAcJ^ma34@{yaGgb5U3@H7!(qV28sy& zKYhpn2_*;;{(or&1Ec-Fpa0*ub|D~(5PDEo0zLRJ=rlnQLO%>t9YzMyO{5102Zbbx qA`C@Qgn*6! delta 15624 zcmY+rb8w(d(6<}gwr$(CZQITUcd)T-Y;0p=dt+;3Z)|_hd%jcW)Z72gbXQGH*GzR^ z{aY&oODG3JSCI#YzyJXOfdOd?VN_g^K+!c#LT3lWZS6}dAu6cpNMXQSMbn!d$<6qi z6c=wEd6e}key_71tlwUZRu|KBgWIIgcSumcBta)vGh}p8{!m)|BJdq$%W4#p)hWV5 zxYLsTBR`oI4qYoZ#|3@BOKwZUoMGeaD9g&_j~tTKPg_FXE(dS5t$FP-@dCXGn*Vc; z6Uzk99c37b`(+#-U(x>O#YQ*qyoht0Cmaky#-0Th=4(EBYlYIlpEZ%%4qh)aWaB-4PBj8*>in^~~~D#9|KUZiUm z{2p!lZ>fE%f*9Mq_4cgVI+QCkT>Nj1zF3IO$_imp_*Dk0uLFbb z6tc?dP;;w%FKBZb=ou&yZe(jgFg^LT8O7D&oWJPy#RYpJ1OxJe{ZLLo7%Y61hA zgoXkdm6VEt;)9c@1ezq(*`;5Z1eJ6}(1~hm8VR=!#ZSKS(y+nD%5W@`6Jz+8bg4=TC}zMJIW_UNCm;&(Cn|arJ#MtP!njXP z94SinnY`)0_iwlU6fg1qOBJtEhdZw7=`5_eCz{}>T2-c9nlXE182V>AU1YZ%AJhVOT* z9w@rQu;=g(4#R8&cUM{yE0eosM@8gt7xhgAwehc#W-+K&dju^epYsL!bPmj)xrTrY z2IH?@xD+FqTWvxsIWm|C<;_@3POO;H#18xk{d#>p9!e?_6o%=;h^%UO9n|C*Rn8;7 zm(2@|%R@<`NPeq!p;W371zOPAL`)@kB94n1yl`m&Zep`SO8qR>SsGX6n8-+&fV}%> zvpqh%I5*uahWYh=lRN*{fFU1z%Q!%87r?%?G3N7FFI|I`EGsV`drm*+jtCtX#npsm z8#x4%lYu2knxAjYV(3^Z&{?3lsi2{_u82gxi#}lS61@SZy$O_5U`b3&B!P;+4jbp= z=a;R$!rBEValvnqkPAyB&1U-kV7(mmJ!_Hl3=9;aO&4nD>`Ou~cFlS_fdDkDU;EQS zZdQ_z;T%({pD-4ruI{zZt1I+ND|0xPLMeH8)|||Lvg#&b>Z6K?^fIjd;F8hK(PIcz zzaa=_K@(AmDy2U;z>sD-pO>&#c0-kW;uiwGJ|tzqIL4l_@o9$al0j(@A~r3vr~KdN zQGqGZ)JyC6>h%OmdG%UtX28D>=TiL1ezH5&{{fUKHI2BnX5RQa?dWh{Z=Z!n9VKUs zF=+RscZFfu(Kw=@w*<3k>ce~oT+v(K?-bJ6Kkk8s@No69zqQtfTlX$Wo<1ofVr(_B zJn;GG;$)cP>$hD(vG!BfRS)CS?=9o=pQ-Nb7)Q)Tj7je;%l^3=|K;CnOSX^#b~joG z=!fW<4b5nwMq>>Tc4lX2!E6MO0LMNhKgf=t6xM43Z-hDp$wK|0t%Yi6V1UqX%GyKl z<_YiDgD#n^ZssB+#Mv|Uxl*;C@TKuED5S(4J@EHQ+DnFg&xTUV$1oRyyWVzu*C86< zDBll!VkobPH5j&B1tx6`g5q`JMfAsB*@bZ{q1>?{$Vmn6F0j-1#I9QceFRVbHoC-L zQ$uN1(LQ8fcW!F4CP)}?-Yfq?53D;QEv&P4rUipsAXeD>2E@UWcB1<=r%(-UYUnVM za!w)Es3MITH$Ok`+~+GD`adsm*c?#@InV=tu&`qi(Ui-)AHg}iA|N)Lwa$1y$T{_5^H zXO~EPKf4%+3IxfE-?}69UO;?d2kA3Q5HN?$Sywmcema_w)a=-I$mCg7eW+bhDiQpY zmqP9%Cw8U@`eR1^<){v=Ecms(i}UtZLus!LcX`39^l*Pk1KEbi{s-Jp_r85Xk5vef zO{J)POYJjha;egYFVy!t0uwbo$J+4UuR6Cc!OtdWOyK%zm>nM%&q}D}@%=$CJ8|w% zNq36Zuk(FnaYtxXXLTyn=v_F#+f#L>1YVokGj>AmR|f+>e=A4n@!!8J($ZsbmzOhn zgfrVK_Rf-&l%vwBqgq z{D^12>A~@Id~IPUuaab?>!pf1kwyKd=A5*G*+i9cvr@3YK4Cl)CuPrrm6G^^RaFIE zvnq6IXfKCcg{I~5pG%&$21iq$@z?2LiTOm`2H;_F+j684`#a;up@ke-_LHR-TAODV z#>oALZ~s-@v~dmhDa3iF4ZXGf+ZrrerVWkS+{)Tu$Xr*er&*$^w7*j($-XWoV1GBjD4`V?feigjuAy86Y3WO6BWxwfH&U#iDd&%)p zBhV1`9=iTJdJv31PJnv{M~6?dvB1FaU)rzqKg1E15tD;kZo;X{w!E~_L*4NjjmIlt z#Q7S3bZAX$csUZlm!}k5=vxbDdPoqH1%THJs{jFk#7L*DpI;)$&|BPe(7GQKu}fE# zWlH-i_tx-z-+zMQqM?zfr-%z|7kz$_`T_hIN|9WAxu1J;vjs88df)qg1sawXJqHFv zQ+|u>fiM9ta|w(FFb!4IxxMAj!e~Q!+lRT+>Aypbn_psc_#1!b2`OZUF#Jr?4l#1| z>awr~0sj$lAtnW*aho-}E~41HLyqY>Wl%3ZJ-LId_jJp=YHSF&bbYS}Y)MJSCIGYF z|NdN}5nM3d3o*Ij(FD`l3D@StkOg&E9eZKAAuVx|=jMPOU?_XFg}%Mt|CJ(1C&0sj zag{VZ>@{|f*4NwhYJot*$Z{Ai(-Sl|Po>W6NL*gle8&2i`{&)dbYz8ulv>_LIPqwh z5keQi8J1E>vbU`uLT$a$s=1vv>I|@xLuW0drp)bBl3Ua%s@UR%GUMg8RgQ7JzN%e4 zJUm!nHHB_dbNCm$p1?;~`Kqg17;G8l=dT4Wj(}_5%#f3zHMX=Y;gm)<8Fy^EWwvPx zP%n6}z2Ew;`STDG??-KO`IK26Uwy&2B#Y%H9R(xPrGL34?Oxl&KG#--D+5HrbmZDU z%(asBn01xB=J}9H_*Fyh$vx2yj^zF|WL7B;hN_?|(d|d4rbd8~kx8x+C8XrzW^S)R z%TKxd%|A4wW8_MQvo*!mRgiPLd!AOGV2BRZgr4~s_*HxW%u1kV4*A&Av)Vxx)O)P= z;fzgW6@`Hg%(nO#Dm+sQ_yRU)QJ?pR%j(=bjJlfQ;^LO;v}+~dd1-@#lg-WS3%sF@ zC;L&(r$PHf6&T=>n`@aHe||`Cr>JS(f9d z@_`QCvf0=t1Ft9r_`&OlOHjku_T^P!{xaqe5Pw_yzb9>a$r}DETn85RaB0Y8Wy*<( zK-c#-NZdM)4?=&{GvDfAk6UUy&xQw>Y{^nM3NHw+I4Gxi>a8$rz8TD&`Vo!$fQ_dw zr7X}?)_;uQqd>{7w6wv`ehuong@Pu7#4CVQE}awBCzdDT%y+3dxBrC4W7~fJ5MSI& zgdw4lePSjwG&GFM0|nAM;|d}tv%e^{l&~qmhDKiu_GM*d72AEH`e(U`;OCP&xMl%- zx*taYx`1cqH{3fYN9p9Knl(&>7s6w@tn@q}OU-mB4ANPKoSa*%!_|NC>mVjy(X#jDlV0I)TP*$U$3BK(OQ1g8s= zFFX`3tuq}b&by%kc~x7_hNhE>l+!tr20qlGXWs>7NJD5i5f*XDcvv_J(Mtfu-+|E5 zfj&-8!)cScw)Y4ZwDy5;J?cD$p8kyj0Y^0X*`Pd1eZkdhZ)S$J;n+RB%^m8RKg#p% zrF&ymOy>284p7@i2zsOSY+Y~V3%n1guCxVgnNrT#sz386cDPR{5d=3BBTguLGjj;q zJ7)yB*sR1m-|g%cwri)Q@F|fh`&74N2gP~~OY;H-Ov;ET1@G*R;NU;Bv{v|Q&%>ah zdmiYgwowDDG%?yKu)GF-GyO<`GL_@N2%%sFHK6H$%H2Mb zj8fkP2V0NC51}8`_YPf1#+j`<=wEQa!ops;y}rJdfG1Efuu7Fp+cFb66B7}IR9|;- zi<>sAnpimXNImiy&xh9eDd`Q_Vqwd7g)X6Pzw1j{P8pg$_sent?V)^e!xHfD{9MAa zRK0zwfRKQIz;5yg15lD<#iN}V}ie30vV35!IeFu2;%9Sk9?DCqg=tAAbP_5rYmNZ!E`8MvMt31)WzqU zAd8*${4z5$>xZ*yeyv|}?=aM`J@49!NL+lbL855gTnq_gnetXp7T?nb6hWwDM#siL z(F~tAVgGSUHE}@TyV>Oqb0*r_ZvSAi@n1*6bwE{46uvjLM=q+ZC2_X1pXT*1pH{f-)$q!e>FElhNrbHJ9OuR$OhC5*oHrnh+v1IrZMKn0pyn)>aJnK@^=l6`w(n z#PKUsZVN5WH`(?(+Cf@+uWdx1+Uu~r`}?;-B!z?`iUW7!)!|_(tS0t9@x2G~uomgw z#%$>rG;)Zt%Z@JnvsM>?j0mBeA?`FO=>16drFUN9N8Ij`G4+;We7JJ|c9(#Px#i8myij#wFasY zX|01{R`k~>@ekY=R_2l>sgMAYmLBH#`l6UWifF{hH-k6W3>MBeoO5Fgl3le3Zn3~Ms98aZ4BVy9@j zBoDI+HPWn>w_6sfuA3AggMqOi2pv3_@C#fjzkuKEhYvo$QCHwomZtE1ap#@;=oC@Y zYg}W0egT^TbrXzQF)!8T1Is|%wOR1b1~i2!Oyq;t)iB z$2vrgv-_Kdzyb3RmKs_JWgPxjpR=EzpII&Mt^x}iig^OJyRPdsaR>tz3X01R3%pBq zgOfSQQ93ceM3oSCbl`3tSXx@H{Pxd`#fFGfAwhRpjn=N+<%%GabXQWOjBWXc;VlzO zj$(Wknu+^IXgA_Q=bQ!bTIIVo(oT5&N z0%Kd&RY?2A!!eBA67(7>6gncHbn2|{Y4Qi>ilc4~1%2n^0fvFT#T(K=Y%yBet#?Cx ze4NF%Rmc&xwZ~=S@X>hJ{}+`B3dqsTKSXSQ)Uj56mghJI)Tu}lBSdI62hP?e)L79BfsM;ME9ssYr+|B*h%3 zlSc{Q15m~C?$oj!zhT@ixWoSIakPSWkg{-c{FJy0#PlKygcmkfQjpVYLf??vc7r0y z()}x(qZkI)f$5S>63OFsyH)Mbo5n86cTpkXZoPFOh@C zOv}u)(Yd`k)esLdC5fl>P2}p%M43BFK~SD-BQ9NjX5eR3rF%UvJ^%MY<8fwNI8vDb zfa*<_>=Pzb*xa@%sn+3zyFSACv9kC^=p8$4BF{_jT9EK*fzcGVCsFG?JF^4 z8{LGN6i21&=Po&09K(Du!3_5SG<@BL-d#1&JT@wFNC7-HJOlqs%=!CJyQ6xvDWEde72b!S8=p8#ro)EBUv1JY7tj^-3XuRB#HF zgJ=~)Coo)1d}M=ykJxAS-u7gQ0Ck-8|N5FtQ~O8x8MF*PakZDe?Ry%@D+{-*bUR*4 z=%}f!!9$Zw6>q@7NP=3t2v;C;RvrJfH!#t^k=QhBkctup?_y#xl>k(kL(5MG7P65w zsqX_iDDJHEanuMUdx_N*9R-;{Ro6;IAa0SKQ!=n@qQij2LEiuA5qj{WA$Z7@k{et_ zkB_+Qzm9Is%jKj5motMpo;=D7>Y$*$DsM0Bn|9?M1v}$C{!zr|)IEBB{*vsuWI+dY zs_Ooco-1sG7xT5bAS{EIvX$687Zjp08F6^_$r>ehdPc^3h@~}xazGi-kVojaP^Z7H z=L%(6;#Ba>q*m*wb+oeHFE1trwh?SuDn`g^TIdeCon4ro(0A`uOErK8u6TmTRrlf~ zUocl6NOv>GM;@IL6O+U@dYxYc$Li-chKlQitWhgss9c5PNEJp<{Ss7;WjUZs(6a!A zh_lmeS@dEjjT!8NX{`>}kKb=g#)7O~<7YLKu~U-Vxmz0@?gCD;W}1O|un;vc!TFFH z*}X!XK0C>s5LqpC6F)MG#(%iaF7PLD{@TqkhFjR{@>e`UC|7y#xmd&_U%(jI3^I|( z^=3F!#bYufF8`x$sggoC=bGv-2Gl*9Y)*oLyWC$y#2RD7#>N6JBCp(S?ZS48Hu2n@ zpXwrzH8jm1I^{(8%ll^|V^XTH32WGeBk=Ict3u*VR~r%8s;XZ*ey@wmr&PP+w@>O_ zpecLXSqm#GZ7oZrNauy*((7m~`03aRwiTEvyw_h$Hiqv~bQz^)q{4Uvt<0T;q~hcZ z>X^T@4el`#Wb=WAKEcJ!)AyeUs-7>_v3i~>!VA1rg;0{{aBLjVCWBz(Xmu{Z z#l^)VN_X0Z0Re7&8Y}o3BDX&Ulz1jie*y!RHV|L}d#qs6)f<3MpNxulV!xQhz#ng+gkWRx z+|;SOG<--t4wb_WutxGS{~W4-W$Qov+8n3Jprt=Y4zIV)G<)TJVtQH{0v|3=MY8^T zE$;@!hz?11_wDr{%!;^QhbH}>PQ9;RtIB=sD3RjLpDjO*SxUpO`kSe&En4bimX~?r z5wkT1_M!oTFRuy*HI|Uqxz#7?`kA-w+9vx8slGrtym7xo0Wj{_2`b1d(<%+*-Kq(D zj^7sR>#jAYm14Xwcw}5dnXLwjt>>l0#ylb#4n~!q7B(mK3tXNvI3^=}&q(eJ8gX=y?L&|+F~1iw`)FFO^U_`afEk*= zSa&Llh`28szIx~s>4b*HYK?2h1F?GHqrT2Ua`CaOub=49V)5%KdqCa(q$w-I6tkgHv2Kr|khuA>R@`oHMTLdoi06AN@yr|4Lg>kD#uD?K z3$QVYU6oZ%YZL*`_4Ks z-vCH+579jgTtxiP#_@kJ@8~wAlbQq&wD;Z|jo*2v_wn%ddyAy#C6251qnT1_X|-Q% zN+}cM6~f?6_Ek*MOK!^B>iG^}9g29tOHSS%wms!)QZTauK?nk6p*Y?xK&FrA5Dzr* zyP6er>VLvH%+9y$nkw+iahccl`z8o*cRAM`&+3GqHO3_kVsma%l-glYfCvHl2AhvHA^>drqKEfIQDy!YOcE)t96qQe~D4dHzWNyIOjQ$4|T zuwn||*xK79TuX<3Zn1-Y{)GWpvK@~HB1j~vLrHs$;lGUCr=Gi6kbPMwd`2XY$tR-; z*n2oTMc|1cxif-{orjThM9p-)rs^P`5RRwY88w^^IxJ742yxt@&*H^*!QoO-@iL`l zcUZJR`EC`4vVO7roPf{XJiTZ&FWz0ah%;$1=B7W~q=Krj@w4kH9|Z@TI$XNn(VaVD zAxd25{KdDqo&WLwv~1uW#J&#TmpIZIJITPt#|EzBcHxIKJsbpbm6YNgq4U<5_(ZSR zy3gWxddMSC`2Q}1(}Eh6q#kvEqK7+*Z}i`Bk+}5C z7YI;?V)U6nzyfZ`B>Owr6bNh2uf1{@sMTyw!8C{16KK^Fbe!u?1p`lUgljlmXr2y8 z#Wk6J$Cn-0YrJ^mp2DIri7WqV9fLtAZ`$zTOzbGmY3h@Dx?zFMq0qt*TgE|q-?n21 zMb6>`qLug)a=;<=g+cSaW6g}4%e~9oX8qSm#!vk+6Y0XU$X3j20>)FK(UJ5TmhXv* zjK6OX+oTaILs8Z`YbFZRVMrB=eh(;Q+z69=BT?D0tum^IMOc>E-+R8*<~sg@eM{i{ zDz>#VSTP`BN36hm2m&3wTz6p;o)I(NXX1zAGw)BUxkI!KQahX$tRHigXU}_n??F0wc_Vyoww{e*HUwb_6 zL*3KJ0yPexF{KHsyIc)p@z%STkGB8*Lbhr0;|l3+&EEl0bej~4o1nKpxZr~WAO$>M zaQ`06MuiJMmJ)pv?Dvfc=kE|F3rB!{`BxKvmd##|eJ>b8zU92%>H&{AUt^zZ|Dm_C zg1^M3na~SLJ6|7>@5FtJ%n!GDU*{U!S7Ut_^-hn%4~^A*oePYH*5u<8DBlkHN*Eu$ ztpt5530J@)cVAq>@qA9~i4CFkn#kCl`R#r9OXp z=O|?>LM-FC4|so8o$l%MdjHE$MzeUW*!}q?@V9u|+=j=`&|P_>bT&2}+pPFp=8Q9* zgX>oWxEi`Kodicjepf?AdAL18Qo%6l8BFD+fxZco(I)8H4Kj+i#rFqh`dy9bPB;w8 zC==jd5hqNL67-w;}ZLIv5ax7`19Q+8g|WMzT=RNA$3DKSa?`DsJ`5 zi~zmvBHj^#g25{fxkZ)5oRvS`UwA)f6XJE{#!vv((Zus&21(%i1}dw4ddaB6O9z=# zG`Pou8Nqvf*m_nU+M4NCboi;qnb-qM8{;PQQTQuM?lEj`Fiwbn@<aoh}Aw5#Q|aLE3zM4DswGn%9Ze=a&O0J zg(I>NT!`vR&SSzzeV#Iif-{oj91+k94qeLb?HQ$*79;j3*R0-u;Na|h5U_SoDJAqP z%_)Q;QPPpPPudXdOeA(TV9&S6kRb7x`)#i>xk-;Lakm6x$0xMS^DDp?Ts6}3RJ10| zSTjOXegG79B64nSH4=kDZlPmMfiEpODH3Z2ij00?cJ~s|?ciazG@{wvBL{HkCkvh? zp3~`&`s&o3P6FC8!(gc6IwoJ;;yhgr!Vad*K%2qVKuGkVc_A(DuVuh~p)uj6Q6L_0 zi3aH+AFxOlN>5|wpTG+F2GtU9YU_q#nGz^_3>pYD!x z2D_;Y^eswQoYF+}T$GrJ)dA|{hm4C@vU$OciGIe*I=V;otAi{k#p$=KwRN=34T6C7 zAZe#XCb#{b{?)FBTepl~yPIGy++DwzoqX*JFs1=q@_ z+L#uwB_yXZ+B4TwPmkcfUHp$uqtv(uhj{C~ZI246c)i-cdLJq3fDw2EVWw|1mv=#^ zv*%Ly?cY{8~#7B@-kX_P`y@&ZlOG+ zD!mn#Jp%sQLLR0HJp~R1V|?~MDU8@>P_?nUP1<{Jjd29@Q9tVBXc^s|bP1}N-lEA| z15`OWzP9#9dg*MVD`f$KK9!?T(zBAX1D^mS=jb9op*#bCV1}#uSjZ$lrj|GCJL+H= z%aQv&dOm`aafWN}2fEp}H=O;!kSFcM1mfDUt8x*8o+7tHjuTm9|BiPdhz%5c7x_B+ zT%8TQw{o;{9Q3|1H0GtQk;9^73z-BxLxltk(-@q>etvfj)EKbFNPrL8^z@I+t*)_- z)1|4RT6M@~>a9~rAAVAfPad<4M{ZtaUp{T4r+x-Mx=5$W;gNI1;0^dbuOIprzvcFOZ9cC(3gqtZ7e2L0yrPx) z-u?$4Hy~R?HB&mq^Yo%C6K=jPDCl=a zZ?q$Q@p8m7>XoJS;pc=^xY_nD78?|VfU4yWAq7(ZWPTm)%N1CI(DBaIW( zYW{oA6Id&Xp~Elq;wo$g%|kqG3N&1)mk0HT-nqO?+w`3L`#%Sz-0#Ey%J1?fGkEUO zDg~d$& z@RlMP=-*}K$ot_JcqMV{kixf*ji}k&s-+l!T^2DIP8yu-6rWnv)Np@BZx6*0i!}Yv zwnm2WL=cyz-9(TXWJxONga{#T$qSATx09u84Z^9#f%NiE(JJZ^@Btq+2n)nm`S;?TV{U$icEDPSl<+RHl$L@VH3@$tkiNK49rAbf~rqKBLDktc1B6G=3h--x6rr*9C&Q# zP|Rs0u}zCXACTj&v|Esa`duNKhv#KW4;%4w(V(CswZnN|_3x=r$VloKDBW=ZrBnE`wvI(#$3SOLRolGW*J_^Q z1vBL%!=qO9L$XSJ1^!4>!?J40S1(P7PCOxIZb(=c-Iyi2s<(crkkhTwb4z8d2x13| zLeG(Ud_x0!d=lf7s4`7{08$Sa+KBFZmJ+QI`uvTs(r+A*rlUk%e|a>m4?Y#Uw_`?a zG!BMFKA1Q=ZY>ul$_z`6uzFswM$dA*(*Er;XWTnNO~BIk!wZCk-AxCyUXzy53$lOr zBTD|WuRcA$CHLWUR;YXrR&}pg$55(usMv(BRji8g~=Ts^=tRp89lU3o!|7VaQVwdkg@iow7bslosih3PmC;&&7= z9FUoJIN-`l2Tnz>xtm+vbyOB(KU?sAIwSR*RG@q8P&YuB{{$+{C=Eox|3$2Td; z;Tn2DS01GX6TeLzS8b0iMUEzZsD=u81q|Bck&6l=GEKfDT8Prp71A;}A=zoI6ZfKO zF4u@MVew)F(IrE(R}HX*5z;HjVTB1Q!lIXEqQw+VDz|zhC)k`FzDiY-sbF;()=ZZF z_qTQ1w}~Gsd7T1TOvS@^DvvEv^Vt49xa;;u{$!90LPV{7Kav)w(3Z11O;pYD)Xq7o z!uVhrv?y?MRSC#$M${dmtAem&eeod4$wqPHMxhfvU;Bf9?tfAhq;D@{WvNDfi#oGJ-20e3nTjxZ^Dbf#w_~Io&Wa}f zVXU3zk}e1>j7J?pm)KyVsB+7*hQWfNiYA=;_aBY`Wr}eefkvco{3cBNuROf2yzC9# z-3P-}s=Ee+UIIg?;5uT;N5v^wgiUz0)ZAA8=Ha;Ka?k}GzsEocbEbDp{R_q{>{4NUNdgWfh-yj~L-_>pgh}_fo>loff6JzQNMVHLX`fzo}h-NDlWFbO{t0pp^j{ zd1Xl^3ut9Urj!odAmRmuv*YTO0b;5xI|T%U0x}b3mF#VVi#U&s9#?0MGO9m*#o~5i zRovNgSD@%lw$CCk87+P!7$}lkzRi-(EuTuEO2kN)acr|uYNpiTvLJiU#>t7wRPS%( z>+y1&C0ze@M8$_S#7nU7iu`r)+&+QZ93|OaFiPGr48e4lOZ_OO>P)jGn;ppHBn zk!T#>Uyy2~x!b0dZ-4Q};$$_L(jPe(4V74x@hI4eH1iKen><_!PDwtwSQ@=C+}K%; zo(UTr+d=Y)Rw$3y!lrO)mj>`}jNmiY5IR<@T7ciAzMa3AIrI?X%MQ2i5tYCE#h(+! zx)o^=AmY-Fv^G;gs-q3^|0*{t>V{cfa+R(+s(Wpw&i*3}JM5u9Uh^=g{DFi_=Ydt& z43Y%q(gNil0RGW6SvIp!QIcl(i3pmn z*3Td*%K~{3a!P&1PLi-Yw!&I@&-TH1@4L@=7|O?{lMeha0D;4 ze{6j$!8o-8gJC|pmStrp|JBtvA6qeU*5=Y(+=GIF+=)j0-wG7G^}X=Dv#=|1C!55k z4ny;V4N|yq@qt;H>W!@5U-$-C4;uqdbSrQPOWHB|ECy@{rs~<+6*ksgn$tyOSeJOW zvgF$c>~Top0IRZ_|6nY+WQB|5=E8luE26i>y6G>ocTk<(WxhXxh4_XOp*^g*h?8QH z*iAGD?mQf3&r_X{dWWyRd_pfz#i&ulOim-kr)LLxP^gJx?qg0vgkAak^D*EFV@ESa z3>$@^dv>g51u!|j@(#xZAUWwZwz~>%NMi@AY#XK)fPZp_9j|7AN!k9%yMCqb1DHr{ zHhSWdj3l?vW;M7%W2$#uy92Lt5;U}gdO!NSi>Y&bnvUJ_?ndV}@A&A(Ht(s!r&FB@ z8VkRFB?&h#V+LGoz~PVGZ+hVC?29SUkS#nv`4HTM1 zQ<%W=w(Em4X3Yl%8d1E`N2ze}_#~JZGtD{1n;hclEx1p{xc^4x&|rI1hhNzA)`>_} zxJVL*VSNeIpiADYJ9vI@JW?JkNa0xJaE#6~0`%li(IVKPH2q0(ph}<&X;`}ba{c=O z1Rr)uDqCX{t@0(3boVaRt0+U-eO=iFI{$RRtAq!Y`#>F_jHxpa2}7=xkv2e9rp5YM zoX2CcS@^pfYo2#**A<={gp?O=;s zbGPB<@rp&~qHn*dV0CFV5JP1ucN*^h$n9-DBb4G-W_7>3|&!C zBN;B7&D+)Z8+rR+T_9Z04UF`awx$2hx zjeqZ~?v5&hgz=NNT&qnt$-osUbE&&3LchEtx-~7Rm=CutGEhzNl(!g-ma37JA;P}+ zXDh^$u5(YAZBa#uH6er~THr!70*dLPzWyb~lH`BYop3ssOdb_;{RW|KiT_f4R6=p^ zCF4{?xOCWskM{0mb1&tJhQz9@%vOnt;Jv|pC!e|VK#uoX+(63GB;-tAK$efRVCat+@m@6qOcoS)d%#&{0ZRSU}qyks7{(cGikf-|$89hYpIzf%e2h@v0Z zYt8J5w!nZoI`8g@m>Ij#36Sr|s8jsAd~qdXB-{vov@Y@_kbSg0xmPMCIhph>>Lo4y ztA>gSt*@L^%Epo??P*{8Occ|V9|tSm;%|qC*%fHLqpjL&C7nL{by6x7WEGzAf^W#= zN&-`<&aRB%=7~=8v#njM+I4aSAoWlcedwCo`TpBiGQ=C+6a}Uj3IG$9N*16LD}In? z{>uYD&$q1lqBQ6iZ@FlYFwBYJ@GqspD4GuFZy!yLYgtUx7ZAqNG$T;KP&EFT1TE_a z&8$A%2qNYO#!>A7j5Kpxl=I6s)2WWZr><{~^t$oD4knGc`KYC&I5yRzO{1wTx?QTd zVV?aawg<7P{9P_L1YD+s_|a8X@)QPg`^wsl66Tz(7}UiYbc(5;pE!?%66i5hN1>v; zjJlOzd_tV7KGgKv9EGhqsgkOlKz96{^lDGUn%szyTBtYh+k35!iSBH@Wjqei>Ps&t zATHpYz$F$Ucfo6mps=Ok?B!rkLvnLAMRn(5w8qiT3`AL92R7FDY-I=HmVPW@#WOqS z&W(hoym%z-E=j$M(vIo#c~kYM^LbiL4ee*K>$B43^Q1Hq1@-I&p<9kvG<31FiM5i= zgm30uQ?a5&(bLl~B#3P9dcsHO`V1&?R z(?>PMEUs4b0@#VQmQyswD%8qUOP1^7$p))XdI(mNKm3BvqYZ3qUBccT>?c)v$^O;< zA+X<&FJsm+Q9ESv;W{ep%D3Q8*4{9r5&c)2jG-wJQ=M6)8xY)s1{F@!>Z+hYX@WlUFJlCzWzliide1fbn7MoT*%C@eQ(_7g{3f$-uqzUV? z53nM&z5KwaA1NG z%6wdLbU|BuB%Kyd%xE&tv4e}+tx(o2n>B9Vxa z5|V%i%8-(mAOQ-WvY#LUmLHU|9YpZo5*?^+igTg@XkJQpqAZx-&y>!e^#9|5|GzWm z|Hl2lVdYW`lQcoUQ<{?$K~qvdlNCUDQk0TKK)+HFl6k=VV^U0Gs8X(z=|FE%@KSi7 KZsY#P`Tqd-EBE*S diff --git a/python/MXTuning.py b/python/MXTuning.py index 27c368b..c2984d3 100755 --- a/python/MXTuning.py +++ b/python/MXTuning.py @@ -76,6 +76,11 @@ class MXTuning(Tuning): num1=np.poly1d([mag1]) den1 = np.poly1d([T1**2,2*T1*d1,1]) + #reiner integrator: 30Hz=0dB -> k=30*2*pi=180 + #num1=np.poly1d([120*120]) + #den1 = np.poly1d([1,0,0]) + + #first resonance frequency f2=np.array([197,199]) d2=np.array([.02,.02])#daempfung @@ -101,9 +106,9 @@ class MXTuning(Tuning): num=num1*num2*numc#*num3 den=den1*den2*denc#*den3 mdl= signal.lti(num, den) #num denum - print num - print den - print mdl + print(num) + print(den) + print(mdl) d={'num':num.coeffs,'num1':num1.coeffs,'num2':num2.coeffs,'numc':numc.coeffs, 'den':den.coeffs,'den1':den1.coeffs,'den2':den2.coeffs,'denc':denc.coeffs} fn=os.path.join(base,'model%d.mat'%mot) @@ -175,9 +180,9 @@ class MXTuning(Tuning): num=num1*num2*num3*num4*num5*numc den=den1*den2*den3*den4*den5*denc mdl= signal.lti(num, den) #num denum - print num - print den - print mdl + print(num) + print(den) + print(mdl) d={'num':num.coeffs,'num1':num1.coeffs,'num2':num2.coeffs,'num3':num3.coeffs,'num4':num4.coeffs,'num5':num5.coeffs,'numc':numc.coeffs, 'den':den.coeffs,'den1':den1.coeffs,'den2':den2.coeffs,'den3':den3.coeffs,'den4':den4.coeffs,'den5':den5.coeffs,'denc':denc.coeffs} fn=os.path.join(base,'model%d.mat'%mot) @@ -196,6 +201,100 @@ class MXTuning(Tuning): # tp print see also: print(np.poly1d([1,2,3], variable='s')), print(np.poly1d([1,2,3], r=True, variable='s')) + def custom_chirp(self): + motor = 1 + amp, minFrq, maxFrq, tSec = (10, 10, 300, 30) + file='/tmp/gather.npz' + # if not os.path.isfile(f): tune.init_stage();plt.close('all') + # tune.bode_chirp(openloop=True, file=f, motor=mot, amp=amp, minFrq=minFrq, maxFrq=maxFrq, tSec=tSec) + prog = ''' + &0 //cout works only in coord 0 + open prog 999 + L11=0 + L10=0 + L12=0 + L13=0 + L14=0 + Gather.Enable=2 + while(L10<300005) + { + L12=10*sin(31.415926535897931*(pow(1.058324104020218,(L10*0.000199996614513))-1)/log(1.058324104020218)) + cout%d:(L12) + L10=L10+1 + } + Gather.Enable=0 + close + b999r + '''%motor + gpascii = self.comm.gpascii + gt = self.gather + print(gpascii.servo_period) + gt.set_phasemode(False) + address=("Motor[1].IqCmd", "Motor[1].ActPos",) + gt.set_address(*address) + #Gather.Enable=1 + gt.set_property(MaxSamples=300000, Period=1) + + # gt.enable(2) + gpascii.send_line(prog) + gpascii.sync() + + gt.wait_stopped() + self.data=data=gt.upload() + meta={'motor':motor,'date':time.asctime(),'minFrq':minFrq,'maxFrq':maxFrq,'tSec':tSec,'amp':amp,'address':address} + np.savez_compressed(file, data=data, meta=meta) + meta['file'] = file + self.bode_chirp_plot(data, meta,True) + pass + + + def bode_sine(self,openloop=True,motor=1,minFrq=1,maxFrq=20,numFrq=15,amp=10,file='/tmp/gather.npz'): + '''calculates phase and amplitude at different frequencies and + saves:#loads and plots the bode diagram''' + if False:# os.path.isfile(file): + f=np.load(file) + bode=f['bode'] + meta=f['meta'].item() + meta['file']=file + else: + gpascii=self.comm.gpascii + #motor 1 maximum: 13750 + #amp= percentage of maximum amplitude + fig = plt.figure() + ax = fig.add_subplot(1, 1, 1) + frqLst=np.logspace(np.log10(minFrq),np.log10(maxFrq),numFrq) + n=len(frqLst) + #frqLst=(10,15,20,25,30) + bode=np.ndarray((n,3)) + bode[:, 0]=frqLst + #for i in range(n): + for i in range(n-1,-1,-1): + frq=frqLst[i] + t=1 + rep=max(1,frq*t) + if openloop: + data=self.do_command('openloopsine',motor,amp,frq,rep,0) + else: + data=self.do_command('sinusoidal',motor,amp,frq,rep,0) + data=data[:,(1,2)] + gpascii.send_line('#1j=0') + time.sleep(1) + ax.clear() + avg=data.mean(0) + print(avg) + ax.plot(data[:, 0]-avg[0] , 'b-', label='input') + ax.plot(data[:, 1]-avg[1], 'g-', label='output') + #plt.pause(.05) + bode[i,1:]=self.phase_amp(frq, rep) + print('frq %g ampl %g phase %g'%tuple(bode[i,:])) + plt.show(block=False);plt.pause(.05) + + meta={'motor':motor,'date':time.asctime()} + np.savez_compressed(file, bode=bode, meta=meta) + meta['file']=file + self.bode_sine_plot(bode, meta) + + def bode(mdl): w,mag,phase = signal.bode(mdl,1000) f=w/(2*np.pi) @@ -214,6 +313,16 @@ def bode(mdl): if __name__=='__main__': from argparse import ArgumentParser,RawDescriptionHelpFormatter + import logging + logger = logging.getLogger(__name__) + logger = logging.getLogger('pbtools.misc.pp_comm') + logger.setLevel(logging.DEBUG) + logging.basicConfig(format=('%(asctime)s %(name)-12s ' + '%(levelname)-8s %(message)s'), + datefmt='%m-%d %H:%M', + ) + + def parse_args(): 'main command line interpreter function' #usage: gpasciiCommunicator.py --host=PPMACZT84 myPowerBRICK.cfg @@ -315,7 +424,13 @@ Examples:'''+''.join(map(lambda s:cmd+s, exampleCmd))+'\n ' tune.bode_sine(openloop=False, file=fn) if os.path.basename(fn).startswith('chirp_cl_mot'): tune.bode_chirp(openloop=False, file=fn) - print 'done' + print('done') + elif mode==7: #further tests + tune.init_stage(); + plt.close('all') + tune.bode_sine() + tune.custom_chirp() + plt.show() #------------------ Main Code ---------------------------------- #ssh_test() diff --git a/python/helicalscan.py b/python/helicalscan.py index 3e681e3..c1d5a7a 100755 --- a/python/helicalscan.py +++ b/python/helicalscan.py @@ -151,19 +151,19 @@ class HelicalScan: param = self.param cx, cz, w, fy, = (0.2,0.3,0.1,0.4) #cx, cz, w, fy, = (10.,20,3.,40) - print 'input : cx:%.6g cz:%.6g w:%.6g fy:%.6g' % (cx,cz,w/d2r*1000.,fy) + print('input : cx:%.6g cz:%.6g w:%.6g fy:%.6g' % (cx,cz,w/d2r*1000.,fy)) (dx,dz,w,y) = self.fwd_transform(cx,cz,w,fy) - print 'fwd_trf: dx:%.6g dz:%.6g w:%.6g fy:%.6g' % (dx,dz,w/d2r*1000.,y) + print('fwd_trf: dx:%.6g dz:%.6g w:%.6g fy:%.6g' % (dx,dz,w/d2r*1000.,y)) (cx,cz,w,fy) = self.inv_transform(dx,dz,w,y) - print 'inv_trf: cx:%.6g cz:%.6g w:%.6g fy:%.6g' % (cx,cz,w/d2r*1000.,fy) + print('inv_trf: cx:%.6g cz:%.6g w:%.6g fy:%.6g' % (cx,cz,w/d2r*1000.,fy)) dx, dz, w, y, = (0.2,0.3,0.1,0.4) #dx, dz, w, y, = (10.,20,3.,40) - print 'input : dx:%.6g dz:%.6g w:%.6g fy:%.6g' % (dx,dz,w/d2r*1000.,y) + print('input : dx:%.6g dz:%.6g w:%.6g fy:%.6g' % (dx,dz,w/d2r*1000.,y)) (cx,cz,w,fy) = self.inv_transform(dx,dz,w,y) - print 'inv_trf: cx:%.6g cz:%.6g w:%.6g fy:%.6g' % (cx,cz,w/d2r*1000.,fy) + print('inv_trf: cx:%.6g cz:%.6g w:%.6g fy:%.6g' % (cx,cz,w/d2r*1000.,fy)) (dx,dz,w,y) = self.fwd_transform(cx,cz,w,fy) - print 'fwd_trf: dx:%.6g dz:%.6g w:%.6g fy:%.6g' % (dx,dz,w/d2r*1000.,y) + print('fwd_trf: dx:%.6g dz:%.6g w:%.6g fy:%.6g' % (dx,dz,w/d2r*1000.,y)) @@ -306,7 +306,7 @@ class HelicalScan: p[i,0]=x_i+r_i*np.sin(phi_i) # x= x_i+r_i*cos(phi_i+w)+cx p[i,1]=y_i # y= y_i p[i,2]=z_i+r_i*np.cos(phi_i) # z= z_i+r_i*sin(phi_i*w) - print p + print(p) ofs=(p[1]+p[0])/2. # = center of the cristal m=Trf.trans(*ofs); self.hOrig=self.pltOrig(m) @@ -327,7 +327,7 @@ class HelicalScan: p[i, 0] = x_i + r_i * np.cos(phi_i) # x= x_i+r_i*cos(phi_i+w)+cx p[i, 1] = y_i # y= y_i p[i, 2] = z_i + r_i * np.sin(phi_i) # z= z_i+r_i*sin(phi_i*w) - print p + print(p) ofs = (p[1] + p[0]) / 2. # = center of the cristal m = Trf.trans(cx,fy,cz) @@ -375,7 +375,7 @@ class HelicalScan: p[i,1]=y_i # y= y_i p[i,2]=z_i+r_i*np.cos(phi_i) # z= z_i+r_i*sin(phi_i*w) ofs=(p[1]+p[0])/2. # = center of the cristal - print 'p, ofs',p,ofs + print('p, ofs',p,ofs) m=Trf.trans(0,0,0); self.hOrig=self.pltOrig(m) hCrist,pt=self.pltCrist(cx=-ofs[0],fy=-ofs[1],cz=-ofs[2]) @@ -458,7 +458,7 @@ class HelicalScan: for k,v in fh.iteritems(): s+=' '+k+': '+str(v.dtype)+' '+str(v.shape)+'\n' setattr(self,k,v) - print s + print(s) def fwd_transform(self,cx,cz,w,fy): #cx,cy: coarse stage @@ -530,7 +530,7 @@ class HelicalScan: param[i, 1] = y param[i, 2:] = HelicalScan.meas_rot_ctr(x) # (bias,ampl,phase) (bias, ampl, phase) = param[i][2:] - print param + print(param) def calcParam(self,x=((-241.,96.,-53.),(-162.,-293.,246.)), @@ -569,15 +569,15 @@ class HelicalScan: x_ = ampl * np.cos(w + phase) + bias print(x_) (dx,dz,w,y_) = (0,0,0,y[0]) - print 'input : dx:%.6g dz:%.6g w:%.6g fy:%.6g' % (dx,dz,w/d2r*1000.,y_) + print('input : dx:%.6g dz:%.6g w:%.6g fy:%.6g' % (dx,dz,w/d2r*1000.,y_)) (cx,cz,w,fy) = self.inv_transform(dx,dz,w,y_) - print 'inv_trf: cx:%.6g cz:%.6g w:%.6g fy:%.6g' % (cx,cz,w/d2r*1000.,fy) + print('inv_trf: cx:%.6g cz:%.6g w:%.6g fy:%.6g' % (cx,cz,w/d2r*1000.,fy)) (dx, dz, w, y_) = (0,0,0,y[1]) - print 'input : dx:%.6g dz:%.6g w:%.6g fy:%.6g' % (dx, dz, w / d2r * 1000., y_) + print('input : dx:%.6g dz:%.6g w:%.6g fy:%.6g' % (dx, dz, w / d2r * 1000., y_)) (cx, cz, w, fy) = self.inv_transform(dx, dz, w, y_) - print 'inv_trf: cx:%.6g cz:%.6g w:%.6g fy:%.6g' % (cx, cz, w / d2r * 1000., fy) + print('inv_trf: cx:%.6g cz:%.6g w:%.6g fy:%.6g' % (cx, cz, w / d2r * 1000., fy)) - print param + print(param) def pltOrig(self,m,h=None): diff --git a/python/shapepath.py b/python/shapepath.py index 1dfc3f7..5c624bb 100755 --- a/python/shapepath.py +++ b/python/shapepath.py @@ -56,7 +56,6 @@ import matplotlib as mpl import matplotlib.pyplot as plt import subprocess as sprc import telnetlib -from utilities import * class ShapePath: def __init__(self,args): diff --git a/python/utilities.py b/python/utilities.py new file mode 100755 index 0000000..77f83ea --- /dev/null +++ b/python/utilities.py @@ -0,0 +1,91 @@ +#!/usr/bin/env python +#*-----------------------------------------------------------------------* +#| | +#| Copyright (c) 2016 by Paul Scherrer Institute (http://www.psi.ch) | +#| | +#| Author Thierry Zamofing (thierry.zamofing@psi.ch) | +#*-----------------------------------------------------------------------* +''' +utilities classes +''' +import json +import numpy as np +import time,os + + +class dotdict(dict): + """dot.notation access to dictionary attributes""" + def __init__(self,arg=None,**kwargs): + if arg!=None: + self.__fill__(arg) + self.__fill__(kwargs) + + def __fill__(self,kw): + for k,v in kw.iteritems(): + if type(v)==dict: + self[k]=dotdict(v) + else: + self[k]=v + if type(v)==list: + for i,w in enumerate(v): + if type(w)==dict: + v[i]=dotdict(w) + pass + + def __dir__(self): + l=dir(object) + #l.extend(self.keys()) + l.extend(map(str,self.keys())) + return l + + def __getattr__(self, attr): + #return self.get(attr) + try: + return self[attr] + except KeyError as e: + raise AttributeError("%r instance has no attribute %r" % (self.__class__, attr)) + + def __repr__(self): + return '<' + dict.__repr__(self)[1:-1] + '>' + + def PrettyPrint(self,indent=0): + for k,v in self.iteritems(): + if type(v)==dotdict: + print(' '*indent,str(k)+':') + v.PrettyPrint(indent+2) + else: + print(' '*indent+str(k)+'\t'+str(v)) + + __setattr__= dict.__setitem__ + __delattr__= dict.__delitem__ + #__getattr__= dict.__getattr__ + + +def ConvUtf8(s): + 'convert unicoded json object to ASCII encoded' + #http://stackoverflow.com/questions/956867/how-to-get-string-objects-instead-of-unicode-ones-from-json-in-python + if isinstance(s, dict): + return {ConvUtf8(key): ConvUtf8(value) for key, value in s.items()} + elif isinstance(s, list): + return [ConvUtf8(element) for element in s] + elif isinstance(s, str): + return s.encode('utf-8') + else: + return s + +class GpasciiCommunicator(): + '''Communicates with the Delta Tau gpascii programm + ''' + gpascii_ack="\x06\r\n" + gpascii_inp='Input\r\n' + + def connect(self, host, username='root', password='deltatau',prompt='ppmac# '): + p=telnetlib.Telnet(host) + print(p.read_until('login: ')) + p.write(username+'\n') + print(p.read_until('Password: ')) + p.write(password+'\n') + print(p.read_until(prompt)) # command prompt + p.write('gpascii -2\n') # execute gpascii command + print(p.read_until(self.gpascii_inp)) + return p