diff --git a/Examples/fit/gisas2d/fit_gisas.py b/Examples/fit/gisas2d/fit_gisas.py
index e546facdceda6dd11f6542a4c5b34210553363c8..abd933db1e17ffbccbd38552499c07d6a625558e 100755
--- a/Examples/fit/gisas2d/fit_gisas.py
+++ b/Examples/fit/gisas2d/fit_gisas.py
@@ -7,7 +7,7 @@ Fake experimental data are generated by gisas_fake1.
 
 import gisas_model1 as model
 import bornagain as ba
-import ba_fitmonitor
+from bornagain import ba_fitmonitor
 import numpy as np
 import os
 from matplotlib import pyplot as plt
diff --git a/Examples/fit51_Basic/consecutive_fitting.py b/Examples/fit51_Basic/consecutive_fitting.py
index dff41607936a0347f5efedc676f4d8179a74d44f..ac78f42047b9c85ca599607cc7ad808e7276743e 100755
--- a/Examples/fit51_Basic/consecutive_fitting.py
+++ b/Examples/fit51_Basic/consecutive_fitting.py
@@ -10,8 +10,7 @@ after that to find precise minimum location.
 import numpy as np
 from matplotlib import pyplot as plt
 import bornagain as ba
-import ba_fitmonitor
-from bornagain import deg, angstrom, nm
+from bornagain import ba_fitmonitor, deg, angstrom, nm
 
 
 def get_sample(params):
@@ -78,7 +77,7 @@ def run_fitting():
     fit_objective = ba.FitObjective()
     fit_objective.addSimulationAndData(get_simulation, real_data, 1)
     fit_objective.initPrint(10)
-    observer = ba_fitmonitor.PlotterGISAS()
+    observer = bornagain.ba_fitmonitor.PlotterGISAS()
     fit_objective.initPlot(10, observer)
     """
     Setting fitting parameters with starting values.
diff --git a/Examples/fit52_Advanced/find_background.py b/Examples/fit52_Advanced/find_background.py
index 193d908131326e4cbea91dec173930769428242a..83a7c1e17db655eb98dbcc3cf98596162b242712 100755
--- a/Examples/fit52_Advanced/find_background.py
+++ b/Examples/fit52_Advanced/find_background.py
@@ -10,8 +10,7 @@ scale and background factors.
 import numpy as np
 from matplotlib import pyplot as plt
 import bornagain as ba
-import ba_fitmonitor
-from bornagain import deg, angstrom, nm
+from bornagain import deg, angstrom, nm, ba_fitmonitor
 
 
 def get_sample(params):
diff --git a/Examples/fit52_Advanced/fit_along_slices.py b/Examples/fit52_Advanced/fit_along_slices.py
index 0c7a786fb8dc96777353fb72843b53f2bc4736c1..29d734681f1bdecb7f1515cd11ac1d2178b25581 100755
--- a/Examples/fit52_Advanced/fit_along_slices.py
+++ b/Examples/fit52_Advanced/fit_along_slices.py
@@ -5,8 +5,7 @@ Fitting example: fit along slices
 
 from matplotlib import pyplot as plt
 import bornagain as ba
-import ba_plot
-from bornagain import deg, angstrom, nm
+from bornagain import deg, angstrom, nm, ba_plot
 
 phi_slice_value = 0.0  # position of vertical slice in deg
 alpha_slice_value = 0.2  # position of horizontal slice in deg
diff --git a/Examples/fit52_Advanced/fit_with_masks.py b/Examples/fit52_Advanced/fit_with_masks.py
index 6a11e39f14494aab4af75953f0a8dad0d0d424aa..5f3d8a3e15fe30e717bcba185326eb5448016f4a 100755
--- a/Examples/fit52_Advanced/fit_with_masks.py
+++ b/Examples/fit52_Advanced/fit_with_masks.py
@@ -6,8 +6,7 @@ Fitting example: fit with masks
 import numpy as np
 from matplotlib import pyplot as plt
 import bornagain as ba
-from bornagain import deg, angstrom, nm
-import ba_fitmonitor
+from bornagain import deg, angstrom, nm, ba_fitmonitor
 
 
 def get_sample(params):
diff --git a/Examples/fit52_Advanced/multiple_datasets.py b/Examples/fit52_Advanced/multiple_datasets.py
index 932cfbf74f190f914cfb0320158f402fc54f9ac9..9eb9272e1dc242ac665735bdff5e8cdb161a788f 100755
--- a/Examples/fit52_Advanced/multiple_datasets.py
+++ b/Examples/fit52_Advanced/multiple_datasets.py
@@ -8,7 +8,7 @@ import matplotlib
 from matplotlib import pyplot as plt
 import bornagain as ba
 from bornagain import deg, nm
-import ba_plot as bp
+import bornagain.ba_plot as bp
 
 def get_sample(params):
     """
