Alexander Ziegler / gammalib

Commit 8945de38f330c9c0723e2a42aed0a30e939cc346

Authored by Jürgen Knödlseder
1 parent bb549118

Implement new variable naming

Inline Side-by-side
Showing 1 changed file with 79 additions and 87 deletions
inst/cta/test/test_irf_elliptical.py
  View file @8945de3
... ... @@ -82,7 +82,7 @@ def cta_roi_arclength(rad, dist, roi):
82 82 # ======================= #
83 83 # Show ellptical response #
84 84 # ======================= #
85   -def show_response(centre, semimajor, semiminor, posang, obsDir, pntDir, rhos):
  85 +def show_response(centre, semimajor, semiminor, posang, obsDir, pntDir, delta_max, rhos, log=False):
86 86 """
87 87 Show ellptical response
88 88 """
... ... @@ -90,55 +90,39 @@ def show_response(centre, semimajor, semiminor, posang, obsDir, pntDir, rhos):
90 90 plt.figure(1)
91 91 plt.title("Elliptical response")
92 92  
93   - # Set some constants
94   - delta_max = 0.5 * gammalib.deg2rad
95   -
96   - # Put posang in [-pi,pi]
97   - #if (posang > gammalib.pi):
98   - # posang -= gammalib.pi
99   - #elif (posang < -gammalib.pi):
100   - # posang += gammalib.pi
  93 + # Convert delta_max to radians
  94 + delta_max *= gammalib.deg2rad
101 95  
102 96 # Compute some stuff
103   - cos_pa = math.cos(posang * gammalib.deg2rad)
104   - sin_pa = math.sin(posang * gammalib.deg2rad)
105   - zeta = centre.dist(obsDir)
106   - obsOmega = centre.posang(obsDir)
107   - eta = pntDir.dist(obsDir)
108   - lamda = centre.dist(pntDir)
109   - src_max = semimajor * gammalib.deg2rad
110   -
111   - # Omega_0
112   - omega0 = 0.0
113   - denom = math.sin(lamda*gammalib.deg2rad) * math.sin(zeta*gammalib.deg2rad)
114   - if (denom != 0.0):
115   - arg = (math.cos(eta*gammalib.deg2rad) - math.cos(lamda*gammalib.deg2rad) * \
116   - math.cos(zeta*gammalib.deg2rad))/denom;
117   - omega0 = gammalib.acos(arg)
118   - pa_pnt = centre.posang(pntDir)
  97 + cos_pa = math.cos(posang * gammalib.deg2rad)
  98 + sin_pa = math.sin(posang * gammalib.deg2rad)
  99 + rho_obs = centre.dist(obsDir)
  100 + posangle_obs = centre.posang(obsDir)
  101 + rho_pnt = centre.dist(pntDir)
  102 + posangle_pnt = centre.posang(pntDir)
  103 + omega_pnt = posangle_pnt - posangle_obs
  104 + src_max = semimajor * gammalib.deg2rad
119 105  
120 106 # Rho limits
121 107 rho_min = 0.0
122   - rho_max = zeta + delta_max
123   - if (zeta > delta_max):
124   - rho_min = zeta - delta_max
  108 + rho_max = rho_obs + delta_max
  109 + if (rho_obs > delta_max):
  110 + rho_min = rho_obs - delta_max
125 111 if rho_max > src_max:
126 112 rho_max = src_max
127 113  
128 114 # Setup ellipse as vector
129   - ellipse_x = []
130   - ellipse_y = []
  115 + x = []
  116 + y = []
131 117 for angle in range(0, 360):
132 118 radians = angle * gammalib.deg2rad
133 119 cos_radians = math.cos(radians)
134 120 sin_radians = math.sin(radians)
135   - x = semimajor * cos_radians * sin_pa - semiminor * sin_radians * cos_pa
136   - y = semimajor * cos_radians * cos_pa + semiminor * sin_radians * sin_pa
137   - ellipse_x.append(x)
138   - ellipse_y.append(y)
  121 + x.append(semimajor * cos_radians * sin_pa - semiminor * sin_radians * cos_pa)
  122 + y.append(semimajor * cos_radians * cos_pa + semiminor * sin_radians * sin_pa)
