Simulating the scanning strategy

Stripeline includes a set of tools to simulate the observation of the sky by the telescope. The way a telescope observes the sky is called scanning strategy, and it is obviously one of the most basic and important task in any simulation of data taking.

It should be noted that in this section only the case of an ideal telescope will be analyzed. A more realistic case is described by a telescope with a fixed set of configuration angles; tools for analyzing a realistic telescope are discussed in Pointing Reconstruction Model.

Strip's location and movements

The Strip telescope, unlike most of the optical telescopes, performs a regular, uninterrupted movement using the so-called azimuth motor, a wheel that makes the whole structure spin around its gravity axis. The height of the telescope can be varied by another motor (the altitude motor), but we foresee that it will be always kept at the same angle during spinning. The following animation shows how the two motor operates: first, the altitude motor moves by roughly 20°, and then the azimuth motor starts spinning.

The red laser beam shows the so-called boresight direction, which points toward the position being observed on the celestial sphere. The motion of the telescope, coupled with the spinning of the Earth, produces a characteristic pattern in the sequence of points on the celestial sphere that are observed by Strip. To understand this pattern, you must be aware that Strip will be deployed to Tenerife, one of the Canarian Islands (Spain), a location in the Northern Emisphere, at latitude ≈28°N:

(The coordinates of the site, as well as its elevation, are available in the constants TENERIFE_LATITUDE_DEG, TENERIFE_LONGITUDE_DEG, and TENERIFE_HEIGHT_M.) Now, consider the animation of the spinning telescope as seen from space:

As above, the red vector shows the boresight direction, while the green vector shows the zenith at Tenerife (the direction along the gravitational force vector, pointings towards the sky); the yellow direction represents the direction of the beam polarizer, and it is a parameter used to convert the Stokes parameters $Q$ and $U$ to celestial coordinates. Remembering that the red beam shows where the telescope is looking, this scanning strategy observes a circumference in the sky. However, because Earth rotates once every 24 hours (much faster in the animation than in reality!), the observation pattern is a spiral that covers a «strip» on the sky, whose apparent height is roughly equal to twice the angle of the boresight wheel.

Simulating the motion of the telescope wheels

We'll now use the tools provided by Stripeline to simulate the region of the sky that is observed by the telescope. You can type the following commands in a Jupyter notebook, if you have installed the awesome IJulia package; otherwise, you can write them directly on Julia's command line. Let's load the packages we'll need to run the simulation: Plots enables the plot function, Healpix implements the Healpix tessellation algorithm on the sphere, and yours truly Stripeline:

using Plots
using Healpix
using Stripeline

Since Stripeline is a simulator, it has been designed to be extremely versatile: the angles of the two wheels discussed above (the elevation wheel and the azimuth wheel) can be specified as arbitrary functions of time, expressed in seconds:

telescope_motors(time_s) = (0.0, deg2rad(20.0), timetorotang(time_s, 1))

Here we use the handy function timetorotang to compute $2\pi\nu t$, the angle as a function of time for a continuous circular motion of 1 rotation per minute. This function must be passed as the first parameter to the function genpointings, which computes the direction and the observing angle for a number of points in time. The direction is encoded as a $N\times 2$ matrix containing the colatitude and the longitude in the two columns; we iterate on the $N$ rows of this matrix and set each pixel in a Healpix map to this value, taking advantage of Healpix.jl's function ang2pix:

function project_to_map(time_range, map)
    # We discard the second return value of "genpointings" with "_":
    # it's the orientation, but for this simple example it is not necessary
    dirs, _ = genpointings(
        telescope_motors,  # Angle of the three wheels as a function of time
        Float64[0, 0, 1],  # Observing direction in the focal plane reference frame
        time_range,        # Array of time values
    )

    # For each sample, set the corresponding pixel in the sky map to 1
    for idx in 1:length(time_range)
        # Extract the colatitude and the longitude from "dirs"
        colat, long = dirs[idx, :]
        pixel_index = ang2pix(map, colat, long)
        map[pixel_index] = 1
    end
end

It's now a matter of creating a Healpix map, wrapping the code in a call to project_to_map, and plotting the result:

# We create a Healpix map that represents the whole sky sphere,
# tessellated up to some resolution NSIDE=128
map = Healpix.HealpixMap{Float64, RingOrder}(128)

# Set the sampling time, the number of seconds used to create one
# measurement. The real value is 0.01, but we use a larger
# sampling time to make the simulation faster
sampling_time_s = 0.05