diff --git a/Examples/fit54_ExternalMinimizer/lmfit_with_plotting.py b/Examples/fit54_ExternalMinimizer/lmfit_with_plotting.py
index 06e1c0c78f39eef48e334f07a628f9b54731e17c..2a42a0d0b2dd1d8cba98b55a2c850afad4e0f508 100755
--- a/Examples/fit54_ExternalMinimizer/lmfit_with_plotting.py
+++ b/Examples/fit54_ExternalMinimizer/lmfit_with_plotting.py
@@ -6,8 +6,7 @@ Fit progress is plotted using lmfit iteration calbback function.
 import numpy as np
 from matplotlib import pyplot as plt
 import bornagain as ba
-import ba_fitmonitor
-from bornagain import deg, angstrom, nm
+from bornagain import deg, angstrom, nm, ba_fitmonitor
 import lmfit
 
 
diff --git a/Examples/fit55_SpecularIntro/FitWithUncertainties.py b/Examples/fit55_SpecularIntro/FitWithUncertainties.py
index c2da703ee6a38e70e6c1de1fd1f84afaec6f1182..2720e256462f5d0111c4ad47eb0ce5b4369cfb29 100755
--- a/Examples/fit55_SpecularIntro/FitWithUncertainties.py
+++ b/Examples/fit55_SpecularIntro/FitWithUncertainties.py
@@ -13,8 +13,8 @@ and use genetic algorithm as the minimizer.
 """
 
 import bornagain as ba
-import ba_fitmonitor as bafim
-import ba_plot
+import bornagain.ba_fitmonitor as bafim
+from bornagain import ba_plot
 import numpy as np
 from matplotlib import pyplot as plt
 import os
diff --git a/Examples/fit55_SpecularIntro/RealLifeReflectometryFitting.py b/Examples/fit55_SpecularIntro/RealLifeReflectometryFitting.py
index 843ebfbef8a49ce43fc14056bba9e71652d0e61b..1882f8ecb1139908a34951dd11d5331814a15fc9 100755
--- a/Examples/fit55_SpecularIntro/RealLifeReflectometryFitting.py
+++ b/Examples/fit55_SpecularIntro/RealLifeReflectometryFitting.py
@@ -35,7 +35,7 @@ from matplotlib import pyplot as plt
 import numpy as np
 import os
 import bornagain as ba
-import ba_plot
+from bornagain import ba_plot
 from scipy.optimize import differential_evolution
 
 
diff --git a/Examples/fit61_Galaxi/fit_galaxi_data.py b/Examples/fit61_Galaxi/fit_galaxi_data.py
index c2dec2da8d6a409b0d201262b3d1e0447714ee42..ef24d65c633ee761d73a1c0aced3adf1761c8fcc 100755
--- a/Examples/fit61_Galaxi/fit_galaxi_data.py
+++ b/Examples/fit61_Galaxi/fit_galaxi_data.py
@@ -4,8 +4,7 @@ Fitting experimental data: spherical nanoparticles with size distribution
 in 3 layers system (experiment at GALAXI).
 """
 import bornagain as ba
-from bornagain import nm
-import ba_fitmonitor
+from bornagain import nm, ba_fitmonitor
 import os
 from matplotlib import pyplot as plt
 
