Josh Cardenzana / gammalib

Commit ab4461b693f3ba185685cc4fe39847e4c729bbc1

Authored by Jürgen Knödlseder
1 parent c238a196

Add test script to illustrate the Psf-based computation for radial models.

Inline Side-by-side
Showing 1 changed file with 216 additions and 0 deletions
inst/cta/test/test_psf_radial.py 0 → 100755
  View file @ab4461b
  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, theta_max, obsDir, psf_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 + delta_mod = obsDir.dist_deg(centre)
  95 + posangle_mod = obsDir.posang(centre)
  96 +
  97 + # Delta limits
  98 + delta_min = 0.0
  99 + delta_max = delta_mod + theta_max
  100 + if (delta_mod > theta_max):
  101 + delta_min = delta_mod - theta_max
  102 + if delta_max > psf_max:
  103 + delta_max = psf_max
  104 +
  105 + # Convert to radians
  106 + delta_mod *= gammalib.deg2rad
  107 + theta_max *= gammalib.deg2rad
  108 +
  109 + # Setup deltas
  110 + deltas = []
  111 + for i in range(51):
  112 + delta = i*0.02
  113 + if (delta >= delta_min and delta <= delta_max):
  114 + deltas.append(delta)
  115 +
  116 + # Plot Psf validity circle
  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(psf_max * cos_radians)
  124 + y.append(psf_max * sin_radians)
  125 + plt.plot(x, y, 'b-')
  126 + plt.xlim([ 2.0*psf_max,-2.0*psf_max])
  127 + plt.ylim([-2.0*psf_max, 2.0*psf_max])
  128 + plt.plot([0.0], [0.0], 'bo')
  129 +
  130 + # Plot model direction
  131 + x_mod = centre.ra_deg() - obsDir.ra_deg()
  132 + y_mod = centre.dec_deg() - obsDir.dec_deg()
  133 + plt.plot([x_mod], [y_mod], 'ro')
  134 +
  135 + # Plot radial model
  136 + x = []
  137 + y = []
  138 + for angle in range(0, 361):
  139 + radians = angle * gammalib.deg2rad
  140 + cos_radians = math.cos(radians)
  141 + sin_radians = math.sin(radians)
  142 + x.append(theta_max * cos_radians * gammalib.rad2deg + x_mod)
  143 + y.append(theta_max * sin_radians * gammalib.rad2deg + y_mod)
  144 + plt.plot(x, y, 'r-')
  145 +
  146 + # Plot deltas
  147 + for delta in deltas:
  148 +
  149 + # Compute some stuff
  150 + delta_rad = delta * gammalib.deg2rad
  151 +
  152 + # Compute half length of arc that lies within model
  153 + dphi = 0.5 * cta_roi_arclength(delta_rad, delta_mod, theta_max)
  154 +
  155 + # Continue only if arc length is positive
  156 + if (dphi > 0.0):
  157 +
  158 + # Set Phi range
  159 + phi_min = -dphi
  160 + phi_max = +dphi
  161 +
  162 + # Plot segment
  163 + n = 100
  164 + x = []
  165 + y = []
  166 + step = (phi_max - phi_min)/n
  167 + for i in range(n+1):
  168 + radians = i*step + phi_min
  169 + cos_radians = math.sin(radians+posangle_mod)
  170 + sin_radians = math.cos(radians+posangle_mod)
  171 + x.append(delta * cos_radians)
  172 + y.append(delta * sin_radians)
  173 + plt.plot(x, y, 'k-')
  174 +
  175 + # Show plot
  176 + plt.show()
  177 +
  178 + # Return
  179 + return
  180 +
  181 +
  182 +# ======================== #
  183 +# Main routine entry point #
  184 +# ======================== #
  185 +if __name__ == '__main__':
  186 + """
  187 + Test radial response computation
  188 + """
  189 + # Dump header
  190 + print("")
  191 + print("*******************************")
  192 + print("* Test radial Psf computation *")
  193 + print("*******************************")
  194 +
  195 + # Check if matplotlib is available
  196 + try:
  197 + import matplotlib.pyplot as plt
  198 + except ImportError:
  199 + print("Matplotlib is not installed on your system. Abort.")
  200 + sys.exit()
  201 +
  202 + # Set Psf size
  203 + psf_max = 1.0
  204 +
  205 + # Set observed photon direction
  206 + obsDir = gammalib.GSkyDir()
  207 + obsDir.radec_deg(0.0, 0.0)
  208 +
  209 + # Set model parameters
  210 + centre = gammalib.GSkyDir()
  211 + centre.radec_deg(0.0, 0.6)
  212 + theta_max = 0.5
  213 +
  214 + # Show elliptical response geometry
  215 + show_response(centre, theta_max, obsDir, psf_max)
  216 +
... ...