Some formulas for burst network analysis

Andrea Viceré (avicere@ligo.caltech.edu)

Contents

This notebook can be read with MathReader or executed with Mathematica, at least version 4.0 (it uses the standard packages listed here)
The basic mathematics needed is the decomposition of rank-2 tensors in spin-weighed spherical harmonics. Some examples of these decomposition are given, maybe not interesting for the general reader: they show how to decompose any tensor in symmetric trace free tensors, and in particular the detector tensor.
The same operation is possible for the signal tensor: each detector is sensitive therefore to an appropriate product of components of detector and signal tensor, which define an extended beam pattern (extended, because it includes the orientation of the binary orbit and the wave polarization).
As the transformation of the spherical harmonics under rotations are easy to write down in terms of Gel'fand functions, these extended beam patterns are obtained by writing the detector tensor in a basis in which it is simple, and then rotating it to the wave basis in which the signal is simple. Being interested in writing the expressions in a reference (or network) frame, common to all detectors, one has that the final result contains the product of the two rotations needed to transform from the network frame to the detector frame, and from the network frame to the wave frame.
Using the positions and orientations of the different detectors, it is then possible to compute the Euler angles of the transformations which specify the detector frames; the network (or Earth, or fide frame) is assumed to be centered on the Earth, have Z along the north pole and X along the Greenwich meridian.
To visually cross check the transformations, please give a look at the global view of the various detectors and the frames attached to each of them.
There is another ingredient which we need to specify the network analysis: the relative sensitivity of the detectors. To this end we have collected the projected spectral densities of the different instruments: we plot also the resulting strain sensitivity. From the sensitivity curves it is immediate to compute the relative sensitivity to coalescing binary events (normalized to one of the detectors, Hanford 4K in our case), which enters the definition of the weighed extended beam patterns  together with the Euler angles specifying the detector frame.
Several possible groups of interferometers can be defined. The detection procedure, and in particular the weight to be applied to the signals extracted from each detector, depend on the set of interferometers chosen, as well as on the relative sensitivities. Having fixed these characteristics of the network, each detector produces in phase and in quadrature signals which need to be combined using a projection tensor, acting on the space of detectors and depending on the position in the sky. We plot examples of entries, for the particular case in which all the interferometers are included in the network. It is of course also needed to introduce proper delays, also depending on the source direction.
A useful quantity to plot is then the energy available to the network: it is just the sum of the squares of the extended beam patterns for each detector. We define this energy also as an average over source orientation and polarization, as a function of the source direction, we take the square root to have a quantity proportional to the SNR, and we plot it for several different network configurations.
When comparing the results given here with those cited in the paper, there are essentially two differences
a) the spectral noise densities have been updated;
b) here the network frame coincides with the Earth frame, and with the fide frame.

Defined functions

By having just the "initialization" cells evaluated, several functions result defined: type ?name for help on each of them.
Ysh defines the spherical harmonics, as a basis of rotations in 3D
Y defines the spin weighed spherical harmonics
STFDecompose allows to express a 3x3 matrix as a sum over spin weighed spherical harmonics. 
T is the Gel'fand function which allows to write down the representation of rotations on the basis of spherical harmonic functions.
EBPCoefficientD is just the D coefficient (see the above cited paper) which embodies that part of the extended beam-pattern which does not depend on the source orientation and polarization.
EBP is instead the true extended beam pattern which depends also on angles defining source orientation and polarization.
latitude, longitude define the position of detectors on the Earth, depending on mnemonics for the detector indices
armXazimuth, armYazimuth define analogously the orientation of the detector arms.
positionEarth returns a versor pointing to a detector location, expressed in the Earth frame.
euler returns the Euler angles needed to rotate coordinates in the Earth frame into coordinates in the detector frame.
fSeism, Snoise define respectively the lower frequency cutoff (the seismic wall) and the noise spectrum (one-sided PSD) of each interferometer.
gAbs, g define the absolute and relative sensitivity of each detector to a coalescing binary signal.
DhatP2, DhatM2 define the components of a network vector having g[itf]*D[itf] as components, and normalized.
Glow, Ghigh express the metric in the helicity subspace relative to some direction, with lower or upper indices.
p is the projection tensor which allows to combine signals from different detectors.

Needed packages

This package contains the rotation matrices and is mandatory

Needs["Geometry`Rotations`"];

These packages are needed only to do the plots

Needs["Graphics`Graphics`"];
Needs["Graphics`Legend`"];
Needs["Graphics`ParametricPlot3D`"];
Needs["Graphics`Shapes`"];

Definition of some special functions

Spherical harmonics

Here the spherical harmonics are defined as sums over the tensors obtained by products (Subscript[n, i]) (Subscript[n, j]) of components of a versor n = {(Subscript[n, 1]),(Subscript[n, 2]),(Subscript[n, 3])} pointing toward a direction specified by angles φ, ϑ.

