diff --git a/.gitignore b/.gitignore
index cd7c206..d670df7 100644
--- a/.gitignore
+++ b/.gitignore
@@ -6,6 +6,7 @@ pyext/build/
 .libs
 .deps
 test/*.fits
+examples/python/*.fits
 inst/cta/test/*.xml
 inst/cta/test/*.fits
 inst/cta/test/*.fits.gz
diff --git a/examples/python/example_model_xml_io.py b/examples/python/example_model_xml_io.py
index 01c6ac6..6fafa2e 100755
--- a/examples/python/example_model_xml_io.py
+++ b/examples/python/example_model_xml_io.py
@@ -18,13 +18,14 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 # =====================================================================
-from gammalib import *
+import gammalib
+from gammalib import GEnergy
 
 #
 # Set the location / centre of the source model to
 # Right Ascension = 42.0 deg and Declination = 43 deg
 # ===================================================
-centre = GSkyDir()
+centre = gammalib.GSkyDir()
 centre.radec_deg(42.0, 43.0)
 
 #
@@ -32,9 +33,9 @@ centre.radec_deg(42.0, 43.0)
 # Normalization:   1.0e-7 ph/cm2/s/MeV @ 100 MeV
 # Spectral index: -2.1
 # ===================================================
-point_spatial  = GModelSpatialPointSource(centre)
-point_spectrum = GModelSpectralPlaw(1.0e-7, -2.1, GEnergy(100.0, "MeV"))
-point = GModelSky(point_spatial, point_spectrum)
+point_spatial  = gammalib.GModelSpatialPointSource(centre)
+point_spectrum = gammalib.GModelSpectralPlaw(1.0e-7, -2.1, GEnergy(100.0, "MeV"))
+point = gammalib.GModelSky(point_spatial, point_spectrum)
 point.name('My point source')
 
 #
@@ -43,9 +44,9 @@ point.name('My point source')
 # Normalization:   4.2e-7 ph/cm2/s/MeV @ 100 MeV
 # Spectral index: -2.4
 # ===================================================
-gauss_spatial  = GModelSpatialRadialGauss(centre, 3.0)
-gauss_spectrum = GModelSpectralPlaw(4.2e-7, -2.4, GEnergy(100.0, "MeV"))
-gauss = GModelSky(gauss_spatial, gauss_spectrum)
+gauss_spatial  = gammalib.GModelSpatialRadialGauss(centre, 3.0)
+gauss_spectrum = gammalib.GModelSpectralPlaw(4.2e-7, -2.4, GEnergy(100.0, "MeV"))
+gauss = gammalib.GModelSky(gauss_spatial, gauss_spectrum)
 gauss.name('My Gaussian source')
 
 #
@@ -53,7 +54,7 @@ gauss.name('My Gaussian source')
 # models to the model container. Print the content of the
 # container
 # =========================================================
-models = GModels()
+models = gammalib.GModels()
 models.append(point)
 models.append(gauss)
 print(models)
@@ -70,6 +71,6 @@ del models
 # Now read the models back from the XML file into the container.
 # Print again the content of the container.
 # ==============================================================
-models2 = GModels(filename)
+models2 = gammalib.GModels(filename)
 print("")
 print(models2)
diff --git a/examples/python/example_radial_models.py b/examples/python/example_radial_models.py
index 4a24e91..6483b23 100755
--- a/examples/python/example_radial_models.py
+++ b/examples/python/example_radial_models.py
@@ -18,7 +18,7 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 # =====================================================================
-from gammalib import *
+import gammalib
 
 
 # ============================ #
@@ -30,11 +30,11 @@ def test_radial_model(name, model, clobber=True):
     """
     # Make a 2D FITS image (cartesian projection, celestial coordinates)
     ra, dec, binsz, npix = 0, 0, 0.01, 300
-    image = GSkymap("CAR", "CEL", ra, dec, -binsz, binsz, npix, npix, 1)
+    image = gammalib.GSkymap("CAR", "CEL", ra, dec, -binsz, binsz, npix, npix, 1)
 
     # Fill the image
-    energy = GEnergy()
-    time   = GTime()
+    energy = gammalib.GEnergy()
+    time   = gammalib.GTime()
     for inx in range(image.npix()):
         dir        = image.inx2dir(inx)
         theta      = centre.dist(dir)
@@ -62,10 +62,15 @@ if __name__ == '__main__':
     Test radial model integration.
     """
     # Set model location / centre
-    centre = GSkyDir()
+    centre = gammalib.GSkyDir()
     centre.radec_deg(0.3, 0.1)
 
     # Test models
-    test_radial_model(name='gauss', model=GModelSpatialRadialGauss(centre, 0.3))
-    test_radial_model(name='disk',  model=GModelSpatialRadialDisk(centre, 0.8))
-    test_radial_model(name='shell', model=GModelSpatialRadialShell(centre, 0.5, 0.1))
+    model = gammalib.GModelSpatialRadialGauss(centre, 0.3)
+    test_radial_model(name='gauss', model=model)
+
+    model = gammalib.GModelSpatialRadialDisk(centre, 0.8)
+    test_radial_model(name='disk',  model=model)
+
+    model = gammalib.GModelSpatialRadialShell(centre, 0.5, 0.1)
+    test_radial_model(name='shell', model=model)
diff --git a/test/test_GApplication.py b/test/test_GApplication.py
index 8e89dcc..9efecf7 100755
--- a/test/test_GApplication.py
+++ b/test/test_GApplication.py
@@ -17,10 +17,7 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 # ==========================================================================
-from gammalib import *
-from math import *
-import os
-import sys
+from gammalib import GPythonTestSuite, GLog, GApplicationPars
 
 
 # ========================================== #
diff --git a/test/test_GFits.py b/test/test_GFits.py
index 9c359bd..f09d070 100755
--- a/test/test_GFits.py
+++ b/test/test_GFits.py
@@ -17,16 +17,15 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 # ==========================================================================
-from gammalib import *
-from math import *
+import math
 import os
 import sys
-
+import gammalib
 
 # =================================== #
 # Test class for GammaLib FITS module #
 # =================================== #
-class Test(GPythonTestSuite):
+class Test(gammalib.GPythonTestSuite):
     """
     Test class for GammaLib FITS module.
     """
@@ -35,7 +34,7 @@ class Test(GPythonTestSuite):
         Constructor.
         """
         # Call base class constructor
-        GPythonTestSuite.__init__(self)
+        gammalib.GPythonTestSuite.__init__(self)
 
         # Return
         return
@@ -69,21 +68,21 @@ class Test(GPythonTestSuite):
             pass
 
         # Create FITS file
-        fits = GFits(file1, True)
+        fits = gammalib.GFits(file1, True)
         sys.stdout.write(".")
 
         # Create images
         nx = 10
         ny = 10
-        img1 = GFitsImageByte(nx, ny)
-        img2 = GFitsImageDouble(nx, ny)
-        img3 = GFitsImageFloat(nx, ny)
-        img4 = GFitsImageLong(nx, ny)
-        img5 = GFitsImageLongLong(nx, ny)
-        img6 = GFitsImageSByte(nx, ny)
-        img7 = GFitsImageShort(nx, ny)
-        img8 = GFitsImageULong(nx, ny)
-        img9 = GFitsImageUShort(nx, ny)
+        img1 = gammalib.GFitsImageByte(nx, ny)
+        img2 = gammalib.GFitsImageDouble(nx, ny)
+        img3 = gammalib.GFitsImageFloat(nx, ny)
+        img4 = gammalib.GFitsImageLong(nx, ny)
+        img5 = gammalib.GFitsImageLongLong(nx, ny)
+        img6 = gammalib.GFitsImageSByte(nx, ny)
+        img7 = gammalib.GFitsImageShort(nx, ny)
+        img8 = gammalib.GFitsImageULong(nx, ny)
+        img9 = gammalib.GFitsImageUShort(nx, ny)
         for x in range(nx):
             for y in range(ny):
                 img1[x, y] = x + y * nx
@@ -127,17 +126,17 @@ class Test(GPythonTestSuite):
 
         # Create table columns
         nrows = 10
-        col1 = GFitsTableBitCol("BIT", nrows)
-        col2 = GFitsTableBoolCol("BOOLEAN", nrows)
-        col3 = GFitsTableByteCol("BYTE", nrows)
-        col4 = GFitsTableDoubleCol("DOUBLE", nrows)
-        col5 = GFitsTableFloatCol("FLOAT", nrows)
-        col6 = GFitsTableLongCol("LONG", nrows)
-        col7 = GFitsTableLongLongCol("LONGLONG", nrows)
-        col8 = GFitsTableShortCol("SHORT", nrows)
-        col9 = GFitsTableStringCol("STRING", nrows, 20)
-        col10 = GFitsTableULongCol("ULONG", nrows)
-        col11 = GFitsTableUShortCol("USHORT", nrows)
+        col1 = gammalib.GFitsTableBitCol("BIT", nrows)
+        col2 = gammalib.GFitsTableBoolCol("BOOLEAN", nrows)
+        col3 = gammalib.GFitsTableByteCol("BYTE", nrows)
+        col4 = gammalib.GFitsTableDoubleCol("DOUBLE", nrows)
+        col5 = gammalib.GFitsTableFloatCol("FLOAT", nrows)
+        col6 = gammalib.GFitsTableLongCol("LONG", nrows)
+        col7 = gammalib.GFitsTableLongLongCol("LONGLONG", nrows)
+        col8 = gammalib.GFitsTableShortCol("SHORT", nrows)
+        col9 = gammalib.GFitsTableStringCol("STRING", nrows, 20)
+        col10 = gammalib.GFitsTableULongCol("ULONG", nrows)
+        col11 = gammalib.GFitsTableUShortCol("USHORT", nrows)
         for i in range(nrows):
             col1[i] = i % 2
             col2[i] = i % 2
@@ -153,7 +152,7 @@ class Test(GPythonTestSuite):
         sys.stdout.write(".")
 
         # Set ASCII table
-        tbl_ascii = GFitsAsciiTable()
+        tbl_ascii = gammalib.GFitsAsciiTable()
         # tbl_ascii.append(col1) # Need to implement ?/!
         # tbl_ascii.append(col2) # Need to implement ?/!
         tbl_ascii.append(col3)
@@ -170,7 +169,7 @@ class Test(GPythonTestSuite):
         sys.stdout.write(".")
 
         # Set binary table
-        tbl_bin = GFitsBinTable()
+        tbl_bin = gammalib.GFitsBinTable()
         tbl_bin.append(col1)
         tbl_bin.append(col2)
         tbl_bin.append(col3)
@@ -196,7 +195,7 @@ class Test(GPythonTestSuite):
         sys.stdout.write(".")
 
         # Re-open FITS file
-        fits = GFits(file1)
+        fits = gammalib.GFits(file1)
         sys.stdout.write(".")
 
         # Get double precision image, take square root of pixel and save in
@@ -204,7 +203,7 @@ class Test(GPythonTestSuite):
         img_double = fits.image("Double")
         for x in range(nx):
             for y in range(ny):
-                img_double[x, y] = sqrt(img_double[x, y])
+                img_double[x, y] = math.sqrt(img_double[x, y])
         sys.stdout.write(".")
 
         # Save into another FITS file
diff --git a/test/test_column_insert_remove.py b/test/test_column_insert_remove.py
index 771bd5d..b879ea1 100644
--- a/test/test_column_insert_remove.py
+++ b/test/test_column_insert_remove.py
@@ -1,13 +1,12 @@
 #! /usr/bin/env python
 
-from gammalib import *
-
+import gammalib
 
 def test_one(a):
     a[0] = 1
     a[1] = 2
     a[2] = 3
-    print a
+    print(a)
     return a
 
 
@@ -18,18 +17,18 @@ def test_two(a):
     a.remove(1, 2)
     a.remove(1, 1)
     a.remove(0, 1)
-    print a
+    print(a)
     return a
 
 
 def test_three(a):
     a.insert(0, 10)
-    print a
+    print(a)
     a[7] = 99
     a.remove(8, 2)
     a[0] = 99
     a.remove(1, 5)
-    print a
+    print(a)
     return a
 
 
@@ -38,7 +37,7 @@ def set(a):
     a[1] = 65536
     a[2] = True
     # a[3] = "Hallo"
-    print a[1]
+    print(a[1])
 
 
 #==========================#
@@ -49,17 +48,17 @@ if __name__ == '__main__':
     Perform testing.
     """
     # Perform tests
-    # a = GFitsTableBitCol("Test",4)
-    # a = GFitsTableBoolCol("Test",4)
-    # a = GFitsTableByteCol("Test",4)
-    # a = GFitsTableDoubleCol("Test",4)
-    # a = GFitsTableFloatCol("Test",4)
-    # a = GFitsTableLongCol("Test",4)
-    # a = GFitsTableLongLongCol("Test",4)
-    # a = GFitsTableShortCol("Test",4)
-    # a = GFitsTableStringCol("Test",4,10)
-    # a = GFitsTableULongCol("Test",4)
-    a = GFitsTableUShortCol("Test", 4)
+    # a = gammalib.GFitsTableBitCol("Test",4)
+    # a = gammalib.GFitsTableBoolCol("Test",4)
+    # a = gammalib.GFitsTableByteCol("Test",4)
+    # a = gammalib.GFitsTableDoubleCol("Test",4)
+    # a = gammalib.GFitsTableFloatCol("Test",4)
+    # a = gammalib.GFitsTableLongCol("Test",4)
+    # a = gammalib.GFitsTableLongLongCol("Test",4)
+    # a = gammalib.GFitsTableShortCol("Test",4)
+    # a = gammalib.GFitsTableStringCol("Test",4,10)
+    # a = gammalib.GFitsTableULongCol("Test",4)
+    a = gammalib.GFitsTableUShortCol("Test", 4)
     a = test_one(a)
     a = test_two(a)
     a = test_three(a)
diff --git a/test/test_diffuse_cube.py b/test/test_diffuse_cube.py
index 9268aab..9aa6cc6 100755
--- a/test/test_diffuse_cube.py
+++ b/test/test_diffuse_cube.py
@@ -22,9 +22,8 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 # ==========================================================================
-from gammalib import *
-from math import *
-
+import math
+from gammalib import GEnergy, GTime, GRan, GModels, GEbounds, GSkyDir
 
 # ======================== #
 # Simulate CTA observation #
@@ -94,7 +93,7 @@ def show_photons(photons, xmlname, e_min, e_max, area, duration, dir, radius, eb
             counts[index] = counts[index] + 1.0
 
         # Create error bars
-        error = [sqrt(c) for c in counts]
+        error = [math.sqrt(c) for c in counts]
 
         # Get model values
         models = GModels(xmlname)
diff --git a/test/test_file_function.py b/test/test_file_function.py
index f7efee7..ec7189c 100755
--- a/test/test_file_function.py
+++ b/test/test_file_function.py
@@ -18,9 +18,8 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 # ==========================================================================
-from gammalib import *
 import math
-
+from gammalib import GEnergy, GModelSpectralFunc, GRan, GModelSpectralPlaw
 
 # ================================= #
 # Generate file function from model #
@@ -74,7 +73,7 @@ def compare_photon_flux(model, file_function, emin=8.7, emax=1123.0, steps=13):
         flux_model = model.flux(e_min, e_max)
         flux_ff = file_function.flux(e_min, e_max)
         if (abs(flux_model - flux_ff) > eps):
-            print e_min, e_max, flux_model, flux_ff, flux_model - flux_ff
+            print(e_min, e_max, flux_model, flux_ff, flux_model - flux_ff)
             n_violate += 1
 
     # Ramp up emax from above first node
@@ -87,7 +86,7 @@ def compare_photon_flux(model, file_function, emin=8.7, emax=1123.0, steps=13):
         flux_model = model.flux(e_min, e_max)
         flux_ff = file_function.flux(e_min, e_max)
         if (abs(flux_model - flux_ff) > eps):
-            print e_min, e_max, flux_model, flux_ff, flux_model - flux_ff
+            print(e_min, e_max, flux_model, flux_ff, flux_model - flux_ff)
             n_violate += 1
 
     # Ramp down emin from above last node
@@ -99,7 +98,7 @@ def compare_photon_flux(model, file_function, emin=8.7, emax=1123.0, steps=13):
         flux_model = model.flux(e_min, e_max)
         flux_ff = file_function.flux(e_min, e_max)
         if (abs(flux_model - flux_ff) > eps):
-            print e_min, e_max, flux_model, flux_ff, flux_model - flux_ff
+            print(e_min, e_max, flux_model, flux_ff, flux_model - flux_ff)
             n_violate += 1
 
     # Ramp down emin from below last node
@@ -112,14 +111,14 @@ def compare_photon_flux(model, file_function, emin=8.7, emax=1123.0, steps=13):
         flux_model = model.flux(e_min, e_max)
         flux_ff = file_function.flux(e_min, e_max)
         if (abs(flux_model - flux_ff) > eps):
-            print e_min, e_max, flux_model, flux_ff, flux_model - flux_ff
+            print(e_min, e_max, flux_model, flux_ff, flux_model - flux_ff)
             n_violate += 1
 
     # Print summary
     if n_violate == 0:
-        print "compare_photon_flux: All tests satisfied accuracy."
+        print("compare_photon_flux: All tests satisfied accuracy.")
     else:
-        print "compare_photon_flux: " + str(n_violate) + " tests exceeded accuracy."
+        print("compare_photon_flux: " + str(n_violate) + " tests exceeded accuracy.")
 
     # Return
     return
@@ -149,7 +148,7 @@ def compare_energy_flux(model, file_function, emin=8.7, emax=1123.0, steps=13):
         flux_model = model.eflux(e_min, e_max)
         flux_ff = file_function.eflux(e_min, e_max)
         if (abs(flux_model - flux_ff) > eps):
-            print e_min, e_max, flux_model, flux_ff, flux_model - flux_ff
+            print(e_min, e_max, flux_model, flux_ff, flux_model - flux_ff)
             n_violate += 1
 
     # Ramp up emax from above first node
@@ -162,7 +161,7 @@ def compare_energy_flux(model, file_function, emin=8.7, emax=1123.0, steps=13):
         flux_model = model.eflux(e_min, e_max)
         flux_ff = file_function.eflux(e_min, e_max)
         if (abs(flux_model - flux_ff) > eps):
-            print e_min, e_max, flux_model, flux_ff, flux_model - flux_ff
+            print(e_min, e_max, flux_model, flux_ff, flux_model - flux_ff)
             n_violate += 1
 
     # Ramp down emin from above last node
@@ -174,7 +173,7 @@ def compare_energy_flux(model, file_function, emin=8.7, emax=1123.0, steps=13):
         flux_model = model.eflux(e_min, e_max)
         flux_ff = file_function.eflux(e_min, e_max)
         if (abs(flux_model - flux_ff) > eps):
-            print e_min, e_max, flux_model, flux_ff, flux_model - flux_ff
+            print(e_min, e_max, flux_model, flux_ff, flux_model - flux_ff)
             n_violate += 1
 
     # Ramp down emin from below last node
@@ -187,14 +186,14 @@ def compare_energy_flux(model, file_function, emin=8.7, emax=1123.0, steps=13):
         flux_model = model.eflux(e_min, e_max)
         flux_ff = file_function.eflux(e_min, e_max)
         if (abs(flux_model - flux_ff) > eps):
-            print e_min, e_max, flux_model, flux_ff, flux_model - flux_ff
+            print(e_min, e_max, flux_model, flux_ff, flux_model - flux_ff)
             n_violate += 1
 
     # Print summary
     if n_violate == 0:
-        print "compare_energy_flux: All tests satisfied accuracy."
+        print("compare_energy_flux: All tests satisfied accuracy.")
     else:
-        print "compare_energy_flux: " + str(n_violate) + " tests exceeded accuracy."
+        print("compare_energy_flux: " + str(n_violate) + " tests exceeded accuracy.")
 
     # Return
     return
@@ -218,7 +217,7 @@ def test_mc(file_function, trials=10):
 
 # Loop over trials
     for i in range(trials):
-        print file_function.mc(e_min, e_max, ran)
+        print(file_function.mc(e_min, e_max, ran))
 
     # Return
     return