Michael Mayer / ctools

Commit d8896b9b0ac75955cef5a1c1c04edd18d36ebfe9

Authored by Jürgen Knödlseder
1 parent b2ed3f25

Create initial version of make_pointing.py script for GPS KSP.

Inline Side-by-side
Showing 1 changed file with 164 additions and 0 deletions
examples/make_pointings.py 0 → 100755
  View file @d8896b9
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This script simulates pointing patterns for CTA surveys.
  4 +#
  5 +# Usage:
  6 +# ./make_pointings.py
  7 +#
  8 +# ==========================================================================
  9 +#
  10 +# Copyright (C) 2015 Juergen Knoedlseder
  11 +#
  12 +# This program is free software: you can redistribute it and/or modify
  13 +# it under the terms of the GNU General Public License as published by
  14 +# the Free Software Foundation, either version 3 of the License, or
  15 +# (at your option) any later version.
  16 +#
  17 +# This program is distributed in the hope that it will be useful,
  18 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  19 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  20 +# GNU General Public License for more details.
  21 +#
  22 +# You should have received a copy of the GNU General Public License
  23 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  24 +#
  25 +# ==========================================================================
  26 +import gammalib
  27 +import sys
  28 +import math
  29 +
  30 +
  31 +# ======================= #
  32 +# Setup double row scheme #
  33 +# ======================= #
  34 +def set_double_row(lmin=-30.0, lmax=30.0, separation=3.0, hours=100.0):
  35 + """
  36 + Setup double row scheme.
  37 + """
  38 + # Initialise observation definition
  39 + obsdef = []
  40 +
  41 + # Compute number of pointings in row
  42 + n_pnt_per_row = int((lmax-lmin)/separation+0.5)
  43 + n_pnt = 2*n_pnt_per_row
  44 +
  45 + # Compute offset from Galactic plane (degrees)
  46 + offset = separation * 0.5 * math.sqrt(0.75)
  47 +
  48 + # Compute duration per pointing (seconds)
  49 + duration = hours*3600.0/float(n_pnt)
  50 +
  51 + # Set upper row
  52 + lon = lmin
  53 + for i in range(n_pnt_per_row):
  54 + obs = {'lon': lon, 'lat': offset, 'duration': duration}
  55 + obsdef.append(obs)
  56 + lon += separation
  57 +
  58 + # Set lower row
  59 + lon = lmin + 0.5*separation
  60 + for i in range(n_pnt_per_row):
  61 + obs = {'lon': lon, 'lat': -offset, 'duration': duration}
  62 + obsdef.append(obs)
  63 + lon += separation
  64 +
  65 + # Return observation definition
  66 + return obsdef
  67 +
  68 +
  69 +# =========================== #
  70 +# Setup Galactic plane survey #
  71 +# =========================== #
  72 +def set_gps():
  73 + """
  74 + Setup Galactic plane survey.
  75 + """
  76 + # Initialise observation definition
  77 + obsdef = []
  78 +
  79 + # Add inner region South
  80 + obsdef.extend(set_double_row(lmin=-60.0, lmax=60.0, hours=780))
  81 +
  82 + # Add Vela & Carina region
  83 + obsdef.extend(set_double_row(lmin=240.0, lmax=300.0, hours=180))
  84 +
  85 + # Add 210-240 region
  86 + obsdef.extend(set_double_row(lmin=210.0, lmax=240.0, hours=60))
  87 +
  88 + # Add Cygnus, Perseus
  89 + obsdef.extend(set_double_row(lmin=60.0, lmax=150.0, hours=450))
  90 +
  91 + # Add Anticentre
  92 + obsdef.extend(set_double_row(lmin=150.0, lmax=210.0, hours=150))
  93 +
  94 + # Return observation definition
  95 + return obsdef
  96 +
  97 +
  98 +# ======================================= #
  99 +# Write observation definition dictionary #
  100 +# ======================================= #
  101 +def write_obsdef(filename, obsdef):
  102 + """
  103 + Write observation definition file.
  104 + """
  105 + # Open file
  106 + file = open(filename, 'w')
  107 +
  108 + # Write header
  109 + file.write("ra,dec,duration\n")
  110 +
  111 + # Loop over pointings
  112 + for obs in obsdef:
  113 +
  114 + # If we have lon,lat then convert into RA,Dec
  115 + if "lon" in obs and "lat" in obs:
  116 + lon = obs["lon"]
  117 + lat = obs["lat"]
  118 + dir = gammalib.GSkyDir()
  119 + dir.lb_deg(lon,lat)
  120 + ra = dir.ra_deg()
  121 + dec = dir.dec_deg()
  122 + else:
  123 + ra = obs["ra"]
  124 + dec = obs["dec"]
  125 +
  126 + # Write information
  127 + file.write("%8.4f,%8.4f,%.4f\n" % (ra, dec, obs["duration"]))
  128 +
  129 + # Close file
  130 + file.close()
  131 +
  132 + # Return
  133 + return
  134 +
  135 +
  136 +# ======================== #
  137 +# Main routine entry point #
  138 +# ======================== #
  139 +if __name__ == '__main__':
  140 + """
  141 + CTA survey simulation.
  142 + """
  143 + # Initialise flags
  144 + need_help = False
  145 +
  146 + # Test for command line arguments
  147 + print(sys.argv[0])
  148 + if (len(sys.argv) > 1):
  149 + if sys.argv[1] == "-h":
  150 + need_help = True
  151 + else:
  152 + need_help = True
  153 +
  154 + # Print help if needed and exit
  155 + if need_help:
  156 + print("Usage: make_pointing.py [OPTIONS]")
  157 + print(" -h Display this usage message")
  158 + sys.exit()
  159 +
  160 + # Test
  161 + obsdef = set_gps()
  162 +
  163 + # Write observation definition file
  164 + write_obsdef("test.dat", obsdef)
... ...