Clear[Ysh];
\(\(Ysh::usage = 
\*""Ysh[l,m] returns the spherical harmonic with quantum numbers l,m, defined in terms of components of a three dimensional versor n={(Subscript[n, 1]), (Subscript[n, 2]), (Subscript[n, 3])}"";\)\)
\(\(Ysh[l_, \ 
          m_Integer?\((#\  ≥ \ 
                  0\  &)\)]\  := 
\((\(-1\))\)^
            m\ Sqrt[\((2*l + 1)\)/\((4\ π)\)*
              Factorial[l - m]/Factorial[l + m]]*\((n\_1\  + \ I\ n\_2)\)^
            m\ *\ 
Module[{j, \ 
              ris}, 
ris\  = \ 
              Expand[Sum[\((\(-1\))\)^
                        j/\((2^l\ Factorial[j]\ Factorial[l - j])\)*
                    Factorial[2\ l - 2\ j]/Factorial[l - m - 2\ j]*
                    n\_3^\((l - m - 2*j)\), {j, 0, 
                    Floor[\((l - m)\)/2]}]]; 
If[
              m == 0, \ 
Expand[
                If[Head[ris]\  === \ Plus, 
                  Map[If[FreeQ[#, \ 
                          n\__], \ #*\((n\_1^2\  + \ n\_2^2\  + \ 
                              n\_3^2)\), #] &, \ ris], 
If[
                    FreeQ[ris, n\__], 
                    ris*\((n\_1^2\  + \ n\_2^2\  + \ n\_3^2)\), \ ris]]], 
              ris]];\)
\)
Ysh[l_,  
          m_Integer?(#  <  0  &)]  :=  (-1)^m *
            Ysh[l,  -m] /. {Complex[a_,  b_]  :>  
              Complex[a,  -b]};

Spin weighed spherical harmonics

They are just the coefficients of the expansion of the spherical  harmonics in terms of the (Subscript[n, i])(Subscript[n, j]) products.

Clear[Y];
\(\(Y::usage = 
\*""Y[2,m] returns the 3x3 matrix defining the (m^th) spin weighed spherical harmonic"";\)\)
Y[2,  m_]  :=  
        Module[{k1, k2}, 
          Table[If[k1 == k2, 1, 1/2]*
              Coefficient[Expand[Ysh[2, m]], (Subscript[n, k1])*(Subscript[n, k2])], {k1, 1, 3}, {k2, 
              1, 3}]];

Decomposition of a tensor in spin weighed spherical harmonics

Clear[STFDecompose];
STFDecompose::usage = 
"STFDecompose[matrix] \
perform the decomposition of a 3x3 matrix into spin weighed spherical \
harmonics (displayed as Hold[Y[2,m]] to prevent their expansion)";
STFDecompose[a_?MatrixQ] := 
        Module[{m}, 
          Sum[Tr[a . Y[2, m]]/(15/(8 π))*
                Conjugate[holdY[2, m]], {m, -2, 2, 1}] /. 
            holdY[x__] :> Hold[Y[x]]];

Gel'fand function T

It is the representation of rotations over the basis of spherical harmonic functions

Clear[T];
T::usage = 
"T[l, m, n, α, β, \
γ] is the Gel'fand function: at fixed order l, it is a l x l matrix in \
the indices m, n, which represent the rotation group over the basis of \
spherical harmonic functions. The α, β, γ arguments \
represent the Euler angles of the 3D rotation.";
T[l_,  m_,  
            n_,  α_,  β_,  γ_]  :=  
Exp[-I m γ] Exp[-I n α]*
            Module[{ris, 
                mu}, 
ris = (-1)^(l - m)*
                  I^(n - m)/(2^l*Factorial[l - m])*
                  Sqrt[Factorial[l - m]*
                      Factorial[
                          l + n]/(Factorial[l + m]*
                            Factorial[l - n])]*(1 - 
                        mu)^(-(n - m)/2)*(1 + 
                        mu)^(-(n + m)/2)*
                  D[(1 - mu)^(l - m)*(1 + mu)^(l + m), {mu, 
                      l - n}]; 
Factor[Simplify[
                      ris] //. {Sqrt[1 - mu^2] :> Sin[β], 
                      Sqrt[1 - mu] :> 
                        1/Sqrt[1 + mu]*
                          Sin[β], (1 - mu)^(3/2) :> 
                        1/Sqrt[1 + mu]*
                          Sin[β]*(1 - mu)} /. {mu :> 
                      Cos[β]}]];;

Examples

Detector tensor in the detector frame

n1  =  1/Sqrt[2] {1, -1, 0};  
    n2  =  1/Sqrt[2] {1, 1, 0};
d  =  
        Outer[Times,  n1,  n1]  -  
          Outer[Times,  n2,  n2];

Decomposition of the tensor d in spin-weighted spherical harmonics

d

STFDecompose[d]

Release[%]

Polarization tensor in the wave frame

eX  =  {1, 0, 0};  eY  =  {0, 1, 0};
eL  =  1/2 Outer[Times,  eX  +  I eY,  eX  +  I eY];  
    eR  =  Conjugate[eL];

STF decomposition of the polarization tensor

eR

STFDecompose[eR]

Release[%]

Extended beam pattern definitions

D coefficients of the extended beam pattern

These are the cofficients of the expansion of the detector tensor in the basis of symmetric trace free tensors. The angles φ, ϑ, ψ  are the Euler angles which allow to rotate a vector in the network frame into a vector in the wave frame. Notice that φ and ϑ alone allow to define  the direction of propagation of the wave. 
On the other hand α,  β, γ are the Euler angles which allow to rotate from the network frame to the detector frame.

Clear[EBPCoefficientD];
\(\(EBPCoefficientD::usage = 
\*""EBPCoefficientD[n, φ, ϑ, ψ, α, β, γ] returns the (n^th) coefficient of the expansion of the detector tensor in the basis of symmetric trace-free tensors. The φ,ϑ,ψ are the Euler angles allowing to rotate a vector in the network frame into a vector in the wave frame; analogously, α,β,γ are Euler angles allowing to rotate from the network frame to the detector frame."";\)\)
EBPCoefficientD[
          n_,  φ_, ϑ_, ψ_, α_, β_, γ_]  :=  -I Module[{s},  
            Sum[T[2, n, 
                  s, φ, ϑ, ψ]*((T[2, 2, 
                          s, α, β, γ]  -  
                        T[2, -2, s, α, β, γ]) /. 
                    Complex[r_, i_] :> 
                      Complex[r, -i]), {s, -2, 2, 
                1}]];
Clear[EBPCoefficientDP2];
EBPCoefficientDP2::usage = 
"EBPCoefficientDP2[\
φ,ϑ,ψ,α,β,γ] is a simplified \
expression for EBPCoefficientD[2,φ,ϑ,ψ,α,\
β,γ]. It can also be used the form \
EBPCoefficientDP2[φ,ϑ,ψ,itf] where itf is a \
mnemonic for the detector name; in this case the coefficient includes also \
the sensitivity (to coalescing binary signals) relative to the LIGO 4km \
detector";
EBPCoefficientDP2[φ_, ϑ_, ψ_, α_, β_, γ_] = 
        Expand[Expand[
              Expand[EBPCoefficientD[
                            2, φ, ϑ, ψ, α, β, γ] //. (-1 + Cos[a_]) (1 + 
                                  Cos[a_]) :>  -Sin[
                                  a]^2 //. 
                        Exp[a_] :>  
                          Cos[Expand[-I a]] + 
                            I Sin[Expand[-I a]] //. {Cos[
                            a_ + b_]  :>  
                          Cos[a] Cos[b] - Sin[a] Sin[b],  
                        Sin[a_ + b_]  :>  
                          Sin[a] Cos[b] + Cos[a] Sin[b]} //. {Cos[
                          2 a_]  :>  1  -  2 Sin[a]^2,  
                      Sin[2 a_]  :>  
                        2 Sin[a] Cos[a]}] //. {Cos[a_]^2  :>  
                    1  -  Sin[a]^2}] /. {Cos[a_]^2  :>  
                1  -  Sin[a]^2}];