diff --git a/Examples/scatter2d/ApproximationDA.py b/Examples/scatter2d/ApproximationDA.py
index 34bf608e39411bb720e440b68c03f38b87a4346e..c26515ae49e6499e889a3d767f082399a6a12278 100755
--- a/Examples/scatter2d/ApproximationDA.py
+++ b/Examples/scatter2d/ApproximationDA.py
@@ -58,7 +58,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/ApproximationLMA.py b/Examples/scatter2d/ApproximationLMA.py
index 8f9982be9f54dcfa4bae7b0c8893309821f1f569..4b0a9ccf393f0591c05263abb766101f5fcc6f8e 100755
--- a/Examples/scatter2d/ApproximationLMA.py
+++ b/Examples/scatter2d/ApproximationLMA.py
@@ -65,7 +65,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/ApproximationSSCA.py b/Examples/scatter2d/ApproximationSSCA.py
index a08caa88421399fc7a642519e5da4885f80464d3..d705710252a7a0ee82a1b6679cc244fd1d561722 100755
--- a/Examples/scatter2d/ApproximationSSCA.py
+++ b/Examples/scatter2d/ApproximationSSCA.py
@@ -59,7 +59,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/BeamDivergence.py b/Examples/scatter2d/BeamDivergence.py
index 08aa1a32b45d2791c4460152a904786e4ed0f3b2..a0f09254b9a892db628f4f5f5c2c1862d55b8fbc 100755
--- a/Examples/scatter2d/BeamDivergence.py
+++ b/Examples/scatter2d/BeamDivergence.py
@@ -57,7 +57,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/BiMaterialCylinders.py b/Examples/scatter2d/BiMaterialCylinders.py
index 98e31ac65ed655d7c1c174e9fa406930b7949cce..23630c243cc185d67845f6a605a9394c2662082a 100755
--- a/Examples/scatter2d/BiMaterialCylinders.py
+++ b/Examples/scatter2d/BiMaterialCylinders.py
@@ -62,7 +62,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/BoxesWithSpecularPeak.py b/Examples/scatter2d/BoxesWithSpecularPeak.py
index 57041e7b5e63a9f704a7c218d77b164bd6a828e7..58381fedd5cb70fca93fe3a3061d68cd85764fed 100755
--- a/Examples/scatter2d/BoxesWithSpecularPeak.py
+++ b/Examples/scatter2d/BoxesWithSpecularPeak.py
@@ -59,7 +59,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/BuriedParticles.py b/Examples/scatter2d/BuriedParticles.py
index 59c187d9935cdbfecdfc2f7fbfc6e5e15a478c53..a644e60d251bd6ba95ac2be198111f200562ff68 100755
--- a/Examples/scatter2d/BuriedParticles.py
+++ b/Examples/scatter2d/BuriedParticles.py
@@ -56,7 +56,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/ConstantBackground.py b/Examples/scatter2d/ConstantBackground.py
index d0e307e557f81df9da93971c61f7c60a15a3ee03..30f359f723aa61dcd77cda9424058fc4496b1f66 100755
--- a/Examples/scatter2d/ConstantBackground.py
+++ b/Examples/scatter2d/ConstantBackground.py
@@ -50,7 +50,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/CoreShellNanoparticles.py b/Examples/scatter2d/CoreShellNanoparticles.py
index 98b4942e2619407968c8418dbb827e6b364ac0b6..ca99f57f55a9a53a2f5490f91fed1bdb40552c1d 100755
--- a/Examples/scatter2d/CoreShellNanoparticles.py
+++ b/Examples/scatter2d/CoreShellNanoparticles.py
@@ -55,7 +55,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/CorrelatedRoughness.py b/Examples/scatter2d/CorrelatedRoughness.py
index fb911718ed38bfb287d7692799b3d18ea5156ee1..b208badd2e87fd88947037b4c1a18feef7bc88a7 100755
--- a/Examples/scatter2d/CorrelatedRoughness.py
+++ b/Examples/scatter2d/CorrelatedRoughness.py
@@ -51,7 +51,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/CosineRipplesAtRectLattice.py b/Examples/scatter2d/CosineRipplesAtRectLattice.py
index f8e2f64db5484f525785779e51e35ae2f854dbd0..01be4cfcb4c970380982fc4eaa8218074fd7f41b 100755
--- a/Examples/scatter2d/CosineRipplesAtRectLattice.py
+++ b/Examples/scatter2d/CosineRipplesAtRectLattice.py
@@ -59,7 +59,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/Cylinders.py b/Examples/scatter2d/Cylinders.py
index 4ebc23eadc386cc8dba31470a41ed47d96def09d..2f83107f473a486a9fe4a4af4d51b0639383d16e 100755
--- a/Examples/scatter2d/Cylinders.py
+++ b/Examples/scatter2d/Cylinders.py
@@ -49,7 +49,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/CylindersAndPrisms.py b/Examples/scatter2d/CylindersAndPrisms.py
index 8ef465bb0f4ec7ce3c31e8a6d9736dd3bfea9a80..4869cd978830b0eb2d8afc18a56b8d4c15c5db95 100755
--- a/Examples/scatter2d/CylindersAndPrisms.py
+++ b/Examples/scatter2d/CylindersAndPrisms.py
@@ -55,7 +55,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/CylindersInAverageLayer.py b/Examples/scatter2d/CylindersInAverageLayer.py
index c2775c51e0596d58726d9ef24569d584f322dddb..1c1bdfb5c8a4131867684bcb3fa6d5c830e169b2 100755
--- a/Examples/scatter2d/CylindersInAverageLayer.py
+++ b/Examples/scatter2d/CylindersInAverageLayer.py
@@ -51,7 +51,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/CylindersInBA.py b/Examples/scatter2d/CylindersInBA.py
index 649859b101339a2e29a34b8a22005e595561f77c..99464e08c7859ba8f6cdbdce024d71f4dba92854 100755
--- a/Examples/scatter2d/CylindersInBA.py
+++ b/Examples/scatter2d/CylindersInBA.py
@@ -46,7 +46,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/CylindersWithSizeDistribution.py b/Examples/scatter2d/CylindersWithSizeDistribution.py
index bd44f7119fa859b249d92190c2f71d0e154cc388..5b7606e798cb1a6c0f2653f5aa383eaef0b6e527 100755
--- a/Examples/scatter2d/CylindersWithSizeDistribution.py
+++ b/Examples/scatter2d/CylindersWithSizeDistribution.py
@@ -52,7 +52,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/DetectorResolutionFunction.py b/Examples/scatter2d/DetectorResolutionFunction.py
index 894b5143fd83e0e99b923b62b72917207aeff929..45589161e828bda2bf1233fccb1bdd68f7c85938 100755
--- a/Examples/scatter2d/DetectorResolutionFunction.py
+++ b/Examples/scatter2d/DetectorResolutionFunction.py
@@ -49,7 +49,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/DodecahedraSAS.py b/Examples/scatter2d/DodecahedraSAS.py
index da5cc69801a42b24599c6d8a2840a02ba81ea7c4..486ee10370528707cd35cc2f19d45fcdcee4d9c4 100755
--- a/Examples/scatter2d/DodecahedraSAS.py
+++ b/Examples/scatter2d/DodecahedraSAS.py
@@ -46,7 +46,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/HalfSpheresInAverageTopLayer.py b/Examples/scatter2d/HalfSpheresInAverageTopLayer.py
index 3e9344901724f6c977c48a05e2299e16348736a5..ff5e97c051ba4f2a33e04d59e7e15d10f2a503ad 100755
--- a/Examples/scatter2d/HalfSpheresInAverageTopLayer.py
+++ b/Examples/scatter2d/HalfSpheresInAverageTopLayer.py
@@ -60,7 +60,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/HexagonalLatticesWithBasis.py b/Examples/scatter2d/HexagonalLatticesWithBasis.py
index b9ba7c19e51e1179f831588b25425ca63cd3e18f..cce34241ea5537f5258d8bb972d7179c68839164 100755
--- a/Examples/scatter2d/HexagonalLatticesWithBasis.py
+++ b/Examples/scatter2d/HexagonalLatticesWithBasis.py
@@ -67,7 +67,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/Interference1DRadialParaCrystal.py b/Examples/scatter2d/Interference1DRadialParaCrystal.py
index a295456cf3119835b6b743d20bebb6ba2e6a40bc..6a4005248d60aa098e6d53795a0b9990ccc7da06 100755
--- a/Examples/scatter2d/Interference1DRadialParaCrystal.py
+++ b/Examples/scatter2d/Interference1DRadialParaCrystal.py
@@ -54,7 +54,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/Interference2DCenteredSquareLattice.py b/Examples/scatter2d/Interference2DCenteredSquareLattice.py
index 3d4d174f276106d30212bfc489fa72dc8b72d367..6da1f63d585d754616dcd2a0837cec7151d37455 100755
--- a/Examples/scatter2d/Interference2DCenteredSquareLattice.py
+++ b/Examples/scatter2d/Interference2DCenteredSquareLattice.py
@@ -65,7 +65,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/Interference2DLatticeSumOfRotated.py b/Examples/scatter2d/Interference2DLatticeSumOfRotated.py
index 20b12bd1822227a0a25ecf93d3db31a14a01e4dc..9bc14a70ed617ae8c68bb3fdeaae2e9a9c1c11cc 100755
--- a/Examples/scatter2d/Interference2DLatticeSumOfRotated.py
+++ b/Examples/scatter2d/Interference2DLatticeSumOfRotated.py
@@ -47,7 +47,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/Interference2DParaCrystal.py b/Examples/scatter2d/Interference2DParaCrystal.py
index d9193a04fe7aad732c8897ae60f599ccc85d529f..c117d174c5e904decc82118f0ce24017e57738e7 100755
--- a/Examples/scatter2d/Interference2DParaCrystal.py
+++ b/Examples/scatter2d/Interference2DParaCrystal.py
@@ -60,7 +60,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/Interference2DRotatedSquareLattice.py b/Examples/scatter2d/Interference2DRotatedSquareLattice.py
index da0bc3ed6a1226410758c482c5977a0c0116fe8e..1054b0f9279e1b3703f6feb7b30905e3d9159789 100755
--- a/Examples/scatter2d/Interference2DRotatedSquareLattice.py
+++ b/Examples/scatter2d/Interference2DRotatedSquareLattice.py
@@ -58,7 +58,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/Interference2DSquareFiniteLattice.py b/Examples/scatter2d/Interference2DSquareFiniteLattice.py
index 40edd4f2260bd3780c6e7a496b27252fed1357e8..3819c3271f7f94a04fec597b525044dffb1347db 100755
--- a/Examples/scatter2d/Interference2DSquareFiniteLattice.py
+++ b/Examples/scatter2d/Interference2DSquareFiniteLattice.py
@@ -56,7 +56,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/MagneticSpheres.py b/Examples/scatter2d/MagneticSpheres.py
index 9b3fdb868c5a82edde88c0f1de333823b51e0be8..eb0b28055abd49735dbcfbf218691855dc96cbf1 100755
--- a/Examples/scatter2d/MagneticSpheres.py
+++ b/Examples/scatter2d/MagneticSpheres.py
@@ -58,7 +58,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/MesoCrystal.py b/Examples/scatter2d/MesoCrystal.py
index bceecc20210ee6b4201d539bf151238214b2e03f..1aae130b97bb896f69770e41665be49a138dc8db 100755
--- a/Examples/scatter2d/MesoCrystal.py
+++ b/Examples/scatter2d/MesoCrystal.py
@@ -61,7 +61,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/ParticlesCrossingInterface.py b/Examples/scatter2d/ParticlesCrossingInterface.py
index 1aa31fd09bf2e6db061c79708d07273c01a79e8e..1a0c7c669fa52ba01574b2abcf423433b8cbdaf1 100755
--- a/Examples/scatter2d/ParticlesCrossingInterface.py
+++ b/Examples/scatter2d/ParticlesCrossingInterface.py
@@ -74,7 +74,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/PolarizedSANS.py b/Examples/scatter2d/PolarizedSANS.py
index 93f758a26dfa36115f0a4e6481e2f3f6ef242a6f..02d9138ed604b1bbd29c4763c8d3236d6dacc4c2 100755
--- a/Examples/scatter2d/PolarizedSANS.py
+++ b/Examples/scatter2d/PolarizedSANS.py
@@ -68,7 +68,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation) # TODO: restore units=ba.Axes.QSPACE)
diff --git a/Examples/scatter2d/PositionVariance.py b/Examples/scatter2d/PositionVariance.py
index 3612292b323b2ee7806d3a9216f5964d39a2e1c4..5ad816efd9b24e935aa4542952512c7cd65fdecb 100755
--- a/Examples/scatter2d/PositionVariance.py
+++ b/Examples/scatter2d/PositionVariance.py
@@ -69,7 +69,7 @@ def run_one(hasVariance, xi, nPlot, title):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     fig, axs = plt.subplots(3, 2, figsize=(10, 13))
 
     xi1 =  5*deg
