diff --git a/inst/cta/test/test_gauss.py b/inst/cta/test/test_gauss.py
index 6bea8e2..bde91d3 100755
--- a/inst/cta/test/test_gauss.py
+++ b/inst/cta/test/test_gauss.py
@@ -35,11 +35,13 @@ def gauss(dir, sigma):
     # Compute derivative
     dh = 0.0001
     s  = GSkyDir()
+    e  = GEnergy()
+    t  = GTime()
     f1 = GModelSpatialRadialGauss(s, sigma)
     f2 = GModelSpatialRadialGauss(s, sigma + dh)
-    v1 = f1.eval_gradients(dir)
+    v1 = f1.eval_gradients(GPhoton(dir, e, t))
     g1 = f1[2].gradient()
-    v2 = f2.eval_gradients(dir)
+    v2 = f2.eval_gradients(GPhoton(dir, e, t))
     g2 = f2[2].gradient()
     g  = (v2 - v1) / dh
 
@@ -79,10 +81,12 @@ def show_gaussian(sigma):
         f_expected = []
         sigma_rad  = sigma * (pi / 180.0)
         norm       = 1.0 / (2.0 * pi * sigma_rad * sigma_rad)
+        eng        = GEnergy()
+        time       = GTime()
         for t in theta:
             s = GSkyDir()
             s.radec_deg(0.0, t)
-            f = gauss.eval(s)
+            f = gauss.eval(GPhoton(s, eng, time))
             e = norm * exp(-0.5 * t * t / sigma / sigma)
             f_gauss.append(f)
             f_expected.append(e)
diff --git a/inst/cta/test/test_model.py b/inst/cta/test/test_model.py
index 0ef4cf6..76e5f73 100755
--- a/inst/cta/test/test_model.py
+++ b/inst/cta/test/test_model.py
@@ -36,7 +36,7 @@ def ptsrc_model(ra=0.0, dec=0.0):
     spatial = GModelSpatialPointSource(pos)
 
     # Set spectral model
-    spectral = GModelSpectralPlaw(1.0, -2.0, 100.0)
+    spectral = GModelSpectralPlaw(1.0, -2.0, GEnergy(100.0, "MeV"))
 
     # Set sky model
     model = GModelSky(spatial, spectral)
@@ -60,7 +60,7 @@ def shell_model(ra=0.3, dec=0.3, radius=0.3, width=0.1):
     radial = GModelSpatialRadialShell(center, radius, width, False)
 
     # Set spectral model
-    spectral = GModelSpectralPlaw(1.0, -2.0, 100.0)
+    spectral = GModelSpectralPlaw(1.0, -2.0, GEnergy(100.0, "MeV"))
 
     # Set sky model
     model = GModelSky(radial, spectral)
@@ -84,7 +84,7 @@ def disk_model(ra=359.6, dec=-0.2, radius=0.4):
     radial = GModelSpatialRadialDisk(center, radius)
 
     # Set spectral model
-    spectral = GModelSpectralPlaw(1.0, -2.0, 100.0)
+    spectral = GModelSpectralPlaw(1.0, -2.0, GEnergy(100.0, "MeV"))
 
     # Set sky model
     model = GModelSky(radial, spectral)
@@ -108,7 +108,7 @@ def gauss_model(ra=359.6, dec=+0.1, sigma=0.2):
     radial = GModelSpatialRadialGauss(center, sigma)
 
     # Set spectral model
-    spectral = GModelSpectralPlaw(1.0, -2.0, 100.0)
+    spectral = GModelSpectralPlaw(1.0, -2.0, GEnergy(100.0, "MeV"))
 
     # Set sky model
     model = GModelSky(radial, spectral)
@@ -128,7 +128,7 @@ def diffuse_model(file="../../../test/data/cena_lobes_parkes.fits"):
     spatial = GModelSpatialDiffuseMap(file)
 
     # Set spectral model
-    spectral = GModelSpectralPlaw(1.0, -2.0, 100.0)
+    spectral = GModelSpectralPlaw(1.0, -2.0, GEnergy(100.0, "MeV"))
 
     # Set sky model
     model = GModelSky(spatial, spectral)
diff --git a/inst/mwl/test/test_python.py b/inst/mwl/test/test_python.py
index aab1509..8df5190 100755
--- a/inst/mwl/test/test_python.py
+++ b/inst/mwl/test/test_python.py
@@ -93,7 +93,7 @@ def fit_spectrum(filename, xmlname):
         yval = 0.0
         for model in obs.models():
             if model.type() == "PointSource":