139 123  
140 124 # Plot ellipse
141   - plt.plot(ellipse_x, ellipse_y, 'r-')
  125 + plt.plot(x, y, 'r-')
142 126 plt.xlim([ 1.5*semimajor,-1.5*semimajor])
143 127 plt.ylim([-1.5*semimajor, 1.5*semimajor])
144 128 plt.plot([0.0], [0.0], 'ro')
... ... @@ -150,26 +134,24 @@ def show_response(centre, semimajor, semiminor, posang, obsDir, pntDir, rhos):
150 134  
151 135 # Draw connecting line model-photon'
152 136 n = 10
153   - step = zeta/n * gammalib.rad2deg
  137 + step = rho_obs/n * gammalib.rad2deg
154 138 x = []
155 139 y = []
156 140 for i in range(n+1):
157   - x.append(i*step * math.sin(obsOmega))
158   - y.append(i*step * math.cos(obsOmega))
  141 + x.append(i*step * math.sin(posangle_obs))
  142 + y.append(i*step * math.cos(posangle_obs))
159 143 plt.plot(x, y, 'b-')
160 144  
161 145 # Plot Psf validity circle
162   - circle_x = []
163   - circle_y = []
  146 + x = []
  147 + y = []
164 148 for angle in range(0, 360):
165 149 radians = angle * gammalib.deg2rad
166 150 cos_radians = math.cos(radians)
167 151 sin_radians = math.sin(radians)
168   - x = delta_max * cos_radians * gammalib.rad2deg + x_obs
169   - y = delta_max * sin_radians * gammalib.rad2deg + y_obs
170   - circle_x.append(x)
171   - circle_y.append(y)
172   - plt.plot(circle_x, circle_y, 'b-')
  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-')
173 155  
174 156 # Plot pointing direction
175 157 x = pntDir.ra_deg() - centre.ra_deg()
... ... @@ -177,28 +159,29 @@ def show_response(centre, semimajor, semiminor, posang, obsDir, pntDir, rhos):
177 159 plt.plot([x], [y], 'go')
178 160  
179 161 # Print some stuff
180   - print("Distance p'-model (zeta) .............: %.2f deg" % (zeta*gammalib.rad2deg))
181   - print("PA of p' from model (obsOmega)........: %.2f deg" % (obsOmega*gammalib.rad2deg))
182   - print("Distance p'-pointing (eta)............: %.2f deg" % (eta*gammalib.rad2deg))
183   - print("Distance model-pointing (lambda) .....: %.2f deg" % (lamda*gammalib.rad2deg))
184   - print("Omega_0 ..............................: %.2f deg" % (omega0*gammalib.rad2deg))
185   - print("PA of pointing from model ............: %.2f deg" % (pa_pnt*gammalib.rad2deg))
186   - print("Rho interval .........................: %.2f-%.2f deg" % (rho_min*gammalib.rad2deg,rho_max*gammalib.rad2deg))
  162 + print("Distance model-event (rho_obs) ..........: %.2f deg" % (rho_obs*gammalib.rad2deg))
  163 + print("PA of event from model (posangle_obs) ...: %.2f deg" % (posangle_obs*gammalib.rad2deg))
  164 + print("Distance model-pointing (rho_pnt) .......: %.2f deg" % (rho_pnt*gammalib.rad2deg))
  165 + print("PA of pointing from model (posangle_pnt) : %.2f deg" % (posangle_pnt*gammalib.rad2deg))
  166 + print("Azimuth of pointing (omega_pnt) .........: %.2f deg" % (omega_pnt*gammalib.rad2deg))
  167 + print("Rho interval ............................: %.2f-%.2f deg" % (rho_min*gammalib.rad2deg,rho_max*gammalib.rad2deg))
187 168  
188 169 # Plot rhos
189 170 for rho in rhos:
190 171  
191 172 # Print header
192   - print("Rho ..................................: %.2f deg" % (rho))
  173 + if log:
  174 + print("Rho ..................................: %.2f deg" % (rho))