Clear[EBPCoefficientDM2];
EBPCoefficientDM2::usage = "EBPCoefficientDM2[φ,\
ϑ,ψ,α,β,γ] is a simplified expression \
for EBPCoefficientD[-2,φ,ϑ,ψ,α,β,\
γ]";
EBPCoefficientDM2[φ_, ϑ_, ψ_, α_, β_, γ_] = 
        Expand[Expand[
              Expand[EBPCoefficientD[-2, φ, ϑ, ψ, α, β, γ] //. (-1 + 
                                  Cos[a_]) (1 + 
                                  Cos[a_]) :>  -Sin[
                                  a]^2 //. 
                        Exp[a_] :>  
                          Cos[Expand[-I a]] + 
                            I Sin[Expand[-I a]] //. {Cos[
                            a_ + b_]  :>  
                          Cos[a] Cos[b] - Sin[a] Sin[b],  
                        Sin[a_ + b_]  :>  
                          Sin[a] Cos[b] + Cos[a] Sin[b]} //. {Cos[
                          2 a_]  :>  1  -  2 Sin[a]^2,  
                      Sin[2 a_]  :>  
                        2 Sin[a] Cos[a]}] //. {Cos[a_]^2  :>  
                    1  -  Sin[a]^2}] /. {Cos[a_]^2  :>  
                1  -  Sin[a]^2}];

Specialized for the interferometer

