renamed from test_ to example_

knodlseder
1 parent 24a7b892
Showing 3 changed files with 88 additions and 42 deletions
test/example_model_xml_io.py
test/test_radial_models.py → test/example_radial_models.py
test/test_model_xml_io.py
+#!/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"""
+
+from gammalib import *
+
+#
+# Set the location / centre of the source model to
+# Right Ascension = 42.0 deg and Declination = 43 deg
+# ===================================================
+center = GSkyDir()
+center.radec_deg(42.0, 43.0)
+
+#
+# Create a point source with power law spectral shape
+# Normalization:   1.0e-7 ph/cm2/s/MeV @ 100 MeV
+# Spectral index: -2.1
+# ===================================================
+point_spatial  = GModelSpatialPtsrc(center)
+point_spectrum = GModelSpectralPlaw(1.0e-7, -2.1)
+point          = GModelPointSource(point_spatial, point_spectrum)
+point.name('My point source')
+
+#
+# Create a 2D Gaussian source with power law spectral shape
+# Gaussian sigma:  3.0 deg
+# Normalization:   4.2e-7 ph/cm2/s/MeV @ 100 MeV
+# Spectral index: -2.4
+# ===================================================
+gauss_spatial  = GModelRadialGauss(center, 3.0)
+gauss_spectrum = GModelSpectralPlaw(4.2e-7, -2.4)
+gauss          = GModelExtendedSource(gauss_spatial, gauss_spectrum)
+gauss.name('My Gaussian source')
+
+#
+# Allocate a GModels model container and append both source
+# models to the model container. Print the content of the
+# container
+# =========================================================
+models = GModels()
+models.append(point)
+models.append(gauss)
+print models
+
+#
+# Write the model container to XML file and delete the container
+# in-memory
+# ==============================================================
+filename = 'models.xml'
+models.save(filename)
+del models
+
+#
+# Now read the models back from the XML file into the container.
+# Print again the content of the container.
+# ==============================================================
+models2 = GModels(filename)
+print
+print models2
 #! /usr/bin/env python
+# ===========================================================================================#
+# This script makes an image from a radial model and verifies that the model
+# is correctly normalized.
+# ===========================================================================================#
 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 a 2D FITS image
-	ra, dec, binsz, npix = 0, 0, 0.01, 200
+	# 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	
 	for pix in range(image.npix()):
-		dir = image.pix2dir(pix)
-		theta = center.dist(dir)
-		image[pix, 0] = model.eval(theta);
+		dir        = image.pix2dir(pix)
+		theta      = center.dist(dir)
+		image[pix] = model.eval(theta);
  
 	# Write it to file
 	filename = name + '.fits'
@@ -20,18 +28,25 @@ def test_radial_model(name, model, clobber=True):
 	image.save(filename, clobber)
  
 	# Check that it integrates to 1
-	integral = 0
+	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
  
+# ================ #
+# Main entry point #
+# ================ #
 if __name__ == '__main__':
-	
+	"""
+	Test radial model integration.
+	"""
+	# Set model location / centre
 	center = GSkyDir()
-	center.radec_deg(0, 0)
+	center.radec_deg(0.3, 0.1)
  
+	# Test models
 	test_radial_model(name='gauss', model=GModelRadialGauss(center, 0.3))
-	test_radial_model(name='disk', model=GModelRadialDisk(center, 0.8))
+	test_radial_model(name='disk',  model=GModelRadialDisk(center, 0.8))
 	test_radial_model(name='shell', model=GModelRadialShell(center, 0.5, 0.1))
-#!/usr/bin/env python
-"""IO example:  GModels <--> XML file"""
-
-from gammalib import *
-
-# Make a GModels object containing two models
-center = GSkyDir()
-center.radec_deg(42, 43)
-
-point_spatial = GModelSpatialPtsrc(center)
-point_spectrum = GModelSpectralPlaw(1.0e-7, -2.1)
-point = GModelPointSource(point_spatial, point_spectrum)
-point.name('point')
-
-gauss_spatial = GModelRadialGauss(center, 3)
-gauss_spectrum = GModelSpectralPlaw(4.2e-7, -2.4)
-gauss = GModelExtendedSource(gauss_spatial, gauss_spectrum)
-gauss.name('gauss')
-
-models = GModels()
-models.append(point)
-models.append(gauss)
-print models
-
-# Write GModels to file and delete it in-memory
-filename = 'models.xml'
-models.save(filename)
-del models
-
-# Now read it back
-models2 = GModels(filename)
-print models2