pylal.sky

1 #!/usr/bin/env python 2 3 # Copyright (C) 2011 Duncan M. Macleod 4 # 5 # This program is free software; you can redistribute it and/or modify it 6 # under the terms of the GNU General Public License as published by the 7 # Free Software Foundation; either version 3 of the License, or (at your 8 # option) any later version. 9 # 10 # This program is distributed in the hope that it will be useful, but 11 # WITHOUT ANY WARRANTY; without even the implied warranty of 12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 13 # Public License for more details. 14 # 15 # You should have received a copy of the GNU General Public License along 16 # with this program; if not, write to the Free Software Foundation, Inc., 17 # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 18 19 # ============================================================================= 20 # Preamble 21 # ============================================================================= 22 23 from __future__ import division 24 import numpy as np,re,sys 25 26 import lal 27 28 from glue import git_version 29 from glue.ligolw import table 30 31 from pylal import inject,date 32 from pylal import inject 33 from pylal.xlal.datatypes.ligotimegps import LIGOTimeGPS 34 35 __author__ = "Duncan M. Macleod <duncan.macleod@astro.cf.ac.uk>" 36 __version__ = git_version.id 37 __date__ = git_version.date 38 39 # useful regex 40 _comment_regex = re.compile(r"\s*([#;].*)?\Z", re.DOTALL) 41 _delim_regex = re.compile('[,\s]') 42 43 # ============================================================================= 44 # Sky Positions structure 45 # ============================================================================= 46

47 -class SkyPositionTable(table.Table):

48 """ 49 glue.ligolw.table.Table holding SkyPosition objects. 50 """ 51 52 tableName = "sky_positions" 53 valid_columns = {"process_id": "ilwd:char",\ 54 "latitude": "real_4",\ 55 "longitude": "real_4",\ 56 "probability": "real_4",\ 57 "solid_angle": "real_4",\ 58 "system": "lstring"} 59

60 - def rotate(self, R):

61 rotated = table.new_from_template(self) 62 for row in self: 63 rotated.append(row.rotate(R)) 64 return rotated

65

66 - def normalize(self):

67 for row in self: 68 row.normalize()

69

70 - def parseTimeDelayDegeneracy(self, ifos, gpstime=lal.GPSTimeNow(),\ 71 dt=0.0005):

72 73 # get detectors 74 detectors = [inject.cached_detector.get(inject.prefix_to_name[ifo])\ 75 for ifo in ifos] 76 77 timeDelays = [] 78 degenerate = [] 79 80 new = table.new_from_template(self) 81 for n,row in enumerate(self): 82 # get all time delays for this point 83 timeDelays.append([]) 84 for i in xrange(len(ifos)): 85 for j in xrange(i+1, len(ifos)): 86 timeDelays[n].append(date.XLALArrivalTimeDiff(\ 87 detectors[i].location,\ 88 detectors[j].location,\ 89 row.longitude,\ 90 row.latitude,\ 91 LIGOTimeGPS(gpstime))) 92 # skip the first point 93 if n==0: 94 degenerate.append(False) 95 new.append(row) 96 continue 97 else: 98 degenerate.append(True) 99 # test this point against all others 100 for i in xrange(0,n): 101 # if check point is degenerate, skip 102 if degenerate[i]: 103 continue 104 # check each time delay individually 105 for j in xrange(0, len(timeDelays[n])): 106 # if this time delay is non-degenerate the point is valid 107 if np.fabs(timeDelays[i][j]-timeDelays[n][j]) >= dt: 108 degenerate[n] = False 109 break 110 else: 111 degenerate[n] = True 112 if degenerate[n]: 113 break 114 115 if not degenerate[n]: 116 new.append(row) 117 118 return new

119 120

121 -class SkyPosition(object):

122 123 __slots__ = SkyPositionTable.valid_columns.keys() 124

125 - def __iter__(self):

126 return iter(tuple((self.longitude, self.latitude)))

127

128 - def __repr__(self):

129 return "SkyPosition(%s, %s)" % (self.longitude, self.latitude)

130

131 - def __str__(self):

132 return "(%s, %s)" % (self.longitude, self.latitude)

133

134 - def rotate(self, R):

135 """ 136 Apply the 3x3 rotation matrix R to this SkyPosition. 137 """ 138 139 cart = SphericalToCartesian(self) 140 rot = np.dot(np.asarray(R), cart) 141 new = CartesianToSpherical(rot, system=self.system) 142 new.normalize() 143 try: 144 new.probability = self.probability 145 except AttributeError: 146 new.probability = None 147 try: 148 new.solid_angle = self.solid_angle 149 except AttributeError: 150 new.solid_angle = None 151 152 return new

153

154 - def normalize(self):

