Pierrick Martin / ctools

Commit 844130f6ff66b7335f01dec9e40ca78d75c2e917

Authored by Jürgen Knödlseder
1 parent 2469a2e5

Add GC survey

Inline Side-by-side
Showing 1 changed file with 204 additions and 47 deletions
examples/make_pointings.py
  View file @844130f
1 1 #! /usr/bin/env python
2 2 # ==========================================================================
3   -# This script simulates pointing patterns for CTA surveys.
  3 +# This script generated pointing patterns for CTA observations.
4 4 #
5 5 # Usage:
6 6 # ./make_pointings.py
... ... @@ -28,52 +28,133 @@ import sys
28 28 import math
29 29  
30 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, site=None):
  31 +# ============================================================ #
  32 +# Get positions of pointings on sky for an interleaved pattern #
  33 +# ============================================================ #
  34 +def get_positions(xmin, xmax, ymin, ymax, step):
35 35 """
36   - Setup double row scheme.
  36 + Get positions for a patch of interleaved pointings on the sky.
  37 + The function does respect the latitude dependence of the
  38 + pointing distance. The step parameter is for a latitude of
  39 + zero.
  40 + """
  41 + # Initialise positions
  42 + positions = []
  43 +
  44 + # Determine ystep and number of rows
  45 + ystep = step * math.sqrt(0.75)
  46 + ny = int((ymax-ymin)/ystep+1.5)
  47 + rescale = (ymax-ymin)/(ystep*(ny-1))
  48 + ystep *= rescale
  49 + print("Number of rows ....: %d" % (ny))
  50 +
  51 + # Loop over rows
  52 + y = ymin
  53 + for row in range(ny):
  54 +
  55 + # Determine xstep and number of pointings in row
  56 + xstep = step * math.cos(gammalib.deg2rad * y)
  57 + nx = int((xmax-xmin)/xstep+1.5)
  58 + rescale = (xmax-xmin)/(xstep*(nx-1))
  59 + xstep *= rescale
  60 +
  61 + # Set x offset
  62 + if row % 2 == 0:
  63 + x = xmin
  64 + else:
  65 + x = xmin + 0.5 * xstep
  66 + #print(row, x)
  67 +
  68 + # Append pointings
  69 + for pnt in range(nx):
  70 + positions.append({'x': x, 'y': y})
  71 + x += xstep
  72 + if x >= 360.0:
  73 + x -= 360.0
  74 +
  75 + # Increment y
  76 + y += ystep
  77 +
  78 + # Return positions
  79 + return positions
  80 +
  81 +
  82 +# ====================== #
  83 +# Setup patch on the sky #
  84 +# ====================== #
  85 +def set_patch(lmin=-30.0, lmax=30.0, bmin=-1.5, bmax=+1.5, \
  86 + separation=3.0, hours=100.0, \
  87 + site=None, lst=True, autodec=0.0):
  88 + """
  89 + Setup pointing patch on the sky.
37 90 """
38 91 # Initialise observation definition
39 92 obsdef = []
40 93  
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
  94 + # Get pointing positions
  95 + pointings = get_positions(lmin, lmax, bmin, bmax, separation)
44 96  
45   - # Compute offset from Galactic plane (degrees)
46   - offset = separation * 0.5 * math.sqrt(0.75)
47   -
48   - # Compute duration per pointing (seconds)
  97 + # Determine number of pointings and pointing duration
  98 + n_pnt = len(pointings)
49 99 duration = hours*3600.0/float(n_pnt)
50 100  
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
  101 + # Observing time collectors
  102 + exposure_south = 0.0
  103 + exposure_north = 0.0
57 104  
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
  105 + # Set observations
  106 + for pnt in pointings:
