From ce3385e65f67511a9fec2f9b56641af2c6939e19 Mon Sep 17 00:00:00 2001
From: sayanmitracode <sayan.mitra@gmail.com>
Date: Mon, 4 Jul 2022 12:18:32 -0500
Subject: [PATCH] made F16 stanalone
---
demo/F16/README.md | 3 +-
demo/F16/aerobench/highlevel/autopilot.py | 97 +++++
.../F16/aerobench/highlevel/controlled_f16.py | 60 +++
demo/F16/aerobench/lowlevel/adc.py | 37 ++
demo/F16/aerobench/lowlevel/cl.py | 53 +++
demo/F16/aerobench/lowlevel/clf16.py | 85 ++++
demo/F16/aerobench/lowlevel/cm.py | 50 +++
demo/F16/aerobench/lowlevel/cn.py | 54 +++
demo/F16/aerobench/lowlevel/conf16.py | 58 +++
demo/F16/aerobench/lowlevel/cx.py | 50 +++
demo/F16/aerobench/lowlevel/cy.py | 9 +
demo/F16/aerobench/lowlevel/cz.py | 31 ++
demo/F16/aerobench/lowlevel/dampp.py | 42 ++
demo/F16/aerobench/lowlevel/dlda.py | 50 +++
demo/F16/aerobench/lowlevel/dldr.py | 52 +++
demo/F16/aerobench/lowlevel/dnda.py | 51 +++
demo/F16/aerobench/lowlevel/dndr.py | 49 +++
.../lowlevel/low_level_controller.py | 109 +++++
demo/F16/aerobench/lowlevel/morellif16.py | 189 ++++++++
demo/F16/aerobench/lowlevel/pdot.py | 29 ++
demo/F16/aerobench/lowlevel/rtau.py | 18 +
demo/F16/aerobench/lowlevel/subf16_model.py | 201 +++++++++
demo/F16/aerobench/lowlevel/tgear.py | 14 +
demo/F16/aerobench/lowlevel/thrust.py | 76 ++++
demo/F16/aerobench/util.py | 287 +++++++++++++
demo/F16/aerobench/visualize/anim3d.py | 404 ++++++++++++++++++
.../visualize/bak_matplotlib.mlpstyle | 8 +
demo/F16/aerobench/visualize/f-16.mat | Bin 0 -> 32222 bytes
demo/F16/aerobench/visualize/plot.py | 294 +++++++++++++
29 files changed, 2459 insertions(+), 1 deletion(-)
create mode 100644 demo/F16/aerobench/highlevel/autopilot.py
create mode 100644 demo/F16/aerobench/highlevel/controlled_f16.py
create mode 100644 demo/F16/aerobench/lowlevel/adc.py
create mode 100644 demo/F16/aerobench/lowlevel/cl.py
create mode 100644 demo/F16/aerobench/lowlevel/clf16.py
create mode 100644 demo/F16/aerobench/lowlevel/cm.py
create mode 100644 demo/F16/aerobench/lowlevel/cn.py
create mode 100644 demo/F16/aerobench/lowlevel/conf16.py
create mode 100644 demo/F16/aerobench/lowlevel/cx.py
create mode 100644 demo/F16/aerobench/lowlevel/cy.py
create mode 100644 demo/F16/aerobench/lowlevel/cz.py
create mode 100644 demo/F16/aerobench/lowlevel/dampp.py
create mode 100644 demo/F16/aerobench/lowlevel/dlda.py
create mode 100644 demo/F16/aerobench/lowlevel/dldr.py
create mode 100644 demo/F16/aerobench/lowlevel/dnda.py
create mode 100644 demo/F16/aerobench/lowlevel/dndr.py
create mode 100644 demo/F16/aerobench/lowlevel/low_level_controller.py
create mode 100644 demo/F16/aerobench/lowlevel/morellif16.py
create mode 100644 demo/F16/aerobench/lowlevel/pdot.py
create mode 100644 demo/F16/aerobench/lowlevel/rtau.py
create mode 100644 demo/F16/aerobench/lowlevel/subf16_model.py
create mode 100644 demo/F16/aerobench/lowlevel/tgear.py
create mode 100644 demo/F16/aerobench/lowlevel/thrust.py
create mode 100644 demo/F16/aerobench/util.py
create mode 100644 demo/F16/aerobench/visualize/anim3d.py
create mode 100644 demo/F16/aerobench/visualize/bak_matplotlib.mlpstyle
create mode 100644 demo/F16/aerobench/visualize/f-16.mat
create mode 100644 demo/F16/aerobench/visualize/plot.py
diff --git a/demo/F16/README.md b/demo/F16/README.md
index 654c1430..d1cc352b 100644
--- a/demo/F16/README.md
+++ b/demo/F16/README.md
@@ -1,3 +1,4 @@
# Adaptation of F16 model
-This is an adaptation of the Python F16 model from (Stanley Bak)[https://github.com/stanleybak/AeroBenchVVPython].
+This is an adaptation of the Python F16 model from [Stanley Bak](https://github.com/stanleybak/AeroBenchVVPython).
+All the files are copied to make this example a standalone executable.
diff --git a/demo/F16/aerobench/highlevel/autopilot.py b/demo/F16/aerobench/highlevel/autopilot.py
new file mode 100644
index 00000000..d369cdc8
--- /dev/null
+++ b/demo/F16/aerobench/highlevel/autopilot.py
@@ -0,0 +1,97 @@
+'''
+Stanley Bak
+Autopilot State-Machine Logic
+
+There is a high-level advance_discrete_state() function, which checks if we should change the current discrete state,
+and a get_u_ref(f16_state) function, which gets the reference inputs at the current discrete state.
+'''
+
+import abc
+from math import pi
+
+import numpy as np
+from numpy import deg2rad
+
+from aerobench.lowlevel.low_level_controller import LowLevelController
+from aerobench.util import Freezable
+
+class Autopilot(Freezable):
+ '''A container object for the hybrid automaton logic for a particular autopilot instance'''
+
+ def __init__(self, init_mode, llc=None):
+
+ assert isinstance(init_mode, str), 'init_mode should be a string'
+
+ if llc is None:
+ # use default
+ llc = LowLevelController()
+
+ self.llc = llc
+ self.xequil = llc.xequil
+ self.uequil = llc.uequil
+
+ self.mode = init_mode # discrete state, this should be overwritten by subclasses
+
+ self.freeze_attrs()
+
+ def advance_discrete_mode(self, t, x_f16):
+ '''
+ advance the discrete mode based on the current aircraft state. Returns True iff the discrete mode
+ has changed. It's also suggested to update self.mode to the current mode name.
+ '''
+
+ return False
+
+ def is_finished(self, t, x_f16):
+ '''
+ returns True if the simulation should stop (for example, after maneuver completes)
+
+ this is called after advance_discrete_state
+ '''
+
+ return False
+
+ @abc.abstractmethod
+ def get_u_ref(self, t, x_f16):
+ '''
+ for the current discrete state, get the reference inputs signals. Override this one
+ in subclasses.
+
+ returns four values per aircraft: Nz, ps, Ny_r, throttle
+ '''
+
+ return
+
+ def get_checked_u_ref(self, t, x_f16):
+ '''
+ for the current discrete state, get the reference inputs signals and check them against ctrl limits
+ '''
+
+ rv = np.array(self.get_u_ref(t, x_f16), dtype=float)
+
+ assert rv.size % 4 == 0, "get_u_ref should return Nz, ps, Ny_r, throttle for each aircraft"
+
+ for i in range(rv.size //4):
+ Nz, _ps, _Ny_r, _throttle = rv[4*i:4*(i+1)]
+
+ l, u = self.llc.ctrlLimits.NzMin, self.llc.ctrlLimits.NzMax
+ assert l <= Nz <= u, f"autopilot commanded invalid Nz ({Nz}). Not in range [{l}, {u}]"
+
+ return rv
+
+class FixedSpeedAutopilot(Autopilot):
+ '''Simple Autopilot that gives a fixed speed command using proportional control'''
+
+ def __init__(self, setpoint, p_gain):
+ self.setpoint = setpoint
+ self.p_gain = p_gain
+
+ init_mode = 'tracking speed'
+ Autopilot.__init__(self, init_mode)
+
+ def get_u_ref(self, t, x_f16):
+ '''for the current discrete state, get the reference inputs signals'''
+
+ x_dif = self.setpoint - x_f16[0]
+
+ return 0, 0, 0, self.p_gain * x_dif
diff --git a/demo/F16/aerobench/highlevel/controlled_f16.py b/demo/F16/aerobench/highlevel/controlled_f16.py
new file mode 100644
index 00000000..c7ee1159
--- /dev/null
+++ b/demo/F16/aerobench/highlevel/controlled_f16.py
@@ -0,0 +1,60 @@
+'''
+Stanley Bak
+Python Version of F-16 GCAS
+ODE derivative code (controlled F16)
+'''
+
+from math import sin, cos
+
+import numpy as np
+from numpy import deg2rad
+
+from aerobench.lowlevel.subf16_model import subf16_model
+from aerobench.lowlevel.low_level_controller import LowLevelController
+
+def controlled_f16(t, x_f16, u_ref, llc, f16_model='morelli', v2_integrators=False):
+ 'returns the LQR-controlled F-16 state derivatives and more'
+
+ assert isinstance(x_f16, np.ndarray)
+ assert isinstance(llc, LowLevelController)
+ assert u_ref.size == 4
+
+ assert f16_model in ['stevens', 'morelli'], 'Unknown F16_model: {}'.format(f16_model)
+
+ x_ctrl, u_deg = llc.get_u_deg(u_ref, x_f16)
+
+ # Note: Control vector (u) for subF16 is in units of degrees
+ xd_model, Nz, Ny, _, _ = subf16_model(x_f16[0:13], u_deg, f16_model)
+
+ if v2_integrators:
+ # integrators from matlab v2 model
+ ps = xd_model[6] * cos(xd_model[1]) + xd_model[8] * sin(xd_model[1])
+
+ Ny_r = Ny + xd_model[8]
+ else:
+ # Nonlinear (Actual): ps = p * cos(alpha) + r * sin(alpha)
+ ps = x_ctrl[4] * cos(x_ctrl[0]) + x_ctrl[5] * sin(x_ctrl[0])
+
+ # Calculate (side force + yaw rate) term
+ Ny_r = Ny + x_ctrl[5]
+
+ xd = np.zeros((x_f16.shape[0],))
+ xd[:len(xd_model)] = xd_model
+
+ # integrators from low-level controller
+ start = len(xd_model)
+ end = start + llc.get_num_integrators()
+ int_der = llc.get_integrator_derivatives(t, x_f16, u_ref, Nz, ps, Ny_r)
+ xd[start:end] = int_der
+
+ # Convert all degree values to radians for output
+ u_rad = np.zeros((7,)) # throt, ele, ail, rud, Nz_ref, ps_ref, Ny_r_ref
+
+ u_rad[0] = u_deg[0] # throttle
+
+ for i in range(1, 4):
+ u_rad[i] = deg2rad(u_deg[i])
+
+ u_rad[4:7] = u_ref[0:3] # inner-loop commands are 4-7
+
+ return xd, u_rad, Nz, ps, Ny_r
diff --git a/demo/F16/aerobench/lowlevel/adc.py b/demo/F16/aerobench/lowlevel/adc.py
new file mode 100644
index 00000000..a70390eb
--- /dev/null
+++ b/demo/F16/aerobench/lowlevel/adc.py
@@ -0,0 +1,37 @@
+'''
+Stanley Bak
+adc.py for F-16 model
+'''
+
+from math import sqrt
+
+def adc(vt, alt):
+ '''converts velocity (vt) and altitude (alt) to mach number (amach) and dynamic pressure (qbar)
+
+ See pages 63-65 of Stevens & Lewis, "Aircraft Control and Simulation", 2nd edition
+ '''
+
+ # vt = freestream air speed
+
+ ro = 2.377e-3
+ tfac = 1 - .703e-5 * alt
+
+ if alt >= 35000: # in stratosphere
+ t = 390
+ else:
+ t = 519 * tfac # 3 rankine per atmosphere (3 rankine per 1000 ft)
+
+ # rho = freestream mass density
+ rho = ro * tfac**4.14
+
+ # a = speed of sound at the ambient conditions
+ # speed of sound in a fluid is the sqrt of the quotient of the modulus of elasticity over the mass density
+ a = sqrt(1.4 * 1716.3 * t)
+
+ # amach = mach number
+ amach = vt / a
+
+ # qbar = dynamic pressure
+ qbar = .5 * rho * vt * vt
+
+ return amach, qbar
diff --git a/demo/F16/aerobench/lowlevel/cl.py b/demo/F16/aerobench/lowlevel/cl.py
new file mode 100644
index 00000000..6bff3591
--- /dev/null
+++ b/demo/F16/aerobench/lowlevel/cl.py
@@ -0,0 +1,53 @@
+'''
+Stanley Bak
+F-16 GCAS in Python
+cl function
+'''
+
+import numpy as np
+
+from aerobench.util import fix, sign
+
+def cl(alpha, beta):
+ 'cl function'
+
+ a = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], \
+ [-.001, -.004, -.008, -.012, -.016, -.022, -.022, -.021, -.015, -.008, -.013, -.015], \
+ [-.003, -.009, -.017, -.024, -.030, -.041, -.045, -.040, -.016, -.002, -.010, -.019], \
+ [-.001, -.010, -.020, -.030, -.039, -.054, -.057, -.054, -.023, -.006, -.014, -.027], \
+ [.000, -.010, -.022, -.034, -.047, -.060, -.069, -.067, -.033, -.036, -.035, -.035], \
+ [.007, -.010, -.023, -.034, -.049, -.063, -.081, -.079, -.060, -.058, -.062, -.059], \
+ [.009, -.011, -.023, -.037, -.050, -.068, -.089, -.088, -.091, -.076, -.077, -.076]], dtype=float).T
+
+ s = .2 * alpha
+ k = fix(s)
+
+ if k <= -2:
+ k = -1
+
+ if k >= 9:
+ k = 8
+
+ da = s - k
+ l = k + fix(1.1 * sign(da))
+ s = .2 * abs(beta)
+ m = fix(s)
+ if m == 0:
+ m = 1
+
+ if m >= 6:
+ m = 5
+
+ db = s - m
+ n = m + fix(1.1 * sign(db))
+ l = l + 3
+ k = k + 3
+ m = m + 1
+ n = n + 1
+ t = a[k-1, m-1]
+ u = a[k-1, n-1]
+ v = t + abs(da) * (a[l-1, m-1] - t)
+ w = u + abs(da) * (a[l-1, n-1] - u)
+ dum = v + (w - v) * abs(db)
+
+ return dum * sign(beta)
diff --git a/demo/F16/aerobench/lowlevel/clf16.py b/demo/F16/aerobench/lowlevel/clf16.py
new file mode 100644
index 00000000..a09addfe
--- /dev/null
+++ b/demo/F16/aerobench/lowlevel/clf16.py
@@ -0,0 +1,85 @@
+'''
+Stanley Bak
+clf16.py for F-16 model
+
+This is the objective function for finding the trim condition of the initial states
+'''
+
+from math import asin, sin
+
+from tgear import tgear
+from conf16 import conf16
+from subf16_model import subf16_model
+
+def clf16(s, x, u, const, model='stevens', adjust_cy=True):
+ '''
+ objective function of the optimization to find the trim conditions
+
+ x and u get modified in-place
+ returns the cost
+ '''
+
+ _, singam, _, _, tr, _, _, _, thetadot, _, _, orient = const
+ gamm = asin(singam)
+
+ if len(s) == 3:
+ u[0] = s[0]
+ u[1] = s[1]
+ x[1] = s[2]
+ else:
+ u[0] = s[0]
+ u[1] = s[1]
+ u[2] = s[2]
+ u[3] = s[3]
+ x[1] = s[4]
+ x[3] = s[5]
+ x[4] = s[6]
+
+ #
+ # Get the current power and constraints
+ #
+ x[12] = tgear(u[0])
+ [x, u] = conf16(x, u, const)
+
+ # we just want the derivative
+ subf16 = lambda x, u: subf16_model(x, u, model, adjust_cy)[0]
+
+ xd = subf16(x, u)
+
+ #
+ # Steady Level flight
+ #
+ if orient == 1:
+ r = 100.0*(xd[0]**2 + xd[1]**2 + xd[2]**2 + xd[6]**2 + xd[7]**2 + xd[8]**2)
+
+ #
+ # Steady Climb
+ #
+ if orient == 2:
+ r = 500.0*(xd[11]-x[0]*sin(gamm))**2 + xd[0]**2 + 100.0*(xd[1]**2 + xd[2]**2) + \
+ 10.0*(xd[6]**2 + xd[7]**2 + xd[8]**2)
+
+
+
+ #
+ # Coord Turn
+ #
+ if orient == 3:
+ r = xd[0]*xd[0] + 100.0 * (xd[1] * xd[1] + xd[2]*xd[2] + xd[11]*xd[11]) + 10.0*(xd[6]*xd[6] + \
+ xd[7]*xd[7]+xd[8]*xd[8]) + 500.0*(xd[5] - tr)**2
+
+ #
+ # Pitch Pull Up
+ #
+
+
+ if orient == 4:
+ r = 500.0*(xd[4]-thetadot)**2 + xd[0]**2 + 100.0*(xd[1]**2 + xd[2]**2) + 10.0*(xd[6]**2 + xd[7]**2 + xd[8]**2)
+
+ #
+ # Scale r if it is less than 1
+ #
+ if r < 1.0:
+ r = r**0.5
+
+ return r
diff --git a/demo/F16/aerobench/lowlevel/cm.py b/demo/F16/aerobench/lowlevel/cm.py
new file mode 100644
index 00000000..4c00a175
--- /dev/null
+++ b/demo/F16/aerobench/lowlevel/cm.py
@@ -0,0 +1,50 @@
+'''
+Stanley Bak
+F-16 GCAS Python
+'''
+
+import numpy as np
+from aerobench.util import fix, sign
+
+def cm(alpha, el):
+ 'cm function'
+
+ a = np.array([[.205, .168, .186, .196, .213, .251, .245, .238, .252, .231, .198, .192], \
+ [.081, .077, .107, .110, .110, .141, .127, .119, .133, .108, .081, .093], \
+ [-.046, -.020, -.009, -.005, -.006, .010, .006, -.001, .014, .000, -.013, .032], \
+ [-.174, -.145, -.121, -.127, -.129, -.102, -.097, -.113, -.087, -.084, -.069, -.006], \
+ [-.259, -.202, -.184, -.193, -.199, -.150, -.160, -.167, -.104, -.076, -.041, -.005]], dtype=float).T
+
+ s = .2 * alpha
+ k = fix(s)
+
+ if k <= -2:
+ k = -1
+
+ if k >= 9:
+ k = 8
+
+ da = s - k
+ l = k + fix(1.1 * sign(da))
+ s = el / 12
+ m = fix(s)
+
+ if m <= -2:
+ m = -1
+
+ if m >= 2:
+ m = 1
+
+ de = s - m
+ n = m + fix(1.1 * sign(de))
+ k = k + 3
+ l = l + 3
+ m = m + 3
+ n = n + 3
+ t = a[k-1, m-1]
+ u = a[k-1, n-1]
+ v = t + abs(da) * (a[l-1, m-1] - t)
+ w = u + abs(da) * (a[l-1, n-1] - u)
+
+ return v + (w - v) * abs(de)
+
diff --git a/demo/F16/aerobench/lowlevel/cn.py b/demo/F16/aerobench/lowlevel/cn.py
new file mode 100644
index 00000000..b06e3729
--- /dev/null
+++ b/demo/F16/aerobench/lowlevel/cn.py
@@ -0,0 +1,54 @@
+'''
+Stanley Bak
+F16 GCAS in Python
+cn function
+'''
+
+import numpy as np
+from aerobench.util import fix, sign
+
+def cn(alpha, beta):
+ 'cn function'
+
+ a = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], \
+ [.018, .019, .018, .019, .019, .018, .013, .007, .004, -.014, -.017, -.033], \
+ [.038, .042, .042, .042, .043, .039, .030, .017, .004, -.035, -.047, -.057], \
+ [.056, .057, .059, .058, .058, .053, .032, .012, .002, -.046, -.071, -.073], \
+ [.064, .077, .076, .074, .073, .057, .029, .007, .012, -.034, -.065, -.041], \
+ [.074, .086, .093, .089, .080, .062, .049, .022, .028, -.012, -.002, -.013], \
+ [.079, .090, .106, .106, .096, .080, .068, .030, .064, .015, .011, -.001]], dtype=float).T
+
+ s = .2 * alpha
+ k = fix(s)
+
+ if k <= -2:
+ k = -1
+
+ if k >= 9:
+ k = 8
+
+ da = s - k
+ l = k + fix(1.1 * sign(da))
+ s = .2 * abs(beta)
+ m = fix(s)
+
+ if m == 0:
+ m = 1
+
+ if m >= 6:
+ m = 5
+
+ db = s - m
+ n = m + fix(1.1 * sign(db))
+ l = l + 3
+ k = k + 3
+ m = m + 1
+ n = n + 1
+ t = a[k-1, m-1]
+ u = a[k-1, n-1]
+
+ v = t + abs(da) * (a[l-1, m-1] - t)
+ w = u + abs(da) * (a[l-1, n-1] - u)
+ dum = v + (w - v) * abs(db)
+
+ return dum * sign(beta)
diff --git a/demo/F16/aerobench/lowlevel/conf16.py b/demo/F16/aerobench/lowlevel/conf16.py
new file mode 100644
index 00000000..f322a6de
--- /dev/null
+++ b/demo/F16/aerobench/lowlevel/conf16.py
@@ -0,0 +1,58 @@
+'''
+Stanley Bak
+Python F-16
+
+Apply constraints to x variable
+used when finding trim conditions
+'''
+
+from math import sin, cos, asin
+from tgear import tgear
+
+def conf16(x, u, const):
+ 'apply constraints to x'
+
+ radgam, singam, rr, pr, tr, phi, cphi, sphi, thetadot, coord, stab, orient = const
+ gamm = asin(singam)
+
+ #
+ # Steady Level Flight
+ #
+ if orient == 1:
+ x[3] = phi # Phi
+ x[4] = x[1] # Theta
+ x[6] = rr # Roll Rate
+ x[7] = pr # Pitch Rate
+ x[8] = 0.0 # Yaw Rate
+
+ #
+ # Steady Climb
+ #
+ if orient == 2:
+ x[3] = phi # Phi
+ x[4] = x[1] + radgam # Theta
+ x[6] = rr # Roll Rate
+ x[7] = pr # Pitch Rate
+ x[8] = 0.0 # Yaw Rate
+
+ #
+ # orient=3 implies coordinated turn
+ #
+ if orient == 3:
+ x[6] = -tr * sin(x[4]) # Roll Rate
+ x[7] = tr * cos(x[4]) * sin(x[3]) # Pitch Rate
+ x[8] = tr * cos(x[4]) * cos(x[3]) # Yaw Rate
+
+ #
+ # Pitch Pull Up
+ #
+ if orient == 4:
+ x[4] = x[1] # Theta = alpha
+ x[3] = phi # Phi
+ x[6] = rr # Roll Rate
+ x[7] = thetadot # Pitch Rate
+ x[8] = 0.0 # Yaw Rate
+
+ x[12] = tgear(u[0])
+
+ return x, u
diff --git a/demo/F16/aerobench/lowlevel/cx.py b/demo/F16/aerobench/lowlevel/cx.py
new file mode 100644
index 00000000..c23a96cb
--- /dev/null
+++ b/demo/F16/aerobench/lowlevel/cx.py
@@ -0,0 +1,50 @@
+'''
+Stanley Bak
+Python F-16 GCAS
+cx
+'''
+
+import numpy as np
+
+from aerobench.util import fix, sign
+
+def cx(alpha, el):
+ 'cx definition'
+
+ a = np.array([[-.099, -.081, -.081, -.063, -.025, .044, .097, .113, .145, .167, .174, .166], \
+ [-.048, -.038, -.040, -.021, .016, .083, .127, .137, .162, .177, .179, .167], \
+ [-.022, -.020, -.021, -.004, .032, .094, .128, .130, .154, .161, .155, .138], \
+ [-.040, -.038, -.039, -.025, .006, .062, .087, .085, .100, .110, .104, .091], \
+ [-.083, -.073, -.076, -.072, -.046, .012, .024, .025, .043, .053, .047, .040]], dtype=float).T
+
+ s = .2 * alpha
+ k = fix(s)
+ if k <= -2:
+ k = -1
+
+ if k >= 9:
+ k = 8
+
+ da = s - k
+ l = k + fix(1.1 * sign(da))
+ s = el / 12
+ m = fix(s)
+ if m <= -2:
+ m = -1
+
+ if m >= 2:
+ m = 1
+
+ de = s - m
+ n = m + fix(1.1 * sign(de))
+ k = k + 3
+ l = l + 3
+ m = m + 3
+ n = n + 3
+ t = a[k-1, m-1]
+ u = a[k-1, n-1]
+ v = t + abs(da) * (a[l-1, m-1] - t)
+ w = u + abs(da) * (a[l-1, n-1] - u)
+ cxx = v + (w - v) * abs(de)
+
+ return cxx
diff --git a/demo/F16/aerobench/lowlevel/cy.py b/demo/F16/aerobench/lowlevel/cy.py
new file mode 100644
index 00000000..fac8c8c3
--- /dev/null
+++ b/demo/F16/aerobench/lowlevel/cy.py
@@ -0,0 +1,9 @@
+'''
+Stanley Bak
+Python F-16 GCAS
+'''
+
+def cy(beta, ail, rdr):
+ 'cy function'
+
+ return -.02 * beta + .021 * (ail / 20) + .086 * (rdr / 30)
diff --git a/demo/F16/aerobench/lowlevel/cz.py b/demo/F16/aerobench/lowlevel/cz.py
new file mode 100644
index 00000000..66f57695
--- /dev/null
+++ b/demo/F16/aerobench/lowlevel/cz.py
@@ -0,0 +1,31 @@
+'''
+Stanley Bak
+Python F-16 GCAS
+Cz function
+'''
+
+import numpy as np
+from aerobench.util import fix, sign
+
+def cz(alpha, beta, el):
+ 'cz function'
+
+ a = np.array([.770, .241, -.100, -.415, -.731, -1.053, -1.355, -1.646, -1.917, -2.120, -2.248, -2.229], \
+ dtype=float).T
+
+ s = .2 * alpha
+ k = fix(s)
+
+ if k <= -2:
+ k = -1
+
+ if k >= 9:
+ k = 8
+
+ da = s - k
+ l = k + fix(1.1 * sign(da))
+ l = l + 3
+ k = k + 3
+ s = a[k-1] + abs(da) * (a[l-1] - a[k-1])
+
+ return s * (1 - (beta / 57.3)**2) - .19 * (el / 25)
diff --git a/demo/F16/aerobench/lowlevel/dampp.py b/demo/F16/aerobench/lowlevel/dampp.py
new file mode 100644
index 00000000..ffca4495
--- /dev/null
+++ b/demo/F16/aerobench/lowlevel/dampp.py
@@ -0,0 +1,42 @@
+'''
+Stanley Bak
+F16 GCAS in Python
+dampp function
+'''
+
+import numpy as np
+from aerobench.util import fix, sign
+
+def dampp(alpha):
+ 'dampp functon'
+
+ a = np.array([[-.267, -.110, .308, 1.34, 2.08, 2.91, 2.76, 2.05, 1.50, 1.49, 1.83, 1.21], \
+ [.882, .852, .876, .958, .962, .974, .819, .483, .590, 1.21, -.493, -1.04], \
+ [-.108, -.108, -.188, .110, .258, .226, .344, .362, .611, .529, .298, -2.27], \
+ [-8.80, -25.8, -28.9, -31.4, -31.2, -30.7, -27.7, -28.2, -29.0, -29.8, -38.3, -35.3], \
+ [-.126, -.026, .063, .113, .208, .230, .319, .437, .680, .100, .447, -.330], \
+ [-.360, -.359, -.443, -.420, -.383, -.375, -.329, -.294, -.230, -.210, -.120, -.100], \
+ [-7.21, -.540, -5.23, -5.26, -6.11, -6.64, -5.69, -6.00, -6.20, -6.40, -6.60, -6.00], \
+ [-.380, -.363, -.378, -.386, -.370, -.453, -.550, -.582, -.595, -.637, -1.02, -.840], \
+ [.061, .052, .052, -.012, -.013, -.024, .050, .150, .130, .158, .240, .150]], dtype=float).T
+
+ s = .2 * alpha
+ k = fix(s)
+
+ if k <= -2:
+ k = -1
+
+ if k >= 9:
+ k = 8
+
+ da = s - k
+ l = k + fix(1.1 * sign(da))
+ k = k + 3
+ l = l + 3
+
+ d = np.zeros((9,))
+
+ for i in range(9):
+ d[i] = a[k-1, i] + abs(da) * (a[l-1, i] - a[k-1, i])
+
+ return d
diff --git a/demo/F16/aerobench/lowlevel/dlda.py b/demo/F16/aerobench/lowlevel/dlda.py
new file mode 100644
index 00000000..1f6b424b
--- /dev/null
+++ b/demo/F16/aerobench/lowlevel/dlda.py
@@ -0,0 +1,50 @@
+'''
+Stanley Bak
+F-16 GCAS in Python
+dlda function
+'''
+
+import numpy as np
+from aerobench.util import fix, sign
+
+def dlda(alpha, beta):
+ 'dlda function'
+
+ a = np.array([[-.041, -.052, -.053, -.056, -.050, -.056, -.082, -.059, -.042, -.038, -.027, -.017], \
+ [-.041, -.053, -.053, -.053, -.050, -.051, -.066, -.043, -.038, -.027, -.023, -.016], \
+ [-.042, -.053, -.052, -.051, -.049, -.049, -.043, -.035, -.026, -.016, -.018, -.014], \
+ [-.040, -.052, -.051, -.052, -.048, -.048, -.042, -.037, -.031, -.026, -.017, -.012], \
+ [-.043, -.049, -.048, -.049, -.043, -.042, -.042, -.036, -.025, -.021, -.016, -.011], \
+ [-.044, -.048, -.048, -.047, -.042, -.041, -.020, -.028, -.013, -.014, -.011, -.010], \
+ [-.043, -.049, -.047, -.045, -.042, -.037, -.003, -.013, -.010, -.003, -.007, -.008]], dtype=float).T
+
+ s = .2 * alpha
+ k = fix(s)
+ if k <= -2:
+ k = -1
+
+ if k >= 9:
+ k = 8
+
+ da = s - k
+ l = k + fix(1.1 * sign(da))
+ s = .1 * beta
+ m = fix(s)
+ if m <= -3:
+ m = -2
+
+ if m >= 3:
+ m = 2
+
+ db = s - m
+ n = m + fix(1.1 * sign(db))
+ l = l + 3
+ k = k + 3
+ m = m + 4
+ n = n + 4
+ t = a[k-1, m-1]
+ u = a[k-1, n-1]
+ v = t + abs(da) * (a[l-1, m-1] - t)
+ w = u + abs(da) * (a[l-1, n-1] - u)
+
+ return v + (w - v) * abs(db)
diff --git a/demo/F16/aerobench/lowlevel/dldr.py b/demo/F16/aerobench/lowlevel/dldr.py
new file mode 100644
index 00000000..8d753484
--- /dev/null
+++ b/demo/F16/aerobench/lowlevel/dldr.py
@@ -0,0 +1,52 @@
+'''
+Stanley Bak
+Python GCAS F - 16
+dldr function
+'''
+
+import numpy as np
+from aerobench.util import sign, fix
+
+def dldr(alpha, beta):
+ 'dldr function'
+
+ a = np.array([[.005, .017, .014, .010, -.005, .009, .019, .005, -.000, -.005, -.011, .008], \
+ [.007, .016, .014, .014, .013, .009, .012, .005, .000, .004, .009, .007], \
+ [.013, .013, .011, .012, .011, .009, .008, .005, -.002, .005, .003, .005], \
+ [.018, .015, .015, .014, .014, .014, .014, .015, .013, .011, .006, .001], \
+ [.015, .014, .013, .013, .012, .011, .011, .010, .008, .008, .007, .003], \
+ [.021, .011, .010, .011, .010, .009, .008, .010, .006, .005, .000, .001], \
+ [.023, .010, .011, .011, .011, .010, .008, .010, .006, .014, .020, .000]], dtype=float).T
+
+ s = .2 * alpha
+ k = fix(s)
+ if k <= -2:
+ k = -1
+
+ if k >= 9:
+ k = 8
+
+ da = s - k
+ l = k + fix(1.1 * sign(da))
+ s = .1 * beta
+ m = fix(s)
+
+ if m <= -3:
+ m = -2
+
+ if m >= 3:
+ m = 2
+
+ db = s - m
+ n = m + fix(1.1 * sign(db))
+ l = l + 3
+ k = k + 3
+ m = m + 4
+ n = n + 4
+ t = a[k-1, m-1]
+ u = a[k-1, n-1]
+
+ v = t + abs(da) * (a[l-1, m-1] - t)
+ w = u + abs(da) * (a[l-1, n-1] - u)
+
+ return v + (w - v) * abs(db)
diff --git a/demo/F16/aerobench/lowlevel/dnda.py b/demo/F16/aerobench/lowlevel/dnda.py
new file mode 100644
index 00000000..d7e62355
--- /dev/null
+++ b/demo/F16/aerobench/lowlevel/dnda.py
@@ -0,0 +1,51 @@
+'''
+Stanley Bak
+F16 GCAS in Python
+dnda function
+'''
+
+import numpy as np
+from aerobench.util import fix, sign
+
+def dnda(alpha, beta):
+ 'dnda function'
+
+ a = np.array([[.001, -.027, -.017, -.013, -.012, -.016, .001, .017, .011, .017, .008, .016], \
+ [.002, -.014, -.016, -.016, -.014, -.019, -.021, .002, .012, .016, .015, .011], \
+ [-.006, -.008, -.006, -.006, -.005, -.008, -.005, .007, .004, .007, .006, .006], \
+ [-.011, -.011, -.010, -.009, -.008, -.006, .000, .004, .007, .010, .004, .010], \
+ [-.015, -.015, -.014, -.012, -.011, -.008, -.002, .002, .006, .012, .011, .011], \
+ [-.024, -.010, -.004, -.002, -.001, .003, .014, .006, -.001, .004, .004, .006], \
+ [-.022, .002, -.003, -.005, -.003, -.001, -.009, -.009, -.001, .003, -.002, .001]], dtype=float).T
+
+ s = .2 * alpha
+ k = fix(s)
+
+ if k <= -2:
+ k = -1
+
+ if k >= 9:
+ k = 8
+
+ da = s - k
+ l = k + fix(1.1 * sign(da))
+ s = .1 * beta
+ m = fix(s)
+ if m <= -3:
+ m = -2
+
+ if m >= 3:
+ m = 2
+
+ db = s - m
+ n = m + fix(1.1 * sign(db))
+ l = l + 3
+ k = k + 3
+ m = m + 4
+ n = n + 4
+ t = a[k-1, m-1]
+ u = a[k-1, n-1]
+ v = t + abs(da) * (a[l-1, m-1] - t)
+ w = u + abs(da) * (a[l-1, n-1] - u)
+
+ return v + (w - v) * abs(db)
diff --git a/demo/F16/aerobench/lowlevel/dndr.py b/demo/F16/aerobench/lowlevel/dndr.py
new file mode 100644
index 00000000..805f1709
--- /dev/null
+++ b/demo/F16/aerobench/lowlevel/dndr.py
@@ -0,0 +1,49 @@
+'''
+Stanley Bak
+F16 GCAS in Python
+dndr function
+'''
+
+import numpy as np
+from aerobench.util import fix, sign
+
+def dndr(alpha, beta):
+ 'dndr function'
+
+ a = np.array([[-.018, -.052, -.052, -.052, -.054, -.049, -.059, -.051, -.030, -.037, -.026, -.013], \
+ [-.028, -.051, -.043, -.046, -.045, -.049, -.057, -.052, -.030, -.033, -.030, -.008], \
+ [-.037, -.041, -.038, -.040, -.040, -.038, -.037, -.030, -.027, -.024, -.019, -.013], \
+ [-.048, -.045, -.045, -.045, -.044, -.045, -.047, -.048, -.049, -.045, -.033, -.016], \
+ [-.043, -.044, -.041, -.041, -.040, -.038, -.034, -.035, -.035, -.029, -.022, -.009], \
+ [-.052, -.034, -.036, -.036, -.035, -.028, -.024, -.023, -.020, -.016, -.010, -.014], \
+ [-.062, -.034, -.027, -.028, -.027, -.027, -.023, -.023, -.019, -.009, -.025, -.010]], dtype=float).T
+
+ s = .2 * alpha
+ k = fix(s)
+ if k <= -2:
+ k = -1
+
+ if k >= 9:
+ k = 8
+
+ da = s - k
+ l = k + fix(1.1 * sign(da))
+ s = .1 * beta
+ m = fix(s)
+ if m <= -3:
+ m = -2
+
+ if m >= 3:
+ m = 2
+
+ db = s - m
+ n = m + fix(1.1 * sign(db))
+ l = l + 3
+ k = k + 3
+ m = m + 4
+ n = n + 4
+ t = a[k-1, m-1]
+ u = a[k-1, n-1]
+ v = t + abs(da) * (a[l-1, m-1] - t)
+ w = u + abs(da) * (a[l-1, n-1] - u)
+ return v + (w - v) * abs(db)
diff --git a/demo/F16/aerobench/lowlevel/low_level_controller.py b/demo/F16/aerobench/lowlevel/low_level_controller.py
new file mode 100644
index 00000000..d472bf30
--- /dev/null
+++ b/demo/F16/aerobench/lowlevel/low_level_controller.py
@@ -0,0 +1,109 @@
+'''
+Stanley Bak
+Low-level flight controller
+'''
+
+import numpy as np
+from aerobench.util import Freezable
+
+class CtrlLimits(Freezable):
+ 'Control Limits'
+
+ def __init__(self):
+ self.ThrottleMax = 1 # Afterburner on for throttle > 0.7
+ self.ThrottleMin = 0
+ self.ElevatorMaxDeg = 25
+ self.ElevatorMinDeg = -25
+ self.AileronMaxDeg = 21.5
+ self.AileronMinDeg = -21.5
+ self.RudderMaxDeg = 30
+ self.RudderMinDeg = -30
+
+ self.NzMax = 6
+ self.NzMin = -1
+
+ self.freeze_attrs()
+
+class LowLevelController(Freezable):
+ '''low level flight controller
+ '''
+
+ old_k_long = np.array([[-156.8801506723475, -31.037008068526642, -38.72983346216317]], dtype=float)
+ old_k_lat = np.array([[37.84483, -25.40956, -6.82876, -332.88343, -17.15997],
+ [-23.91233, 5.69968, -21.63431, 64.49490, -88.36203]], dtype=float)
+
+ old_xequil = np.array([502.0, 0.0389, 0.0, 0.0, 0.0389, 0.0, 0.0, 0.0, \
+ 0.0, 0.0, 0.0, 1000.0, 9.0567], dtype=float).transpose()
+ old_uequil = np.array([0.1395, -0.7496, 0.0, 0.0], dtype=float).transpose()
+
+ def __init__(self, gain_str='old'):
+ # Hard coded LQR gain matrix from matlab version
+
+ assert gain_str == 'old'
+
+ # Longitudinal Gains
+ K_long = LowLevelController.old_k_long
+ K_lat = LowLevelController.old_k_lat
+
+ self.K_lqr = np.zeros((3, 8))
+ self.K_lqr[:1, :3] = K_long
+ self.K_lqr[1:, 3:] = K_lat
+
+ # equilibrium points from BuildLqrControllers.py
+ self.xequil = LowLevelController.old_xequil
+ self.uequil = LowLevelController.old_uequil
+
+ self.ctrlLimits = CtrlLimits()
+
+ self.freeze_attrs()
+
+ def get_u_deg(self, u_ref, f16_state):
+ 'get the reference commands for the control surfaces'
+
+ # Calculate perturbation from trim state
+ x_delta = f16_state.copy()
+ x_delta[:len(self.xequil)] -= self.xequil
+
+ ## Implement LQR Feedback Control
+ # Reorder states to match controller:
+ # [alpha, q, int_e_Nz, beta, p, r, int_e_ps, int_e_Ny_r]
+ x_ctrl = np.array([x_delta[i] for i in [1, 7, 13, 2, 6, 8, 14, 15]], dtype=float)
+
+ # Initialize control vectors
+ u_deg = np.zeros((4,)) # throt, ele, ail, rud
+
+ # Calculate control using LQR gains
+ u_deg[1:4] = np.dot(-self.K_lqr, x_ctrl) # Full Control
+
+ # Set throttle as directed from output of getOuterLoopCtrl(...)
+ u_deg[0] = u_ref[3]
+
+ # Add in equilibrium control
+ u_deg[0:4] += self.uequil
+
+ ## Limit controls to saturation limits
+ ctrlLimits = self.ctrlLimits
+
+ # Limit throttle from 0 to 1
+ u_deg[0] = max(min(u_deg[0], ctrlLimits.ThrottleMax), ctrlLimits.ThrottleMin)
+
+ # Limit elevator from -25 to 25 deg
+ u_deg[1] = max(min(u_deg[1], ctrlLimits.ElevatorMaxDeg), ctrlLimits.ElevatorMinDeg)
+
+ # Limit aileron from -21.5 to 21.5 deg
+ u_deg[2] = max(min(u_deg[2], ctrlLimits.AileronMaxDeg), ctrlLimits.AileronMinDeg)
+
+ # Limit rudder from -30 to 30 deg
+ u_deg[3] = max(min(u_deg[3], ctrlLimits.RudderMaxDeg), ctrlLimits.RudderMinDeg)
+
+ return x_ctrl, u_deg
+
+ def get_num_integrators(self):
+ 'get the number of integrators in the low-level controller'
+
+ return 3
+
+ def get_integrator_derivatives(self, t, x_f16, u_ref, Nz, ps, Ny_r):
+ 'get the derivatives of the integrators in the low-level controller'
+
+ return [Nz - u_ref[0], ps - u_ref[1], Ny_r - u_ref[2]]
diff --git a/demo/F16/aerobench/lowlevel/morellif16.py b/demo/F16/aerobench/lowlevel/morellif16.py
new file mode 100644
index 00000000..4632b78c
--- /dev/null
+++ b/demo/F16/aerobench/lowlevel/morellif16.py
@@ -0,0 +1,189 @@
+'''
+Stanley Bak
+F16 GCAS in Python
+
+Morelli dynamics (Polynomial interpolation)
+'''
+
+def Morellif16(alpha, beta, de, da, dr, p, q, r, cbar, b, V, xcg, xcgref):
+ 'desc'
+
+ #alpha=max(-10*pi/180,min(45*pi/180,alpha)) # bounds alpha between -10 deg and 45 deg
+ #beta = max( - 30 * pi / 180, min(30 * pi / 180, beta)) #bounds beta between -30 deg and 30 deg
+ #de = max( - 25 * pi / 180, min(25 * pi / 180, de)) #bounds elevator deflection between -25 deg and 25 deg
+ #da = max( - 21.5 * pi / 180, min(21.5 * pi / 180, da)) #bounds aileron deflection between -21.5 deg and 21.5 deg
+ #dr = max( - 30 * pi / 180, min(30 * pi / 180, dr)) #bounds rudder deflection between -30 deg and 30 deg
+
+ # xcgref = 0.35
+ #reference longitudinal cg position in Morelli f16 model
+
+
+ phat = p * b / (2 * V)
+ qhat = q * cbar / (2 * V)
+ rhat = r * b / (2 * V)
+ ##
+ a0 = -1.943367e-2
+ a1 = 2.136104e-1
+ a2 = -2.903457e-1
+ a3 = -3.348641e-3
+ a4 = -2.060504e-1
+ a5 = 6.988016e-1
+ a6 = -9.035381e-1
+
+ b0 = 4.833383e-1
+ b1 = 8.644627
+ b2 = 1.131098e1
+ b3 = -7.422961e1
+ b4 = 6.075776e1
+
+ c0 = -1.145916
+ c1 = 6.016057e-2
+ c2 = 1.642479e-1
+
+ d0 = -1.006733e-1
+ d1 = 8.679799e-1
+ d2 = 4.260586
+ d3 = -6.923267
+
+ e0 = 8.071648e-1
+ e1 = 1.189633e-1
+ e2 = 4.177702
+ e3 = -9.162236
+
+ f0 = -1.378278e-1
+ f1 = -4.211369
+ f2 = 4.775187
+ f3 = -1.026225e1
+ f4 = 8.399763
+ f5 = -4.354000e-1
+
+ g0 = -3.054956e1
+ g1 = -4.132305e1
+ g2 = 3.292788e2
+ g3 = -6.848038e2
+ g4 = 4.080244e2
+
+ h0 = -1.05853e-1
+ h1 = -5.776677e-1
+ h2 = -1.672435e-2
+ h3 = 1.357256e-1
+ h4 = 2.172952e-1
+ h5 = 3.464156
+ h6 = -2.835451
+ h7 = -1.098104
+
+ i0 = -4.126806e-1
+ i1 = -1.189974e-1
+ i2 = 1.247721
+ i3 = -7.391132e-1
+
+ j0 = 6.250437e-2
+ j1 = 6.067723e-1
+ j2 = -1.101964
+ j3 = 9.100087
+ j4 = -1.192672e1
+
+ k0 = -1.463144e-1
+ k1 = -4.07391e-2
+ k2 = 3.253159e-2
+ k3 = 4.851209e-1
+ k4 = 2.978850e-1
+ k5 = -3.746393e-1
+ k6 = -3.213068e-1
+
+ l0 = 2.635729e-2
+ l1 = -2.192910e-2
+ l2 = -3.152901e-3
+ l3 = -5.817803e-2
+ l4 = 4.516159e-1
+ l5 = -4.928702e-1
+ l6 = -1.579864e-2
+
+ m0 = -2.029370e-2
+ m1 = 4.660702e-2
+ m2 = -6.012308e-1
+ m3 = -8.062977e-2
+ m4 = 8.320429e-2
+ m5 = 5.018538e-1
+ m6 = 6.378864e-1
+ m7 = 4.226356e-1
+
+ n0 = -5.19153
+ n1 = -3.554716
+ n2 = -3.598636e1
+ n3 = 2.247355e2
+ n4 = -4.120991e2
+ n5 = 2.411750e2
+
+ o0 = 2.993363e-1
+ o1 = 6.594004e-2
+ o2 = -2.003125e-1
+ o3 = -6.233977e-2
+ o4 = -2.107885
+ o5 = 2.141420
+ o6 = 8.476901e-1
+
+ p0 = 2.677652e-2
+ p1 = -3.298246e-1
+ p2 = 1.926178e-1
+ p3 = 4.013325
+ p4 = -4.404302
+
+ q0 = -3.698756e-1
+ q1 = -1.167551e-1
+ q2 = -7.641297e-1
+
+ r0 = -3.348717e-2
+ r1 = 4.276655e-2
+ r2 = 6.573646e-3
+ r3 = 3.535831e-1
+ r4 = -1.373308
+ r5 = 1.237582
+ r6 = 2.302543e-1
+ r7 = -2.512876e-1
+ r8 = 1.588105e-1
+ r9 = -5.199526e-1
+
+ s0 = -8.115894e-2
+ s1 = -1.156580e-2
+ s2 = 2.514167e-2
+ s3 = 2.038748e-1
+ s4 = -3.337476e-1
+ s5 = 1.004297e-1
+
+ ##
+ Cx0 = a0 + a1 * alpha + a2 * de**2 + a3 * de + a4 * alpha * de + a5 * alpha**2 + a6 * alpha**3
+ Cxq = b0 + b1 * alpha + b2 * alpha**2 + b3 * alpha**3 + b4 * alpha**4
+ Cy0 = c0 * beta + c1 * da + c2 * dr
+ Cyp = d0 + d1 * alpha + d2 * alpha**2 + d3 * alpha**3
+ Cyr = e0 + e1 * alpha + e2 * alpha**2 + e3 * alpha**3
+ Cz0 = (f0 + f1 * alpha + f2 * alpha**2 + f3 * alpha**3 + f4 * alpha**4) * (1 - beta**2) + f5 * de
+ Czq = g0 + g1 * alpha + g2 * alpha**2 + g3 * alpha**3 + g4 * alpha**4
+ Cl0 = h0 * beta + h1 * alpha * beta + h2 * alpha**2 * beta + h3 * beta**2 + h4 * alpha * beta**2 + h5 * \
+ alpha**3 * beta + h6 * alpha**4 * beta + h7 * alpha**2 * beta**2
+ Clp = i0 + i1 * alpha + i2 * alpha**2 + i3 * alpha**3
+ Clr = j0 + j1 * alpha + j2 * alpha**2 + j3 * alpha**3 + j4 * alpha**4
+ Clda = k0 + k1 * alpha + k2 * beta + k3 * alpha**2 + k4 * alpha * beta + k5 * alpha**2 * beta + k6 * alpha**3
+ Cldr = l0 + l1 * alpha + l2 * beta + l3 * alpha * beta + l4 * alpha**2 * beta + l5 * alpha**3 * beta + l6 * beta**2
+ Cm0 = m0 + m1 * alpha + m2 * de + m3 * alpha * de + m4 * de**2 + m5 * alpha**2 * de + m6 * de**3 + m7 * \
+ alpha * de**2
+
+
+ Cmq = n0 + n1 * alpha + n2 * alpha**2 + n3 * alpha**3 + n4 * alpha**4 + n5 * alpha**5
+ Cn0 = o0 * beta + o1 * alpha * beta + o2 * beta**2 + o3 * alpha * beta**2 + o4 * alpha**2 * beta + o5 * \
+ alpha**2 * beta**2 + o6 * alpha**3 * beta
+ Cnp = p0 + p1 * alpha + p2 * alpha**2 + p3 * alpha**3 + p4 * alpha**4
+ Cnr = q0 + q1 * alpha + q2 * alpha**2
+ Cnda = r0 + r1 * alpha + r2 * beta + r3 * alpha * beta + r4 * alpha**2 * beta + r5 * alpha**3 * beta + r6 * \
+ alpha**2 + r7 * alpha**3 + r8 * beta**3 + r9 * alpha * beta**3
+ Cndr = s0 + s1 * alpha + s2 * beta + s3 * alpha * beta + s4 * alpha**2 * beta + s5 * alpha**2
+ ##
+
+ Cx = Cx0 + Cxq * qhat
+ Cy = Cy0 + Cyp * phat + Cyr * rhat
+ Cz = Cz0 + Czq * qhat
+ Cl = Cl0 + Clp * phat + Clr * rhat + Clda * da + Cldr * dr
+ Cm = Cm0 + Cmq * qhat + Cz * (xcgref - xcg)
+ Cn = Cn0 + Cnp * phat + Cnr * rhat + Cnda * da + Cndr * dr - Cy * (xcgref - xcg) * (cbar / b)
+
+ return Cx, Cy, Cz, Cl, Cm, Cn
diff --git a/demo/F16/aerobench/lowlevel/pdot.py b/demo/F16/aerobench/lowlevel/pdot.py
new file mode 100644
index 00000000..dbfb1c32
--- /dev/null
+++ b/demo/F16/aerobench/lowlevel/pdot.py
@@ -0,0 +1,29 @@
+'''
+Stanley Bak
+Python F-16
+power derivative (pdot)
+'''
+
+from aerobench.lowlevel.rtau import rtau
+
+def pdot(p3, p1):
+ 'pdot function'
+
+ if p1 >= 50:
+ if p3 >= 50:
+ t = 5
+ p2 = p1
+ else:
+ p2 = 60
+ t = rtau(p2 - p3)
+ else:
+ if p3 >= 50:
+ t = 5
+ p2 = 40
+ else:
+ p2 = p1
+ t = rtau(p2 - p3)
+
+ pd = t * (p2 - p3)
+
+ return pd
diff --git a/demo/F16/aerobench/lowlevel/rtau.py b/demo/F16/aerobench/lowlevel/rtau.py
new file mode 100644
index 00000000..a2f00d6e
--- /dev/null
+++ b/demo/F16/aerobench/lowlevel/rtau.py
@@ -0,0 +1,18 @@
+'''
+Stanley Bak
+Python F-16
+
+Rtau function
+'''
+
+def rtau(dp):
+ 'rtau function'
+
+ if dp <= 25:
+ rt = 1.0
+ elif dp >= 50:
+ rt = .1
+ else:
+ rt = 1.9 - .036 * dp
+
+ return rt
diff --git a/demo/F16/aerobench/lowlevel/subf16_model.py b/demo/F16/aerobench/lowlevel/subf16_model.py
new file mode 100644
index 00000000..b807b165
--- /dev/null
+++ b/demo/F16/aerobench/lowlevel/subf16_model.py
@@ -0,0 +1,201 @@
+'''
+Stanley Bak
+Python F-16 subf16
+outputs aircraft state vector deriative
+'''
+
+# x[0] = air speed, VT (ft/sec)
+# x[1] = angle of attack, alpha (rad)
+# x[2] = angle of sideslip, beta (rad)
+# x[3] = roll angle, phi (rad)
+# x[4] = pitch angle, theta (rad)
+# x[5] = yaw angle, psi (rad)
+# x[6] = roll rate, P (rad/sec)
+# x[7] = pitch rate, Q (rad/sec)
+# x[8] = yaw rate, R (rad/sec)
+# x[9] = northward horizontal displacement, pn (feet)
+# x[10] = eastward horizontal displacement, pe (feet)
+# x[11] = altitude, h (feet)
+# x[12] = engine thrust dynamics lag state, pow
+#
+# u[0] = throttle command 0.0 < u(1) < 1.0
+# u[1] = elevator command in degrees
+# u[2] = aileron command in degrees
+# u[3] = rudder command in degrees
+#
+
+from math import sin, cos, pi
+
+from aerobench.lowlevel.adc import adc
+from aerobench.lowlevel.tgear import tgear
+from aerobench.lowlevel.pdot import pdot
+from aerobench.lowlevel.thrust import thrust
+from aerobench.lowlevel.cx import cx
+from aerobench.lowlevel.cy import cy
+from aerobench.lowlevel.cz import cz
+from aerobench.lowlevel.cl import cl
+from aerobench.lowlevel.dlda import dlda
+from aerobench.lowlevel.dldr import dldr
+from aerobench.lowlevel.cm import cm
+from aerobench.lowlevel.cn import cn
+from aerobench.lowlevel.dnda import dnda
+from aerobench.lowlevel.dndr import dndr
+from aerobench.lowlevel.dampp import dampp
+
+from aerobench.lowlevel.morellif16 import Morellif16
+
+def subf16_model(x, u, model, adjust_cy=True):
+ '''output aircraft state vector derivative for a given input
+
+ The reference for the model is Appendix A of Stevens & Lewis
+ '''
+
+ assert model in ['stevens', 'morelli']
+ assert len(x) == 13
+ assert len(u) == 4
+
+ xcg = 0.35
+
+ thtlc, el, ail, rdr = u
+
+ s = 300
+ b = 30
+ cbar = 11.32
+ rm = 1.57e-3
+ xcgr = .35
+ he = 160.0
+ c1 = -.770
+ c2 = .02755
+ c3 = 1.055e-4
+ c4 = 1.642e-6
+ c5 = .9604
+ c6 = 1.759e-2
+ c7 = 1.792e-5
+ c8 = -.7336
+ c9 = 1.587e-5
+ rtod = 57.29578
+ g = 32.17
+
+ xd = x.copy()
+ vt = x[0]
+ alpha = x[1]*rtod
+ beta = x[2]*rtod
+ phi = x[3]
+ theta = x[4]
+ psi = x[5]
+ p = x[6]
+ q = x[7]
+ r = x[8]
+ alt = x[11]
+ power = x[12]
+
+ # air data computer and engine model
+ amach, qbar = adc(vt, alt)
+ cpow = tgear(thtlc)
+
+ xd[12] = pdot(power, cpow)
+
+ t = thrust(power, alt, amach)
+ dail = ail/20
+ drdr = rdr/30
+
+ # component build up
+
+ if model == 'stevens':
+ # stevens & lewis (look up table version)
+ cxt = cx(alpha, el)
+ cyt = cy(beta, ail, rdr)
+ czt = cz(alpha, beta, el)
+
+ clt = cl(alpha, beta) + dlda(alpha, beta) * dail + dldr(alpha, beta) * drdr
+ cmt = cm(alpha, el)
+ cnt = cn(alpha, beta) + dnda(alpha, beta) * dail + dndr(alpha, beta) * drdr
+ else:
+ # morelli model (polynomial version)
+ cxt, cyt, czt, clt, cmt, cnt = Morellif16(alpha*pi/180, beta*pi/180, el*pi/180, ail*pi/180, rdr*pi/180, \
+ p, q, r, cbar, b, vt, xcg, xcgr)
+
+ # add damping derivatives
+
+ tvt = .5 / vt
+ b2v = b * tvt
+ cq = cbar * q * tvt
+
+ # get ready for state equations
+ d = dampp(alpha)
+ cxt = cxt + cq * d[0]
+ cyt = cyt + b2v * (d[1] * r + d[2] * p)
+ czt = czt + cq * d[3]
+ clt = clt + b2v * (d[4] * r + d[5] * p)
+ cmt = cmt + cq * d[6] + czt * (xcgr-xcg)
+ cnt = cnt + b2v * (d[7] * r + d[8] * p)-cyt * (xcgr-xcg) * cbar/b
+ cbta = cos(x[2])
+ u = vt * cos(x[1]) * cbta
+ v = vt * sin(x[2])
+ w = vt * sin(x[1]) * cbta
+ sth = sin(theta)
+ cth = cos(theta)
+ sph = sin(phi)
+ cph = cos(phi)
+ spsi = sin(psi)
+ cpsi = cos(psi)
+ qs = qbar * s
+ qsb = qs * b
+ rmqs = rm * qs
+ gcth = g * cth
+ qsph = q * sph
+ ay = rmqs * cyt
+ az = rmqs * czt
+
+ # force equations
+ udot = r * v-q * w-g * sth + rm * (qs * cxt + t)
+ vdot = p * w-r * u + gcth * sph + ay
+ wdot = q * u-p * v + gcth * cph + az
+ dum = (u * u + w * w)
+
+ xd[0] = (u * udot + v * vdot + w * wdot)/vt
+ xd[1] = (u * wdot-w * udot)/dum
+ xd[2] = (vt * vdot-v * xd[0]) * cbta/dum
+
+ # kinematics
+ xd[3] = p + (sth/cth) * (qsph + r * cph)
+ xd[4] = q * cph-r * sph
+ xd[5] = (qsph + r * cph)/cth
+
+ # moments
+ xd[6] = (c2 * p + c1 * r + c4 * he) * q + qsb * (c3 * clt + c4 * cnt)
+
+ xd[7] = (c5 * p-c7 * he) * r + c6 * (r * r-p * p) + qs * cbar * c7 * cmt
+ xd[8] = (c8 * p-c2 * r + c9 * he) * q + qsb * (c4 * clt + c9 * cnt)
+
+ # navigation
+ t1 = sph * cpsi
+ t2 = cph * sth
+ t3 = sph * spsi
+ s1 = cth * cpsi
+ s2 = cth * spsi
+ s3 = t1 * sth-cph * spsi
+ s4 = t3 * sth + cph * cpsi
+ s5 = sph * cth
+ s6 = t2 * cpsi + t3
+ s7 = t2 * spsi-t1
+ s8 = cph * cth
+ xd[9] = u * s1 + v * s3 + w * s6 # north speed
+ xd[10] = u * s2 + v * s4 + w * s7 # east speed
+ xd[11] = u * sth-v * s5-w * s8 # vertical speed
+
+ # outputs
+
+ xa = 15.0 # sets distance normal accel is in front of the c.g. (xa = 15.0 at pilot)
+ az = az-xa * xd[7] # moves normal accel in front of c.g.
+
+ ####################################
+ ###### peter additions below ######
+ if adjust_cy:
+ ay = ay+xa*xd[8] # moves side accel in front of c.g.
+
+ # For extraction of Nz
+ Nz = (-az / g) - 1 # zeroed at 1 g, positive g = pulling up
+ Ny = ay / g
+
+ return xd, Nz, Ny, az, ay
diff --git a/demo/F16/aerobench/lowlevel/tgear.py b/demo/F16/aerobench/lowlevel/tgear.py
new file mode 100644
index 00000000..8c0fe473
--- /dev/null
+++ b/demo/F16/aerobench/lowlevel/tgear.py
@@ -0,0 +1,14 @@
+'''
+Stanley Bak
+Python F-16 GCAS
+'''
+
+def tgear(thtl):
+ 'tgear function'
+
+ if thtl <= .77:
+ tg = 64.94 * thtl
+ else:
+ tg = 217.38 * thtl - 117.38
+
+ return tg
diff --git a/demo/F16/aerobench/lowlevel/thrust.py b/demo/F16/aerobench/lowlevel/thrust.py
new file mode 100644
index 00000000..3de5f8b1
--- /dev/null
+++ b/demo/F16/aerobench/lowlevel/thrust.py
@@ -0,0 +1,76 @@
+'''
+Stanle Bak
+Python F-16
+Thrust function
+'''
+
+import numpy as np
+
+from aerobench.util import fix
+
+def thrust(power, alt, rmach):
+ 'thrust lookup-table version'
+
+ a = np.array([[1060, 670, 880, 1140, 1500, 1860], \
+ [635, 425, 690, 1010, 1330, 1700], \
+ [60, 25, 345, 755, 1130, 1525], \
+ [-1020, -170, -300, 350, 910, 1360], \
+ [-2700, -1900, -1300, -247, 600, 1100], \
+ [-3600, -1400, -595, -342, -200, 700]], dtype=float).T
+
+ b = np.array([[12680, 9150, 6200, 3950, 2450, 1400], \
+ [12680, 9150, 6313, 4040, 2470, 1400], \
+ [12610, 9312, 6610, 4290, 2600, 1560], \
+ [12640, 9839, 7090, 4660, 2840, 1660], \
+ [12390, 10176, 7750, 5320, 3250, 1930], \
+ [11680, 9848, 8050, 6100, 3800, 2310]], dtype=float).T
+
+ c = np.array([[20000, 15000, 10800, 7000, 4000, 2500], \
+ [21420, 15700, 11225, 7323, 4435, 2600], \
+ [22700, 16860, 12250, 8154, 5000, 2835], \
+ [24240, 18910, 13760, 9285, 5700, 3215], \
+ [26070, 21075, 15975, 11115, 6860, 3950], \
+ [28886, 23319, 18300, 13484, 8642, 5057]], dtype=float).T
+
+ if alt < 0:
+ alt = 0.01 # uh, why not 0?
+
+ h = .0001 * alt
+
+ i = fix(h)
+
+ if i >= 5:
+ i = 4
+
+ dh = h - i
+ rm = 5 * rmach
+ m = fix(rm)
+
+ if m >= 5:
+ m = 4
+ elif m <= 0:
+ m = 0
+
+ dm = rm - m
+ cdh = 1 - dh
+
+ # do not increment these, since python is 0-indexed while matlab is 1-indexed
+ #i = i + 1
+ #m = m + 1
+
+ s = b[i, m] * cdh + b[i + 1, m] * dh
+ t = b[i, m + 1] * cdh + b[i + 1, m + 1] * dh
+ tmil = s + (t - s) * dm
+
+ if power < 50:
+ s = a[i, m] * cdh + a[i + 1, m] * dh
+ t = a[i, m + 1] * cdh + a[i + 1, m + 1] * dh
+ tidl = s + (t - s) * dm
+ thrst = tidl + (tmil - tidl) * power * .02
+ else:
+ s = c[i, m] * cdh + c[i + 1, m] * dh
+ t = c[i, m + 1] * cdh + c[i + 1, m + 1] * dh
+ tmax = s + (t - s) * dm
+ thrst = tmil + (tmax - tmil) * (power - 50) * .02
+
+ return thrst
diff --git a/demo/F16/aerobench/util.py b/demo/F16/aerobench/util.py
new file mode 100644
index 00000000..cf2c4796
--- /dev/null
+++ b/demo/F16/aerobench/util.py
@@ -0,0 +1,287 @@
+'''
+Utilities for F-16 GCAS
+'''
+
+from math import floor, ceil
+import numpy as np
+
+class StateIndex:
+ 'list of static state indices'
+
+ VT = 0
+ VEL = 0 # alias
+
+ ALPHA = 1
+ BETA = 2
+ PHI = 3 # roll angle
+ THETA = 4 # pitch angle
+ PSI = 5 # yaw angle
+
+ P = 6
+ Q = 7
+ R = 8
+
+ POSN = 9
+ POS_N = 9
+
+ POSE = 10
+ POS_E = 10
+
+ ALT = 11
+ H = 11
+
+ POW = 12
+
+class Freezable():
+ 'a class where you can freeze the fields (prevent new fields from being created)'
+
+ _frozen = False
+
+ def freeze_attrs(self):
+ 'prevents any new attributes from being created in the object'
+ self._frozen = True
+
+ def __setattr__(self, key, value):
+ if self._frozen and not hasattr(self, key):
+ raise TypeError("{} does not contain attribute '{}' (object was frozen)".format(self, key))
+
+ object.__setattr__(self, key, value)
+
+class Euler(Freezable):
+ '''fixed step euler integration
+
+ loosely based on scipy.integrate.RK45
+ '''
+
+ def __init__(self, der_func, tstart, ystart, tend, step=0, time_tol=1e-9):
+ assert step > 0, "arg step > 0 required in Euler integrator"
+ assert tend > tstart
+
+ self.der_func = der_func # signature (t, x)
+ self.tstep = step
+ self.t = tstart
+ self.y = ystart.copy()
+ self.yprev = None
+ self.tprev = None
+ self.tend = tend
+
+ self.status = 'running'
+
+ self.time_tol = time_tol
+
+ self.freeze_attrs()
+
+ def step(self):
+ 'take one step'
+
+ if self.status == 'running':
+ self.yprev = self.y.copy()
+ self.tprev = self.t
+ yd = self.der_func(self.t, self.y)
+
+ self.t += self.tstep
+
+ if self.t + self.time_tol >= self.tend:
+ self.t = self.tend
+
+ dt = self.t - self.tprev
+ self.y += dt * yd
+
+ if self.t == self.tend:
+ self.status = 'finished'
+
+ def dense_output(self):
+ 'return a function of time'
+
+ assert self.tprev is not None
+
+ dy = self.y - self.yprev
+ dt = self.t - self.tprev
+
+ dydt = dy / dt
+
+ def fun(t):
+ 'return state at time t (linear interpolation)'
+
+ deltat = t - self.tprev
+
+ return self.yprev + dydt * deltat
+
+ return fun
+
+def get_state_names():
+ 'returns a list of state variable names'
+
+ return ['vt', 'alpha', 'beta', 'phi', 'theta', 'psi', 'P', 'Q', 'R', 'pos_n', 'pos_e', 'alt', 'pow']
+
+def printmat(mat, main_label, row_label_str, col_label_str):
+ 'print a matrix'
+
+ if isinstance(row_label_str, list) and len(row_label_str) == 0:
+ row_label_str = None
+
+ assert isinstance(main_label, str)
+ assert row_label_str is None or isinstance(row_label_str, str)
+ assert isinstance(col_label_str, str)
+
+ mat = np.array(mat)
+ if len(mat.shape) == 1:
+ mat.shape = (1, mat.shape[0]) # one-row matrix
+
+ print("{main_label} =")
+
+ row_labels = None if row_label_str is None else row_label_str.split(' ')
+ col_labels = col_label_str.split(' ')
+
+ width = 7
+
+ width = max(width, max([len(l) for l in col_labels]))
+
+ if row_labels is not None:
+ width = max(width, max([len(l) for l in row_labels]))
+
+ width += 1
+
+ # add blank space for row labels
+ if row_labels is not None:
+ print("{: <{}}".format('', width), end='')
+
+ # print col lables
+ for col_label in col_labels:
+ if len(col_label) > width:
+ col_label = col_label[:width]
+
+ print("{: >{}}".format(col_label, width), end='')
+
+ print('')
+
+ if row_labels is not None:
+ assert len(row_labels) == mat.shape[0], \
+ "row labels (len={}) expected one element for each row of the matrix ({})".format( \
+ len(row_labels), mat.shape[0])
+
+ for r in range(mat.shape[0]):
+ row = mat[r]
+
+ if row_labels is not None:
+ label = row_labels[r]
+
+ if len(label) > width:
+ label = label[:width]
+
+ print("{:<{}}".format(label, width), end='')
+
+ for num in row:
+ #print("{:#<{}}".format(num, width), end='')
+ print("{:{}.{}g}".format(num, width, width-3), end='')
+
+ print('')
+
+
+def fix(ele):
+ 'round towards zero'
+
+ assert isinstance(ele, float)
+
+ if ele > 0:
+ rv = int(floor(ele))
+ else:
+ rv = int(ceil(ele))
+
+ return rv
+
+def sign(ele):
+ 'sign of a number'
+
+ if ele < 0:
+ rv = -1
+ elif ele == 0:
+ rv = 0
+ else:
+ rv = 1
+
+ return rv
+
+def extract_single_result(res, index, llc):
+ 'extract a res object for a sinlge aircraft from a multi-aircraft simulation'
+
+ num_vars = len(get_state_names()) + llc.get_num_integrators()
+ num_aircraft = res['states'][0].size // num_vars
+
+ if num_aircraft == 1:
+ assert index == 0
+ rv = res
+ else:
+ rv = {}
+ rv['status'] = res['status']
+ rv['times'] = res['times']
+ rv['modes'] = res['modes']
+
+ full_states = res['states']
+ rv['states'] = full_states[:, num_vars*index:num_vars*(index+1)]
+
+ if 'xd_list' in res:
+ # extended states
+ key_list = ['xd_list', 'ps_list', 'Nz_list', 'Ny_r_list', 'u_list']
+
+ for key in key_list:
+ rv[key] = [tup[index] for tup in res[key]]
+
+ return rv
+
+class SafetyLimits(Freezable):
+ 'a class for holding a set of safety limits.'
+
+ def __init__(self, **kwargs):
+ self.altitude = kwargs['altitude'] if 'altitude' in kwargs and kwargs['altitude'] is not None else None
+ self.Nz = kwargs['Nz'] if 'Nz' in kwargs and kwargs['Nz'] is not None else None
+ self.v = kwargs['v'] if 'v' in kwargs and kwargs['v'] is not None else None
+ self.alpha = kwargs['alpha'] if 'alpha' in kwargs and kwargs['alpha'] is not None else None
+
+ self.psMaxAccelDeg = kwargs['psMaxAccelDeg'] if 'psMaxAccelDeg' in kwargs and kwargs['psMaxAccelDeg'] is not None else None
+ self.betaMaxDeg = kwargs['betaMaxDeg'] if 'betaMaxDeg' in kwargs and kwargs['betaMaxDeg'] is not None else None
+
+ self.freeze_attrs()
+
+class SafetyLimitsVerifier(Freezable):
+ 'given some limits (in a SafetyLimits) and optional low-level controller (in a LowLevelController), verify whether the simulation results are safe.'
+
+ def __init__(self, safety_limits, llc=None):
+ self.safety_limits = safety_limits
+ self.llc = llc
+
+ def verify(self, results):
+ # Determine the number of state variables per tick of the simulation.
+ if self.llc is not None:
+ num_state_vars = len(get_state_names()) + \
+ self.llc.get_num_integrators()
+ else:
+ num_state_vars = len(get_state_names())
+ # Check whether the results are sane.
+ assert (results['states'].size % num_state_vars) == 0, \
+ "Wrong number of state variables."
+
+ # Go through each tick of the simulation and determine whether
+ # the object(s) was (were) in a safe state.
+ for i in range(results['states'].size // num_state_vars):
+ _vt, alpha, beta, _phi, \
+ _theta, _psi, _p, _q, \
+ _r, _pos_n, _pos_e, alt, \
+ _, _, _, _ = results['states'][i]
+ nz = results['Nz_list'][i]
+ ps = results['ps_list'][i]
+
+ if self.safety_limits.altitude is not None:
+ assert self.safety_limits.altitude[0] <= alt <= self.safety_limits.altitude[1], "Altitude ({}) is not within the specified limits ({}, {}).".format(alt, self.safety_limits.altitude[0], self.safety_limits.altitude[1])
+
+ if self.safety_limits.Nz is not None:
+ assert self.safety_limits.Nz[0] <= nz <= self.safety_limits.Nz[1], "Nz ({}) is not within the specified limits ({}, {}).".format(nz, self.safety_limits.Nz[0], self.safety_limits.Nz[1])
+
+ if self.safety_limits.alpha is not None:
+ assert self.safety_limits.alpha[0] <= alpha <= self.safety_limits.alpha[1], "alpha ({}) is not within the specified limits ({}, {}).".format(nz, self.safety_limits.alpha[0], self.safety_limits.alpha[1])
+
+ if self.safety_limits.psMaxAccelDeg is not None:
+ assert ps <= self.safety_limits.psMaxAccelDeg, "Ps is not less than the specified max."
+
+ if self.safety_limits.betaMaxDeg is not None:
+ assert beta <= self.safety_limits.betaMaxDeg, "Beta is not less than the specified max."
diff --git a/demo/F16/aerobench/visualize/anim3d.py b/demo/F16/aerobench/visualize/anim3d.py
new file mode 100644
index 00000000..9042af44
--- /dev/null
+++ b/demo/F16/aerobench/visualize/anim3d.py
@@ -0,0 +1,404 @@
+'''
+3d plotting utilities for aerobench
+'''
+
+import math
+import time
+import os
+import traceback
+
+from scipy.io import loadmat
+
+import numpy as np
+from numpy import rad2deg
+
+# these imports are needed for 3d plotting
+from mpl_toolkits.mplot3d import Axes3D
+from mpl_toolkits.mplot3d.art3d import Line3D, Poly3DCollection
+
+import matplotlib
+import matplotlib.animation as animation
+import matplotlib.pyplot as plt
+
+from aerobench.visualize import plot
+from aerobench.util import StateIndex
+
+def get_script_path(filename=__file__):
+ '''get the path this script'''
+ return os.path.dirname(os.path.realpath(filename))
+
+def make_anim(res, filename, viewsize=1000, viewsize_z=1000, f16_scale=30, trail_pts=60,
+ elev=30, azim=45, skip_frames=None, chase=False, fixed_floor=False,
+ init_extra=None, update_extra=None):
+ '''
+ make a 3d plot of the GCAS maneuver.
+
+ see examples/anim3d folder for examples on usage
+ '''
+
+ plot.init_plot()
+ start = time.time()
+
+ if not isinstance(res, list):
+ res = [res]
+
+ if not isinstance(viewsize, list):
+ viewsize = [viewsize]
+
+ if not isinstance(viewsize_z, list):
+ viewsize_z = [viewsize_z]
+
+ if not isinstance(f16_scale, list):
+ f16_scale = [f16_scale]
+
+ if not isinstance(trail_pts, list):
+ trail_pts = [trail_pts]
+
+ if not isinstance(elev, list):
+ elev = [elev]
+
+ if not isinstance(azim, list):
+ azim = [azim]
+
+ if not isinstance(skip_frames, list):
+ skip_frames = [skip_frames]
+
+ if not isinstance(chase, list):
+ chase = [chase]
+
+ if not isinstance(fixed_floor, list):
+ fixed_floor = [fixed_floor]
+
+ if not isinstance(init_extra, list):
+ init_extra = [init_extra]
+
+ if not isinstance(update_extra, list):
+ update_extra = [update_extra]
+
+ #####
+ # fill in defaults
+ if filename == '':
+ full_plot = False
+ else:
+ full_plot = True
+
+ for i, skip in enumerate(skip_frames):
+ if skip is not None:
+ continue
+
+ if filename == '': # plot to the screen
+ skip_frames[i] = 5
+ elif filename.endswith('.gif'):
+ skip_frames[i] = 2
+ else:
+ skip_frames[i] = 1 # plot every frame
+
+ if filename == '':
+ filename = None
+
+ ##
+ all_times = []
+ all_states = []
+ all_modes = []
+ all_ps_list = []
+ all_Nz_list = []
+
+ for r, skip in zip(res, skip_frames):
+ t = r['times']
+ s = r['states']
+ m = r['modes']
+ ps = r['ps_list']
+ Nz = r['Nz_list']
+
+ t = t[0::skip]
+ s = s[0::skip]
+ m = m[0::skip]
+ ps = ps[0::skip]
+ Nz = Nz[0::skip]
+
+ all_times.append(t)
+ all_states.append(s)
+ all_modes.append(m)
+ all_ps_list.append(ps)
+ all_Nz_list.append(Nz)
+
+ ##
+
+ fig = plt.figure(figsize=(8, 7))
+ ax = fig.add_subplot(111, projection='3d')
+
+ ##
+
+ parent = get_script_path()
+ plane_point_data = os.path.join(parent, 'f-16.mat')
+
+ data = loadmat(plane_point_data)
+ f16_pts = data['V']
+ f16_faces = data['F']
+
+ plane_polys = Poly3DCollection([], color=None if full_plot else 'k')
+ ax.add_collection3d(plane_polys)
+
+ ax.set_xlabel('X [ft]', fontsize=14)
+ ax.set_ylabel('Y [ft]', fontsize=14)
+ ax.set_zlabel('Altitude [ft]', fontsize=14)
+
+ # text
+ fontsize = 14
+ time_text = ax.text2D(0.05, 0.97, "", transform=ax.transAxes, fontsize=fontsize)
+ mode_text = ax.text2D(0.95, 0.97, "", transform=ax.transAxes, fontsize=fontsize, horizontalalignment='right')
+
+ alt_text = ax.text2D(0.05, 0.93, "", transform=ax.transAxes, fontsize=fontsize)
+ v_text = ax.text2D(0.95, 0.93, "", transform=ax.transAxes, fontsize=fontsize, horizontalalignment='right')
+
+ alpha_text = ax.text2D(0.05, 0.89, "", transform=ax.transAxes, fontsize=fontsize)
+ beta_text = ax.text2D(0.95, 0.89, "", transform=ax.transAxes, fontsize=fontsize, horizontalalignment='right')
+
+ nz_text = ax.text2D(0.05, 0.85, "", transform=ax.transAxes, fontsize=fontsize)
+ ps_text = ax.text2D(0.95, 0.85, "", transform=ax.transAxes, fontsize=fontsize, horizontalalignment='right')
+
+ ang_text = ax.text2D(0.5, 0.81, "", transform=ax.transAxes, fontsize=fontsize, horizontalalignment='center')
+
+ trail_line, = ax.plot([], [], [], color='r', lw=2, zorder=50)
+
+ extra_lines = []
+
+ for func in init_extra:
+ if func is not None:
+ extra_lines.append(func(ax))
+ else:
+ extra_lines.append([])
+
+ first_frames = []
+ frames = 0
+
+ for t in all_times:
+ first_frames.append(frames)
+ frames += len(t)
+
+ def anim_func(global_frame):
+ 'updates for the animation frame'
+
+ index = 0
+ first_frame = False
+
+ for i, f in enumerate(first_frames):
+ if global_frame >= f:
+ index = i
+
+ if global_frame == f:
+ first_frame = True
+ break
+
+ frame = global_frame - first_frames[index]
+ states = all_states[index]
+ times = all_times[index]
+ modes = all_modes[index]
+ Nz_list = all_Nz_list[index]
+ ps_list = all_ps_list[index]
+
+ print(f"Frame: {global_frame}/{frames} - Index {index} frame {frame}/{len(times)}")
+
+ speed = states[frame][0]
+ alpha = states[frame][1]
+ beta = states[frame][2]
+ alt = states[frame][11]
+
+ phi = states[frame][StateIndex.PHI]
+ theta = states[frame][StateIndex.THETA]
+ psi = states[frame][StateIndex.PSI]
+
+ dx = states[frame][StateIndex.POS_E]
+ dy = states[frame][StateIndex.POS_N]
+ dz = states[frame][StateIndex.ALT]
+
+ if first_frame:
+ ax.view_init(elev[index], azim[index])
+
+ for i, lines in enumerate(extra_lines):
+ for line in lines:
+ line.set_visible(i == index)
+
+ time_text.set_text('t = {:.2f} sec'.format(times[frame]))
+
+ if chase[index]:
+ ax.view_init(elev[index], rad2deg(-psi) - 90.0)
+
+ colors = ['red', 'blue', 'green', 'magenta']
+
+ mode_names = []
+
+ for mode in modes:
+ if not mode in mode_names:
+ mode_names.append(mode)
+
+ mode = modes[frame]
+ mode_index = modes.index(mode)
+ col = colors[mode_index % len(colors)]
+ mode_text.set_color(col)
+ mode_text.set_text('Mode: {}'.format(mode.capitalize()))
+
+ alt_text.set_text('h = {:.2f} ft'.format(alt))
+ v_text.set_text('V = {:.2f} ft/sec'.format(speed))
+
+ alpha_text.set_text('$\\alpha$ = {:.2f} deg'.format(rad2deg(alpha)))
+ beta_text.set_text('$\\beta$ = {:.2f} deg'.format(rad2deg(beta)))
+
+ nz_text.set_text('$N_z$ = {:.2f} g'.format(Nz_list[frame]))
+ ps_text.set_text('$p_s$ = {:.2f} deg/sec'.format(rad2deg(ps_list[frame])))
+
+ ang_text.set_text('[$\\phi$, $\\theta$, $\\psi$] = [{:.2f}, {:.2f}, {:.2f}] deg'.format(\
+ rad2deg(phi), rad2deg(theta), rad2deg(psi)))
+
+ s = f16_scale[index]
+ s = 25 if s is None else s
+ pts = scale3d(f16_pts, [-s, s, s])
+
+ pts = rotate3d(pts, theta, psi - math.pi/2, -phi)
+
+ size = viewsize[index]
+ size = 1000 if size is None else size
+ minx = dx - size
+ maxx = dx + size
+ miny = dy - size
+ maxy = dy + size
+
+ vz = viewsize_z[index]
+ vz = 1000 if vz is None else vz
+
+ if fixed_floor[index]:
+ minz = 0
+ maxz = vz
+ else:
+ minz = dz - vz
+ maxz = dz + vz
+
+ ax.set_xlim([minx, maxx])
+ ax.set_ylim([miny, maxy])
+ ax.set_zlim([minz, maxz])
+
+ verts = []
+ fc = []
+ ec = []
+ count = 0
+
+ # draw ground
+ if minz <= 0 <= maxz:
+ z = 0
+ verts.append([(minx, miny, z), (maxx, miny, z), (maxx, maxy, z), (minx, maxy, z)])
+ fc.append('0.8')
+ ec.append('0.8')
+
+ # draw f16
+ for face in f16_faces:
+ face_pts = []
+
+ count = count + 1
+
+ if not full_plot and count % 10 != 0:
+ continue
+
+ for findex in face:
+ face_pts.append((pts[findex-1][0] + dx, \
+ pts[findex-1][1] + dy, \
+ pts[findex-1][2] + dz))
+
+ verts.append(face_pts)
+ fc.append('0.2')
+ ec.append('0.2')
+
+ plane_polys.set_verts(verts)
+ plane_polys.set_facecolor(fc)
+ plane_polys.set_edgecolor(ec)
+
+ # do trail
+ t = trail_pts[index]
+ t = 200 if t is None else t
+ trail_len = t // skip_frames[index]
+ start_index = max(0, frame-trail_len)
+
+ pos_xs = [pt[StateIndex.POS_E] for pt in states]
+ pos_ys = [pt[StateIndex.POS_N] for pt in states]
+ pos_zs = [pt[StateIndex.ALT] for pt in states]
+
+ trail_line.set_data(np.asarray(pos_xs[start_index:frame]), np.asarray(pos_ys[start_index:frame]))
+ trail_line.set_3d_properties(np.asarray(pos_zs[start_index:frame]))
+
+ if update_extra[index] is not None:
+ update_extra[index](frame)
+
+ plt.tight_layout()
+
+ interval = 30
+
+ if filename.endswith('.gif'):
+ interval = 60
+
+ anim_obj = animation.FuncAnimation(fig, anim_func, frames, interval=interval, \
+ blit=False, repeat=True)
+
+ if filename is not None:
+
+ if filename.endswith('.gif'):
+ print("\nSaving animation to '{}' using 'imagemagick'...".format(filename))
+ anim_obj.save(filename, dpi=60, writer='imagemagick') # dpi was 80
+ print("Finished saving to {} in {:.1f} sec".format(filename, time.time() - start))
+ else:
+ fps = 40
+ codec = 'libx264'
+
+ print("\nSaving '{}' at {:.2f} fps using ffmpeg with codec '{}'.".format(filename, fps, codec))
+
+ # if this fails do: 'sudo apt-get install ffmpeg'
+ try:
+ extra_args = []
+
+ if codec is not None:
+ extra_args += ['-vcodec', str(codec)]
+
+ anim_obj.save(filename, fps=fps, extra_args=extra_args)
+ print("Finished saving to {} in {:.1f} sec".format(filename, time.time() - start))
+ except AttributeError:
+ traceback.print_exc()
+ print("\nSaving video file failed! Is ffmpeg installed? Can you run 'ffmpeg' in the terminal?")
+ else:
+ plt.show()
+
+def scale3d(pts, scale_list):
+ 'scale a 3d ndarray of points, and return the new ndarray'
+
+ assert len(scale_list) == 3
+
+ rv = np.zeros(pts.shape)
+
+ for i in range(pts.shape[0]):
+ for d in range(3):
+ rv[i, d] = scale_list[d] * pts[i, d]
+
+ return rv
+
+def rotate3d(pts, theta, psi, phi):
+ 'rotates an ndarray of 3d points, returns new list'
+
+ sinTheta = math.sin(theta)
+ cosTheta = math.cos(theta)
+ sinPsi = math.sin(psi)
+ cosPsi = math.cos(psi)
+ sinPhi = math.sin(phi)
+ cosPhi = math.cos(phi)
+
+ transform_matrix = np.array([ \
+ [cosPsi * cosTheta, -sinPsi * cosTheta, sinTheta], \
+ [cosPsi * sinTheta * sinPhi + sinPsi * cosPhi, \
+ -sinPsi * sinTheta * sinPhi + cosPsi * cosPhi, \
+ -cosTheta * sinPhi], \
+ [-cosPsi * sinTheta * cosPhi + sinPsi * sinPhi, \
+ sinPsi * sinTheta * cosPhi + cosPsi * sinPhi, \
+ cosTheta * cosPhi]], dtype=float)
+
+ rv = np.zeros(pts.shape)
+
+ for i in range(pts.shape[0]):
+ rv[i] = np.dot(pts[i], transform_matrix)
+
+ return rv
diff --git a/demo/F16/aerobench/visualize/bak_matplotlib.mlpstyle b/demo/F16/aerobench/visualize/bak_matplotlib.mlpstyle
new file mode 100644
index 00000000..58c7faed
--- /dev/null
+++ b/demo/F16/aerobench/visualize/bak_matplotlib.mlpstyle
@@ -0,0 +1,8 @@
+font.family : serif
+xtick.labelsize : 14
+ytick.labelsize : 14
+
+axes.labelsize : 20
+axes.titlesize : 28
+
+path.simplify : False
diff --git a/demo/F16/aerobench/visualize/f-16.mat b/demo/F16/aerobench/visualize/f-16.mat
new file mode 100644
index 0000000000000000000000000000000000000000..dd42bdfb56ac815ef1c1cb6259d92a7d7b8c484d
GIT binary patch
literal 32222
zcma&McTf{f)HbRDq9Pz5B2pqCBGN>pm#BzH7m;2fDgx3wfshCYh%^-usi8=30xCVB
zNbk~XLQSZFkdTmifA4$m{pb7ZJ7;ERcXnsamf2_bJm;7^cy9dQ;dN#CJJ(GfJePap
z?&%_X{i&yO;2R%*FOBO@b*&9dRqn}N*Y$UC4s>~a-N#$w`g6CS>kooluisU@en(wH
zQB^}x@%r66iptmjFXhz#FQ*Mmxc{q<>zz6kn!TP5GyWPR+za@mxt|Qs4Y}SO<W^Ns
z9J}%9!vbk1+i6(Vf3<RUHSXJ#!+qe3p%22Ix84_R{8cRI_`)U4Y>!FPydLeANxPJv
z5y)!H!nz*U{&H?fn*?3d*%t(ktT<|udc#q!At=ld9fGPMk<?)iE<KNq%KKGzy?pPP
zS`Vy@&H1yXb|u`)>Ni(RvPaa@yLhMIon@7J-%{vIM6Hx9Wcf&Rr5iJ}cti=w_)Q$D
z19ugk$aoBODO?FJX`3E;|K)ciC<dU?`-t#1_Su3U)RpILQ0&m4>>{J($`xN}H`l+T
z+MeV=D$7o^tgMISiG55}gfV1U`HIr_E5%^fuW2_9g#_K#JgH`4F+D@l-w-u5_^z}G
zSwvsuNPk<uY{|_pLKa$~#ImW__m9mG5p+1kBtIcEtV_b2T@`#-e>|7~%RUKtj<2}@
zbzW$Hg#uQtcKCRT%DGwk?e#Su)%+;~2Op|WoY-)3Gyd3O7k9c6UIJ@McMn#Ywa=&Y
z#!~&J7-K&=(rM(npqbPWr{i?Oll;AFF1CX*PAbAeuS8nXMfYm1xhk(Y`>}6(|HP9D
z31;_9epV<v-wB|U*pPL3^7jVZ|Hzi0!t~M3d}i+a=D9CLYMr#^^4t|ntngS&CNO@G
zuxZjHm-9`FALaZJ)x7aN*6jEO#IJTB4qiKo8$3$udSIryR6bZ+;Xr;+`)KfTGPv?>
z*72{XBMQpwjhPHtbZ@4o`l4iwDGjaK_4P~2T(J;o(<4?~^{W1JTIh0-W9O$#?%AE`
zgP&g_q?(;bD>6~QmOR9E*I@ay=`@jyujAhdqQ^rQzc^=&rR-b_x1Ez=G)zV7<jWH{
zi1|LHe01G5ti<+XpgEB8F(?md!%ayH%17A%D9J$uOdEyc(>p!T7xY;H!f9$WbS&Zv
zWev@u6)ADR;@!@nb{c>_qxFu33_!%8RszWOBQpoup{JzK><*EtzUOp0<U-G4(14xC
zPS^E6DfEd905($bMkmFz2}s__gQJ*D?qR7+&^?qJV;u71#Ndb@5eY{tk+~7?K`XDw
zLI`$lG%M|$%)wc%4k;Q3s0)Wa1$ifQz&g?nXk@KHUxh@6ub=jh!t0(C(9CtDkzaiV
zfWz#WfO{Y}#u%jFMDHkR6ELh#7KNsdNXH>pJjjANv9J|;GQUm?aK(tsqw}7L)+b+P
z#MG>W_@2^f=L|jLisp8Rvi7}*=!l@?iOE}*RIAs$c?D}Sus8w%5+d;SzWh7qY<zzs
zMPJJ~`3^feSX0bNU)W~evDYPuJ+7Ltn-Xe3c3#H*dTfTZ*4{!+x@#A1zpXYqA8ZG2
zYalFDd6Cqr+TsztEfTuf)00yJ3DSBRC6*BGY0au0nTh$>uTSR`zv{rN^}UU$Pe?KQ
z>Dho%<T_2Cu9b>&9TeQi+2o>&?oT*8lAmae9ndy2b>AQUnsQHdy)j45Y=jbV6ZY@u
zbv5SF9?3IpKVdQAR#gB3^L~-0G^K<(0`1ej+V?t)=FjEw9T9~fm(^_>x?Ue*6j6|7
zrG|xcCaQOgPtixzWWxOhKQk`W=vqYB;K~M@)0cTtRfq1g&t0&q&k4gnekl8J2Z`*8
ze!hp0Y^<A{l#(9h9>Y!3#oVXUMgRcM(zc?}8<zavp^8Rtkti30eMBMMU$uz3<84UC
zdK;ejsTQ@iH2>4`4=Mt0Z5tjfPNn6lR#!LV#O)Ul7M<LJTBR0i!43|3qU?)S-L0Z6
z9?hPG+xH)JeBC9R4&%IUZ;lPC?;!t!R3_6O%kyb4?F5O6e5kQTF#UY+w3y}nq;#Vn
zA&W@AkvV$mEa;Bg_XrlfHN0vQ<mJ}KkUxrgcPqP+_b+MvY9uLA^YX>^SAS$Q>E_(L
z>_pfs{8;~QI48PmF>(~S{*T2md7|7YPaw$&^lwndZURfB-=EVQBSfx}xTUnAEG6zP
z-d#uGah55-M`OAD_7tm}hhr}eqmS=g`y#sEwe{YJ*F&HtjIBg$6B;$a9Sl`avn6+l
zcIn(Z8XuPUA?gT7Ofdp9tA(kJCEQ%+CiuRO0o-6P<~VLOG#UY11Aop5-7~38>3*=D
zx_0vGX|P(G*!6I$ri%jETsA<--)o8~(cLM!(iNfy9qx&tB{C*WFuxOS2|91u$elgC
z2EdZYk67S?kj{TB@+6ClVq0=yk%!&*WUy*j(_@89yM29yICm3^Z6`ull+(Sk{r4xS
zPWQr2bUiYyHiVj+eie@|nRr~&c5u@xj8m<0BZU4?PFfbt=N~2pJfOUhJ0A!3%Pe6_
zPDsML53%?9)L}mrif!cgeJboY2}vINRG50jNeS#DrnZh%ZKYm(Q@slCzt=WkADH<m
zypO*`8EzB8=fNcw#FA5OP{iI$yPhVB_0Mz&%iGcOcee?DQ9Azfc<5>Vk@2-Zns@ru
zE6)L40)IUH=BD$%C%<dy{`9&2r7QfGFAn$vI0|CKL~n>bS?Nx^$?KE&QQ)fLSuI}i
zTVF!;x;cJ7GQ9Ql)v3W~IkQWKCW0dO@87v}>GJG`Z@QzyrJCz`;^U`kT%V0xfr<0T
zW&<6MHiJkrD8h_HMo2(tIZ~O?;WsF`ICa)Y1IzGQ`=y+Fs-<rxUw)M|GNJ|87peRO
z!R<o&`*sz@uHcfU!g=atoj;vuHrxYlBrR+1E&Lgom<|^Og4jW39vV!fA7duMgY@qj
zgDGSkzSFF{upCgQJV<iXuZIABH}F%86DtODWK&gUS1D#%q>{v6oiWKNIa2crFVhXT
zGGci-S{+e)5OO@C@ijQz(!L>U?Z8p<%BJf43(sHcVT>0_t^U3dueUsf9Ii}5UPMW)
zlNZh5bv-*tqY-5`K8eFKm#<B`*{)5}I+m}{lMV#*%-)vBUV2rdUcdY5uAWV+U%xL&
zO(rQ|8J;qV!5kivPblxh6<_{2X*7NL0rMEHV;?wk3I6u%ErFE)mgfrv#X_qEub^v$
zp+~b5s(Vv0GUAa)!Lw>L#>+CEdFQ0$OQ|NsGwePmYA&vxKTacP!*8}4gD5t^#T^sY
z3#WvhJ8GazmN-0pJROvDdV3wYOQw}wr)ShK^b;}!Dhfz}*VNR{0NvZ?tK*`nKeg#8
zl!?ZBsH`Bun{Ut(raAChT#py7^*K%s$|sj0s<qZ1UHVw5I2}E|+F5n}rm>kE?w_GV
zaf3E29B4u0Q`8<1dH5=~cSos>vUlEmiZJC8CK?`O^$7WC{~8Jm)RBj$_ja#s?i76M
zDLT2u^>!<S0w5Z>?++HM3)X9zn!CJ=jx*IP!`O-S1bmU~b-iJ`mRV$*Q2fT`w5f+M
z&%A5f#8xp)BZ1pgsrC(yESL@!J$#fXm#XF2OU-&!b6-_xBJJgOdy^xNsZ^1So|jtd
z0e%y=RJv7a(^7RLt;pk@=jGF`-0W<t{`v3!*aV-Q`2eIvpOX0O%z77NnEpV2=54=p
z+mGazEx%;flHVwzy8y>N`}rwPd+SCuGWYI;;XIEv<bqIEf4Z<lN1P^A3)@*G>zJ|K
zzb^cR@O0*9<GrlcuF5sMZ))vTbNa_`*WK)MpX||++?tk9yuVQRg?yH_wW(m0O7KA4
zPBj61`q9`UAkbFck3B?wa4etvu~PG7qPhQV9kxPR8`Et4DaJy)^NmVrUp;l8Fa?m-
z)%0a|;cEtB_yL#S#$PLvOJ4pH!`<O;VYaAGyoC=&kPFA@JRF$se9;N7{|1eRq`0=A
z#O%MdhHhs$oMH&BKkeuzR}|a6j^MZsc(@a-tpO<oFP%Z&U2PD%rF@_-xn6eN=?LUS
z_)sCPSZlYC+_EH*zo42lCsHY{8|R7GeEOX0aQ%iN<O6-&PI}Vbc%8Tu0)#xj|DK%j
z9%H=|mej0vE9;3&%A<76^wU0uJqhg~$#G>KnAi34cVW-h%kWBZRA=U9Q7bkf;K&$}
z&Pz<{=I~iiJ-k$OG=4hW7~%A#4SvIL(YuqLc4yh$QNr$9CYIpjR1{*am;JcnRgUi0
z_UNefL1Bd}XVcCkpXmvoU3T^p8*qr?+yTrJ&GI1<nIZ21xf%-cNgG6W+l{9sA0G*q
z{^STLJJos{ZGN66{T=>vMs7NtkQ<R)bplSo*FAO$GsWq(*SrX_5PnGa?*RbxgTt8B
z9@dO^8VT<J>L<!?!f;qSHi9BiL2~C|Aq96~WD2W!+V$+3+1RM#m6YX*rzKx;pmx3|
zXKr%ss-m<x387gVMG@J)9+wS|6izQ~uppj4q)UrN&OCn@bbqIiQV`ZWvlK!+T2s$!
zxSWLGC43nQ<t2}>3peB=93Jk>%-enrFrZ#oi<_T$MI7E!alHObfqQpn<@2MLFno5s
zs#-L$$``ctcY8kh_rSl3*WKIB8f*;Ld!agNO1AlO8iuV^x+5EDxiJXS>w{<Izboe7
zj(puzE|lc+H+99eL}3YiJyK2o-X;sVR_9&lxmAhwghuFtcS5JR&mBmhSx%3h3ysA|
zo+F>Nk3~E25Y9*GM1tj9F=GnnD3>09>X`XR9)cu!=|LVlYuNq&DWXqG-)VgPWF2$A
zrak2vzVy<N9BTl_@$J5vhj&nyhOAkS_&ZuogUfkZrO?cGD^Xxp=!JIBz(I8Nr6E1m
zqw^h*M|bpBkp4XqL;Qb<5)08>WsG;2=d6^6c^s@7&Jm{~lJjtpB@V${^n^Tk{zxDh
z%K88%j<sf`b7=9`hrcJC?a+LWZhP6<72dra6LQ&ia18jrkJ@`%4`0hE{%LWWN5buA
zs`O%ur*E;1*ZzK&#>po0MA^mdYwtO=KsTPB@9x=r`d0AadCUbaIk_uWZ(Wh&wg1bW
zB69AM#3v_SUVQU9!XMN`+FUleuJvo;GcjW4BF40sG@2W*KS!->n{Gf5_kz9wBeoV>
zb(J^$H_akOQWCyPr+C{iCkDkF*iA}=NV3~nQ?B2HHOWWleA`W;%SNx0Fs*FPrhbji
z5@Grs08}msdg=VTDF{7AXTfg|zpJ+X7)GpAg3XaPY8GXg_UTzONw+W%&AOSGbrEpf
zFvS`t<HSe6(rp)Vse6DZz!pWazq;*OtCZiK(%vC@dir_0mh#GYjgI-OXKV_^1IbhS
zLC|)rY_M@y>}c^Ldu-cN!Fy6>kQ0x|!$&ier(uh*z_t@3Yhh_S;d_HM`&8Lot)k!|
z`{(9A^FU*hda7Pr7qt_JLc_$OxhGv8Q#aS`S}*`Mco~q)q#x-g!|KST%L>vA+P`H%
z&St=Jdc;g7C1j>rqY>A*3=|N)tA;@Id;sXg&`X93;&6=V9}lbTe;f~lr(iWsB3wrr
z+jc)S?F}@~_+;C?tnNyo@-n~9X~I;PC+4q`Z8Bi<`?s|%$^D4iS(WlaOY6k>-hFke
zyAYG;s&eM?vXTt#&!8MKap~8LK62J|b7}ex=4tx{`G*<<N5dN4U@AqRLEwGR*At$W
zwdlEkkdS1XTGiq;0Lbnby=~2_@sVj4bgzbp1=fIilA)_zk|`k5g^VoEnY5bMRk*d|
zlMRl#U_;WB5xZQ#<2kyXNV+j(tvF-Jzh)tLov1S{MPjBiPj*jCy!O<04iHqBQ_IXs
z=+SE*h)v}hc>hFrCUKFuOjr$G5+G^rY%6%6=JB&~Sqoo+4?>btT0CqA#1$@~Zsl^c
zk3y3a4&upSvc9k5({~(z8P;#57zeo_qM`xw;>gO9w<o0o5ag=vsz6=T#2cALt_5lh
zL;Oz<Km<$dc>(IYF#2mgCF{@fM*k?5vFNt}bv0#Lzj@s7HR46+`d+?=xD~lQg*0Wm
zE=|e@-={^_yzoiQQQSf*7|R{oADr2}&`~q|aYjKt8Aik(f8cle@QTT0%DpSxf~SuN
z{+>xPOI}!}=0#vMr5$~`nm09<g;yHTX5@;8dl)tF)SnOmIr;j}bz2su80)pB`UCy@
zdjin&9U$jffd^WJ-vt^R!%11wg|588zvdkk0ARPS4~dkyJ(yDh@=GQr5g3XpkiP}}
zAa(Ga0@@7cQ3b9Fz6WM@C^;W|i*}9I=h&1y%>6Z7qaYma?jR3cl<?E1c7qUHdl+u5
zWN<Pw26g5*zn0A>QIzitA}ey|+`j%{?k#*C+sXGwxthllD7}ic7AR0Hf%wI08+6t0
z=vA9tcfx|xs@{QJVW<?~66UYLR$TY!K1$$pb>W|IpqUG<NpaDqmW4{H8W}X23tan`
z@=$l;ctJb3fsnH^E7O`ih1nM&v_waHbZvv?sW^+Sp>eh*GPysyeBOO2ITBeeyREb;
zG8uY)MB=!<5Cnhz=lVxligZn1`T1Q<7lHY4zCV|L?Aa^nJK5B!6@r^2><Q1NJ_a0A
zbV*f)uLy0*wKmd?<3h>ZOQ`^+dFAE9%4ba8-oK}hOam=1HztYstQ~J3Q`L?Q7bk)-
zj0Jl58@13}auz%HsLo{Ui6B`$Suu~C7Ib@&5mT*MKocW-c2ae}3ruSjf+SHQ;cTPl
zt5oBirYtl~f_5cg=X6FE01)$rV-IbmQYG^tj#;0JOL!7)dfL8U#w;0@FI0uiu)|&i
ztcpH^80e@4KQkM1QWec@Y6#zwc;&yBM)hhB$-;}DSZuO?Ey<S&*p?(@&C3zT*mFIt
zmZiAdY8IgVEK;}+31ao2ij1e<`^PmG<0i)gH13+No3A9C!neFaNWl6mtMJLQbBUpk
zro8v#)Lb3@!kQpb;M508T5;g>>81x=j4%#*?kg!m<I>S)9<9nNsGR*fY~68QnkZT9
zcOBOszDh=w$^1%XMED2hG0DRM?EpcMZov5s_i46m*?F_p#-o2BE3s4k&yy&&tD@6d
z>wRjxW@_B-IXqC$2Xx~({j<iJN<UOMMTf*e^<!y**zKc%xXGCfP*VWZ=ynq6w;RHo
znXr`ffhv7KbXxEB-U#SF>Z)!8zo;?WxwPjB={uTKnmjMu1jAnd8{mtlpbUEd0tHi)
zkDIUy5d`1~<sj|72&iA`bic<=8#{5VB<y`qoL}C^e{__6kR7hk3@^KLntq$=8KYhG
zQKCnG+x%a!j|R=)re!|GDwOuF?VHa>#c!uTdg`Gqycv#num4##q|X;q)4-y`n&V~8
zztd;deGXltx4Dy7N0nyb7NEay^C|)VoTtk#a`CNwv^3iJx2-nptgz(jJng6k!&%!h
zVhngwN8lXXdb_#PA{nX+d_+9ALs#XlrUi|Kw=+(rwKo_<cIP;K@P}rn&_9Bh`0v{=
z!L<QhNU}z*R!_i{!~V~~A9OTwEpU3=*=KNI<<U^l)DqnfW`()T8=^vcyD_AoGTWd1
z&qxCrQ`_wQ_18BoP#TAUy=kqc5eltBldU{Z?TVmO@FS<)56p4bXWAVR=A(V#(cAi~
zi)~z^4U=#!_<@6w$o<8;fDhcOsnn)p$kUpzFU=iG>75A9TzypK$cx-az!>mj0rzF4
zpIz<9d1N%4bNB9n>gGw%`<%0kkSZa4=Iuz*&5ORyXoWqW4xh7XiU;uHxmjx#?=e8a
zRWv(p=|cqTXSmqj+Ub&mRDG(`2c=`4Y49$`7v^lrBd4oqw&O5P{JEVe><2BSK8t9r
zDIIbA^Swt<p@XPn^$v~?380zN3Re%}ds3il)HygCFdt=CMz=l4+85ArYmfer+WQXX
z=$CZuplX$k_K~F8j^sZ}ONM&zt@3t;na`$<VdBjuIUTu!^q8&==CF&~ImeE30xqyr
ztvHSIt&;r6Qi=AAiQYF%e&Q8B5XaQn}G%5j0W2eg3Uq4PlCDwdU~8TdXj`OcVr
z`-*dPSTdd811UgZgL*swXjS?}LCsqq_6lMFEZVMlm`PKt6bt=5{pKfl6-7AkD(#_D
z{K`q`reI(WB?%hk{TIQ~ZdrI-n=dq8g7V5UDz&&gI#xE%#8C4?YM}%7ao4bpXl)MY
z;5BT>5^#+zqDlfa*ecF^2p7Tzl%TjpF&_dr@U2`Ym+ohu8GsRhBmD<PJ?Pe&!^=#k
zdsx0565A#RDKem<<76QZmk2rMmebxehHmxnWekJ44`YuqI+~sx31d`eGsbx^w_ByS
z7H`z*9A=w4Nml!s4z5VzWx-!j&Sa5~M8Z7oxhT2iAzYVfa1u{a$Lp+_F;<j!xG0<Q
zqHo!CUVbp+AZxh$?oUPVrKl;$`1mo^K=n8d<kW43FV_styV8}$c)WC9XRVK^?yp0^
ze3XGmjfJj-zfb<luD=?~Z1G!%j%edykd%H~4AZfYyFh}T0KE`e#Rks0i!HfIuy7z2
zP3QR(YCXulYwD1G_%cEl$VKeqMpEy2w?0NW5KkQ@;qLtez8LH6)ZuFnxNwrGf_=3W
zVl_V+76;A1-#IjBY2qgEnOlofw2lCDYmqH(ofJeqG`(o4Qf6Ifn|pp(t?A@SgdtW`
z-}Nc(Fk;<fpPlZIq(s4l>O^9`lusOJ8|-zF?)?VPZHbzGy3}#F-qLN$wpi%TX?Oh9
zp!b@%abPdS)^SGrA4nJspyPzI^O#;^Ocg9l#)tl2kcG)yw>=3Tnz>eTZRT0^AMBf!
z+PEk7RgzV$U&ZYTRNc&tDoQW$U!L;cGx}WqTgY&?w#a1d8S&$%Cvq7E53d-!+0<0B
zl~gkL{jDO`egxBc^c_rN%K?o9AY8$yBX9Vq`qAck;_i4*s&wN_gULg~yKWuf4bTK*
z*}Ns7Nv8kosHtF+g5+yC|9`*C1adNzTRuRwsRT`F^X_}wRi}T)a>|MA)Cx^JD=c?=
z)RqPh&bHo|AbIa7D;B9{K6_}OpJJ?i7$C@YFMPWW31e$I(xKNBW@Xu_s5u-@el^Sh
zLcqqmLvzw!Eoe(Bg$w*2v%%7`<?%5$7ny>SN|ld?$Go@*9JM|?j;+2pYvUd3Hs2fb
zPR2Pn=0NG5#afSwiF8TR3TRr;dd9>nAnq8+zx}!4jniBFYI}eh$%Qg3dfBEK@eYBq
zh&UX%6CFn2rq+c3*d1&6+1D|eTfrNfJ;=B#E8(@-paw$)P)i~;*+aq6d_#%pd^D+^
zu#$e4vbhtE^e+!)V$^#V*HZzrd{SdWGYia}oP5_ZTxoyxL4LYlE~K*FJRM!2y;QLG
z%RQ%lD)y3;-XopxApe(Q!#0yeRM#!N*@mkag`)YS%tY^lP<1-^D9N%Rq#28@Q_@z6
z`Jgooy=^F{Bwk<D>%G=Pd#e7YI+9%51rUYhX!v8SI%maq?OB28ov@iZ`vst5BIp{%
z`aJ#Fx`T#j4<dk6|5@KtT)(&b8;*S%6})u5g<p*8BO1v&)N&TsjQ#Esu4%Dgcx_ht
z<hgGkZwzwYow*ZnH{VP&<L(uLp^rzjgJzZEU0hn4jPTD2K8*Zj@a<*8JEH-)%U#JM
zB7qBXjJHa!=d+z6M#egB9F463UGFn2%WaumUo_5Jk$dy5F+ae!RR4*A%U0_?xdlbP
z`D}9A1tw=UXid4&4{3A4m5df57D{RUhCW76o!wOF;F?(YG4=0p7kiu4m%lkBM9A8n
zgzHKFJ8vmFa!LI^V6=c|)yZu0C*fq&Snb~Nj`~GdZCyi58o)EL@o$Lc=A7V3Px-!@
z%@lXZuiH)+1TilsB_A{D2QVIqtt&!B$9vTuFAT|GAT#LL?wzH?u=edhOGrcmtZ==@
zG7<86T_O&KAJ{U8x12>@2ulzDZEAIhef<~6{znHe@15H*8c2m&$L%`AMu?66z3nYu
z`={cqs{!a_<stx8o8J*h->Ssl%EEE0{`n}p0n^;ne8lYqNb8kMqzslK)k^F~sz$Ta
zx?JpCYhI?+@P7yOjfT)h``bNPP#kFP`jqAEaPbPYdtHDzUiJ^FkC!{ExrRbIy<=ue
z9?=8?uq8s&<8T7v2!HVX4MP;>-e%r#!|5x^>!-fFfX+2_`sT>Jar<6XJm7fv5zT=Z
zZ||b7Nn3d_DTK8?M@|1ZQ14-bTb<OOH<qdJ4fuOR{6Zlf&f52;l@D&X<$M7QAgv|6
zTZxZFA`yZ%>CpjpzTVON%M0MY?@V^sQP9Ta(v-yJQMz5?$c2I6qu)~UDV2)PdiD5O
zN?u_6<5l?;VrOxy0>VhuUF#Y(cze0LlmVoW@c_Yp;*8GbB9BWSmg1W~Sw<ZHZAebd
zo~EWu8~Y~Y1XZl*VNQo8&)V-q<Of%GiDcfS+@|%Xhd#q|!kbxzFNK$DY{s)pdI8#g
z!Jj*1{6%I0+?IcUGd88Jq4HG^rXP)g4mO=>#tqkyFRalX=MQRFwGJ{5)RN#o8_+zv
z^8PEXQSivMlc<oNQog+xPR{yHonol9OXPgVvHIM4?|30bRZNtPbSdQ};BO%|uJOup
zEacOX2O<xA3IE_g2+R&#zl1lXMYl?<=N#~@Ym;r74({%mt!t0?L?P<W5p`*sRSMh5
zWcDC^GUs;IEDMp7)_0Z}7Me9JK(q{9iDvv<X>C%2WnvRF<g^5cj{9t-I9~ly+$sQl
zJ~T=2s6K}MEu#Q@{ebmIJB(%YhSNJ&vePN`H{(K$w}Ph+S+$5wqoF&Tv~P{<&@P76
zN`Jf0N=>_udxyk(c%<rwT#aM`9qk;pL;Z{}d&k}1_5^QvsPWh6Y-F!NHkLd|<Iyr%
zeSj-myF$|+sQ}kc3(yQ#Ss;n+6HmEs-!LU;KgS@yRp8vwsl1~artL=OA!+Tew(3a@
ztTq}<`|Z6=vCoctEyeG8i}5L&`++LAE#`IUo&iFl|J2fY{i0QxbIZK=7d4i^;aF;!
z3c4amA}aVRENJn?0hsG!Bai*?>XsC$U&+hO?&3+emUQq^)2D!<PFQSBDe|Wa(ccQv
zn%FWVSYNyUSm2-_%7ZhBT<98Sxbj;{q6f}JX*4{r`wtJZGXfN+yEjdg<^l7#%~#XC
z`y(kbYznvdw4Xz+94+s|A-2?Fu~z-H=8wH+upS%OPo39NJe$ry2G8Kk*tl{9T6C@o
zOqG{40=XHf{oBiF><k}%s{3-;UJ)3FPgx#>D=nPK_0?Oq3P)SZ(yW(YgLyEvHJzIf
z<0HYT5Dxel6)WPWa3*EoJp+su%cWl2kN@{pG&)>n-=xW5)SnXS(k%&SQ--y-U4`^Z
zwm#nQ;zuqmxE=TWxxLa|p8So1t-{Ok`=zlsy-W*F-;AN?6rMEJ!7)ySgvJX_-n`Hs
zYmcWUo~>wai{d)fwMcS7$+r(eRcZFOsSyz_A&Ux03h#Q~IJfbGU-W8scD1zHTdEp%
z;2atOp}WPnpiJ)VX-%KNiq-KfOxU1uyKBHn>Ti2`dd<PC{#tU6q$BbLh`?EvI?GPm
zy6p^{e2_Y6yVS|RCl{im13^eBEcl-_*sL$4gM+NQtwdwfdCKJw!UEI^TA(^|9QTGN
zbA;}m6>aL~`XDq(eYTPn<^b<%(bG-~a0|Z<eWz50mb7h85bkSun9Pq839bp@DLGK7
zitABW7ZhV^k-5X&oZ?&gohreZP~1P9Y^jtE&h*--D4|4m8TnNIT+pM?JKNt6E)NU+
zb}Eg=%NKz;w{L{@gR!acA+=iQN`Ai=e4<5bSAqn!HZG8L%lp_E?*!8OLzdn;oom(A
zvN{JiCC_(owN=_|pW9AB1eM8Tsd+Ifdysd8&n+Xe=62|?bxOQTXBC!2%btdvZ!~WA
zIFrnG0B`3w2Z?N#SQkc%WTu=m-9~*$4q90BB5{g%W^0eYW(+?PscRLC)L;o8)qnHZ
zD>LjTB6(9yu(R2jbh)+zN^Sb)0@uHK&8j^SmOtdg@EX!GQw%S-mPkaO7ZRi@9lQ3q
z_K^fhPPjJYK*JgN@YjU*yKQGldpW{N>w%|9Owds+djKc2iiWKesgf1;RnY)b9ix>5
zL}!&e6Jy{Uj`xOZ@^kFD-ey^7EOo}60dc-1WG}fchPROz%Xo=GSBU1GuTD+b(bCsz
z^_<E6vRWA69r&y8v(q%!#*S}4m9t88K9BaA1;SMS#C1qXWvK`(4zGRUDCGMScqsBe
z^8&8(t@HJqUwcpdPwD5W|Fzbs|GN^Dz0R)gPwoyVE0t&0=%+#lgO&RK_Auwy&3TMA
zHOW=211gYd^dUT?hM$%W3O?MNNZ{afTG}e4*cRm?L}!lDe=&osH*T|#eGZNy?5ECW
zdmhXliPTk_9uSC)x!g_fbq)rv4@p&~|Ic7CcpmoyU=-rt%vu;?M9C5b%~x50D(*zB
zxor=+R5(`fmgNjGoC%Ska7N@dOIQtf@ABSx+whl<TL99Y6X&dF7GDv%vKBFad*iS+
ztLTlS;vOCoc(!Phztry!KEe3m*ai2Ow6QVnOrJ|lbP%?wFt}?k5l2gN`Im6qr`Dsv
zl?Xpf0_`jIATBojr)^imP6&ST(0k=`d9@vW&62xDpY+FcclyWRe)6%7?`w*H5qI9e
z_UuX+#SssY7*thFffX_}b5?O{x^vcgxYgb2QAxE~)=67LPUTv9Mw4G*DSFcR4TDZu
zQ%(Ihq~&YS_OnZ9Xt&QV*%R75WNxloyg&Q;&3;RDcgxg-EuHyh^Aq)c&Pa$%CSq|j
zF++c-;q3e0PcTrEEW{GtRSCeI`F(7qIm24Pe#r;LtAw&|3#`tF&hqh!sQ#95v^`@_
zwE7l5`q);%4DPo2b;w*DnqWfPRZ2E35z16^j{GmJLcl~T%$`(tfW4RdHoi$S#Qhq#
zp$&Net({IL6Z6In;U$N4f<z7EFCuS^$Rz>3->ONN-*1=JYNON1%R8?gd0%4&md~Z`
zDLC<JzCNBi?oNm;K^0N%y3U*PHE?CLVMmdKlmCzPw*lLH-MEs%8e+_=x_C}$@S%mp
z=YQ6G0Y>?g_Y8^#EhqjyGWl3t{9*mum9X08YsJO!LN}d1```DC%@GV{eI0vx?sZna
zL%jzALmHV_?8S@#QA#XID0SqK{RlNjK&Qgcd%G03ekOg3Q)AngI#zcSq%r>jhMap=
zU3e;$ls#S&NN`$;)8yOSaT{s#FWEPryVt5~7C^PMUEgR8F@E^?rb(@;pnt6EKR))y
zqRE1(55%Osl@$K;=Pd!w9{b(<VX?258i1S(yf&i+fxdXIp?!Q%RH!I_9#PE~^y~Op
zd#yESe0(N9N9PZqu+W-RQ|*VVDO5@Lc17DBZekicoj-H@+G52yjs#Jvzm-SXkN=i!
zeEWzhR<;&%@dC`y{?6TdIZ)ck=+_`+_scEaTh(o5*Oo?AyD{O#&x*6^aeo&t+xZ=>
zzuGu0O;|z5bsx7q>+`m<=mr2iEVcgxCH>*37(l$J@);xNeF9#n3d(h*Jyr01wA9W=
zWU-!NUrYE8tSM)fvD}S?X6O{CB0C`ooCh0XQ+uJ)&Ju7&{?E}z8{|vucc-xJaqTfB
zqMY{d#UO%Ci<wzuWPp`O@|eYfJZ>sDacY9z11-!on4b6wrEhZ5EX`UvgH6<V%B4n6
zqqq#+Yf5AyV`<8Rf$RTT{IM>sB61&v?SbJTD#IU_Xmz89vKQ=k)bK?<^=@4g(LF&r
zIeI`z0`7Q`r4`>QPD(88dq%8exE0KNte7Ko69V+l!1{_Cr=2Ih?e-O0Zs_l``wSP)
znlz*BB?5m8+uQ8^)YC6=6b?ZgzJRV&<7T}$r!5wLy3~#~Y`b?n?z-eNBYX|!J&nu?
z!(3`Ro`J1Db@5BXy*r~#pWY6Xk+wI(<z9Vj!&_W3u+g9jhv?O8hsuytYdWhkeyN2R
z^<^mvPD(PmYvKroC3+JA$+L8qYg;O8DGogFhvW2Me$M5u4_p4|{N6PE8LJ`Bgsh^U
za#{{)x&*d>Gm;&Ewp)O&kE9W<D%BUwe$WFFG-`4~Gk8bj2A%4}BaeZ>KEYqD#T`MH
zj%?cdjh;_P?Dr1`J$i8XS-xvP{)<OfNfy-Z!6dk2^8K{{EtCD@U10VX-dp&DXn_%a
zsb>y0J_0_}8=1bCj0}4&Df0jW>xt6(e!`SV%Hmk8OyF|O%h45C&SA4D6&;-?M(mFh
zrBTzRMTJr`Y$AciML+L-{2p>M?VgN9;e3ig3Eh$?u`gm)5*a72DH5TSbhGt!8Ey0M
zIsT#<Gw`;%d%r0Bvwn!18~PwVc_TfwA7cV+nF*J&7%59*8R9Xx9Ne=rmJ@6cDDS-7
zOfffD_r0+=;r-@r&fX0zac$*$1<j<TgVEpj{KI`$4{OWrVUF~Yb>!WD)xvx4)@*D|
z9!{6~hQB)Ns4+dVLM8no?0(W=3kY-DDiJk8YJRWCcb3={O?;+X+r>rngy`$KlH}$#
zkxzAFmaZHxHrnZ?!%&(KNNnu9oIkhSwE{^z6Wdc4ApJb-Mc9Ig?mui_s|vN#;h!|L
zCtfgfMQ~C7toqCpmnTCu^;@vW)5$tDqa}b^9v3<?>>bDU8J$vYwR}@_OjxAQ_E~5t
zK&^lq9S?mcLpg)Ma;yD!f0Ce)nv}B)v^JS5HyVXjCSMedMxy1&7jmQEXi4%p(I^;N
zBJ=_yY6h(tcq-z)>KI?(S^9(YF`2-N5xOJ})c0jZRnuet-HO32#J?wetqZ?$uDIv1
zO-LM578RfBB7&cRZ3S(OyB2qjoeJLt*)B&$aJb_}R>oUQ^e%^lf@Unt)dMvQfU5g5
z4*uXfX~75$MoY32CE_Vu>Lp)kLvKC~uFz~eICIg*{!*h6tH&zS-Oxcnr7e7c_fshk
zE#&ZEO%%L3`=3Le8`Q^-4>5W)gb+9jo3WF99Bu_5&3f^w?iiz#m90LUc{Tj0uFg(u
zr-q-Mdw7?3ZvMJ|4dSx>)qm|w)D7MO+)^^pgEu;yLM2XIvqXj19h1mjcIr*shVbLv
z*ox;>j;+^ye07?F8(SC#Qrc^ujU5z7>DN^v7XAc~L)$)~6fQFoV4jCN&&7I7hz$QP
z)=7Os<eea0*$DF6bq}_FO5$X(<B$eBh-J;HlhnVaj&+|4@EHyJtrEZ6RS3MgOOTY)
zt_xn{niySeg;qB5e(#`fa4<ynD?ZlzUJPjqJ&ZA1TS~ayh`gX|9UoA(P0FfahHV)f
zrx}7yZd?<|My!p(c<{`-Ck*ggd90>*>-JW}TE!#3;mE!VDLXIMyVB#=g-s^&!%=#U
zIf*xrpz=@&OnFHPl03av@1+se_(zTQk($r|)WqKxv_pWZ$-B`ZgE|*(U)EAUpp<iU
zo5*X^i8&$3nW5D;m0ja~q(wBm8#YQqxYn|txStQ+=iux54l0Y<;k^rxaVJWN_3uO*
zCq|nOb0A;&iKsRuScq@l_^vS25O}DEjZJUxskkv3&)m{POh3}4kcL76Dqc7khRUlG
zhO&_=QLn>R;<M+{991j8M>r1k2E$IAb8e~rYL3y)O(1EK)?__Ndky0eZ=sU%p#ew%
ze=f@*q~UZ>OlKTXAMW#H%G1v*e@OfCQFMr~ORH7HsBL(%pIQg?%Wnl*^7@G~6m*|^
zNUgd^F@b8}%#fTI%FkpST{8F^cJqKu<1<P1(|Um7$uN?sfWkt$g|MtT(I@S_KNuHV
zF;~e{mwVK-eDl;cL$3qmfR-!+j$6^LCPz#rP$Jp4&vKOlLT6VK!KZzOcmg6Wo3*PR
zuxYU>97Y@nwX%)^irBH7lxWd>Gwjr~f@cIrKoqha<w@snid0_>@8cnN7t+q@oCQTG
zbi&aoDR^`8g+T75haqGbmObpD_842>8K-CrT08KaPUL@lJF-K6NCd(LTWL55U+474
zGB^1BFQ3FOKGKm$|J>W>yk6>s5QlvpY!^9+(|Ra>U22RSBBa9wiW_|eXwQ8|KOG*k
zr}MI<W_DkxU7#ri((ZI>={o#SdtC)0(DczkcceprfdG9!*Tf2cSbg2~6UTwYB0H3u
z88<plQAMyb9S#d#{kyE7gjBUS*j3P^gxswdP05&v+l>A7bS5(2g5qXh+C)f&%V=jx
z`+21*E>0TzS930!@sZ3&5kb-P^4NdXgv=Vwzf+lJUjWnEM)>@%PWEKz#BP(z1y`&<
zz~|PM0m`j6tshUIdP5lTFOXYLk9xsZ5P?fNFNUiM=e7r+`GKX71XPNHwsTJ_R%V;&
zf-)G0L)RQsC30EIWR-mW=LaWO0Wt^V6;lWh;7+(LC@PGJkx6dLUG}w&0H~jWeS*~K
zAZPilefL72U&CkbTPc|NRxKGxn2r*df9dQSW)A7-FjkUOBX@@y9%^x8)|z!k??kei
z!Cu-m?e*p^%|9UMZC*yFX%%Jfz6XbipW}12%F>KOpMoD8;}!A@yV2Il6`k8tw0xs=
zuFXWqUk6U)f{%}z>{In~z!=z|mPS7!t<SpcQ-HLUj(EiQbpnqET$9h`CP9p<l1gjM
z2wE%+l>^)9&ZfQwOw&(;-&4F!&O3?Y)M_SM??}ysh=Yg^_qp{}iKKvpD?xSB3VtdU
zhk2zOUay9KLBp4}hWFoO5~gF&wjq^7Bl^+{c<Bn&{t!0Fh^EZMCW=yPvL9L^v#oVT
zB{N68*q_rY1=3{x5RJ4r#fgWmRx_%Z=OZff4Sa<0Dic39S}rVqU+6?k8aFSR@CBSz
z{VDqPBqFUGbZq&y9tLs|5vcd(n8RIog3mieET^BQIk!nO`sh97?D6_dms*;Sxv?J3
ztyLOVAwm0-3lnzVNhZ7Ds#il~7%-qUuo;R@(H?*Gcso<7M1QrMEbsSc+lFGC@rj<E
zFyY+4Cf;N3HU#Jn&(4wk)$T`?Qt-Sp5uZR$RD{Z9^O;mm2pQHK_^b<tiU-*g9KV{*
z=oBIvf;5Upgl6GS6=<<^(2C4X-KFibW)2_^q$^?nN>BH;g8yi#6}VT+qETGG-D6V{
zlm6RjE7jF+`59CiIBEs?rf$_k6nZ~o!oj&4_6LD6lu-yaYYctk7s!Eu@Xc78`NJ!%
z!cUWi!{JlRpU8IoPaf#nNtepb43epH`s-@>^Zn1s1y1CK%E8W1yPc(41T{-2ppx7P
zd{<MsU{~Xn0|BmTn{07~%i>dx;$LoZhQ8DtWBc9k3K<Q0(#rX8+C>Z#E2fVV`F)Aj
zrL#Ay*>kr7G&`N3ZN%_y%GBIu&HWC9^UYf0pXaXKK03KGw|Q&Faw~~AOm5)aj`Mke
zi`xHDgL>gmk%3Ix4S5_R1+iIuRb{<tA-bbvmLmZ+>qXS$w&Zi`^w#=dCv!(g-AB^n
zVMA3d#riTq(3VPX&#IV+kph1JG_H01C~jWAmK53WN3GC1K29065&{ZJFi_r4vc@Qm
z*z1H`Wy1|W+^?-l^r^se+ahfbS%U^nY<6hWBNLhGQX@PtJwQi4Cshzg)Zz;l_kTPc
z{&`eUJ&Gh6wEPn8(rORJdAlI%N+Il^*NB-zNLBOSqjQrf5xT4HRV@pbQWC1Yr!kU^
zFH`h@vjMentF(_8Z8Y~<9z;0e5b-yx^cl`?KTEZ&fSLHj1KB6J`}r<LXp+&5c@$J#
zH^m+r16~nMqT6jhF09U`Zqr7?Rp2BeRygvJZE8xu@Um62>h_)leaTmOqKMC;AveHn
zHEes%th{O2(2fWF3=e(DMunSiT*C1W2pqPu1au%$z%O*Za}YYfy_mh*)h(F;zcPay
z3rE&k;|?I(HbBemrpRoMw8q2N2p=Jol@7ZdHjnup6&-E&9hU@=uX?2ovc}j9v0==P
z{eG*WX=efNa>oR6EogN04Mi>e6^9R&(3nNEg8-B5f8u79w(fz?EqSl7UcO}TwVUlu
z!MWGxuvoRy_88~%$PdP^L0TDpXHUE6KltkSG+3+g;)M&bdL8e?uG>qp@^)p$?7Dw3
zE!o=L+GoA|H$(jyaY_@Bo!ym3$K=l7+7Ol?%LwAso(F>4-a+)c^J8tKUl+#=emE+<
z#>ynV&CRLlQ2Nfw8tLTbJ;dv<*yZIH=uL4xJ}~_);r+9TpWq6=(~%pS)t{5UxPGi~
zBh~ROsUDXW<UvEK!C}8s{_M{&yG9)O1CNIk;Nq*8B?6&7chAwf<<|cFbq{^67Pp>)
z5;$0<Tdx7dqg>CxLYjlSa<%X`5&=Q^ZPD`$7vW<XVhNO_oagGYvPt`#`(Zaqy)IN`
z!mlFLZp4dH7baE1i0Thui`r@sr#>!<Ps2@;wgwD-vc3k=iGN9AQ)&!=4AoZz7)Q4-
zEf;TTpl%e&W7GaE%Rw(^?Xx_bA?{neDR*#uXsGi1`|+GJOEp(zL1l0Y*S!zN6C)B;
zIp%zj8>INd+XcUkHh~Q<!D1#k4cNcC>hd3f$06392@elTK7|5G(4Q0-n>-6?GBbjg
zYhkPFKHDY#Sb$ZapkC|^%Il$ln0h<OvA&Q7sB(bm89#^BEx3q2*NM?#!zxA~wr0&Z
z)yhD6_&>5M1O%Hq?8&~^6{?<y_L{90?+X5EAS`3T%<~g;yi31~h`%Q9=zVgs{gvJS
zX1CQ?4{=`~h&4AZF<n*_cecy#Nja>et<!U&j1Y^XC6l`)1xi=c-8x)3_w7&;3AuYj
zoR#)R3uP#;0x@o8?eV^mXmJ+6&;f33b-5HEv<wsX|6O5hO=d$Lwl2Xc+xf_yT$+gB
zmUpVVEc-wC$P&;#Gjal1T;ZJno~2RD>>U8lrY7$5P72R{c)#@%{DTnw@}Vw5e5B9g
zFW@9O?0)X-&Ka5Q(mxt4a~f<;Oa+B^W^dPo4tU7@Blf}2iX?Vt00^FEp&0bAr|0n)
zNUG>fvWKwO>`47)=7CIT8V>Zsc}WAey&z@1eKVo%ql%Y2HdCHjGC=B&f<-PO9|4(a
zSh~l81Fw$&TZ4PX^dGPDUUN>ahkwQ{TnK$_z+~qu$4(sNT5paV8bBM3o9l1QuJ{V|
z^h(ThKa-N&=3WEn1KKX_R07EAIgYyKnh&|%5dp1rzd)?YH1}ie$o{`U-ybz{xrc<k
zh&!1hY|f9)BU-dlZzhsQ?F%I)+(I1uXY;$tt|>4?X*RdO*7q&m<Z|Rp4@f}2{UZXC
z7Uzq-epf;@ngiPZwGm}_TqVfoYQTS&#ghcqM9qS}HH`vFnfz^aK|Kq&Caw#x8~x$x
zv4^sEN-bJd>wS9T^xx-vrKml4hIs~ew~LEIkgub69Hqa&`|U6eAnzn0yzlN91_zO@
z%>EwNwmzH_^7E&Ma?;znmrf-o(3h<tAD7}dd)151j-F;85o-9`?={QSa1uC(Zc(*v
zQ6AZFKi_sF?EVEQaAUtL5H=$RVu&f}twn47y-k&YM>bpqt{26L;s_(PpN4o>DP>m3
z+3MSwZ3VAMotKbktMiCwu-X3%shkJ~_}f*lwD#AAs^Z9Epl<aebZt?+LvrEGw`5%t
z%&2-EEh=oiIiXmIv!0`k1<AuuGsNej(RX-NS+ZkD#FGVY$jciw7DM`^$C4VV!V?%W
z$-sI`AX9a>muDjA=_(;S7iDbb=uUbufxW7Mq76#&{r4GZ4#)RMnXdXeto>C3jcKu%
zU3bDtpkaR+o{-xTH@WO*Znf5?mrKm6r1-2@FR$S$K9_7hYm=UUEq;>%ud{*#w{)Au
zneS6LXKHofy|?Mc+Y`-ioLWk8=SPM9Xd0YcB6c$zFD<t|EyNuLm=fWnT3Vc`9>D&X
z=#j2!N&TVlyrAYS^Nr@+q7;{HWelj8$7zCUlfkTet2GMYl{n6|LOgPeRM?f-kGB2{
zLYI{1?ML9a5+#!tz{7wt?lAwiqYt_y%(j)O8Kt8}@P8gXgPnqn(YVIke=!18<m*6G
z))K9|!J@pc@}iqk5kG#gJYbytUx{`(uss!<6|_pr&DL|}q=q{knEc~!PM{lH$u*kB
zu>m&D$EQr{6%*1gZq<AP3=c#%IHIid;)4KLc6KWzVI|@{Y9cySfX%YI<NTuIlJXr=
zDI<g<2cK_0t-%JX62aZqR#+tbcv`Bcp=PU|c>!u5St<|rc+Nb?nxs|d3b{cr<|jh8
z+84!UXuJY&w3BTY_kn7lYe}(m*9Dd<*EYwK&~huiHHin`s*iCoaWU5Q{rho~F+Hb=
zj3@1?MRHR1!Sz6$aqYA(`b95m^zln|Lap4nG07Qf^!T$|6Ds3ox$%k*x`j%YDSa~w
z?9~$<hZY2b<ZXd*&G$4O?GkWT^Gqg8(d>cQ8K2FJXWjOJ8S5~_@#&F;CvzmtJq{V}
zBG4jY=W-+<A27!h{z%EQppI75zn7}}kt?Ukqy<c8zKOIfm&uzkPan#1sk^d8Yj5r>
zJ7Pa@7Ww4(O0)y|0mKLa8k=Nq3M~I7f5Uvad$J^MTcyr>@{V>Emj_tAMv-J1w()Q9
zZ$;^sZ}P8wuN+l^Pwlh&>e3~0C(1};a{18hyM8W|u{#Oe-x~(H_g^*zCsYR*I%?c=
z8e*zvTDNmjUsOlwycBLLD7pR1Te0M{DivCg^$1`CPHgTH?mbyASx-{Z*q?LRi$=Hu
zkt>wpIPZV0?K-6aBWuP<K#E+Td~a_#$GQj;v5|tkSX)vIoQ*sBYgVJPe4`o$Kn^|7
z`IuwA<$%K*Iz2B3B^Cxjqv`LsH}4*t!p9u;S++%Ld<N^1^S_27ZVN2!6rE%4N|q%+
zlqq5}h1TTCXYgJv_7*|R-NP1#{q#O9P&xf_fYt)$Tc}gr;SVFcmf4OFF-oIdnarxT
zX16KX*7K1R38+?7coD<kZ>z${e;;i}b&lSR>-At>ZgYOb6zu0vNA+9M`&)g-eYL~g
zXR}5{;8%CjDEMvR{-9Hk1U3aD_;nmuz79XSFXc=Az1%J=SVccaf3|6)PI_WV)ib+#
zaE~MapXX+*sHGdV8W!bhWzPI(MiZ?`;AmO;KC_1jXy|9hoF20r8~o1l1D=tDp)YZX
z9-lYW*-*>{R?;>RPPd6@9;bQ=*EUxWH#Jhz_c63{ANI?sO{%+AC_<ECznL5mCLgUH
ziyBldHOm#K43N_VMm5}VZTGGQR)7Y;S01p>5n0KAU|sePPM_7q9}|CSDG*-_TcEYG
zQ+_5B#}Ag-yv}+ZZ-WK?BU~cV_kQj^hOwf=ecE?^3Q;TwP1`cnYgKPm1JyRgh8Bjh
zPKtMPe~mY!9P^G+HwNHn%bvlqR*xCiJ{Igt?5Ww9&{v?LjNfw!7?1h|Z{TtLS|oq)
zLJ<2w44_?*amt6xw&;<d<ol>KQfigUVegQAq>lxi1ZkXt{K622r@MTrTI{RBHDQag
zIGPo-#S;4ah%k(&;-<{Jf5YdftHb|gF;~EJ&F<^p2ezAH*|5A=o|fd=1?w!XC`5F1
zbZ&HVG&I^f+KEbJUNw!%jU0(&M9N1wMI}ceqT1T0s!^XvW8pUW-1$QJQu#{xz<i_6
zxUpKcC6(@>A^GtmYQ*WPo#PsN@z&?YG6hx{`L?1}35cxNU84Ydru`vPjxDUZNv&YX
zh%V~acce~1a*vc~fbg;jx2R4udM;NSPDGe?zkq>M?GKl2bUAR>W^Sv!89p2o8Trsw
zB7LWzEAG&7<nzHbZAiJx<UAXO=pd_lqjqG)>=Nm<wpmxWmQT)j1eqCwllF3eXG^`2
z+-5ehWCvvJZNZ<s0xjZ7)Bm3UTRo)0_ch%=YxmL6u>F+q&cdLc%MO_#zJohj1nxD{
zx;M($I`ypiI6Y$T(H3V+-6c)ijq;~zU6q{OaW(jAnw+bq8{jPe>nN_654oQ>weR43
z?|%RQ0RR8xRcTODR}`*nDk_LApiHDlOM%pAsoJU)AFZNLi?|f)3Mdw}KoJCW+}Om8
zqE&2BPyrF@hKj~zKsX?pkgx>^fv`$~1SkX$6Gfy`CHI7x$2d&=!8p_V=gxiie&>Am
zJKwo4X6<hwfebDz?!0%=$RHJ_^^IGUjO(^qXjwp|<P^K4rI%9TrZaJpnIdY{!Ne;8
zJGgLKl$@Gjln*+3o11m-2!Lgg=V-OJ1UQ)+rP)S8P~pp<^2)J;ndy13E%d8ddrG(v
zU~fHXrOP9TA$g-BY;GKtwh~+`Finli%*|T}4;g&WHyQ}O4uaQrcll(C)Of-dzD#lQ
z*v+~}ne?299KJ?CJyvQRZj*~B_Q}bX3D3lo3daGnkv-eT85F~<;#t%qKOy)86gvAc
z%Ard-mDLg@f<`SL{}&g^)qFQuuUh<Fa}L;zQ?`9-{y?4Avv0^D=JcD<2474L%4Aiv
zK<|r15ic*u2wqC~)&k)J1;>DW(U}ji9|vt&8pp=C)P$MwNuO7gP{HSPMxU4~rtF+G
zbsXcWDOG$P%>Fi(F}JxC4*IUI*}0>f@UMmsnKSR3nODLXhM>gITdd}rKeDo+!x2aw
zG6=pMGA=oUr)syuNJEC}Sc5iL{r#U&OI|jCt%LLNdGi~{y=viQXS9*Ru^f_bZ}Gma
zou}rnx>q@5-FVe~$Gr*&U!sE;y3B0t^$fCs)PdAbK=>7tbH(a?)Ir5NQJt^IFloQL
zOA@Gbcu1*u3aPa7C&ymBTuJ;`OW}Fc6>(~P2_ALEIH)(qL)|ehU4N|mU~}YN9w&^I
z3uUhJt84;!;84HoqiYL!>ic0mP$#Sx^2NHLo>*79o~S?i0PBE$z<Qx?ux{ud)D8Or
z`x@sdjj{vDt+|lyk?O=W;SKN)`W}6R{zpGyAE2+WAJAXe7w9wOi+x7#KkQ5N5xsBe
zef?kGlfmARx4mQd&Y`dA?;!o16xBB|YigQ_o-*Rk2Cy{taG%JnC%9@@K4;2h`}x8_
z-{moN57Ij2h2-A(;F-8>-w|mp@k4)~sn$H=n>-kn^Hf$9T0nFyg0t~XdU4PCdGQ~O
zU&-}rAoYAqzE_d85GDksWH`IrMt&5zC*?GSlQsKj|NlZ`zrMxlr=^sit!{bZ)&hzf
z6k=v~AdOnL>~@V--Bs}F$YdHb^S~h4Dw;F35UN?cVd*gfn8Y_`{~lbV{%*bY*aL-C
zzWV%vbzbeuUCJOV*<!_lPinx!?PjO@C<)k~S$wnXNdp{aM^a10GPpT*qQ*jv7Ld$2
zlNC_eLiBH^cFJwn`ZILx{%`AFd`zmmvgM=?lyGEX<h<T%a<H<jbaCibfbY1lEPjF<
zY}{Of_dZiV)IYAe`Z02d@Rs#l>Qq3}p*_*dqUFFjV)8>y$H09$o*ti`-!bDa|H@K1
zR8&}8zN@2zpL0_cqwD2R6rT_jI!XyY3D^($O>zi~8n3lNPYJC{@0o`GEr;2GPj1c8
zAGlA))8oHy`G;En)M?>`Uf$hcS9@i}XRj2nFQA3L>GgNfES-(lpDW<pZvAM+^ZqsN
z<N2&fKVA}RDj_t&aK@?E_hE9nW0#Cj!iN1Fo24WA*SL@8G4H$L-<7ZGov5L3#(o&8
z^-$-1TmHdj#5F%{+?t{kKySFI%`UkF4E>YNX(trIa6x;|*ViS`87~g>wa<s=Z+?p*
z5;z~*8{%KZ2kCg8X@aK&3XcEbrsKhfQ-|&r!n%R`bUZyiJ^#_UjhY%~OQ7o;C(rvV
zDV+NCjloMnb2vLddyN$8W;8AG4HANJ*S2}fw@bl#`*b};x)3&mg-1;f8n{o#)8oHy
z`G;Ene*wzOeMS!e0028T004NLEZGO#jpZN5@#k9ieIsPcp7AGpCo3X6*|L)nDY8c-
zD_i(yM`Vv|vUkRl$ci#E%SahdWc2s{U7zE;U)Q<M^W69E_r1=2uCrbzMO<9mKON)Z
z0@+4!v8}(j=YPkiii<O{Rie6Vt&D_@{%?yANs$~cA`Q~xC1gQ%<bsZS84>=26nF(`
zQ3RP#963-5dE?3=Kgy#ZDxxr|pcty5Bx;~$Ts-Qc0UDzzTHqbDL3?z>2dEYI5$d5k
z-bOF9L|?Q;f4q-@81#HQ6tVro5vv=8SiiA|<%~xxKX=?0n1bn;iP@No1z3z_Sb^2Z
z6SoFau^zLq5%aJGORya)u?wZ+c4Hm(Vgn9fGY(=K4r3>d;={P(IEi0y78h^{SMfV;
z;x_J}bKG5=LL$!LAui)FuH#?)36{N@_<U;++aD+|R+mJ1v3@ToFP4*1d9nOdDlIY~
z6Edr8$ca42tMa2D3aiwrC^Di1vZzwXg|c{Ay^2>*K^0MzP#mwLq^c%WZy;9nrYfgu
zqP(h&imDE(pgyXqhEmxGHPA#Vo2q6~S4%a=TdEc6;a#a}qZ+8TQrQk~qk~j-l*;$f
zLbOz!(HdQ)vW--A!+WBqj1C;#Nqwk#<0JG#cMMRUs2-}9`V@UJ1pUQlG6tyOYM>g4
z!5D*K_#7kDM2yBHj8#)GUQAFk@C9bc-mlbTH5=2!bTtn%u~3d!ggIEM=BjV704wE)
zRak<xYMELtR*KbPty-@(VvG7#ZB|>=4z*3~QoCh%d)0m<r~^1864Vj7K86HwLLJ5_
zaZLS!vvTV^&fpTRs-x<h9Q7Mch--4hO?6wHR)69;Zpr^nsJoJLPyMLwtG|(`9>UM#
zfqI19Vw!r2J=me1D0$Hr|5M&}^-LM<#V!HXDyyw_dc9;P(YI8n&#D)+cSm_i^%<2+
zC)F<^nNEWgNU1NX)H=P+h)hVS3#%;JtDv%KFJ5KW-fc<9rgPyB<=s^|b#8H2W!CAC
z52<lQ70_NGq?WT*M0<r%T>q*{=nLuwN^7sIzM;x#udKMHUejKAaa2{(UPW<8y{^3~
zVxOv}y@9e?_o^D&t1fn`n%bKnt8kmDt-UR>CN_v7x+IE=cwH9_&=^h80`H&=+M^@h
z*B|N+bQj%OcSR5MMn4R|rx=3IL}xuh>PG32dJI0tL`)KI>8WUl8EA%CXpK2|5A)DT
zFOV5^)gS4FlCVe;7Gpe?VX~N_S7JKWU?$dMHa21|wqU;ACJEhjH@#gFc1XfbEX5uy
z7c2CBtj70PheOzaBiM}N*s6b$gr0H^PfEfmN%$GNaRz(EK79cRxP%{Y6+hy49K%hV
z(0@omZ{17ZmV`eg;V+!VJ)9Hg^#feQBV55#{Dx<^0iKZNmh|td`)F%!OLu3yyV{G4
zm#7o<-}(jP{iD5<#(Sc@)W-Xt_R<+I_U4%peWi0hoym9svLXp`m}DjwQXnr<AwSZh
zATppZKGq-W{<^5~vY-UAqqNCs$|4WSBOfZF0IHx6s==@OPjq$T6-P~!LcA$s>flw>
zM+G!OWi&xmG{-<aP@YRIjrS&6qc+-_w@iD~Lq{~k2Y4G@&<x!$NDtD3bq~`+`u0Ms
zZ(p<(?aTmlz^CYhA?S?H&=n)m-HeejOb^v#%{b}b(|nFTn23)>e=`{aF%5(9C5GWE
zjKJ3zZHCKh%3R~S&%*>P#1~?cS&FGxju}{mSy+oX_!h&px5>;ETd)Azu~;lIyYUV7
zV!v5o4w#iDK{~HC-<vh&AlBnBj^MaCW;U3O<^*C1rw~gxjkDNge#LjVjD6;+xh8qP
zNnV2aU9zs@5N_d^`2#<hKk+l};*7Xp5^)I+%|Fujig{!n)BB0^ekz^+#ZAP%#2qz8
zhPB6yv)(xqS}!S*+ZSv~oHyrX_Wss=liGR@%mb57uG3qemcc$U8LjukWU^09X6rpO
zUN-A#o84NQ1G(`s{)2*&{g=sU@0i^7uE}ffnS3&rb0)B_AQ_6-RJIt>qNGh{OWE|c
zw9R14;3ZqmX0op$D=Hv6DkBf7q5!H(ehm~xEnC#aqXg<oQawqkkFsbaNsaLunxYa~
zNJ=Xi?;w`c2G!9XHPO-5ww+K%be7RY^6J@-&`@-j(a82d6ZA%N(Mm>3+Yj$z0NP<7
zI*3B5gB>hcne}JZ8!SfJG58!4ky&TfBkg26#ZHpk8JLARm?FB_X?DI{gl=M~U0}bl
z%k2Ah6}}Ne<QxvQ>oG)(vKufSn=w>Om8@xYJHEm$OcnF&HtD<<^Taayon3DC+Z8qe
z%f(vzgI#A2W3AX~k76TEU<(dOx1X^SXRuZ5v8SZ-1?&;~?Rk5^UX(R(83)86dkqJ1
z-5#+waU8dC5?Ae6dqc)O>6(bsc!+cMzN9?HdHYCGp2(X17nkr<Iw<G;ZawR6+Q50S
zlqAl(Do5S0FE}rSyJM5N`}RfWJ+NLH=lx?-x<@v>^Pbpr?y1e-p4p7<zRBdgXZD$O
zE(;=LM;vm2w||${C3DGLN|zt0Q4obt#ND$myR<H!Bo%e(Trroy6?Yd*Dd)X}vdHUl
zi(KwiR6u1^MRk$Wy@`sbjcUl@>Yx_tqoK&|8lxVvxu$61vbyHaB}YEj3I*^kin^9E
z-!?MewyvEkfet8*CG!4I&UM0v=z`M5>*hTFnJ7{BbY36y6V2qjG?VV1pfv`eh-l*m
zNa7I59O~M;VUjc4d9mzKDDOt2GsdASCZGqtKwmU*{oNE9)6oz!Q3G=@4-2tG40h9G
zESDLr#HUz;NoK9{M&erx!zK(672Fn##deIxF8mkYVY1AxzMJYAxMecK6>dLP<9qZs
z2c5Sbhp`byu|iaG$FU73u@k>w56+5x?mQB35&qn?Ggq8<2*2S7ZeX9N;%?z5{E44&
z7iaL6JLmq!ulNTu#AWvw*YGc{BQEgHNo7?R`}%v^c|6sFJI;#&@4icP$pY`8^HK!f
zm+nhf#-s|or_M_oc+Z@dA@CBV@(ni|F9qHL=Vb}JrEaN98stE7<c41>F9vy$2KkX5
z1(8%_4+<j}ieZ+U<&v2aftL@ZQ2^zTU6>#fVwbE)5)=>SxZ*(tR7O=~3kn4_P!zTB
zw~i0IQmBivXn;bZd@$El4;rH;nxb~l9CgqN_3<vMi;BT~*C=R*Cg^}>K}WR22WX8h
zXe25J3tiiw8{R`ttkAszuM_&BGy0>gcs*F`x&{N$1A{R|4-LFN7><uI3SC7tc~swc
ze*8GZ`iu`o<G&bz$ry%s9j~VaU*ao#jrmxFrC5$tSc`SR25bzrU~8}gyRkQz9_+_#
ze2)b<gk?B_)i{pz!B5x}oWi!?7wo~=U|(<^3Al(Ka0Nf&HypzaoCt2=WN;fl2Y=xV
z?gi(9MEr_}xQxfRhJSG#v7bV220HZq2>1vc#=3@F?n&SKL9)<$5O^s<?@{2T3caU+
zmp1gC1zv{GQ{hW6VP*u#h9n|sm<!2~7cU||(x4#HqcAdtMUg2ij;tt!93p2}M|*if
zuN?BBJPM#93ZV*$gjG>29Id^Qp*KN$Wx_IHlHkqItA%(}7xlu1XcRU<bF@OdZXJ4U
z&>pe1+$r=Li5g)S)J8Yd74^ejXdL!I3w(@s!v1KBf#`_AXe>I1!_XBY&_nbL$DnsO
z4*f6zNyCZw1d}ibQ_))t4QF6DW?_^V9e$0m;XI7TLfj7)hu#-hhACKqv0{3-8Z)sD
zv&EcnBj$#iu>jk!DBOXi*p21bi@9P=cmV5hFkBTL!iMliYzmKJD^6g0coMtt3%<iy
zY!c_g?{N|P!%NTmUU}a4+Vfu55zDw4p2KaN#vPo(v|w6rKYV~kc#3Dxk%}xrBt>$h
zLR#o3eUuTIkqwDquE={FdU+!+_CscVWI{pY5IQO>NktG#Dvp#<DWpUGC}&g-dGH$Y
zp%Mz<brcf$qc@~;breHQl#Jq029=_+Q9YDLLsZ1usDfswCMrcOrE@F1iFZ*uYKNB5
zTTvaUs~dGd19U=Tbe6qc<hWMRM`#^&LyxEz+D5%mJL-dyQD3BxbJ8~IhyGD}>F{<m
z0I}J9ioww^jKF9Njy}gkOcGO~anZEsOU#UBV@@;=^P@#r94(MB9!v3!SP`wlT6~L5
z7!pmDjOo!f%!+p4>u47iM0+qU+AZhiJFLTgY(^|?IDWubu{HV;{yuDvj$vnX0!yNw
zuv3(cPGKKTBLT5{eh?SnAIrnhB}@=UqpMgNUBmwY00960q}XMYRM)lz;O)B_cXxM}
zKp+s@-95OwLvRTa+}+*XA-KD{dvMn~&%A%``~8^XJ6+VNHP^0c26S>ZbQXEcBn>@b
zE;G-WYs@(22J@Q9h$IfZVLmV~nem~!%opT#=pL^B?f%X@V3vnoFfuI2vrrHd7UZ3G
zVgLOVfr-pSV`4I~!s0Rsn8ZvnCMEO#olg~(mdU_mX0kCknOtFcn0!nDrVvw<DHc|e
zDZ`X!Dlt`=YGE~*I!t}0A=89O9o8IaiL^#?Bkhn5NGBw9Sn;rKNDrhJQXT1o^g{+9
z#lr@L4Pk~cBbm|6IA(m<L}m&zoteeVWd?=KWBzx30rP+7moWc3zm!=Pwvt)HtY<bc
zTbXTPJDEMqKIR~EggF{^f;r8cV=gjRn5$tonA^-f<{|ThQ9=tV9EC)f2rnW~BoSFe
z711b$h$&)=I22FB7YRioN+Obq<RXPgDGa7lQyLKt)9EO^h=S>ilu5+GbQa1g5@0$z
z<q*lxkW1tisW6>~@`|*W&PVx0222;Af+91f3sGT_4bw%bn8=BS5~7qSEy{{=R6$e}
zl|>b*CaQ~?q88N=bwz#AfEtQMqKRlq%|#2*O0*VjL>WxCqxPZ_raMw6Q3KOmsH>=l
z>F(4+G{$r<>MdHLp%2ne^vCZ4y~IE<2(Jos7emBQOb5D%VPZI110BUk@i(>w+KJI(
z3{H_kj1}X>1ezo!i>YE7%@8xiY%z!CiTProSVT+2Qn6gDpjBeEkV4@z8Hee0v|dcX
z^hVkwW?^~@Z50bJy`6T5WtiSYyTuwbq{OGX7SnrapV)xu19VVq#`Ix2BDP`r7#$Zo
zF@2Iwi9Kk@jnCwaI4922MR7@75m)KDxFK$d+jLjl6A#2gdMuuZXW}`%5HCeZd}sTJ
zvzUHEZ^dOyzo!r4CZ<2pXK^3XU+J58is>KpQ@lb$8$?Pae~Q2lVWgGcMBuCNG9*8X
zz$Xz-hL;~i;Jt__BgwZS@J2+DeQ}CtGP;Z@V^JI#SH_nKkc2X^OhU<Ia+y-5LQ>1L
zG96`*8D(Y}0~1*(n~aB<97rNe<f7a%1!nRfX)y6C<&&A@Z!!xe3exW~J7x+aIb{*f
z6r<vpDk1Y?q7;%(=9gt;Sy=!R<*9-!gqcc65m{Ll#Y9!AhN<eZBxY*JnzA&eYRfvZ
z946{feOXa9z|<d@YACB<qA@kWR8!0}!&EKVf?CPCm}o<7G1ZP~gjaW<rm~~#go#c_
zXW1IB>Pp@4s`j!6_2g@M$v<Tm-04gGWOv*dK)vKZIY{=FLu4P^`HP0h{<t%OM&j;Z
zIZBR}Lvd#;jg!N1X96-({w*ib6gdWWrqOgcUe3VXnQ|iT%%-_=D(=ju1-QEqcW24P
zv_#I4%j7)VSwSmtca>a(J8NjITq@Vg<**V~!&<pfZjxK%R#+xC%N=s3+%0#>y>g%2
zk4Fda^&fdi9+5}oF?o<q$dmFk9vziu>6AP#FW_t7lDsId$g6N2k1xxcbW7frcjY~J
zAkX051x!4mCwSd6`CML;5AnN~^h!RJujOO;MxK*z(f$tm4&wwL<a=B{(Pa4%*Dv%{
ze#Z4X{g8X)Px%dxqzWjtORmFTH&Pk0YNy;TH{e}swGmgRyb7ss>R%aNy}<5>DiBFM
z!-=oS$haG(qN->rf{LhOP)xk~E?%8U##V80FAzn=SNUZeJPK1~WCC16Dv?TvYZ6MT
z662bjQmChRXUWtXnL<TYDOD<T2BNDp>OMZ9J2I_GuQE_36&0UdIy}lkSyeQAnwixN
znO)^jIVrdL4WDZ8H7}i(d2r27$7McT3(^5u0N28_Ulzi(DD9R-a4kVQWN}<e(RNu9
z*DbP?+8|5gYj*C>MP*evRY6sx%BqU0rm9m-RZG=Tb*a8;pc<-1)I>E^%~cDerD~;G
zt2U~wYOgv_C)HVXRo$qE>Zy9GKdG<krv|8j)E09?)KK-88m305ku*w;R%6vTnyAX-
zB$H{1s*00Lqv@(PPBN2bsXuU%IW$)_!&&Adt<(aXrajshqJ5DXj^-t_RE<INa$2b-
zp?Ni}QB%>pj@GLgXx>Ph)NC|wp{;5jy0;?>)DG27bwvA4wC_@j(Y%NLQOnT0j}E9+
zXg)-T)mk(krDJLXnorP4wHeK)>5STj?sLcvbsnE#XS82H`$e@I&6nwl+Kc9EbVD6P
z^DVlqj-dH2-BZWW{D2;+Q)qrnPt;j-KSR!|=lBkFMf<;Kf1xg+`4zoZSJC{I-m9Bv
z{z#wH9W;NTuj)RUzta!(2+cwVq<)HS&0k$qcXS(V^h>ll9q{@MTElc8y#9dJh&m8S
ze@1H*9f+#Gp*6Y=#Lz#{8cPRa>!3I;^+d1L@hJgXLpl*9*5Pzg9YIIb$ti`7tfS~u
zlv+pA(REr%r(>cq17*~4(3qKet9WS4O4)QmH0Ge3I*HD$lc6yW<<%+Cn2+-7G-xb9
z1$BBf7NWvB6B>(BAC(1-#i@kOj>b|{TIbSb^)F~FPZjj9Xskq)^>1jbO4ao5Xskgs
zbrCeyrhcjz8tYO$T@sBA=nq{+H`3+M*o2zuifC+3Ep!z$wxZU$IvU$jJ6#Kn9cX~6
zqdOr%YZvOu*6!3(H%4b~`cpSUXJ6{4TcUFS4b*MWIhcm%_UQZz8HV&%gH%U!kI*A^
zXEcwZ(YhO&$I^J+3(XU0lJ0}%DKu60NAq-=p$DOP7R}Z}(LEd)s0OPMdM@%eI_J{@
zJqDeNXo;SH&SkV*Pe$iTTBWC<a}BN4Gts%8Ht0EM-9($wx<zl*^U%4ScIbuZ+(moz
zQgrU6eR>5t51{pcUab$J^N>D_-XnS)I*-wDy%C)!>9jtB_EWe9{paZd`!CUD^zYVJ
z>AJqD|G|#ibVu*Uj(f;K?6{9D5A;Lqd8ChG#}j(0PjJUseGa|<A{X>U{gPgx`-*-;
zZ}lxSzNZiRE;>KzPw0H0KkF}O{fgFa`Y{@R&`<pgt<p$?e}2tbZOn6YztGOy(82#G
z2!uBgOe7NykE2i&vs*_pQO!H-h;9NObqw=K$276bXC2#o!IrotuKA|pneW(<zyuN+
z{Ma1)MPiIeY^+IYl2Hni(xf(ND4lV5?uAZnUh0(Ql}>G5>$K(#q&FTX$w-+@M3cot
zGto^p%5Gws9N3-H#5TE19BlptiD&Yd_$Gn*mGYTH*!&w6FiEiYcPeC3VQ&#CiY>)V
z25cxnB~4~i+GNFsvQ!RR%416fY{_jZQDtnYV)9}`HL7m%V?#}<WeS?wrVbM9sZR~K
zry(^m#j&LcH8rJ7b5jOeT2gt_3VT|cHl`xBw4*Afz3E`8VoN9LY-(UjSE_BgnYyMP
zHuRvL<_~P>O(|t#Z0JMHOkdOA#8LyyKr`43p})*9Gs29dQD(FmYsS%d+?iyCn#pF0
znP#TbOf$>OF>`4iPP4!)G>gp=T4v^(6=tbfX;zsvW-YD5#`R_pwrr#cW;nKNrY&X+
zwrr#A*tWw=#fDw9+swv>e`qgGu+J>Qh68jETMuDd7PA@~j?htje#gxQY&c1$uxXPy
zLuavRt2s{>aGHx|CpKKBD`pQiT%+qa!40z?8*b5UY`ueRQO#j&xJUQR3A8_?N9LG$
zVoqbjGkR{$VZ#f0iPOBoDPH6H)?C7-_w>PB#irNh6Hf6NJFer=S98;R!{)d8E;jt2
zpV%(!18h+?pzUL0ZLq`JK*$C=!r4H0`xN&gT5I3nJdy1??1^dv(d<X;dShbPw<f0j
zf?Y9eZ2Q^7v0qJG`_0_c-?1yc4MeiedYh0E*)SW<CZVKOT4j?{3TtdC8^NZwQEXJ3
zmeSei*qebe;ciBpCyR}2v)ULon~i0&+n6?ojc4QAT$J0!!5OlnE00Zv_B{4in~#38
z1?=y(5EZdSZE<YRi|0#HDO(WFm!YyYu`Pk;%2Nef8qZatM7A8BsYccDPO9Q-O-NvC
z;A;ZB%i1<GPFcsM!P)CkeVg6}Xa585t0B(b2xre`8`~!K7s_p$+2%GU{c0EMak`}q
zw6+cLFLYZQXlIk)r<M*j(9vF3ooyEz=w>IV9=4Y4iM_pSU2OZ4`q&28*pK?-Gy|}G
zAhr*}_VKzoHV&boI6-4Po8PvuE$wg`VcXcY_HP<x+uIIy42`v&uyH(1z`LCUlX0zL
zr_wZhZpG{jnrUa*Id(41w+rkdyO@^Yq|0$#X?x?8t7(nxi&HMQJ@sImXgzJPf8kV{
zXfsaV-)^M>b{I~$op#ufIN>hZZAasT|Il7L4kz3XmGC{PVkg=aIOhSo3cC)`VY?Rl
zR@g>*9d;d~<8~u<oupHC3wE8M-|aT+I!EX2PVBl!m+T(wx<Xg&KJ2@052E`f-Lgl}
ze~0ecqv*d+59|r_KcdI>H2R-XVS5(+&*@)#0sSxOmA#DqH}n>N0^iviX#YST?QOJw
zrri1-+P~5_`w;Cv=%;;xcIg5|?K8A%7clMx+MNq{_ZsbC?u`xS-lI2y3q*9E&>PtW
zqPQ>UjphQ;-FNiHbb(k-IO*b0af=__-5eX3XXCl}ZjlWvu?bxwm)M1O5nNJA<|4Vs
zE(L|xQC&2bic-56E~ZOM>0E4dW}u8N38per7MC1T+33AWiK!fv)1|>wZu-Tg$5dX5
zi+5^qmOui$(-6)Qa=*HnHXp9PxdQaNE98o}qEO6b!j=+L(q+Y#e6F<1jtyn0oXdp`
z6{w=igAJAGH#gf>am6uJjjFp+n9T2Lx-yumO?6y(Ox2_Mt`esHpn`6$ZRo0EsxdWj
zH85GgtyHxz)tp+mx|nK3tz835wWUICzHR3kVX6albWQnwCaPwb>Ox&zOH6gA9<B|h
zdQlO#(DruiG1Z6qx=xrZ?E1UTm>NifTsKS&p`orPriM`>x7ZeT!`%o>1^T$Z-6+=^
z8STcnzi6l%=f=CiG{{YKliUEB;-<RkZU)VAv)x=bj~2LvZn0ZJ%iMCe(ygL3ZmnDI
zHXs|_X14{|>bAQb$WFK0?V-JHpF7~Dq45wMc5~2pl#aPYXgonD-3m0GrZa9GI?o}S
z(0KvbhR#dKE;L@Dt8Tx$?hc~yCf#yJ(0GULy5nfPPY>KFG(Muo?kqZ=BInWh9Jz?j
z7szEazM|LentSVRpz%F@aJSL;i9Wk~X#7gw+(R_}pr7suI;EHX89MPFwEkaoS|4ye
z=nVNlIR6ff5quz`|A@xOJ`lx!L1Q!@i0;3mF{Tg1@&cW45QWZo2z~(Y2`HiWKI9Wq
zQXk$&@X0BKkK`l!RFv9B_0fD<O6Ozvm_7q#^s&*I8HtO|tVnz`W~UrJ5gKz*8I=r;
z6;$9CpAwCEsfJ2}#(b3Dr$=J}D(Ewzvk;O6okP`7l?{y}RG^5@iN>N-+`qNgY)Ky|
z?I+l>zMKzK^p$*=EAFfMYQBcANws|)UjcXO`APVb+Q0`I`bNHqZ%WO4bKlapqSn5I
zFX`J-J6{Hk9jKG9jGe807hfHX-Ke{-i`Gu~p49gZd@t(l8>6ug_4O^#*q;XYHfS6~
zgME864y9qfGg=4u5xy%L|E5vC2U>^ux2`7|$I>|eCmJWvMBfjMlWB?{h{kC&-48+I
zOq%V7qjid(>qnw-J}vO0(K_2N@?+4rgqHg8Xk1P!{3JB4qSbyX8rRY~KLd>$Xp^6V
z*42KCpNGb6wB0X6>n6X`FGk~T+T)j@aWC!jE75p>4x(+nKTJpbMt{t&!G?qWgx`$T
zWB!!iipDc^*6%>;34h-2LgPic<o`kA6}sy8qwzW&^*8-3bOsLjJ7~P?@1gIQzwaNQ
z<Ai_cAE7aD#y{~-vFVh5=AZj(bk)D`Fa2e@<X`(Y{sO)8@BK&riN5%+{=5G{A|yj9
zq(e63C=?0{g%3rb$e}2qXrbs7GZYK2iyev^ibn}T2}6lPNkU0OAJCXQ6i5;JhQ?H(
zKx(`?PADy<!>bd8GSGhj00960q}OG%6-Tf(;O?sKo*vxY-CcsayF0-l5FjK#fItWq
zG{GSuxCagH4#C~s-7PrWch~xUeQVucclCPa^q!e|tEcxq9Cjved=Brt_qjxFAI-<`
zd3`<~+sE|<L_tsfE1$-v^@T+dpWbKi#YAzR$!GQ@MJb=vXY*x538a?uIgwaFRP=d}
zTG?0e`H@&nR7Y(MUkHh{L~UObiFHLiB$V_GL_=R1iH${lq&D^Ck=R^(?JFU*rEle{
zBC!oJ+u(W$-`3YeVtdiS*Fj=O(FqCP`ex8ss=H#eKz%nP{D5Oitgt(dm9YaoaLnst
z`kpw(@IU(K{(IjV`MuGhx9@`ySzBMx4{K}U`}+ZYpdajqh@pO%AK^#(Q8@d%A1%iE
zaejiI<bUyh_&@zb%%=KjellLo^t1e5cs1A0^E2>jl%IpTQGTHxBgXkfehKz#ng7SH
z5Uc!Zzt*o48~g&l5qEF$Q;@mQZ}HQSxyf(yvyr*Q@9^`HxeZBk@m?-R^=`4puR--*
zu@C(Q`TgPmRzKVy5{La}BpwyV{8pqM_9y&yq#pBsBYCGkEmHekNIWaf`G1jkL0rU>
z@AsF)Wq%NfSH(4d1c^7qO@ADz*YJ#|kaS1f^(WEcrhgTjLE;1P(4R-*WAOyLa>+jx
z&-@i6z7Q|{EC1R*7tj0~?DshA`&<9cUl+Ihd!&BwAJHM=j{l6EU4*;7cojYtk35AO
zKJXFuy$yrVhf%_)B6=7jj1|Tfal?3Ff-oUQqA*FAR3r~mgsH;RB5jy1bg2Hn)@KYe
zg;~O^B72x4%oXMqdBc2Rfv_OPS7G6>h$t2o4@-unM47N`m;u!l!ir(#u!^V_Ru5~2
zwM3n;ZdgBTfYC5)95xZn!scO%u%&1nwh7Clx?R{FzJX5A8M?rC@I7=7dxSlq7xaO?
z&=2~<02qX47!nQ*hlvs4$ndxDcQGa$8;-{gBnT&jlfud3&u~gOHJm19gfqjzsQ&-!
z=b(0u)X)2W_4C68s9o^C)GrQ~gv-L^Vnw(zTpg|v>%#Tn#&8qH=5TAcP3#DFhP%T(
zVsE%FT!iX_;i2$IcvKt@PlTt!)8cG+F1!$4#JCh*39pLl;f?TCcw5{J?}Z0Y&59H&
ztDM5eOGH)CR16hU#8z=sJQZIgREbm)l~g2GDO4(zTBKD*r3a@ns!S@1h^DfrtSXj>
zqq3_UDuGC(a;e-ZnMk4Xs(dPqNT)KWf{;-aQiWA!kwq0%#Z)$tU6oKJRZfvhl~!d`
z9+6j-Q{`2DQ9xBum7$udu4<}UqK>Mo>Z=B#k!q}(s%GMA)k3vWtwmeaPIXY<h)$}i
z>I~IYSM{B$E$XW8)eownXsmjuo~pTMp?ay_s*Px;eo}o^NAa!tMg0r|@aiWuNDRhZ
zz0|K_D5`s?;bH_*zgMHgZ|K%ljTU3nSnSDIHC|0nlhkDKr<$Uss%c_|nyF^1Ibxog
zuNJDm#bUKYEmO;dR^zZI<B+gYtx|u8DQb;ctEP*YYQ5T^=8E}hliI8ni6v^Q+NS;y
zq*kbkuu|<*yVPp2M(t7ms&!($+Nbubjbf8Js1B(uVyil$j;ie<y*iFPIj&BsQ|gR5
zE6%G6>XN!FuBvP5hPo+kt2^qRx-TB8N9u`sDk`ZH*wYh8c%fdZbK-(}t^QM2#5MI+
zy;HZu9rZzdR1d@>^;vyU&qQUVb!Q;`1!s+pu=*qJ4s?Xq?~oEjM?}^CAw9Z|h@oGi
zYb+hnU&YmNbX*-@C(wy>Vv$rQ*U2EIPNmc6v?9IEs53xjokeHU*&&|J2?=y=oku4T
z$#e>x4^rs@x}Z)cGU!aY5M<FsbWxoHa_QobTbI<ObY789=htOWQ&#HA>k6o-pbP7Y
z5UHvR#dHZ>4NB@7x~48I%AlgQE~o22d0kJ}M@4;INjK09QBy@Xf~vZ@ZVEMYbN#ih
z1+{exRJ7D}bt|Z++vv8aXp5S5sAvz3^*6esZl*gybKO~YK}8o-w9>70S7@WV>F+Vx
z;qD*g-fwkx=&XC{A8|*-ce=OkgM0eupKwnP-4A-|{`wc)OZ3(Q^&pI&bYDGK|Ezz7
z{(6`mt_O&LdZZqOD-lEVP(2!k>9Km89w8!EC%|udlAf%`=s#ero}&NK<6(lHswe7c
zFiFqQGxeW(7EIA|^jut-i!1YBx?Z3c;_5;@TQ7n+dWl|&D@*l!y$lxUfAk9dH!Rkx
zaCVhm4l7`-UIXj(2E7iITl8kVRd3U~^iJ3fdtiq?p!dT;{V!f0(T6cR1_$&B%unf~
za2oUD`iwpc=k*2ryaea;6@67-6F2lt+;JOso`+NVuD*wR?(0YTp}wOp!d?ATKhqE3
zo_?-R>K90UiGHtOpMI_XgZ=s~*7F8C^#-%|Sm!(a0p9CR`m^2+AN3dbgkL72tvD!J
zsTFB5WNC|z*rb`(VCHc26wqTGg`ViB@Ki^q81zubpo>^>KrwMOq-YceqEkGIPf;K$
z-Ovds2PCH06pNBVA}Xd6;Z-uKs*|B6Ii&zkspyeTNvSC%W@#xMJ<#bW4ZXneg-(zB
z^ppX)aVe9yr8CoY$V%BL2gTMoC<|THS@0?s<)%0~H|3-&Ixk&@e3YLGQgK}nKX2>E
zyf77^LUc?Q!mJn-r=oO77sae3m7)@~N0-2?4E58cX*ZOia#Wtm(l%Wdvx+oYSD?*M
zkt)+hU5Pe9Wm>PR&}zs<d1Va+(6KsIqv}+XYEd1kE9z4NYDA4iQ))(EQw!0GT2otU
zCpu6^`i8!xPSk}u!*|q;exUB6C;dphsgLMO{ir|vA_meR8iKkmG?ad&;WP|J(kS|!
zMvJjDjwVnInh3RMGW|jIL<5>ae^C?BjHc0aYAIUNOqxaQ#Wz^V9IRw6lD|XpJQ_js
zVH7Q-ziEsZM~i6*O`@eRnU>Q(^rx6YD`^!?71L-9t)&@aCatFpG+WF;@+Rs^o00rI
zlDE(T+6oJ4JMExFVlnNa-L#DMz;fD4`)Gw&NeAd4trlzOFdd<FVm%$B<Frw1Lh?!M
z;VC3{NAhXfMrU9<oul)#Q|zLPbcz0@%dnTO(ly#I4$uv{Nr%K?x=nZJs5nOV=sukg
zCz1RR@6ID6_eAnzI!jOB96h7wbU|FCm-LFR&}+C#Z|E&u7dPlVeV|+7HhrScbXVNN
zf9NznpooV^?m^6tsTXpelHsRFb{y~vig<}H6H)j-iVANiI>+F56!D&7acurb5uYe7
z$Kx*)p*R62WD?AYxi=E~faRo|45<;GQ*cV8rsQax8lrPrPRB7tEY83gIgW_SnK=u`
z7YR5UXXivBG3Vr;C<*t4q@0`caB`7?^KpJo1!=e-q~$_fnA3|4T$GD(CXty-a7oT8
zvT<oH!#PAwF30^S7yk^oxdK<@ydoc0<|<qOzT#?7h-+|7E+UF@ZLY(`MG3CQ^|_QN
z%?-H`mlfr>3HPV+{0mgzX55@BiOSr9TXI#X&aI&ax8-(ROVs8L{0-L?^|%v%%MC<B
z?!sNUv1r2GcmOr!fzXV9kmMfR6RABVxfl25Hrxl=a$oMp?L`Of&%ba-(TNB0Anq)>
z@DTo$zZ2bfI1fYiaQ>bL(@11@=TXpuf9KKsqv*wBc^v=5<DoB4<VpOq=+A%fpFBVe
z<iB_-4;DjsI?v#tVi*si5j>KI&?x?mXTk40hv)JbF_!1^0-nGNVInW$#XMR3!Ap4=
zPZ59dKfHpciRrwGSMy9Ug4gmoB(CEPJX_4+O}v@siTS*hxAEV+9TxFU-o;DAQr?5q
ze|b6YmE3*2lJ~<ZKFEi7jabV^_!u9Dqex$m*?Q!k;?u}KEj`ZiIpm*{9vAo$U*;Wr
z1$Oc^zRtVF9(1|Mx6tP{y6opWe1Px5LB7uq_^>#_C;0?KvLC}4e#*~~8*!drAo(i4
zglqho|Kl6*Ml#>>JEY#>_elJ}_xK~+=g<6wABtm4Mw=5jK4Jrpu|j7qvp1oMYNA0D
z69c@7WiH|wVwr9HiesBs=o8n(!+68-O{7Z#^PUsJ2lPm6K5-KG%*jl0ql7lbu(2j3
zq&8_xI+I>xG?`2m<KPLWHcvUNdB*9@bIxd9aAuqjOjZa@c9X+Ihd3q|#5H+LUK1Y@
zn0zLo$!`*w0+84wF<+UarVu1EMNCnXLa-@rN|;n4jVWbHn+zh8DQn81LphVpl!xr5
zqN!wZLT*z9@|bF-IyzKGhZ^Wm6AGHzrVcvPF@;TCC}QfH2Bw%OZW@`!rX(~qP0**Q
z^!nPgK%W-UtCeYODu{}vt!ZZ}LsioOs+o?a6Gjd5E&9|louRhrYQDp$YwDS9P~ZGu
zx|@cgk?CoE#Asrgn`ZE}>23O$mZFvEYtrjBqOJMa^v7s#63_rM$P9!b=2tVo3^&7I
zq#0#?H>1T^GY;2>nek?*nPmPjlVOVa%S<!V#Y{5`=ewENW}cY~3(P{y=a@z2Z?nWK
zGfUwgvjX$>W~CWyR>7}k4dx@v1Q=!3ne}Fj7-u%(`gpU+{Ao7BUuLV>W@d<4W(Uqs
zHapFHvkMlQJ($lm|Dw}gvmBp+<>;~>^UCIcS#1u%8gm5m3FauQHOI{fbcxttPMOo_
zblPk&XJD&2XU?1LVu!hi`!}0QX1BQvd(2hLcbRMGblvPTv0=ZtiTPo33yzpO=B_yi
z$IN{=ZXTLP=A<}fo|vat$x|dg!+37anwN0Syf*)t3*sVH@D@GZNuMj`1G;`RH_S`E
zfi9oT7j*t&ZX0E9qZ3)S_wa}Ep|Lj7C$Qc=F%eHq6szn56BVqDZtt5H=9P(QU!hlQ
z8^^ve5zo;-o_%lP!&~%BXg`=l_M=H`Kbb513H_4V#Mapu5ZL54g$+d%o64rPT98d^
z(^&)2Yz7<4W`x)_v&~}TiWerE&4DX9Y*d@w#<aO?Je%7luz75Jo7X0>`5>t+U<=wr
zBB3pWtU?&6p@=PNi`x>Slr3$`+H#^WR!|t%D%r}mq$q8x+G@7Et!Pt1C0oPRv{gh^
zTie#P)v@yGI9uN~u(fPGJWV~!8`*|78N{$nL{6K@HigW#x&7K^71{7iEwS!aSa)mN
z1}kr4+uE<N+V-M@&5x&8z(rsjceLNyRvb~4yV$Na6-PwjZuWb7i6Xk&9`;8&fg+lq
zXCK?l{shf!Kl`(7A#&JXaK5P>fX)Nac@R1eMyK)I0Udv}L$QL0Huf7^3_98ocBJhr
zy4c_B@Af;<&5p5S?GK`d9d9StUUnj$dJ@KDsB8bUQ?PGSWbdZhX?D7pX=m9vcCMIj
z7udhCnnf6kZC|?-`q}07AKPF2Vu#p1JOry+Wmnr_Vz6Cn*J0%m1MCLW497b-9IM=9
zH``I-H@nqtv!lfryTk6Z<HUHo8}CtVtg|lGxWq2AdtkZUYxmg|VzFImTk}fvJ7^Et
zHDawjVvpMOVuL-7nvM9xY(&44_LSWsw%RlHtlchl*z@*+-6j6Dmte2GVz1i$Vu`(u
zojG7{z(IS<-nNIu5qsC(v&Y17`v7$(Y#BIdAKAzDv^Zm*+GqBhIB#Frm-eEqV=voR
zaK-*--`H#7x_yUzziA8dO?3ERKiWIuuKjGk*!$vv)$Wyjh)?`Obikip_ryj#wax|Z
zxs7;XL-*Q7f&XkY7u~(J5$|kF_m_QdW4ZUpj^pCGPd4JSjqegTCA3TAqBwFD@Rgd~
zOgq~qaS=(~e7n#lcM(xtN{HrCyEHC_i0RU~sW!HY%CTJ>m%(Lp@kD%=*=2DFMIx6C
zb*Wr-NbPdETrRCh=kmB|HobdI=~0=_<#(AxW>?UC<+6%wt}yD1OAHokNroql@Gb>v
zi?|}LC=_=kTq##tly&7?4p$y>x{9un%PsP_ysir5bJbjRS3neWHC-)NNECK;P&308
zcNJV+sOaju2ClLw>Z-a%P|Y=QO<fI9(=~TryV|0TYl)g!wxp})T0wo+#<g_~MKRad
zEv3e&?BKp}%|vt8$$jfuh?cGkYUbF|uC?n5ZCp3^y=x~*xb|)$wMS(S*VA<rom?;1
z+jSOQ+)t>PXUn?pTwmzses=xc52BRo?gl^)H^>ckKZ;)NS2xu45kI-%s99iBxDl=l
zj6~JX?l<wf8;xtD-556v_l^_e-4J9<6qDRQWdG^@fWO>SH{Hz;v)pVq*Ub|P+(Ng=
zEf!1NGWU;LAy&CHZZ)iP8{B%><hHoYu+43EJKZjr;`YE)x7Y1+GsSFoz#Vk+#X@)3
z9dS#<GIz`!cPqtex7MA6^=_j(4V&FocNVt0^X`J%4gb1Ju-9F2SKWSbz+HDY+#zw;
z-Ez0xQE|-Ob@$u}anhY~58$*r;~v3Tcg{V5^X{2@?k>V*_Y$tS*X}=eO<Z?x-8*+v
z+;Sh>M|Vfub)VfAcV9ek51kGkBJ(l6EI&pjetm;yF5<cKK^VMpQQ)<U7DNx;xQMqd
zW)LfQ?;<|9I6>UtlZ*K5;s*%=C3HYRBH+LTNx%llg5*IUe2_9o6+{hEL$n}mkS>TJ
zVg?z4j6rM>C&(OR3F3+PLAD@!kWeHF5(haUNsu(i4atJMLB1e`NEs9e3I=I{s+0z)
zHR-D$eGri$C>%7TOd@koG$<Bi71@FkLCGM8$Qk4cN<;1-PcW47Aax{_4e|#O1%h%x
z`QW9k5JdjftQ1rUCfJB-LG_?UP&=pvp{pO%3+e|AgGND<psA=COmfX(qH7Vf3R;V{
zLA#)R&_Q$zI-%x&00030|E1GefX-zY1#l(Kd5B0MMTzW0gvb^#)`$?<B9eU@Oa_B#
z><q)$XBf*^nz1iq-}f~;kv)61Oo-0+T-Q|>zl;BJ&hx(S|NDZ_5kaVq&gg<#t`54P
zJL;o{R6`^5lCRMlebCI+LkpylmiQ6<(8jeze+)o-*8yS5biyF%j3Ed|SJw>@7>b@4
zCLIxs5z-r@FdBVb4}{=B>JYU?VjMzUe?(yd2D&grD>E3AWC*5UDt>Yin2s43ju;t)
zk(eo?5Qo_q<Az}@=1L^uF(2by6c%C;Cc0=OC^H#}G6hSp6w};vEXNAOVx>&NEUcC|
zBw;P)xERdCdWpwIY{CM!5L>VnKf45MS7s@8$TIB0Znf>kV(gVAScUzPh&4DUNjQun
zSm#z^1CGi@9LEW4c3W@?r?Jg#$6000VUL_w?@nBhi<)y0mvCAAmvIGG)qfS&)OZ-z
z-3{DS?j~;GjC$_4yEv)-d+t7tY3?IDl*dTM6FhZ)<2h3BC8l6XzAUMjnxwh3e3j{$
zfv?Hy%)~dDS+eN8EPRXE_zttXcljrBN-}aWH=iNkIr1>C-VDgf{4Ah119Gqs3-cA`
z)0`qct0;?UPBEXAktHM(OR|)5rT9KeOD>jSS$^R1upG-P9gv?DS&0Q*Ay#G;HCEt<
z{9KB%8iV+eE6$HuQ%bTn>nK-;pYRJQ!+QLZpSp5vz^{}JsK`ca%+Fk9HepjWHe^+P
zBh}cPEm+;vU@d+twb@3U74*C=zms}w$M0F+HQ*2IpmYa4?Z{?Ick+go>>{n$Lk&IH
zjXkx(p1SU(6X>PpV4ZrfR-TT%*_z#@4g0F2q1G1S>+8o*tu55oH;7>}gyD?fP&Zt6
zM#^Z8VI;@7ft<+6oFo%Cozo<mGZ`yWIft_)hH*-b^lUuGDzm^P=)4lRm`k{f%iT(@
z;u_AAkxb$gu9b<L!SynUv$#>FaV|GYth(cwsBf95hHc!=6>b%Ga+hWW4CZd`(aHmc
zai6~BI_{Sw9@Gvt@Q|$KCLWRXYMn1zcudz@nNzm$1W$5@+sV^Bqn!ur=2`yAy>1`R
z>$i7+zsW&f)R`RSCE3rT+T&5RUF9|Y;tugRAL6)LZ}28hxzoJOJ35tsv%Jecc+Q>Y
zeLmnr_lWm&@3D(wvSuXnFJ92hC%nw3n)j5?G~p_ryG4AV6Mn%IcAK~OvIV5H8@$d`
zc3po$7%h#ZwRD!=Ww4Bv$=-08EvsdbZ1%Qgw;b+0%W1hRw@YpRH=ob)N&zcqg{_Dy
zYQ?O$m2f4kw3U)FR@TZ{c~{XYS!Ju@@>q~nlNwgjYFiyw*Xmh)Yv3AMV{0T$tf@7#
z=B}l+vewqdRke23*4j%4=_H+{t8|y15^TM!kAz4+36%j7CW9qhdkD7(8*0OCgd1g}
zZLCGQ@fM{sh_Z<mZIf+^n`YB3#$w%&f62#revTV&bCsWK^Zqlx$QDY1CE5~O>XzFI
zTV<<VlC86~vcWdmX4~Sn*>>A$yW9fXXM5#<9kj!C#2vL`cEV1&({|R*$T>T27wn?D
zY**}>{qFYIKV7zmqYnT806G8w0C=3^V_;xV;9_841JWD}5E=wnq5PMe3=E76Ak4_%
e45XQ$VxwRbjDk@x3P!;w7zLwXcmn`KnF5A0qgu-V
literal 0
HcmV?d00001
diff --git a/demo/F16/aerobench/visualize/plot.py b/demo/F16/aerobench/visualize/plot.py
new file mode 100644
index 00000000..c6901452
--- /dev/null
+++ b/demo/F16/aerobench/visualize/plot.py
@@ -0,0 +1,294 @@
+'''
+Stanley Bak
+Python code for F-16 animation video output
+'''
+
+import math
+import os
+
+import numpy as np
+from numpy import rad2deg
+
+import matplotlib
+import matplotlib.animation as animation
+import matplotlib.pyplot as plt
+
+from aerobench.util import get_state_names, StateIndex
+
+def get_script_path(filename=__file__):
+ '''get the path this script'''
+ return os.path.dirname(os.path.realpath(filename))
+
+def init_plot():
+ 'initialize plotting style'
+
+ matplotlib.use('TkAgg') # set backend
+
+ parent = get_script_path()
+ p = os.path.join(parent, 'bak_matplotlib.mlpstyle')
+
+ plt.style.use(['bmh', p])
+
+def plot_overhead(run_sim_result, waypoints=None, llc=None):
+ '''altitude over time plot from run_f16_sum result object
+
+ note: call plt.show() afterwards to have plot show up
+ '''
+
+ init_plot()
+
+ res = run_sim_result
+ fig = plt.figure(figsize=(7, 5))
+
+ ax = fig.add_subplot(1, 1, 1)
+
+ full_states = res['states']
+
+ if llc is not None:
+ num_vars = len(get_state_names()) + llc.get_num_integrators()
+ num_aircraft = full_states[0, :].size // num_vars
+ else:
+ num_vars = full_states[0, :].size
+ num_aircraft = 1
+
+ for i in range(num_aircraft):
+ states = full_states[:, i*num_vars:(i+1)*num_vars]
+
+ ys = states[:, StateIndex.POSN] # 9: n/s position (ft)
+ xs = states[:, StateIndex.POSE] # 10: e/w position (ft)
+
+ ax.plot(xs, ys, '-')
+
+ label = 'Start' if i == 0 else None
+ ax.plot([xs[0]], [ys[1]], 'k*', ms=8, label=label)
+
+ if waypoints is not None:
+ xs = [wp[0] for wp in waypoints]
+ ys = [wp[1] for wp in waypoints]
+
+ ax.plot(xs, ys, 'ro', label='Waypoints')
+
+ ax.set_ylabel('North / South Position (ft)')
+ ax.set_xlabel('East / West Position (ft)')
+
+ ax.set_title('Overhead Plot')
+
+ ax.axis('equal')
+
+ ax.legend()
+
+ plt.tight_layout()
+
+def plot_attitude(run_sim_result, title='Attitude History', skip_yaw=True, figsize=(7, 5)):
+ 'plot a single variable over time'
+
+ init_plot()
+
+ res = run_sim_result
+ fig = plt.figure(figsize=figsize)
+
+ ax = fig.add_subplot(1, 1, 1)
+ ax.ticklabel_format(useOffset=False)
+
+ times = res['times']
+ states = res['states']
+
+ indices = [StateIndex.PHI, StateIndex.THETA, StateIndex.PSI, StateIndex.P, StateIndex.Q, StateIndex.R]
+ labels = ['Roll (Phi)', 'Pitch (Theta)', 'Yaw (Psi)', 'Roll Rate (P)', 'Pitch Rate (Q)', 'Yaw Rate (R)']
+ colors = ['r-', 'g-', 'b-', 'r--', 'g--', 'b--']
+
+ rad_to_deg_factor = 180 / math.pi
+
+ for index, label, color in zip(indices, labels, colors):
+ if skip_yaw and index == StateIndex.PSI:
+ continue
+
+ ys = states[:, index] # 11: altitude (ft)
+
+ ax.plot(times, ys * rad_to_deg_factor, color, label=label)
+
+ ax.set_ylabel('Attitudes & Rates (deg, deg/s)')
+ ax.set_xlabel('Time (sec)')
+
+ if title is not None:
+ ax.set_title(title)
+
+ ax.legend()
+ plt.tight_layout()
+
+def plot_outer_loop(run_sim_result, title='Outer Loop Controls'):
+ 'plot a single variable over time'
+
+ init_plot()
+
+ res = run_sim_result
+ assert 'u_list' in res, "Simulation must be run with extended_states=True"
+ fig = plt.figure(figsize=(7, 5))
+
+ ax = fig.add_subplot(1, 1, 1)
+ ax.ticklabel_format(useOffset=False)
+
+ times = res['times']
+ u_list = res['u_list']
+ ps_list = res['ps_list']
+ nz_list = res['Nz_list']
+ ny_r_list = res['Ny_r_list']
+
+ # u is: throt, ele, ail, rud, Nz_ref, ps_ref, Ny_r_ref
+ # u_ref is: Nz, ps, Ny + r, throttle
+ ys_list = []
+
+ ys_list.append(nz_list)
+ ys_list.append([u[4] for u in u_list])
+
+ ys_list.append(ps_list)
+ ys_list.append([u[5] for u in u_list])
+
+ ys_list.append(ny_r_list)
+ ys_list.append([u[6] for u in u_list])
+
+ # throttle reference is not included... although it's just a small offset so probably less important
+ ys_list.append([u[0] for u in u_list])
+
+ labels = ['N_z', 'N_z,ref', 'P_s', 'P_s,ref', 'N_yr', 'N_yr,ref', 'Throttle']
+ colors = ['r', 'r', 'lime', 'lime', 'b', 'b', 'c']
+
+ for i, (ys, label, color) in enumerate(zip(ys_list, labels, colors)):
+ lt = '-' if i % 2 == 0 else ':'
+ lw = 1 if i % 2 == 0 else 3
+
+ ax.plot(times, ys, lt, lw=lw, color=color, label=label)
+
+ ax.set_ylabel('Autopilot (deg & percent)')
+ ax.set_xlabel('Time (sec)')
+
+ if title is not None:
+ ax.set_title(title)
+
+ ax.legend()
+ plt.tight_layout()
+
+def plot_inner_loop(run_sim_result, title='Inner Loop Controls'):
+ 'plot inner loop controls over time'
+
+ init_plot()
+
+ res = run_sim_result
+ assert 'u_list' in res, "Simulation must be run with extended_states=True"
+ fig = plt.figure(figsize=(7, 5))
+
+ ax = fig.add_subplot(1, 1, 1)
+ ax.ticklabel_format(useOffset=False)
+
+ times = res['times']
+ u_list = res['u_list']
+
+ # u is throt, ele, ail, rud, Nz_ref, ps_ref, Ny_r_ref
+ ys_list = []
+
+ rad_to_deg_factor = 180 / math.pi
+
+ for i in range(4):
+ factor = 1.0 if i == 0 else rad_to_deg_factor
+ ys_list.append([u[i] * factor for u in u_list])
+
+ labels = ['Throttle', 'Elevator', 'Aileron', 'Rudder']
+ colors = ['b-', 'r-', '#FFA500', 'm-']
+
+ for ys, label, color in zip(ys_list, labels, colors):
+ ax.plot(times, ys, color, label=label)
+
+ ax.set_ylabel('Controls (deg & percent)')
+ ax.set_xlabel('Time (sec)')
+
+ if title is not None:
+ ax.set_title(title)
+
+ ax.legend()
+ plt.tight_layout()
+
+def plot_single(run_sim_result, state_name, title=None):
+ 'plot a single variable over time'
+
+ init_plot()
+
+ res = run_sim_result
+ fig = plt.figure(figsize=(7, 5))
+
+ ax = fig.add_subplot(1, 1, 1)
+ ax.ticklabel_format(useOffset=False)
+
+ times = res['times']
+ states = res['states']
+
+ index = get_state_names().index(state_name)
+ ys = states[:, index] # 11: altitude (ft)
+
+ ax.plot(times, ys, '-')
+
+ ax.set_ylabel(state_name)
+ ax.set_xlabel('Time')
+
+ if title is not None:
+ ax.set_title(title)
+
+ plt.tight_layout()
+
+def plot_altitude(run_sim_result):
+ '''altitude over time plot from run_f16_sum result object
+
+ note: call plt.show() afterwards to have plot show up
+ '''
+
+ plot_single(run_sim_result, 'alt')
+
+
+def plot2d(filename, times, plot_data_list):
+ '''plot state variables in 2d
+
+ plot data list of is a list of (values_list, var_data),
+ where values_list is an 2-d array, the first is time step, the second is a state vector
+ and each var_data is a list of tuples: (state_index, label)
+ '''
+
+ num_plots = sum([len(var_data) for _, var_data in plot_data_list])
+
+ fig = plt.figure(figsize=(7, 5))
+
+ for plot_index in range(num_plots):
+ ax = fig.add_subplot(num_plots, 1, plot_index + 1)
+ ax.tick_params(axis='both', which='major', labelsize=16)
+
+ sum_plots = 0
+ states = None
+ state_var_data = None
+
+ for values_list, var_data in plot_data_list:
+ if plot_index < sum_plots + len(var_data):
+ states = values_list
+ state_var_data = var_data
+ break
+
+ sum_plots += len(var_data)
+
+ state_index, label = state_var_data[plot_index - sum_plots]
+
+ if state_index == 0 and isinstance(states[0], float): # state is just a single number
+ ys = states
+ else:
+ ys = [state[state_index] for state in states]
+
+ ax.plot(times, ys, '-')
+
+ ax.set_ylabel(label, fontsize=16)
+
+ # last one gets an x axis label
+ if plot_index == num_plots - 1:
+ ax.set_xlabel('Time', fontsize=16)
+
+ plt.tight_layout()
+
+ if filename is not None:
+ plt.savefig(filename, bbox_inches='tight')
+ else:
+ plt.show()
--
GitLab