193 175  
194 176 # Compute some stuff
195   - rho_rad = rho * gammalib.deg2rad
  177 + rho_rad = rho * gammalib.deg2rad
196 178  
197 179 # Compute half length of arc that lies within Psf validity circle
198   - domega = 0.5 * cta_roi_arclength(rho_rad, zeta, delta_max)
  180 + domega = 0.5 * cta_roi_arclength(rho_rad, rho_obs, delta_max)
199 181  
200 182 # Print domega
201   - print(" domega .............................: %.2f deg" % (domega*gammalib.rad2deg))
  183 + if log:
  184 + print(" domega .............................: %.2f deg" % (domega*gammalib.rad2deg))
202 185  
203 186 # Continue only if arc length is positive
204 187 if (domega > 0.0):
... ... @@ -206,10 +189,10 @@ def show_response(centre, semimajor, semiminor, posang, obsDir, pntDir, rhos):
206 189 # Precompute cosine and sine terms for azimuthal integration
207 190 cos_rho = math.cos(rho_rad);
208 191 sin_rho = math.sin(rho_rad);
209   - cos_psf = cos_rho * math.cos(zeta)
210   - sin_psf = sin_rho * math.sin(zeta)
211   - cos_ph = cos_rho * math.cos(lamda)
212   - sin_ph = sin_rho * math.sin(lamda)
  192 + cos_psf = cos_rho * math.cos(rho_obs)
  193 + sin_psf = sin_rho * math.sin(rho_obs)
  194 + cos_ph = cos_rho * math.cos(rho_pnt)
  195 + sin_ph = sin_rho * math.sin(rho_pnt)
213 196  
214 197 # Full containment within ellipse
215 198 if (rho < semiminor):
... ... @@ -219,7 +202,8 @@ def show_response(centre, semimajor, semiminor, posang, obsDir, pntDir, rhos):
219 202 omega_max = +domega
220 203  
221 204 # Print
222   - print(" Omega range ........................: %.2f - %.2f deg" % (omega_min*gammalib.rad2deg,omega_max*gammalib.rad2deg))
  205 + if log:
  206 + print(" Omega range ........................: %.2f - %.2f deg" % (omega_min*gammalib.rad2deg,omega_max*gammalib.rad2deg))
223 207  
224 208 # Plot segment
225 209 n = 100
... ... @@ -228,8 +212,8 @@ def show_response(centre, semimajor, semiminor, posang, obsDir, pntDir, rhos):
228 212 step = (omega_max - omega_min)/n
229 213 for i in range(n+1):
230 214 radians = i*step + omega_min
231   - cos_radians = math.sin(radians+obsOmega)
232   - sin_radians = math.cos(radians+obsOmega)
  215 + cos_radians = math.sin(radians+posangle_obs)
  216 + sin_radians = math.cos(radians+posangle_obs)
233 217 x.append(rho * cos_radians)
234 218 y.append(rho * sin_radians)
235 219 plt.plot(x, y, 'k-')
... ... @@ -244,7 +228,8 @@ def show_response(centre, semimajor, semiminor, posang, obsDir, pntDir, rhos):
244 228 omega_width = math.acos(math.sqrt(arg1/arg2))
245 229  
246 230 # Print Omega width
247   - print(" omega_width ........................: %.2f deg" % (omega_width*gammalib.rad2deg))
  231 + if log:
  232 + print(" omega_width ........................: %.2f deg" % (omega_width*gammalib.rad2deg))
248 233  
249 234 #
250 235 if (omega_width > 0.0):
... ... @@ -253,7 +238,7 @@ def show_response(centre, semimajor, semiminor, posang, obsDir, pntDir, rhos):
253 238 # to connecting line between model and observed photon. Make
254 239 # sure that omega_0 is within [-pi,pi] thus that the omega
255 240 # intervals are between [-2pi,2pi]
256   - omega_0 = posang * gammalib.deg2rad - obsOmega
  241 + omega_0 = posang * gammalib.deg2rad - posangle_obs