-                yval += model.value(d, e, t) * energy * energy * 1.6021765e-6
+                yval += model.value(GPhoton(d, e, t)) * energy * energy * 1.6021765e-6
         y.append(yval)
     plt.plot(x, y)
 
diff --git a/test/example_model_xml_io.py b/test/example_model_xml_io.py
index 6fe4ab7..cd202aa 100755
--- a/test/example_model_xml_io.py
+++ b/test/example_model_xml_io.py
@@ -1,11 +1,23 @@
 #!/usr/bin/env python
-# ==========================================================================
-# This script builds a set of models, saves the models into an XML file and
-# loads the models back into memory.
-# ==========================================================================
-
-"""IO example:  GModels <--> XML file"""
-
+# =====================================================================
+# 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
+#
+# 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
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# 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 *
 
 #
@@ -21,7 +33,7 @@ center.radec_deg(42.0, 43.0)
 # Spectral index: -2.1
 # ===================================================
 point_spatial  = GModelSpatialPointSource(center)
-point_spectrum = GModelSpectralPlaw(1.0e-7, -2.1, 100.0)
+point_spectrum = GModelSpectralPlaw(1.0e-7, -2.1, GEnergy(100.0, "MeV"))
 point = GModelSky(point_spatial, point_spectrum)
 point.name('My point source')
 
@@ -32,7 +44,7 @@ point.name('My point source')
 # Spectral index: -2.4
 # ===================================================
 gauss_spatial  = GModelSpatialRadialGauss(center, 3.0)
-gauss_spectrum = GModelSpectralPlaw(4.2e-7, -2.4, 100.0)
+gauss_spectrum = GModelSpectralPlaw(4.2e-7, -2.4, GEnergy(100.0, "MeV"))
 gauss = GModelSky(gauss_spatial, gauss_spectrum)
 gauss.name('My Gaussian source')
 
diff --git a/test/example_radial_models.py b/test/example_radial_models.py
index 0a30d4c..d12c605 100755
--- a/test/example_radial_models.py
+++ b/test/example_radial_models.py
@@ -1,8 +1,23 @@
 #! /usr/bin/env python
-# ==========================================================================
-# This script makes an image from a radial model and verifies that the model
-# is correctly normalized.
-# ==========================================================================
+# =====================================================================
+# This script makes an image from a radial model and verifies that the
+# model is correctly normalized.
+#
+# Copyright (C) 2011-2013 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
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# 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 *
 
 
@@ -10,30 +25,34 @@ from gammalib import *
 # Radial model testing routine #
 # ============================ #
 def test_radial_model(name, model, clobber=True):
-    """Make an image of a spatial model
-    and check that it integrates to 1"""
+    """
+    Make an image of a spatial model and check that it integrates to 1.
+    """
     # 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)
 
     # Fill the image
+    energy = GEnergy()
+    time   = GTime()
     for pix in range(image.npix()):
-        dir = image.pix2dir(pix)
-        theta = center.dist(dir)
-        image[pix] = model.eval(theta)
+        dir        = image.pix2dir(pix)
+        theta      = center.dist(dir)
+        image[pix] = model.eval(theta, energy, time)
 
     # Write it to file
     filename = name + '.fits'
-    print 'Writing', filename
+    print('Writing', filename)
     image.save(filename, clobber)
 
     # Check that it integrates to 1
     integral = 0.0
     for pix in range(image.npix()):
         integral += image.omega(pix) * image[pix]
-    print 'integral = %g (should be 1)' % integral
-    print 'integral error = %g (should be 0)' % (integral - 1)
-    print
+    print('integral = %g (should be 1)' % integral)
+    print('integral error = %g (should be 0)' % (integral - 1))
+    print('')
+
 
 # ================ #
 # Main entry point #
diff --git a/testall.sh b/testall.sh
index b384b33..caf9ea1 100755
--- a/testall.sh
+++ b/testall.sh
@@ -2,7 +2,7 @@
 #
 # Run all test and example scripts that come with the package
 #
-# Copyright (C) 2011 Juergen Knoedlseder
+# Copyright (C) 2011-2013 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
@@ -20,6 +20,17 @@
 base=$PWD
 
 
+# examples
+# ========
+echo
+echo "=====> examples"
+cd examples/python
+./models_howto.py
+./xml_howto.py
+./xml_html_create.py
+cd $base
+
+
 # test
 # ====
 echo