64 107  
65   - # Optionally add-in site dependent IRF
66   - if site != None:
67   - if site == "North":
68   - caldb = "tenerife"
69   - irf = "DESY20140105_50h"
70   - else:
71   - caldb = "aar"
72   - irf = "DESY20140105_50h"
73   - for obs in obsdef:
  108 + # Set positions and duration
  109 + obs = {'lon': pnt['x'], 'lat': pnt['y'], 'duration': duration}
  110 +
  111 + # Optionally add-in site dependent IRF
  112 + if site != None:
  113 +
  114 + # Handle automatic site switching
  115 + if site == "Automatic":
  116 + pos = gammalib.GSkyDir()
  117 + pos.lb_deg(pnt['x'], pnt['y'])
  118 + dec = pos.dec_deg()
  119 + if (dec >= autodec):
  120 + obs_site = "North"
  121 + else:
  122 + obs_site = "South"
  123 + else:
  124 + obs_site = site
  125 +
  126 + # Set site IRF
  127 + if obs_site == "North":
  128 + caldb = "tenerife"
  129 + if lst:
  130 + irf = "DESY20140105_50h"
  131 + else:
  132 + irf = "DESY20140105_50h_noLST"
  133 + exposure_north += duration
  134 + else:
  135 + caldb = "aar"
  136 + if lst:
  137 + irf = "DESY20140105_50h"
  138 + else:
  139 + irf = "DESY20140105_50h_noLST"
  140 + exposure_south += duration
  141 +
  142 + # Set IRF information
74 143 obs['caldb'] = caldb
75 144 obs['irf'] = irf
76 145  
  146 + # Append observation
  147 + obsdef.append(obs)
  148 +
  149 + # Dump statistics
  150 + print("Number of pointings: %d (%.2f sec)" % (n_pnt,duration))
  151 + print("South array .......: %.2f hours (%.2f%%)" % \
  152 + (exposure_south/3600.0, \
  153 + 100.0*exposure_south/(exposure_south+exposure_north)))
  154 + print("North array .......: %.2f hours (%.2f%%)" % \
  155 + (exposure_north/3600.0, \
  156 + 100.0*exposure_north/(exposure_south+exposure_north)))
  157 +