# Run a simulation lasting one minute (which corresponds to
# exactly one rotation of the telescope around the azimuth axis)
project_to_map(0.0:sampling_time_s:60.0, map)

# Plot the map, using the orthographic projection
plot(map, orthographic)

The result should match your expectations: considering the animation of the spinning telescope as seen from space, the sequence of points set to 1 in the map corresponds to the directions visited by the red vector during one spin. Let's re-run the animation over a longer time span, one hour; no need to re-allocate the map, as we'll simply overwrite it:

# Run a simulation lasting one hour (60 rotations
# of the telescope, while the Earth slowly spins)
project_to_map(0.0:sampling_time_s:3600.0, map)

plot(map, orthographic)

Running the simulation over one hour shows the effect of the Earth's rotation, which makes the circle move Eastward.

A more complex example

In the examples above, we set to 1 the pixels in the map that were observed by the Strip telescope at least one time. A much more informative way of plotting these maps is to count the number of times a pixel has been observed, as typically one wants to make several observations and then average them together to produce one estimate with its own error bar. The number stored in each pixel is the so-called hit count for that pixel, and it is related to the overall amount of time spent by the telescope observing that direction.

Let's modify the function project_to_map so that it increments the value of a pixel:

function project_to_map(time_range, map)
    dirs, _ = genpointings(telescope_motors, Float64[0, 0, 1], time_range)

    for idx in 1:length(time_range)
        colat, long = dirs[idx, :]
        pixel_index = ang2pix(map, colat, long)

        # Increment the "hit count"
        map[pixel_index] += 1
    end
end

Now let's recreate the simulation above, using the @animate macro to produce a movie of the hit count as time passes:

# Reset the map, so that each pixel is set to zero
map[:] .= 0

# Create one frame per each minute of observation
anim = @animate for minute in 0:60
    # Start and end times of this minute; note that we drop
    # the last sample from "end_time_s", as it will be included
    # in the next iteration
    start_time_s = minute * 60
    end_time_s = (minute + 1) * 60 - sampling_time_s

    project_to_map(start_time_s:sampling_time_s:end_time_s, map)

    # The keyword "clim" fixes the lower and upper limits of the
    # color bar. Avoiding to do so results in the color scale flickering
    # between frames in the animation (try it!)
    plot(map, orthographic, clim=(0, 100))
end

# Save the result into an animated GIF file, with 10 frames per second.
# As one frame in our animation lasts one minute, this means that each
# second in the animation corresponds to 10 minutes
gif(anim, "scanning-animation.gif", fps = 10)

Note that most of the pixels are observed a few times, but those on the uppermost and lowermost part of the strip have a much higher hit count. We clipped the maximum value shown in the color bar to 100, but we can make Julia compute the maximum value in the map:

println("The maximum hit count in the map is ", maximum(map))

The maximum hit count in the map is 168.0

A more complex example

So far, we have simulated the behavior of the Strip instrument in quite ideal conditions: we started each simulation from time $t = 0$, without specifying an absolute date.

The function genpointings accepts a starting time expressed as a Julian date; in this case, it uses a slower algorithm to find the direction and orientation of the telescope that considers several other effects:

1. The position of the Earth with respect to the Sun;
2. The precession of the Earth's axis;
3. The nutation of the Earth's axis due to other bodies in the Solar System;
4. Stellar aberration;
5. Atmospheric refraction, although this computation is valid only for visible wavelengths and should therefore not be used when simulating observations done by Strip detectors (which operates at microwave lengths).

To specify times, you can use the functions jdcnv and daycnv from AstroLib:

using AstroLib, Dates

# Assume that the simulation starts on January, 1st 2022, 15:00:00
start_day = DateTime(2022, 1, 1, 15, 0, 0)

println("The simulation starts from JD $start_day")

dirs, orientations = genpointings(
    telescope_motors,
    Float64[0, 0, 1],
    0:sampling_time_s:60.0,
    start_day,   # We specify here the JD when the simulation starts
)

([0.7290778648544944 5.40295952303128; 0.7290886554394445 5.405651780783916; … ; 0.7290869258580213 5.40463561231746; 0.7290852961325847 5.4073279088281945], [1.5681503232122127, 1.561224985902629, 1.5542998523707596, 1.547375068780534, 1.5404507812319306, 1.5335271356890912, 1.5266042779704168, 1.5196823537015656, 1.5127615083062913, 1.505841886934236  …  1.6304723842806206, 1.6235518176778878, 1.6166301411637374, 1.6097075000857792, 1.6027840400149536, 1.5958599067354513, 1.5889352461736386, 1.58201020438792, 1.5750849275215972, 1.5681595617930688])

