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/test/example_model_xml_io.py b/examples/python/example_model_xml_io.py
index 01c6ac6..ef3577f 100755
--- a/test/example_model_xml_io.py
+++ b/examples/python/example_model_xml_io.py
@@ -3,7 +3,7 @@
 # This script builds a set of models, saves the models into an XML file
 # and loads the models back into memory.
 #
-# Copyright (C) 2011-2013 Juergen Knoedlseder
+# Copyright (C) 2011-2015 Juergen Knoedlseder
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -18,13 +18,13 @@
 # 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
 
 #
 # 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 +32,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, gammalib.GEnergy(100.0, "MeV"))
+point = gammalib.GModelSky(point_spatial, point_spectrum)
 point.name('My point source')
 
 #
@@ -43,9 +43,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, gammalib.GEnergy(100.0, "MeV"))
+gauss = gammalib.GModelSky(gauss_spatial, gauss_spectrum)
 gauss.name('My Gaussian source')
 
 #
@@ -53,7 +53,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 +70,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/test/example_radial_models.py b/examples/python/example_radial_models.py
index 4a24e91..6483b23 100755
--- a/test/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..6fe69b4 100755
--- a/test/test_GApplication.py
+++ b/test/test_GApplication.py
@@ -1,7 +1,7 @@
 # ==========================================================================
 # This module performs unit tests for the GammaLib application module.
 #
-# Copyright (C) 2012 Juergen Knoedlseder
+# Copyright (C) 2012-2015 Juergen Knoedlseder
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -17,16 +17,13 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 # ==========================================================================
-from gammalib import *
-from math import *
-import os
-import sys
+import gammalib
 
 
 # ========================================== #
 # Test class for GammaLib application module #
 # ========================================== #
-class Test(GPythonTestSuite):
+class Test(gammalib.GPythonTestSuite):
     """
     Test class for GammaLib application module.
     """
@@ -36,7 +33,7 @@ class Test(GPythonTestSuite):
         Constructor.
         """
         # Call base class constructor
-        GPythonTestSuite.__init__(self)
+        gammalib.GPythonTestSuite.__init__(self)
 
         # Set members
         self.logfile = "test_logger.log"
@@ -65,7 +62,7 @@ class Test(GPythonTestSuite):
         Test GLog.
         """
         # Allocate logger
-        log = GLog(self.logfile, True)
+        log = gammalib.GLog(self.logfile, True)
 
         # Test print methods
         log("Test __call__(std::string)")
@@ -161,7 +158,7 @@ class Test(GPythonTestSuite):
         # Test GApplicationPars constructor with bad filename.
         self.test_try("Test GApplicationPars constructor with bad filename")
         try:
-            pars = GApplicationPars("testme.par")
+            pars = gammalib.GApplicationPars("testme.par")
             self.test_try_failure("Exception not thrown.")
         except:
             self.test_try_success()
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..45091ed 100755
--- a/test/test_diffuse_cube.py
+++ b/test/test_diffuse_cube.py
@@ -6,7 +6,7 @@
 #
 # --------------------------------------------------------------------------
 #
-# Copyright (C) 2013 Juergen Knoedlseder
+# Copyright (C) 2013-2015 Juergen Knoedlseder
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -22,9 +22,8 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 # ==========================================================================
-from gammalib import *
-from math import *
-
+import math
+import gammalib
 
 # ======================== #
 # Simulate CTA observation #