155 """ 156 Normalize this SkyPosition to be within standard limits 157 [0 <= longitude < 2pi] and [-pi < latitude <= pi] 158 """ 159 160 # first unwind positions, i.e. make sure that: 161 # [0 <= alpha < 2pi] and [-pi < delta <= pi] 162 163 # normalise longitude 164 while self.longitude < 0: 165 self.longitude += lal.TWOPI 166 while self.longitude >= lal.TWOPI: 167 self.longitude -= lal.TWOPI 168 169 # normalise latitude 170 while self.latitude <= -lal.PI: 171 self.latitude += lal.TWOPI 172 while self.latitude > lal.TWOPI: 173 self.latitude -= lal.TWOPI 174 175 # 176 # now get latitude into canonical interval [-pi/2, pi/2) 177 # 178 179 if self.latitude > lal.PI_2: 180 self.latitude = lal.PI - self.latitude 181 if self.latitude < lal.PI: 182 self.longitude += lal.PI 183 else: 184 self.longitude -= lal.PI 185 186 if self.latitude < -lal.PI_2: 187 self.latitude = -lal.PI - self.latitude 188 if self.longitude < lal.PI: 189 self.longitude += lal.PI 190 else: 191 self.longitude -= lal.PI

192

193 - def opening_angle(self, other):

194 """ 195 Calcalate the opening angle between this SkyPosition and the other 196 SkyPosition 197 """ 198 199 if self == other: 200 return 0.0 201 202 s = np.sin 203 c = np.cos 204 205 return np.arccos(s(self.latitude) * s(other.latitude) +\ 206 c(self.latitude) * c(other.latitude) *\ 207 c(self.longitude - other.longitude))

208

209 - def get_time_delay(self, ifo1, ifo2, gpstime):

210 """ 211 Return the time delay for this SkyPosition between arrival at ifo1 212 relative to ifo2, for the given gpstime. 213 """ 214 det1 = inject.cached_detector.get(inject.prefix_to_name[ifo1]) 215 det2 = inject.cached_detector.get(inject.prefix_to_name[ifo2]) 216 217 return date.XLALArrivalTimeDiff(det1.location, det2.location,\ 218 self.longitude, self.latitude,\ 219 LIGOTimeGPS(gpstime))

220 221 222 223 # ============================================================================= 224 # Convert between coordinate systems 225 # ============================================================================= 226 227 # a few last definitions (taken from SkyCoordinates.c) 228 lal.LAL_ALPHAGAL = 3.366032942 229 lal.LAL_DELTAGAL = 0.473477302 230 lal.LAL_LGAL = 0.576 231

232 -def ConvertSkyPosition(skyPoint, system, zenith=None, gpstime=None):

233 """ 234 Convert the SkyPosition object skyPoint from it's current system to the 235 new system. 236 237 Valid systems are : 'horizon', 'geographic', 'equatorial', 'ecliptic',\ 238 'galactic' 239 240 SkyPosition zenith and gpstime should be given as appropriate for 'horizon' 241 and 'geographic' conversions respectively. 242 """ 243 244 valid_systems = ['horizon', 'geographic', 'equatorial', 'ecliptic',\ 245 'galactic'] 246 247 system = system.lower() 248 skyPoint.system = skyPoint.system.lower() 249 250 if system not in valid_systems: 251 raise AttributeError("Unrecognised system='%s'" % system) 252 253 while skyPoint.system != system: 254 255 if skyPoint.system == 'horizon': 256 skyPoint = HorizonToSystem(skyPoint, zenith) 257 258 elif skyPoint.system == 'geographic': 259 if system == 'horizon': 260 if zenith.system == 'geographic': 261 skyPoint = SystemToHorizon(skyPoint, zenith) 262 elif zenith.system == 'equatorial': 263 skyPoint = GeographicToEquatorial(skyPoint, gpstime) 264 else: 265 raise AttibuteError("Cannot convert from geographic to "+\ 266 "horizon with zenith point not in "+\ 267 "geogrphic of equatorial systems.") 268 else: 269 skyPoint = GeographicToEquatorial(skyPoint, gpstime) 270 271 if system == 'horizon' and zenith.system == 'equatorial': 272 skyPoint = SystemToHorizon(skyPoint, zenith) 273 elif system == 'ecliptic': 274 skyPoint = EquatorialToEcliptic(skyPoint) 275 elif system == 'galactic': 276 skyPoint = EquatorialToGalactic(skyPoint) 277 else: 278 skyPoint = EquatorialToGeographic(skyPoint, gpstime) 279 280 elif skyPoint.system == 'ecliptic': 281 skyPoint = EclipticToEquatorial(skyPoint) 282 283 elif skyPoint.system == 'galactic': 284 skyPoint = GalacticToEquatorial(skyPoint) 285 286 return skyPoint

287

288 -def HorizonToSystem(input, zenith):