Passing start_day will make genpointings use a much slower algorithm to compute the pointings; you should use this syntax only if your simulation really needs the increased precision. Typical examples where you really want to do this are the following:

1. You want to estimate when and how long an object in the sky (e.g., the Crab Nebula, Jupiter) will be visible by Strip;
2. A variation of the previous point is to compute when bright objects (Sun, Moon, etc.) might be dangerously close to the boresight direction of the telescope, in order to decide when the telescope will need to be shut down to prevent overheating or saturations in the detectors.
3. You want to produce a sky map to be compared with those of other experiments running at the same time as Strip;
4. You want to assess how much effects like precession, nutation, and aberration affect the measurements.

Reference documentation

telescopetoground(wheelanglesfn, time_s)

telescopetoground(3600.0) do
    # Boresight motor keeps a constant angle equal to 0°
    # Altitude motor remains at 20° from the Zenith
    # Ground motor spins at 1 RPM
    (0.0, deg2rad(20.0), timetorotang(time_s, 1))
end

telescopetoground(wheelanglesfn, time_s, config_ang)

groundtoearth(groundq, time_s, latitude_deg; day_duration_s=86400.0)

genpointings!(wheelanglesfn, beam_dir, timerange_s, dirs, psi;
              polaxis = Float64[1.0, 0.0, 0.0],
              latitude_deg = TENERIFE_LATITUDE_DEG,
              ground = false,
              config_ang::Union{ConfigAngles, Nothing} = nothing)
genpointings(wheelanglesfn, beam_dir, timerange_s;
             polaxis = Float64[1.0, 0.0, 0.0],
             latitude_deg = TENERIFE_LATITUDE_DEG,
             ground = false,
             config_ang::Union{ConfigAngles, Nothing} = nothing)
genpointings!(wheelanglesfn, beam_dir, timerange_s, t_start, dirs, psi;
              polaxis = Float64[1.0, 0.0, 0.0],
              latitude_deg = TENERIFE_LATITUDE_DEG,
              longitude_deg = TENERIFE_LONGITUDE_DEG,
              height_m = TENERIFE_HEIGHT_M,
              precession = true,
              nutation = true,
              aberration = true,
              refraction = true,
              config_ang::Union{ConfigAngles, Nothing} = nothing)
genpointings(wheelanglesfn, beam_dir, timerange_s, t_start;
             polaxis=Float64[1.0, 0.0, 0.0],
             latitude_deg=TENERIFE_LATITUDE_DEG,
             longitude_deg=TENERIFE_LONGITUDE_DEG,
             height_m=TENERIFE_HEIGHT_M,
             precession = true,
             nutation = true,
             aberration = true,
             refraction = true,
             config_ang::Union{ConfigAngles, Nothing} = nothing)

dir, psi = genpointings([0, 0, 1], 0:0.1:1) do time_s
    # Boresight motor keeps a constant angle equal to 0°
    # Altitude motor remains at 20° from the Zenith
    # Ground motor spins at 1 RPM
    (0.0, deg2rad(20.0), timetorotang(time_s, 1))
end

import Dates

dirs, psi = genpointings(time_s -> (0, deg2rad(20),
                                    timetorotang(time_s, 1)),
                         [0, 0, 1],
                         0:0.1:1,
                         Dates.DateTime(2022, 01, 01, 0, 0, 0),
                         latitude_deg=10.0,
                         longitude_deg=20.0,
                         height_m = 1000) do time_s

timetorotang(time, rpm)

northdir(θ, ϕ)
eastdir(θ, ϕ)

julia> northdir(π/2, 0) ≈ [0, 0, 1]
true
julia> eastdir(π/2, 0) ≈ [0, 1, 0]
true

northdir(θ, ϕ)
eastdir(θ, ϕ)

julia> northdir(π/2, 0) ≈ [0, 0, 1]
true
julia> eastdir(π/2, 0) ≈ [0, 1, 0]
true

polarizationangle(northdir, eastdir, poldir)

julia> polarizationangle([0, 0, 1], [0, 1, 0], [0, 1, 0])
0.0
julia> polarizationangle([0, 0, 1], [0, 1, 0], [0, 0, 1]) |> rad2deg
90.0
julia> polarizationangle([0, 0, 1], [0, 1, 0], [0, 0, -1]) |> rad2deg
-90.0