diff --git a/Examples/scatter2d/RectangularGrating.py b/Examples/scatter2d/RectangularGrating.py
index 68dd1a5c2e640bd527e9e46521fb90c0583a14dc..c64f4317d006eb69d14339e550be0e93758e7309 100755
--- a/Examples/scatter2d/RectangularGrating.py
+++ b/Examples/scatter2d/RectangularGrating.py
@@ -63,7 +63,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/RotatedPyramids.py b/Examples/scatter2d/RotatedPyramids.py
index 2059fb009df29daac672ad084f415f54eb4d04d0..b70951f3ccf7d1d51cc6eecbad7342de9abdae84 100755
--- a/Examples/scatter2d/RotatedPyramids.py
+++ b/Examples/scatter2d/RotatedPyramids.py
@@ -50,7 +50,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/SpheresAtHexLattice.py b/Examples/scatter2d/SpheresAtHexLattice.py
index 1bcc4e12fe0fff2d5adee19ae2346a1131c57dfb..eb01ae4d3116a09abc389e4b5893c9f3eeab2b08 100755
--- a/Examples/scatter2d/SpheresAtHexLattice.py
+++ b/Examples/scatter2d/SpheresAtHexLattice.py
@@ -58,7 +58,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/TriangularRipple.py b/Examples/scatter2d/TriangularRipple.py
index 8ed97467ca50a599500fc0b303f8c88560bfaa4d..479f9328034fd5e7e80323a34cddeec31cef95db 100755
--- a/Examples/scatter2d/TriangularRipple.py
+++ b/Examples/scatter2d/TriangularRipple.py
@@ -59,7 +59,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/scatter2d/TwoTypesOfCylindersWithSizeDistribution.py b/Examples/scatter2d/TwoTypesOfCylindersWithSizeDistribution.py
index b22792a589ea5f59ee9644d5fd8d745cbdea7123..6d395efc44f42210a92483281e8c91aedacae0d6 100755
--- a/Examples/scatter2d/TwoTypesOfCylindersWithSizeDistribution.py
+++ b/Examples/scatter2d/TwoTypesOfCylindersWithSizeDistribution.py
@@ -61,7 +61,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/specular/AlternatingLayers.py b/Examples/specular/AlternatingLayers.py
index 9b957efe501cf00317797ce852b41b93a66d27ed..dfc717193a065b10d2c810cf20c22dc64dc0e3f2 100755
--- a/Examples/specular/AlternatingLayers.py
+++ b/Examples/specular/AlternatingLayers.py
@@ -40,7 +40,7 @@ def get_simulation(sample, scan_size=500):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/specular/BeamAngularDivergence.py b/Examples/specular/BeamAngularDivergence.py
index 815095463709989b3faf51e4537cdb9611aadcab..3f2e347d567a827f315d5f28c43459ad7b7d9725 100755
--- a/Examples/specular/BeamAngularDivergence.py
+++ b/Examples/specular/BeamAngularDivergence.py
@@ -90,7 +90,7 @@ if __name__ == '__main__':
     datadir = os.getenv('BORNAGAIN_EXAMPLE_DATA_DIR', '')
     data_fname = os.path.join(datadir, "genx_angular_divergence.dat.gz")
 
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/specular/BeamFullDivergence.py b/Examples/specular/BeamFullDivergence.py
index 79a342412cc2e4f4c64d5cb3932f3036203a6ac7..9b585dba78aaef7b852f989b983911a4c18c157c 100755
--- a/Examples/specular/BeamFullDivergence.py
+++ b/Examples/specular/BeamFullDivergence.py
@@ -72,7 +72,7 @@ def get_simulation(sample, scan_size=500):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/specular/TOFRWithResolution.py b/Examples/specular/TOFRWithResolution.py
index 654cd03059ff871ea9ed0f056c21f7dcb9573ef1..b698900634254c411169949a6e918ef19b30d746 100755
--- a/Examples/specular/TOFRWithResolution.py
+++ b/Examples/specular/TOFRWithResolution.py
@@ -66,7 +66,7 @@ def get_simulation(sample, scan_size=500):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/specular/TimeOfFlightReflectometry.py b/Examples/specular/TimeOfFlightReflectometry.py
index 7335778a2807c9fa7478aec2d400848a29537cd5..bb10942d71c6faabebf02ff5efddd3cc2f632690 100755
--- a/Examples/specular/TimeOfFlightReflectometry.py
+++ b/Examples/specular/TimeOfFlightReflectometry.py
@@ -55,7 +55,7 @@ def get_simulation(sample, scan_size=500):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/varia/AccessingSimulationResults.py b/Examples/varia/AccessingSimulationResults.py
index 5ae4b9e66cde9c40544a3c58236624e4cbe91578..928339d37e3ae239dbbc0cd8aa12ceddc39805a4 100755
--- a/Examples/varia/AccessingSimulationResults.py
+++ b/Examples/varia/AccessingSimulationResults.py
@@ -7,7 +7,7 @@ import math
 import random
 import bornagain as ba
 from bornagain import angstrom, deg, nm, nm2, kvector_t