77 158 # Return observation definition
78 159 return obsdef
79 160  
... ... @@ -81,7 +162,7 @@ def set_double_row(lmin=-30.0, lmax=30.0, separation=3.0, hours=100.0, site=None
81 162 # =========================== #
82 163 # Setup Galactic plane survey #
83 164 # =========================== #
84   -def set_gps():
  165 +def set_gps(separation=3.0, bmin=-1.3, bmax=1.3, lst=True):
85 166 """
86 167 Setup Galactic plane survey.
87 168 """
... ... @@ -89,19 +170,72 @@ def set_gps():
89 170 obsdef = []
90 171  
91 172 # Add inner region South
92   - obsdef.extend(set_double_row(lmin=-60.0, lmax=60.0, hours=780, site="South"))
  173 + obsdef.extend(set_patch(lmin=-60.0, lmax=60.0, bmin=bmin, bmax=bmax, \
  174 + separation=separation, hours=780, \
  175 + site="South", lst=lst))
93 176  
94 177 # Add Vela & Carina region
95   - obsdef.extend(set_double_row(lmin=240.0, lmax=300.0, hours=180, site="South"))
  178 + obsdef.extend(set_patch(lmin=240.0, lmax=300.0, bmin=bmin, bmax=bmax, \
  179 + separation=separation, hours=180, \
  180 + site="South", lst=lst))
96 181  
97 182 # Add 210-240 region
98   - obsdef.extend(set_double_row(lmin=210.0, lmax=240.0, hours=60, site="South"))
  183 + obsdef.extend(set_patch(lmin=210.0, lmax=240.0, bmin=bmin, bmax=bmax, \
  184 + separation=separation, hours=60, \
  185 + site="South", lst=lst))
99 186  
100 187 # Add Cygnus, Perseus
101   - obsdef.extend(set_double_row(lmin=60.0, lmax=150.0, hours=450, site="North"))
  188 + obsdef.extend(set_patch(lmin=60.0, lmax=150.0, bmin=bmin, bmax=bmax, \
  189 + separation=separation, hours=450, \
  190 + site="North", lst=lst))
102 191  
103 192 # Add Anticentre
104   - obsdef.extend(set_double_row(lmin=150.0, lmax=210.0, hours=150, site="North"))
  193 + obsdef.extend(set_patch(lmin=150.0, lmax=210.0, bmin=bmin, bmax=bmax, \
  194 + separation=separation, hours=150, \
  195 + site="North", lst=lst))
  196 +
  197 + # Return observation definition
  198 + return obsdef
  199 +
  200 +
  201 +# ========================== #
  202 +# Setup Extragalactic survey #
  203 +# ========================== #
  204 +def set_extgal(separation=3.0, lst=True):
  205 + """
  206 + Setup Extragalactic survey.
  207 + """
  208 + # Initialise observation definition
  209 + obsdef = []
  210 +
  211 + # Set patch
  212 + obsdef.extend(set_patch(lmin=-90.0, lmax=90.0, bmin=+5.0, bmax=+90.0, \
  213 + separation=separation, hours=1000, \
  214 + site="Automatic", lst=lst, autodec=25.0))
  215 +
  216 + # Return observation definition
  217 + return obsdef
  218 +
  219 +
  220 +# ========================= #
  221 +# Setup Galactic centre KSP #
  222 +# ========================= #
  223 +def set_gc(lst=True):
  224 + """
  225 + Setup Galactic centre KSP.
  226 + """
  227 + # Initialise observation definition
  228 + obsdef = []
  229 +
  230 + # Central wobble
  231 + obsdef.extend(set_patch(lmin=-0.5, lmax=0.5, bmin=-0.5, bmax=0.5, \
  232 + separation=0.1, hours=525, \
  233 + site="South", lst=lst))
  234 +
  235 + # Extended region
  236 + obsdef.extend(set_patch(lmin=-10.0, lmax=10.0, bmin=-10.0, bmax=10.0, \
  237 + separation=1.5, hours=300, \
  238 + site="South", lst=lst))
105 239  
106 240 # Return observation definition
107 241 return obsdef
... ... @@ -151,27 +285,50 @@ def write_obsdef(filename, obsdef):
151 285 # ======================== #
152 286 if __name__ == '__main__':
153 287 """
154   - CTA survey simulation.
  288 + CTA pointing pattern generation.
155 289 """
156 290 # Initialise flags
157 291 need_help = False
158 292  
159 293 # Test for command line arguments
160   - print(sys.argv[0])
  294 + #print(sys.argv[0])
161 295 if (len(sys.argv) > 1):
162 296 if sys.argv[1] == "-h":
163 297 need_help = True
164 298 else:
165   - need_help = True
  299 + obsname = sys.argv[1]
  300 + else:
  301 + need_help = True
166 302  
167 303 # Print help if needed and exit
168 304 if need_help:
169 305 print("Usage: make_pointing.py [OPTIONS]")
170 306 print(" -h Display this usage message")
  307 + print(" gps Galactic plane scan (2 row scheme)")
  308 + print(" gps3 Galactic plane scan (3 row scheme)")
  309 + print(" extgal Extragalactic survey")
  310 + print(" gc Galactic centre")
171 311 sys.exit()
172 312  
173   - # Test
174   - obsdef = set_gps()
  313 + # Galactic plane scan
  314 + if obsname == "gps":
  315 + obsdef = set_gps(lst=True)
  316 + write_obsdef("gps.dat", obsdef)
  317 + elif obsname == "gps3":
  318 + obsdef = set_gps(separation=1.5, lst=True)
  319 + write_obsdef("gps3.dat", obsdef)
  320 +
  321 + # Extragalactic survey
  322 + elif obsname == "extgal":
  323 + obsdef = set_extgal(lst=True)
  324 + write_obsdef("extgal.dat", obsdef)
  325 +
  326 + # Galactic centre
  327 + elif obsname == "gc":
  328 + obsdef = set_gc(lst=True)
  329 + write_obsdef("gc.dat", obsdef)
175 330  
176   - # Write observation definition file
177   - write_obsdef("gps.dat", obsdef)
  331 + # Invalid pattern
  332 + else:
  333 + print("Unknown option \""+obsname+"\"")
  334 +
... ...