289 """ 290 Convert the SkyPosition object input from 'horizon' to the inherited system 291 using the SkyPosition zenith 292 """ 293 294 if input.system.lower() != 'horizon': 295 raise AttributeError("Cannot convert from horizon for point in "+\ 296 "system='%s'" % input.system.lower()) 297 298 if zenith.system.lower() != 'equatorial'\ 299 and zenith.system.lower() != 'geographic': 300 raise AttributeError("zenith must have coordinate system = "+\ 301 "'equatorial' or 'geographic'") 302 303 # intermediates 304 sinP = np.sin(zenith.latitude) 305 cosP = np.cos(zenith.latitude) 306 sina = np.sin(input.latitude) 307 cosa = np.cos(input.latitude) 308 sinA = np.sin(input.longitude) 309 cosA = np.cos(input.longitude) 310 311 # final components 312 sinD = sina*sinP + cosa*cosA*cosP 313 sinH = cosa*sinA 314 cosH = sina*cosP - cosa*cosA*sinP 315 316 # output 317 output = SkyPosition() 318 output.system = zenith.system 319 output.latitude = np.arcsin(sinD) 320 output.longitude = zenith.longitude - np.arctan2(sinH, cosH) 321 output.probability = input.probability 322 output.solid_angle = input.solid_angle 323 output.normalize() 324 325 return output

326

327 -def SystemToHorizon(input, zenith):

328 """ 329 Convert the SkyPosition object input from the inherited system to 'horizon' 330 using the SkyPosition zenith 331 """ 332 333 if input.system.lower() != zenith.system.lower(): 334 raise AttributeError("input coordinate system must equal zenith system") 335 336 if zenith.system.lower() != 'equatorial'\ 337 and zenith.system.lower() != 'geographic': 338 raise AttributeError("zenith must have coordinate system = "+\ 339 "'equatorial' or 'geographic'") 340 341 # intermediates 342 h = zenith.longitude - input.longitude 343 sinH = np.sin(h) 344 cosH = np.cos(h) 345 sinP = np.sin(zenith.latitude) 346 cosP = np.cos(zenith.latitude) 347 sinD = np.sin(input.latitude) 348 cosD = np.cos(input.latitude) 349 350 # final components 351 sina = sinD*sinP + cosD*cosP*cosH 352 sinA = cosD*sinH 353 cosA = sinD*cosP - cosD*sinP*cosH 354 355 # output 356 output = SkyPosition() 357 output.system = 'horizon' 358 output.latitude = np.arcsin(sina) 359 output.longitude = np.arctan2(sinA, cosA) 360 output.probability = input.probability 361 output.solid_angle = input.solid_angle 362 output.normalize() 363 364 return output

365

366 -def GeographicToEquatorial(input, gpstime):

367 """ 368 Convert the SkyPosition object input from the inherited 'geographical' 369 system to to 'equatorial'. 370 """ 371 372 if input.system.lower() != 'geographic': 373 raise AttributeError("Input system!='geographic'") 374 375 # get GMST 376 gmst = np.fmod(date.XLALGreenwichSiderealTime(gpstime, 0),\ 377 lal.TWOPI) 378 379 # output 380 output = SkyPosition() 381 output.system = 'equatorial' 382 output.longitude = np.fmod(input.longitude + gmst, lal.TWOPI) 383 output.latitude = input.latitude 384 output.probability = input.probability 385 output.solid_angle = input.solid_angle 386 output.normalize() 387 388 return output

389

390 -def EquatorialToEcliptic(input):

391 """ 392 Convert the SkyPosition object input from the inherited 'equatorial' 393 system to to 'ecliptic'. 394 """ 395 396 # intermediates 397 sinD = np.sin(input.latitude) 398 cosD = np.cos(input.latitude) 399 sinA = np.sin(input.longitude) 400 cosA = np.cos(input.longitude) 401 402 # components 403 sinB = sinD*np.cos(lal.IEARTH)\ 404 - cosD*sinA*np.sin(lal.IEARTH) 405 sinL = cosD*sinA*np.cos(lal.IEARTH)\ 406 + sinD*np.sin(lal.IEARTH) 407 cosL = cosD*cosA 408 409 # output 410 output.system = 'ecliptic' 411 output.latitude = np.arcsin(sinB) 412 output.longitude = np.arctan2(sinL, cosL) 413 output.normalize() 414 415 return output

416

417 -def EquatorialToGalactic(input):

418 """ 419 Convert the SkyPosition object input from the inherited 'equatorial' 420 system to to 'galactic'. 421 """ 422 423 # intermediates. */ 424 a = -lal.LAL_ALPHAGAL + input.longitude 425 sinD = np.sin(input.latitude) 426 cosD = np.cos(input.latitude) 427 sinA = np.sin(a) 428 cosA = np.cos(a) 429 430 # components. */ 431 sinB = cosD*np.cos(lal.LAL_DELTAGAL)*cosA\ 432 + sinD*np.sin(lal.LAL_DELTAGAL) 433 sinL = sinD*np.cos(lal.LAL_DELTAGAL)\ 434 - cosD*cosA*np.sin(lal.LAL_DELTAGAL) 435 cosL = cosD*sinA 436 437 # output 438 output = SkyPosition() 439 output.system = 'galactic' 440 output.latitude = np.arcsin(sinB) 441 output.longitude = np.arctan2(sinL, cosL) + lal.LAL_LGAL 442 output.probability = input.probability 443 output.solid_angle = input.solid_angle 444 output.normalize() 445 446 return output

