Skip to content
GitLab
Menu
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
einsteinathome
pulsatingscience
Commits
da5a166a
Commit
da5a166a
authored
Sep 12, 2017
by
Oliver Bock
Browse files
Merge branch 'antenna-patterns' of gitlab:einsteinathome/pulsatingscience into antenna-patterns
parents
08406ff4
8c3f7e60
Changes
1
Hide whitespace changes
Inline
Side-by-side
src/antenna.c
0 → 100644
View file @
da5a166a
#include
<math.h>
#include
<stdio.h>
#include
<stdlib.h>
#include
<strings.h>
// For converting degrees to radians
const
double
deg_to_rad
=
M_PI
/
180
.
0
;
// Event time
// GPS 1187008882
// UTC: Thu 2017 Aug 17 12:41:04
// Julian Day 2457983.02852
// Modified Julian Day 57982.52852
// The convention in this code is that Lattitude is measured north of
// the equator. So Lattitude -10 means 10 degrees SOUTH of the
// equator. Longitude is measured East of Greenwich. So longitude
// -90 means 90 degrees West of Greenwich (somewhere in the USA).
// Galaxy NGC 4993 (from Wikipedia) Right ascension 13h 09m 47.2s
// Declination −23° 23′ 4″ DEC and RA in radians. We will need to
// convert this to location over the Earth at the event time. It will
// turn out to be Lattitude −23° 23′ 4″ and Longitude 41.092981
// degrees
// From Chapter 11 of "Astronomical Algorithms" by Jean Meeus, published 1991
void
print_galaxy_coordinates
()
{
// Here is the Julian day of the event, including a fractional part
double
JD
=
2457983
.
02
852
;
// This is the Julian day (note previous integer part!)
double
JD0
=
2457982
.
5
;
// Compute Julian century t
double
day
=
JD0
-
2451545
.
0
;
double
t
=
day
/
36525
.
0
;
// The following is equation 11.4 of the previous reference
double
degrees
=
280
.
46061837
+
360
.
98564736629
*
(
JD
-
2451545
.
0
)
+
t
*
t
*
0
.
0003
87933
-
t
*
t
*
t
/
38710000
.
0
;
degrees
=
fmod
(
degrees
,
360
.
0
);
printf
(
"Greenwich Mean Sidereal time is %f degrees
\n
"
,
degrees
);
// Now the longitude of our source (with sign conventions above) is given by L = RA - MST, where
// MST (in degrees) is given above and RA is the RA of our source
// Here a positive number means "East of Greenwich"
double
longitude
=
(
13
*
3600
.
0
+
9
*
60
+
47
.
2
)
/
240
-
degrees
;
printf
(
"Longitude of source is %f degrees
\n
"
,
longitude
);
}
// Note these are in the order lattitude,longitude. See routine
// "print_galaxy_coordinates" to convert location on sky to Earth
// location at event time.
const
double
galaxy
[
2
]
=
{
-
1
.
0
*
deg_to_rad
*
(
23
.
0
+
23
.
0
/
60
.
0
+
4
.
0
/
3600
.
0
),
deg_to_rad
*
41
.
092981
};
struct
Detector
{
// null terminated char string
char
name
[
8
];
// lattitude north of equator, radians
// longitude east from Greenwich, radians
double
location
[
2
];
// orientation of Y arm CCW from North, radians
double
orientation
;
// Unit vector from center of Earth to detector
double
vec
[
3
];
double
north
[
3
];
double
east
[
3
];
// Unit vectors along the two arms
double
lx
[
3
];
double
ly
[
3
];
};
// LLO, LHO, Virgo, in that order
struct
Detector
detectors
[
3
];
struct
Source
{
// lattitude north of equator, radians
// longitude east from Greenwich, radians
double
location
[
2
];
// Unit vector from Earth to source
double
vec
[
3
];
// Unit vectors defining plane perpendicular to the source
// direction
double
u
[
3
];
double
v
[
3
];
// Angle defining line of source ellipse
double
orientation
;
};
struct
Source
source
;
// input is a lattitude/longitude set in radians
// output is unit vectors
void
make_unit_vectors
(
double
*
out
,
double
*
in
)
{
double
lat
=
in
[
0
];
double
lon
=
in
[
1
];
out
[
0
]
=
cos
(
lon
)
*
cos
(
lat
);
out
[
1
]
=
sin
(
lon
)
*
cos
(
lat
);
out
[
2
]
=
sin
(
lat
);
// sanity check
double
len
=
out
[
0
]
*
out
[
0
]
+
out
[
1
]
*
out
[
1
]
+
out
[
2
]
*
out
[
2
];
if
(
fabs
(
len
-
1
.
0
)
>
0
.
0001
)
{
fprintf
(
stderr
,
"Squares check failed!
\n
"
);
exit
(
1
);
}
}
void
make_u_v_vectors
(
struct
Source
*
src
)
{
// lattitude and longitude
double
lat
=
src
->
location
[
0
];
double
lon
=
src
->
location
[
1
];
// construct unit vectors perpendicular to line of sight
double
u
[
3
],
v
[
3
];
u
[
0
]
=
-
sin
(
lat
)
*
cos
(
lon
);
u
[
1
]
=
-
sin
(
lat
)
*
sin
(
lon
);
u
[
2
]
=
cos
(
lat
);
v
[
0
]
=
-
sin
(
lon
);
v
[
1
]
=
cos
(
lon
);
v
[
2
]
=
0
.
0
;
// now make a rotated basis, rotated by oriention angle
double
cpsi
=
cos
(
src
->
orientation
);
double
spsi
=
sin
(
src
->
orientation
);
int
i
;
for
(
i
=
0
;
i
<
3
;
i
++
)
{
src
->
u
[
i
]
=
cpsi
*
u
[
i
]
+
spsi
*
v
[
i
];
src
->
v
[
i
]
=
-
spsi
*
u
[
i
]
+
cpsi
*
v
[
i
];
}
}
// input is a lattitude/longitude set in radians
// output is unit vectors in the north and east directions
void
make_north_east_vectors
(
struct
Detector
*
det
)
{
double
lat
=
det
->
location
[
0
];
double
lon
=
det
->
location
[
1
];
det
->
north
[
0
]
=
-
sin
(
lat
)
*
cos
(
lon
);
det
->
north
[
1
]
=
-
sin
(
lat
)
*
sin
(
lon
);
det
->
north
[
2
]
=
cos
(
lat
);
det
->
east
[
0
]
=
-
sin
(
lon
);
det
->
east
[
1
]
=
cos
(
lon
);
det
->
east
[
2
]
=
0
.
0
;
double
cpsi
=
cos
(
det
->
orientation
);
double
spsi
=
sin
(
det
->
orientation
);
int
i
;
for
(
i
=
0
;
i
<
3
;
i
++
)
{
det
->
lx
[
i
]
=
cpsi
*
det
->
east
[
i
]
+
spsi
*
det
->
north
[
i
];
det
->
ly
[
i
]
=
-
spsi
*
det
->
east
[
i
]
+
cpsi
*
det
->
north
[
i
];
}
}
// sets out = a X b
void
cross_prod
(
double
*
out
,
double
*
a
,
double
*
b
)
{
out
[
0
]
=
a
[
1
]
*
b
[
2
]
-
a
[
2
]
*
b
[
1
];
out
[
1
]
=
a
[
2
]
*
b
[
0
]
-
a
[
0
]
*
b
[
2
];
out
[
2
]
=
a
[
0
]
*
b
[
1
]
-
a
[
1
]
*
b
[
0
];
}
void
populate_source
()
{
source
.
location
[
0
]
=
-
1
.
0
*
deg_to_rad
*
(
23
.
0
+
23
.
0
/
60
.
0
+
4
.
0
/
3600
.
0
);
source
.
location
[
1
]
=
deg_to_rad
*
41
.
092981
;
// Initially set to zero, fix later
source
.
orientation
=
0
.
0
;
make_unit_vectors
(
source
.
vec
,
source
.
location
);
make_u_v_vectors
(
&
source
);
//
}
// initially take detector locations and orientations from
// https://arxiv.org/pdf/gr-qc/9607075.pdf
// Be sure to check this later!!
void
populate_detectors
(){
// LLO
// Sanity checked using Google Earth!
strcpy
(
detectors
[
0
].
name
,
"LLO"
);
detectors
[
0
].
location
[
0
]
=
30
.
56
*
deg_to_rad
;
detectors
[
0
].
location
[
1
]
=
-
90
.
77
*
deg_to_rad
;
detectors
[
0
].
orientation
=
198
.
0
*
deg_to_rad
;
// LHO
// Sanity checked using Google Earth!
strcpy
(
detectors
[
1
].
name
,
"LHO"
);
detectors
[
1
].
location
[
0
]
=
46
.
45
*
deg_to_rad
;
detectors
[
1
].
location
[
1
]
=
-
119
.
41
*
deg_to_rad
;
detectors
[
1
].
orientation
=
126
.
8
*
deg_to_rad
;
// VIRGO
// Sanity checked using Google Earth!
strcpy
(
detectors
[
2
].
name
,
"Virgo"
);
detectors
[
2
].
location
[
0
]
=
43
.
63
*
deg_to_rad
;
detectors
[
2
].
location
[
1
]
=
10
.
5
*
deg_to_rad
;
detectors
[
2
].
orientation
=
71
.
5
*
deg_to_rad
;
int
i
;
// Coordinate system has x/y plane through the equator, north pole
// along positive z axis, and Greenwich passing through x axis (ie
// y=0).
for
(
i
=
0
;
i
<
3
;
i
++
)
{
make_unit_vectors
(
detectors
[
i
].
vec
,
detectors
[
i
].
location
);
make_north_east_vectors
(
detectors
+
i
);
}
}
void
print_detector
(
int
det
)
{
printf
(
"name: %s
\n
"
"lattitude (north): %f
\n
"
"longitude (east): %f
\n
"
"orientation: (CCW from North): %f
\n
"
"earth center to detector %f %f %f
\n\n
"
,
detectors
[
det
].
name
,
detectors
[
det
].
location
[
0
],
detectors
[
det
].
location
[
1
],
detectors
[
det
].
orientation
,
detectors
[
det
].
vec
[
0
],
detectors
[
det
].
vec
[
1
],
detectors
[
det
].
vec
[
2
]
);
}
void
print_source
()
{
printf
(
"source at:
\n
"
"lattitude: %f
\n
"
"longitude: %f
\n
"
"earth center to source %f %f %f
\n\n
"
,
source
.
location
[
0
],
source
.
location
[
1
],
source
.
vec
[
0
],
source
.
vec
[
1
],
source
.
vec
[
2
]
);
}
// This dots the mass quadrupole with the antenna functions
void
get_UUUVVV
(
int
det
,
double
*
uu
,
double
*
uv
,
double
*
vv
)
{
double
UU
=
0
,
VV
=
0
,
UV
=
0
;
int
i
,
j
;
for
(
i
=
0
;
i
<
3
;
i
++
)
for
(
j
=
0
;
j
<
3
;
j
++
)
UU
+=
source
.
u
[
i
]
*
source
.
u
[
j
]
*
(
detectors
[
det
].
lx
[
i
]
*
detectors
[
det
].
lx
[
j
]
-
detectors
[
det
].
ly
[
i
]
*
detectors
[
det
].
ly
[
j
]);
for
(
i
=
0
;
i
<
3
;
i
++
)
for
(
j
=
0
;
j
<
3
;
j
++
)
UV
+=
(
source
.
u
[
i
]
*
source
.
v
[
j
]
+
source
.
v
[
i
]
*
source
.
u
[
j
])
*
(
detectors
[
det
].
lx
[
i
]
*
detectors
[
det
].
lx
[
j
]
-
detectors
[
det
].
ly
[
i
]
*
detectors
[
det
].
ly
[
j
]);
for
(
i
=
0
;
i
<
3
;
i
++
)
for
(
j
=
0
;
j
<
3
;
j
++
)
VV
+=
source
.
v
[
i
]
*
source
.
v
[
j
]
*
(
detectors
[
det
].
lx
[
i
]
*
detectors
[
det
].
lx
[
j
]
-
detectors
[
det
].
ly
[
i
]
*
detectors
[
det
].
ly
[
j
]);
*
uu
=
UU
;
*
uv
=
UV
;
*
vv
=
VV
;
}
int
main
(
int
argc
,
char
*
argv
[])
{
double
iota
=
atof
(
argv
[
1
]);
printf
(
"Iota = %f degrees
\n
"
,
iota
);
iota
*=
deg_to_rad
;
double
ci
=
cos
(
iota
);
double
si
=
sin
(
iota
);
// setup and test
// print_galaxy_coordinates();
populate_detectors
();
populate_source
();
// for (i=0; i<3; i++) print_detector(i);
// print_source();
// Now find waveforms. At each site, waveform is A sin^2(wt) + B cos^2(wt) + C sin(wt)cos(wt)
double
UU
,
UV
,
VV
;
int
i
;
for
(
i
=
0
;
i
<
3
;
i
++
)
{
double
A
=
0
,
B
=
0
,
C
=
0
;
get_UUUVVV
(
i
,
&
UU
,
&
UV
,
&
VV
);
A
=
ci
*
ci
*
(
2
.
0
*
VV
/
3
.
0
-
UU
/
3
.
0
);
B
=
2
.
0
*
(
UU
+
VV
)
/
3
.
0
;
C
=
ci
*
UV
;
printf
(
"For detector %s the waveform is %f sin^2(wt) + %f cos^2(wt) + %f sin(wt)cos(wt)
\n
"
,
detectors
[
i
].
name
,
A
,
B
,
C
);
}
return
0
;
}
Write
Preview
Supports
Markdown
0%
Try again
or
attach a new file
.
Attach a file
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment