Leonardo Di Venere / gammalib

Commit 813fe75129fec24a971d62c52eaa6d5d61c8d9db

Authored by Jürgen Knödlseder
1 parent e9e98866

Add test script for radial response computation

Inline Side-by-side
Showing 1 changed file with 228 additions and 0 deletions
inst/cta/test/test_irf_radial.py 0 → 100755
  View file @813fe75
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This script tests the radial response computation. The script requires
  4 +# matplotlib to be installed.
  5 +#
  6 +# --------------------------------------------------------------------------
  7 +#
  8 +# Copyright (C) 2014 Juergen Knoedlseder
  9 +#
  10 +# This program is free software: you can redistribute it and/or modify
  11 +# it under the terms of the GNU General Public License as published by
  12 +# the Free Software Foundation, either version 3 of the License, or
  13 +# (at your option) any later version.
  14 +#
  15 +# This program is distributed in the hope that it will be useful,
  16 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  17 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  18 +# GNU General Public License for more details.
  19 +#
  20 +# You should have received a copy of the GNU General Public License
  21 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  22 +#
  23 +# ==========================================================================
  24 +import gammalib
  25 +import math
  26 +
  27 +
  28 +# ================= #
  29 +# Compute arclength #
  30 +# ================= #
  31 +def cta_roi_arclength(rad, dist, roi):
  32 + """
  33 + Returns length of circular arc within circular ROI
  34 +
  35 + Parameters:
  36 + rad - Circle radius in radians (<pi)
  37 + dist - Circle centre distance to ROI centre (<pi)
  38 + roi - Radius of ROI in radians
  39 + """
  40 + # Initialise
  41 + arclength = 0.0
  42 +
  43 + # Handle special case that circle centre matches ROI centre
  44 + if dist == 0.0:
  45 + if rad > roi:
  46 + arclength = 0.0
  47 + else:
  48 + arclength = 2.0 * math.pi
  49 +
  50 + # ... otherwise circle and ROI centres are not identical
  51 + else:
  52 +
  53 + # Handle special case that we evaluate exactly at the circle
  54 + # centre. In this case we have in fact a point, and if this point
  55 + # falls within the ROI it has a formal arclength of 2pi.
  56 + if rad == 0.0:
  57 + if dist > roi:
  58 + arclength = 0.0
  59 + else:
  60 + arclength = 2.0 * math.pi
  61 +
  62 + # ... otherwise we have to handle the general case
  63 + else:
  64 + d = roi - dist
  65 + if -rad >= d:
  66 + arclength = 0.0
  67 + elif rad <= d:
  68 + arclength = 2.0 * math.pi
  69 + else:
  70 + cosrad = math.cos(rad)
  71 + sinrad = math.sin(rad)
  72 + cosroi = math.cos(roi)
  73 + cosdist = math.cos(dist)
  74 + sindist = math.sin(dist)
  75 + cosang = (cosroi - cosdist * cosrad) / (sindist * sinrad)
  76 + arclength = 2.0 * gammalib.acos(cosang)
  77 +
  78 + # Return
  79 + return arclength
  80 +
  81 +
  82 +# ==================== #
  83 +# Show radial response #
  84 +# ==================== #
  85 +def show_response(centre, radius, obsDir, delta_max):
  86 + """
  87 + Show radial response
  88 + """
  89 + # Create figure
  90 + plt.figure(1)
  91 + plt.title("Radial response")
  92 +
  93 + # Compute some stuff
  94 + rho_obs = centre.dist_deg(obsDir)
  95 + posangle_obs = centre.posang(obsDir)
  96 +
  97 + # Rho limits
  98 + rho_min = 0.0
  99 + rho_max = rho_obs + delta_max
  100 + if (rho_obs > delta_max):
  101 + rho_min = rho_obs - delta_max
  102 + if rho_max > radius:
  103 + rho_max = radius
  104 +
  105 + # Convert to radians
  106 + delta_max *= gammalib.deg2rad
  107 + rho_obs *= gammalib.deg2rad
  108 +
  109 + # Setup rhos
  110 + rhos = []
  111 + for i in range(51):
  112 + rho = i*0.02
  113 + if (rho >= rho_min and rho <= rho_max):
  114 + rhos.append(rho)
  115 +
  116 + # Plot radial model
  117 + x = []
  118 + y = []
  119 + for angle in range(0, 361):
  120 + radians = angle * gammalib.deg2rad
  121 + cos_radians = math.cos(radians)
  122 + sin_radians = math.sin(radians)
  123 + x.append(radius * cos_radians)
  124 + y.append(radius * sin_radians)
  125 + plt.plot(x, y, 'r-')
  126 + plt.xlim([ 2.0*radius,-2.0*radius])
  127 + plt.ylim([-2.0*radius, 2.0*radius])
  128 + plt.plot([0.0], [0.0], 'ro')
  129 +
  130 + # Plot observed direction
  131 + x_obs = obsDir.ra_deg() - centre.ra_deg()
  132 + y_obs = obsDir.dec_deg() - centre.dec_deg()
  133 + plt.plot([x_obs], [y_obs], 'bo')
  134 +
  135 + # Draw connecting line model-photon'
  136 + n = 10
  137 + step = rho_obs/n * gammalib.rad2deg
  138 + x = []
  139 + y = []
  140 + for i in range(n+1):
  141 + x.append(i*step * math.sin(posangle_obs))
  142 + y.append(i*step * math.cos(posangle_obs))
  143 + plt.plot(x, y, 'b-')
  144 +
  145 + # Plot Psf validity circle
  146 + x = []
  147 + y = []
  148 + for angle in range(0, 361):
  149 + radians = angle * gammalib.deg2rad
  150 + cos_radians = math.cos(radians)
  151 + sin_radians = math.sin(radians)
  152 + x.append(delta_max * cos_radians * gammalib.rad2deg + x_obs)
  153 + y.append(delta_max * sin_radians * gammalib.rad2deg + y_obs)
  154 + plt.plot(x, y, 'b-')
  155 +
  156 + # Plot rhos
  157 + for rho in rhos:
  158 +
  159 + # Compute some stuff
  160 + rho_rad = rho * gammalib.deg2rad
  161 +
  162 + # Compute half length of arc that lies within Psf validity circle
  163 + domega = 0.5 * cta_roi_arclength(rho_rad, rho_obs, delta_max)
  164 +
  165 + # Continue only if arc length is positive
  166 + if (domega > 0.0):
  167 +
  168 + # Set Omega range
  169 + omega_min = -domega
  170 + omega_max = +domega
  171 +
  172 + # Plot segment
  173 + n = 100
  174 + x = []
  175 + y = []
  176 + step = (omega_max - omega_min)/n
  177 + for i in range(n+1):
  178 + radians = i*step + omega_min
  179 + cos_radians = math.sin(radians+posangle_obs)
  180 + sin_radians = math.cos(radians+posangle_obs)
  181 + x.append(rho * cos_radians)
  182 + y.append(rho * sin_radians)
  183 + plt.plot(x, y, 'k-')
  184 +
  185 + # Show plot
  186 + plt.show()
  187 +
  188 + # Return
  189 + return
  190 +
  191 +
  192 +# ======================== #
  193 +# Main routine entry point #
  194 +# ======================== #
  195 +if __name__ == '__main__':
  196 + """
  197 + Test radial response computation
  198 + """
  199 + # Dump header
  200 + print("")
  201 + print("************************************")
  202 + print("* Test radial response computation *")
  203 + print("************************************")
  204 +
  205 + # Check if matplotlib is available
  206 + try:
  207 + import matplotlib.pyplot as plt
  208 + except ImportError:
  209 + print("Matplotlib is not installed on your system. Abort.")
  210 + sys.exit()
  211 +
  212 + # Set Psf size
  213 + delta_max = 0.6
  214 +
  215 + # Set model parameters
  216 + centre = gammalib.GSkyDir()
  217 + centre.radec_deg(0.0, 0.0)
  218 + radius = 0.8
  219 +
  220 + # Set observed photon direction
  221 + obsDir = gammalib.GSkyDir()
  222 + obsDir.radec_deg(0.4, -0.3)
  223 + #obsDir.radec_deg(-0.4, -0.3)
  224 + #obsDir.radec_deg(0.4, 0.4)
  225 +
  226 + # Show elliptical response geometry
  227 + show_response(centre, radius, obsDir, delta_max)
  228 +
... ...