447

448 -def EquatorialToGeographic(input, gpstime):

449 """ 450 Convert the SkyPosition object input from the inherited 'equatorial' 451 system to to 'geographic'. 452 """ 453 454 if input.system.lower() != 'equatorial': 455 raise AttributeError("Input system is not 'equatorial'") 456 457 # get GMST 458 gmst = np.fmod(date.XLALGreenwichSiderealTime(gpstime, 0),\ 459 lal.TWOPI) 460 461 # output 462 output = SkyPosition() 463 output.system = 'geographic' 464 output.longitude = np.fmod(input.longitude - gmst, lal.TWOPI) 465 output.latitude = input.latitude 466 output.probability = input.probability 467 output.solid_angle = input.solid_angle 468 output.normalize() 469 470 return output

471

472 -def EclipticToEquatorial(input):

473 """ 474 Convert the SkyPosition object input from the inherited 'eliptic' 475 system to to 'equatorial'. 476 """ 477 478 # intermediates 479 sinB = np.sin(input.latitude) 480 cosB = np.cos(input.latitude) 481 sinL = np.sin(input.longitude) 482 cosL = np.cos(input.longitude) 483 484 # components 485 sinD = cosB*sinL*np.sin(lal.IEARTH)\ 486 + sinB*np.cos(lal.IEARTH) 487 sinA = cosB*sinL*np.cos(lal.IEARTH)\ 488 - sinB*np.sin(lal.IEARTH) 489 cosA = cosB*cosL 490 491 # output 492 output.system = 'equatorial' 493 output.latitude = np.arcsin(sinD) 494 output.longitude = np.arctan2(sinA, cosA) 495 output.normalize() 496 497 return output

498 499 # ============================================================================= 500 # Read grid from file 501 # ============================================================================= 502

503 -def fromfile(fobj, loncol=0, latcol=1, probcol=None, angcol=None, system=None,\ 504 degrees=False):

505 """ 506 Returns a SkyPositionTable as read from the file object fobj. 507 """ 508 509 l = SkyPositionTable() 510 511 if degrees: 512 func = np.radians 513 else: 514 func = float 515 516 for line in fobj: 517 line = _comment_regex.split(line)[0] 518 if not line: continue 519 line = re.split(_delim_regex, line) 520 p = SkyPosition() 521 p.longitude = func(float(line[loncol])) 522 p.latitude = func(float(line[latcol])) 523 if probcol: 524 p.probability = line[probcol] 525 else: 526 p.probability = None 527 if angcol: 528 p.solid_angle = line[angcol] 529 else: 530 p.solid_angle = None 531 p.system = system 532 l.append(p) 533 534 return l

535

536 -def tofile(fobj, grid, delimiter=" ", degrees=False):

537 """ 538 Writes the SkyPositionTable object grid to the file object fobj 539 """ 540 541 for p in grid: 542 if degrees: 543 lon = np.degrees(p.longitude) 544 lat = np.degrees(p.latitude) 545 else: 546 lon = p.longitude 547 lat = p.latitude 548 fobj.write("%s%s%s\n" % (lon, delimiter, lat))

549 550 # ============================================================================= 551 # Generate grids 552 # ============================================================================= 553

554 -def SkyPatch(ifos, ra, dec, radius, gpstime, dt=0.0005, sigma=1.65,\ 555 grid='circular'):

556 """ 557 Returns a SkyPositionTable of circular rings emanating from a given 558 central ra and dec. out to the maximal radius. 559 """ 560 561 # form centre point 562 p = SkyPosition() 563 p.longitude = ra 564 p.latitude = dec 565 566 # get detectors 567 ifos.sort() 568 detectors = [] 569 for ifo in ifos: 570 if ifo not in inject.prefix_to_name.keys(): 571 raise ValueError("Interferometer '%s' not recognised." % ifo) 572 detectors.append(inject.cached_detector.get(\ 573 inject.prefix_to_name[ifo])) 574 575 alpha = 0 576 for i in xrange(len(ifos)): 577 for j in xrange(i+1,len(ifos)): 578 # get opening angle 579 baseline = date.XLALArrivalTimeDiff(detectors[i].location,\ 580 detectors[j].location,\ 581 ra, dec, LIGOTimeGPS(gpstime)) 582 ltt = inject.light_travel_time(ifos[i], ifos[j]) 583 angle = np.arccos(baseline/ltt) 584 585 # get window 586 lmin = angle-radius 587 lmax = angle+radius 588 if lmin < lal.PI_2 and lmax > lal.PI_2: 589 l = lal.PI_2 590 elif np.fabs(lal.PI_2-lmin) <\ 591 np.fabs(lal.PI_2-lmax): 592 l = lmin 593 else: 594 l = lmax 595 596 # get alpha 597 dalpha = ltt * np.sin(l) 598 if dalpha > alpha: 599 alpha = dalpha 600 601 # get angular resolution 602 angRes = 2 * dt/alpha 603 604 # generate grid 605 if grid.lower() == 'circular': 606 grid = CircularGrid(angRes, radius) 607 else: 608 raise RuntimeError("Must use grid='circular', others not coded yet") 609 610 # 611 # Need to rotate grid onto (ra, dec) 612 # 613 614 # calculate opening angle from north pole 615 north = [0, 0, 1] 616 angle = np.arccos(np.dot(north, SphericalToCartesian(p))) 617 #angle = north.opening_angle(p) 618 619 # get rotation axis 620 axis = np.cross(north, SphericalToCartesian(p)) 621 axis = axis / np.linalg.norm(axis) 622 623 # rotate grid 624 R = _rotation(axis, angle) 625 grid = grid.rotate(R) 626 627 # 628 # calculate probability density function 629 # 630 631 # assume Fisher distribution in angle from centre 632 kappa = (0.66*radius/sigma)**(-2) 633 634 # compute probability 635 for p in grid: 636 overlap = np.cos(p.opening_angle(grid[0])) 637 p.probability = np.exp(kappa*(overlap-1)) 638 639 probs = [p.probability for p in grid] 640 for p in grid: 641 p.probability = p.probability/max(probs) 642 643 return grid