-import ba_plot
+from bornagain import ba_plot
 from matplotlib import pyplot as plt
 from matplotlib import rcParams
 
diff --git a/Examples/varia/AllFormFactorsAvailable.py b/Examples/varia/AllFormFactorsAvailable.py
index fa762cbe53d622fb8290f3b5c34eff6b071e5425..adc5394cfa175ddcfdb0cde44367cebb1fce72ee 100755
--- a/Examples/varia/AllFormFactorsAvailable.py
+++ b/Examples/varia/AllFormFactorsAvailable.py
@@ -4,8 +4,7 @@ All form factors available in BornAgain in the Born Approximation
 """
 import numpy
 import bornagain as ba
-import ba_plot
-from bornagain import deg, angstrom
+from bornagain import deg, angstrom, ba_plot
 from matplotlib import pyplot as plt
 
 phi_min, phi_max = -2, 2.0
diff --git a/Examples/varia/AxesInDifferentUnits.py b/Examples/varia/AxesInDifferentUnits.py
index ec39f315c7042cb5d6ca8a77513bd7d89105dc57..2a2af592994f48e820538048b6a02f011d14036b 100755
--- a/Examples/varia/AxesInDifferentUnits.py
+++ b/Examples/varia/AxesInDifferentUnits.py
@@ -4,8 +4,7 @@ In this example we demonstrate how to plot simulation results with
 axes in different units (nbins, mm, degs and QyQz).
 """
 import bornagain as ba
