Commit 553eff5b authored by Whitfield, Ross's avatar Whitfield, Ross
Browse files

Update IntegratePeaksMD_ellipsoid_test.py PlotPeakByLofValue_test.py...

Update IntegratePeaksMD_ellipsoid_test.py PlotPeakByLofValue_test.py WANDPowderReduction_group_test.py calgoniotest.py centroidmd_leanelasticpeak.py compare_LeakElasticPeaks_test.py ellipsoid_shapes.py ellipsoid_shell.py ellipsoid_shell3.py etest.py etest_fake.py etest_fake2.py etest_fake3.py event_test.py findleanelasticpeakstest.py goniometer_test.py goniometer_test_long.py indexpeakstest.py leanpeaksworkpsacetest.py multig.py peak_ellipsoid_integrate_test.py peaks_test.py peaksworkspacesaveload.py predictpeakstest.py quick_peak_test.py
parent 631038a0
# import mantid algorithms, numpy and matplotlib
from mantid.simpleapi import *
import matplotlib.pyplot as plt
import numpy as np
ws = CreateMDWorkspace(Dimensions='3', Extents='0.5,1.5,0.5,1.5,0.5,1.5',
Names='Q_lab_x,Q_lab_y,Q_lab_z', Units='U,U,U',
Frames='QLab,QLab,QLab',
SplitInto='2', SplitThreshold='50')
expt_info = CreateSampleWorkspace()
ws.addExperimentInfo(expt_info)
nevents=10000000
FakeMDEventData(ws, UniformParams=f'{nevents},0.98,1.02,0.98,1.02,0.98,1.02')
# add peak
p = CreatePeaksWorkspace(InstrumentWorkspace='ws', NumberOfPeaks=0, OutputWorkspace='peaks')
SetUB('peaks', 1,1,1,90,90,90)
pk = p.createPeak([1,1,1]) # axes aligned with viewing axes
p.addPeak(pk)
a = 0.005
b = 0.013
c = 0.013
#a=b=c=0.005
p=IntegratePeaksMD(InputWorkspace='ws', PeakRadius=[a,b,c], Ellipsoid=True, PeaksWorkspace='peaks', OutputWorkspace='ellipsoid')
print("Integrate Intensity =", p.getPeak(0).getIntensity())
print("Approx. Expected Intensity=", 4/3*np.pi*(a*b*c)*nevents*25**3)
from mantid.simpleapi import *
import numpy as np
ws = CreateSampleWorkspace()
#create string of workspaces to fit (ws,i0; ws,i1, ws,i2 ...)
workspaces = [ws.name() + ',i%d' % i for i in range(ws.getNumberHistograms())]
workspaces = ';'.join(workspaces)
function = "name=Gaussian,Height=10.0041,PeakCentre=10098.6,Sigma=48.8581"
peaks = PlotPeakByLogValue(workspaces, function, Spectrum=1)
function = "name=Gaussian,Height=10.0041,PeakCentre=10098.6,Sigma=48.8581;name=FlatBackground,A0=0.3"
peaks = PlotPeakByLogValue(workspaces, function, Spectrum=1)
# import mantid algorithms, numpy and matplotlib
from mantid.simpleapi import *
import matplotlib.pyplot as plt
import numpy as np
from mantid.simpleapi import (
WANDPowderReduction,
CreateSampleWorkspace,CloneWorkspace, GroupWorkspaces
)
from mantid.api import WorkspaceGroup
import numpy as np
event_data = CreateSampleWorkspace(
NumBanks=1,
BinWidth=20000,
PixelSpacing=0.1,
BankPixelWidth=100,
WorkspaceType="Event",
)
event_cal = CreateSampleWorkspace(
NumBanks=1,
BinWidth=20000,
PixelSpacing=0.1,
BankPixelWidth=100,
WorkspaceType="Event",
Function="Flat background",
)
event_bkg = CreateSampleWorkspace(
NumBanks=1,
BinWidth=20000,
PixelSpacing=0.1,
BankPixelWidth=100,
WorkspaceType="Event",
Function="Flat background",
)
pd_out = WANDPowderReduction(
InputWorkspace=[event_data, event_data],
CalibrationWorkspace=event_cal,
BackgroundWorkspace=event_bkg,
Target="Theta",
NumberBins=1000,
NormaliseBy="None",
Sum=True,
)
pd_out2 = WANDPowderReduction(
InputWorkspace=[event_data, event_data],
CalibrationWorkspace=event_cal,
BackgroundWorkspace=event_bkg,
Target="Theta",
NumberBins=1000,
NormaliseBy="None",
Sum=True,
)
event_data2 = CloneWorkspace(event_data)
event_data_group2 = WorkspaceGroup()
event_data_group2.addWorkspace(event_data)
event_data_group2.addWorkspace(event_data2)
pd_out = WANDPowderReduction(
InputWorkspace=event_data_group2,
CalibrationWorkspace=event_cal,
BackgroundWorkspace=event_bkg,
Target="Theta",
NumberBins=1000,
NormaliseBy="None",
Sum=False,
)
event_data_group = GroupWorkspaces([event_data,event_data2])
pd_out = WANDPowderReduction(
InputWorkspace=event_data_group,
CalibrationWorkspace=event_cal,
BackgroundWorkspace=event_bkg,
Target="Theta",
NumberBins=1000,
NormaliseBy="None",
Sum=False,
)
\ No newline at end of file
from mantid.geometry import Goniometer
from mantid.kernel import V3D, FloatTimeSeriesProperty
import numpy as np
from mantid.simpleapi import CreatePeaksWorkspace, HFIRCalculateGoniometer, AddTimeSeriesLog, LoadEmptyInstrument
omega = np.deg2rad(42)
R = np.array([[np.cos(omega), 0, np.sin(omega)],
[0, 1, 0],
[-np.sin(omega), 0, np.cos(omega)]])
wl = 1.54
k = 2*np.pi/wl
theta = np.deg2rad(47)
phi = np.deg2rad(13)
q_lab = np.array([-np.sin(theta)*np.cos(phi),
-np.sin(theta)*np.sin(phi),
1 - np.cos(theta)]) * k
q_sample = np.dot(np.linalg.inv(R), q_lab)
peaks = CreatePeaksWorkspace(OutputType="LeanElasticPeak", NumberOfPeaks=0)
p = peaks.createPeakQSample(q_sample)
peaks.addPeak(p)
HFIRCalculateGoniometer(peaks, wl, OverrideProperty=True, InnerGoniometer=True)
g = Goniometer()
g.setR(peaks.getPeak(0).getGoniometerMatrix())
print(g.getEulerAngles('YZY'))
assert np.isclose(g.getEulerAngles('YZY')[0], 42)
chi = np.deg2rad(-3)
phi = np.deg2rad(23)
R1 = np.array([[np.cos(omega), 0, -np.sin(omega)], # omega 0,1,0,-1
[ 0, 1, 0],
[np.sin(omega), 0, np.cos(omega)]])
R2 = np.array([[ np.cos(chi), np.sin(chi), 0], # chi 0,0,1,-1
[-np.sin(chi), np.cos(chi), 0],
[ 0, 0, 1]])
R3 = np.array([[np.cos(phi), 0, -np.sin(phi)], # phi 0,1,0,-1
[ 0, 1, 0],
[np.sin(phi), 0, np.cos(phi)]])
R = np.dot(np.dot(R1, R2), R3)
q_sample = np.dot(np.linalg.inv(R), q_lab)
peaks = CreatePeaksWorkspace(OutputType="LeanElasticPeak", NumberOfPeaks=0)
p = peaks.createPeakQSample(q_sample)
p.setGoniometerMatrix(np.dot(R2, R3))
peaks.addPeak(p)
HFIRCalculateGoniometer(peaks, wl)
g = Goniometer()
g.setR(peaks.getPeak(0).getGoniometerMatrix())
print(g.getEulerAngles('YZY'))
assert np.isclose(g.getEulerAngles('YZY')[0], -42)
assert np.isclose(g.getEulerAngles('YZY')[1], 3)
assert np.isclose(g.getEulerAngles('YZY')[2], -23)
peaks = CreatePeaksWorkspace(NumberOfPeaks=0, OutputType="LeanElasticPeak")
#peaks.run().getGoniometer().setR(np.dot(R1, R2))
peaks.run().getGoniometer().setR(R)
p = peaks.createPeakQSample(q_sample)
peaks.addPeak(p)
HFIRCalculateGoniometer(peaks, wl, OverrideProperty=True, InnerGoniometer=True)
g = Goniometer()
g.setR(peaks.getPeak(0).getGoniometerMatrix())
print(g.getEulerAngles('YZY'))
assert np.isclose(g.getEulerAngles('YZY')[0], -42)
assert np.isclose(g.getEulerAngles('YZY')[1], 3)
assert np.isclose(g.getEulerAngles('YZY')[2], -23)
from mantid.simpleapi import *
peaks = CreatePeaksWorkspace(NumberOfPeaks=0)
p = peaks.createPeakQSample((1,1,1)) # starting at [1,1,1]
p.setWavelength(1.54)
peaks.addPeak(p)
ws = CreateMDWorkspace(Dimensions='3', Extents='0.5,1.5,0.5,1.5,0.5,1.5',
Names='x,y,z', Units='rlu,rlu,rlu',
Frames='QSample,QSample,QSample')
# Create a fake peak at [1.1, 0.9, 1.05]
FakeMDEventData(ws, PeakParams='1000000,1.1,0.9,1.05,0.2')
centroid_peaks = CentroidPeaksMD(ws, peaks)
print("Qsample (should be approx [1.1, 0.9, 1.05]) =",centroid_peaks.getPeak(0).getQSampleFrame())
print("Wavelength =",centroid_peaks.getPeak(0).getWavelength())
from mantid.simpleapi import *
from mantid.geometry import OrientedLattice
from itertools import permutations
import numpy as np
ol=OrientedLattice(5,7,12,90,90,120)
ub=ol.getUB()
print(ub)
peaks = CreatePeaksWorkspace(NumberOfPeaks=0)
peaks.addPeak(peaks.createPeakQSample(2*np.pi*np.dot(ub,[1,1,1])))
peaks.addPeak(peaks.createPeakQSample(2*np.pi*np.dot(ub,[1,1,0])))
peaks.addPeak(peaks.createPeakQSample(2*np.pi*np.dot(ub,[1,2,0])))
peaks2 = CreatePeaksWorkspace(NumberOfPeaks=0)
peaks2.addPeak(peaks2.createPeakQSample(2*np.pi*np.dot(ub,[1,1,1])))
peaks2.addPeak(peaks2.createPeakQSample(2*np.pi*np.dot(ub,[1,1,0])))
peaks2.addPeak(peaks2.createPeakQSample(2*np.pi*np.dot(ub,[1,2,0])))
CompareWorkspaces(peaks, peaks2, Tolerance=1e-4)
ws = CreateSampleWorkspace()
peaks = CreatePeaksWorkspace(InstrumentWorkspace='ws')
peaks2 = CreatePeaksWorkspace(InstrumentWorkspace='ws')
peaks.getPeak(0).setRunNumber(123)
CompareWorkspaces(peaks, peaks2, Tolerance=1e-4)
# import mantid algorithms, numpy and matplotlib
from mantid.simpleapi import *
import matplotlib.pyplot as plt
import numpy as np
ws = CreateMDWorkspace(Dimensions='3', Extents='-2,2,-2,2,-2,2',
Names='Q_lab_x,Q_lab_y,Q_lab_z', Units='U,U,U',
Frames='QLab,QLab,QLab',
SplitInto='2', SplitThreshold='50')
expt_info = CreateSampleWorkspace()
ws.addExperimentInfo(expt_info)
# add peak
p = CreatePeaksWorkspace(InstrumentWorkspace='ws', NumberOfPeaks=0, OutputWorkspace='peaks')
SetUB('peaks', 1,1,1,90,90,90)
pk = p.createPeak([1,1,1]) # axes aligned with viewing axes
p.addPeak(pk)
a=0.3
b=0.4
c=0.5
IntegratePeaksMD(InputWorkspace='ws', PeakRadius=[a,b,c], Ellipsoid=True, PeaksWorkspace='peaks', OutputWorkspace='ellipsoid')
a=0.4
b=0.3
c=0.3
b_in_a = 0.4
b_in_b = 0.4
b_in_c = 0.4
b_out_a = 0.4
b_out_b = 0.5
b_out_c = 0.5
IntegratePeaksMD(InputWorkspace='ws', PeakRadius=[a,b,c], BackgroundInnerRadius=[b_in_a, b_in_b, b_in_c], BackgroundOuterRadius=[b_out_a, b_out_b, b_out_c], Ellipsoid=True, PeaksWorkspace='peaks', OutputWorkspace='ellipsoid_bg')
IntegratePeaksMD(InputWorkspace='ws',
PeakRadius=0.3,
BackgroundInnerRadius=0.4,
BackgroundOuterRadius=0.5,
PeaksWorkspace='peaks', OutputWorkspace='sphere_bg')
import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.patches import Circle, Ellipse, Patch, PathPatch, Wedge
from matplotlib.transforms import Affine2D, IdentityTransform
import numpy as np
class EllipticalShell(Patch):
"""
Elliptical shell patch.
"""
def __str__(self):
return f"EllipticalShell(center={self.center}, width={self.width}, height={self.height}, thick={self.thick}, angle={self.angle})"
def __init__(self, center, width, height, thick, angle=0.0, **kwargs):
"""
Draw an elliptical ring centered at *x*, *y* center with outer width (horizontal diameter)
*width* and outer height (vertical diameter) *height* with a fractional ring thickness of *thick*
( [r_outer - r_inner]/r_outer where r is the major radius).
Valid kwargs are:
%(Patch)s
"""
super().__init__(**kwargs)
self.center = center
self.height, self.width = height, width
self.thick = thick
self.angle = angle
self._recompute_path()
# Note: This cannot be calculated until this is added to an Axes
self._patch_transform = IdentityTransform()
def _recompute_path(self):
# Form the outer ring
arc = Path.arc(theta1=0.0, theta2=360.0)
print(f'{arc=}')
# Draw the outer unit circle followed by a reversed and scaled inner circle
v1 = arc.vertices
v2 = arc.vertices[::-1] * float(1.0 - self.thick) # self.thick is fractional thickness
print(f'{v1=} {v2=}')
v = np.vstack([v1, v2, v1[0, :], (0, 0)])
print(f'{v=}')
c = np.hstack([arc.codes, arc.codes, Path.MOVETO, Path.CLOSEPOLY])
print(f'{c=}')
c[len(arc.codes)] = Path.MOVETO
# Final shape acheieved through axis transformation. See _recompute_transform
self._path = Path(v, c)
def _recompute_transform(self):
"""NOTE: This cannot be called until after this has been added
to an Axes, otherwise unit conversion will fail. This
makes it very important to call the accessor method and
not directly access the transformation member variable.
"""
center = (self.convert_xunits(self.center[0]), self.convert_yunits(self.center[1]))
width = self.convert_xunits(self.width)
height = self.convert_yunits(self.height)
self._patch_transform = Affine2D() \
.scale(width * 0.5, height * 0.5) \
.rotate_deg(self.angle) \
.translate(*center)
def get_patch_transform(self):
self._recompute_transform()
return self._patch_transform
def get_path(self):
if self._path is None:
self._recompute_path()
return self._path
fig, ax = plt.subplots(subplot_kw={'aspect': 'equal'})
e3 = EllipticalShell(center=(0, 0),
width=1,
height=1,
thick=0.5,
angle=0)
ax.add_patch(e3)
e3.set_facecolor((1, 0, 0))
e1 = Ellipse(xy=(0, 0),
width=0.1,
height=1,
angle=0)
ax.add_patch(e1)
e1.set_facecolor((0, 0, 1))
e2 = Ellipse(xy=(0, 0),
width=0.1,
height=1,
angle=90)
ax.add_patch(e2)
e2.set_facecolor((0, 1, 0))
ax.set_xlim(-1, 1)
ax.set_ylim(-1, 1)
arc = Path.arc(theta1=0.0, theta2=360.0)
v1 = arc.vertices
v2 = arc.vertices[::-1] * float(1.0 - 0.5) # self.thick is fractional thickness
ax.scatter(v1[:,0], v1[:,1])
ax.scatter(v2[:,0], v2[:,1])
plt.show()
import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.patches import Circle, Ellipse, Patch, PathPatch, Wedge
from matplotlib.transforms import Affine2D, IdentityTransform
import numpy as np
class EllipticalShell(Patch):
"""
Elliptical shell patch.
"""
def __str__(self):
return f"EllipticalShell(center={self.center}, width={self.width}, height={self.height}, thick={self.thick}, angle={self.angle})"
def __init__(self, center, width, height, thick, angle=0.0, **kwargs):
"""
Draw an elliptical ring centered at *x*, *y* center with outer width (horizontal diameter)
*width* and outer height (vertical diameter) *height* with a fractional ring thickness of *thick*
( [r_outer - r_inner]/r_outer where r is the major radius).
Valid kwargs are:
%(Patch)s
"""
super().__init__(**kwargs)
self.center = center
self.height, self.width = height, width
self.thick = thick
self.angle = angle
self._recompute_path()
# Note: This cannot be calculated until this is added to an Axes
self._patch_transform = IdentityTransform()
def _recompute_path(self):
# Form the outer ring
arc = Path.arc(theta1=0.0, theta2=360.0)
print(f'{arc=}')
# Draw the outer unit circle followed by a reversed and scaled inner circle
v1 = arc.vertices
v2 = np.zeros_like(v1)
v2[:, 0] = v1[::-1, 0] * float(1.0 - self.thick[0])
v2[:, 1] = v1[::-1, 1] * float(1.0 - self.thick[1])
print(f'{v1=} {v2=}')
v = np.vstack([v1, v2, v1[0, :], (0, 0)])
print(f'{v=}')
c = np.hstack([arc.codes, arc.codes, Path.MOVETO, Path.CLOSEPOLY])
print(f'{c=}')
c[len(arc.codes)] = Path.MOVETO
# Final shape acheieved through axis transformation. See _recompute_transform
self._path = Path(v, c)
def _recompute_transform(self):
"""NOTE: This cannot be called until after this has been added
to an Axes, otherwise unit conversion will fail. This
makes it very important to call the accessor method and
not directly access the transformation member variable.
"""
center = (self.convert_xunits(self.center[0]), self.convert_yunits(self.center[1]))
width = self.convert_xunits(self.width)
height = self.convert_yunits(self.height)
self._patch_transform = Affine2D() \
.scale(width * 0.5, height * 0.5) \
.rotate_deg(self.angle) \
.translate(*center)
def get_patch_transform(self):
self._recompute_transform()
return self._patch_transform
def get_path(self):
if self._path is None:
self._recompute_path()
return self._path
fig, ax = plt.subplots(subplot_kw={'aspect': 'equal'})
e3 = EllipticalShell(center=(0, 0),
width=1,
height=1,
thick=(0.5, 0.0),
angle=0)
ax.add_patch(e3)
e3.set_facecolor((1, 0, 0))
e1 = Ellipse(xy=(0, 0),
width=0.1,
height=1,
angle=0)
ax.add_patch(e1)
e1.set_facecolor((0, 0, 1))
e2 = Ellipse(xy=(0, 0),
width=0.1,
height=1,
angle=90)
ax.add_patch(e2)
e2.set_facecolor((0, 1, 0))
ax.set_xlim(-1, 1)
ax.set_ylim(-1, 1)
arc = Path.arc(theta1=0.0, theta2=360.0)
v1 = arc.vertices
v2 = arc.vertices[::-1] * float(1.0 - 0.5) # self.thick is fractional thickness
ax.scatter(v1[:,0], v1[:,1])
ax.scatter(v2[:,0], v2[:,1])
plt.show()
# import mantid algorithms, numpy and matplotlib
from mantid.simpleapi import *
import matplotlib.pyplot as plt
import numpy as np
Load(Filename='SXD23767.raw', OutputWorkspace='SXD23767', LoaderName='LoadRaw', LoaderVersion='3')
ConvertToDiffractionMDWorkspace(InputWorkspace='SXD23767', OutputWorkspace='SXD23767MD', OneEventPerBin='0')
FindSXPeaks(InputWorkspace='SXD23767', PeakFindingStrategy='AllPeaks', AbsoluteBackground=2000,
ResolutionStrategy='AbsoluteResolution', XResolution=200, PhiResolution=3, TwoThetaResolution=3,
OutputWorkspace='peaks')
IntegratePeaksMD(InputWorkspace='SXD23767MD', PeakRadius=0.1, Ellipsoid=True, PeaksWorkspace='peaks', OutputWorkspace='peaks_int',
UseOnePercentBackgroundCorrection=False)
IntegratePeaksMD(InputWorkspace='SXD23767MD', PeakRadius='0.1,0.2,0.3', Ellipsoid=True, PeaksWorkspace='peaks', OutputWorkspace='peaks_int2',
UseOnePercentBackgroundCorrection=False)
IntegratePeaksMD(InputWorkspace='SXD23767MD', PeakRadius=0.1, PeaksWorkspace='peaks', OutputWorkspace='peaks_ints',
UseOnePercentBackgroundCorrection=False)
IntegratePeaksMD(InputWorkspace='SXD23767MD', PeakRadius='0.2,0.1,0.1', Ellipsoid=True, PeaksWorkspace='peaks', OutputWorkspace='peaks_intbg',
UseOnePercentBackgroundCorrection=False, BackgroundInnerRadius='0.2,0.2,0.1', BackgroundOuterRadius='0.2,0.3,0.2')
IntegratePeaksMD(InputWorkspace='SXD23767MD', PeakRadius='0.1,0.2,0.3', Ellipsoid=True, PeaksWorkspace='peaks', OutputWorkspace='peaks_intbg2',
UseOnePercentBackgroundCorrection=False, BackgroundInnerRadius='0.1,0.2,0.3', BackgroundOuterRadius='0.2,0.3,0.4')
IntegratePeaksMD(InputWorkspace='SXD23767MD', PeakRadius=0.1, Ellipsoid=False, PeaksWorkspace='peaks', OutputWorkspace='peaks_int01',
UseOnePercentBackgroundCorrection=False)
IntegratePeaksMD(InputWorkspace='SXD23767MD', PeakRadius='0.1,0.1,0.1', Ellipsoid=True, PeaksWorkspace='peaks', OutputWorkspace='peaks_int02',
UseOnePercentBackgroundCorrection=False)
IntegratePeaksMD(InputWorkspace='SXD23767MD', PeakRadius=0.1, Ellipsoid=False, PeaksWorkspace='peaks', OutputWorkspace='peaks_int11',
UseOnePercentBackgroundCorrection=False, BackgroundInnerRadius=0.15, BackgroundOuterRadius=0.2)
IntegratePeaksMD(InputWorkspace='SXD23767MD', PeakRadius='0.1,0.1,0.1', Ellipsoid=True, PeaksWorkspace='peaks', OutputWorkspace='peaks_int12',
UseOnePercentBackgroundCorrection=False, BackgroundInnerRadius='0.15,0.15,0.15', BackgroundOuterRadius='0.2,0.2,0.2')
IntegratePeaksMD(InputWorkspace='SXD23767MD', PeakRadius='0.1,0.1,0.1', Ellipsoid=True, PeaksWorkspace='peaks', OutputWorkspace='peaks_int03',
UseOnePercentBackgroundCorrection=False, BackgroundInnerRadius='0.15,0.15,0.15', BackgroundOuterRadius='0.15,0.15,0.15')
IntegratePeaksMD(InputWorkspace='SXD23767MD', PeakRadius='0.1,0.1,0.1', Ellipsoid=True, PeaksWorkspace='peaks', OutputWorkspace='peaks_int04',
UseOnePercentBackgroundCorrection=False, BackgroundOuterRadius='0.1,0.1,0.1')
IntegratePeaksMD(InputWorkspace='SXD23767MD', PeakRadius=0.1, Ellipsoid=False, PeaksWorkspace='peaks', OutputWorkspace='peaks_int05',
UseOnePercentBackgroundCorrection=False, BackgroundInnerRadius=0.15, BackgroundOuterRadius=0.15)
IntegratePeaksMD(InputWorkspace='SXD23767MD', PeakRadius='0.1,0.1,0.1', Ellipsoid=True, PeaksWorkspace='peaks', OutputWorkspace='peaks_int06',
UseOnePercentBackgroundCorrection=False, BackgroundOuterRadius='0.05,0.05,0.05')
IntegratePeaksMD(InputWorkspace='SXD23767MD', PeakRadius='0.1,0.1,0.1', Ellipsoid=True, PeaksWorkspace='peaks', OutputWorkspace='peaks_int21',
UseOnePercentBackgroundCorrection=False, BackgroundOuterRadius='0.2,0.05,0.2')
# import mantid algorithms, numpy and matplotlib
from mantid.simpleapi import *
import matplotlib.pyplot as plt
import numpy as np