1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
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
40 _comment_regex = re.compile(r"\s*([#;].*)?\Z", re.DOTALL)
41 _delim_regex = re.compile('[,\s]')
42
43
44
45
46
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
61 rotated = table.new_from_template(self)
62 for row in self:
63 rotated.append(row.rotate(R))
64 return rotated
65
69
72
73
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
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
93 if n==0:
94 degenerate.append(False)
95 new.append(row)
96 continue
97 else:
98 degenerate.append(True)
99
100 for i in xrange(0,n):
101
102 if degenerate[i]:
103 continue
104
105 for j in xrange(0, len(timeDelays[n])):
106
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
122
123 __slots__ = SkyPositionTable.valid_columns.keys()
124
126 return iter(tuple((self.longitude, self.latitude)))
127
129 return "SkyPosition(%s, %s)" % (self.longitude, self.latitude)
130
132 return "(%s, %s)" % (self.longitude, self.latitude)
133
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
155 """
156 Normalize this SkyPosition to be within standard limits
157 [0 <= longitude < 2pi] and [-pi < latitude <= pi]
158 """
159
160
161
162
163
164 while self.longitude < 0:
165 self.longitude += lal.TWOPI
166 while self.longitude >= lal.TWOPI:
167 self.longitude -= lal.TWOPI
168
169
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
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
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
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
225
226
227
228 lal.LAL_ALPHAGAL = 3.366032942
229 lal.LAL_DELTAGAL = 0.473477302
230 lal.LAL_LGAL = 0.576
231
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
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
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
312 sinD = sina*sinP + cosa*cosA*cosP
313 sinH = cosa*sinA
314 cosH = sina*cosP - cosa*cosA*sinP
315
316
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
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
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
351 sina = sinD*sinP + cosD*cosP*cosH
352 sinA = cosD*sinH
353 cosA = sinD*cosP - cosD*sinP*cosH
354
355
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
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
376 gmst = np.fmod(date.XLALGreenwichSiderealTime(gpstime, 0),\
377 lal.TWOPI)
378
379
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
391 """
392 Convert the SkyPosition object input from the inherited 'equatorial'
393 system to to 'ecliptic'.
394 """
395
396
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
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
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
418 """
419 Convert the SkyPosition object input from the inherited 'equatorial'
420 system to to 'galactic'.
421 """
422
423
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
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
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
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
458 gmst = np.fmod(date.XLALGreenwichSiderealTime(gpstime, 0),\
459 lal.TWOPI)
460
461
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
473 """
474 Convert the SkyPosition object input from the inherited 'eliptic'
475 system to to 'equatorial'.
476 """
477
478
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
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
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
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
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
562 p = SkyPosition()
563 p.longitude = ra
564 p.latitude = dec
565
566
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
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
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
597 dalpha = ltt * np.sin(l)
598 if dalpha > alpha:
599 alpha = dalpha
600
601
602 angRes = 2 * dt/alpha
603
604
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
612
613
614
615 north = [0, 0, 1]
616 angle = np.arccos(np.dot(north, SphericalToCartesian(p)))
617
618
619
620 axis = np.cross(north, SphericalToCartesian(p))
621 axis = axis / np.linalg.norm(axis)
622
623
624 R = _rotation(axis, angle)
625 grid = grid.rotate(R)
626
627
628
629
630
631
632 kappa = (0.66*radius/sigma)**(-2)
633
634
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
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
656 n = max(1, int(np.ceil(lal.TWOPI\
657 * np.sin(theta)/resolution)))
658
659 dp = lal.TWOPI / n
660
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
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
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
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
695 ifos = parse_sites(ifos)
696 assert len(ifos)==2, "More than two sites in given list of detectors"
697
698
699 detectors = [inject.cached_detector.get(inject.prefix_to_name[ifo])\
700 for ifo in ifos]
701
702
703 ltt = inject.light_travel_time(ifos[0], ifos[1])
704
705
706 max = int(np.floor(ltt/dt))
707
708
709 baseline = detectors[1].location - detectors[0].location
710 baseline /= np.linalg.norm(baseline)
711
712
713
714
715
716
717
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
734 lineofnodes = np.cross([0,0,1], zaxis)
735 lineofnodes /= np.linalg.norm(lineofnodes)
736
737
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
743 R = _rotation_euler(alpha, beta, gamma)
744
745
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
755 t = i * dt
756
757 if abs(t) >= ltt: continue
758
759 e = SkyPosition()
760 e.system = 'geographic'
761
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
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
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
799 for i,ifo in enumerate(ifos):
800 detectors.append(inject.cached_detector.get(inject.prefix_to_name[ifo]))
801
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
808 angle = np.arccos(np.dot(baseline[0], baseline[1]))
809
810
811
812
813
814 tiling = tiling.lower()
815 t = np.zeros(len(baseline))
816 tdgrid = []
817
818
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
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
834 elif re.match('hex', tiling):
835 dm = (2*dx+dy)
836 dn = (-dx+dy)
837 dx = dm - dn
838
839 grid = []
840 orig = np.array([0,0])
841
842
843 while orig[1] >= -ny: orig -= dn
844
845
846 for y in xrange(-ny,ny+1):
847 orig[0] = orig[0] % dx[0]
848
849 minx = -nx
850 while minx % 3 != orig[0]: minx+=1
851
852 x = np.arange(minx, nx+1, dx[0])
853
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
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876 zaxis = baseline[0]
877
878
879 yaxis = np.cross(zaxis, baseline[1])
880 yaxis /= np.linalg.norm(yaxis)
881
882
883 xaxis = np.cross(yaxis, zaxis)
884 xaxis /= np.linalg.norm(xaxis)
885
886
887 north = np.array([0, 0, 1])
888 lineofnodes = np.cross(baseline[0], north)
889 lineofnodes /= np.linalg.norm(lineofnodes)
890
891
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
897 R = _rotation_euler(alpha, beta, gamma)
898
899
900
901
902
903 l = SkyPositionTable()
904
905 k = 0
906 for i in grid:
907 t = dt*i
908
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
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
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
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
963 angle = np.arccos(np.dot(cart, baseline))
964
965
966 longitudes = np.linspace(0, lal.TWOPI, n, endpoint=False)
967 latitudes = [lal.PI_2-angle]*len(longitudes)
968
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
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
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
1000 p = SkyPosition()
1001 p.longitude = ra
1002 p.latitude = dec
1003 p.system = 'equatorial'
1004
1005
1006 ifos = parse_sites(ifos)
1007 assert len(ifos)==2, "More than two sites in given list of detectors"
1008
1009
1010 ltt = inject.light_travel_time(ifos[0], ifos[1])
1011
1012 dt = ltt/n
1013
1014 return TwoSiteSkyGrid(ifos, gpstime, dt=dt, point=p, sky='full')
1015
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
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
1041 dtheta = lal.TWOPI/n
1042
1043
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
1050 l = SkyPositionTable()
1051
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
1064
1065
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
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
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
1130 """
1131 Form rotation matrix from Euler angles.
1132 """
1133 c = np.cos
1134 s = np.sin
1135
1136
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
1151 """
1152 Returns a new list of interferometers containing one only per site.
1153 I.e. this removes H2 if included.
1154 """
1155
1156
1157 ifos = list(ifos)
1158
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