644

645 -def CircularGrid(resolution, radius):

646 """ 647 Generates a SkyPositionTable of circular grids around the North Pole with 648 given angular resolution and a maximal radius. 649 """ 650 651 l = SkyPositionTable() 652 653 theta = 0 654 while theta < radius+resolution: 655 # get number of points in ring 656 n = max(1, int(np.ceil(lal.TWOPI\ 657 * np.sin(theta)/resolution))) 658 # get phi step 659 dp = lal.TWOPI / n 660 # get solid angle 661 if theta == 0 or theta == lal.PI: 662 dO = 4*lal.PI * np.sin(0.25*resolution)**2 663 else: 664 dO = 4*lal.PI/n * np.sin(theta) * np.sin(0.5*resolution) 665 666 # assign points in ring 667 for j in xrange(n): 668 p = SkyPosition() 669 p.system = 'equatorial' 670 p.longitude = (-lal.PI + dp/2) + dp*j 671 p.latitude = lal.PI_2 - theta 672 p.solid_angle = dO 673 p.probability = 0 674 p.normalize() 675 l.append(p) 676 theta += resolution 677 678 # assign uniform probability distribution 679 totSolidAngle = sum([p.solid_angle for p in l]) 680 for i,p in enumerate(l): 681 l[i].probability = p.solid_angle/totSolidAngle 682 683 return l

684

685 -def TwoSiteSkyGrid(ifos, gpstime, dt=0.0005, sky='half', point=None):

686 """ 687 Generates a SkyPositionTable for a two-site all-sky grid, covering either 688 a hemisphere (sky='half') or full sphere (sky='full'). 689 690 The grid can be forced to pass through the given SkyPosition point if 691 required. 692 """ 693 694 # remove duplicate site references 695 ifos = parse_sites(ifos) 696 assert len(ifos)==2, "More than two sites in given list of detectors" 697 698 # get detectors 699 detectors = [inject.cached_detector.get(inject.prefix_to_name[ifo])\ 700 for ifo in ifos] 701 702 # get light travel time 703 ltt = inject.light_travel_time(ifos[0], ifos[1]) 704 705 # get number of points 706 max = int(np.floor(ltt/dt)) 707 708 # get baseline 709 baseline = detectors[1].location - detectors[0].location 710 baseline /= np.linalg.norm(baseline) 711 712 # 713 # construct rotation 714 # 715 716 # xaxis points overhead of detector 1 or given point 717 # zaxis is baseline, or something else 718 if point: 719 xaxis = SphericalToCartesian(point) 720 xaxis /= np.linalg.norm(xaxis) 721 zaxis = np.cross(xaxis, baseline) 722 zaxis /= np.linalg.norm(zaxis) 723 else: 724 xaxis = detectors[0].location 725 xaxis /= np.linalg.norm(xaxis) 726 zaxis = baseline 727 yaxis = np.cross(zaxis, xaxis) 728 yaxis /= np.linalg.norm(yaxis) 729 730 yaxis = np.cross(xaxis, zaxis) 731 yaxis /= np.linalg.norm(yaxis) 732 733 # construct Euler rotation 734 lineofnodes = np.cross([0,0,1], zaxis) 735 lineofnodes /= np.linalg.norm(lineofnodes) 736 737 # work out Euler angles 738 alpha = np.arccos(np.dot([1,0,0], lineofnodes)) 739 beta = np.arccos(np.dot([0,0,1], zaxis)) 740 gamma = np.arccos(np.dot(lineofnodes, xaxis)) 741 742 # get rotation 743 R = _rotation_euler(alpha, beta, gamma) 744 745 # construct sky table 746 l = SkyPositionTable() 747 if sky=='half': longs = [0] 748 elif sky=='full': longs = [0,lal.PI] 749 else: raise AttributeError("Sky type \"%s\" not recognised, please use " 750 "'half' or 'full'" % sky) 751 752 for long in longs: 753 for i in xrange(-max, max+1): 754 # get time-delay 755 t = i * dt 756 # if time-delay greater that light travel time, skip 757 if abs(t) >= ltt: continue 758 # construct object 759 e = SkyPosition() 760 e.system = 'geographic' 761 # project time-delay onto prime meridian 762 e.longitude = long 763 e.latitude = lal.PI_2-np.arccos(-t/ltt) 764 e.probability = None 765 e.solid_angle = None 766 e.normalize() 767 e = e.rotate(R) 768 e = GeographicToEquatorial(e, LIGOTimeGPS(gpstime)) 769 if e not in l: 770 l.append(e) 771 772 return l

