/
flat_fftgrid.jl
68 lines (50 loc) · 2.37 KB
/
flat_fftgrid.jl
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
# derivatives can either be applied in fourier space by multipliying by im*k or
# in map space by finite differences. this tracks which option to use for a
# given field
abstract type ∂modes end
struct fourier∂ <: ∂modes end
struct map∂ <: ∂modes end
promote_rule(::Type{map∂}, ::Type{fourier∂}) = fourier∂
# Flat{Nside,θpix,∂mode} is a flat sky pixelization with `Nside` pixels per side
# and pixels of width `θpix` arcmins, where derivatives are done according to ∂mode
abstract type Flat{Nside,θpix,∂mode<:∂modes} <: Pix end
# for convenience
Flat(;Nside, θpix=θpix₀, ∂mode=fourier∂) = Flat{Nside,θpix,∂mode}
Nside(::Type{P}) where {N,P<:Flat{N}} = N
# default angular resolution used by a number of convenience constructors
θpix₀ = 1
@doc """
The number of FFTW threads to use. This must be set via e.g.:
CMBLensing.FFTW_NUM_THREADS[] = 4
*before* creating any `FlatField` objects; subsequent changes to this variable
will be ignored. The default value is the environment variable `FFTW_NUM_THREADS`,
or if that is not specified its `Sys.CPU_THREADS÷2`.
"""
const FFTW_NUM_THREADS = Ref{Int}()
@init FFTW_NUM_THREADS[] = parse(Int,get(ENV,"FFTW_NUM_THREADS","$(Sys.CPU_THREADS÷2)"))
@generated function FlatInfo(::Type{T}, ::Type{Arr}, ::Val{θpix}, ::Val{Nside}) where {T<:Real, Arr<:AbstractArray, θpix, Nside}
FFTW.set_num_threads(FFTW_NUM_THREADS[])
Δx = T(deg2rad(θpix/60))
FFT = plan_rfft(Arr{T}(undef,Nside,Nside))
Δℓ = T(2π/(Nside*Δx))
nyq = T(2π/(2Δx))
Ωpix = T(Δx^2)
x,k = (ifftshift(-Nside÷2:(Nside-1)÷2),) .* (Δx,Δℓ)
kmag = @. sqrt(k'^2 + k^2)
ϕ = @. angle(k' + im*k)[1:Nside÷2+1,:]
sin2ϕ, cos2ϕ = @. sin(2ϕ), cos(2ϕ)
if iseven(Nside)
sin2ϕ[end, end:-1:(Nside÷2+2)] .= sin2ϕ[end, 2:Nside÷2]
end
@namedtuple(T, θpix, Nside, Δx, Δℓ, nyq, Ωpix, x, k, kmag, sin2ϕ=Arr(sin2ϕ), cos2ϕ=Arr(cos2ϕ), FFT)
end
function Cℓ_to_2D(::Type{P}, ::Type{T}, Cℓ) where {T,N,P<:Flat{N}}
Complex{T}.(nan2zero.(Cℓ.(fieldinfo(P,T).kmag[1:N÷2+1,:])))
end
@doc doc"""
pixwin(θpix, ℓ)
Returns the pixel window function for square flat-sky pixels of width `θpix` (in
arcmin) evaluated at some `ℓ`s. This is the scaling of k-modes, the scaling of
the power spectrum will be pixwin^2.
"""
pixwin(θpix, ℓ) = @. sinc(ℓ*deg2rad(θpix/60)/2π)