257 242 if (omega_0 > gammalib.pi):
258 243 omega_0 -= gammalib.pi
259 244 elif (omega_0 < -gammalib.pi):
... ... @@ -265,8 +250,10 @@ def show_response(centre, semimajor, semiminor, posang, obsDir, pntDir, rhos):
265 250 omega2_min = omega1_min + gammalib.pi
266 251 omega2_max = omega1_max + gammalib.pi
267 252  
268   - print(" Omega range 1 (initial) ............: %.2f - %.2f deg" % (omega1_min*gammalib.rad2deg,omega1_max*gammalib.rad2deg))
269   - print(" Omega range 2 (initial) ............: %.2f - %.2f deg" % (omega2_min*gammalib.rad2deg,omega2_max*gammalib.rad2deg))
  253 + # Logging
  254 + if log:
  255 + print(" Omega range 1 (initial) ............: %.2f - %.2f deg" % (omega1_min*gammalib.rad2deg,omega1_max*gammalib.rad2deg))
  256 + print(" Omega range 2 (initial) ............: %.2f - %.2f deg" % (omega2_min*gammalib.rad2deg,omega2_max*gammalib.rad2deg))
270 257  
271 258 # Limit each of the intervals to [omega_min, omega_max]
272 259 omega_min = -domega
... ... @@ -317,10 +304,11 @@ def show_response(centre, semimajor, semiminor, posang, obsDir, pntDir, rhos):
317 304 omega2_max = 0.0
318 305  
319 306 # Print results
320   - print(" omega_0 ............................: %.2f deg" % (omega_0*gammalib.rad2deg))
321   - print(" Omega range ........................: %.2f - %.2f deg" % (omega_min*gammalib.rad2deg,omega_max*gammalib.rad2deg))
322   - print(" Omega range 1 ......................: %.2f - %.2f deg" % (omega1_min*gammalib.rad2deg,omega1_max*gammalib.rad2deg))
323   - print(" Omega range 2 ......................: %.2f - %.2f deg" % (omega2_min*gammalib.rad2deg,omega2_max*gammalib.rad2deg))
  307 + if log:
  308 + print(" omega_0 ............................: %.2f deg" % (omega_0*gammalib.rad2deg))
  309 + print(" Omega range ........................: %.2f - %.2f deg" % (omega_min*gammalib.rad2deg,omega_max*gammalib.rad2deg))
  310 + print(" Omega range 1 ......................: %.2f - %.2f deg" % (omega1_min*gammalib.rad2deg,omega1_max*gammalib.rad2deg))
  311 + print(" Omega range 2 ......................: %.2f - %.2f deg" % (omega2_min*gammalib.rad2deg,omega2_max*gammalib.rad2deg))
324 312  
325 313 # Draw Omega range lines and circle segments
326 314 if (omega1_min < omega1_max):
... ... @@ -329,14 +317,14 @@ def show_response(centre, semimajor, semiminor, posang, obsDir, pntDir, rhos):
329 317 #x = []
330 318 #y = []
331 319 #for i in range(n+1):
332   - # x.append(i*step * math.sin(omega1_min+obsOmega))
333   - # y.append(i*step * math.cos(omega1_min+obsOmega))
  320 + # x.append(i*step * math.sin(omega1_min+posangle_obs))
  321 + # y.append(i*step * math.cos(omega1_min+posangle_obs))
334 322 #plt.plot(x, y, 'k-.')
335 323 #x = []
336 324 #y = []
337 325 #for i in range(n+1):
338   - # x.append(i*step * math.sin(omega1_max+obsOmega))
339   - # y.append(i*step * math.cos(omega1_max+obsOmega))
  326 + # x.append(i*step * math.sin(omega1_max+posangle_obs))
  327 + # y.append(i*step * math.cos(omega1_max+posangle_obs))
340 328 #plt.plot(x, y, 'k-.')
341 329 n = 100
342 330 x = []
... ... @@ -344,8 +332,8 @@ def show_response(centre, semimajor, semiminor, posang, obsDir, pntDir, rhos):
344 332 step = (omega1_max - omega1_min)/n
345 333 for i in range(n+1):
346 334 radians = i*step + omega1_min
347   - cos_radians = math.sin(radians+obsOmega)
348   - sin_radians = math.cos(radians+obsOmega)
  335 + cos_radians = math.sin(radians+posangle_obs)
  336 + sin_radians = math.cos(radians+posangle_obs)