773

774 -def ThreeSiteSkyGrid(ifos, gpstime, dt=0.0005, tiling='square', sky='half'):

775 """ 776 Generates a SkyPositionTable for a three-site all-sky grid, covering either 777 a hemisphere (sky='half') or full sphere (sky='full'), using either 778 'square' or 'hexagonal tiling. 779 """ 780 781 # remove duplicate site references 782 pop = [] 783 for i,ifo in enumerate(ifos): 784 if i==0: continue 785 site = ifo[0] 786 if site in [x[0] for x in ifos[0:i]]: 787 sys.stderr.write("Duplicate site reference, ignoring %s.\n" % ifo) 788 pop.append(i) 789 for p in pop[::-1]: 790 ifos.pop(p) 791 792 assert len(ifos)==3 793 794 detectors = [] 795 T = [] 796 baseline = [] 797 798 # load detectors 799 for i,ifo in enumerate(ifos): 800 detectors.append(inject.cached_detector.get(inject.prefix_to_name[ifo])) 801 # get light travel time and baseline for other detectors to first 802 if i>=1: 803 T.append(inject.light_travel_time(ifos[0], ifos[i])) 804 baseline.append(detectors[i].location - detectors[0].location) 805 baseline[i-1] /= np.linalg.norm(baseline[i-1]) 806 807 # get angle between baselines 808 angle = np.arccos(np.dot(baseline[0], baseline[1])) 809 810 # 811 # construct tiling vectors 812 # 813 814 tiling = tiling.lower() 815 t = np.zeros(len(baseline)) 816 tdgrid = [] 817 818 # square grid 819 dx = np.array([1,0]) 820 dy = np.array([0,1]) 821 nx = int(np.floor(T[0]/dt)) 822 ny = int(np.floor(T[1]/dt)) 823 824 # hexagonal grid 825 if tiling == 'square': 826 dm = dx 827 dn = dy 828 nm = nx 829 nn = ny 830 x = np.linspace(-nx, nx, 2*nx+1) 831 y = np.linspace(-ny, ny, 2*ny+1) 832 grid = np.array([[i,j] for i in x for j in y]) 833 # hexagonal grid 834 elif re.match('hex', tiling): 835 dm = (2*dx+dy) 836 dn = (-dx+dy) 837 dx = dm - dn 838 # tile grid so that origin gets tiled 839 grid = [] 840 orig = np.array([0,0]) 841 842 # find bottom left point on grid 843 while orig[1] >= -ny: orig -= dn 844 845 # loop over y 846 for y in xrange(-ny,ny+1): 847 orig[0] = orig[0] % dx[0] 848 # work out minimal x value 849 minx = -nx 850 while minx % 3 != orig[0]: minx+=1 851 # tile x 852 x = np.arange(minx, nx+1, dx[0]) 853 # append to grid 854 grid.extend([[i,y] for i in x]) 855 orig += dn 856 orig[0] = orig[0] % dx[0] 857 grid = np.asarray(grid) 858 859 # 860 # Following Rabaste, Chassande-Motin, Piran (2009), construct an ellipse 861 # of physical values in 2D time-delay space for 3 detectors 862 # 863 864 # 865 # (theta, phi) coordinate system is as follows: 866 # detector 1 is origin 867 # z-axis joins positively to detector 2 868 # x-axis is orthogonal to z-axis in plane joining D1, D2, and D3 869 # y-axis forms a right-handed coordinate system 870 # 871 872 # so need to rotate from Rabaste network coordinates into 873 # earth-fixed coordinates at the end 874 875 # set z-axis in Rabaste frame 876 zaxis = baseline[0] 877 878 # work out y-axis in Rabaste frame 879 yaxis = np.cross(zaxis, baseline[1]) 880 yaxis /= np.linalg.norm(yaxis) 881 882 # work out x-axis in Rabaste frame 883 xaxis = np.cross(yaxis, zaxis) 884 xaxis /= np.linalg.norm(xaxis) 885 886 # work out lines of nodes 887 north = np.array([0, 0, 1]) 888 lineofnodes = np.cross(baseline[0], north) 889 lineofnodes /= np.linalg.norm(lineofnodes) 890 891 # work out Euler angles 892 alpha = np.arccos(np.dot(xaxis, lineofnodes)) 893 beta = np.arccos(np.dot(baseline[0], north)) 894 gamma = np.arccos(np.dot(lineofnodes, [1,0,0])) 895 896 # construct rotation 897 R = _rotation_euler(alpha, beta, gamma) 898 899 # 900 # construct sky position table 901 # 902 903 l = SkyPositionTable() 904 # loop over time-delay between D1 and D2 905 k = 0 906 for i in grid: 907 t = dt*i 908 # construct coefficients of elliptical equation in quadratic form 909 A = -T[1]/T[0] * t[0] * np.cos(angle) 910 B = T[1]**2 * ((t[0]/T[0])**2 - np.sin(angle)**2) 911 # test condition for inside ellipse and place point 912 condition = t[1]**2 + 2*A*t[1] + B 913 if condition <= 0: 914 ntheta = np.arccos(-t[0]/T[0]) 915 nphi = np.arccos(-(T[0]*t[1]-T[1]*t[0]*np.cos(angle))/\ 916 (T[1]*np.sqrt(T[0]**2-t[0]**2)*np.sin(angle))) 917 p = SkyPosition() 918 p.system = 'geographic' 919 p.longitude = nphi 920 p.latitude = lal.PI_2 - ntheta 921 p.probability = None 922 p.solid_angle = None 923 p.normalize() 924 p = p.rotate(R) 925 p = GeographicToEquatorial(p, LIGOTimeGPS(gpstime)) 926 l.append(p) 927 if sky=='full': 928 p2 = SkyPosition() 929 p2.longitude = p.longitude 930 p2.latitude = p.latitude + lal.PI 931 p2.system = 'equatorial' 932 p2.normalize() 933 l.append(p2) 934 935 return l

