From 0524bfa6e8e125a58db45ad4f99b81189baf1a62 Mon Sep 17 00:00:00 2001 From: Luke Campagnola Date: Thu, 22 May 2014 01:22:12 -0400 Subject: [PATCH] Added new demos: - Relativity simulator - Optics simulator - Mechanical chain simulator --- examples/__main__.py | 5 + examples/optics/__init__.py | 1 + examples/optics/pyoptic.py | 582 ++++++++++++++++++++++++++ examples/optics/schott_glasses.csv.gz | Bin 0 -> 37232 bytes examples/optics_demos.py | 170 ++++++++ examples/relativity | 1 + examples/relativity_demo.py | 23 + examples/verlet_chain/__init__.py | 1 + examples/verlet_chain/chain.py | 110 +++++ examples/verlet_chain/make | 3 + examples/verlet_chain/maths.so | Bin 0 -> 8017 bytes examples/verlet_chain/relax.c | 48 +++ examples/verlet_chain/relax.py | 23 + examples/verlet_chain_demo.py | 111 +++++ 14 files changed, 1078 insertions(+) create mode 100644 examples/optics/__init__.py create mode 100644 examples/optics/pyoptic.py create mode 100644 examples/optics/schott_glasses.csv.gz create mode 100644 examples/optics_demos.py create mode 160000 examples/relativity create mode 100644 examples/relativity_demo.py create mode 100644 examples/verlet_chain/__init__.py create mode 100644 examples/verlet_chain/chain.py create mode 100755 examples/verlet_chain/make create mode 100755 examples/verlet_chain/maths.so create mode 100644 examples/verlet_chain/relax.c create mode 100644 examples/verlet_chain/relax.py create mode 100644 examples/verlet_chain_demo.py diff --git a/examples/__main__.py b/examples/__main__.py index e972c60a..6c2e8b97 100644 --- a/examples/__main__.py +++ b/examples/__main__.py @@ -32,6 +32,11 @@ examples = OrderedDict([ ('Auto-range', 'PlotAutoRange.py'), ('Remote Plotting', 'RemoteSpeedTest.py'), ('HDF5 big data', 'hdf5.py'), + ('Demos', OrderedDict([ + ('Optics', 'optics_demos.py'), + ('Special relativity', 'relativity_demo.py'), + ('Verlet chain', 'verlet_chain_demo.py'), + ])), ('GraphicsItems', OrderedDict([ ('Scatter Plot', 'ScatterPlot.py'), #('PlotItem', 'PlotItem.py'), diff --git a/examples/optics/__init__.py b/examples/optics/__init__.py new file mode 100644 index 00000000..577c24da --- /dev/null +++ b/examples/optics/__init__.py @@ -0,0 +1 @@ +from pyoptic import * \ No newline at end of file diff --git a/examples/optics/pyoptic.py b/examples/optics/pyoptic.py new file mode 100644 index 00000000..486f653d --- /dev/null +++ b/examples/optics/pyoptic.py @@ -0,0 +1,582 @@ +# -*- coding: utf-8 -*- +from PyQt4 import QtGui, QtCore +import pyqtgraph as pg +#from pyqtgraph.canvas import Canvas, CanvasItem +import numpy as np +import csv, gzip, os +from pyqtgraph import Point + +class GlassDB: + """ + Database of dispersion coefficients for Schott glasses + + Corning 7980 + """ + def __init__(self, fileName='schott_glasses.csv'): + path = os.path.dirname(__file__) + fh = gzip.open(os.path.join(path, 'schott_glasses.csv.gz'), 'rb') + r = csv.reader(fh.readlines()) + lines = [x for x in r] + self.data = {} + header = lines[0] + for l in lines[1:]: + info = {} + for i in range(1, len(l)): + info[header[i]] = l[i] + self.data[l[0]] = info + self.data['Corning7980'] = { ## Thorlabs UV fused silica--not in schott catalog. + 'B1': 0.68374049400, + 'B2': 0.42032361300, + 'B3': 0.58502748000, + 'C1': 0.00460352869, + 'C2': 0.01339688560, + 'C3': 64.49327320000, + 'TAUI25/250': 0.95, ## transmission data is fabricated, but close. + 'TAUI25/1400': 0.98, + } + + for k in self.data: + self.data[k]['ior_cache'] = {} + + + def ior(self, glass, wl): + """ + Return the index of refraction for *glass* at wavelength *wl*. + + The *glass* argument must be a key in self.data. + """ + info = self.data[glass] + cache = info['ior_cache'] + if wl not in cache: + B = map(float, [info['B1'], info['B2'], info['B3']]) + C = map(float, [info['C1'], info['C2'], info['C3']]) + w2 = (wl/1000.)**2 + n = np.sqrt(1.0 + (B[0]*w2 / (w2-C[0])) + (B[1]*w2 / (w2-C[1])) + (B[2]*w2 / (w2-C[2]))) + cache[wl] = n + return cache[wl] + + def transmissionCurve(self, glass): + data = self.data[glass] + keys = [int(x[7:]) for x in data.keys() if 'TAUI25' in x] + keys.sort() + curve = np.empty((2,len(keys))) + for i in range(len(keys)): + curve[0][i] = keys[i] + key = 'TAUI25/%d' % keys[i] + val = data[key] + if val == '': + val = 0 + else: + val = float(val) + curve[1][i] = val + return curve + + +GLASSDB = GlassDB() + + +def wlPen(wl): + """Return a pen representing the given wavelength""" + l1 = 400 + l2 = 700 + hue = np.clip(((l2-l1) - (wl-l1)) * 0.8 / (l2-l1), 0, 0.8) + val = 1.0 + if wl > 700: + val = 1.0 * (((700-wl)/700.) + 1) + elif wl < 400: + val = wl * 1.0/400. + #print hue, val + color = pg.hsvColor(hue, 1.0, val) + pen = pg.mkPen(color) + return pen + + +class ParamObj: + # Just a helper for tracking parameters and responding to changes + def __init__(self): + self.__params = {} + + def __setitem__(self, item, val): + self.setParam(item, val) + + def setParam(self, param, val): + self.setParams(**{param:val}) + + def setParams(self, **params): + """Set parameters for this optic. This is a good function to override for subclasses.""" + self.__params.update(params) + self.paramStateChanged() + + def paramStateChanged(self): + pass + + def __getitem__(self, item): + return self.getParam(item) + + def getParam(self, param): + return self.__params[param] + + +class Optic(pg.GraphicsObject, ParamObj): + + sigStateChanged = QtCore.Signal() + + + def __init__(self, gitem, **params): + ParamObj.__init__(self) + pg.GraphicsObject.__init__(self) #, [0,0], [1,1]) + + self.gitem = gitem + self.surfaces = gitem.surfaces + gitem.setParentItem(self) + + self.roi = pg.ROI([0,0], [1,1]) + self.roi.addRotateHandle([1, 1], [0.5, 0.5]) + self.roi.setParentItem(self) + + defaults = { + 'pos': Point(0,0), + 'angle': 0, + } + defaults.update(params) + self._ior_cache = {} + self.roi.sigRegionChanged.connect(self.roiChanged) + self.setParams(**defaults) + + def updateTransform(self): + self.resetTransform() + self.setPos(0, 0) + self.translate(Point(self['pos'])) + self.rotate(self['angle']) + + def setParam(self, param, val): + ParamObj.setParam(self, param, val) + + def paramStateChanged(self): + """Some parameters of the optic have changed.""" + # Move graphics item + self.gitem.setPos(Point(self['pos'])) + self.gitem.resetTransform() + self.gitem.rotate(self['angle']) + + # Move ROI to match + try: + self.roi.sigRegionChanged.disconnect(self.roiChanged) + br = self.gitem.boundingRect() + o = self.gitem.mapToParent(br.topLeft()) + self.roi.setAngle(self['angle']) + self.roi.setPos(o) + self.roi.setSize([br.width(), br.height()]) + finally: + self.roi.sigRegionChanged.connect(self.roiChanged) + + self.sigStateChanged.emit() + + def roiChanged(self, *args): + pos = self.roi.pos() + # rotate gitem temporarily so we can decide where it will need to move + self.gitem.resetTransform() + self.gitem.rotate(self.roi.angle()) + br = self.gitem.boundingRect() + o1 = self.gitem.mapToParent(br.topLeft()) + self.setParams(angle=self.roi.angle(), pos=pos + (self.gitem.pos() - o1)) + + def boundingRect(self): + return QtCore.QRectF() + + def paint(self, p, *args): + pass + + def ior(self, wavelength): + return GLASSDB.ior(self['glass'], wavelength) + + + +class Lens(Optic): + def __init__(self, **params): + defaults = { + 'dia': 25.4, ## diameter of lens + 'r1': 50., ## positive means convex, use 0 for planar + 'r2': 0, ## negative means convex + 'd': 4.0, + 'glass': 'N-BK7', + 'reflect': False, + } + defaults.update(params) + d = defaults.pop('d') + defaults['x1'] = -d/2. + defaults['x2'] = d/2. + + gitem = CircularSolid(brush=(100, 100, 130, 100), **defaults) + Optic.__init__(self, gitem, **defaults) + + def propagateRay(self, ray): + """Refract, reflect, absorb, and/or scatter ray. This function may create and return new rays""" + + """ + NOTE:: We can probably use this to compute refractions faster: (from GLSL 120 docs) + + For the incident vector I and surface normal N, and the + ratio of indices of refraction eta, return the refraction + vector. The result is computed by + k = 1.0 - eta * eta * (1.0 - dot(N, I) * dot(N, I)) + if (k < 0.0) + return genType(0.0) + else + return eta * I - (eta * dot(N, I) + sqrt(k)) * N + The input parameters for the incident vector I and the + surface normal N must already be normalized to get the + desired results. eta == ratio of IORs + + + For reflection: + For the incident vector I and surface orientation N, + returns the reflection direction: + I – 2 ∗ dot(N, I) ∗ N + N must already be normalized in order to achieve the + desired result. + """ + + + + iors = [self.ior(ray['wl']), 1.0] + for i in [0,1]: + surface = self.surfaces[i] + ior = iors[i] + p1, ai = surface.intersectRay(ray) + #print "surface intersection:", p1, ai*180/3.14159 + #trans = self.sceneTransform().inverted()[0] * surface.sceneTransform() + #p1 = trans.map(p1) + if p1 is None: + ray.setEnd(None) + break + p1 = surface.mapToItem(ray, p1) + + #print "adjusted position:", p1 + #ior = self.ior(ray['wl']) + rd = ray['dir'] + a1 = np.arctan2(rd[1], rd[0]) + ar = a1 - ai + np.arcsin((np.sin(ai) * ray['ior'] / ior)) + #print [x for x in [a1, ai, (np.sin(ai) * ray['ior'] / ior), ar]] + #print ai, np.sin(ai), ray['ior'], ior + ray.setEnd(p1) + dp = Point(np.cos(ar), np.sin(ar)) + #p2 = p1+dp + #p1p = self.mapToScene(p1) + #p2p = self.mapToScene(p2) + #dpp = Point(p2p-p1p) + ray = Ray(parent=ray, ior=ior, dir=dp) + return [ray] + + +class Mirror(Optic): + def __init__(self, **params): + defaults = { + 'r1': 0, + 'r2': 0, + 'd': 0.01, + } + defaults.update(params) + d = defaults.pop('d') + defaults['x1'] = -d/2. + defaults['x2'] = d/2. + gitem = CircularSolid(brush=(100,100,100,255), **defaults) + Optic.__init__(self, gitem, **defaults) + + def propagateRay(self, ray): + """Refract, reflect, absorb, and/or scatter ray. This function may create and return new rays""" + + surface = self.surfaces[0] + p1, ai = surface.intersectRay(ray) + if p1 is not None: + p1 = surface.mapToItem(ray, p1) + rd = ray['dir'] + a1 = np.arctan2(rd[1], rd[0]) + ar = a1 + np.pi - 2*ai + ray.setEnd(p1) + dp = Point(np.cos(ar), np.sin(ar)) + ray = Ray(parent=ray, dir=dp) + else: + ray.setEnd(None) + return [ray] + + +class CircularSolid(pg.GraphicsObject, ParamObj): + """GraphicsObject with two circular or flat surfaces.""" + def __init__(self, pen=None, brush=None, **opts): + """ + Arguments for each surface are: + x1,x2 - position of center of _physical surface_ + r1,r2 - radius of curvature + d1,d2 - diameter of optic + """ + defaults = dict(x1=-2, r1=100, d1=25.4, x2=2, r2=100, d2=25.4) + defaults.update(opts) + ParamObj.__init__(self) + self.surfaces = [CircleSurface(defaults['r1'], defaults['d1']), CircleSurface(-defaults['r2'], defaults['d2'])] + pg.GraphicsObject.__init__(self) + for s in self.surfaces: + s.setParentItem(self) + + if pen is None: + self.pen = pg.mkPen((220,220,255,200), width=1, cosmetic=True) + else: + self.pen = pg.mkPen(pen) + + if brush is None: + self.brush = pg.mkBrush((230, 230, 255, 30)) + else: + self.brush = pg.mkBrush(brush) + + self.setParams(**defaults) + + def paramStateChanged(self): + self.updateSurfaces() + + def updateSurfaces(self): + self.surfaces[0].setParams(self['r1'], self['d1']) + self.surfaces[1].setParams(-self['r2'], self['d2']) + self.surfaces[0].setPos(self['x1'], 0) + self.surfaces[1].setPos(self['x2'], 0) + + self.path = QtGui.QPainterPath() + self.path.connectPath(self.surfaces[0].path.translated(self.surfaces[0].pos())) + self.path.connectPath(self.surfaces[1].path.translated(self.surfaces[1].pos()).toReversed()) + self.path.closeSubpath() + + def boundingRect(self): + return self.path.boundingRect() + + def shape(self): + return self.path + + def paint(self, p, *args): + p.setRenderHints(p.renderHints() | p.Antialiasing) + p.setPen(self.pen) + p.fillPath(self.path, self.brush) + p.drawPath(self.path) + + +class CircleSurface(pg.GraphicsObject): + def __init__(self, radius=None, diameter=None): + """center of physical surface is at 0,0 + radius is the radius of the surface. If radius is None, the surface is flat. + diameter is of the optic's edge.""" + pg.GraphicsObject.__init__(self) + + self.r = radius + self.d = diameter + self.mkPath() + + def setParams(self, r, d): + self.r = r + self.d = d + self.mkPath() + + def mkPath(self): + self.prepareGeometryChange() + r = self.r + d = self.d + h2 = d/2. + self.path = QtGui.QPainterPath() + if r == 0: ## flat surface + self.path.moveTo(0, h2) + self.path.lineTo(0, -h2) + else: + ## half-height of surface can't be larger than radius + h2 = min(h2, abs(r)) + + #dx = abs(r) - (abs(r)**2 - abs(h2)**2)**0.5 + #p.moveTo(-d*w/2.+ d*dx, d*h2) + arc = QtCore.QRectF(0, -r, r*2, r*2) + #self.surfaces.append((arc.center(), r, h2)) + a1 = np.arcsin(h2/r) * 180. / np.pi + a2 = -2*a1 + a1 += 180. + self.path.arcMoveTo(arc, a1) + self.path.arcTo(arc, a1, a2) + #if d == -1: + #p1 = QtGui.QPainterPath() + #p1.addRect(arc) + #self.paths.append(p1) + self.h2 = h2 + + def boundingRect(self): + return self.path.boundingRect() + + def paint(self, p, *args): + return ## usually we let the optic draw. + #p.setPen(pg.mkPen('r')) + #p.drawPath(self.path) + + def intersectRay(self, ray): + ## return the point of intersection and the angle of incidence + #print "intersect ray" + h = self.h2 + r = self.r + p, dir = ray.currentState(relativeTo=self) # position and angle of ray in local coords. + #print " ray: ", p, dir + p = p - Point(r, 0) ## move position so center of circle is at 0,0 + #print " adj: ", p, r + + if r == 0: + #print " flat" + if dir[0] == 0: + y = 0 + else: + y = p[1] - p[0] * dir[1]/dir[0] + if abs(y) > h: + return None, None + else: + return (Point(0, y), np.arctan2(dir[1], dir[0])) + else: + #print " curve" + ## find intersection of circle and line (quadratic formula) + dx = dir[0] + dy = dir[1] + dr = (dx**2 + dy**2) ** 0.5 + D = p[0] * (p[1]+dy) - (p[0]+dx) * p[1] + idr2 = 1.0 / dr**2 + disc = r**2 * dr**2 - D**2 + if disc < 0: + return None, None + disc2 = disc**0.5 + if dy < 0: + sgn = -1 + else: + sgn = 1 + + + br = self.path.boundingRect() + x1 = (D*dy + sgn*dx*disc2) * idr2 + y1 = (-D*dx + abs(dy)*disc2) * idr2 + if br.contains(x1+r, y1): + pt = Point(x1, y1) + else: + x2 = (D*dy - sgn*dx*disc2) * idr2 + y2 = (-D*dx - abs(dy)*disc2) * idr2 + pt = Point(x2, y2) + if not br.contains(x2+r, y2): + return None, None + raise Exception("No intersection!") + + norm = np.arctan2(pt[1], pt[0]) + if r < 0: + norm += np.pi + #print " norm:", norm*180/3.1415 + dp = p - pt + #print " dp:", dp + ang = np.arctan2(dp[1], dp[0]) + #print " ang:", ang*180/3.1415 + #print " ai:", (ang-norm)*180/3.1415 + + #print " intersection:", pt + return pt + Point(r, 0), ang-norm + + +class Ray(pg.GraphicsObject, ParamObj): + """Represents a single straight segment of a ray""" + + sigStateChanged = QtCore.Signal() + + def __init__(self, **params): + ParamObj.__init__(self) + defaults = { + 'ior': 1.0, + 'wl': 500, + 'end': None, + 'dir': Point(1,0), + } + self.params = {} + pg.GraphicsObject.__init__(self) + self.children = [] + parent = params.get('parent', None) + if parent is not None: + defaults['start'] = parent['end'] + defaults['wl'] = parent['wl'] + self['ior'] = parent['ior'] + self['dir'] = parent['dir'] + parent.addChild(self) + + defaults.update(params) + defaults['dir'] = Point(defaults['dir']) + self.setParams(**defaults) + self.mkPath() + + def clearChildren(self): + for c in self.children: + c.clearChildren() + c.setParentItem(None) + self.scene().removeItem(c) + self.children = [] + + def paramStateChanged(self): + pass + + def addChild(self, ch): + self.children.append(ch) + ch.setParentItem(self) + + def currentState(self, relativeTo=None): + pos = self['start'] + dir = self['dir'] + if relativeTo is None: + return pos, dir + else: + trans = self.itemTransform(relativeTo)[0] + p1 = trans.map(pos) + p2 = trans.map(pos + dir) + return Point(p1), Point(p2-p1) + + + def setEnd(self, end): + self['end'] = end + self.mkPath() + + def boundingRect(self): + return self.path.boundingRect() + + def paint(self, p, *args): + #p.setPen(pg.mkPen((255,0,0, 150))) + p.setRenderHints(p.renderHints() | p.Antialiasing) + p.setCompositionMode(p.CompositionMode_Plus) + p.setPen(wlPen(self['wl'])) + p.drawPath(self.path) + + def mkPath(self): + self.prepareGeometryChange() + self.path = QtGui.QPainterPath() + self.path.moveTo(self['start']) + if self['end'] is not None: + self.path.lineTo(self['end']) + else: + self.path.lineTo(self['start']+500*self['dir']) + + +def trace(rays, optics): + if len(optics) < 1 or len(rays) < 1: + return + for r in rays: + r.clearChildren() + o = optics[0] + r2 = o.propagateRay(r) + trace(r2, optics[1:]) + +class Tracer(QtCore.QObject): + """ + Simple ray tracer. + + Initialize with a list of rays and optics; + calling trace() will cause rays to be extended by propagating them through + each optic in sequence. + """ + def __init__(self, rays, optics): + QtCore.QObject.__init__(self) + self.optics = optics + self.rays = rays + for o in self.optics: + o.sigStateChanged.connect(self.trace) + self.trace() + + def trace(self): + trace(self.rays, self.optics) + diff --git a/examples/optics/schott_glasses.csv.gz b/examples/optics/schott_glasses.csv.gz new file mode 100644 index 0000000000000000000000000000000000000000..8df4ae14a38b916d5591749d049804097f8e86b8 GIT binary patch literal 37232 zcmV)HK)t^oiwFpOZhcb#19M|&Z*+8DXKZ0}b7gZbV{>)@?Y-NwB-eE$_^z)g<4a3z zc~0&7<&hUbYJt|H5jZwB{RR;w3N|3X20+R5Uo$_@UpC9#*N)h6qE4ZpR!coCfUL|s zdC!RLd)<8b`irkV|L(iT*Izuoe)sr?FCPE!?(xfSzWU~m-@f_nH($Ja{31SnksiOu zk6*^eFVo|f`SGjp_*H!TDn0%@JpMdBzWV(4fBVJfj~H6``uC4_PhWj$@AErd{fplD z=0i{Y?E7!uefQm)Z~pN8m%sh|EBos2{_gL-{OySo3+?K4#tS=S$mpZ%fC`I?kJ=9XHRHmGfmKA5SRfJCsn*o5#I*o;S~r$mciC z?RZ(P=e@ZcH_wmA@w~Ymi}NGW`4Q>-h;)AS@CvP72uEU9&=SSr88?VEOo9Ay`hZ8r? zk64EjH@9PPendJyBAp+R+I|clzxdT1_s<@me*Miie|YovpMU$s*TC4H{^h5C z_W1eVe*5`%U)uAwyZ^$T`grr^Yb=N7f?HS&&ylxkn)c0c%ND|Oz^(6b@f|L{#lfs1#z_=FdH%I|RT_qg~H7j}yDK{+OUPMYnx^Tr7JHmQu~{P_Bl zC;KRSmnVGK&)&TG>o?y$zWV5OTYE*@rB>|-QyF#_2uL5b?ZOi#c$Sa zd&BS6ZF__Nux{I19M7>gIqtW&(fz-|M}L3b{BGU+?(wV7zkm1m_T5*{Z~pq*FMs>)@qhl!(_h&_vWuVp=JA(*@$^@J@=1Ec zXTw=2ye?rpa(m7lmoYvvwhMbCmzr~~^AVm?n4#y`o86wd=h|aSb~}VQdoOwPh`&k? zrH0=6BgXhVM(Yh9QRAOKg@zl=Zj8rMdhX-TpX{TieTdy<&$k!bjrmB=HQ0@4H%fTq zXTwpoyX@U<6R~fL`wV>c4&hPol%@wV9tgN^?6(>&BHkt2GN|^?mK45g9r#^sRJ$Dd zUA$+`v0e7tJ^THImzkH5@4)-{o_1Jp+4yL|_n!&R=>xY@Hpc62>AHH!##_> z3cuEzAAPxhIlZ=SUo>uCj+c$k54c@3UKRM|&aY*Mj{lMY#tyL^Lp%P(w!vbDyYg*o z?fg}b@to{*rU(B{>49B1>1f-DAH}v*T+cdiJ!{N7opI{4rapPyI-cKRN}bm=WL*y~ z+SqF2HO|({MQqo$G3y%JvR&J^;OlDON^GqqzTYUd20!v4_P`ZpH`2h5N;%+Ql`*i$ z)69c+_P~Xcc-Ug$V@Y@FDwUG?OunE-oy-u2j;dc(UHPFFj0_(3TXHMU*$(kn?>_$m z+xYjNfBE&7U;pmS?`Rgj|EE8^d!+O>y(M6p;rK6q^AnPL#3}1q>crSu@E?!l9VWU= z8)W!QE``aXGr7jK4pO;7hT93>BBm2#{$W2B$@^ElqdIS2j<*w@Q(OHB$^U;s^3$Ju zQhC5~;SUVmJ|&N5J6yRxn;AJ!978EIC~@-HTnx>Qu#AmCxDGpBq44XK+FQc^4Z_uq zNF26G(sRrX;Z_CVX3K{VHsmHg)(GziY=+?wu7F`UAY5=xh|d$Q4IoWBwjILpfNYJ7-=-T=u5z=mz1=j~`qZv&8Mg zUIXdsG&_%v2ar1-bvSD6Xl5N? zY}Ye$|L%ootbG!tSTEepM^8le3^ViE@X71ZiFW()L}M^yX|bkpTeLuPX+nJ!x>M$^ z(P(ji;-sP%TP z=J@SN_|1eZti72en&4u`!|~f>f*qz4K&)hk$m~z-DO>FM)rc-aqa-4-B8SX>Nru%PtUv^`5i(E2a>s+53Yon z`}X0xNxmD}eardq#{i)2>H*2K*%*`DmL!rO#DHpt&rg641rGJY8#H7PEYl*PZ6NOh69TlMC3ct1sj< zZwngNoq(>kM0tNqg&Rd|dOT&l-U9ES8>8r_a$S!hlW7|qfjY8a%S}25o*3IJxR?r~ zldhMf5oN0d@R|j1`m9B^qK_!_@PwEfv4GumZ7Afm?O^t^s#)e_2m1p6&u`~j8NEOG zR6jQIU}HnE^ALb;yWG*8%#*G>LFPF`^E!C;Q;+5(WW{>6Omi^@g}YE{woK~8CL8Pk zm=j~a_Q^KJPIM0g8J%czumN~Ox0y?!yPG3Mbo&-S!|hf#!PY+U? z9z);ug0-i7*fih0iq|gd%kQu6sc%UCPbW0hCf2e7E^AhI(7Yxb8w~ z?6mDQ!OO*UEe>UF83t1A5!LplenJ|So;MtmgD0sX5erB&g^>jCS`to+g_8n5Ba5my z$v!O!(v^@_G}ofAd>CoZmxOh-?zLal0uXISHJZ3wc<|K*d<_;Z`NeZ-0(u_7+E z{mm{7WCht>{H|(2$LR9+?bKQv#nW4C3+j!FD?~FX26CEgJ1zEN!M;5LuN*oBIIxS#pRokeciLP)r=a~U&P#`?IT<2*{Ca*0r z=5>-wY%Wc{tyf%IgpKPMxDD75Cro>mD{g+FZ4-wpJtKd<-6(RN3o{U~=lDyhi^G-f z$^g=dOG2)Jf9iESo0-v%p$yM7SDI+ZL7!pahg@ezX%#N9Z`jYsEyJ=se)a)V`{p0N z{OYSWfBWvY-~9gFn=ilquYdgV+kf)MKmYE{H-Go$uRj0!4WRdU^Yu61zxftV{QbMH zzj*h>o4@S7xk?l^IJ5 z{An5PP#yd^IQWZ4)a|g&9yn)QbUZo%@93enzoNOT&Um{YUb?MWy9-N{{%FWbWe})c zIq2{armvpcfv1bPOq9y1+a|&+wjBOm&JaG%Z~;HVH|3>fP_U~n@G)WZ9y{zw_Hb?@ zz8!E2*aG*;&=3*9xdR|{uIaJ^*te@V7oA&$?Ws&I^4li!M1g~Ts{y3_Qanp*o zj6B;jUj8ou`uM}3w^+QsWkMh475XQ5UJ+$ac{(tSYm>%u!QJW!8%D^;czYI#u zh4-A8o9iC%u|M5SxjURZqM#a|5@CRr$4*>y}fXieT(=#JijN=%zk-P;H04 zTO?BLbM56HN_V%n{Vfy7hSX0~LWtZP>nH=*v;~Gs^H6em%L4>lnsZL;Y+I~3K3!ST zJk9%wUbg~-?FbXM1qqL8KKvsBK{KPm3tKF#8cyO7S7*m+ZjRNt%ySY}!_^~yihFpx z0%y->7e1nVgx~28Cr9~{7rqb^!apBed#!0V{+QtPwL;JrBLPg&XYsajQVD;_8JjfO1qnf)$_k&N9aa?*=xRysuHvdbp(MB=OQ9OgYGo@L-m zQzXsORSk5P2DYjvFHF6NoiFKcgFAjtfXVt`b}UYrTtwNw{`DvPp$yw`V8cWT^BI0q zMTG(hG9o#eOB9nFvfNtnrI@oVjpn&QJH8f5d-clff64fEF?2hU*}lH5f$BBpURLoc zCik)#P@^hb$MUkXr^o7%=-yzOk39H1(%l4g3+p9kxrzlYWge$PC?Uo3!J#PwfLU%yi;8>J-^> zSJa$DuBnmXWKxnd$R0CR__0*+*^J-?#;dwzTY>!? z^DE#8e9$hV3SF4&{?AXp`1vQQw)dboq$B<>9b}Tyu$$)TNi*9D3m1UF)HYNk`V1WQ zLt<*4j!(AjnQxth=VF`c!1)Unzrsl4eD>)+bQ&vXR+8ttNAa*3MJ(CnRn6*0bbNxf z9M$M?M91NHpnAH@i`ewFi{l=g!Ac*IC#t9#J+^@PlIt!T%D38;=ZClWmV2Cs%3<@P zwP(m#D`)NMg$JzwO}Ww)FyH-y*J60Vv7cYbJ!nE)esq)vP1c7p{;X4>QoBXWz-n>U z2oHUI;H$Tq9x4P>yr?5rW~;s1nSRjBCo@gp@)jLd#>h#raG1L8L|bDrh1W`)nK@NL z&5VP`!b7oeQSjCkgSR>FYn2H%l|x2QLd}jk@1F+qtTpWz`3QmLd!@jQla62Tn&6z1pNLe9zCHUkmRCx;l zh3k?TxSHiUv1hi?wj#O37KRIP@tHqF=NYvyGoV#b7)nsV908OMBF@Yh6;4y0ZcUS0 z(`hd>sG;T1}A|&hlz! zaE^ytyIHENnR}_We#G@iN>qm==2mcS$7P+8^HH4HC7hdED{t4{V8fkR3U_yd0@YgN zNeN7|%ZcTc8Hwz;Qypn4VD(}RixT$2=!VxT4_+%$c!;xW4tU9Soa_@;aknhU(WAHB zOyAjl8uI+M;o{PEU0_btGK+FSIYZsoZbZ@F@SANLT=nN}(9H*d@d zj|)jsHKIU!XFF$Hkn0Y)INgEm2k&Cgfvu{>%Jlie^FU-iY+^rXTC1X@@buv+5mnKN zxV*)oZtixf#_rng+(;pg_t$nJE}r-uLbTFSd}h9qG#qvWhZ8l*X2f3gnD={j*rgM8 z$_xv{^qU4lq|TkbsqlC$5h3OTGi{kWTYp#*C?|C$9ENI4$QvMMV5UlI5)PNAiUpU5 zk|JRsGnJ;OQMoCFo4d?XQ9Z|uo0ekfN-2}wpwzLJrjk(CrAjJQQOejUSzUFa)_ZFk z?8?d!7q7yKN?1cx8nk<1NM0%GyetU`cJjlI?6CY%`x@?Ry17f~8d_pQv#)EP{0zr= zLGQ43RolsEB8Eorn^*`b%!3~0&JLi(wA?f%%%|-%L^W|;9boVv$IL+PeYl?O)W^25fyEZo5)!n!K7?#A{&8(-n zySSN5)M``DO4W_`nLDA*94Dk^s(nglETq~~V0Ip1oE4z z?PxM;A|l$mHeI_@XmOOBt|GLWuOv2Zo%(ycXO-f*i{|WWtEHNp=x9cBGj>5NiI~;$ zks42k#ZH9Vbd#;hOfx?he^zC6L?|*RQ_m@? znq)b8vG=ImJ_whI$r_TsBFGqxJ7o@FLkv$kMRR{uK?djx$fT7YY#u@a-zza&gUy6{ znX?`Yar3Zt-J?5C<=9@^oXs@SZtaTuj2E`jIB-|JQtDLXMWm5FUXXeuyD*KxF{Ll}w;$V($_>1Lqe-EE& ztES3jPHLo}ZKO++>6yjganiDow2O32$pAI$Xom&=WWyWZ6Et0PDGvK5cg>axIzf|8 z1k*1{%h;mk-?fTUXkeq`Y96J?%SGhm*1PVX;vBX)MP?4g7k0TXaHeVWz8)ANUCZ|w zQ_r70xTPHdZ!i77WeNU!(bR^7)ReYJeN+Yj} z=LYBPDy`A(xK;^@pn}AI+N*8=0yno-okIh3GH3aWBdTUqw@16y`{ zeje&F%Tq|btZZC2C1mN#V&{la*VPlX$#)DpGujGTbL^+V{Z>9BD<$AnD*ftQV`WWp zaAKlU7m# zR6~WVjuJ|(c2;o-E7#sKC7tXEMFEX%UiP(mQu5IbT~+DGHWAy~?ZEQuf!|=xO2qLq zCwaDtP7ed^A~BY7@R;%2S^mw1-f~VV<)|&`u+PFVKY4JQcQ^Odeu58o z2JPQ>FZIYED{CaDK0j7Rk__^eDe|TuWQh*P8N>#IYNxqr^4#pCO>kCFu2=L&T$A5T zAO%%2OO<;JvM{fk09d;!A!zuKNfw-1^-8L8KB%PX{?YFX7uQuB=I(6+Bl>kyRu z9Nu&A-c<8IgW}dUw@*S@su+&XngZCUSvM#B-N1V^N7@xvbKeGRm%qf{!|LtAZkFIO=me< z@X%90X*x`EvwDVk&8vT}^yO36;Q;N`0qP~tNJ_(E85dE8ninZCl_|s`v2m?TY&Jm`xR%cvVg2Ev-b48I2Re zo;9DrV|p3uI>)H=tPKt(!vNVb`C#o}0Hqy2ECap8;BVBaVy*d1V0RwXxV zr;N79VYm_O2H+BJ6w6949l>r{tvexkusdF)cGRsYQs3)VyFN2E%>`wpL_%9je6ibV5NBa0d^+_^8JYpvB(CxF3Q7~qE_4DzG>-x|GqMT@gj1O*5XgLyn79rMV{9JY(5(}fMC*teot4{+kmnj=+^ z?P(%PkAp98^w_lM^^@Kgk8gFaQ>7xVwwdVG&lx_`eNV5Pl%spS9)zEM?OtERKICsJ z^SCPF$@=Zx=IXQKJX>v z-3umR`yy)dsfTUzI%Q_gK>V%E&fP4X*)+qqow-M6+l-oZ3pyAQSUmgQgh97SnP3uGNEK{iS~d(%D7rYLZ#T9U*SNLAzjC*{(8uh>EeujC0^U3tk8Ib(cf zxR~*txBcuB@q$C~49PHy1G_hH88>1(a%hP@f9wXENwN3@x& z3PzjWr`A|0fsl&mGG2!qvP%~Dh<-DN7RQ|4m!>|p9`d*^PLSlA`qX0gzO|-Yd=`=o zhohW;>qE{}9R9rvsc=O{F{!wGuma>j>#Q$|PVRc0(`CH{NJ+v`l}`2KYEFk~+o`d3 z`0Q60BDmSU0O6MN zb~NlrDZFhl+FxeC=d|NurZzLho?a+DBqd((SnZx4z`N?J{oPG9KSSWTL`u=9t8L2`|j10dqb8x<^;T73-jg z^#@yd<+O!NenHWhK=>cH!vTDx{d2ni__6zy&+ z9Nth;?$zL~Hgnt@GPIP{R8Yr4Pnd|TOoY08p)NXii) zllU>Th;n{GQCdVRl_rrQGw(n93DN%(AUcJ&P+)&N8J-W(V#c{k>a;M|WN!+xAh9>i zfmZKotP@RSY!0xP7|rPq5jyG&7VGjai?ODG9xhRwAX}k%z1Rmhb1C-ad*|vQa9JVn zq2ul94UM&hZ_TRhyibdy;9bOQGN(4C~)&(~kL@gHt?+bgkh@D)a!+<3*n-RtPwHq_5B zZ#^l?cQ(I(1YN_hQyGsO?m%Na)U>4UIXs8UqsiLi(mpJ!_Nd)8Hs`&Iup2x;jwy-q zS0&gIt5Q?3RqyU}qVK*&<9n4%`92jW*dpNIee^v^!!h$v9U_RU$D8;Y;=HF^1xHtK z5mPQRTZd$;lMq$JD7x}sx1twCLZztNJd#umtGmOWd@3y}iQCu~hHZ|@XLO$*Ie55> zH1|~NxYJ=?2-$fdZy`@(l6wEN{fp4?;Cofwdg&-FNX=Nm=e%~g?fwo%{8C9ld0e~h2} z8k6DH5-N?OzyW$fQKpu!mQLX` z+Se9`t~>2u>De_B2QyZsPyIby6LARe&`fmUwWh`1SpCG}&o>ajjwg?(Q*(i~>6(fz z$X2}TcxDfkOInuN`4Oo=#O-beEA$2<*qo5DmE(3KPB)QyUva=%c6f$A4 zM1@(UnOD=AZyHc=b5V6PCAP_h4mcvLKXw)9A&MWCz9A70I7S zg?N&cW>y+1$h>U4EWE_FK$g`Lg;t)rUj{zfd0Bax_?*EtdNmjq6P0J;N-2}-Lw6fHbq6khzVt-tLH_+vy+KRB%&2$kr2{mxp zBI`=0H3_Kp$+WdyPk}x*Z4SBX?=pPTI0RRVFZy^}{8VwQF38GsJ0wYTf1@Uh1p}gF zgc^K_&SM)#VW@e9NQJm0dQ&Yxjn?D`osgW?==nEA3hZoG`V0bG;csp=4DD0eWANN zpf;DKb3UuwDcU?}J&KlWE?I;hA4F0vkmbT&93x6F-MnIhbXA5`2LV;XZsN7B8wV<4#N7@qOnB#+R2D zcBA}!aiCPIqF21RElKA$o!!u#n7``3+|`=?d1-<>aYH_WzFRo(?<879;j6}c_4RGl zq9z88RQHsjA{n1quyW;KH)gLc%Q81tFLr?zTR$$1I$h0J!Vs~#y73UVBRV%Dc?6v_ zd;p6$-WD^i)w2h_vohdW5xG~hAyXu-2ioy?W26)80GVzGyJ<+tr4af#6hYWK4d`;S ze}E(qsXVDeq57DE#Yg5d+u~C#=_c1tI;B+?+sb}co#@ax+hr9LeENxNJ?khYh!$7; zEL^6OQx&Hj2#HIEVOFKIy7IjCMTgH#urC`CeybY{PE+UuAWm%Mu+M-oZ!U;j%Chfa z=_<1K07dC`4w!o5y2N|sI;K0b z%8?oUFEjs%2J8dqR$JWyh2h3*Kum(}jwx5;CaB-rFNy2NPl0B~%yfO3fZ(R3&9! zYZG;(2G@+h+ZB)G)~tuDLRB2s`gAO37ZB7YD(fRJunECqwIFuao`gzONr>5tvP(@C zHk-$vn0HOXf{(~FM>NMf_Us!_v7lUCf93`h$Rel503j|0aLTe0+H3%7FB;jK3fB;B zNh?7Xn)eCRoloqH#k3c!>62Hb#^OQ0N=F~V?^vOJb|<_osul2yUh+d4DJm+ zo1JwVD*nbe(dg4Tcuzh11iZst_~c0xLnYRiwpD%CU2G6f6~r4jKe{S)*)_IKhD*b{ zz8eeBbc-dimDG1%HrpD>vxx_jgh_5Pl|-?iuogOQap9a@x2@&bM}{=#2p7q<#H7JI zNi78}McwbR;^u;=Z#|1P3zZix*IdNN%6(mL%b5M}>TfTgr;!2XOfU5qwXRVn&{yu; zXc5O$ge@l;o04Om)ebv_q)Aj#(RFIlLv9jEEz@z!?maDwp#VG-r!UTObo!8?tc;Rv zU?@33xrxM*AM&6R6LEG%?4|e_y*Rw(It*nk8-|J00-`v&4f5eXL-Dc8-AB;eW_GoD zXq%_*d4SS6D@kQ@lymTz2dn{wyywjk$zLJBu?gcCt zZCW(27h_3^G>oe=G1uS`(!k2IG8oEm%6$%7IPKtgjh(SqsRWv>TU{mIJt@A;2-^9W zb*hr3@aTlNk}Op~j9aU~h2Q+)0@b~}cw>p~6wO`JBAzO6(pQ$&^3fm^nONHdk4f7# ziwi2mP(mCFUQjRS6@~akbFTga<-@FFulkZ0&mnrZo|w<2%Ca6^c~c4X>Xz4b70Ibp zS@7DtTTe`t>ypN`1|Zbdt9Om&_}4>>P}Bb}lg+!#0gC1bprg$VI3 zN;*9Oit#9qR%Oet0|_Jr+}!UO?k$(%y$I?}Kmy)=j_q`Kn~#l^pC0BVJrI2NWyw`o zZ~>i&J*yOGjo2j`^_x!wgL;Pq;Uel@Rh0>gj9b!LLQv7taL#xJ8LuNJh=2=aDy}FcNlC0sF-MfB- z^XV2BNmyR18L)-vf>3j)Fcu`(7Tuy>)-4rqB*;ehm>xiPM^&Bza%QKeo%@hA&CarJ z8c&$kC=0JHGrLt?Vo0fPWsa*H?xgIHS!~u$r~Bl>vA4v=zSK{39X$#zkQ%Z00XeVm z!xnqxydrbrkfF8Xk{6cWJBa}Bn9!a_%&aO!!SpkBlFeL4kw}3>55rcH-13Pw$vl?S zaV^-uppc3}ctc=f=JIgN)rOR8sXLa(bb+ZuQ1O^UklZZqBl?!-?T=O;QR1}i;7+`| z^TS7!csnnma`9Ag&cv-wC*l^CH~#c{GTXbRTRE|(HX>!IHfFbTx)qVnbN2{i$Ffsl z(&jCXd#L;xmT%#)``&xreXEW8e93U^9^x-U$@4u?Gg~20#6dYpMAq)a5XT#h_qSCWDa9zC+e_kI zg60@Y;9V4Tn?=%!6zOdr5I5V^6N9UtETuOmf+{5{%!?#iUFfZda+h5wPTRGQC2+mq zxTbuqStQEq1T@uiECE($`B)^IS>h>anCg&{1DOe+6iB!2l#9&j5n_oe5B2$CCPZ!; zL};j-Lz52^CM8Ik&iw1)t&_IX6f@2s{-VnX%pkkJ6SGWhk!&)lRV8W0E(vqywWVdA zgn2i~d!JKr+nNt9nd5eNlFR3Mo{zEzI3Db?$$3H)>0QJSm-P28-@HNR-qybN_wLW(k`tIN5Bbt2QqZYTt0oq>D*-)&!muRSd9*ZYX z--~7JnRar};@X&z;4&aDrBgT#`?jlyK5AEohJ;bO3#FQ|Qb)#A1U-jzYUWHDOHlZC zj}maSUiZOTexiA_K7?)JV%W(tonBm7cdi|p!X-9K&G?S=a zf#zTToOH5Y)K%C)4S;Wr2?PauDjKTj8k@~6DTz*Ids8KoAfQdy=Ofrxjy%vb$%c6K zy9DYWL9P8MO8^IotI4POwdX@{Ktf_t?tE?FW`OYtUik{3~eIM%y z$$>783%zkT986?LF=uX3uQN|xD?M*~74N-GCz=(H@{Fq3t@5SHPD+-$*70RbhciPY z&bYc5?|hZl7%BS(BOq&oIn^`NAE`Q%lI2oC0i$cnU_;owre`$=C-tkLYyv@s?hQo` zSAPr|J1|w19ZZTf9acb7Q?9Wz1WiHLphp?Go>Y+PIrzj7Q8PQF4t1?2sS+qphr#Um zwqOC(%4Auu-fRU+PUt2Y@NRqRLAb@7AV*mf*WR0@YLLV1y;dt3R=p;uXR1b5n%Xzj z7M)#(RKl2fIn!-5>i26J+b0Hri_(kzH@PIvKC9};h$lPA^tkYqu!P0mcK3$kXWx;>3uky5n7H8N5-1qI_w)X81WN0Zf(sPiE{Z4L2> z42$_y$C0fXN7=XR2`!h?NJrSaa@%f!PhWeS0_312V5d>56@S|MFmZ3RxaF}O2YjZY zUiU?LPP;-OMvBY~cz>_T|CgB{Ug3%|iW|%1M>JGze2=j)(!4iVY8q;pk_S(Gs@N&^ zQhpKpR_R1WrL=}HiY9)IwNVlmDMm=%pqeE47bwzEt%}E4b9soG8il5&mZokNm#9Lo zBhgee_h$q>-`AR5^X!8gp~pK@(bu($>*>v~Dmj_+$K3(c5;nExH5$QDIq|md;B6nqT(RQ)TfoK5Pc!#YUEKvX-zNr|Su(RA? ztlWH_4c-Is3JMM&@NkJo0yXM51E@)6u1t0`n1`y>)kMMWt+D2dd0BkiyZjBXMSzVC z;N%0fD&yDd`o;Cox}K9BF(wFf0d)wdv2U5GF1kgCLq#SR*04B9=2Rk(o~e?I+ZU)JPOv=+kKb`+#Cm&~b5b`%Z6W)bat zU|-Ltj7KnnXB)GDht1FE zn(W<8upHkYlbKyL08R2-Xi~C5A(N)M&gO_`^|v==xG>W|e@i7D6T8}NDbrLd#Dh`VyQYlD%4l(Fa}xR(2{?wJAsM%s9}?=SZ5ltpfQ~)Jl5IQFsX53rER)W&8DZ zM9lg%9_9m~ECAW#xk@AufvPJO}4lDG;8_8N<6tKrJ+qnmAAWQ4yV`BrcM=y+hJ#e zXN?DoL^Q)b(iW7yL3NRbbX1os>bkC;K;9wvT6glEM!uzJ?+N61#;p>7Op1WQHJQ95 zRNXc95tko^<7?2|i$h`io&vobcY;!%D7LRmE+37q!o6~@Y^Y$>S5nOd8ZngeBP5od@)e6sSk%s?(K2JY(v`EEdb-XnfF224 z{b`l~74#(#P*!A`eaoPL!l2+u;dCRrbQK;p0z-#SVJFgJlSyqKe?ZRQ?!{eAQ z-5k1ptFoAt{e^>FH*74H0+hh9f#!V!vhG!#?U{ShC~Y*0zCkdaH5?6~h~>h4EE_p; zH6y;(v~UkfU)vSrI`f`*YKf?P(;1o;wXwBR2J7&1Pr#u<(jImAZUe zX*Lg*SK#AHM`WPW;RNFtKK^+IC&5VjDoAeUW$K7WuJ{0hytOPOZ|gjU+A)X_KtpKW zsd&!Idh)yJ>tfJRz|}L{-k*Ks80qrEW28L=)t1qHIYye5o`D_(cr>M{*Vs0?G*pR5 zwmC(RF_4si(tBNUTisk^{W{zBEC@hu3I$tqBsX@nt#Xb_sXs;PxDw^9?~kDp(^jdc z9*HI&lK6;+96(Z{DIa=Tt6!mrX|beT+aUrZwJnS#QN?y_rJnbUwY({#j?E}=zRCb+S=Oa)H@t;GStThr1}FNP`}x3EE|{5g&H<`!R;8@flg9X&0}q}Gk1YlV z&0_?Y#g)jwu@+-aZS8{_W}Jqr=>lH^y_-b^w~)<@#)c}~Kp(JeTo~>pEYZcyi|#b4 z4=)JAd5yhY=}ul&4M)oN6fu9Hpd63jv%N<_xgNo1_LMJ|vfA>>ds_SVQ1jkNyl0a3 z>g7hY2sb$9(|lis->vRD_13(vzyHJC!yC8xVpz8n*N#b5W7LpVE)J~wu%}7I!JH3} z?;vcvQeoU6ux&vn3$Cduk(dxP6k}>IT_0-9_vTv+f*e9+a^*UqD-8LqQIoVfDt^3f z#ZepdxElShD;C$DF3v&q)2tamhF)c*DwOmNi&~_Bh53^zzlb7UHQ4bKa+Q0l%34rt zDV;ECvNP-C#TXAba7^qoA z?kyMu0qzoe2P_n*YB;t~!M0+BMiX#_o*toD?NmcvM;VTQejf}xnicy_#suxBDuai1 zpCsOchIJ8|X^ETXL-H;y34m}J7jpFew_A6ASx+HEpQ6hHr0pp}$5wD9v81mwk)lCy zAyzE9b9eV<0>NT~(by=7o*I%PyrDLoj&S#G@683SFSrmZ5R@CTCX4JP+}-JYC0xEt zM2P2RvOC7^5^ICvC(0RaTL7Au4F_b&OV92OXXlnsNVyTgk0YBt)WS9>aMxC92s6SZ z`5Im7wMQRx*Cgp$eDIx&DTuA4Np!*4t|eRvJw+{JjqJKCXKo8c0+^lQYMU~;(8p9Q zcc`XoK5#tg#%2%Y5Zhu`?L`S~O+U`Tb7pHc0U00~R(Hgx^8(*Bn}A$X2(3JYP0FsU zT57xiIbm!z^bOGF^HQ@C0OpaL5(;V<{{aBQz@qh)0G8uu-_2Vh>}R&7{zvE*^rhhe zej-9y;kTeAav)5_B5tH{NG3KMv-TvMQM_fKFAa1IY8TOAbp@Sd{}uCW2*sUieW?+*c7aq~a5j0~RIw;BP=ZXOdd zi}RG)u1xo|E)O{!Zw(0R3BkK@&mmr`Lm*xrdC2OD`x(#R(c9+xJ>|4D$X%3kc6zSb z@Z^$;PP}oSyd8$$NhjIaFINCRRs(Ez^WDym`KT$;_nv^QCw}h<@L3~#H)&8#Rrs~^ z__Yp6cyb6~b(34%Oc9g^QV!r3Ni0>l8M@5K+57lfXd@!H(E-Pa(r6L0*wOi`L)Ts7 zYSWNFdycZgxk_F-cW>+S*Jnj(Sj6GBm?OvC>i$thL3Zm720V8UaU@{H)7I?MT%R4n zvO=$N44xbI79Jnn6{RXB@dJ%ltYb>7u57YeMIsx$E|q50O;)XDjdNeTCL~si9e(IM z77j+JZ`JAM=oNQ5%5AB6((ED)U$L96c3eaP<|<1A0VC83{T@4B5T{oNT}MT&Rjzcj zehd{^W0B<=4kU*QA=q=IF8x}{RZR6ZsfAUniXh+LVxsDh)yLHNKf>89y zB-D#WLh-Q-Z`-8Ty*&D!B=Ig;w-NKl=3~2jod^l#M6dhIr22_DOl*6aNtK(eCD&!8 zR8;jLUVE6ko7vaEyQj404ajmdm7ocmEA)-dLnMeHZa2x9 zc#z;EVd&@bP%^fOS5b7>fm&~pO$v1}RzVE`UB0fNQVju7K6a@-Y`XikQ;o~ejQ1MV z0B~8%xSnZ^U(gf~@e!NPoHx@dTRvv5?o{_-?D4#Uv^ewd98sOOOjSRogAw$V1PSaW2Zl7Q`Wp4_A?5B&@Au0w z-N_M~1AVW+?rPAg7L2Q*JK!V5$7|I4@g>;12QA*Cou zFM|2w0%0S#t>~B0T6MmjXc26zRs6^Cxz3CbWjOU~!-+E;@|vU2RfR3luQkI9>rRY` z;Ps@S=T|2yUe+;k#p?Cx!2E+8sv-9OC?Yc|7Y*)XRpBU$TQLx`hLz1KE1Ts?Q!bY6 z%8((uL6(m#nwl7>Nur`Y)ZqqYfnM`UqRV8=YY@Chj2IPzRCu!GZpZzRV2>f}dLsE3 z{N+xPri)RT<^hsB{Ap(XP`n*F$90-(;0=w@Zqx;2+g%3u#{_b;7zIR=!Oj)N3F+*g zcg~`Q7Lnr==bEIsn8l;q6f)Am-4FYWup$N7a`fp1x6UWqscyjU?8XB-?1b^yVBGF1 zQ?~^3eNf-UN8Kc1_7O%~C!R;cbyaZ#IRHaIyuZC$BI^Kj-_pcBMd`JW_%$}gUYyeH zuzguRXKT#!kXW;@WQDj$lL`${aa`J=&5-y2_h5n7wR6ewLy?z7A)I~EIY~XXFRAmu z6th5(x~g}(k5L>6GU2df~1By^gb0c zPqo|+I~}4ZM+_k%xJ%!_E$)|rThaiXE?89y&uie$0=Mk6B(`4xxA)>YBl$(DdejZv z^#pf)JJdW6i33BULtG*MRJ%n^P?Q#8)H&}6N2_pD)xmt08K)W;V}dT9ynaLoAS`;o z%zYH{x5Lbo!JSdkY~I)dGIE$iE-iH-pOKb2Elb^^Lp&Z>bsdXG+=k_bBN^6?2hKr% z@koZvV`rFzIRqL`i1!`hd&29rr10Ju*Efhurj57j-s!{l-y7#u$njy`1b4p!{lt0N zaK5L$UL!Z|_w32Jjjnt#-?wIZetWUCRQ?oML22A5+MlOu%3}$@UpO7X;coxbpC|=; zPHF*p04}jeN%OGfnrb%OO`c5^qTn!%*(**(wi{!RxPHaW=t+w2cHOnE+gjh@zRtOA z6mKI+Y8TRYUd&n&TQX@I>8?fHN3PBqNwY1u=iHSJbW(m2$-C;*opZ6}mJ_{Q4Amq} zS8xA&ITt&tP@y@4$>(6tmzAmMT#t-)4zuP4sU&%7jvYLz#;OTK#VeS)_^1Pi+b$tR zS2f%LGuq~ml6I-V4n$V#tR0)-)M~R%;BI)z&VIJV-5H@3;v`VJT{%)tdH1$@r(zXt z0-fq9j`YwyGnNi~??TgciBWsBz0zN8vCdz4mYFz zs@kOz0`fy7cEmb1-(Nj3l4x?iU0pLCvhrnTl?uv64~Xg^m)!Y7G~6M1q>V`2^JLV8 zwdZz^v?Y>(x#}KaEApP6=SoVEl;$#4!;`!I`m!SF=y|YxWi!H7zO7ER<~KF*+O}I1 zX%COIM+#@2%))hLxhY;x-$p4hS&2hZQg>?7c=LK3?I5*z!X_!HEO*i~)%Lok7#4HpfzNn(}n(V)Tr8$h0RM*GQE(6Jqu9v=&C=5v%DaOI4`7nK<-g_ zI%KV)cl^=cdK~oToln_!#k$wph_-t!CR@SvR#;Z;r=MZ97#?cvj7l+=Cx8D(cb|UiA&>4~j0RE;kx;*MasJ@PggbJevztALv&A`A4d#VrP}D zyGZW+X`MAw)a;DfVp9 z0vA43JrGrqwpnXM=51^Wc3D!5Fw5dv4$aWkfW*@p%lUSmk4TBaZdD2Bw8pI459!CY zDXZ69H>_L>RU^AKHkSttKVWdZoKI>9R$b0>j1%(r=5BXhmT@51w!R#7s?@>r4^(D5MY)7gYvlMObC~e?>z+}88$Plf#*X@$;FxabmsH`e;c#qqC&eA?- z-k04U_Xpx$p?0^sbKOPBr^S92s65dnyEg_;z zCyYA9$<(ojO>_XaAgRcv%5UpFp#`F60=SqUqs+zlPz73o7Z?UCzJw%XuPd4rXA!zE zztcAYx$P3hv#%{8Sxy(GNXyO?=`x~2UH8R1E`-j*omNb&<2N+PWq6kR1%Rb0A`)#{ z!zGvP&6hSoOe$`WTkg;3GCV1m0`O;_WR7mY>{2C_L4_+6Eh&t1sBZHg8SfaOQeT4`((lMK@e_jeV$y@#){>f2j@Fuh;6>oBAC?fY~pe}ji9H$T`9ckcJapx$B5V9pj z@`gb}9k)%w7o&-TIb6_UbO~n|V%0itK?wJOk!YK*3mEg_2W{?{N?zMp)pZKYg3mnD zk%ENBxHN92v$CcXt2BJ!P^?my1AWm^c2w4)-e1l9J-Qs!Fr9r6lP8!C9rk%ibRlHN zMij;~hr%v*Iramrv`}~o4^c*gN0$0pXF+Pd&6QXk6=PdfFe{;YtYGd55_zsgH4&{# z&&WIwgdO=Y0g(YYd$#r69OkZ;K&NRkpdQ|!kcYP+#j_1*5uuxUlUCQ?WON(1w)~De z$I*s_z&NQXkT0xkT7n0lCTiwj=9!Q3eD8;+lKNQBlNvR)s+POQ60F3PxX}aj8bU#o2WuZwi-c6D}vUQ{{*y23fYRv9i+Xkiza4)G-Ue+F|X5I+>~XU=;<7 z+1qhf=!%E!%jTZhV)Q6)$6JuDPUuDLUP0H$jTCbaGg<(LT;22^`HOOKTov=&&97ym zKijt~!=qJOw7HEW)?K!psG!C2G0kS#S*f;6vSyb^JOAhbTsrG*>8y{m%;Y?yRX25~ zH+we3hoc=~9AWJ6SCm7$xgY*7M3|}2tKN$f7OVU3V`@kn%kkCjKS?BP5EL1EYr}n^ zP9wPPRFEG4;LW{#wFQvROydNj`W58fiaa#a^&15jTvkIodTr^{4p(#3%uHBE7|mz+ zqApBh#sR*O)3qZdKXhRkB=2st8wc%T$sl7F8It9 z=6W)z?e2_8cTrcnt+ScO zkwOaIC^Qsymc5DhbYRD(CM7Nq+?vOGyXCzCw5N@>lI|Yi^^Fs`LsQgdu{HBq(k-5( zRt$3KC3;_b{D(QtdyKqA&^s!f%Pb0`sxL&PRIEg=Px1U3@m8{JYC?c4<2==(U>XaM z$th9-CH>l@ql&D|R3C0@-iPa5M58QC8CRAFQT?}#6e%;INwjrQrv`doc&({{*M$lU zYka8vFnO+7(+V-Hust>~jx}i_xF)Pb=$*LyLR>YdJHCi&A=*N6byeBaN(KwlwF=_; zizZxq%>()^5F19LRy3~^$lHKde6P~)L%pYPA=@R!yach8ZFfUvjcK5P7@ZM7Mi)|Y zvhqOXt0mE-Amontz}(WrYB7SMtSKXVcRZ>;PLnz`WMY1=Wjtk0+mWt)pEg@Pl0x&a zj7z=3GrveLl3yV~s57s+Xv(kAq;8k4+_vRLHyG{`scJ@77}?ap%AFKaL7Q!)t0>E} zW_pEJ@b?;Xd23qQvA9d5(ve(+V16&j1NE&hNKR|ID90=-X^RGAP07YcCX{c54Khme zcv;DOYR9--TXzi@?~wOeCL7O*NM!+L5DCJzL0MHt@|;5mbnb&~_w}}d>Fj19RCBsv zZHqg#yLJrI$h=WpeL@+!*yU>nPt?j18b0|P?mtz+jHw21u zt7e#*?OlinmMz>6so$dpDg$q&Sn=!>E;e+Qw{pi(fL`-ay8Snl=gI*CQp8P-$$p0% z9+XkDT!w-;uT_o>?!z-W9InTH0@r98OFG6I3#eIKPbtMJc?uS2Jz{PVtEHWWkMFEimz}D!wbP5rkQE?$Cr|09$w7MP zp>c$f_8V@*Ktjx0Y0Z7j8nGK(vPA?1X#P9B0?Y(p)wP>~SynxI!)J4HQx21VsC zf&|I0OO#MkIw+1Ns0Ef>EFjL?%v1Yulw3x?V+l!b77e zi0X`x9rIegH1wonHOIdo;OhA*pE5uT$il|Zba_&TBo0>9sjUIUIsYiY5pOsJIDhqj z9l(YEv4{(RTzEV!z%^7YW>C+N`ADcXLV*RB;N!)&x&KFGw;GFz+#n5iF?9yO4I*e| z_>P(DtW`8-^BnAmCy*lpa;jiOqeiu|F{cda=#IENS1PnP;`)+~D2I^JaPoat4I?`2 zD4;e)@r~|+j)PP{%_gakW~r)0O+uFj!agcz%#Fh%&^cYEo$aHoU1l z-N$4-5Fn!<;GQH%eC;9r4Oe;Xsy=8qGiFP8I#ShzUX%_UYId`utmn-t)JzxhWg8b4 zWQn(hi9oVYzh#AM(CgI+BM>gWIGosZe<$2Il2q9 zq~ePC3CT^BI}K;hKtpxX0z!jk1JCHS|8OGC^R8jB?)(R-Oa0u=c+H; z04h~(&KVDvML)L0#HD}ip3cbIcTx-a7DBmKmyV^cs{N{=leiYk%kQy_%v$?MHNqQ- z7$e4|KrBJyZ!#2)-fCh_#-mf;?pd~=`0A}i$$j*8YV|kO-m2oFq3Ow3{qzIwwIapoh2xf*}?83@t-0ZmlV zp-Glx${*?ALu`|hC*WD=AksUl!@V+<+uRcL#A)|l>`*Jd2{zA@%L13)(!~y`tvpL7O!>pAHlM=|x#kpXw|Xzr zkzFJfmEpiqwJ&SVCAuje9Dxd5${2p}ILqHtQP%5t3qRrQqwyvy5tKz&{!V@b8|h_s zJi-b8gFwf7d~}S2BgHPpZ%`R2FU?7=lKK4?MHnFfL-IDY?pgtMC}cmT1h#paf;5InR2bb9eYCATB0uKyo(#$ z>v8ui$?ftoHC?^fOy5gAUR~0=&_Q#S;$zXu;@M`k^k_ymQDf6^;~;?(Oh2AW`5?66QnzhN7~lIxHR~>~32P!76z$F$Y@ejS zRVAdC>KePo{?)F=oZ_{+>5bCGC#|U8DOEB20tmPkI6{qZJ!-uO{1~1WjnDGM!ozQ=-%Kn7x_JDJ~ZtZ5c&q=9_a@AjO25dpenO z#~@^~$=QHiuDb-RiaYJ?==Sgo=)>DDPd-3!S3Pc`N2loTDwBANSYo3{de+>gth_jX zX5aUR-QaFq)yppNTJdDRM_rp+3yF#<-iInrkn5ves+DW7)Clq#c1ff=;9mpz2Rqd6 zR3o`c+9s9jm1GEiMAoDMbu~yA7*@v*Dd$r_PRqU>hHZD|N1M}sQLK_AH(dl!ezHbd z6q)=i@Fj^Lm|dkWZm`ndB}M&RD2`y@kt>KU6(elhvutODG*7} z%$T4y&B3TWla6gWK#|X2PQmFSjP4!ML*fZJBxUs^ZKY&Z$~^H^eVM&_ONqh3Vuic-Lf(48T97DZIiR+f@*lSN_g=6YkAA1m-o*+< zNpmG@?KpcciAw~h!>-p(d0>2BXXi2fRxpyJk^B#)>_Y*({C*zuh&BK=eAHRhUrhoV zt|i_Lt`K{@f_k{1>BC;HSfy-+E<{~V^>u5)>UD9m&s(a(sk;85M#r^Wiu$;t>lw}| zf?4fJ%7=J6b8S$7_OPgRN( zc9-o?s$xaMJj&So!FHZ16s2tAvt6Yo47mh?eAdlSvW@ScpfP zmmq$Kj%w8Lb2zhSOL`$F|JMECq$2F#!|LRq_*qh(O|si{S4`@=(;hao_LuiasBjld z5r0ojSGzf`FEA_Is)i}au|=WX&GEje+9qTkC%AMW&Ls|^a zsUj=*z;?|sSGQnAXURh`J+E9Y+tX_hed?eFHrCxvg?$H^LR4u**iXRH97tQ*rELu0 zs!Mg%@HE|=H_XA1=tBW3T_-Lk`<<~!7S;AJI}XUifiP6sK;_u1Ra@qhdZ0wvn4ZZl zE0bN)LYX+T66B3)vn-o-est$pBgdlz58@$o666)v^9gJ2%u{q-V9uRMqOJRyQZg;? z$u!;0=W^a%XQ9!>Np59bshN7K<+6@1ibrA_2s%DzdDmUooNxDtZ+Yx{F;}gDgskUl z2AFWIu-y_8Za!r(Z5)yj!lH<*cc?$dE5XmyEmL+rA2@?d>p zN1c-82kRbQrwZH#xF||egg(AVDWN5>q#~e*sTdix3OIvH6fA{B%ZtH5*6IS;GG()z zrYwnv*A6c;95VR+_?luWsh}X!=A}^?wdU@rf>zX+_#uy_!s|nw+zYIyHQ5Jr`{gFC zBdyMDHidHEQvGV1I;dxGLVAI+bebKfv^vUmrIS6qwTllZlE#;bNag|C! zuh6RiGtUFvTcC=!p&A;0--F^&6J1-b4Bgo!Qgrwz*O8Zj5fXBD?LDXJRn3Vor%DYX zQ3eQWiq}507)_R1jvDZGS(Yk6sc~c!r<-f4U;%;dKp3(={INMIubWh+D^2e9a1d4R z5)i^PuC3_Q`f!=mFbO%lhx>;L0%ghMJhJx{O)k|UYY6Zywz*_{3}a>g@z{qRsnL`* z5QPp?Hrdv>c7e{1!@g8%9Qq{tX%;Waj+&p2Vy$gKWrExU+9NW*%?(+Y%aDa-$i~80 zo{q86mmza~4mV_Z_X_7})8c8!+7hHsF54BoTB{H;$<0h|$V5Wvr_s5NN8BT(ao)cV zZrIJ=rR75o7vFgsE#X<+P%dB-YeDAcD&4b#ZLdr1G9*I%jisb-9&!lXa7?WQT|x}-7RrKwgkyyL&B z6D_H|BdA!fDy!FXQjV%1Mm1a^*c0Ad7oa!_3~vpc&(Ar&4{uTwbfM4yfsHIPFySu% zZ%jR7K1S2Dwg^iWrs`x!po;cD9qz7CjJ83bTuWu641+Q;u^d&8g{p9NIwlX$YV*g^ zfjDVh#v`9v-6ia9wUSb+Yk(Scv%;ZvTSQIo2jSHVo_HC`J#ds51;-6)P6gb`5f{;# zy|^j}s%)!{LRu1ep!zI+A0Jj3_R`6TGEE@MI}32b`6ZisFPzGGl{r*WhsYXRZ$~+^ zCyKByeg7A`MlTRiD!j1f6QOvr5K<$5Vr$$<6#0Q2B> zCxHeTHF%0yS(eyxLpcI3Fw=!{v%9F1!&CXG6%fl1+OrouV<`teTX}a#Tp)Nbs95rBX0MVW4O|E2D%SiB4rext-9GgEYmtxK4R~l z+VT!A6Ri}#^dprWlpDzOZ{C0wRy~5L;H&E>H}5?G+T3vVTG?Uy#(N^Qs4g-Nw29;A zWrys#hVKEJPi5e|6oVSNOr0r+gA1t$$it$e*xTWe-aE;lc0}wRETr+5<`y(9dlqD&E(MNf3DqO*2PlH=S zHe@Ip(nQ%~pw`n0)RwSKS-ifHmJ_*Zf7!WcLTJZ) zyd5l#gzw1IpWLIR_x+V%*M7KzV zZO@H#l{jq$d!Aoq&eEEfHux(v6T2}Z3m-&P{@@)}icz!p0EoPK21~ellFAHXI4{zK z-Q22JPerNfpsU&+y7q~V(L@|Zp5=?VjRP|wN=exQwSbv}$`(S6I=Xg8szuUufrCJS zM>BEr8B?S>t7LWAvCx+vS;87>zkJN2X|WIYBUfuwaRuah_>ncNhl`VD!y(dGlm?lH z9~sV<;jQ2|Nui5q%Uj6SHzRTx@H5U;z7G~tg)eV;$gMJ`)XN#NU4h0Etvb+()^r<2A8Op<`7iurhLw7NhBzXn?6t9&blghk25LaF%yLsTF>Ti za{Ta1C$f5*O;eK@y4mcRc9iYXel|ljibmPalC=;ra*-$!A`ys@WBS}!I9N!&tHnJv zOo$EIk2EW3yKaSDtMNr&aOBu{CPBjztZ|!X=dWJu+ETl4yXI)wl99=A^yGQTHEMpP zNtQl_O6`$zuDKP9`%;dzIOcBF`f0w}HeAF#NL~t_>VB6J4KH0uv2BGkUGDyydV{4e zHBP)X4&52d^k)eIxX5p%krfKrdbiLZtOIpO%HxkAOwltVe`sLHkTS4~N zzJd>yr8O%$gp{=iD_`g>y*#rcMi+XO!H?6pCV)2cDzv4QBPpwxJ(k8yJ&M_30j(c7lF33_ICCWVOo7xD&?U)v zuinS?_>!s-qr?s_2WKSBmQ{*G`oG?T)whW0LKxYqZGBtcVN-XCXExHc4OMn;?@#`F zc0|dK?e+l65Lv%Br#juGP|9hFNpuGPkyd1H}7SD3jjmsQV zK8vSWTbfVHtc@dIvUEnrnai~2aOMIt0?Xg?N}KzNd%jOAdh6=lS&v)3&%H-cR($>p zPe{7*vL4IJ^ZHuZ*^lk~U=PM*3R^x3x@9b5Kr?f zB2p4^%Zeyh=o8ua4nF;dm+1`CTsk|wv|$36#-NwGvv|+&@6tNF^qFfLL0RG%>P3@m zhEmn8@AZ=n``5E!<%HW2p}mf|{O<>Ib;H~ylyNVuy~@ub5WK)9Rd<;O(n2koRA~n-Luy3SEISyff);KFM1jHthD7ZN=*GN*=u*s{q9XV^9TqHFkDRr`$HeGs7WOLm7 zatqs>cqju$Z)==E+O7^I@!im0nl|;i#7#T*B1F1#vv!8oE(0oWa<{beG;r!)Cluv4 z#aBuFI!Mqi2A>AnRTgXFAh#x!3i~4J zV(UI_bR?DGFe-u9{yw%zl-t@A`_=~&`{=^neJ**~E7zXX0E1Vo$0batc@WbZW4x?3 zUH+3fqZ(awUUR{7a+ZS&6At&nrv*BAz$yunnM~0loR03)7^o)#Z#=1YoD!g3 zzwWlrdS7{6Grvn*edWXZ%9%G@RqF8U$+U-uTJ737zBHfGK|IcfP4rH5S~GUAKUM`y za~WJNe%!ml?iJ*{voo9+@#N%I2KVOEd}3i^H}`wC2eN3|&-lr1v+U32??JlhI z+E;GJ*Yq0QmP<(@A~{$w6GGNd`Vwr3#L|fz$npR;2tzs^T^3@Jf|7LH**QMFb3Cg< zyiFH5XxiXa<*rptId;Ve4Khhuwo+7tJkroI?(Mdm$>B}v7}tyl>;in>t}#n>G;Ry% zwNXf}r<@ z-q9s97Z=6v&fEaa`acu&ml2`{V{S)Ooxn8&D;g8;jZTQhc>}L#)b*(rwg_v|w2e?Q z26`(YZAQPxCgN;BpHpQ4A!PJ!hz3{3Lh&&?0=#B6l9#0~m6A7-*Wf6==s!4$hx?{J zM=3hjBp$t~>CnpdVT7liL~<0Xe`Zr|xsDLsvHaTQt|qk=7n~xtp;!v2TZe&sS8qg3 zS91y*HxGOF&z&RsiVVaBz=~_Pu7l0gR~j3)vu$N<1mJ=(z$n2a$1ny;pnJ~TL>JrE zB*iKf;Ox9i^PgML6k$p+$Rq2McC>l{qD~EMv1>^|Ov$ww?A84pBn4sDn))50v)(23 za9z+~dj{_|52Cw7^$fK*MUOG0YHE0G<(MEE=JZZ#SME(J$=GWp3{mC2Dte97hr@+eZpt)-{1Dz^5V0l_fW423)cer3&D-7M$o7 zK?Qs`H#(^ofdc>u({Tiy@3F{g=kNw~P5+ih^@>IIS=vkCv0rfS~DcHm+Tb|!Na*#(Tg zpy1ysdz9!}O=c%bey&K{5+(j6b2Dz6;wpz3z+G{Vw*2DWzL!ItX z?|!+la`vH!zOuIJiAo6}ZJSMyJN|1=${&Ls>(6`zL!$CjC0N;>$k!$Xel z$Wmz>lQ6@4!XS7{vT2SYcqgGO6n7+xx*dWT+P09>ye(+JTv;k?11hb#Y7A5s1nbpp z5l(v?Lw9C%P_Z^OccF%)nS7|APxMwUTon8yj(62qyi_IrsDQcnUfz8uX>3jo5Ofaw z28p|%LC}6rZ5Zh-)&BmLX&DDFi!2|2aQR;y!j-3^E6eH--u22pEMK^(Fe7DRW7GAvniXF&XmQ6~}^spMOwCqb%20i`+ z;eT@or~0-XN#sxU>ynx&W3{hY2Sy3mCE#oRbfW)jdwG|}hJE>xj^0;0e#v8TsD0RT zBwLXRyc0HO2DqsW-Cf;Q@M@b6%cFMgL3z%4teuj}t__ov2ySQnW%<%U-g3O=y;f(G zzb+r~bU*T5lTE%p%T$^Ys?~T*?6W@dw%}iji3Ne3YNwc3mRHR3)4n?oXq#@ji^VL7 z2`}B$T_W-BG{5Fp@2>yeQnlsm`bG;9i$MY?>P$R!H6o|paPBw$m%&QF>_zmW_a-3% z{vPbWjC_X+BDTm@;r#H?rgf;-(40yM7S5s`!r;y?99BiM)1jy=vUlhss;wT?YO$u4 zZi|;ut~ZN#^qZ=BATvkMoY<69kKL%n#nHs~()=b>Z4PV})pELKPD%QJk{Mz3Ef)MN z0dEa8^{QcrM~(`EQ!Gq=M#6mg!p>LtwEF(#sQ~-34f$VO+JRAunX93+tUVH)Vwm9S z92GlmQ|YJx;c4hQtBhoZHD{c$398>L_WhY$x@M=^nm3SEatexs#qa9NnAXE8wrj{* zjs$W)dh6#N>AIdlyxmXZ)QMe~S|@MK?hevtqepD-Q_0 zD^S=wHxcLRf{?|evyh5T#wkjXWc8}_i)*V&BzNU|D_Rd_bftngT;*XDN%xPYuy}|7 zTSs@rQPqpCi<`z-#(Ro5iyc75S&3&LJa|~FI}(Slu6$Kn>SxUdrAw`$xVrh$Ef$wb zfmWR^2@`cZ`Ipb`IMmG>pXu$YKDpqqsY|tHagr)(+lFfP-7QE*Wve{CO)2?CWN%h8 zFkRWvmit?FVbRh*8I;NxmbtHao*fydMTb@aR2TwOcEUx8W0SW8SVh=V7Z(DIIV zKlh%idd;+Mg(5Ezl@lMX6v*vE($xF>wxj!T*DUrg@R}0oYgTo)q??wN%SJIFCG{Fb zv8suTlN_-cju$xs>~;RDkKxlT@oJnijA=w9r6dgOX%|X`wiz`uV74LjI9kH+#_y>r zZcJ)^j7^%Nad`4cauCQ0S1r>xZ_8D7Q%#rAlTR8s-inBW4Z{ykvuZw|a|*f$i)y2; z8mp3(E^9*-hf!akFyi&r)|Eh*u}pU@m}RJ`nmn1Oa~ zyKw2f$g;KO#z+(ap2MS(Z(At_o5$O%C?Y`_M_uClTlt10b-pyR8i`P1B|EUe8pab?g7wZFbloxe1$=`|RSSTPStJt6jav%KPP(WV(#s z+sL`~oK50wUEE9ciEDsy-+-AoaDwHs?7p9VtsP@>}it3p~s?D4P%t#y4U}#H%4hTtt(%3s#5pAV_(iM0|3W&N=4N&Tb zNh~sO)7?)#(3R(np0K4(sXxlW^NX|FIOq=3mPJC|t)bl1%Q&n3_Pq|Xvm*AwYCjjh zqjIyA`LvR9SzXz*Zy_Bc&$?O2zC-_?M~%4Wmv8pdlS8K)>$;~h2pU$K(+1K+hgPQ3 zHeE_lb#<6Y%r7xx4oy9I=MlA76jXWqp??$zsWKWNY9FP@$PT4ni?b9o%eaQ(BXFvc zOQA{I%vNCZC^=Lk*D{E@UXwJEk@)d{<^_-UByQSUGDHAE8R{!fd zEkhj=)vK-A<|Kxgykp`RmBT|-d~%j&ZKAuj(XD4*GW&Lb&BDb*L8IM@{`J$ZfBuR1 z^;LrgY(p~>T`GEF?b6NzVop)1k#|#Jw+v91d4QtXNf)>SlWQ2Hv&P{-@jOGA`7sPo z_EMDGU;8gTYd_yc9WQq7ZZ+;1kvmf%Zy=IA+Q_mSjcI7w$EMYo%eAoiSoAzq7*{jO z$Af9FY~G+Vwo;1^8y9We259wOw#Gz(MS7XAx6u2jAWU>Tc zB~G`@s|u9DqqKF#me_QS#=4f3)>QYd&wCfyYjNUh2^{S$N*3QTukGY8so%v#QpIGY z&*k0Hni)k4Sj8@i=X}ZfV8Ra(1-}QEcu#OuGg@siCsyf5sTb%pAl7sLvK8mudtS7! zSZ}&ks%LhFG{4)JiXo@ys&?I7rC-I&%j(8bNm7~%no`p-FoqjS)y9$P^o?YWv#6Pa z$IdzW?08o|4K7S1>~iC-Mao4ZD~B~zKG9zi?{}H*BDSGzD)L#Y%RVf+uN%XrkP^EG zeojuTFaHqaBid%{l@=2K^7bR+B*QKk&2@ErN&>2;8;TT8| zn$gvEWL>BCSyLl~lVf90ws!f#13LgGNwVWQdM`F#*fHz~0TMT3l^1qly0Oa)-0Pp4 znI0SWcjxYJeF9D*%x5ZCmN4H2(yR((-{|$P7h3L@J1=VAiMAYV-85OH@}*$1$}v+k z?BUZEiLwY<$XB2y zWsHm&ciaRs1T<%Ihzhj``Z^nW4+nm z>6+#@``Zc7a@b#@c6`wOg7+X#=qHHY?C(p{n@`g_2vkf)Y*GVfaYw{4$kQKSHSD4@ zo<|uVEv`Kfi@-XQ;#bib`d|Z>tr4}W4piL*&gf0!((?+{o5-h{#xF~kAeEk2s$6OV z2^B72we<{;32{A0Yd8(>be6MhWs)7mQ0K>c46jliqBBSvz547i@d!`DD_yneL_{~d z>1ue@eH}21!G4Ym@AyHt8eHJh@FLxcVEV5=RhuBmj!dbcPqx9}whAto2DUOB%$Ur` zdR#$sU7{=1F}$dX0E$2$23TdU*O?!kDYQ-I0~lD)$5TofBi9SD3tT79DMwKs%4*a} zix~VQeLfo8r{eR~11GLGUSgNiAS_ZKe4wjD`Iw7-YCfku1Wk&^W!T=wz7usVm<-eD zQKN~FRd1a5iyXsVTG^4+H+*fEJ7%{Uta=PIfN-nxzzCcZmE&Hp^BP1o$tZ42Xyl$PzBN$YB&SnZUNN%zhL0w%NdH`f-0S z39M1d(c3bk$h|1IFnOE;Ypylrpk%PDF|^2~NnVql*T#KiBUzj{t}nH!Q&x7bmJn3SWl->#i^c9*$vfvBPO_?P{5jJT_lgONiMNmy`l zPU4W%pqY(|99E)|$j^QN@NfR{%dftA^SAGQ`_1p)z4`L%|N6%-zx}5-zx(R*@4kEU z`FC%=`MWoN_4(It68h@jeErS$Z@$G7fB)|5FW!B@H~jkDKfd|m-FIL9?(3huaN41Q zJ$tv!1mY!lfILF09 z^ra%hEh6UUr5#_9hl?EJF1AA1(r|4zsmgzge6c06w|+fuXcf_8pZ$8Y`t0HDu_4Ty z&$k;lIA55Qiu$)^JSkSFpO#j}g{s@x$K!Ox)w;cZ__f$u%!;v!gtT8%YM7vBGMt@x za}63AhaQ^3WDC;d8YfX=YLW}xDnUA$dhkXU!D0*%EXFK){&e-uslj-VO)Iw*30*Xe zb(mT7soBb@?j|abxL>v)X{DpK%;LFi&xpxNHzeTr92HDPk)!TgKu~UciNi0g7>BZ0 z4q$QGA?Ycaf0}+&&>Sg6V0)kNr&jT*DOnMNIFpoEK$LwHdtrnYVa}P|eIODfFd%vO zlWnT5Ljdc+WCOumT`-g4JzsKtHQ`4NMuP70?J+BsEeH7}ymW(TZ{+Y*G2P`NYf9V0dbs%&t4@`P$)y)6H)E1#piU5=x%W% z;3fmxs0l<+V-)5E1|sQ6q;~nod4s`g zuIbjVb5T~xr(IK5D-WOu4YEzH^|0HcxDJZFA!|4s!bOc4`UR-rAORCu$&nLgyJ$>) zAbD1#kutndY;>QP+=I`rQ}H(Brh(&Y0J#@mV_SpBHBHD_x!W@FxpqN{Xhl&CFT#%Vke&!9PY!2Q zR46ZDjI@~pQZ7iDDk1?kvdEUaSy)PfVY!Ium|1>k;x5g=>_fQ98jBoU(QI{`W&!IE z6z?zEir(X7!1Z`vqhNjDa-}~*BT`(LfbUqhx~gk;?yRGocqbyyCdRN+B~;o}CwM$- zukp~Ty%;71<-(ujH|JgzeyvA+@5kA)Ky`Z`|ChQI@3h->XXn#euE;|j>{B8mPj7x5$2=>@voO!XHmGE9ze?jmgvQyb>?#@mpp^R zQ42mW=m+1N?3S+=vOzpCE&YQG|V*nS> z<=1Qn8;H`BxUL0bWbv%f!AyZV8=@^2A1tT%VDR+L4Bpu?7lHn0RSQQmnzwd7oQ!+2 zblcIYwq`P)WRlPmi)FYmMdaB}pXDBW2Ids3%tppj7x2apvJwnGT* zSWL~8thbo#g3}LQ-YLCzT&L~pZeLL(i&Bzq&i;h+PhKJ#d2N4 z;0C)W+LZkU=F~Wz7dJS%mXddpr`EQfw4rDW(?o^DxUj+#6DN zc^37(yks=$v&!Kn9~gHUl=p3KR{sao7w>`&A0}OHk`vyB$|NzYcO0qd-8FVM^z9cBWm!ZRbVi#MGHC0~1YRl@E-_WM1IMANC1ADV0K z_23#!x298FpG~JrOzpL%)4GYHHScMe-G}>0<+(Dukwa59$C5hxUjq*|Ty zB=>vxah0R$6XjNwHF}jb<)cc&W>YxT{o_eb3BX2`H6#Mtn;C*Ak1=r=BR#54%2O8A zu>|P6GPg{9V2n%xuq~3znv?+Nzn~XgF@SAh-HZ@UL#Pri9#9zGOPO8K1ip;vN8hSh?lHifj89ztNa?38vB>FM@(YM=d`UC69`?6uP7- z?68H34>hyU7)NAk8CBVyw~)b3Jl1ra?J=yMPTf-%)98J7&i9U zOE&65_hP2VQ5vsYB?W($Y?z@M+(dfww5d0-)p_*tENN4k_l`cOQ2V6B>5rs7b_qt2jRKR&$%-OGIwZLJC zaNKMKP1INxFRFWzw{S0DRXraXsm+7yx*F40vrqK$N^(oxHFuuur^(!J_T7M;K2tP{ z4B7=wa+2Z%+kYrGFUHMdu^{?o7uap*Yc_RnKHX`N-WeB?D?7Xgf@(~pQIa&V z8cQYw1Ldp$xet+$-9EROs*{h^plVwq1L3d|TtWRPF^W_Dal=a~C5F|rOcnZ^wiGlufg zj;M1%6cY{!(X=3%kNj6~hz=(~iF%xlhP-~L_79d-qQsk67Cp+Kq*jQz>@xz$oGMD+ zL<#&$G0PPa(`K$-)!5(gdYL$Q%BeUaCBCYJT)UQJz?~o{h4r?7)`#=ZdmOzpU zRnJCTVobfhP#Sg}L#9bq9keWoh!w&5%=tA~L-o-^_pr%(UTNxjh#yv2#dZv`r;s8t zU0kADf_Tvq?`H7c4tsBk-MYZNw`S9B*^H#fiMe0RV-mk_%Vqy)^LTf8@7gV!krcUz zpBT5!_~6a1p3cGmHn3)~szkCL)ZS`8sd(p*sz;nOv!EF*Qm(pDJpDna zmLVyYxkLzw#i1lYUY%p=@$*s?XJxFy35dFtaOatoYMXpa8rM>himSup$z440u3Re|{peXj(%xLxs#9|jMim=y_gi;>*00nXU zir%1?_VFmK+Jh8W+s`ea?<7vFviIhnkO)PaQ&Z<+@CJL#}WXQ$P?un zXe3tY$7%mpo9dPH5nY_{ltL?L+%T{60R@A~@t0=S4$Qz_s_?Y^i*)@nbnIL_8x-XQ7HI~K{OYa zje!4*7kFx91*-3}zW0`tR}xABuO!f_1L?bj!4#UiKtwxVN6^250M1AP@VdmNar0)i~cp|^_puTg-oQ~O7 zqb&}2)LAh1RXN-gbTq08F{7)tAbrFr+mOPNAwB2Wm6cO_0#`upHy7UykMm5qjk;?>Quz* zZ~1jwNAnsNWXVWmDRD9E_}dcq#(A^Zpka%25p%MG3NHkjOv9E{XS-nY2Gbf>j;VB? zq#lziLmceGCSn8oNT+bvJ-OSDg3MHgSCTHpi#p>-OOb!1N8~X3yP+UB)`%*p!uGlG zfGf~_tv-0UNE8pBr!aYHaN(0Tsp-l5>Vq0?s`++hpEWL_8H72Ljz+!8GZUnET_rgvVU9-Gn*JLnbg7~C2GH2o2)V`XHnJflgQyf&XN8+ 4*np.pi: + m1['angle'] = 315 + 1.5*np.sin(phase) + m1a['angle'] = 315 + 1.5*np.sin(phase) + else: + m2['angle'] = 135 + 1.5*np.sin(phase) + m2a['angle'] = 135 + 1.5*np.sin(phase) + phase += 0.2 + +timer = QtCore.QTimer() +timer.timeout.connect(update) +timer.start(40) + + + + + +## Start Qt event loop unless running in interactive mode or using pyside. +if __name__ == '__main__': + import sys + if (sys.flags.interactive != 1) or not hasattr(QtCore, 'PYQT_VERSION'): + QtGui.QApplication.instance().exec_() diff --git a/examples/relativity b/examples/relativity new file mode 160000 index 00000000..876a0a80 --- /dev/null +++ b/examples/relativity @@ -0,0 +1 @@ +Subproject commit 876a0a80b705dad71e5b1addab9b859cfc292f20 diff --git a/examples/relativity_demo.py b/examples/relativity_demo.py new file mode 100644 index 00000000..24a1f476 --- /dev/null +++ b/examples/relativity_demo.py @@ -0,0 +1,23 @@ +# -*- coding: utf-8 -*- +""" +Special relativity simulation + + + +""" +import initExample ## Add path to library (just for examples; you do not need this) +import pyqtgraph as pg +from relativity import RelativityGUI + +pg.mkQApp() +win = RelativityGUI() +win.setWindowTitle("Relativity!") +win.resize(1100,700) +win.show() +win.loadPreset(None, 'Twin Paradox (grid)') + +## Start Qt event loop unless running in interactive mode or using pyside. +if __name__ == '__main__': + import sys + if (sys.flags.interactive != 1) or not hasattr(pg.QtCore, 'PYQT_VERSION'): + pg.QtGui.QApplication.instance().exec_() diff --git a/examples/verlet_chain/__init__.py b/examples/verlet_chain/__init__.py new file mode 100644 index 00000000..abd9e103 --- /dev/null +++ b/examples/verlet_chain/__init__.py @@ -0,0 +1 @@ +from chain import ChainSim \ No newline at end of file diff --git a/examples/verlet_chain/chain.py b/examples/verlet_chain/chain.py new file mode 100644 index 00000000..9671e0a5 --- /dev/null +++ b/examples/verlet_chain/chain.py @@ -0,0 +1,110 @@ +import pyqtgraph as pg +import numpy as np +import time +from relax import relax + + +class ChainSim(pg.QtCore.QObject): + + stepped = pg.QtCore.Signal() + relaxed = pg.QtCore.Signal() + + def __init__(self): + pg.QtCore.QObject.__init__(self) + + self.damping = 0.1 # 0=full damping, 1=no damping + self.relaxPerStep = 10 + self.maxTimeStep = 0.01 + + self.pos = None # (Npts, 2) float + self.mass = None # (Npts) float + self.fixed = None # (Npts) bool + self.links = None # (Nlinks, 2), uint + self.lengths = None # (Nlinks), float + self.push = None # (Nlinks), bool + self.pull = None # (Nlinks), bool + + self.initialized = False + self.lasttime = None + self.lastpos = None + + def init(self): + if self.initialized: + return + + assert None not in [self.pos, self.mass, self.links, self.lengths] + + if self.fixed is None: + self.fixed = np.zeros(self.pos.shape[0], dtype=bool) + if self.push is None: + self.push = np.ones(self.links.shape[0], dtype=bool) + if self.pull is None: + self.pull = np.ones(self.links.shape[0], dtype=bool) + + + # precompute relative masses across links + l1 = self.links[:,0] + l2 = self.links[:,1] + m1 = self.mass[l1] + m2 = self.mass[l2] + self.mrel1 = (m1 / (m1+m2))[:,np.newaxis] + self.mrel1[self.fixed[l1]] = 1 # fixed point constraint + self.mrel1[self.fixed[l2]] = 0 + self.mrel2 = 1.0 - self.mrel1 + + for i in range(100): + self.relax(n=10) + + self.initialized = True + + def makeGraph(self): + #g1 = pg.GraphItem(pos=self.pos, adj=self.links[self.rope], pen=0.2, symbol=None) + brushes = np.where(self.fixed, pg.mkBrush(0,0,0,255), pg.mkBrush(50,50,200,255)) + g2 = pg.GraphItem(pos=self.pos, adj=self.links[self.push & self.pull], pen=0.5, brush=brushes, symbol='o', size=(self.mass**0.33), pxMode=False) + p = pg.ItemGroup() + #p.addItem(g1) + p.addItem(g2) + return p + + def update(self): + # approximate physics with verlet integration + + now = pg.ptime.time() + if self.lasttime is None: + dt = 0 + else: + dt = now - self.lasttime + self.lasttime = now + + if self.lastpos is None: + self.lastpos = self.pos + + # remember fixed positions + fixedpos = self.pos[self.fixed] + + while dt > 0: + dt1 = min(self.maxTimeStep, dt) + dt -= dt1 + + # compute motion since last timestep + dx = self.pos - self.lastpos + self.lastpos = self.pos + + # update positions for gravity and inertia + acc = np.array([[0, -5]]) * dt1 + inertia = dx * (self.damping**(dt1/self.mass))[:,np.newaxis] # with mass-dependent damping + self.pos = self.pos + inertia + acc + + self.pos[self.fixed] = fixedpos # fixed point constraint + + # correct for link constraints + self.relax(self.relaxPerStep) + self.stepped.emit() + + + def relax(self, n=50): + # speed up with C magic + relax(self.pos, self.links, self.mrel1, self.mrel2, self.lengths, self.push, self.pull, n) + self.relaxed.emit() + + diff --git a/examples/verlet_chain/make b/examples/verlet_chain/make new file mode 100755 index 00000000..78503f2a --- /dev/null +++ b/examples/verlet_chain/make @@ -0,0 +1,3 @@ +gcc -fPIC -c relax.c +gcc -shared -o maths.so relax.o + diff --git a/examples/verlet_chain/maths.so b/examples/verlet_chain/maths.so new file mode 100755 index 0000000000000000000000000000000000000000..62aff3214ec02e8d87bef57c290d33d7efc7f472 GIT binary patch literal 8017 zcmeHMeN0=|6~D$1AP{VlCS##fyeL&kH6D!RgC!${z=M~_BpJn8GBtZOwgFGZCiZg* z-4_{Y6vw5aRW;?GN+{B%sGX)ri?+0lHj`{^R;?24)Co=fgLSF~QKVA#A-a!v=iGab zdGFb7y8W@ga>4K1^E)5+-23k5yWdm2-6akOqvT`<7;-aZ0%?~5tyX4$w6j)L4$pd4 z$91LZnu00!-g?0hWz53?EMpz!YB&qjBQlaUmk731QnEu9?dqgmozy3qkyRmDA>6Q1 zp!mBb<#xJ5>JddtUx4axKP)lHFIqj@M7h??ouiK3QI|c45>WlFI7v zx;+4eIN{fG#K(Gny7Sb8x2o#EhtK`+)nDv=INY`HCdPnrEQ{LzLT0;zm9|$RhE=SF z-$C`=JFore`ESpkI{x4*Qyg*TW7#_@f5Ja@(bbqBKWre zFBCv)5&b^Ex5Lk#E&$+Wn_08lV-dZ`@hz;?hCe59yFSLUu|R#daJ-s%Wq$#drzAW# zvMM$sJH^ZRA~5Ot&`2z*Ck%hw&~>JVqhW*TgFu*msJ~YahT@^2aKZ@1`+GYhv1q8@ zKM)BCSz(DD81th8e}no_el?(5q}~PO0ak+;vZv)Q*nbu!UF*%5mWXsJrwjC zeu!rvkr3ek6b-T-@1cX8dW+Jc>=qH{o+Z$W*8T*H`~+m_T_v}MD;ad!OJpU-EA{tL z*&Y=(yjkK6@_mp#@$)VZ_lRmVBoJ6I;nc*4FPd-~qlhn?a9Je6Y}JI9b3{DqWITwO zySm4OoAtHHI~9w2L0OypRmDxlvb&#O?_t@8UVx`-TRY^CA4ca(3t31HT|gd($I=|I zXs@Nqb_1wAoiR$XbKKFOYuj10VcJyLD9WbV27vgqovT{v18s7(=E;(iH^K0)mBMi4 zWOf0|1N|=x{T7q{?5a~s-Oy%lKdL$AwAAa`+jo=Pe)Dg+{W}KOzmN74Z65=|k`HT> zZ9m7H56UyDwRGDbfLm;XkQayHaq{)DIRG4gxjBeQ$;CU_cCj4HjBOCy5NKenHu)g_ z?*k0JvU4Ywz6K7K`rt7*jqbHGcc!tbsqb9YQpp)D<-4e*dZ)c9^}ILLJMo5k1*B{# z3h(<3^(xQzK|ZZs)j$!d^?s?AR%ftkY39hJ)N4XCyHKi49fRgI%dV%@YhX6@z^~B} z$Sw;zEv4QPqREV-pm;8=UN2%fFGR&G7gk(ub$-S5xO!{FRjV!{3)ti89J0&EE)KdH zKzO2;3jsQT`0?4r*Z!T&g4WYx&|F$tkd&Ii8U~=g>I8)E`WZW$N$wN1UaQ%85P@$t z=u*_o>H0!qXfUYwnm2dO+kryr;H?6q&4A0^ zyzZP`F7&F>zH3_G9c`-mUGe^W@c!NeH|x_rXc4nuS_n$eJvQ&#Q1%?q8&DK<9_$7A zFHE|hya9dSzbj%nzli+qlJ;G<23a)vi{|LDCy1!gWud6K+fNy_*){KE=z3DU>VQ|i zU#)XJbLn5%?4)`H_$&KWkL#Ip!2jmvG@q-)G{+qE&i%OdYGzNKoj#5Iy`e58Y64Hl*U_+eH)cx_ee_Uu%{+j_% z1F!y(Z~%&ofg^Y*+&e!cDR4$&N+32e5{er3ru(0G952}CsGk=5K0(f9HzJlPko!CI zDYB<=MD#0CllmL=XL5k+Gmas$r*TH~nC!?{F6xjy_5XIj&^OuBxFt&C6jXSkaY4KW z85}cYPve9r%@-6u=@IpTJ&r-*X&e&u$bv9_ESut&FbhJ4>V_DNd!iWxAyN5cPxX(2 z%xq8dNRD!AVUKf{-F^%(jEm+unrDgDNP8k!mN_MWG24&IaZ8luJ+g7j4AJk}>}RAs z(I49iBs=nV)@D!Z0#TYTN#Ev<;ddE~pWa`w?`*~FFWT(K-3#8gU%sQ zyHok}9&drnTt0nYna4^&Y7iv%BzvN7fy8W2>n?o=x|jX$ZT7U@G{^v{916hBzXt-u zsQ&c5uDOUkwFCM4BV;iCW&$K7`$udz>S7{V3wbJ=3*_VLvi#d-5b|V4F!I#*2}4>T zCzy6Q9&sP_1;kMZiRK?1UxmuoKF8HqrOmm$e4$nN>4a!$Ju+)JA!1rzthkfWx?#mj z*+hXLove(}Ja5%s$7uex;^n#d*@{=>_HR}^|GPBGo$$q*=0~eOe5aM|ZN>Ay>k2;` za`yvP{TfE|j}?bOCC9fFcQcwdtavS>`M`>AgnJN=8RB^-+r()6nuSi>^~(LQ72nKg zyjk%rdEd!V?qsxEUem_@RLVZet|21V1*8POd)wNxzlgpLxP#@_HzNIHY<|v2eA$M- zAaVP-<|XcbtzI+^Ug3Ct`!7g8bWQ@6r#tX;;EZo;|0-};s2z&ccN6G4V3C=Z>y>oR zKF{&|c0hR{|BYNf|2^*Gcz!?Wl=?JJQn^Ptp5ISGfO|;EM9%58$N`-sy;2A`_(R!I(Z2i4FK8deDf)6S_Y+#^BU4 z8VMPppt5K0u3g}{fD*rd5~m@!12W>{Oq`#B$&rz9Ffnob!pfq`Z +#include + +void relax( + double* pos, + long* links, + double* mrel1, + double* mrel2, + double* lengths, + char* push, + char* pull, + int nlinks, + int iters) + { + int i, l, p1, p2; + double x1, x2, y1, y2, dx, dy, dist, change; +// printf("%d, %d\n", iters, nlinks); + for( i=0; i lengths[l] ) + dist = lengths[l]; + + change = (lengths[l]-dist) / dist; + dx *= change; + dy *= change; + + pos[p1] -= mrel2[l] * dx; + pos[p1+1] -= mrel2[l] * dy; + pos[p2] += mrel1[l] * dx; + pos[p2+1] += mrel1[l] * dy; + } + } +} \ No newline at end of file diff --git a/examples/verlet_chain/relax.py b/examples/verlet_chain/relax.py new file mode 100644 index 00000000..95a2b7b3 --- /dev/null +++ b/examples/verlet_chain/relax.py @@ -0,0 +1,23 @@ +import ctypes +import os + +so = os.path.join(os.path.dirname(__file__), 'maths.so') +lib = ctypes.CDLL(so) + +lib.relax.argtypes = [ + ctypes.c_void_p, + ctypes.c_void_p, + ctypes.c_void_p, + ctypes.c_void_p, + ctypes.c_void_p, + ctypes.c_void_p, + ctypes.c_void_p, + ctypes.c_int, + ctypes.c_int, + ] + +def relax(pos, links, mrel1, mrel2, lengths, push, pull, iters): + nlinks = links.shape[0] + lib.relax(pos.ctypes, links.ctypes, mrel1.ctypes, mrel2.ctypes, lengths.ctypes, push.ctypes, pull.ctypes, nlinks, iters) + + diff --git a/examples/verlet_chain_demo.py b/examples/verlet_chain_demo.py new file mode 100644 index 00000000..6ed97d48 --- /dev/null +++ b/examples/verlet_chain_demo.py @@ -0,0 +1,111 @@ +""" +Mechanical simulation of a chain using verlet integration. + + + +""" +import initExample ## Add path to library (just for examples; you do not need this) + +import pyqtgraph as pg +from pyqtgraph.Qt import QtCore, QtGui +import numpy as np + +from verlet_chain import ChainSim + +sim = ChainSim() + + +chlen1 = 80 +chlen2 = 60 +npts = chlen1 + chlen2 + +sim.mass = np.ones(npts) +sim.mass[chlen1-15] = 100 +sim.mass[chlen1-1] = 500 +sim.mass[npts-1] = 200 + +sim.fixed = np.zeros(npts, dtype=bool) +sim.fixed[0] = True +sim.fixed[chlen1] = True + +sim.pos = np.empty((npts, 2)) +sim.pos[:chlen1, 0] = 0 +sim.pos[chlen1:, 0] = 10 +sim.pos[:chlen1, 1] = np.arange(chlen1) +sim.pos[chlen1:, 1] = np.arange(chlen2) + +links1 = [(j, i+j+1) for i in range(chlen1) for j in range(chlen1-i-1)] +links2 = [(j, i+j+1) for i in range(chlen2) for j in range(chlen2-i-1)] +sim.links = np.concatenate([np.array(links1), np.array(links2)+chlen1, np.array([[chlen1-1, npts-1]])]) + +p1 = sim.pos[sim.links[:,0]] +p2 = sim.pos[sim.links[:,1]] +dif = p2-p1 +sim.lengths = (dif**2).sum(axis=1) ** 0.5 +sim.lengths[(chlen1-1):len(links1)] *= 1.05 # let auxiliary links stretch a little +sim.lengths[(len(links1)+chlen2-1):] *= 1.05 +sim.lengths[-1] = 7 + +push1 = np.ones(len(links1), dtype=bool) +push1[chlen1:] = False +push2 = np.ones(len(links2), dtype=bool) +push2[chlen2:] = False +sim.push = np.concatenate([push1, push2, np.array([True], dtype=bool)]) + +sim.pull = np.ones(sim.links.shape[0], dtype=bool) +sim.pull[-1] = False + +mousepos = sim.pos[0] + + +def display(): + global view, sim + view.clear() + view.addItem(sim.makeGraph()) + +def relaxed(): + global app + display() + app.processEvents() + +def mouse(pos): + global mousepos + pos = view.mapSceneToView(pos) + mousepos = np.array([pos.x(), pos.y()]) + +def update(): + global mousepos + #sim.pos[0] = sim.pos[0] * 0.9 + mousepos * 0.1 + s = 0.9 + sim.pos[0] = sim.pos[0] * s + mousepos * (1.0-s) + sim.update() + +app = pg.mkQApp() +win = pg.GraphicsLayoutWidget() +win.show() +view = win.addViewBox() +view.setAspectLocked(True) +view.setXRange(-100, 100) +#view.autoRange() + +view.scene().sigMouseMoved.connect(mouse) + +#display() +#app.processEvents() + +sim.relaxed.connect(relaxed) +sim.init() +sim.relaxed.disconnect(relaxed) + +sim.stepped.connect(display) + +timer = pg.QtCore.QTimer() +timer.timeout.connect(update) +timer.start(16) + + +## Start Qt event loop unless running in interactive mode or using pyside. +if __name__ == '__main__': + import sys + if (sys.flags.interactive != 1) or not hasattr(QtCore, 'PYQT_VERSION'): + QtGui.QApplication.instance().exec_()