349 337 x.append(rho * cos_radians)
350 338 y.append(rho * sin_radians)
351 339 plt.plot(x, y, 'k-')
... ... @@ -356,14 +344,14 @@ def show_response(centre, semimajor, semiminor, posang, obsDir, pntDir, rhos):
356 344 #x = []
357 345 #y = []
358 346 #for i in range(n+1):
359   - # x.append(i*step * math.sin(omega2_min+obsOmega))
360   - # y.append(i*step * math.cos(omega2_min+obsOmega))
  347 + # x.append(i*step * math.sin(omega2_min+posangle_obs))
  348 + # y.append(i*step * math.cos(omega2_min+posangle_obs))
361 349 #plt.plot(x, y, 'k-.')
362 350 #x = []
363 351 #y = []
364 352 #for i in range(n+1):
365   - # x.append(i*step * math.sin(omega2_max+obsOmega))
366   - # y.append(i*step * math.cos(omega2_max+obsOmega))
  353 + # x.append(i*step * math.sin(omega2_max+posangle_obs))
  354 + # y.append(i*step * math.cos(omega2_max+posangle_obs))
367 355 #plt.plot(x, y, 'k-.')
368 356 n = 100
369 357 x = []
... ... @@ -371,8 +359,8 @@ def show_response(centre, semimajor, semiminor, posang, obsDir, pntDir, rhos):
371 359 step = (omega2_max - omega2_min)/n
372 360 for i in range(n+1):
373 361 radians = i*step + omega2_min
374   - cos_radians = math.sin(radians+obsOmega)
375   - sin_radians = math.cos(radians+obsOmega)
  362 + cos_radians = math.sin(radians+posangle_obs)
  363 + sin_radians = math.cos(radians+posangle_obs)
376 364 x.append(rho * cos_radians)
377 365 y.append(rho * sin_radians)
378 366 plt.plot(x, y, 'k-')
... ... @@ -381,8 +369,8 @@ def show_response(centre, semimajor, semiminor, posang, obsDir, pntDir, rhos):
381 369 #x = []
382 370 #y = []
383 371 #for i in range(10):
384   - # x.append(i*0.1 * math.sin(0.0+obsOmega))
385   - # y.append(i*0.1 * math.cos(0.0+obsOmega))
  372 + # x.append(i*0.1 * math.sin(0.0+posangle_obs))
  373 + # y.append(i*0.1 * math.cos(0.0+posangle_obs))
386 374 #plt.plot(x, y, 'k-')
387 375  
388 376 # Build circle
... ... @@ -425,6 +413,9 @@ if __name__ == &#39;__main__&#39;:
425 413 print("Matplotlib is not installed on your system. Abort.")
426 414 sys.exit()
427 415  
  416 + # Set Psf size
  417 + delta_max = 0.4
  418 +
428 419 # Set model parameters
429 420 centre = gammalib.GSkyDir()
430 421 centre.radec_deg(0.0, 0.0)
... ... @@ -434,7 +425,8 @@ if __name__ == &#39;__main__&#39;:
434 425  
435 426 # Set observed photon direction
436 427 obsDir = gammalib.GSkyDir()
437   - obsDir.radec_deg(0.4, -0.3)
  428 + #obsDir.radec_deg(0.4, -0.3)
  429 + obsDir.radec_deg(-0.4, -0.3)
438 430  
439 431 # Set pointing direction
440 432 pntDir = gammalib.GSkyDir()
... ... @@ -444,5 +436,5 @@ if __name__ == &#39;__main__&#39;:
444 436 rhos = [i*0.02 for i in range(51)]
445 437  
446 438 # Show elliptical response geometry
447   - show_response(centre, semimajor, semiminor, posang, obsDir, pntDir, rhos)
  439 + show_response(centre, semimajor, semiminor, posang, obsDir, pntDir, delta_max, rhos)
448 440  
... ...