EBP[φ_, ϑ_, ψ_, ε_, 
          itf_Integer] := EBP[φ, ϑ, ψ, ε, ##] & @@ euler[itf];

As with EBP, they are specialized for the interferometer

EBPCoefficientDP2[φ_, ϑ_, ψ_, 
          itf_Integer] := EBPCoefficientDP2[φ, ϑ, ψ, ##] & @@ euler[itf];
EBPCoefficientDM2[φ_, ϑ_, ψ_, 
          itf_Integer] := EBPCoefficientDM2[φ, ϑ, ψ, ##] & @@ euler[itf];

Complex antenna patterns of the detectors

The real part gives the sensitivity to a wave polarized "+" in the wave frame, and the imaginary part the same quantity for a wave polarized "x" in the wave frame

Clear[F];
F::usage = "F[ϕ,θ,α,β,γ] returns a \
complex number. The real part is the sensitivity to a wave with a \
polarization + in the wave frame, and the imaginary part the sensitivity to a \
wave with a polarization x in the wave frame; they depend on the direction of \
propagation of the wave, specified by the angles ϕ,θ,0 of \
orientation of the wave frame with respect to the network frame, and on the \
orientation of the detector frame, dictated by the angles α,β,\
γ.Also the form F[ϕ,θ,itf] can be used, where itf is one \
of GEO,VIRGO,TAMA,L2K.L4K,H4K; in this case, a factor Sin[armOpening[itf]] is \
included.";
F[ϕ_, θ_, α_, β_, γ_] := 
        EBPCoefficientDM2[ϕ, θ, 
          0, α, β, γ];
F[ϕ_, θ_, itf_] := 
        Sin[armOpening[itf]]*
          EBPCoefficientDM2[ϕ, θ, 0, itf];

Interferometer locations and orientations

Data taken from the paper "Gravitational wave detector sites", B.Allen, gr-qc/9607075, 1997,
formatted to comply with the specifications adopted in the Frame format (LIGO-T970130-D-E and VIR-SPE-LAP-5400-102).
More precisely, the longitude is given east of Greenwich and the arm azimuths are measured from east.
Data are input in degrees and converted in radians.

Initials and mnemonics for the various detectors

Clear[namesITF];
\(\(namesITF = {"G", "H", "H", "L", "T", "V"};\)\)
Clear[GEO, H2K, H4K, L4K, TAMA, VIRGO];
\(\(Map[\((\(Evaluate[#[\([1]\)]]\)::"usage" = 
              StringJoin[
                ToString[#[\([1]\)]], " is the index referring to the ", 
                ToString[#[\([2]\)]], " interferometric observatory."])\) &, {{GEO, "German-British 600m"}, {H2K, "Hanford 2km LIGO"}, {H4K, "Hanford 4km LIGO"}, {L4K, "Livingston 4km LIGO"}, {TAMA, "Japanese 300m"}, {VIRGO, "French-Italian 3km"}}];\)\)
GEO  = 1;
H2K = 2;
H4K = 3;
L4K = 4;
TAMA = 5;
VIRGO = 6;

Grouping of the interferometers

We define several lists with different combinations of interferometers

itfAll = {GEO, H4K, H2K, L4K, TAMA, 
          VIRGO};
itfLIGO = {H4K, H2K, L4K};
itfLIGOandVIRGO = 
        Append[itfLIGO, VIRGO];
itfLIGOandGEO = Append[itfLIGO, GEO];
itfGEOandVIRGO = {GEO, VIRGO};

Latitudes and longitudes

Latitude north of the equator, longitude east of Greenwich, measured in radians

deg2rad = π/180;

Clear[latitude];
latitude::usage = "latitude[detector] gives the latitude of the \
corner station, in radians North of the equator";
latitude[GEO] = (52.25)*deg2rad;
latitude[
          H2K] = latitude[H4K] = (46.45)*
            deg2rad;
latitude[L4K] = (30.56)*deg2rad;
latitude[TAMA] = (35.68)*deg2rad;
latitude[VIRGO] = (43.63)*
          deg2rad;

Clear[longitude];
longitude::usage = "longitude[detector] gives the longitude of the \
corner station, in radians East of Greenwich";
longitude[GEO] = (9.81)*deg2rad;
longitude[
          H2K] = longitude[H4K] = (-119.41)*
            deg2rad;
longitude[L4K] = (-90.77)*
          deg2rad;
longitude[TAMA] = (139.54)*deg2rad;
longitude[VIRGO] = (10.5)*deg2rad;

Orientations of the arms

measured counter-clock-wise from east

Clear[armXazimuth];
\(\(armXazimuth::usage = \*""armXazimuth[detector] gives the angle in radians between the X arm of the detector and the local East, measured counter-clock-wise."";\)\)
armXazimuth[GEO] = 
        Mod[291.61 + 90, 360]*deg2rad;
armXazimuth[
          H2K] = armXazimuth[H4K] = (36.8 + 90)*
            deg2rad;
armXazimuth[L4K] = (108.0 + 90)*
          deg2rad;
armXazimuth[TAMA] = (90.0 + 90)*
          deg2rad;
armXazimuth[VIRGO] = 
        Mod[341.5 + 90, 360]*
          deg2rad;

Clear[armYazimuth];
armYazimuth::usage = "armYazimuth[detector] gives the angle in \
radians between the Y arm of the detector and the local East, measured \
counter-clock-wise.";
armYazimuth[GEO] = (25.94 + 90)*
          deg2rad;
armYazimuth[
          H2K] = armYazimuth[H4K] = (126.8 + 90)*
            deg2rad;
armYazimuth[L4K] = (198.0 + 90)*
          deg2rad;
armYazimuth[TAMA] = (180.0 + 90)*
          deg2rad;
armYazimuth[VIRGO] = (71.5 + 90)*
          deg2rad;

Angle between the arms

Clear[armOpening];
armOpening::usage = "armOpening[detector] gives the angle in \
radians between the X and Y arm of the detector";
armOpening[detector_] = 
        armYazimuth[detector] - armXazimuth[detector];

Network and detector frames

Earth (or network or fide) frame

This frame is defined by having the X axis joining Earth's center and the equatorial point having 0 longitude,
the Z axis pointing to the true north. The Y axis follows from being a right handed system. For simplicity, we shall use this as the reference (or fide) frame.

Clear[positionEarth];
\(\(positionEarth::usage = "positionEarth[detector] returns a unit versor, in the Earth frame, locating the corner station of the detector";\)\)
positionEarth[itf_] := {Cos[latitude[itf]]*Cos[longitude[itf]], 
          Cos[latitude[itf]]*Sin[longitude[itf]], 
          Sin[latitude[itf]]};

Detector frames

These frames are defined by having the interferometer arms in the X,Y plane, and the Z axis pointing outward from the Earth. Further, the X axis is chosen to bisect the angle between the detector's arms. The important quantities to know are the Euler angles needed to transform coordinates in the fide to coordinates in the detector frame.

Detector frame axes definition, in the Earth frame

We compute here just what are the coordinates of the axes of the detector frame, expressed in the Earth frame

normalize[
          v_?VectorQ] := 
        v/Sqrt[v . v];

  
dfArmsAxes[itf_] := 
        Module[{localEast, localNorth, dfAxisX, dfAxisY, dfAxisZ, armX, 
            armY}, 
  
dfAxisZ = 
            normalize[positionEarth[itf]]; 
  
localEast = 
            normalize[
              Cross[dfAxisZ, dfAxisZ*{1, 1, 0}]]; 
  
localNorth = 
            normalize[Cross[dfAxisZ, localEast]]; 
  
armX = 
            normalize[
              Cos[armXazimuth[itf]]*localEast + 
                Sin[armXazimuth[itf]]*localNorth]; 
armY = 
            normalize[
              Cos[armYazimuth[itf]]*localEast + 
                Sin[armYazimuth[itf]]*localNorth]; 
  
dfAxisX = 
            normalize[armX + armY]; 
  
dfAxisY = 
            normalize[Cross[dfAxisZ, dfAxisX]]; 
{armX, 
            armY, dfAxisX, dfAxisY, dfAxisZ}];

Euler angles to go from one frame to another

These Euler angles are needed to trasform coordinates expressed in the Fide frame into coordinates expressed in the detector frame. We have chosen here as fide frame the Earth frame.

(* first we compute the Euler angles to do the opposite,  
      translating coordinates in the detector frame into coordinates in the Fide frame *) 

eulerInverse[
          itf_] := 
        Module[{rotation, ϕ, ϑ, ψ}, 
 (* this is the rotation matrix which rotates coordinates in the
detector frame into coordinates in the fide frame *) 
rotation = 
            Take[dfArmsAxes[itf], {3, 5}]; 
 (* 
            setting it equal to a rotation matrix we find the angles *) 
ϑ = 
            ArcCos[rotation[[3, 3]]]; 
ϕ = 
            If[Sin[ϑ] ≠ 
                0.0, 
             ArcTan[-rotation[[3, 2]]/Sin[ϑ], 
                rotation[[3, 1]]/
                  Sin[ϑ]], 
             ArcCos[rotation[[1, 1]]]]; 
ψ = 
            If[Sin[ϑ] ≠  
                0.0, 
             ArcTan[
                rotation[[2, 3]]/Sin[ϑ], 
                rotation[[1, 3]]/
                  Sin[ϑ]], 
             0.0]; 
 (* 
            notice that if ϑ == 
              0 only the ϕ + ψ combination is relevant, 
and conventionally we set ψ = 
              0 in this case *) 
 {ϕ, ϑ, ψ}];

 (* 
      inverse,  
      to transform coordinates in the Fide frame into
coordinates in the detector frame *) 
Clear[euler];
euler::usage = "euler[detector] gives the Euler angles which allow \
to transform coordinates in the network (or Earth) frame into coordinates in \
the detector frame.";
euler[itf_]  := -Reverse[eulerInverse[itf]];

Graphical representation

We show the detector frames on the surface of a sphere, just as a check

\( (*\ 
      the\ axes\ of\ the\ network\ frame\ *) 
\(networkFrameTriad\  = 
        Graphics3D[{Thickness[0.01], RGBColor[1, 0, 0], 
            Line[{{0, 0, 0}, {1.5, 0, 0}}], 
            Text["X", {1.5, 0, \(-0.25\)}, 
              TextStyle -> \ {FontFamily -> \ "Times-Roman", 
                  FontSize -> \ 17}], RGBColor[0, 1, 0], 
            Line[{{0, 0, 0}, {0, 1.5, 0}}], 
            Text["Y", {0, 1.5, \(-0.25\)}, 
              TextStyle -> \ {FontFamily -> \ "Times-Roman", 
                  FontSize -> \ 17}], RGBColor[0, 0, 1], 
            Line[{{0, 0, 0}, {0, 0, 1.5}}], 
            Text["Z", {0, 0, 1.75}, 
              TextStyle -> \ {FontFamily -> \ "Times-Roman", 
                  FontSize -> \ 
                    17}]}];\)

 (*\ 
      a\ triad\ of\ axes, \ 
      the\ longest\ is\ the\ Z\ one\ *) \)
AxesTriad  = 
        Graphics3D[{Thickness[0.01], RGBColor[1, 0, 0], 
            Line[{{0, 0, 0}, {0.25, 0, 0}}], RGBColor[0, 1, 0], 
            Line[{{0, 0, 0}, {0, 0.25, 0}}], RGBColor[0, 0, 1], 
            Line[{{0, 0, 0}, {0, 0, 
                  0.25}}]}];

  
ArmsITF[itf_] := 
        Module[{arms}, 
arms = 
            Take[dfArmsAxes[itf], {1, 2}]/
              4. ; 
TranslateShape[
            Graphics3D[{Thickness[0.01], Line[{{0, 0, 0}, arms[[1]]}], 
                Line[{{0, 0, 0}, arms[[2]]}]}],  
            positionEarth[
              itf]]
];
AxesITF[itf_] := 
        TranslateShape[RotateShape[AxesTriad, ##]  &  @@  euler[itf], 
          positionEarth[itf]];
NameITF[itf_]  := 
        Graphics3D[
          Text[namesITF[[itf]], positionEarth[itf], 
            TextStyle -> {FontFamily -> Times, FontSize -> 16, 
                FontColor -> 
                  RGBColor[0, 1, 
                    0]}]];

  
earth = 
        Graphics3D[Sphere[1.0]];

  
detectorGraphics[itf_] := 
        Sequence[NameITF[itf], ArmsITF[itf], 
          AxesITF[itf]];
networkGraphics = 
        Sequence[detectorGraphics[GEO], detectorGraphics[L4K], 
          detectorGraphics[H4K], detectorGraphics[VIRGO], 
          detectorGraphics[TAMA]];

View from above the United States

The black lines represent the detector arms, while the colored lines represent the axes of the detector frame.
The red one is X and bisects the two arms, the blu is Z (vertical), the green is Y in a right-handed axes triad.
We also show the triad of axes of the Earth frame, centered on the Earth.

Show[networkGraphics, earth, networkFrameTriad, 
      ViewPoint -> {0,  -1,  1}]

Export["networkAboveUS.eps", 
        Show[networkGraphics, earth, networkFrameTriad, 
          DisplayFunction -> Identity, 
          ViewPoint -> {0, -1, 1}], "epsi", 
        ImageResolution -> 160];

View from above Europe

Show[networkGraphics, earth, networkFrameTriad, 
      ViewPoint -> {1, 0, 1}]

Export["networkAboveEU.eps", 
        Show[networkGraphics, earth, networkFrameTriad, 
          DisplayFunction -> Identity, 
          ViewPoint -> {1, 0, 1}], "epsi", 
        ImageResolution -> 160];

Noise spectra for the different interferometers (simplified)

Here we collect the projected sensitivities of the different instruments, disregarding all the resonances and stopping at the "seismic noise" wall. The "model" includes only the pendular and internal thermal noise, and the shot noise. 
The spectrum is measured in units of 1/Hz, and the seismic wall is chosen as the frequency at which the thermal and seismic noises cross.

fSeism::usage = "fSeism[detector] returns the frequency at which \
the seismic noise wall is presumed to limit the detector sensitivity";
Snoise::usage = "Snoise[detector][f] returns the one-sided spectral \
density of the strain noise presumed to limit the detector sensitivity at \
frequency f";

GEO

The GEO sensitivity curve depends on the choice of the tuning for the signal recycling, which affects the shot noise.
Presuming that a wideband operation gives a better sensitivity to coalescing binaries signal, we have chosen to remove the signal recycling altogether in these examples
Source: http://www.geo600.uni-hannover.de/geo600/project/wideband.html

fSeism[GEO] = 50;
fNyquist[GEO] = 8192;
Snoise[GEO][f_] = (4.1*10^(-36))/f^5 + 
          9. *10^(-43)/f + 
          10^(-44)*(1 + (f/577)^2);

LIGO I

The LIGO I sensitivity considered here is obtained by a simple fit to the curves computed for the 2K and 4K interferometers by Kent Blackburn, and accessible through the ASIS site http://www.lsc-group.phys.uwm.edu/lscasis, following the documents link.
Note that the different length of the FP cavities in the 2K and 4K interferometers leads to a different knee frequency.
Source: http: // www . ligo . caltech . 
                edu/~kent/ASIS_NM/noise_models . html

fSeism[H4K  |  L4K] = 40;
fNyquist[H2K | H4K | L4K] = 8192;
Snoise[L4K  |  H4K][f_] = 
        5.6*10^(-36)/f^5 + 3.9*10^(-44)/f + 
          1.1*10^(-46)*(1 + (f/
                      83)^2);
fSeism[H2K] = 40;
Snoise[H2K][f_] = 
        2.1*10^(-35)/f^5 + 2.25*10^(-43)/f + 
          4.35*10^(-46)*(1 + (f/182)^2);

TAMA

Source: http://tamago.mtk.nao.ac.jp/tama/recom/recom3/index.html

fSeism[TAMA] = 50;
fNyquist[TAMA] = 10^4;
Snoise[TAMA][f_] = (6.6*10^(-31))/f^5 + 
          3.2*10^(-40)/f + 
          1.78*10^(-42)*(1 + (f/500)^2);

VIRGO

Source: http://www.virgo.infn.it/senscurve

fSeism[VIRGO] = 4. ;
fNyquist[VIRGO] = 10^4;
Snoise[VIRGO][f_]  = (9*10^(-37))/f^5 + 
          4.5*10^(-43)/f + 
          3.24*10^(-46)*(1 + (f/500)^2);

Plot of the sensitivities

This plot does not include the seismic wall

\(strainPlot = 
      ShowLegend[
        LogLogPlot[{Sqrt[\(Snoise[GEO]\)[f]], Sqrt[\(Snoise[H4K]\)[f]], \ 
            Sqrt[\(Snoise[H2K]\)[f]], Sqrt[\(Snoise[TAMA]\)[f]], 
            Sqrt[\(Snoise[VIRGO]\)[f]]}, {f, fSeism[VIRGO], 10^4}, 
          AxesLabel -> {Hz, ""}, 
PlotRange -> {{fSeism[VIRGO], 10^4}, {10^\(-23\), 
                10^\(-18\)}}, 
TextStyle -> {FontFamily -> "Times-Roman", 
              FontSize -> 
                17}, 
PlotStyle -> {{Dashing[{0.05, 
                    0.05}], Thickness[0.005], 
                RGBColor[1, 0, 0]}, {Dashing[{0.075, 0.025}], 
                Thickness[0.005], 
                RGBColor[0, 1, 0]}, {Dashing[{0.008, 0.012, 0.05, 0.012}], 
                Thickness[0.005], RGBColor[0, 0, 1]}, {Dashing[{0.01, 0.01}], 
                Thickness[0.005], RGBColor[1, 0.5, 0]}, {Thickness[0.005], 
                GrayLevel[0.25]}}, 
          GridLines -> {{Automatic, GrayLevel[0.75]}, {Automatic, 
                GrayLevel[0.75]}}, 
          DisplayFunction -> 
            Identity], 
{{{Graphics[{Dashing[{0.25, 0.25}], 
                  RGBColor[1, 0, 0], Thickness[0.06], 
                  Line[{{0, 0}, {1, 0}}]}], 
              FontForm["GEO600", {"Times-Roman", 
                  17}]}, {Graphics[{Dashing[{0.6, 0.2}], RGBColor[0, 1, 0], 
                  Thickness[0.06], Line[{{0, 0}, {1, 0}}]}], 
              FontForm["H4K, L4K", {"Times-Roman", 
                  17}]}, {Graphics[{Dashing[10*{0.004, 0.016, 0.05, 0.016}], 
                  RGBColor[0, 0, 1], Thickness[0.06], 
                  Line[{{0, 0}, {1, 0}}]}], 
              FontForm["H2K", {"Times-Roman", 
                  17}]}, {Graphics[{Dashing[{0.15, 0.25}*0.75], 
                  RGBColor[1, 0.5, 0], Thickness[0.06], 
                  Line[{{0, 0}, {1, 0}}]}], 
              FontForm["TAMA300", {"Times-Roman", 
                  17}]}, {Graphics[{GrayLevel[0.25], Thickness[0.06], 
                  Line[{{0, 0}, {1, 0}}]}], 
              FontForm["VIRGO", {"Times-Roman", 17}]}}, 
          LegendPosition -> {0.0, 0.13}, LegendSize -> {0.6, 0.4}, 
          LegendSpacing -> 0, LegendTextSpace -> 1.1, 
          LegendShadow -> {0.01, \(-0.01\)}}]\)

Export["sens.eps", strainPlot, "epsi", 
        ImageResolution -> 160, ImageSize -> 768, 
        ImageRotated -> True];

Inverse RMS of the different ITFs

This quantity appears in the calculation of the sensitivity of the detectors to a delta signal

Clear[rmsInverse];
\(\(rmsInverse::usage = \*""rmsInverse[detector] returns simply (∫_(fSeism[detector])^(fNyquist[detector])) (1/(Snoise[detector]))[f] \[DifferentialD]f"";\)
\)
rmsInverse[itf_]  :=  rmsInverse[itf] = 
          Module[{f}, 
NIntegrate[
              1/Snoise[itf][f], {f, fSeism[itf], 
                fNyquist[itf]}]];

Response to a δ signal

This is the response of the detector to a δ signal having certain amplitudes \(TraditionalForm\`\(\(A\_+\)\(,\)\)\) (Subscript[A, ×]) in the wave frame

deltaResponse::usage = "deltaResponse[{Aplus_,Atimes_},ϕ,\
θ,itf] gives the response of the interferometer itf to a δ \
function signal, with amplitudes Aplus and Atimes in the wave frame, \
propagating along the directions specified by the Euler angles \
ϕ,θ,0 which dictate the rotation of the wave frame with respect \
to the network frame";
deltaResponse[{Aplus_, Atimes_}, ϕ_, θ_, itf_] := 
        Re[(Aplus - I Atimes)*F[ϕ, θ, itf]]*
          rmsInverse[itf];

deltaResponse[{10^-21, 0}, 0.1, 0.5, GEO]

Network signal y

\(\(yNetwork::usage = \*""yNetwork[{Aplus,Atimes},ϕ,θ,itfGroup] gives a complex number, containing the response of the network, specified by the list itfGroup, to a  δ function signal, with amplitudes Aplus and Atimes in the wave frame, propagating along the directions specified by the Euler angles ϕ,θ,0 which dictate the rotation of the wave frame with respect to the network frame. It computes (∑_itf)F[ϕ,θ,itf]*deltaResponse[{Aplus,Atimes},ϕ,θ,itf], hence the real part corresponds to the + polarization, the imaginary part to the x polarization."";\)\)
yNetwork[{Aplus_, Atimes_}, ϕ_, θ_,  itfGroup_List] := 
        Plus @@ ((F[ϕ, θ, #]*
                    deltaResponse[{Aplus, 
                        Atimes}, ϕ, θ, #]) &  /@  
              itfGroup);

Matrix Θ

This is the matrix which specifies the inner product: its inverse, sandwiched with yNetwork, gives the likelyhood

Clear[Θ];
Θ[ϕ_?NumberQ, θ_?NumberQ, 
          itfGroup_List] := 
Plus  @@ ((Outer[
                      Times, #, #] &  /@  ((F[ϕ, θ, #] & /@  itfGroup) /. 
                      Complex[a_, b_] :>  {a, 
                          b}))* (rmsInverse[#] &  /@  
                  itfGroup));

Likelyhood

The likelyhood here expressed is valid only for δ function signals.

Clear[L]; 
    L::usage = "L[{Aplus,Atimes},ϕ,θ,itfGroup] is the \
log-likelyhood for a signal with polarizations Aplus and Atimes in the wave \
frame, specified by Euler angles ϕ,θ, for a group of \
interferometers specified by the list itfGroup. Analogously L[A,ϕ,\
θ,itfGroup] gives the same quantity, but averaged over the possible \
polarizations, with a scale amplitude A.";
L[{Aplus_?NumberQ, Atimes_?NumberQ}, ϕ_, θ_, 
          itfGroup_List] :=  
        Module[{y}, 
          y = yNetwork[{Aplus, Atimes}, ϕ, θ, itfGroup] /. 
              Complex[a_, b_] :> {a, b};  
          y . Inverse[Θ[ϕ, θ, itfGroup]] . 
            y];
L[A_?NumberQ, ϕ_, θ_, itfGroup_List] :=  
      Module[{yPlus, yTimes, iTheta}, 
        yPlus = yNetwork[{A, 0}, ϕ, θ, itfGroup] /. 
            Complex[a_, b_] :>  {a, b}; 
        yTimes = 
          yNetwork[{0, A}, ϕ, θ, itfGroup] /. 
            Complex[a_, b_] :>  {a, 
                b}; 
iTheta = 
          Inverse[Θ[ϕ, θ, 
              itfGroup]]; 
1/
            2 (yPlus . iTheta . yPlus + 
              yTimes . iTheta . yTimes)]

Plots of the sensitivity to unit length burst for the different combinations

The viewpoint is either automatic, above US or above the EU; we include a factor 1/Sqrt[2 * (2 burstLength)] = 1/2

plotFrameTriad  = 
        Graphics3D[{Thickness[0.01], RGBColor[1, 0, 0], 
            Line[{{0, 0, 0}, {35, 0, 0}}], 
            Text["X", {35, 0, -10}, 
              TextStyle ->  {FontFamily ->  "Times-Roman", 
                  FontSize ->  17}], RGBColor[0, 1, 0], 
            Line[{{0, 0, 0}, {0, 35, 0}}], 
            Text["Y", {0, 35, -10}, 
              TextStyle ->  {FontFamily ->  "Times-Roman", 
                  FontSize ->  17}], RGBColor[0, 0, 1], 
            Line[{{0, 0, 0}, {0, 0, 35}}], 
            Text["Z", {5, 5, 35}, 
              TextStyle ->  {FontFamily ->  "Times-Roman", 
                  FontSize ->  17}]}];

All the ITF, passionately

plotSnrAll = 
        SphericalPlot3D[
          1/2*L[10^(-23), φ - (π/2), π - ϑ, itfAll], {ϑ, 0, π}, {φ, 0, 
            2 π}, DisplayFunction -> Identity, 
          PlotRange -> All];
Show[plotSnrAll, plotFrameTriad, 
        DisplayFunction :> $DisplayFunction];
Show[plotSnrAll, plotFrameTriad, 
        DisplayFunction :>  $DisplayFunction, 
        ViewPoint -> {0, -1, 1}];
Show[plotSnrAll, plotFrameTriad, 
      DisplayFunction :> $DisplayFunction, 
      ViewPoint -> {1, 0, 1}]; 
    Show[plotSnrAll, plotFrameTriad, 
      DisplayFunction :> $DisplayFunction, 
      ViewPoint -> {3.057,  0.766,  1.231}];

Export["snrAll.eps", 
      Show[plotSnrAll, plotFrameTriad, 
        ViewPoint -> {3.057,  0.766,  1.231}, 
        DisplayFunction -> Identity], "epsi",  
      ImageResolution -> 160]

Just LIGO going alone

Show[plotSnrLIGO, 
        plotFrameTriad, 
        DisplayFunction :>  $DisplayFunction];
Show[plotSnrLIGO, plotFrameTriad, 
        DisplayFunction :> $DisplayFunction, 
        ViewPoint -> {0, -1, 1}];
Show[plotSnrLIGO, plotFrameTriad, 
        DisplayFunction :> $DisplayFunction, 
        ViewPoint -> {1, 0, 1}];
Show[plotSnrLIGO, plotFrameTriad, 
        DisplayFunction :> $DisplayFunction, 
        ViewPoint -> {3.057,  0.766,  1.231}];
Show[plotSnrLIGO, plotFrameTriad, 
        DisplayFunction :> $DisplayFunction, 
        ViewPoint -> {1, 1, 2}];

Export["snrLIGO.eps", 
      Show[plotSnrLIGO, plotFrameTriad, 
        ViewPoint -> {3.057,  0.766,  1.231}, 
        DisplayFunction -> Identity], "epsi",  
      ImageResolution -> 160]

Show[plotSnrAll, plotFrameTriad, 
        DisplayFunction :> $DisplayFunction, 
        ViewPoint -> {1, 1, 2}];

Show[plotSnrLIGO, plotFrameTriad, 
        DisplayFunction :> $DisplayFunction, 
        ViewPoint -> {0, 0, 2}];

Show[plotSnrGEOandVIRGO, plotFrameTriad, 
        DisplayFunction :> $DisplayFunction, 
        ViewPoint -> {0, 0, 2}];

Show[plotSnrAll, plotFrameTriad, 
        DisplayFunction :> $DisplayFunction, 
        ViewPoint -> {0, 0, 2}];

The 4km LIGO detectors

plotSnrLIGO4 = 
        SphericalPlot3D[
          1/2*L[10^(-23), φ - (π/2), π - ϑ, {H4K, L4K}], {ϑ, 0, π}, {φ, 0, 
            2  π},  PlotRange -> All, 
          DisplayFunction -> Identity];
Show[plotSnrLIGO4, plotFrameTriad, 
        DisplayFunction :>  $DisplayFunction];
Show[plotSnrLIGO4, plotFrameTriad, 
        DisplayFunction :> $DisplayFunction, 
        ViewPoint -> {0, -1, 1}];
Show[plotSnrLIGO4, plotFrameTriad, 
        DisplayFunction :> $DisplayFunction, 
        ViewPoint -> {1, 0, 1}];
Show[plotSnrLIGO4, plotFrameTriad, 
        DisplayFunction :> $DisplayFunction, 
        ViewPoint -> {3.057,  0.766,  1.231}];

Export["snrLIGO4.eps", 
      Show[plotSnrLIGO4, plotFrameTriad, DisplayFunction -> Identity, 
        ViewPoint -> {3.057,  0.766,  1.231}], "epsi",  
      ImageResolution -> 160]

LIGO and GEO

plotSnrLIGOandGEO = 
        SphericalPlot3D[
          1/2*L[10^(-23), φ - (π/2), π - ϑ, itfLIGOandGEO], {ϑ, 0, π}, {φ, 0, 
            2  π}, PlotRange -> All, 
          DisplayFunction -> Identity];
Show[plotSnrLIGOandGEO, plotFrameTriad, 
        DisplayFunction :> $DisplayFunction];
Show[plotSnrLIGOandGEO, plotFrameTriad, 
        DisplayFunction :> $DisplayFunction, 
        ViewPoint -> {0, -1, 1}];
Show[plotSnrLIGOandGEO, plotFrameTriad, 
        DisplayFunction :> $DisplayFunction, 
        ViewPoint -> {1, 0, 1}];
Show[plotSnrLIGOandGEO, plotFrameTriad, 
        DisplayFunction :> $DisplayFunction, 
        ViewPoint -> {3.057,  0.766,  1.231}];

Export["snrLIGOandGEO.eps", 
      Show[plotSnrLIGOandGEO, plotFrameTriad, 
        DisplayFunction -> Identity, 
        ViewPoint -> {3.057,  0.766,  1.231}], "epsi",  
      ImageResolution -> 160]

LIGO and VIRGO

plotSnrLIGOandVIRGO = 
        SphericalPlot3D[
          1/2*L[10^(-23), φ - (π/2), π - ϑ, itfLIGOandVIRGO], {ϑ, 0, π}, {φ, 0, 
            2  π}, PlotRange -> All, 
          DisplayFunction -> Identity];
Show[plotSnrLIGOandVIRGO, plotFrameTriad, 
        DisplayFunction :> $DisplayFunction];
Show[plotSnrLIGOandVIRGO, plotFrameTriad, 
        DisplayFunction :> $DisplayFunction, 
        ViewPoint -> {0, -1, 1}];
Show[plotSnrLIGOandVIRGO, plotFrameTriad, 
        DisplayFunction :> $DisplayFunction, 
        ViewPoint -> {1, 0, 1}];
Show[plotSnrLIGOandVIRGO, plotFrameTriad, 
        DisplayFunction :> $DisplayFunction, 
        ViewPoint -> {3.057,  0.766,  1.231}];

Export["snrLIGOandVIRGO.eps", 
      Show[plotSnrLIGOandVIRGO, plotFrameTriad, 
        DisplayFunction -> Identity, 
        ViewPoint -> {3.057,  0.766,  1.231}], "epsi",  
      ImageResolution -> 160]

GEO and VIRGO

plotSnrGEOandVIRGO = 
        SphericalPlot3D[
          1/2*L[10^(-23), φ - (π/2), π - ϑ, itfGEOandVIRGO], {ϑ, 0, π}, {φ, 0, 
            2  π}, PlotRange -> All, 
          DisplayFunction -> Identity];
Show[plotSnrGEOandVIRGO, plotFrameTriad, 
        DisplayFunction :> $DisplayFunction];
Show[plotSnrGEOandVIRGO, plotFrameTriad, 
        DisplayFunction :> $DisplayFunction, 
        ViewPoint -> {0, -1, 1}];
Show[plotSnrGEOandVIRGO, plotFrameTriad, 
        DisplayFunction :> $DisplayFunction, 
        ViewPoint -> {1, 0, 1}];
Show[plotSnrGEOandVIRGO, plotFrameTriad, 
        DisplayFunction :> $DisplayFunction, 
        ViewPoint -> {3.057,  0.766,  1.231}];

Export["snrGEOandVIRGO.eps", 
      Show[plotSnrGEOandVIRGO, plotFrameTriad, 
        DisplayFunction -> Identity, 
        ViewPoint -> {3.057,  0.766,  1.231}], "epsi",  
      ImageResolution -> 160]