Module s2dw_core_mod

module s2dw_core_mod

        ! Uses
    use s2dw_types_mod
    use s2dw_error_mod
    use s2dw_dl_mod
    use omp_lib

        ! Types
    public type s2dw_wav_abg

        ! Variables
    integer, private, parameter :: FFTW_ESTIMATE = 64
    integer, private, parameter :: FFTW_MEASURE = 0
    integer, private, parameter :: FFTW_FORWARD = -1
    integer, private, parameter :: FFTW_BACKWARD = 1

        ! Subroutines and functions
    public subroutine s2dw_core_init_kernels (K_gamma, Phi2, bl_scoeff, J, B, alpha)
    private function s2dw_core_phi2_integrand (k, bw, alpha) result (integrand)
    public subroutine s2dw_core_init_directionality (Slm, B, N)
    public function s2dw_core_admiss (admiss, K_gamma, Phi2, B, J) result (pass)
    public subroutine s2dw_core_analysis_flm2wav (wav, scoeff, flm, K_gamma, Slm, J, B, N, bl_scoeff, alpha)
    public subroutine s2dw_core_analysis_flm2wav_dynamic (wavdyn, scoeff, flm, K_gamma, Slm, J, B, N, bl_scoeff, alpha)
    public subroutine s2dw_core_synthesis_wav2flm (flm, wav, scoeff, K_gamma, Phi2, Slm, J, B, N, bl_scoeff, alpha)
    public subroutine s2dw_core_synthesis_wav2flm_dynamic (flm, wavdyn, scoeff, K_gamma, Phi2, Slm, J, B, N, bl_scoeff, alpha)
    public subroutine s2dw_core_free_wavdyn (wavdyn)
    private function s2dw_core_params_valid (J, B, alpha, N) result (valid)
    public function s2dw_core_comp_Jmax (B, alpha) result (J_max)
    public function s2dw_core_assimilate_int (x, y) result (assimilated)
    private function quad_weights (theta_t, B) result (w)
    private function binomial (n, r) result (b)
    private function logfact (n) result (logfactn)
    private subroutine trapzd (func, B, alpha, a_lim, b_lim, s, n)
    private function qtrap (func, B, alpha, a_lim, b_lim) result (integral)
    private function qsimp (func, B, alpha, a_lim, b_lim) result (integral)

end module s2dw_core_mod

Functionality to perform scale discretised wavelet transform on the sphere.

Author: J. D. McEwen (mcewen@mrao.cam.ac.uk)

Version: 0.1 November 2007

Description of Types

s2dw_wav_abg

public type s2dw_wav_abg
    real (kind=dp), allocatable, dimension (:,:,:) :: coeff
end type s2dw_wav_abg

Structure to store wavelet coefficients for given scale. Required for dynamic memory allocation of wavelet coefficients.

Description of Variables

FFTW_ESTIMATE

integer, private, parameter :: FFTW_ESTIMATE = 64

FFTW_MEASURE

integer, private, parameter :: FFTW_MEASURE = 0

FFTW_FORWARD

integer, private, parameter :: FFTW_FORWARD = -1

FFTW_BACKWARD

integer, private, parameter :: FFTW_BACKWARD = 1

Description of Subroutines and Functions

s2dw_core_init_kernels

public subroutine s2dw_core_init_kernels (K_gamma, Phi2, bl_scoeff, J, B, alpha)
    real (kind=dp), intent(out), dimension (0:J,0:B-1) :: K_gamma
    real (kind=dp), intent(out), dimension (0:B-1) :: Phi2
    integer, intent(out) :: bl_scoeff
    integer, intent(in) :: J
    integer, intent(in) :: B
    real (kind=dp), intent(in) :: alpha
    ! Calls: s2dw_error
end subroutine s2dw_core_init_kernels

Compute kernels, scaling functions and band limits.

Notes:

Kernels are computed from differences of scaling functions, which are computed directly. It is important to compute kernels in this manner to ensure that any errors in the numerical integration collapse when computing the resolution of the identity.

Variables:

K_gamma(0:J, 0:B-1): Kernel functions computed for each scale [output].
Phi2(0:B-1): Scaling function squared at maximum j=J [output].
bl_scoeff: Upper band limit for scaling coefficients [output].
J: Maximum analysis scale depth [input].
B: Harmonic band limit [input].

Author: J. D. McEwen

Version: 0.1 October 2007

s2dw_core_phi2_integrand

private function s2dw_core_phi2_integrand (k, bw, alpha) result (integrand)
    real (kind=dp), intent(in) :: k
    integer, intent(in) :: bw
    real (kind=dp), intent(in) :: alpha
    real (kind=dp) :: integrand
end function s2dw_core_phi2_integrand

Integrand of definite integral to be calculated to compute the scaling function.

Notes:

Integrand is given by K_Phi^2(k')/k'.

Variables:

k: Dummy (scale) variable of integration [input].
B: Harmonic band limit [input].
alpha: Basis dilation factor [input].
integrand: Value of integrand computed [output].

Author: J. D. McEwen

Version: 0.1 October 2007

s2dw_core_init_directionality

public subroutine s2dw_core_init_directionality (Slm, B, N)
    complex (kind=dpc), intent(out), dimension (0:B-1,0:N-1) :: Slm
    integer, intent(in) :: B
    integer, intent(in) :: N
    ! Calls: s2dw_error
end subroutine s2dw_core_init_directionality

Compute wavelet directionality coefficients.

Variables:

Slm(0:B-1, 0:N-1): Directionality coefficients [output].
B: Harmonic band limit [input].
N: Azimuthal band limit [input].

Author: J. D. McEwen

Version: 0.1 October 2007

s2dw_core_admiss

public function s2dw_core_admiss (admiss, K_gamma, Phi2, B, J) result (pass)
    real (kind=dp), intent(out), dimension (0:B-1) :: admiss
    real (kind=dp), intent(in), dimension (0:J,0:B-1) :: K_gamma
    real (kind=dp), intent(in), dimension (0:B-1) :: Phi2
    integer, intent(in) :: B
    integer, intent(in) :: J
    logical :: pass
end function s2dw_core_admiss

Compute resolution of the identity from scaled kernels and scaling function and check equal to unity for all el.

Notes:

Resolution of identity must be unity for all el to ensure admissibility satisfied, i.e. to ensure exact reconstruction is possible.
Resolution of identity computed by Phi2(el) + sum_{jj=0}^{J} K_gamma^2(jj,el)

Variables:

admiss(0:B-1): Resolution of the identity computed for each el (all values should be close to one if passed admissibility test) [output].
K_gamma(0:J, 0:B-1): Kernel functions computed for each scale [input].
Phi2(0:B-1): Scaling function squared at maximum j=J [input].
J: Maximum analysis scale depth [input].
B: Harmonic band limit [input].
pass: Logical specifying whether kernels and scaling function passed admissibilty test [output].

Author: J. D. McEwen

Version: 0.1 October 2007

s2dw_core_analysis_flm2wav

public subroutine s2dw_core_analysis_flm2wav (wav, scoeff, flm, K_gamma, Slm, J, B, N, bl_scoeff, alpha)
    real (kind=dp), intent(out), dimension (0:J, 0:2*B-2, 0:2*B-1, 0:N-1) :: wav
    complex (kind=dpc), intent(out), dimension (0:bl_scoeff-1, 0:bl_scoeff-1) :: scoeff
    complex (kind=dpc), intent(in), dimension (0:B-1, 0:B-1) :: flm
    real (kind=dp), intent(in), dimension (0:J, 0:B-1) :: K_gamma
    complex (kind=dpc), intent(in), dimension (0:B-1, 0:N-1) :: Slm
    integer, intent(in) :: J
    integer, intent(in) :: B
    integer, intent(in) :: N
    integer, intent(in) :: bl_scoeff
    real (kind=dp), intent(in) :: alpha
    ! Calls: dfftw_destroy_plan, dfftw_execute_dft_c2r, dfftw_plan_dft_c2r_1d, s2dw_dl_beta_operator, s2dw_error
end subroutine s2dw_core_analysis_flm2wav

Compute wavelet and scaling coefficients from harmonic coefficients of original signal using fast algorithm.

Notes:

Wavelet and scaling coefficients contain all information required to reconstruct the original signal exactly.

Variables:

wav(0:J, 0:2*B-2, 0:2*B-1, 0:N-1): Wavelet coefficients [output].
scoeff(0:bl_scoeff-1, 0:bl_scoeff-1): Scaling coefficients [output].
flm(0:B-1, 0:B-1): Harmonic coefficients of signal [input].
K_gamma(0:J, 0:B-1): Kernel functions computed for each scale [input].
Slm(0:B-1, 0:N-1): Directionality coefficients [input].
J: Maximum analysis scale depth [input].
B: Harmonic band limit [input].
N: Azimuthal band limit [input].
bl_scoeff: Upper band limit for scaling coefficients [input].
alpha: Basis dilation factor [input].

Author: J. D. McEwen

Version: 0.1 October 2007

s2dw_core_analysis_flm2wav_dynamic

public subroutine s2dw_core_analysis_flm2wav_dynamic (wavdyn, scoeff, flm, K_gamma, Slm, J, B, N, bl_scoeff, alpha)
    type (s2dw_wav_abg), intent(out), allocatable, dimension (:) :: wavdyn
    complex (kind=dpc), intent(out), dimension (0:bl_scoeff-1, 0:bl_scoeff-1) :: scoeff
    complex (kind=dpc), intent(in), dimension (0:B-1, 0:B-1) :: flm
    real (kind=dp), intent(in), dimension (0:J, 0:B-1) :: K_gamma
    complex (kind=dpc), intent(in), dimension (0:B-1, 0:N-1) :: Slm
    integer, intent(in) :: J
    integer, intent(in) :: B
    integer, intent(in) :: N
    integer, intent(in) :: bl_scoeff
    real (kind=dp), intent(in) :: alpha
    ! Calls: dfftw_destroy_plan, dfftw_execute_dft_c2r, dfftw_plan_dft_c2r_1d, s2dw_dl_beta_operator, s2dw_error
end subroutine s2dw_core_analysis_flm2wav_dynamic

Compute wavelet and scaling coefficients from harmonic coefficients of original signal using fast algorithm, dynamically allocating the memory required for the wavelet coefficients dynamically.

Notes:

The memory required to store the wavelet coefficients is allocated dynamically herein.
Wavelet and scaling coefficients contain all information required to reconstruct the original signal exactly.

Variables:

wavdyn(0:J)%coeff: Wavelet coefficients for each scale (memory allocated herein) [output].
scoeff(0:bl_scoeff-1, 0:bl_scoeff-1): Scaling coefficients [output].
flm(0:B-1, 0:B-1): Harmonic coefficients of signal [input].
K_gamma(0:J, 0:B-1): Kernel functions computed for each scale [input].
Slm(0:B-1, 0:N-1): Directionality coefficients [input].
J: Maximum analysis scale depth [input].
B: Harmonic band limit [input].
N: Azimuthal band limit [input].
bl_scoeff: Upper band limit for scaling coefficients [input].
alpha: Basis dilation factor [input].

Author: J. D. McEwen

Version: 0.1 February 2008

s2dw_core_synthesis_wav2flm

public subroutine s2dw_core_synthesis_wav2flm (flm, wav, scoeff, K_gamma, Phi2, Slm, J, B, N, bl_scoeff, alpha)
    complex (kind=dpc), intent(out), dimension (0:B-1, 0:B-1) :: flm
    real (kind=dp), intent(in), dimension (0:J, 0:2*B-2, 0:2*B-1, 0:N-1) :: wav
    complex (kind=dpc), intent(in), dimension (0:bl_scoeff-1, 0:bl_scoeff-1) :: scoeff
    real (kind=dp), intent(in), dimension (0:J, 0:B-1) :: K_gamma
    real (kind=dp), intent(in), dimension (0:B-1) :: Phi2
    complex (kind=dpc), intent(in), dimension (0:B-1, 0:N-1) :: Slm
    integer, intent(in) :: J
    integer, intent(in) :: B
    integer, intent(in) :: N
    integer, intent(in) :: bl_scoeff
    real (kind=dp), intent(in) :: alpha
    ! Calls: dfftw_destroy_plan, dfftw_execute_dft_r2c, dfftw_plan_dft_r2c_1d, s2dw_dl_beta_operator, s2dw_error
end subroutine s2dw_core_synthesis_wav2flm

Synthesis harmonic coefficients of signal from wavelet and scaling coefficients.

Notes:

This routines uses less memory than using s2dw_core_synthesis_wav2wig followed by s2dw_core_synthesis_wig2flm since the Wigner coefficients do not need to be stored for all j.

Variables:

flm(0:B-1, 0:B-1): Synthesised harmonic coefficients of signal [output].
wav(0:J, 0:2*B-2, 0:2*B-1, 0:N-1): Wavelet coefficients [input].
scoeff(0:bl_scoeff-1, 0:bl_scoeff-1): Scaling coefficients [input].
K_gamma(0:J, 0:B-1): Kernel functions computed for each scale [input].
Phi2(0:B-1): Scaling function squared at maximum j=J [input].
Slm(0:B-1, 0:N-1): Directionality coefficients [input].
J: Maximum analysis scale depth [input].
B: Harmonic band limit [input].
N: Azimuthal band limit [input].
bl_scoeff: Upper band limit for scaling coefficients [input].
alpha: Basis dilation factor [input].

Author: J. D. McEwen

Version: 0.1 October 2007

s2dw_core_synthesis_wav2flm_dynamic

public subroutine s2dw_core_synthesis_wav2flm_dynamic (flm, wavdyn, scoeff, K_gamma, Phi2, Slm, J, B, N, bl_scoeff, alpha)
    complex (kind=dpc), intent(out), dimension (0:B-1, 0:B-1) :: flm
    type (s2dw_wav_abg), intent(inout), allocatable, dimension (:) :: wavdyn
    complex (kind=dpc), intent(in), dimension (0:bl_scoeff-1, 0:bl_scoeff-1) :: scoeff
    real (kind=dp), intent(in), dimension (0:J, 0:B-1) :: K_gamma
    real (kind=dp), intent(in), dimension (0:B-1) :: Phi2
    complex (kind=dpc), intent(in), dimension (0:B-1, 0:N-1) :: Slm
    integer, intent(in) :: J
    integer, intent(in) :: B
    integer, intent(in) :: N
    integer, intent(in) :: bl_scoeff
    real (kind=dp), intent(in) :: alpha
    ! Calls: dfftw_destroy_plan, dfftw_execute_dft_r2c, dfftw_plan_dft_r2c_1d, s2dw_dl_beta_operator, s2dw_error
end subroutine s2dw_core_synthesis_wav2flm_dynamic

Synthesise harmonic coefficients of signal from wavelet and scaling coefficients, where the memory of the wavelet coefficients has been allocated dynamically (destroying wavelet coefficients in the process).

Notes:

The wavelet coefficients are destroyed by this routine and their memory is freed.
This routines uses less memory than using s2dw_core_synthesis_wav2wig followed by s2dw_core_synthesis_wig2flm since the Wigner coefficients do not need to be stored for all j.

Variables:

flm(0:B-1, 0:B-1): Synthesised harmonic coefficients of signal [output].
wavdyn(0:J)%coeff: Wavelet coefficients for each scale (destroyed and freed on output) [input/output].
scoeff(0:bl_scoeff-1, 0:bl_scoeff-1): Scaling coefficients [input].
K_gamma(0:J, 0:B-1): Kernel functions computed for each scale [input].
Phi2(0:B-1): Scaling function squared at maximum j=J [input].
Slm(0:B-1, 0:N-1): Directionality coefficients [input].
J: Maximum analysis scale depth [input].
B: Harmonic band limit [input].
N: Azimuthal band limit [input].
bl_scoeff: Upper band limit for scaling coefficients [input].
alpha: Basis dilation factor [input].

Author: J. D. McEwen

Version: 0.1 February 2008

s2dw_core_free_wavdyn

public subroutine s2dw_core_free_wavdyn (wavdyn)
    type (s2dw_wav_abg), intent(inout), allocatable, dimension (:) :: wavdyn
end subroutine s2dw_core_free_wavdyn

Free memory corresponding to dynamically allocated wavelet coefficients.

Variables:

wavdyn(0:J)%coeff: Wavelet coefficients for each scale (destroyed and freed on output) [input/output].

Author: J. D. McEwen

Version: 0.1 February 2008

s2dw_core_params_valid

private function s2dw_core_params_valid (J, B, alpha, N) result (valid)
    integer, intent(in) :: J
    integer, intent(in) :: B
    real (kind=dp), intent(in) :: alpha
    integer, optional, intent(in) :: N
    logical :: valid
end function s2dw_core_params_valid

Check size parameters valid.

Notes:

Only needs to be called once when initialising Kernels, then array sizes ensure other objects of correct size.

Variables:

J: Maximum analysis scale depth [input].
B: Harmonic band limit [input].
N: Azimuthal band limit [input].
valid: Logical specifying whether passed parameter check [output].

Author: J. D. McEwen

Version: 0.1 October 2007

s2dw_core_comp_Jmax

public function s2dw_core_comp_Jmax (B, alpha) result (J_max)
    integer, intent(in) :: B
    real (kind=dp), intent(in) :: alpha
    integer :: J_max
end function s2dw_core_comp_Jmax

Compute maximum analysis depth allowed.

Variables:

B: Harmonic band limit [input].
J_max: Maximum analysis scale depth allowed [output].

Author: J. D. McEwen

Version: 0.1 November 2007

s2dw_core_assimilate_int

public function s2dw_core_assimilate_int (x, y) result (assimilated)
    real (kind=dp), intent(in) :: x
    integer, intent(out) :: y
    logical :: assimilated
end function s2dw_core_assimilate_int

If a real value is very close to an integer then set to integer, otherwise leave unaltered. Necessary to avoid numerical precision problems when using ceiling and floor functions.

Variables:

x: Input value to (possibly) assimilate to integer [input].
y: Value x is assimilated to (or 0 if not assimilated) [output].
assimilated: Logical indicating where assimilation occured [output].

Author: J. D. McEwen

Version: 0.1 November 2007

quad_weights

private function quad_weights (theta_t, B) result (w)
    real (kind=dp), intent(in) :: theta_t
    integer, intent(in) :: B
    real (kind=dp) :: w
end function quad_weights

Compute quadrature weights for exact quadrature with measure dcos(theta). Weights are derived by Driscoll and Healy.

Variables:

theta_t: Theta value to compute weight for [input].
B: Harmonic band limit [input].
w: Corresponding weight [output]

Author: J. D. McEwen

Version: 0.1 October 2007

binomial

private function binomial (n, r) result (b)
    integer, intent(in) :: n
    integer, intent(in) :: r
    integer :: b
end function binomial

Compute Binomial coefficient nCr.

Variables:

n: Upper Binomial coefficient index [input].
r: Lower Binomial coefficient index [input].
b: Binomial coefficient value [output].

Author: J. D. McEwen

Version: 0.1 October 2007

logfact

private function logfact (n) result (logfactn)
    integer, intent(in) :: n
    real (kind=dp) :: logfactn
    ! Calls: s2dw_error
end function logfact

Computes the natural logarithm of an (integer) factorial.

Variables:

n: Integer to compute factorial of [input].
logfactn: Natural logarithm of factorial value computed [output].

Author: J. D. McEwen

Version: 0.1 October 2007

trapzd

private subroutine trapzd (func, B, alpha, a_lim, b_lim, s, n)
    interface func
        function func (k, B, alpha) result (integrand)
            real (kind=dp), intent(in) :: k
            integer, intent(in) :: B
            real (kind=dp), intent(in) :: alpha
            real (kind=dp) :: integrand
        end function func
    end interface func
    integer, intent(in) :: B
    real (kind=dp), intent(in) :: alpha
    real (kind=dp), intent(IN) :: a_lim
    real (kind=dp), intent(IN) :: b_lim
    real (kind=dp), intent(INOUT) :: s
    integer, intent(IN) :: n
end subroutine trapzd

Computes nth stage of refinement of extended trapezoidal rule. Adapted from numerical recipes for computing scaling functions squared.

Notes:

Numerical recipies comment: This routine computes the nth stage of refinement of an extended trapezoidal rule. func is input as the name of the function to be integrated between limits a and b, also input. When called with n=1, the routine returns as s the crudest estimate of int_b^a f(x)dx. Subsequent calls with n=2,3,... (in that sequential order) will improve the accuracy of s by adding 2n-2 additional interior points. s should not be modified between sequential calls.
Adapted for use in evaluating integrals required to compute scaling functions squared (i.e. also takes parameters required to compute these integrals).

Variables:

func: "Pointer" to integrand function [input].
B: Harmonic band limit parameter for integrand function [input].
alpha:Basis dilation factor for integrand function [input].
a_lim: Lower limit to evalue definite integral for [input].
b_lim: Upper limit to evalue definite integral for [input].
s: Value of integral computed to current order [input].
n: Refinement order [output].

Author: J. D. McEwen

Version: 0.1 June 2005

qtrap

private function qtrap (func, B, alpha, a_lim, b_lim) result (integral)
    interface func
        function func (k, B, alpha) result (integrand)
            real (kind=dp), intent(in) :: k
            integer, intent(in) :: B
            real (kind=dp), intent(in) :: alpha
            real (kind=dp) :: integrand
        end function func
    end interface func
    integer, intent(in) :: B
    real (kind=dp), intent(in) :: alpha
    real (kind=dp), intent(in) :: a_lim
    real (kind=dp), intent(in) :: b_lim
    real (kind=dp) :: integral
    ! Calls: s2dw_error, trapzd
end function qtrap

Computes the integral of the function func from a to b using the trapezoid method. Adapted from numerical recipes for computing scaling functions squared.

Notes:

Numerical recipies comment: Returns the integral of the function func from a to b. The parameter EPS should be set to the desired fractional accuracy and JMAX so that 2 to the power JMAX-1 is the maximum allowed number of steps. Integration is performed by the trapezoidal rule.
Adapted for use in evaluating integrals required to compute scaling functions squared (i.e. also takes parameters required to compute these integrals).

Variables:

func: "Pointer" to integrand function [input].
B: Harmonic band limit parameter for integrand function [input].
alpha:Basis dilation factor for integrand function [input].
a_lim: Lower limit to evalue definite integral for [input].
b_lim: Upper limit to evalue definite integral for [input].
integral: Value of the evaluated integral [output].

Author: J. D. McEwen

Version: 0.1 June 2005

qsimp

private function qsimp (func, B, alpha, a_lim, b_lim) result (integral)
    interface func
        function func (k, B, alpha) result (integrand)
            real (kind=dp), intent(in) :: k
            integer, intent(in) :: B
            real (kind=dp), intent(in) :: alpha
            real (kind=dp) :: integrand
        end function func
    end interface func
    integer, intent(in) :: B
    real (kind=dp), intent(in) :: alpha
    real (kind=dp), intent(IN) :: a_lim
    real (kind=dp), intent(IN) :: b_lim
    real (kind=dp) :: integral
    ! Calls: s2dw_error, trapzd
end function qsimp

Computes the integral of the function func from a to b using Simpson's rule. Adapted from numerical recipes for computing scaling functions squared.

Notes:

Numerical recipies comment: Returns the integral of the function func from a to b. The parameter EPS should be set to the desired fractional accuracy and JMAX so that 2 to the power JMAX-1 is the maximum allowed number of steps. Integration is performed by Simpson's rule.
Adapted for use in evaluating integrals required to compute scaling functions squared (i.e. also takes parameters required to compute these integrals).

Variables:

func: "Pointer" to integrand function [input].
B: Harmonic band limit parameter for integrand function [input].
alpha:Basis dilation factor for integrand function [input].
a_lim: Lower limit to evalue definite integral for [input].
b_lim: Upper limit to evalue definite integral for [input].
integral: Value of the evaluated integral [output].

Author: J. D. McEwen

Version: 0.1 June 2005