@@ -34,20 +33,20 @@ def simulate(xmlname, e_min, e_max, area, duration, dir, radius):
     Simulate CTA observation.
     """
     # Allocate MC parameters
-    emin = GEnergy()
-    emax = GEnergy()
-    tmin = GTime(0.0)
-    tmax = GTime(duration)
+    emin = gammalib.GEnergy()
+    emax = gammalib.GEnergy()
+    tmin = gammalib.GTime(0.0)
+    tmax = gammalib.GTime(duration)
 
     # Define MC parameters
     emin.TeV(e_min)
     emax.TeV(e_max)
 
     # Allocate random number generator
-    ran = GRan()
+    ran = gammalib.GRan()
 
     # Load models and extract first model
-    models = GModels(xmlname)
+    models = gammalib.GModels(xmlname)
     model = models[0]
     #print(model)
     #print(model.spatial())
@@ -80,7 +79,8 @@ def show_photons(photons, xmlname, e_min, e_max, area, duration, dir, radius, eb
         plt.title("MC simulated photon spectrum (" + str(e_min) + '-' + str(e_max) + " TeV)")
 
         # Setup energy range covered by data
-        ebds = GEbounds(ebins, GEnergy(e_min, "TeV"), GEnergy(e_max, "TeV"))
+        ebds = gammalib.GEbounds(ebins, gammalib.GEnergy(e_min, "TeV"),
+                                        gammalib.GEnergy(e_max, "TeV"))
 
         # Create energy axis
         energy = []
@@ -94,13 +94,13 @@ 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)
+        models = gammalib.GModels(xmlname)
         mod    = models[0]
         model  = []
-        t      = GTime()
+        t      = gammalib.GTime()
         mod.spatial().set_mc_cone(dir, radius)
         for i in range(ebds.size()):
             eng    = ebds.elogmean(i)
@@ -173,7 +173,7 @@ if __name__ == '__main__':
     e_max    = 100.0              # 100 TeV
     area     = 1.0                # 1 cm2
     duration = 1.0                # 1 sec
-    dir      = GSkyDir()
+    dir      = gammalib.GSkyDir()
     dir.radec_deg(84.17263, 22.01444)
     radius   = 1.0
 
diff --git a/test/test_file_function.py b/test/test_file_function.py
index f7efee7..6b1d043 100755
--- a/test/test_file_function.py
+++ b/test/test_file_function.py
@@ -2,7 +2,7 @@
 # ==========================================================================
 # This script tests the GModelSpectralFunc file function spectral model.
 #
-# Copyright (C) 2011 Jurgen Knodlseder
+# Copyright (C) 2011-2015 Jurgen Knodlseder
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -18,9 +18,8 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 # ==========================================================================
-from gammalib import *
 import math
-
+import gammalib
 
 # ================================= #
 # Generate file function from model #
@@ -37,14 +36,14 @@ def make_file_function(model, emin=0.1e6, emax=100.0e6, nodes=10):
     file = open(filename, "w")
     for i in range(nodes):
         energy = math.pow(10.0, math.log10(emin) + i * dloge)
-        eng = GEnergy()
+        eng    = gammalib.GEnergy()
         eng.MeV(energy)
         value = model.eval(eng)
         file.write(str(energy) + " " + str(value) + "\n")
     file.close()
 
     # Generate file function
-    file_function = GModelSpectralFunc(filename)
+    file_function = gammalib.GModelSpectralFunc(filename)
 
     # Return
     return file_function
@@ -62,8 +61,8 @@ def compare_photon_flux(model, file_function, emin=8.7, emax=1123.0, steps=13):
     n_violate = 0
 
     # Allocate boundaries
-    e_min = GEnergy()
-    e_max = GEnergy()
+    e_min = gammalib.GEnergy()
+    e_max = gammalib.GEnergy()
 
     # Ramp up emax from below first node
     e_min.MeV(emin)
@@ -72,9 +71,9 @@ def compare_photon_flux(model, file_function, emin=8.7, emax=1123.0, steps=13):
         energy = math.pow(10.0, math.log10(emin) + i * dloge)
         e_max.MeV(energy)
         flux_model = model.flux(e_min, e_max)
-        flux_ff = file_function.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
@@ -137,8 +136,8 @@ def compare_energy_flux(model, file_function, emin=8.7, emax=1123.0, steps=13):
     n_violate = 0
 
     # Allocate boundaries
-    e_min = GEnergy()
-    e_max = GEnergy()
+    e_min = gammalib.GEnergy()
+    e_max = gammalib.GEnergy()
 
     # Ramp up emax from below first node
     e_min.MeV(emin)
@@ -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
@@ -208,17 +207,17 @@ def test_mc(file_function, trials=10):
     Print a bunch of Monte Carlo energies.
     """
     # Allocate boundaries
-    e_min = GEnergy()
-    e_max = GEnergy()
-    ran = GRan()
+    e_min = gammalib.GEnergy()
+    e_max = gammalib.GEnergy()
+    ran   = GRan()
 
     # Set energy range
     e_min.TeV(0.1)
     e_max.TeV(100.0)
 
-# Loop over trials
+    # 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
@@ -232,7 +231,7 @@ if __name__ == '__main__':
     Main entry point
     """
     # Allocate power law
-    model = GModelSpectralPlaw(5.7, -2.48)
+    model = gammalib.GModelSpectralPlaw(5.7, -2.48)
     model["Prefactor"].scale(1.0e-16)
     model["PivotEnergy"].value(0.3)
     model["PivotEnergy"].scale(1.0e6)