Domenico Tiziani / gammalib

Commit f8f5bf0c86be20b332f10b408aa893e34348dca3

Authored by Jürgen Knödlseder
1 parent 86da4f6c

Add shell model visualization.

Inline Side-by-side
Showing 1 changed file with 132 additions and 0 deletions
inst/cta/test/test_shell.py 0 → 100755
  View file @f8f5bf0
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This script shows the shell model radial distribution.
  4 +#
  5 +# Copyright (C) 2014 Juergen Knoedlseder
  6 +#
  7 +# This program is free software: you can redistribute it and/or modify
  8 +# it under the terms of the GNU General Public License as published by
  9 +# the Free Software Foundation, either version 3 of the License, or
  10 +# (at your option) any later version.
  11 +#
  12 +# This program is distributed in the hope that it will be useful,
  13 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  15 +# GNU General Public License for more details.
  16 +#
  17 +# You should have received a copy of the GNU General Public License
  18 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  19 +#
  20 +# ==========================================================================
  21 +import gammalib
  22 +import math
  23 +# Check if matplotlib is available
  24 +try:
  25 + import matplotlib.pyplot as plt
  26 + has_matplotlib = True
  27 +except ImportError:
  28 + print("Matplotlib is not (correctly) installed on your system.")
  29 + has_matplotlib = False
  30 +
  31 +
  32 +# ================ #
  33 +# Compute function #
  34 +# ================ #
  35 +def function(theta, radius=0.30, width=0.10):
  36 + """
  37 + """
  38 + # Initialise
  39 + value = 0.0
  40 + sin_theta = math.sin(theta*gammalib.deg2rad)
  41 + sin_theta2 = sin_theta*sin_theta
  42 + sin_theta_in = math.sin(radius*gammalib.deg2rad)
  43 + sin_theta2_in = sin_theta_in*sin_theta_in
  44 + sin_theta_out = math.sin((radius+width)*gammalib.deg2rad)
  45 + sin_theta2_out = sin_theta_out*sin_theta_out
  46 +
  47 + # Compute function value
  48 + if (sin_theta2 < sin_theta2_out):
  49 + value = math.sqrt(sin_theta2_out - sin_theta2)
  50 + if (sin_theta2 < sin_theta2_in):
  51 + value -= math.sqrt(sin_theta2_in - sin_theta2)
  52 +
  53 + # Return value
  54 + return value
  55 +
  56 +
  57 +# ============= #
  58 +# Show function #
  59 +# ============= #
  60 +def show(theta_max=0.5, theta_bins=100):
  61 + """
  62 + """
  63 +
  64 + # Continue only if matplotlib is available
  65 + if has_matplotlib:
  66 +
  67 + # Create figure
  68 + plt.figure(1)
  69 + plt.title("Shell function")
  70 +
  71 + # Create theta axis
  72 + theta = []
  73 + dtheta = theta_max/theta_bins
  74 + for i in range(theta_bins):
  75 + theta.append((i)*dtheta)
  76 +
  77 + # Get function
  78 + value = []
  79 + for i in range(theta_bins):
  80 + value.append(function(theta[i]))
  81 +
  82 + # Plot function
  83 + plt.plot(theta, value, 'r-')
  84 +
  85 + # Set axes
  86 + plt.xlabel("Theta angle (degrees)")
  87 + plt.ylabel("Function")
  88 +
  89 + # Create figure
  90 + plt.figure(2)
  91 + plt.title("MC function")
  92 +
  93 + # Get functions
  94 + value = []
  95 + envelope = []
  96 + for i in range(theta_bins):
  97 + value.append(function(theta[i])*math.sin(theta[i]*gammalib.deg2rad))
  98 + envelope.append(math.sin(0.4*gammalib.deg2rad)*math.sin(0.4*gammalib.deg2rad))
  99 +
  100 + # Plot function
  101 + plt.plot(theta, value, 'r-')
  102 + plt.plot(theta, envelope, 'b-')
  103 +
  104 + # Set axes
  105 + plt.xlabel("Theta angle (degrees)")
  106 + plt.ylabel("Function*sin(theta)")
  107 +
  108 + # Notify
  109 + print("PLEASE CLOSE WINDOW TO CONTINUE ...")
  110 +
  111 + # Show plot
  112 + plt.show()
  113 +
  114 + # Return
  115 + return
  116 +
  117 +
  118 +#==========================#
  119 +# Main routine entry point #
  120 +#==========================#
  121 +if __name__ == '__main__':
  122 + """
  123 + Show Shell function profile
  124 + """
  125 + # Dump header
  126 + print("")
  127 + print("*******************************")
  128 + print("* Show Shell function profile *")
  129 + print("*******************************")
  130 +
  131 + # Show
  132 + show()
... ...