-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathmainmatter.tex
703 lines (593 loc) · 27.3 KB
/
mainmatter.tex
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
\section{Introduction}
Since we have to work with many coordinate systems it can be hard to
keep track of the definitions of each. This document is meant as a
reference to easily find the relations between the different systems.
First coordinate systems used by \hisparc will be discussed. Including
the units that are used and where it is used. Then other coordinate
systems that we have to deal with are discussed, including ways to
convert from those to our usual coordinate systems.
In \cref{sec:geographic} we discuss geographic coordinates to define the
coordinates of the observer (detection station and scintillator detectors). In
\cref{sec:time} we discuss different time keeping systems and how we can
relate a time of observation to the rotation of the Earth w.r.t the celestial
sphere. The systems to define a point on the celestial sphere are described in
\cref{sec:celestial}. In \cref{sec:corsika} the coordinate systems which are
used in \corsika simulations are described. Finally in \cref{sec:example} an
example shows the process to go from a \hisparc detection to equatorial
coordinates.
Make sure that angles are converted to the correct units for the
cosine and sine functions when converting between coordinate systems.
The use of radians or degrees as units depends on the programming language.
\section{Geographic}
\label{sec:geographic}
Geographic coordinate systems define a point on the Earth. For \hisparc
two systems are used: WGS84 and ENU. \cref{fig:wgs84_ecef_enu}
illustrates the relationships between these systems. Additionally a
compass based system is used to determine detector locations relative to
the \gps antenna. These systems are described in the following sections.
The conversion formulae are taken from \cite[sec.~K]{usno:2014aa}.
\begin{figure}
\centering
\input{figures/WGS84_ECEF_ENU}
\caption{Relationship between the WGS84 (orange), ECEF (blue) and ENU
(green) coordinate systems. The Prime meridian is the IERS
Reference Meridian.}
\label{fig:wgs84_ecef_enu}
\end{figure}
\subsection{World Geodetic System 1984 (WGS84)}
The location of a \hisparc station is determined by means of a \gps
antenna (one antenna for each station). The \gps returns coordinates in
the WGS84 coordinate system. This defines a position with a latitude
($\varphi$), longitude ($\lambda$) and altitude ($h$). The latitude and
longitude are defined in decimal degrees, the altitude in meters.
Latitude is the angular distance between a location and the equator,
angles towards north are positive. The longitude is the angular distance
of a location to the Reference Meridian (from now on referred to as the
prime meridian) defined by the International Earth Rotation and
Reference Systems Service (IERS), angles towards east are positive. The
altitude is the height above an ellipsoid that approximates the shape of
the Earth. The parameters which define this ellipsoid are given in
\cref{eq:wgs84} \cite{nga:2014aa}.
\subsection{Earth-Centered, Earth-Fixed (ECEF)}
The ECEF coordinate system is used as intermediate between WGS84 and
ENU. ECEF has its origin as the center of the Earth (Earth-Centered),
and its axes are aligned with the WGS84 system. The z-axis points to
North. The x-axis crosses the Earth surface where the equator
(\SI{0}{\degree} latitude) intersects the prime meridan (\SI{0}{\degree}
longitude). The y-axis intersects the equator at \SI{90}{\degree}
longitude. The axes rotate with the Earth (Earth-Fixed). A position is
defined by the $X$, $Y$, and $Z$ coordinates, each given in meters.
The formulae to convert WGS84 coordinates to ECEF coordinates are:
%
\begin{equation}
\begin{array}{l}
X = \left(N + h\right) \cos{\varphi} \cos{\lambda} \ , \\
Y = \left(N + h\right) \cos{\varphi} \sin{\lambda} \ , \\
Z = \left(\frac{b^2}{a^2} N + h\right) \sin{\varphi} \ .
\end{array}
\end{equation}
%
For the paramaters that define the ellipsoidal approximation of the
Earth surface we have
%
\begin{equation}
\label{eq:wgs84}
\begin{array}{l}
a = \SI{6378137.0}{\meter} \ , \\
f = \frac{1}{298.257223563} \ , \\
b = a (1 - f) \ , \\
e = \sqrt{2 f - f^2} \ , \\
N = \frac{a}{\sqrt{1 - e^2 \sin^2 \varphi}} \ ,
\end{array}
\end{equation}
%
here $a$ is the semi-major axis, $f$ the flattening, $b$ the
semi-minor axis and $e$ the eccentricity of the ellipsoid. $N$ is the
Normal, which is the distance between a location on the ellipsoid and
the intersection of its normal and the z-axis of the ellipsoid. For a
position on the equator (i.e. $\varphi = \SI{0}{\degree}$) $N = a$.
\begin{figure}
\centering
\input{figures/Ellipsoid}
\caption{Ellipsoidal approximation of Earth. Showing the definitions
of the semi-major axis, the semi-minor axis, and the Normal.}
\label{fig:ellipsoid}
\end{figure}
\subsection{East, North, Up (ENU)}
East, North, Up is used to obtain relative locations and angles with
respect to a reference position. It is a local coordinate system in a
plane tangential to the WGS84 elipsoid. From the reference position the
positive x-axis (East) is in the eastern direction, the positive y-axis
(North) is towards north, and the positive z-axis (Up) is towards the
Zenith. All distances are in meters.
The formulae to convert ECEF to ENU with a given reference position are:
%
\begin{equation}
\begin{bmatrix}
x \\
y \\
z
\end{bmatrix}
=
\begin{bmatrix}
-\sin(\lambda_r) & \cos(\lambda_r) & 0. \\
-\sin(\varphi_r) \cos(\lambda_r) & -\sin(\varphi_r) \sin(\lambda_r) & \cos(\varphi_r) \\
\cos(\varphi_r) \cos(\lambda_r) & \cos(\varphi_r) \sin(\lambda_r) & \sin(\varphi_r)
\end{bmatrix}
\begin{bmatrix}
X - X_r \\
Y - Y_r \\
Z - Z_r
\end{bmatrix}
\end{equation}
%
where $(X, Y, Z)$ and $(x, y, z)$ are the location in ECEF and
ENU respectively. $(\varphi_r, \lambda_r)$ and $(X_r, Y_r, Z_r)$ are the
reference position in WGS84 and ECEF respectively.
\subsection{`Compass'}
For reconstruction of events the precise location of the individual
detectors in each station is required. The location of detectors are
taken relative to the \gps antenna, because that position is known. We
use a system that requires a distance and angle measurement, see
\cref{fig:enu_compass}. The distance ($r$) between the \gps and a
detector (center of the scintillator) is measured (in meters), the angle
($\alpha$) relative to the geographic north pole (not magnetic north!)
is determined using a compass and the magnetic declination. This angle
is in degrees and goes from North to East (NESW). The height difference
($z$) between the \gps and detectors can also be measured, however, this
is only important if the detectors are not all at the same altitude or
are at a very different altitude from the \gps. The value of $z$ is positive
if the detector is higher than the \gps antenna. The rotation angle of
the detector itself ($\beta$) is defined as the angle of the long side
of the detector relative to the North Pole (positive towards East).
When a compass is used the angle needs to be compensated for the
difference between magnetic and geographic north, the difference is
called the magnetic declination. The magnetic declination is date and
location dependent \cite{canada:2013aa}. The magnetic declination can be
up to \SI{3}{\degree} for some \hisparc station locations.
\begin{figure}
\centering
\input{figures/ENU_compass}
\caption{Relationship between the ENU (green) and Compass (magenta)
coordinate system. The detectors of station 503 are shown.
Its \gps is used as the origin for the ENU coordinates. The
$\alpha$ and $r$ coordinate of detector 4 (blue), and the
$\beta$ coordinate of detector 1 (black) are illustrated.
The values are given in \cref{sec:example}.}
\label{fig:enu_compass}
\end{figure}
The formulae to convert the location of a station in Compass coordinates
to ENU are:
%
\begin{equation}
\begin{array}{l}
x = r * \sin{\alpha} \ , \\
y = r * \cos{\alpha} \ , \\
z = z \ .
\end{array}
\end{equation}
\section{Time}
\label{sec:time}
For \hisparc events it is important to know precisely when they occur, a
\gps is used to get synchronized and precise timing. It gives
unique timestamps. GPS time needs to be converted to Local Sidereal Time for
the conversion of celestial coordinates in \cref{sec:celestial}. The
full transformation chain is as follows: GPS $\to$ UTC $\to$ JD $\to$
GMST $\to$ LST.
\subsection{\gps time}
\gps time started counting on 6 January 1980 00:00 UTC (Coordinated
Universal Time), i.e. \num{315964800} UNIX seconds after the UNIX epoch
of 1 January 1970 00:00 UTC. \gps time is a continuously increasing
number. This is different from UTC time in which certain seconds are
repeated when `leap seconds' are added \cite{usno:2012aa}. Since 1
July 2015 the difference between UTC and \gps time is 17 seconds.
From the \hisparc electronics we get a \gps timestamp, adjusted to the
UNIX epoch, i.e. $\mathrm{GPS} + \mathrm{UNIX(GPS=0)}$. The timestamp is
given in nanoseconds. The \gps timestamp is therefore a large number
which requires 64-bits to be represented digitally. For instance the
\gps timestamp for 1 December 2014 at 00:00:29.886222166 (\gps time) is
%
\begin{equation}
\SI{1101427229886222166}{\nano\second} +
\SI{315964800000000000}{\nano\second} =
\SI{1417392029886222166}{\nano\second} \ .
\end{equation}
\subsection{UTC time}
UTC time is similar to \gps but corrected with leap seconds to keep it more in
sync with solar time. A leap second means that a second is counted twice, one
timestamp refers to two real seconds. The IERS decides when to add leap
seconds. \hisparc initially used UTC time but switched to \gps time around
2010-2011. UTC time can cause issues with duplicate timestamps. This makes it
impossible to distinguish between two real seconds.
\subsubsection{Leap seconds}
Leap seconds can be introduced at the end of either 30 June or 31
December, the leap second is represented in UTC time as 23:59:60.
\cref{table:leapseconds} gives a list of the recent leap seconds,
which are important for \hisparc. Leap seconds are announced about 6
months before they take effect.
\begin{table}
\centering
\begin{tabular}{ l c }
\toprule
Date & Leap seconds \\
\midrule
31 December, 1998 & 13 \\
31 December, 2005 & 14 \\
31 December, 2008 & 15 \\
30 June, 2012 & 16 \\
30 June, 2015 & 17 \\
\bottomrule
\end{tabular}
\caption{Leap seconds in effect or introduced after 2002. This is the
amount of seconds that \gps time is ahead of UTC time.}
\label{table:leapseconds}
\end{table}
\subsection{Julian Date (JD)}
JD is a decimal number counting the number of days since 1 January 4713
B.C. 12:00:00. Calculating this is not trivial since one has to account
for the fact that the 5th up to and including the 14th of October in
1582 do not exist. Those ten days were skipped because the Julian
calendar, which was used before that gap, introduced too many leap
years, every year divisible by 4 was a leap year. The Gregorian calendar
used afterwards compensated by having less leap years. Leap years still
occur when years are divisible by 4, except when they are also divisible
by 100 (e.g. 1900), except when they are also divisible by 400 (e.g.
2000) \cite{acf:2014aa}. Consider these examples: 1999 is not divisible
by 4 so not a leap year, 2004 is divisible by 4 and not by 100 and is a
leap year, 1900 is divisible by 4 and 100 but not 400 so is not a leap
year, and 2000 is divisible by 4, 100 and 400 and is therefore a leap
year.
The following conversion is valid for dates after 15 October 1582. The
brackets $\floor{\ }$ indicate that the decimal part of the value should
be discarded. The $\mathit{month}$ starts at 1 for January and
$\mathit{hour}$ is the decimal number of hours.
%
\begin{equation}
\begin{array}{l}
a = \floor{\frac{14 - \mathit{month}}{12}} \ , \\
A = \floor{\frac{\mathit{year - a}}{100}} \ , \\
B = 2 - A + \floor{\frac{A}{4}} \ , \\
C = \floor{365.25 (\mathit{year} - a)} \ , \\
D = \floor{30.6001 (\mathit{month} + 12 a + 1)} \ , \\
\mathit{JD} = B + C + D + \mathit{day} +
\frac{\mathit{hour}}{24} + 1720994.5 \ .
\end{array}
\end{equation}
\subsubsection{Modified Julian Date (MJD)}
Modified Julian Date to keep number shorter uses epoch 17 November 1858
00:00:00. The difference between a Modified Julian Date and a Julian
Date is 2400000.5 days. The .5 days moves the reference time from noon
to midnight.
%
\begin{equation}
\mathit{MJD} = \mathit{JD} - 2400000.5 \ .
\end{equation}
\subsection{Sidereal Time}
In sidereal time the rotation of Earth relative to the `fixed'
background stars is taken. A sidereal day is a full rotation of Earth
around its axis, relative to background stars. In a solar day the Earth
rotates more than \SI{360}{\degree} around its axis because it orbits
around the Sun, causing the position of the Sun relative to the Earth to
change. The Earth orbits the Sun every 365.25 days, so in one day the
orbit proceeds by $\frac{\SI{360}{\degree}}{\SI{365.25}{\day}} \approx
\SI{1}{\degree\per\day}$. So an extra rotation of approximately
\SI{1}{\degree} around its axis is required to make the same part of the
Earth face the Sun (a solar day). A sidereal day is therefore around
\SI{4}{\minute} shorter than a solar day. In \cref{fig:sidereal_time}
the difference between a solar and sidereal day is illustrated.
\begin{figure}
\centering
\input{figures/Sidereal_time}
\caption{This shows the orbit of Earth (green) around the Sun
(yellow). The black line on the Earth represents the prime
meridian. A full rotation of the Earth causes the prime
meridian to point to the same background stars, this is a
sidereal day. It takes longer, due to the orbit of the
Earth, for the prime meridian to point to the Sun again,
which is a solar day. Used sizes and angles are
illustrative, not realistic.}
\label{fig:sidereal_time}
\end{figure}
\subsubsection{Greenwich Mean Sidereal Time (GMST)}
GMST is the hour angle of the vernal equinox (\cref{sec:celestial})
with respect to the prime meridian at Greenwich, expressed in decimal
hours \cite{kaplan:2011aa}.
%
\begin{equation}
\mathit{GMST} = 6.697374558 + 0.06570982441908\ \mathit{D0} +
1.00273790935\ H + 0.000026\ T^2 \ .
\end{equation}
%
Here $D0$ is the Julian Date of the previous midnight using the J2000
epoch (1 January, 2000 12:00 UTC or Julian date 2451545.0), $H$ is the
number of hours since the last midnight, and $T$ is the number of
centuries since J2000.
\begin{figure}
\centering
\input{figures/WGS84_GMST_LST}
\caption{Top-down view of Earth showing the relation between GMST
(teal) and LST (blue) for an observer at a longitude
$\lambda$ (orange). The observer and prime meridian move with
the rotation of the Earth, the vernal equinox does not.}
\label{fig:wgs84_gmst_lst}
\end{figure}
\subsubsection{Local Sidereal Time (LST)}
Similar to GMST but takes observers longitude into account. To get the
Local Sidereal Time the observers longitude has to be added to the GMST.
To convert the longitude from degrees to hours it has to be divided by
15, $\frac{\SI{360}{\degree}}{\SI{24}{\hour}} =
\SI{15}{\degree\per\hour}$.
%
\begin{equation}
\begin{array}{l}
\mathit{LST} = \mathit{GMST} + \frac{\lambda}{15} \ .
\end{array}
\end{equation}
%
The relation between LST and GMST is illustrated in
\cref{fig:wgs84_gmst_lst}.
\section{Celestial}
\label{sec:celestial}
Here we describe the different coordinate systems that are used to
define the direction (origin) of an air shower. The zenith and azimuth
coordinates rotate with the rotation of the Earth, because they are
anchored to the location of the observer. The Equatorial system is
linked to the celestial sphere and is independent from the rotation of
Earth, the time of observation, and the location of the observer. To
transform from zenith and azimuth to Equatorial the position of the
observer and time of observation are required. The Horizontal coordinate
system is used as intermediate, it also defines an azimuth coordinate,
this is different from the azimuth coordinate in the zenith and azimuth
system.
\subsection{Zenith and azimuth coordinates}
When a station detects a shower we try to reconstruct the direction of
its origin, i.e. the position on the sky the shower came from. This
direction is then given by a zenith ($\theta$) and azimuth ($\phi$)
coordinate. These coordinates define a point on the semi-sphere that is
the sky above the detection station. The Zenith is the point directly
above the observer. The zenith angle is the angle between the direction
and the Zenith point, so straight up is \SI{0}{\radian} and the horizon
\SI{\pi / 2}{\radian}. The azimuth is the direction in the horizontal
plane, it starts at East (\SI{0}{\radian}) then goes to North (ENWS).
This is illustrated in pink in \cref{fig:enu_zenazi}.
We neither expect nor consider air showers from below the horizon, so the
zenith angles, defined in radians, are an angle in the range [0,
$\frac{\pi}{2}$). The azimuth is restricted to the range [$-\pi$, $\pi$).
\begin{figure}
\centering
\input{figures/ENU_ZenAzi}
\caption{Relationship between the ENU (green) and Zenith Azimuth (pink)
coordinate systems. The zenith and azimuth coordinates define
a point on the sky.}
\label{fig:enu_zenazi}
\end{figure}
\subsection{Horizontal coordinate system}
This is a system used as intermediary for some coordinate conversions.
It uses azimuth ($A$) and altitude ($a$) to define a direction on the sky. The
altitude is the complement of the zenith, so \SI{0}{\radian} is horizontal
and \SI{\pi / 2}{\radian} is the Zenith. The azimuth definitions between the
zenith and azimuth, and Horizontal systems differ, in Horizontal
coordinates it moves from North to East (NESW). This is illustrated in
gold in \cref{fig:enu_zenazi_horizontal}.
The formulae to convert from zenith ($\theta$) and azimuth ($\phi$) to
altitude ($a$) and azimuth ($A$) are:
%
\begin{equation}
\begin{array}{l}
a = \frac{\pi}{2} - \theta \ , \\
A = \frac{\pi}{2} - \phi \ .
\end{array}
\end{equation}
Conversions can cause values to go beyond the allowed range of values
for angles or times. They may have to be brought back into the range
after the conversion.
\begin{figure}
\centering
\input{figures/ENU_ZenAzi_Horizontal}
\caption{Relationship between the ENU (green), Zenith Azimuth (pink)
and horizontal (gold) coordinate systems.}
\label{fig:enu_zenazi_horizontal}
\end{figure}
\subsection{Equatorial (J2000)}
\cref{fig:equatorial} shows the relation between the celestial sphere
and Equatorial coordinates. The right ascension ($\alpha$) is the angle
between the projection of the position on the celestial sphere ($\star$)
on the plane of the Celestial equator and the vernal equinox
($\vernal$). The declination ($\delta$) is the angle between the plane
of the Celestial equator and the sky position. Due to precession (change
of direction of Earth's rotational axis) the celestial positions slowly
change over time. To compensate for precession coordinates are
calculated as they would have been on a specific date, the J2000 epoch
is commonly used. The equatorial coordinates are fixed to the celestial
sphere and not influenced by the rotation of the Earth.
\begin{figure}
\centering
\input{figures/Equatorial}
\caption{This figures shows the relation between a position on the
sky (black), Equatorial coordinates (red), and the
Celestial sphere.}
\label{fig:equatorial}
\end{figure}
Given an observers position in WGS84, a time of observation in LST
(derived from the \gps timestamp and observer longitude) and horizontal
coordinates (converted from the zenith and azimuth coordinates) pointing
to the source the following formulae can be used to get the
corresponding equatorial coordinates \cite[p. 37]{duffet-smith:1990aa}.
%
\begin{equation}
\label{eq:equatorial}
\begin{array}{l}
\delta = \arcsin{\left((\sin{a} \sin{\varphi}) +
(\cos{a} \cos{\varphi} \cos{A})\right)} \ , \\
\mathit{HA} = \arccos{\left(\frac{\sin{a} - (\sin{\varphi} \sin{\delta})}
{\cos{\varphi} \cos{\delta}}\right)} \ , \\
\alpha = \mathit{LST} - \mathit{HA} \ .
\end{array}
\end{equation}
The output range of $\arccos$ is [0, $\pi$), while that of $\mathit{HA}$ is
[0, $2\pi$). If the azimuth is positive use: $\mathit{HA} = 2 \pi -
\mathit{HA}$.
Where $\varphi$ the geographical latitude, $a$ the altitude, $A$ the
azimuth, $\mathit{LST}$ the Local sidereal time, $\delta$ is the
declination, $\alpha$ the right ascension and the intermediate variable
$\mathit{HA}$ is the hour angle. This is illustrated in
\cref{fig:wgs84_zenazi_lst_equatorial}.
\begin{figure}
\centering
\input{figures/WGS84_ZenAzi_LST_Equatorial}
\caption{Relationship between the WGS84 (orange), ENU (green),
Zenith azimuth (pink), and clock (teal) and Equatorial
(red) coordinate systems. $\vernal$ points to the vernal
equinox, the point where the ecliptic and celestial equator
cross in the Aries zodiac. WGS84 defines the position of
the observer, LST defines the time of observation, ENU is
the local coordinate system (only East is shown because it
is used for the azimuth), and zenith and azimuth give the
direction of shower origin. Combined these give the same
position on the sky as the equatorial coordinates.}
\label{fig:wgs84_zenazi_lst_equatorial}
\end{figure}
\section{\corsika}
\label{sec:corsika}
\corsika is an extensive air shower simulation program which we use. It
uses coordinate systems which are slightly different from the ones we
employ \cite{heck:2013aa}. The \sapphire framework provides a script to
convert \corsika output to \hdf format. During this conversion the
coordinates are transformed to the \hisparc coordinate system. The
following sections describe these transformations.
\begin{figure}
\centering
\input{figures/ENU_ZenAzi_CORSIKA}
\caption{Relationship between the ENU (green), Zenith Azimuth (pink),
and \corsika (brown) coordinates.}
\label{fig:enu_zenazi_corsika}
\end{figure}
\subsection{Geographic}
\corsika defines positions on the ground (or observation level) relative
to the point where the shower axis intersects the observation level.
Positive x axis points to magnetic North, positive y axis to the West,
and the z axis upwards. This is shown by the dashed brown lines in
\cref{fig:enu_zenazi_corsika}.
The formulae to convert \corsika to ENU are:
%
\begin{equation}
\begin{array}{l}
x_{\mathrm{ENU}} = -y_{\mathrm{CORSIKA}} \ , \\
y_{\mathrm{ENU}} = x_{\mathrm{CORSIKA}} \ , \\
z_{\mathrm{ENU}} = z_{\mathrm{CORSIKA}} \ . \\
\end{array}
\end{equation}
\subsection{Celestial}
\corsika looks from the point of view of the shower, so not the
direction it came form, but the direction it goes to. The \corsika
zenith angle ($\Theta$) is identical to the \hisparc zenith ($\theta$),
\SI{0}{\radian} is a shower from the zenith and \SI{\pi / 2}{\radian} is
a horizontal shower. The \corsika azimuth angle ($\Phi$) differs from
the \hisparc definition of azimuth ($\phi$). The \corsika azimuth is the
angle the shower is heading to with respect to the North, while the
\hisparc azimuth is the angle the shower is coming from with respect to
the East. $\Phi = \SI{0}{\radian}$ is a shower heading towards North, so
coming from South, which we would define as $\phi = \SI{-\pi /
2}{\radian}$. The (positive) rotation of the angle is in the same
direction, from North to West. This is shown by the brown arced lines in
\cref{fig:enu_zenazi_corsika}.
The formulae to convert the \corsika direction to \hisparc zenith and
azimuth are:
%
\begin{equation}
\begin{array}{l}
\theta = \Theta \ , \\
\phi = \Phi - \frac{\pi}{2} \ . \\
\end{array}
\end{equation}
\section{Example}
\label{sec:example}
An example will now be given for a real \hisparc event detected by
station 503.
\subsection{Measured}
From the \gps of station 503 we know its geographical location in WGS84,
%
\begin{equation}
(\varphi, \lambda, h) = (\SI{52.3562600}{\degree},
\SI{4.9529440}{\degree},
\SI{51.4}{\meter}) \ .
\end{equation}
The location of its four detectors have been measured in compass
coordinates (see \cref{fig:enu_compass})
\begin{center}
\begin{tabular}{ l c c c c }
\toprule
Detector & $\alpha$ [\si{\degree}] & $r$ [\si{\meter}] &
$\beta$ [\si{\degree}] & $z$ [\si{\meter}] \\
\midrule
1 & 315 & 8.97 & 135 & 0 \\
2 & 315 & 3.15 & 135 & 0 \\
3 & 225 & 5.09 & 225 & 0 \\
4 & 45 & 4.89 & 225 & 0 \\
\bottomrule
\end{tabular}
\end{center}
An event was detected on \gps timestamp
\SI{1333018296870008589}{\nano\second}. The relative arrival times in
each detector are \SIlist{0;2.5;5;12.5}{\nano\second}.
\subsection{Conversions}
To reconstruct the direction of the shower we use \sapphire. First it
converts the detector compass coordinates to ENU (relative to the \gps)
\begin{center}
\begin{tabular}{ l c c c }
\toprule
Detector & $x$ [\si{\meter}] & $y$ [\si{\meter}] & $z$ [\si{\meter}] \\
\midrule
1 & -6.34 & 6.34 & 0 \\
2 & -2.23 & 2.23 & 0 \\
3 & -3.60 & -3.60 & 0 \\
4 & 3.46 & 3.46 & 0 \\
\bottomrule
\end{tabular}
\end{center}
The local zenith and azimuth direction is then reconstructed using the
relative arrival times and the ENU positions. In this case
%
\begin{equation}
\begin{array}{l}
\theta = \SI{0.3818}{\radian} \ , \\
\phi = \SI{3.0030}{\radian} \ ,
\end{array}
\end{equation}
%
which is equal to these horizontal coordinates
%
\begin{equation}
\begin{array}{l}
a = \SI{1.1890}{\radian} \ , \\
A = \SI{-1.4322}{\radian} \ .
\end{array}
\end{equation}
The \gps timestamp represents the date 29 March 2012 at 10:51:21. On
that date 15 leap seconds were in effect. The corresponding UTC
timestamp, Julian Date, and GMST are
%
\begin{equation}
\begin{array}{l}
\mathit{UTC} = \SI{1333018281870008589}{\nano\second} \ , \\
\mathit{JD} = \SI{2456015.952336}{\day} \ , \\
\mathit{GMST} = \SI{23.3389}{\hour} = 23^\mathrm{h}20^\mathrm{m}20^\mathrm{s} \ . \\
\end{array}
\end{equation}
Using the longitude of the station the Local Sidereal Time can be
calculated from the GMST
%
\begin{equation}
\mathit{LST} = \SI{23.6691}{\hour} = 23^\mathrm{h}40^\mathrm{m}9^\mathrm{s} \ .
\end{equation}
Now the equatorial coordinates can be calculated
%
\begin{equation}
\begin{array}{l}
\alpha = \SI{5.5848}{\radian} = 21^\mathrm{h}19^\mathrm{m}57^\mathrm{s} \ , \\
\delta = \SI{0.8730}{\radian} = \ang{50;1;} \ .
\end{array}
\end{equation}
\section{Acknowledgements}
We want to thank Dr. J.J.M. Steijger for his careful checking of our
methods.