Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[WIP] Better scanning #44

Merged
merged 6 commits into from
Apr 12, 2020
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 8 additions & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
@@ -1,5 +1,13 @@
# HEAD

- (**Breaking change**) Drop support for Julia 1.0, 1.1, and 1.2
- (**Breaking change**) Several changes to `genpointings`
([#44](https://github.com/lspestrip/Stripeline.jl/pull/44)):
- In `genpointings`, rename `refract` keyword into `refraction`, for
consistency with `precession`, `nutation`, and `aberration`.
- Fix bug [#42](https://github.com/lspestrip/Stripeline.jl/issues/42)
- Add function `genpointings!`
- Update the coordinates of the site in Tenerife
- (**Breaking change**) Turn `rank` and `comm` parameters into
keywords for `generate_noise_mpi`
[#39](https://github.com/lspestrip/Stripeline.jl/pull/39)
Expand Down
13 changes: 13 additions & 0 deletions docs/src/scanning.md
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,19 @@ DocTestSetup = quote
end
```

```@raw html
<iframe
src="https://www.google.com/maps/embed?pb=!1m14!1m8!1m3!1d28777729.760743093!2d-34.4392714819374!3d28.30115741125523!3m2!1i1024!2i768!4f13.1!3m3!1m2!1s0x0%3A0x0!2zMjjCsDE4JzAwLjkiTiAxNsKwMzAnMzYuNCJX!5e0!3m2!1sen!2sit!4v1586614115629!5m2!1sen!2sit"
width="600"
height="450"
frameborder="0"
style="border:0;"
allowfullscreen=""
aria-hidden="false"
tabindex="0">
</iframe>
```

# Simulating the scanning strategy

```@docs
Expand Down
243 changes: 171 additions & 72 deletions src/scanning.jl
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,7 @@
# +-----+-----------------------------+------------------------------------------+
# | MCS | Mount Coordinate System | Centered at the basis of the telescope |
# +-------+---------------------------+------------------------------------------+
# | ACS | A zimuth Coordinate System | Moveable, centered in the middle of |
# | ACS | Azimuth Coordinate System | Moveable, centered in the middle of |
# | | | ground motor. Positive angles go from |
# | | | North to *East* (beware!) |
# +-----+-----------------------------+------------------------------------------+
Expand Down Expand Up @@ -80,10 +80,15 @@ import Dates

export TENERIFE_LATITUDE_DEG, TENERIFE_LONGITUDE_DEG, TENERIFE_HEIGHT_M
export timetorotang, telescopetoground, groundtoearth
export genpointings, northdir, eastdir, polarizationangle
export genpointings!, genpointings, northdir, eastdir, polarizationangle

const TENERIFE_LATITUDE_DEG = 28.3
const TENERIFE_LONGITUDE_DEG = -16.509722
"Latitude of the LSPE/Strip site in Tenerife, in degrees"
const TENERIFE_LATITUDE_DEG = 28.30026

"Longitude of the LSPE/Strip site in Tenerife, in degrees"
const TENERIFE_LONGITUDE_DEG = -16.51012

"Height of the LSPE/Strip site in Tenerife, in meters"
const TENERIFE_HEIGHT_M = 2390

include("quaternions.jl")
Expand Down Expand Up @@ -275,67 +280,100 @@ function quat_to_angles(boreaxis, polaxis, quat)
end


function genpointings(wheelanglesfn,
dir,
timerange_s;
polaxis = Float64[1.0, 0.0, 0.0],
latitude_deg = 0.0,
ground = false,
day_duration_s = 86400.0)
function genpointings!(wheelanglesfn,
beam_dir,
timerange_s,
dirs,
psi;
polaxis = Float64[1.0, 0.0, 0.0],
latitude_deg = TENERIFE_LATITUDE_DEG,
ground = false,
day_duration_s = 86400.0)

if ground
dirs = Array{Float64}(undef, length(timerange_s), 4)
ψ = Array{Float64}(undef, length(timerange_s), 2)
@assert size(dirs, 2) == 4
@assert size(psi, 2) == 2
else
dirs = Array{Float64}(undef, length(timerange_s), 2)
ψ = Array{Float64}(undef, length(timerange_s))
@assert size(dirs, 2) == 2
@assert size(psi, 2) == 1
end
@assert size(dirs, 1) == size(psi, 1)

for (idx, time_s) = enumerate(timerange_s)

# This converts the RDP into the MCS (ground reference frame)
groundq = telescopetoground(wheelanglesfn, time_s)
# This converts the MCS into the celestial reference frame
quat = groundtoearth(groundq, time_s, latitude_deg; day_duration_s = day_duration_s)

θ, ϕ, curψ = quat_to_angles(dir, polaxis, quat)
θ, ϕ, curpsi = quat_to_angles(beam_dir, polaxis, quat)
(dirs[idx, 1], dirs[idx, 2]) = (θ, ϕ)

if ground
# Re-run the transformation algorithm using the ground quaternion
θ_ground, ϕ_ground, ψ_ground = quat_to_angles(dir, polaxis, groundq)
θ_ground, ϕ_ground, psi_ground = quat_to_angles(beam_dir, polaxis, groundq)

(dirs[idx, 3], dirs[idx, 4]) = (θ_ground, ϕ_ground)
(ψ[idx, 1], ψ[idx, 2]) = (curψ, ψ_ground)
(psi[idx, 1], psi[idx, 2]) = (curpsi, psi_ground)
else
ψ[idx] = curψ
psi[idx] = curpsi
end
end

(dirs, ψ)
end


function genpointings(wheelanglesfn,
dir,
timerange_s,
t_start;
polaxis = Float64[1.0, 0.0, 0.0],
latitude_deg = 0.0,
longitude_deg = 0.0,
height_m = 0,
precession = true,
nutation = true,
aberration = true,
refract = true)
beam_dir,
timerange_s;
polaxis = Float64[1.0, 0.0, 0.0],
latitude_deg = TENERIFE_LATITUDE_DEG,
ground = false,
day_duration_s = 86400.0)

skydirs = Array{Float64}(undef, length(timerange_s), 2)
skyψ = Array{Float64}(undef, length(timerange_s))
if ground
dirs = Array{Float64}(undef, length(timerange_s), 4)
psi = Array{Float64}(undef, length(timerange_s), 2)
else
dirs = Array{Float64}(undef, length(timerange_s), 2)
psi = Array{Float64}(undef, length(timerange_s))
end

genpointings!(
wheelanglesfn,
beam_dir,
timerange_s,
dirs,
psi;
polaxis = polaxis,
latitude_deg = latitude_deg,
ground = ground,
day_duration_s = day_duration_s,
)

(dirs, psi)
end

function genpointings!(wheelanglesfn,
beam_dir,
timerange_s,
t_start,
skydirs,
skypsi;
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)

@assert size(skydirs, 1) == size(skypsi, 1)
@assert size(skydirs, 2) == 2
@assert size(skypsi, 2) == 1

for (idx, time_s) = enumerate(timerange_s)
groundq = telescopetoground(wheelanglesfn, time_s)
rotmatr = rotationmatrix_normalized(groundq)
vector = rotmatr * dir
vector = rotmatr * beam_dir

jd = AstroLib.jdcnv(t_start + Dates.Nanosecond(round(Int64, time_s * 1e9)))
Dec_rad, Ra_rad = vector2equatorial(vector,
Expand All @@ -346,33 +384,92 @@ function genpointings(wheelanglesfn,
precession,
nutation,
aberration,
refract)
refraction)

poldir = rotmatr * polaxis
north = northdir(π / 2 - Dec_rad, Ra_rad)
east = eastdir(π / 2 - Dec_rad, Ra_rad)

skydirs[idx, 1] = Dec_rad
skydirs[idx, 2] = Ra_rad
skyψ[idx] = polarizationangle(north, east, poldir)
skypsi[idx] = polarizationangle(north, east, poldir)
end
end

function 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 = false)

(skydirs, skyψ)
skydirs = Array{Float64}(undef, length(timerange_s), 2)
skypsi = Array{Float64}(undef, length(timerange_s))

genpointings!(
wheelanglesfn,
beam_dir,
timerange_s,
t_start,
skydirs,
skypsi,
polaxis = polaxis,
latitude_deg = latitude_deg,
longitude_deg = longitude_deg,
height_m = height_m,
precession = precession,
nutation = nutation,
aberration = aberration,
refraction = refraction,
)

(skydirs, skypsi)
end


@doc raw"""
genpointings(wheelanglesfn, dir, timerange_s;
polaxis=Float64[1.0, 0.0, 0.0],
latitude_deg=0.0, ground=false)
genpointings(wheelanglesfn, dir, timerange_s, t_start;
genpointings!(wheelanglesfn, beam_dir, timerange_s, dirs, psi;
polaxis = Float64[1.0, 0.0, 0.0],
latitude_deg = TENERIFE_LATITUDE_DEG,
ground = false)
genpointings(wheelanglesfn, beam_dir, timerange_s;
polaxis = Float64[1.0, 0.0, 0.0],
latitude_deg = TENERIFE_LATITUDE_DEG,
ground = false)
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)
genpointings(wheelanglesfn, beam_dir, timerange_s, t_start;
polaxis=Float64[1.0, 0.0, 0.0],
latitude_deg=0.0, longitude_deg=0.0, height_m=0.0)

Generate a set of pointings for a STRIP detector. The parameter
`wheelanglesfn` must be a function which takes as input a time in
seconds and returns a 3-tuple containing the angles (in radians) of
the three motors:
latitude_deg=TENERIFE_LATITUDE_DEG,
longitude_deg=TENERIFE_LONGITUDE_DEG,
height_m=TENERIFE_HEIGHT_M,
precession = true,
nutation = true,
aberration = true,
refraction = true)

Generate a set of pointing directions for a STRIP detector. Each
function is provided in two flavours: the ones ending with `!` save
the results in the last two parameters `dirs` and `psi`, while the
others automatically allocate memory and return their results as a
pair `(dirs, psi)`.

The parameter `wheelanglesfn` must be a function which takes as input
a time in seconds and returns a 3-tuple containing the angles (in
radians) of the three motors:

1. The boresight motor (rotation around the ``z`` axis,
counterclockwise)
Expand All @@ -385,7 +482,7 @@ the three motors:

The meaning of the parameters/keywords is the following:

- `dir` specifies the pointing direction of the mean (the boresight is
- `beam_dir` specifies the pointing direction of the mean (the boresight is
[0, 0, 1]). It must be normalized.

- `timerange_s` is an enumerable type that specifies at which times
Expand All @@ -404,46 +501,48 @@ The meaning of the parameters/keywords is the following:

- `polaxis` is the polarization axis; it must be normalized.

- `precession` is a boolen parameter that allow to include the Earth's precession
effect; default is true.
- `precession`: if `true` (the default), the Earth's precession is
taken into account.

- `nutation` is a boolen parameter that allow to include the Earth's nutation
effect; default is true.
- `nutation`: if `true` (the default), the Earth's nutation is taken
into account.

- `aberration` is a boolen parameter that allow to include the stellar aberration
effect; default is true.
- `aberration`: if `true` (the default), stellar aberration is taken
into account.

- `refract` is a boolen parameter that allow to include the refraction correction
; default is true.
- `refraction`: if `true`, refraction corrections are taken into account.
As these corrections are only valid for optical wavelengths, the
default is `false`.

# Return values

In the first form, the function returns a 2-tuple containing the sky
directions (a N×2 array containing declination and right ascension, in
Equatorial coordinates) and the polarization angle for each time step.
If `t_start` is not provided, the function `genpointings` returns a
2-tuple containing the sky directions (a N×2 array containing
declination and right ascension, in Equatorial coordinates) and the
polarization angle for each time step. The function `genpointings!` sets
the values in the last two parameters `dirs` and `psi`.

In the second form, the function returns a 2-tuple (4-tuple)
containing the directions (a N×2 or Nx4 array containing the
colatitude and the longitude) and the polarization angles at each time
step.
If `t_start` is provided, the function `genpointings` returns a
2-tuple (4-tuple) containing the directions (a N×2 or Nx4 array
containing the colatitude and the longitude) and the polarization
angles at each time step; `genpointings!` works as above.


# Examples

Here is an example of the first kind:
Here is an example using the form without `t_start`:

`````julia
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
return (0.0, deg2rad(20.0), timetorotang(time_s, 1))
(0.0, deg2rad(20.0), timetorotang(time_s, 1))
end
`````

And here is an example of the second kind. Note that, unlike the
previous example, we use a lambda function instead of the `do...end`
notation.
And here is an example using `t_start`; unlike the previous example,
we use a lambda function instead of the `do...end` notation.

`````julia
import Dates
Expand All @@ -452,7 +551,7 @@ dirs, psi = genpointings(time_s -> (0, deg2rad(20),
timetorotang(time_s, 1)),
[0, 0, 1],
0:0.1:1,
Dates.DateTime(2019, 01, 01, 0, 0, 0),
Dates.DateTime(2022, 01, 01, 0, 0, 0),
latitude_deg=10.0,
longitude_deg=20.0,
height_m = 1000) do time_s
Expand Down
Loading