936

937 -def ISOTimeDelayLine(ifos, ra, dec, gpstime=None, n=3):

938 """ 939 Returns the n-point SkyPositionTable describing a line of constant time 940 delay through the given ra and dec. for the given 2-tuple ifos. 941 """ 942 943 if gpstime: 944 gpstime = LIGOTimeGPS(gpstime) 945 946 p = SkyPosition() 947 p.longitude = ra 948 p.latitude = dec 949 p.system = 'equatorial' 950 p.probability = None 951 p.solid_angle = None 952 if gpstime: 953 p = EquatorialToGeographic(p, gpstime) 954 cart = SphericalToCartesian(p) 955 956 # get baseline 957 detectors = [inject.cached_detector.get(inject.prefix_to_name[ifo])\ 958 for ifo in ifos] 959 baseline = detectors[0].location - detectors[1].location 960 baseline = baseline / np.linalg.norm(baseline) 961 962 # get angle to baseline 963 angle = np.arccos(np.dot(cart, baseline)) 964 965 # get evenly spaced ring over north pole 966 longitudes = np.linspace(0, lal.TWOPI, n, endpoint=False) 967 latitudes = [lal.PI_2-angle]*len(longitudes) 968 # get axis and angle of rotation 969 north = np.array([0,0,1]) 970 axis = np.cross(north, baseline) 971 axis = axis / np.linalg.norm(axis) 972 angle = np.arccos(np.dot(north, baseline)) 973 R = _rotation(axis, angle) 974 975 # construct sky table 976 iso = SkyPositionTable() 977 for lon,lat in zip(longitudes, latitudes): 978 e = SkyPosition() 979 e.longitude = lon 980 e.latitude = lat 981 e.probability = None 982 e.solid_angle = None 983 e.system = 'geographic' 984 e.normalize() 985 e = e.rotate(R) 986 if gpstime: 987 e = GeographicToEquatorial(e, gpstime) 988 iso.append(e) 989 990 return iso

991

992 -def MaxTimeDelayLine(ifos, ra, dec, gpstime=None, n=3):

993 """ 994 Return the n-point SkyPositionTable describing the line perpendicular to 995 the line of constant time delay through the given ra and dec for the 996 2-tuple ifos. 997 """ 998 999 # construct sky point 1000 p = SkyPosition() 1001 p.longitude = ra 1002 p.latitude = dec 1003 p.system = 'equatorial' 1004 1005 # remove duplicate site references 1006 ifos = parse_sites(ifos) 1007 assert len(ifos)==2, "More than two sites in given list of detectors" 1008 1009 # get light travel time 1010 ltt = inject.light_travel_time(ifos[0], ifos[1]) 1011 # get number of points 1012 dt = ltt/n 1013 1014 return TwoSiteSkyGrid(ifos, gpstime, dt=dt, point=p, sky='full')

1015

1016 -def MaxTimeDelayLine3(ifo1, ifo2, ra, dec, gpstime=None, n=3):