-from bornagain import angstrom, deg, nm, nm2, kvector_t
-import ba_plot
+from bornagain import angstrom, deg, nm, nm2, kvector_t, ba_plot
 from matplotlib import pyplot as plt
 from matplotlib import rcParams
 
diff --git a/Examples/varia/CustomFormFactor.py b/Examples/varia/CustomFormFactor.py
index 2b0b439cbda626fe4b06ed15f3731184b6a82098..1f40b6ad6446ed68d9aef862315a2318c45e4e6f 100755
--- a/Examples/varia/CustomFormFactor.py
+++ b/Examples/varia/CustomFormFactor.py
@@ -80,7 +80,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Examples/varia/DepthProbe.py b/Examples/varia/DepthProbe.py
index 78e086348acb92c35fe1b1409351b5aaba9c78bf..00b8ce7b94172aab5afca3625e171749c0643f8a 100755
--- a/Examples/varia/DepthProbe.py
+++ b/Examples/varia/DepthProbe.py
@@ -82,7 +82,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation, aspect='auto')
diff --git a/Examples/varia/FindPeaks.py b/Examples/varia/FindPeaks.py
index 6f8effff5f06983eff45ce6bdc2304104ec48703..9aa6405b8f526cc5cd6ae844b0aeb7df6927c01a 100755
--- a/Examples/varia/FindPeaks.py
+++ b/Examples/varia/FindPeaks.py
@@ -5,8 +5,7 @@ Monte-carlo integration is used to get rid of
 large-particle form factor oscillations.
 """
 import bornagain as ba
-from bornagain import deg, angstrom, nm, micrometer
-import ba_plot
+from bornagain import deg, angstrom, nm, micrometer, ba_plot
 from matplotlib import pyplot as plt
 
 
diff --git a/Examples/varia/GratingMC.py b/Examples/varia/GratingMC.py
index 758754f5d5525bafd24b964731286affe43aca34..a45dd34dac35e2183d356f2122da8f751d0a6ebb 100755
--- a/Examples/varia/GratingMC.py
+++ b/Examples/varia/GratingMC.py
@@ -5,8 +5,7 @@ Monte-carlo integration is used to get rid of
 large-particle form factor oscillations.
 """
 import bornagain as ba