1017 """ 1018 Alternative implementation of MaxTimeDelayLine. 1019 """ 1020 1021 if gpstime: 1022 gpstime = LIGOTimeGPS(gpstime) 1023 1024 p = SkyPosition() 1025 p.longitude = ra 1026 p.latitude = dec 1027 p.system = 'equatorial' 1028 p.probability = None 1029 p.solid_angle = None 1030 if gpstime: 1031 p = EquatorialToGeographic(p, gpstime) 1032 cart = SphericalToCartesian(p) 1033 1034 # get baseline 1035 detectors = [inject.cached_detector.get(inject.prefix_to_name[ifo])\ 1036 for ifo in [ifo1,ifo2]] 1037 baseline = detectors[0].location - detectors[1].location 1038 baseline = baseline / np.linalg.norm(baseline) 1039 1040 # get angular spacing 1041 dtheta = lal.TWOPI/n 1042 1043 # get axis and angle of rotation 1044 north = np.array([0,0,1]) 1045 axis = np.cross(cart, baseline) 1046 axis = axis / np.linalg.norm(axis) 1047 R = _rotation(axis, dtheta) 1048 1049 # set up list 1050 l = SkyPositionTable() 1051 # append central point 1052 l.append(p) 1053 1054 for i in xrange(1, n): 1055 l.append(l[i-1].rotate(R)) 1056 1057 if gpstime: 1058 for i,p in enumerate(l): 1059 l[i] = GeographicToEquatorial(p, gpstime) 1060 return l

1061 1062 # ============================================================================= 1063 # Helpers 1064 # ============================================================================= 1065

1066 -def SphericalToCartesian(skyPoint):

1067 """ 1068 Convert SkyPosition object into Cartesian 3-vector 1069 """ 1070 1071 p = np.zeros(3) 1072 phi = skyPoint.longitude 1073 theta = lal.PI_2 - skyPoint.latitude 1074 a = np.sin(phi) 1075 b = np.cos(phi) 1076 c = np.sin(theta) 1077 d = np.cos(theta) 1078 1079 p[0] = c*b 1080 p[1] = c*a 1081 p[2] = d 1082 1083 return p

1084

1085 -def CartesianToSpherical(x, system='equatorial'):

1086 """ 1087 Convert 3-tuple Cartesian sky position x to SkyPosition object in the 1088 given system 1089 """ 1090 1091 assert len(x)==3 1092 1093 p = SkyPosition() 1094 p.system = system 1095 p.longitude = np.arctan2(x[1], x[0]) 1096 p.latitude = lal.PI_2 - np.arccos(x[2]) 1097 p.probability = None 1098 p.solid_angle = None 1099 p.normalize() 1100 1101 return p

1102

1103 -def _rotation(axis, angle):

1104 """ 1105 Form 3x3 rotation matrix to rotate about a given 3-tuple axis by a given 1106 angle 1107 """ 1108 1109 R = np.zeros((3,3)) 1110 1111 ux,uy,uz = axis 1112 c = np.cos(angle) 1113 s = np.sin(angle) 1114 1115 R[0][0] = c + ux**2*(1-c) 1116 R[0][1] = ux*uy*(1-c) - uz*s 1117 R[0][2] = ux*uz*(1-c) + uy*s 1118 1119 R[1][0] = uy*ux*(1-c) + uz*s 1120 R[1][1] = c + uy**2*(1-c) 1121 R[1][2] = uy*uz*(1-c) - ux*s 1122 1123 R[2][0] = uz*ux*(1-c) - uy*s 1124 R[2][1] = uz*uy*(1-c) + ux*s 1125 R[2][2] = c + uz**2*(1-c) 1126 1127 return R

1128

1129 -def _rotation_euler(alpha, beta, gamma):

1130 """ 1131 Form rotation matrix from Euler angles. 1132 """ 1133 c = np.cos 1134 s = np.sin 1135 1136 # Define rotation matrix 1137 R = np.array([[c(alpha) * c(gamma) + s(alpha) * s(beta) * s(gamma), 1138 c(beta) * s(alpha), 1139 c(gamma) * s(alpha) * s(beta) - c(alpha) * s(gamma)], 1140 [c(alpha) * s(beta) * s(gamma) - c(gamma) * s(alpha), 1141 c(alpha) * c(beta), 1142 s(alpha) * s(gamma) + c(alpha) * c(gamma) * s(beta)], 1143 [c(beta) * s(gamma), 1144 -s(beta), 1145 c(beta) * c(gamma)]], 1146 dtype=float) 1147 1148 return R

1149

1150 -def parse_sites(ifos):

1151 """ 1152 Returns a new list of interferometers containing one only per site. 1153 I.e. this removes H2 if included. 1154 """ 1155 1156 # rebind ifos 1157 ifos = list(ifos) 1158 # remove duplicate site references 1159 ifos2 = [] 1160 for ifo in ifos: 1161 if len([i for i in ifos2 if i.startswith(ifo[0])]) == 0: 1162 ifos2.append(ifo) 1163 1164 return ifos2

1165

Source Code for Module pylal.sky