-from bornagain import deg, angstrom, nm, micrometer
-import ba_plot
+from bornagain import deg, angstrom, nm, micrometer, ba_plot
 from matplotlib import pyplot as plt
 
 
diff --git a/Examples/varia/Interference1DLattice.py b/Examples/varia/Interference1DLattice.py
index bb64d51b5828eaf60f936155b3182e96eabe9f5b..bc8a5aa71f2ee1822e6787123a96cdd75cf7d3f8 100755
--- a/Examples/varia/Interference1DLattice.py
+++ b/Examples/varia/Interference1DLattice.py
@@ -62,7 +62,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation, intensity_min=1e-03)
diff --git a/Examples/varia/LargeParticlesFormFactor.py b/Examples/varia/LargeParticlesFormFactor.py
index 7a04c0677fa8bb07877af49e307c33f1fd41c944..7164794cda0fac02d3a94c890adb28eeb29aef47 100755
--- a/Examples/varia/LargeParticlesFormFactor.py
+++ b/Examples/varia/LargeParticlesFormFactor.py
@@ -8,8 +8,7 @@ oscillates rapidly within one detector bin and analytical calculations
 In this case Monte-Carlo integration over detector bin should be used.
 """
 import bornagain as ba
-from bornagain import deg, angstrom, nm
-import ba_plot
+from bornagain import deg, angstrom, nm, ba_plot
 from matplotlib import pyplot as plt
 
 default_cylinder_radius = 10*nm
diff --git a/Examples/varia/OffSpecularSimulation.py b/Examples/varia/OffSpecularSimulation.py
index 5d19cb021a80d176080c1e12b174f8e28f736b49..b500ca643dc51124511b33b0b0633176d3f23a84 100755
--- a/Examples/varia/OffSpecularSimulation.py
+++ b/Examples/varia/OffSpecularSimulation.py
@@ -70,7 +70,7 @@ def get_simulation(sample):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation, intensity_min=1)
diff --git a/Examples/varia/SpecularSimulationWithRoughness.py b/Examples/varia/SpecularSimulationWithRoughness.py
index 5941d025d1a71a15809ad02e3554263cbc1a675a..fe7472671d56af03870eb5af16e8dbc19cbe2533 100755
--- a/Examples/varia/SpecularSimulationWithRoughness.py
+++ b/Examples/varia/SpecularSimulationWithRoughness.py
@@ -52,7 +52,7 @@ def get_simulation(sample, scan_size=500):
 
 
 if __name__ == '__main__':
-    import ba_plot
+    from bornagain import ba_plot
     sample = get_sample()
     simulation = get_simulation(sample)
     ba_plot.run_and_plot(simulation)
diff --git a/Wrap/Python/__init__.py.in b/Wrap/Python/__init__.py.in
index 7dbac77823ede7609380d02511becfe1ce6efbac..01e522db3128b006ada25a2b9c495ed808b01873 100644
--- a/Wrap/Python/__init__.py.in
+++ b/Wrap/Python/__init__.py.in
@@ -3,9 +3,6 @@ import sys, os
 sys.path.append(os.path.abspath(
                 os.path.join(os.path.split(__file__)[0], '@BA_MODULES_IMPORT_PATH@')))
 
-# this line needed to fix python3 import problem
-sys.path.append(os.path.abspath(os.path.dirname(__file__)))
-
 # this is needed to adapt to the changes in Python 3.8 on Windows regarding dll loading
 # see https://docs.python.org/3/whatsnew/3.8.html#ctypes
 if sys.version_info >= (3, 8, 0) and sys.platform == 'win32':
diff --git a/Wrap/Python/ba_fitmonitor.py b/Wrap/Python/ba_fitmonitor.py
index 33b14cd73031073c77b3b19ed9e68b9f5260f347..fb3ab0f4f8ea450626a95fca00947903c9153a83 100644
--- a/Wrap/Python/ba_fitmonitor.py
+++ b/Wrap/Python/ba_fitmonitor.py
@@ -13,7 +13,7 @@
 #  **************************************************************************  #
 
 import bornagain as ba
-import ba_plot as bp
+from bornagain import ba_plot as bp
 try:  # workaround for build servers
     import numpy as np
     from matplotlib import pyplot as plt